diff --git a/src/USER-MISC/pair_e3b.cpp b/src/USER-MISC/pair_e3b.cpp index 1f7e9fa3bb..a904939e42 100644 --- a/src/USER-MISC/pair_e3b.cpp +++ b/src/USER-MISC/pair_e3b.cpp @@ -33,28 +33,28 @@ //these are defined here to avoid confusing hardcoded indices, but //they do not allow flexibility. If they are changed the code will break #define DIM 3 -#define NUMH 2 //number of hydrogen atoms per water molecule -#define NUMO 2 //number of oxygen atoms per pair of water molecules -#define BOND_DELTA 1.01 //buffer for OH bonds that aren't perfectly rigid +#define NUMH 2 //number of hydrogen atoms per water molecule +#define NUMO 2 //number of oxygen atoms per pair of water molecules +#define BOND_DELTA 1.01 //buffer for OH bonds that aren't perfectly rigid -static constexpr double NOT_SET=1.024e300; +static constexpr double NOT_SET = 1.024e300; using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairE3B::PairE3B(LAMMPS *lmp) : Pair(lmp),pairPerAtom(10) +PairE3B::PairE3B(LAMMPS *lmp) : Pair(lmp), pairPerAtom(10) { single_enable = 0; restartinfo = 0; one_coeff = 1; - nextra=4; //store and tally pot energy terms eA, eB, eC, and e2 + nextra = 4; //store and tally pot energy terms eA, eB, eC, and e2 pvector = new double[nextra]; allocatedE3B = false; - pairO = nullptr; - pairH = nullptr; - exps = nullptr; - del3 = nullptr; + pairO = nullptr; + pairH = nullptr; + exps = nullptr; + del3 = nullptr; fpair3 = nullptr; sumExp = nullptr; } @@ -87,23 +87,22 @@ PairE3B::~PairE3B() void PairE3B::compute(int eflag, int vflag) { - int i,j,k,h,ii,jj,hh,kk,inum,jnum,otherO; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair,rsq,tmpexp; - double fxtmp,fytmp,fztmp,fix,fiy,fiz; - double delxh,delyh,delzh,rsqh,tmpr; - double scFact1,scFact2,scEng,scDer; - int *ilist,*jlist,*numneigh,**firstneigh; + int i, j, k, h, ii, jj, hh, kk, inum, jnum, otherO; + double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair, rsq, tmpexp; + double fxtmp, fytmp, fztmp, fix, fiy, fiz; + double delxh, delyh, delzh, rsqh, tmpr; + double scFact1, scFact2, scEng, scDer; + int *ilist, *jlist, *numneigh, **firstneigh; bool addedH; - if (natoms != atom->natoms) - error->all(FLERR,"pair E3B requires a fixed number of atoms"); + if (natoms != atom->natoms) error->all(FLERR, "pair E3B requires a fixed number of atoms"); - ev_init(eflag,vflag); + ev_init(eflag, vflag); //clear sumExp array - memset(sumExp,0.0,nbytes); + memset(sumExp, 0.0, nbytes); evdwl = 0.0; - pvector[0]=pvector[1]=pvector[2]=pvector[3]=0.0; + pvector[0] = pvector[1] = pvector[2] = pvector[3] = 0.0; double **x = atom->x; double **f = atom->f; @@ -120,8 +119,7 @@ void PairE3B::compute(int eflag, int vflag) // loop over half neighbor list of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; - if (type[i]!=typeO) - continue; + if (type[i] != typeO) continue; xtmp = x[i][0]; ytmp = x[i][1]; @@ -130,42 +128,41 @@ void PairE3B::compute(int eflag, int vflag) // two-body interactions jlist = firstneigh[i]; - jnum = numneigh[i]; + jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; //skip unless O-O interaction - if (type[j]!=typeO) - continue; + if (type[j] != typeO) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; //OO distance + rsq = delx * delx + dely * dely + delz * delz; //OO distance //two body interaction //not shifted b/c k2=4.87/A, so at cutoff (5.2A) e^(-kr) = 1e-11 if (rsq < rc2sq) { - tmpr = sqrt(rsq); - tmpexp = e2 * exp(-k2*tmpr); - fpair = k2 * tmpexp / tmpr; + tmpr = sqrt(rsq); + tmpexp = e2 * exp(-k2 * tmpr); + fpair = k2 * tmpexp / tmpr; - fxtmp = delx*fpair; - fytmp = dely*fpair; - fztmp = delz*fpair; - fix += fxtmp; - fiy += fytmp; - fiz += fztmp; + fxtmp = delx * fpair; + fytmp = dely * fpair; + fztmp = delz * fpair; + fix += fxtmp; + fiy += fytmp; + fiz += fztmp; f[j][0] -= fxtmp; f[j][1] -= fytmp; f[j][2] -= fztmp; if (evflag) { - ev_tally(i,j,nlocal,newton_pair,tmpexp,0.0,fpair,delx,dely,delz); + ev_tally(i, j, nlocal, newton_pair, tmpexp, 0.0, fpair, delx, dely, delz); pvector[0] += tmpexp; } - } //end if rsqmap(tag[otherO]+hh+1); + for (kk = 0; kk < NUMO; kk++) { + k = pairO[npair][kk]; + otherO = pairO[npair][(kk + 1) % 2]; + for (hh = 0; hh < NUMH; hh++) { + h = atom->map(tag[otherO] + hh + 1); //if hydrogen atom is missing, bond potential or shake will //catch this, so don't need to check here //if (h<0) // error->one(FLERR,"hydrogen atom missing"); - h = domain->closest_image(otherO,h); + h = domain->closest_image(otherO, h); pairH[npair][kk][hh] = h; delxh = x[k][0] - x[h][0]; delyh = x[k][1] - x[h][1]; delzh = x[k][2] - x[h][2]; - rsqh = delxh*delxh + delyh*delyh + delzh*delzh; + rsqh = delxh * delxh + delyh * delyh + delzh * delzh; if (rsqh < rc3sq) { - tmpr = sqrt(rsqh); - tmpexp = exp(-k3*tmpr); + tmpr = sqrt(rsqh); + tmpexp = exp(-k3 * tmpr); if (tmpr > rs) { - scFact1 = rc3-tmpr; - scFact2 = sc_num + 2*tmpr; - scEng = scFact1*scFact1*scFact2*sc_denom; - scDer = k3*scEng - 6*scFact1*(rs-tmpr)*sc_denom; - } else { + scFact1 = rc3 - tmpr; + scFact2 = sc_num + 2 * tmpr; + scEng = scFact1 * scFact1 * scFact2 * sc_denom; + scDer = k3 * scEng - 6 * scFact1 * (rs - tmpr) * sc_denom; + } else { scDer = k3; scEng = 1.0; } //need to keep fpair3 separate from del3 for virial - fpair3[npair][kk][hh] = scDer*tmpexp/tmpr; - tmpexp *= scEng; - exps[npair][kk][hh] = tmpexp; + fpair3[npair][kk][hh] = scDer * tmpexp / tmpr; + tmpexp *= scEng; + exps[npair][kk][hh] = tmpexp; del3[npair][kk][hh][0] = delxh; del3[npair][kk][hh][1] = delyh; del3[npair][kk][hh][2] = delzh; //accumulate global vector of sum(e^kr) //tags start at 1, so subtract one to index sumExp - sumExp[tag[k]-1] += tmpexp; - sumExp[tag[h]-1] += tmpexp; + sumExp[tag[k] - 1] += tmpexp; + sumExp[tag[h] - 1] += tmpexp; addedH = true; } else { - exps [npair][kk][hh] = 0.0; + exps[npair][kk][hh] = 0.0; fpair3[npair][kk][hh] = 0.0; - } //if < rc3sq - } //end loop through 2 Hs - } //end for kk in NUMO + } //if < rc3sq + } //end loop through 2 Hs + } //end for kk in NUMO //if added a pair, check if array is too small and grow if (addedH) { npair++; - if (npair >= pairmax) - error->one(FLERR,"neigh is too small"); + if (npair >= pairmax) error->one(FLERR, "neigh is too small"); } - } //end if < rc3deltaSq - } //end for jj neigh + } //end if < rc3deltaSq + } //end for jj neigh //add 2-body forces on i f[i][0] += fix; f[i][1] += fiy; f[i][2] += fiz; - } //end for ii + } //end for ii //communicate sumExp array //tested that no change in speed with MPI_IN_PLACE - MPI_Allreduce(MPI_IN_PLACE,sumExp,maxID,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(MPI_IN_PLACE, sumExp, maxID, MPI_DOUBLE, MPI_SUM, world); //now loop through list of pairs, calculating 3body forces - int j2,otherH; - double partA,partB,partC; + int j2, otherH; + double partA, partB, partC; for (ii = 0; ii < npair; ii++) { - for (kk=0; kkntypes; - memory->create(setflag,n+1,n+1,"pair:setflag"); - memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); } void PairE3B::allocateE3B() @@ -349,28 +344,26 @@ void PairE3B::allocateE3B() allocatedE3B = true; //TODO: get memory->grow working for 4d arrays - pairmax = atom->nlocal*pairPerAtom; //initial guess for size of pair lists - memory->create(pairO ,pairmax,NUMO ,"pair:pairO"); - memory->create(pairH ,pairmax,NUMO,NUMH ,"pair:pairH"); - memory->create(exps ,pairmax,NUMO,NUMH ,"pair:exps"); - memory->create(fpair3,pairmax,NUMO,NUMH ,"pair:fpair3"); - memory->create(del3 ,pairmax,NUMO,NUMH,DIM,"pair:del3"); + pairmax = atom->nlocal * pairPerAtom; //initial guess for size of pair lists + memory->create(pairO, pairmax, NUMO, "pair:pairO"); + memory->create(pairH, pairmax, NUMO, NUMH, "pair:pairH"); + memory->create(exps, pairmax, NUMO, NUMH, "pair:exps"); + memory->create(fpair3, pairmax, NUMO, NUMH, "pair:fpair3"); + memory->create(del3, pairmax, NUMO, NUMH, DIM, "pair:del3"); //set del3 to zero to silence valgrind memcheck errors //don't need to do this in every call to compute() because we set //exps and fpair3 to zero, and all uses of del3 are multiplied by one of those - int ii,jj,kk,ll; - for (ii=0; iinatoms; - maxID=find_maxID(); - if (!natoms) - error->all(FLERR,"No atoms found"); - memory->create(sumExp,maxID,"pair:sumExp"); + natoms = atom->natoms; + maxID = find_maxID(); + if (!natoms) error->all(FLERR, "No atoms found"); + memory->create(sumExp, maxID, "pair:sumExp"); nbytes = sizeof(double) * maxID; } @@ -380,8 +373,8 @@ void PairE3B::allocateE3B() void PairE3B::settings(int narg, char **arg) { - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); - typeO=utils::inumeric(FLERR,arg[0],false,lmp); + if (narg != 1) error->all(FLERR, "Illegal pair_style command"); + typeO = utils::inumeric(FLERR, arg[0], false, lmp); } /* ---------------------------------------------------------------------- @@ -393,67 +386,66 @@ void PairE3B::coeff(int narg, char **arg) if (!allocated) allocate(); //1=* 2=* 3/4=1st keyword/value - if (narg < 4) - error->all(FLERR,"There must be at least one keyword given to pair_coeff"); + if (narg < 4) error->all(FLERR, "There must be at least one keyword given to pair_coeff"); - // clear setflag since coeff() called once with I,J = * * + // clear setflag since coeff() called once with I,J = * * int n = atom->ntypes; for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int j = i; j <= n; j++) setflag[i][j] = 0; - setflag[typeO][typeO]=1; + setflag[typeO][typeO] = 1; //parse keyword/value pairs - double bondL=0.0; //OH bond length - bool repeatFlag=false; + double bondL = 0.0; //OH bond length + bool repeatFlag = false; //clear parameters - e2=ea=eb=ec=k3=k2=NOT_SET; - rs=rc3=rc2=0.0; + e2 = ea = eb = ec = k3 = k2 = NOT_SET; + rs = rc3 = rc2 = 0.0; - int iarg = 2; //beginning of keyword/value pairs + int iarg = 2; //beginning of keyword/value pairs while (iarg < narg) { char *keyword = arg[iarg++]; - if (checkKeyword(keyword,"Ea",1,narg-iarg)) - ea=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"Eb",1,narg-iarg)) - eb=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"Ec",1,narg-iarg)) - ec=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"K3",1,narg-iarg)) - k3=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"Rs",1,narg-iarg)) - rs=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"Rc3",1,narg-iarg)) - rc3=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"Rc2",1,narg-iarg)) - rc2=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"bondL",1,narg-iarg)) - bondL=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"E2",1,narg-iarg)) - e2=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"K2",1,narg-iarg)) - k2=utils::numeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"neigh",1,narg-iarg)) - pairPerAtom=utils::inumeric(FLERR,arg[iarg++],false,lmp); - else if (checkKeyword(keyword,"preset",1,narg-iarg)) { - int presetFlag=utils::inumeric(FLERR,arg[iarg++],false,lmp); - presetParam(presetFlag,repeatFlag,bondL); - } else error->all(FLERR,"Keyword {} is unknown",keyword); + if (checkKeyword(keyword, "Ea", 1, narg - iarg)) + ea = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "Eb", 1, narg - iarg)) + eb = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "Ec", 1, narg - iarg)) + ec = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "K3", 1, narg - iarg)) + k3 = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "Rs", 1, narg - iarg)) + rs = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "Rc3", 1, narg - iarg)) + rc3 = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "Rc2", 1, narg - iarg)) + rc2 = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "bondL", 1, narg - iarg)) + bondL = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "E2", 1, narg - iarg)) + e2 = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "K2", 1, narg - iarg)) + k2 = utils::numeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "neigh", 1, narg - iarg)) + pairPerAtom = utils::inumeric(FLERR, arg[iarg++], false, lmp); + else if (checkKeyword(keyword, "preset", 1, narg - iarg)) { + int presetFlag = utils::inumeric(FLERR, arg[iarg++], false, lmp); + presetParam(presetFlag, repeatFlag, bondL); + } else + error->all(FLERR, "Keyword {} is unknown", keyword); } checkInputs(bondL); //cutmax for neighbor listing - cutmax = MAX(rc2,rc3); - rc2sq = rc2*rc2; - rc3sq = rc3*rc3; - rc3deltaSq = (rc3+bondL)*(rc3+bondL); + cutmax = MAX(rc2, rc3); + rc2sq = rc2 * rc2; + rc3sq = rc3 * rc3; + rc3deltaSq = (rc3 + bondL) * (rc3 + bondL); - double tmpfact=1.0/(rc3-rs); - sc_denom=tmpfact*tmpfact*tmpfact; - sc_num=rc3-3*rs; + double tmpfact = 1.0 / (rc3 - rs); + sc_denom = tmpfact * tmpfact * tmpfact; + sc_num = rc3 - 3 * rs; } /* ---------------------------------------------------------------------- @@ -462,132 +454,131 @@ void PairE3B::coeff(int narg, char **arg) void PairE3B::init_style() { - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style E3B requires atom IDs"); - if (force->newton_pair == 0) - error->all(FLERR,"Pair style E3B requires newton pair on"); - if (typeO<1 || typeO>atom->ntypes) - error->all(FLERR,"Invalid Otype: out of bounds"); + if (atom->tag_enable == 0) error->all(FLERR, "Pair style E3B requires atom IDs"); + if (force->newton_pair == 0) error->all(FLERR, "Pair style E3B requires newton pair on"); + if (typeO < 1 || typeO > atom->ntypes) error->all(FLERR, "Invalid Otype: out of bounds"); // need a half neighbor list - neighbor->request(this,instance_me); + neighbor->request(this, instance_me); - if (!force->pair_match("tip4p",false,0)) - if (comm->me==0) error->warning(FLERR,"E3B pair_style is designed for use with hybrid/overlay tip4p style"); + if (!force->pair_match("tip4p", false, 0)) + if (comm->me == 0) + error->warning(FLERR, "E3B pair_style is designed for use with hybrid/overlay tip4p style"); if (!allocatedE3B) allocateE3B(); } - static const char cite_E3B1[] = - "Explicit Three-Body (E3B) potential for water:\n\n" - "@article{kumar_water_2008,\n" - "title = {Water Simulation Model with Explicit Three-Molecule Interactions},\n" - "volume = {112},\n" - "doi = {10.1021/jp8009468},\n" - "number = {28},\n" - "journal = {J Phys. Chem. B},\n" - "author = {Kumar, R. and Skinner, J. L.},\n" - "year = {2008},\n" - "pages = {8311--8318}\n" - "}\n\n"; + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{kumar_water_2008,\n" + "title = {Water Simulation Model with Explicit Three-Molecule Interactions},\n" + "volume = {112},\n" + "doi = {10.1021/jp8009468},\n" + "number = {28},\n" + "journal = {J Phys. Chem. B},\n" + "author = {Kumar, R. and Skinner, J. L.},\n" + "year = {2008},\n" + "pages = {8311--8318}\n" + "}\n\n"; static const char cite_E3B2[] = - "Explicit Three-Body (E3B) potential for water:\n\n" - "@article{tainter_robust_2011,\n" - "title = {Robust three-body water simulation model},\n" - "volume = {134},\n" - "doi = {10.1063/1.3587053},\n" - "number = {18},\n" - "journal = {J. Chem. Phys},\n" - "author = {Tainter, C. J. and Pieniazek, P. A. and Lin, Y.-S. and Skinner, J. L.},\n" - "year = {2011},\n" - "pages = {184501}\n" - "}\n\n"; + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{tainter_robust_2011,\n" + "title = {Robust three-body water simulation model},\n" + "volume = {134},\n" + "doi = {10.1063/1.3587053},\n" + "number = {18},\n" + "journal = {J. Chem. Phys},\n" + "author = {Tainter, C. J. and Pieniazek, P. A. and Lin, Y.-S. and Skinner, J. L.},\n" + "year = {2011},\n" + "pages = {184501}\n" + "}\n\n"; static const char cite_E3B3[] = - "Explicit Three-Body (E3B) potential for water:\n\n" - "@article{tainter_reparametrized_2015,\n" - "title = {Reparametrized {E3B} (Explicit Three-Body) Water Model Using the {TIP4P/2005} Model as a Reference},\n" - "volume = {11},\n" - "doi = {10.1021/acs.jctc.5b00117},\n" - "number = {5},\n" - "journal = {J. Chem. Theory Comput.},\n" - "author = {Tainter, Craig J. and Shi, Liang and Skinner, James L.},\n" - "year = {2015},\n" - "pages = {2268--2277}\n" - "}\n\n"; + "Explicit Three-Body (E3B) potential for water:\n\n" + "@article{tainter_reparametrized_2015,\n" + "title = {Reparametrized {E3B} (Explicit Three-Body) Water Model Using the {TIP4P/2005} Model " + "as a Reference},\n" + "volume = {11},\n" + "doi = {10.1021/acs.jctc.5b00117},\n" + "number = {5},\n" + "journal = {J. Chem. Theory Comput.},\n" + "author = {Tainter, Craig J. and Shi, Liang and Skinner, James L.},\n" + "year = {2015},\n" + "pages = {2268--2277}\n" + "}\n\n"; -void PairE3B::presetParam(const int flag,bool &repeatFlag,double &bondL) { - if (repeatFlag) - error->all(FLERR,"Cannot request two different sets of preset parameters"); - repeatFlag=true; +void PairE3B::presetParam(const int flag, bool &repeatFlag, double &bondL) +{ + if (repeatFlag) error->all(FLERR, "Cannot request two different sets of preset parameters"); + repeatFlag = true; - if (ea!=NOT_SET || eb!=NOT_SET || ec!=NOT_SET || - e2!=NOT_SET || k3!=NOT_SET || k2!=NOT_SET || - bondL!=0.0 || rs!=0.0 || rc3!=0.0 || rc2!=0.0 ) - error->all(FLERR,"Preset keyword will overwrite another keyword setting"); + if (ea != NOT_SET || eb != NOT_SET || ec != NOT_SET || e2 != NOT_SET || k3 != NOT_SET || + k2 != NOT_SET || bondL != 0.0 || rs != 0.0 || rc3 != 0.0 || rc2 != 0.0) + error->all(FLERR, "Preset keyword will overwrite another keyword setting"); - double econv=1.0,lconv=1.0; - if (strcmp(update->unit_style,"real") == 0) { - econv=1.0/4.184; - lconv=1.0; - } else if (strcmp(update->unit_style,"metal") == 0) { - econv=0.103653271; - lconv=1.0; - } else if (strcmp(update->unit_style,"si") == 0) { - econv=1.660578e-21; - lconv=1e-10; - } else if (strcmp(update->unit_style,"cgs") == 0) { - econv=1.660578e-14; - lconv=1e-8; - } else error->all(FLERR, "Pre-defined E3B parameters have not been set for {} units.", - update->unit_style); + double econv = 1.0, lconv = 1.0; + if (strcmp(update->unit_style, "real") == 0) { + econv = 1.0 / 4.184; + lconv = 1.0; + } else if (strcmp(update->unit_style, "metal") == 0) { + econv = 0.103653271; + lconv = 1.0; + } else if (strcmp(update->unit_style, "si") == 0) { + econv = 1.660578e-21; + lconv = 1e-10; + } else if (strcmp(update->unit_style, "cgs") == 0) { + econv = 1.660578e-14; + lconv = 1e-8; + } else + error->all(FLERR, "Pre-defined E3B parameters have not been set for {} units.", + update->unit_style); //here parameters are defined in kJ/mol and A //they will be converted to the lammps units after - if (flag==2008) { - error->all(FLERR,"\"preset 2008\" is not yet supported, because this would require distinct k3 coefficients, " - "use \"preset 2011\" or \"preset 2015\""); + if (flag == 2008) { + error->all(FLERR, + "'preset 2008' is not yet supported, because this would require distinct k3 " + "coefficients, use 'preset 2011' or 'preset 2015'"); if (lmp->citeme) lmp->citeme->add(cite_E3B1); - ea = 4699.6; - eb =-2152.9; - ec = 1312.7; + ea = 4699.6; + eb = -2152.9; + ec = 1312.7; //ka = 1.0/1.88; //kb = 1.0/1.71; //kc = 1.0/1.56; - e2 = 1.925e6; - k2 = 4.67; - rs = 5.0; + e2 = 1.925e6; + k2 = 4.67; + rs = 5.0; rc3 = 5.2; rc2 = 5.2; bondL = 0.9572; - } else if (flag==2011) { + } else if (flag == 2011) { if (lmp->citeme) lmp->citeme->add(cite_E3B2); - ea = 1745.7; - eb =-4565.0; - ec = 7606.8; - k3 = 1.907; - e2 = 2.349e6; - k2 = 4.872; - rs = 5.0; + ea = 1745.7; + eb = -4565.0; + ec = 7606.8; + k3 = 1.907; + e2 = 2.349e6; + k2 = 4.872; + rs = 5.0; rc3 = 5.2; rc2 = 5.2; bondL = 0.9572; - } else if (flag==2015) { + } else if (flag == 2015) { if (lmp->citeme) lmp->citeme->add(cite_E3B3); - ea = 150.0; - eb =-1005.0; - ec = 1880.0; - k3 = 1.907; - e2 = 0.453e6; - k2 = 4.872; - rs = 5.0; + ea = 150.0; + eb = -1005.0; + ec = 1880.0; + k3 = 1.907; + e2 = 0.453e6; + k2 = 4.872; + rs = 5.0; rc3 = 5.2; rc2 = 5.2; bondL = 0.9572; } else - error->all(FLERR,"Unknown argument: preset only takes 2011 or 2015 as arguments"); + error->all(FLERR, "Unknown argument: preset only takes 2011 or 2015 as arguments"); //convert units ea *= econv; @@ -599,7 +590,7 @@ void PairE3B::presetParam(const int flag,bool &repeatFlag,double &bondL) { rs *= lconv; rc2 *= lconv; rc3 *= lconv; - bondL *= lconv*BOND_DELTA; + bondL *= lconv * BOND_DELTA; } /* ---------------------------------------------------------------------- @@ -608,14 +599,15 @@ void PairE3B::presetParam(const int flag,bool &repeatFlag,double &bondL) { //pair.cpp::init uses this to set cutsq array, used for neighboring, etc double PairE3B::init_one(int i, int j) { - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); return cutmax; } -bool PairE3B::checkKeyword(const char *thiskey,const char *test,const int nVal, const int nRem) { - if (strcmp(thiskey,test) == 0) { - if (nRemall(FLERR,"Too few arguments to '{}' keyword.",test); +bool PairE3B::checkKeyword(const char *thiskey, const char *test, const int nVal, const int nRem) +{ + if (strcmp(thiskey, test) == 0) { + if (nRem < nVal) error->all(FLERR, "Too few arguments to '{}' keyword.", test); return true; } return false; @@ -630,46 +622,34 @@ tagint PairE3B::find_maxID() tagint max = 0; tagint maxID; - for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); - MPI_Allreduce(&max,&maxID,1,MPI_LMP_TAGINT,MPI_MAX,world); + for (int i = 0; i < nlocal; i++) max = MAX(max, tag[i]); + MPI_Allreduce(&max, &maxID, 1, MPI_LMP_TAGINT, MPI_MAX, world); return maxID; } -void PairE3B::checkInputs(const double &bondL) { +void PairE3B::checkInputs(const double &bondL) +{ //first check that all necessary values were set - if (rc2==0.0) - error->all(FLERR,"rc2 keyword missing"); - if (rs==0.0) - error->all(FLERR,"Rs keyword missing"); - if (rc3==0.0) - error->all(FLERR,"Rc3 keyword missing"); - if (bondL==0.0) - error->all(FLERR,"bondL keyword missing"); - if (ea==NOT_SET) - error->all(FLERR,"Ea keyword missing"); - if (eb==NOT_SET) - error->all(FLERR,"Eb keyword missing"); - if (ec==NOT_SET) - error->all(FLERR,"Ec keyword missing"); - if (k3==NOT_SET) - error->all(FLERR,"K3 keyword missing"); - if (e2==NOT_SET) - error->all(FLERR,"E2 keyword missing"); - if (k2==NOT_SET) - error->all(FLERR,"K2 keyword missing"); + if (rc2 == 0.0) error->all(FLERR, "rc2 keyword missing"); + if (rs == 0.0) error->all(FLERR, "Rs keyword missing"); + if (rc3 == 0.0) error->all(FLERR, "Rc3 keyword missing"); + if (bondL == 0.0) error->all(FLERR, "bondL keyword missing"); + if (ea == NOT_SET) error->all(FLERR, "Ea keyword missing"); + if (eb == NOT_SET) error->all(FLERR, "Eb keyword missing"); + if (ec == NOT_SET) error->all(FLERR, "Ec keyword missing"); + if (k3 == NOT_SET) error->all(FLERR, "K3 keyword missing"); + if (e2 == NOT_SET) error->all(FLERR, "E2 keyword missing"); + if (k2 == NOT_SET) error->all(FLERR, "K2 keyword missing"); //now test that values are within acceptable ranges - if (k2 < 0.0 or k3 < 0.0) - error->all(FLERR,"exponential decay is negative"); - if (bondL<0.0) - error->all(FLERR,"OH bond length is negative"); - if (rc2 < 0.0 || rc3 < 0.0 || rs < 0.0) - error->all(FLERR,"potential cutoff is negative"); - if (rs > rc3) - error->all(FLERR,"potential switching distance is larger than cutoff"); - if (rs==rc3) - error->warning(FLERR,"potential switching distance is equal to cutoff: this is untested and not conserve energy"); - if (pairPerAtom<0) - error->all(FLERR,"neigh is negative"); + if (k2 < 0.0 or k3 < 0.0) error->all(FLERR, "exponential decay is negative"); + if (bondL < 0.0) error->all(FLERR, "OH bond length is negative"); + if (rc2 < 0.0 || rc3 < 0.0 || rs < 0.0) error->all(FLERR, "potential cutoff is negative"); + if (rs > rc3) error->all(FLERR, "potential switching distance is larger than cutoff"); + if (rs == rc3) + error->warning(FLERR, + "potential switching distance is equal to cutoff: this is untested and not " + "conserve energy"); + if (pairPerAtom < 0) error->all(FLERR, "neigh is negative"); }