jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14449 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2016-01-15 15:53:23 +00:00
parent bfa988207e
commit 984a1afa72
7 changed files with 114 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
src/USER-SMTBQ/pair_smtbq.cpp
View File

@ -38,10 +38,10 @@
<http://www.gnu.org/licenses/>.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_smtbq.h"
#include "atom.h"
#include "comm.h"
@ -147,7 +147,6 @@ PairSMTBQ::PairSMTBQ(LAMMPS *lmp) : Pair(lmp)
fct = NULL;
maxpage = 0;
// set comm size needed by this Pair
@ -258,6 +257,9 @@ void PairSMTBQ::coeff(int narg, char **arg)
if (!allocated) allocate();
if (strstr(force->pair_style,"hybrid"))
error->all(FLERR,"Pair style SMTBQ is not compatible with hybrid styles");
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
@ -460,11 +462,11 @@ void PairSMTBQ::read_file(char *file)
for (i=1; i <= maxintparam; i++)
intparams[i].mode = (char*) malloc(sizeof(char)*6);
QEqMode = (char*) malloc(sizeof(char)*18);
Bavard = (char*) malloc(sizeof(char)*5);
QInitMode = (char*) malloc(sizeof(char)*18);
writepot = (char*) malloc(sizeof(char)*5);
writeenerg = (char*) malloc(sizeof(char)*5);
QEqMode = (char*) malloc(sizeof(char)*19);
Bavard = (char*) malloc(sizeof(char)*6);
QInitMode = (char*) malloc(sizeof(char)*19);
writepot = (char*) malloc(sizeof(char)*6);
writeenerg = (char*) malloc(sizeof(char)*6);
// Little loop for ion's parameters
@ -844,7 +846,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
int *ilist,*jlist,*numneigh,**firstneigh;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double rsq,iq,jq,Eself,natom;
double rsq,iq,jq,Eself;
double ecovtot,ErepOO,ErepMO,Eion,Ecoh;
double **tmp,**tmpAll,*nmol;
double dq,dqcov;
@ -891,13 +893,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
#ifdef LAMMPS_BIGBIG
long int *tag = atom->tag;
#else
tagint *tag = atom->tag;
#endif
int *type = atom->type;
int newton_pair = force->newton_pair;
int nlocal = atom->nlocal;
@ -950,10 +946,6 @@ void PairSMTBQ::compute(int eflag, int vflag)
iq = q[i];
gp = flag_QEq[i];
if (gp == 0 && itype > 0) natom += 1.0;
// if (gp == 0 && itype > 0) nmol[gp] += 1.0;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -999,6 +991,7 @@ void PairSMTBQ::compute(int eflag, int vflag)
for (jj = 0; jj < jnum; jj++) {
// ===============================
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
jtag = tag[j]; jq = q[j];
@ -2345,7 +2338,9 @@ void PairSMTBQ::QForce_charge(int loop)
if (loop == 0) {
// ==================
memset(sbcov,0,sizeof(double)*atom->nmax);
memset(coord,0,sizeof(double)*atom->nmax);
memset(sbmet,0,sizeof(double)*atom->nmax);
for (ii = 0; ii < inum; ii ++) {
//--------------------------------
@ -2354,8 +2349,6 @@ void PairSMTBQ::QForce_charge(int loop)
gp = flag_QEq[i];
sbcov[i] =coord[i]= sbmet[i] = 0.0;
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
@ -2454,6 +2447,7 @@ void PairSMTBQ::QForce_charge(int loop)
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
m = intype[itype][jtype];
jq = q[j];
@ -2897,16 +2891,17 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
{
int ii,i,jj,j,kk,k,itype,jtype,ktype,jnum,m,gp,zz,z,kgp;
int iproc,team_elt[10][nproc],team_QEq[10][nproc][5];
int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp,nboite;
int *ilist,*jlist,*numneigh,**firstneigh,ngp,igp;
double delr[3],xtmp,ytmp,ztmp,rsq;
int **flag_gp, *nelt, **tab_gp;
int QEq,QEqall[nproc];
double **x = atom->x;
int *type = atom->type;
const int nlocal = atom->nlocal;
const int nghost = atom->nghost;
const int nall = nlocal + nghost;
int inum = list->inum;
int nlocal = atom->nlocal;
int nghost = atom->nghost;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -2916,7 +2911,6 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// On declare et initialise nos p'tits tableaux
// +++++++++++++++++++++++++++++++++++++++++++++++++
nboite = nlocal + nghost;
int **tabtemp,**Alltabtemp, *gptmp, *Allgptmp;
memory->create(tabtemp,10*nproc+10,nproc,"pair:tabtemp");
@ -2924,19 +2918,19 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
memory->create(gptmp,10*nproc+10,"pair:gptmp");
memory->create(Allgptmp,10*nproc+10,"pair:Allgptmp");
memory->create(flag_gp,nproc,nboite,"pair:flag_gp");
memory->create(nelt,nboite,"pair:nelt");
memory->create(tab_gp,10,nboite,"pair:flag_gp");
memory->create(flag_gp,nproc,nall,"pair:flag_gp");
memory->create(nelt,nall,"pair:nelt");
memory->create(tab_gp,10,nall,"pair:flag_gp");
for (i = 0; i < nlocal+nghost ; i++) { flag_QEq[i] = 0; }
for (i = 0; i < nall ; i++) { flag_QEq[i] = 0; }
for (i = 0; i < 10*nproc; i++) {
gptmp[i] = 0; Allgptmp[i] = 0;
for (j=0;j<nproc;j++) { tabtemp[i][j] = 0;
Alltabtemp[i][j] = 0;}
}
for (i = 0; i < 10; i++) {
for (k = 0; k < nboite; k++) { tab_gp[i][k] = 0;
for (k = 0; k < nall; k++) { tab_gp[i][k] = 0;
if (i == 0) nelt[k] = 0;
}
for (j = 0; j < nproc; j++) {
@ -2949,7 +2943,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// printf ("groupeQEq me %d - nloc %d nghost %d boite %d\n",
// me,nlocal,nghost,nboite);
// me,nlocal,nghost,nall);
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++
// On identifie les atomes rentrant dans le schema QEq +
@ -2968,7 +2962,9 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++ )
{
j = jlist[jj] ; jtype = map[type[j]];
j = jlist[jj] ;
j &= NEIGHMASK;
jtype = map[type[j]];
if (jtype == itype) continue;
m = intype[itype][jtype];
@ -3007,7 +3003,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
for (m = 0; m < nproc; m++) {
for (i = 0; i < nboite; i++) { flag_gp[m][i] = 0; }
for (i = 0; i < nall; i++) { flag_gp[m][i] = 0; }
}
// OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
@ -3057,7 +3053,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// ---------------------------------------------
jlist = firstneigh[k];
jnum = numneigh[k];
for (j = 0; j < nboite; j++ )
for (j = 0; j < nall; j++ )
{
jtype = map[type[j]];
if (jtype == ktype) continue;
@ -3104,7 +3100,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
}
nelt[ngp] = 0;
for (z = nlocal; z < nboite; z++) {
for (z = nlocal; z < nall; z++) {
if (flag_gp[me][z] == ngp) flag_gp[me][z] = igp;
}
@ -3152,7 +3148,9 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++ )
{
j = jlist[jj] ; jtype = map[type[j]];
j = jlist[jj] ;
j &= NEIGHMASK;
jtype = map[type[j]];
if (jtype != 0) continue;
m = 0;
@ -3192,7 +3190,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
}
nelt[kgp] = 0;
for (k = 0; k < nboite; k++) {
for (k = 0; k < nall; k++) {
if (flag_gp[me][k] == kgp) flag_gp[me][k] = igp;
}
@ -3223,7 +3221,7 @@ void PairSMTBQ::groupQEqAllParallel_QEq()
// =============== End of COMM =================
for (i = 0; i < nboite; i++) {
for (i = 0; i < nall; i++) {
m = 10*me + flag_gp[me][i];
if (m == 10*me) continue; // Pas de groupe zero
@ -3553,7 +3551,7 @@ int PairSMTBQ::Tokenize( char* s, char*** tok )
mot = NULL;
strncpy( test, s, MAXLINE );
strncpy( test, s, MAXLINE-1 );
for( mot = strtok(test, sep); mot; mot = strtok(NULL, sep) ) {
strncpy( (*tok)[count], mot, MAXLINE );

15
src/create_atoms.cpp
View File

@ -416,14 +416,19 @@ void CreateAtoms::command(int narg, char **arg)
}
// for MOLECULE mode:
// set molecule IDs for created atoms if used
// reset new molecule bond,angle,etc and special values
// molecule can mean just a mol ID or bonds/angles/etc or mol templates
// set molecule IDs for created atoms if atom->molecule_flag is set
// reset new molecule bond,angle,etc and special values if defined
// send atoms to new owning procs via irregular comm
// since not all atoms I created will be within my sub-domain
// perform special list build if needed
if (mode == MOLECULE) {
int molecule_flag = atom->molecule_flag;
int molecular = atom->molecular;
tagint *molecule = atom->molecule;
// molcreate = # of molecules I created
int molcreate = (atom->nlocal - nlocal_previous) / onemol->natoms;
@ -444,8 +449,7 @@ void CreateAtoms::command(int narg, char **arg)
// including molecules existing before this creation
tagint moloffset;
tagint *molecule = atom->molecule;
if (molecule) {
if (molecule_flag) {
tagint max = 0;
for (int i = 0; i < nlocal_previous; i++) max = MAX(max,molecule[i]);
tagint maxmol;
@ -480,13 +484,12 @@ void CreateAtoms::command(int narg, char **arg)
tagint **improper_atom4 = atom->improper_atom4;
int **nspecial = atom->nspecial;
tagint **special = atom->special;
int molecular = atom->molecular;
int ilocal = nlocal_previous;
for (int i = 0; i < molcreate; i++) {
if (tag) offset = tag[ilocal]-1;
for (int m = 0; m < natoms; m++) {
if (molecular) molecule[ilocal] = moloffset + i+1;
if (molecule_flag) molecule[ilocal] = moloffset + i+1;
if (molecular == 2) {
atom->molindex[ilocal] = 0;
atom->molatom[ilocal] = m;

13
src/fix_nh.cpp
View File

@ -78,6 +78,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
omega_mass_flag = 0;
etap_mass_flag = 0;
flipflag = 1;
dipole_flag = 0;
// turn on tilt factor scaling, whenever applicable
@ -321,6 +322,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"update") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
else error->all(FLERR,"Illegal fix nvt/npt/nph command");
iarg += 2;
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
@ -417,6 +423,13 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
p_period[0] != p_period[2]))
error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
if (dipole_flag) {
if (!atom->sphere_flag)
error->all(FLERR,"Using update dipole flag requires atom style sphere");
if (!atom->mu_flag)
error->all(FLERR,"Using update dipole flag requires atom attribute mu");
}
if ((tstat_flag && t_period <= 0.0) ||
(p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) ||

1
src/fix_nh.h
View File

@ -109,6 +109,7 @@ class FixNH : public Fix {
int eta_mass_flag; // 1 if eta_mass updated, 0 if not.
int omega_mass_flag; // 1 if omega_mass updated, 0 if not.
int etap_mass_flag; // 1 if etap_mass updated, 0 if not.
int dipole_flag; // 1 if dipole is updated, 0 if not.
int scaleyz; // 1 if yz scaled with lz
int scalexz; // 1 if xz scaled with lz

36
src/fix_nh_sphere.cpp
View File

@ -91,6 +91,42 @@ void FixNHSphere::nve_v()
}
}
/* ----------------------------------------------------------------------
perform full-step update of position with dipole orientation, if requested
-----------------------------------------------------------------------*/
void FixNHSphere::nve_x()
{
// standard nve_x position update
FixNH::nve_x();
// update mu for dipoles
// d_mu/dt = omega cross mu
// renormalize mu to dipole length
if (dipole_flag) {
double msq,scale,g[3];
double **mu = atom->mu;
double **omega = atom->omega;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mu[i][3] > 0.0) {
g[0] = mu[i][0] + dtv * (omega[i][1]*mu[i][2]-omega[i][2]*mu[i][1]);
g[1] = mu[i][1] + dtv * (omega[i][2]*mu[i][0]-omega[i][0]*mu[i][2]);
g[2] = mu[i][2] + dtv * (omega[i][0]*mu[i][1]-omega[i][1]*mu[i][0]);
msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
scale = mu[i][3]/sqrt(msq);
mu[i][0] = g[0]*scale;
mu[i][1] = g[1]*scale;
mu[i][2] = g[2]*scale;
}
}
}
/* ----------------------------------------------------------------------
perform half-step scaling of rotatonal velocities
-----------------------------------------------------------------------*/

1
src/fix_nh_sphere.h
View File

@ -26,6 +26,7 @@ class FixNHSphere : public FixNH {
protected:
void nve_v();
void nve_x();
void nh_v_temp();
};

15
src/fix_property_atom.cpp
View File

@ -222,8 +222,21 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf,
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
if (values[0] == NULL) {
char str[128];
sprintf(str,"Too few lines in %s section of data file",keyword);
error->one(FLERR,str);
}
int format_ok = 1;
for (j = 1; j < nwords; j++) {
values[j] = strtok(NULL," \t\n\r\f");
if (values[j] == NULL) format_ok = 0;
}
if (!format_ok) {
char str[128];
sprintf(str,"Incorrect %s format in data file",keyword);
error->all(FLERR,str);
}
itag = ATOTAGINT(values[0]) + id_offset;
if (itag <= 0 || itag > map_tag_max) {