diff --git a/doc/src/pair_vashishta.txt b/doc/src/pair_vashishta.txt index c396087b30..5c6ccb0a88 100644 --- a/doc/src/pair_vashishta.txt +++ b/doc/src/pair_vashishta.txt @@ -8,26 +8,36 @@ pair_style vashishta command :h3 pair_style vashishta/omp command :h3 +pair_style vashishta/table command :h3 [Syntax:] -pair_style vashishta :pre +pair_style style args :pre + +style = {vashishta} or {vashishta/table} +args = list of arguments for a particular style :ul + {vashishta} args = none + {vashishta/table} args = Ntable cutinner + Ntable = # of tabulation points + cutinner = tablulate from cutinner to cutoff :pre [Examples:] pair_style vashishta pair_coeff * * SiC.vashishta Si C :pre +pair_style vashishta/table 100000 0.2 +pair_coeff * * SiC.vashishta Si C :pre + [Description:] -The {vashishta} style computes the combined 2-body and 3-body -family of potentials developed in the group of Vashishta and -co-workers. By combining repulsive, screened Coulombic, -screened charge-dipole, and dispersion interactions with a -bond-angle energy based on the Stillinger-Weber potential, -this potential has been used to describe a variety of inorganic -compounds, including SiO2 "Vashishta1990"_#Vashishta1990, -SiC "Vashishta2007"_#Vashishta2007, +The {vashishta} and {vashishta/table} styles compute the combined +2-body and 3-body family of potentials developed in the group of +Vashishta and co-workers. By combining repulsive, screened Coulombic, +screened charge-dipole, and dispersion interactions with a bond-angle +energy based on the Stillinger-Weber potential, this potential has +been used to describe a variety of inorganic compounds, including SiO2 +"Vashishta1990"_#Vashishta1990, SiC "Vashishta2007"_#Vashishta2007, and InP "Branicio2009"_#Branicio2009. The potential for the energy U of a system of atoms is @@ -43,10 +53,19 @@ both zero at {rc}. The summation over three-body terms is over all neighbors J and K within a cut-off distance = {r0}, where the exponential screening function becomes zero. -Only a single pair_coeff command is used with the {vashishta} style which -specifies a Vashishta potential file with parameters for all -needed elements. These are mapped to LAMMPS atom types by specifying -N additional arguments after the filename in the pair_coeff command, +The {vashishta} style computes these formulas analytically. The +{vashishta/table} style tabulates the analytic values for {Ntable} +points from cutinner to the cutoff of the potential. The points are +equally spaced in R^2 space from cutinner^2 to cutoff^2. For the +two-body term in the above equation, a linear interpolation for each +pairwise distance between adjacent points in the table. In practice +the tabulated version can run 3-5x faster than the analytic version +with little loss of accuracy for Ntable = ???. + +Only a single pair_coeff command is used with either style which +specifies a Vashishta potential file with parameters for all needed +elements. These are mapped to LAMMPS atom types by specifying N +additional arguments after the filename in the pair_coeff command, where N is the number of LAMMPS atom types: filename @@ -95,59 +114,52 @@ r0 (distance units) C costheta0 :ul -The non-annotated parameters are unitless. -The Vashishta potential file must contain entries for all the -elements listed in the pair_coeff command. It can also contain -entries for additional elements not being used in a particular -simulation; LAMMPS ignores those entries. -For a single-element simulation, only a single entry is required -(e.g. SiSiSi). For a two-element simulation, the file must contain 8 -entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that -specify parameters for all permutations of the two elements -interacting in three-body configurations. Thus for 3 elements, 27 -entries would be required, etc. +The non-annotated parameters are unitless. The Vashishta potential +file must contain entries for all the elements listed in the +pair_coeff command. It can also contain entries for additional +elements not being used in a particular simulation; LAMMPS ignores +those entries. For a single-element simulation, only a single entry +is required (e.g. SiSiSi). For a two-element simulation, the file +must contain 8 entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, +CCSi, CCC), that specify parameters for all permutations of the two +elements interacting in three-body configurations. Thus for 3 +elements, 27 entries would be required, etc. -Depending on the particular version of the Vashishta potential, -the values of these parameters may be keyed to the identities of -zero, one, two, or three elements. -In order to make the input file format unambiguous, general, -and simple to code, -LAMMPS uses a slightly confusing method for specifying parameters. -All parameters are divided into two classes: two-body and three-body. -Two-body and three-body parameters are handled differently, -as described below. -The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma, and r0. -They appear in the above formulae with two subscripts. -The parameters Zi and Zj are also classified as two-body parameters, -even though they only have 1 subscript. -The three-body parameters are B, C, costheta0. -They appear in the above formulae with three subscripts. -Two-body and three-body parameters are handled differently, -as described below. +Depending on the particular version of the Vashishta potential, the +values of these parameters may be keyed to the identities of zero, +one, two, or three elements. In order to make the input file format +unambiguous, general, and simple to code, LAMMPS uses a slightly +confusing method for specifying parameters. All parameters are +divided into two classes: two-body and three-body. Two-body and +three-body parameters are handled differently, as described below. +The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma, +and r0. They appear in the above formulae with two subscripts. The +parameters Zi and Zj are also classified as two-body parameters, even +though they only have 1 subscript. The three-body parameters are B, +C, costheta0. They appear in the above formulae with three +subscripts. Two-body and three-body parameters are handled +differently, as described below. -The first element in each entry is the center atom -in a three-body interaction, while the second and third elements -are two neighbor atoms. Three-body parameters for a central atom I -and two neighbors J and K are taken from the IJK entry. -Note that even though three-body parameters do not depend on the order of -J and K, LAMMPS stores three-body parameters for both IJK and IKJ. -The user must ensure that these values are equal. -Two-body parameters for an atom I interacting with atom J are taken from -the IJJ entry, where the 2nd and 3rd -elements are the same. Thus the two-body parameters -for Si interacting with C come from the SiCC entry. Note that even -though two-body parameters (except possibly gamma and r0 in U3) -do not depend on the order of the two elements, -LAMMPS will get the Si-C value from the SiCC entry -and the C-Si value from the CSiSi entry. The user must ensure -that these values are equal. Two-body parameters appearing -in entries where the 2nd and 3rd elements are different are -stored but never used. It is good practice to enter zero for -these values. Note that the three-body function U3 above -contains the two-body parameters gamma and r0. So U3 for a -central C atom bonded to an Si atom and a second C atom -will take three-body parameters from the CSiC entry, but -two-body parameters from the CCC and CSiSi entries. +The first element in each entry is the center atom in a three-body +interaction, while the second and third elements are two neighbor +atoms. Three-body parameters for a central atom I and two neighbors J +and K are taken from the IJK entry. Note that even though three-body +parameters do not depend on the order of J and K, LAMMPS stores +three-body parameters for both IJK and IKJ. The user must ensure that +these values are equal. Two-body parameters for an atom I interacting +with atom J are taken from the IJJ entry, where the 2nd and 3rd +elements are the same. Thus the two-body parameters for Si interacting +with C come from the SiCC entry. Note that even though two-body +parameters (except possibly gamma and r0 in U3) do not depend on the +order of the two elements, LAMMPS will get the Si-C value from the +SiCC entry and the C-Si value from the CSiSi entry. The user must +ensure that these values are equal. Two-body parameters appearing in +entries where the 2nd and 3rd elements are different are stored but +never used. It is good practice to enter zero for these values. Note +that the three-body function U3 above contains the two-body parameters +gamma and r0. So U3 for a central C atom bonded to an Si atom and a +second C atom will take three-body parameters from the CSiC entry, but +two-body parameters from the CCC and CSiSi entries. :line @@ -181,13 +193,7 @@ two different element types, mixing is performed by LAMMPS as described above from values in the potential file. This pair style does not support the "pair_modify"_pair_modify.html -{shift} and {tail} options, but it supports the {table} option. -Turning on the {table} option will use a linear interpolation table -for force and energy of the two-body term with {2**N} entries. -Tabulation can reduce the computational cost by a factor of 2-3 -with 20-bit to 12-bit table sizes, respectively, at a small to -moderate reduction in precision. It is not recommended to use -smaller than 12-bit tables. +shift, table, and tail options. This pair style does not write its information to "binary restart files"_restart.html, since it is stored in potential files. Thus, you @@ -209,11 +215,11 @@ the "Making LAMMPS"_Section_start.html#start_3 section for more info. This pair style requires the "newton"_newton.html setting to be "on" for pair interactions. -The Vashishta potential files provided with LAMMPS (see the -potentials directory) are parameterized for metal "units"_units.html. -You can use the Vashishta potential with any LAMMPS units, but you would need -to create your own Vashishta potential file with coefficients listed in the -appropriate units if your simulation doesn't use "metal" units. +The Vashishta potential files provided with LAMMPS (see the potentials +directory) are parameterized for metal "units"_units.html. You can +use the Vashishta potential with any LAMMPS units, but you would need +to create your own Vashishta potential file with coefficients listed +in the appropriate units if your simulation doesn't use "metal" units. [Related commands:] @@ -224,11 +230,14 @@ appropriate units if your simulation doesn't use "metal" units. :line :link(Vashishta1990) -[(Vashishta1990)] P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B 41, 12197 (1990). +[(Vashishta1990)] P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B +41, 12197 (1990). :link(Vashishta2007) -[(Vashishta2007)] P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino. J. Appl. Phys. 101, 103515 (2007). +[(Vashishta2007)] P. Vashishta, R. K. Kalia, A. Nakano, +J. P. Rino. J. Appl. Phys. 101, 103515 (2007). :link(Branicio2009) -[(Branicio2009)] Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002 +[(Branicio2009)] Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed +Matter 21 (2009) 095002 diff --git a/examples/vashishta/in.indiumphosphide-tabulated b/examples/vashishta/in.indiumphosphide.table similarity index 96% rename from examples/vashishta/in.indiumphosphide-tabulated rename to examples/vashishta/in.indiumphosphide.table index 3224b9fb75..954e366477 100644 --- a/examples/vashishta/in.indiumphosphide-tabulated +++ b/examples/vashishta/in.indiumphosphide.table @@ -42,9 +42,8 @@ mass 2 30.98 # choose potential -pair_style vashishta +pair_style vashishta/table 100000 0.2 pair_coeff * * InP.vashishta In P -pair_modify table 16 # setup neighbor style diff --git a/examples/vashishta/in.sio2-tabulated b/examples/vashishta/in.sio2.table similarity index 91% rename from examples/vashishta/in.sio2-tabulated rename to examples/vashishta/in.sio2.table index dc1381af2b..5b9016d96a 100644 --- a/examples/vashishta/in.sio2-tabulated +++ b/examples/vashishta/in.sio2.table @@ -11,9 +11,8 @@ replicate 4 4 4 velocity all create 2000.0 277387 mom yes displace_atoms all move 0.05 0.9 0.4 units box -pair_style vashishta +pair_style vashishta/table 100000 0.2 pair_coeff * * SiO.1990.vashishta Si O -pair_modify table 16 neighbor 0.3 bin neigh_modify delay 10 diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 46f057641f..19f6017907 100755 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -14,7 +14,6 @@ /* ---------------------------------------------------------------------- Contributing author: Yongnan Xiong (HNU), xyn@hnu.edu.cn Aidan Thompson (SNL) - Anders Hafreager (UiO), andershaf@gmail.com ------------------------------------------------------------------------- */ #include @@ -32,6 +31,7 @@ #include "neigh_list.h" #include "memory.h" #include "error.h" + using namespace LAMMPS_NS; #define MAXLINE 1024 @@ -52,15 +52,6 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) params = NULL; elem2param = NULL; map = NULL; - - neigh3BodyMax = 0; - neigh3BodyCount = NULL; - neigh3Body = NULL; - - useTable = false; - nTablebits = 12; - forceTable = NULL; - potentialTable = NULL; } /* ---------------------------------------------------------------------- @@ -75,11 +66,6 @@ PairVashishta::~PairVashishta() memory->destroy(params); memory->destroy(elem2param); - memory->destroy(forceTable); - memory->destroy(potentialTable); - memory->destroy(neigh3BodyCount); - memory->destroy(neigh3Body); - if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); @@ -87,86 +73,10 @@ PairVashishta::~PairVashishta() } } -void PairVashishta::modify_params(int narg, char **arg) -{ - if (narg == 0) error->all(FLERR,"Illegal pair_modify command"); - - int iarg = 0; - while (iarg < narg) { - if (strcmp(arg[iarg],"table") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); - - nTablebits = force->inumeric(FLERR,arg[iarg+1]); - if ((nTablebits < 0) || (nTablebits > sizeof(int)*CHAR_BIT-1)) - error->all(FLERR,"Illegal interpolation table size"); - - if (nTablebits == 0) { - useTable = false; - } else { - useTable = true; - } - - iarg += 2; - } else if (strcmp(arg[iarg],"tabinner") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); - - tabinner = force->numeric(FLERR,arg[iarg+1]); - if (tabinner <= 0.0) - error->all(FLERR,"Illegal inner cutoff for tabulation"); - iarg += 2; - } - } - updateTables(); -} - -void PairVashishta::validateNeigh3Body() { - if (atom->nlocal > neigh3BodyMax) { - neigh3BodyMax = atom->nmax; - memory->destroy(neigh3BodyCount); - memory->destroy(neigh3Body); - memory->create(neigh3BodyCount,neigh3BodyMax, "pair:vashishta:neigh3BodyCount"); - memory->create(neigh3Body,neigh3BodyMax, 1000, "pair:vashishta:neigh3Body"); // TODO: pick this number more wisely? - } -} - -void PairVashishta::updateTables() -{ - memory->destroy(forceTable); - memory->destroy(potentialTable); - forceTable = NULL; - potentialTable = NULL; - - if (!useTable) return; - - int ntable = 1<create(forceTable,nelements,nelements,ntable+1,"pair:vashishta:forceTable"); - memory->create(potentialTable,nelements,nelements,ntable+1,"pair:vashishta:potentialTable"); - - for (int i = 0; i < nelements; i++) { - for (int j = 0; j < nelements; j++) { - int ijparam = elem2param[i][j][j]; - for(int idx=0; idx <= ntable; idx++) { - double rsq = tabinnersq + idx*deltaR2; - double fpair, eng; - // call analytical version of two-body term function - twobody(params[ijparam], rsq, fpair, 1, eng, false); - forceTable[i][j][idx] = fpair; - potentialTable[i][j][idx] = eng; - } - } - } -} - /* ---------------------------------------------------------------------- */ void PairVashishta::compute(int eflag, int vflag) { - validateNeigh3Body(); int i,j,k,ii,jj,kk,inum,jnum,jnumm1; int itype,jtype,ktype,ijparam,ikparam,ijkparam; tagint itag,jtag; @@ -201,8 +111,6 @@ void PairVashishta::compute(int eflag, int vflag) ytmp = x[i][1]; ztmp = x[i][2]; - neigh3BodyCount[i] = 0; // Reset the 3-body neighbor list - // two-body interactions, skip half of them jlist = firstneigh[i]; @@ -213,21 +121,6 @@ void PairVashishta::compute(int eflag, int vflag) j &= NEIGHMASK; jtag = tag[j]; - jtype = map[type[j]]; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - ijparam = elem2param[itype][jtype][jtype]; - - if (rsq <= params[ijparam].cutsq2) { - neigh3Body[i][neigh3BodyCount[i]] = j; - neigh3BodyCount[i]++; - } - - if (rsq > params[ijparam].cutsq) continue; - if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { @@ -238,7 +131,17 @@ void PairVashishta::compute(int eflag, int vflag) if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } - twobody(params[ijparam],rsq,fpair,eflag,evdwl, useTable); + jtype = map[type[j]]; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + ijparam = elem2param[itype][jtype][jtype]; + if (rsq > params[ijparam].cutsq) continue; + + twobody(¶ms[ijparam],rsq,fpair,eflag,evdwl); f[i][0] += delx*fpair; f[i][1] += dely*fpair; @@ -251,8 +154,6 @@ void PairVashishta::compute(int eflag, int vflag) evdwl,0.0,fpair,delx,dely,delz); } - jlist = neigh3Body[i]; - jnum = neigh3BodyCount[i]; jnumm1 = jnum - 1; for (jj = 0; jj < jnumm1; jj++) { @@ -635,55 +536,32 @@ void PairVashishta::setup_params() if (params[m].cut > cutmax) cutmax = params[m].cut; if (params[m].r0 > cutmax) cutmax = params[m].r0; } - - updateTables(); } /* ---------------------------------------------------------------------- */ -void PairVashishta::twobody(const Param ¶m, double rsq, double &fforce, - int eflag, double &eng, bool tabulated) +void PairVashishta::twobody(Param *param, double rsq, double &fforce, + int eflag, double &eng) { - if (tabulated) { - if (rsq < tabinnersq) { - sprintf(estr,"Pair distance < table inner cutoff: " - "ijtype %d %d dist %g",param.ielement+1,param.jelement+1,sqrt(rsq)); - error->one(FLERR,estr); - } + double r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3; - const int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2; - // double - int will only keep the 0.xxxx part - const double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; + r = sqrt(rsq); + rinvsq = 1.0/rsq; + r4inv = rinvsq*rinvsq; + r6inv = rinvsq*r4inv; + reta = pow(r,-param->eta); + lam1r = r*param->lam1inv; + lam4r = r*param->lam4inv; + vc2 = param->zizj * exp(-lam1r)/r; + vc3 = param->mbigd * r4inv*exp(-lam4r); - double force0 = forceTable[param.ielement][param.jelement][tableIndex]; - double force1 = forceTable[param.ielement][param.jelement][tableIndex+1]; - fforce = (1.0 - fraction)*force0 + fraction*force1; // force is linearly interpolated between the two values - if (evflag) { - double energy0 = potentialTable[param.ielement][param.jelement][tableIndex]; - double energy1 = potentialTable[param.ielement][param.jelement][tableIndex+1]; - eng = (1.0 - fraction)*energy0 + fraction*energy1; - } - } else { - double r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3; - - r = sqrt(rsq); - rinvsq = 1.0/rsq; - r4inv = rinvsq*rinvsq; - r6inv = rinvsq*r4inv; - reta = pow(r,-param.eta); - lam1r = r*param.lam1inv; - lam4r = r*param.lam4inv; - vc2 = param.zizj * exp(-lam1r)/r; - vc3 = param.mbigd * r4inv*exp(-lam4r); - - fforce = (param.dvrc*r - - (4.0*vc3 + lam4r*vc3+param.big6w*r6inv - - param.heta*reta - vc2 - lam1r*vc2) - ) * rinvsq; - if (eflag) eng = param.bigh*reta - + vc2 - vc3 - param.bigw*r6inv - - r*param.dvrc + param.c0; - } + fforce = (param->dvrc*r + - (4.0*vc3 + lam4r*vc3+param->big6w*r6inv + - param->heta*reta - vc2 - lam1r*vc2) + ) * rinvsq; + if (eflag) eng = param->bigh*reta + + vc2 - vc3 - param->bigw*r6inv + - r*param->dvrc + param->c0; } /* ---------------------------------------------------------------------- */ @@ -742,15 +620,3 @@ void PairVashishta::threebody(Param *paramij, Param *paramik, Param *paramijk, if (eflag) eng = facrad; } - -/* ---------------------------------------------------------------------- - * add memory usage for tabulation - * ---------------------------------------------------------------------- */ - -double PairVashishta::memory_usage() -{ - double bytes = Pair::memory_usage(); - if (useTable) - bytes += 2*nelements*nelements*sizeof(double)*(1< +#include +#include +#include +#include "pair_vashishta_table.h" +#include "atom.h" +#include "neighbor.h" +#include "neigh_request.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairVashishtaTable::PairVashishtaTable(LAMMPS *lmp) : PairVashishta(lmp) +{ + neigh3BodyMax = 0; + neigh3BodyCount = NULL; + neigh3Body = NULL; + forceTable = NULL; + potentialTable = NULL; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairVashishtaTable::~PairVashishtaTable() +{ + memory->destroy(forceTable); + memory->destroy(potentialTable); + memory->destroy(neigh3BodyCount); + memory->destroy(neigh3Body); +} + +/* ---------------------------------------------------------------------- */ + +void PairVashishtaTable::compute(int eflag, int vflag) +{ + int i,j,k,ii,jj,kk,inum,jnum,jnumm1; + int itype,jtype,ktype,ijparam,ikparam,ijkparam; + tagint itag,jtag; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,rsq1,rsq2; + double delr1[3],delr2[3],fj[3],fk[3]; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // reallocate 3-body neighbor list if necessary + // NOTE: using 1000 is inefficient + // could make this a LAMMPS paged neighbor list + + if (nlocal > neigh3BodyMax) { + neigh3BodyMax = atom->nmax; + memory->destroy(neigh3BodyCount); + memory->destroy(neigh3Body); + memory->create(neigh3BodyCount,neigh3BodyMax, + "pair:vashishta:neigh3BodyCount"); + memory->create(neigh3Body,neigh3BodyMax,1000, + "pair:vashishta:neigh3Body"); + } + + // loop over full neighbor list of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itag = tag[i]; + itype = map[type[i]]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + // reset the 3-body neighbor list + + neigh3BodyCount[i] = 0; + + // two-body interactions, skip half of them + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtag = tag[j]; + + jtype = map[type[j]]; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + ijparam = elem2param[itype][jtype][jtype]; + + if (rsq <= params[ijparam].cutsq2) { + neigh3Body[i][neigh3BodyCount[i]] = j; + neigh3BodyCount[i]++; + } + + if (rsq > params[ijparam].cutsq) continue; + + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp && x[j][1] < ytmp) continue; + if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + + twobody_table(params[ijparam],rsq,fpair,eflag,evdwl); + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,fpair,delx,dely,delz); + } + + jlist = neigh3Body[i]; + jnum = neigh3BodyCount[i]; + jnumm1 = jnum - 1; + + for (jj = 0; jj < jnumm1; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = map[type[j]]; + ijparam = elem2param[itype][jtype][jtype]; + delr1[0] = x[j][0] - xtmp; + delr1[1] = x[j][1] - ytmp; + delr1[2] = x[j][2] - ztmp; + rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; + if (rsq1 >= params[ijparam].cutsq2) continue; + + for (kk = jj+1; kk < jnum; kk++) { + k = jlist[kk]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem2param[itype][ktype][ktype]; + ijkparam = elem2param[itype][jtype][ktype]; + + delr2[0] = x[k][0] - xtmp; + delr2[1] = x[k][1] - ytmp; + delr2[2] = x[k][2] - ztmp; + rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; + if (rsq2 >= params[ikparam].cutsq2) continue; + + threebody(¶ms[ijparam],¶ms[ikparam],¶ms[ijkparam], + rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); + + f[i][0] -= fj[0] + fk[0]; + f[i][1] -= fj[1] + fk[1]; + f[i][2] -= fj[2] + fk[2]; + f[j][0] += fj[0]; + f[j][1] += fj[1]; + f[j][2] += fj[2]; + f[k][0] += fk[0]; + f[k][1] += fk[1]; + f[k][2] += fk[2]; + + if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairVashishtaTable::twobody_table(const Param ¶m, double rsq, + double &fforce, int eflag, double &eng) +{ + // use analytic form if rsq is inside inner cutoff + + if (rsq < tabinnersq) { + Param *pparam = const_cast (¶m); + PairVashishta::twobody(pparam,rsq,fforce,eflag,eng); + return; + } + + // double -> int will only keep the 0.xxxx part + + const int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2; + const double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; + + // force/energy are linearly interpolated between two adjacent values + + double force0 = forceTable[param.ielement][param.jelement][tableIndex]; + double force1 = forceTable[param.ielement][param.jelement][tableIndex+1]; + fforce = (1.0 - fraction)*force0 + fraction*force1; + + if (evflag) { + double energy0 = potentialTable[param.ielement][param.jelement][tableIndex]; + double energy1 = potentialTable[param.ielement][param.jelement] + [tableIndex+1]; + eng = (1.0 - fraction)*energy0 + fraction*energy1; + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairVashishtaTable::settings(int narg, char **arg) +{ + if (narg != 2) error->all(FLERR,"Illegal pair_style command"); + + ntable = force->inumeric(FLERR,arg[0]); + tabinner = force->numeric(FLERR,arg[1]); + + if (tabinner <= 0.0) + error->all(FLERR,"Illegal inner cutoff for tabulation"); +} + +/* ---------------------------------------------------------------------- */ + +void PairVashishtaTable::setup_params() +{ + PairVashishta::setup_params(); + + create_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairVashishtaTable::create_tables() +{ + memory->destroy(forceTable); + memory->destroy(potentialTable); + forceTable = NULL; + potentialTable = NULL; + + tabinnersq = tabinner*tabinner; + + deltaR2 = (cutmax*cutmax - tabinnersq) / (ntable-1); + oneOverDeltaR2 = 1.0/deltaR2; + + memory->create(forceTable,nelements,nelements,ntable+1, + "pair:vashishta:forceTable"); + memory->create(potentialTable,nelements,nelements,ntable+1, + "pair:vashishta:potentialTable"); + + // tabulalate energy/force via analytic twobody() in parent + + int i,j,idx; + double rsq,fpair,eng; + + for (i = 0; i < nelements; i++) { + for (j = 0; j < nelements; j++) { + int ijparam = elem2param[i][j][j]; + for (idx = 0; idx <= ntable; idx++) { + rsq = tabinnersq + idx*deltaR2; + PairVashishta::twobody(¶ms[ijparam],rsq,fpair,1,eng); + forceTable[i][j][idx] = fpair; + potentialTable[i][j][idx] = eng; + } + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of tabulation arrays +------------------------------------------------------------------------- */ + +double PairVashishtaTable::memory_usage() +{ + double bytes = 2*nelements*nelements*sizeof(double)*ntable; + return bytes; +} diff --git a/src/MANYBODY/pair_vashishta_table.h b/src/MANYBODY/pair_vashishta_table.h new file mode 100755 index 0000000000..0052159ae1 --- /dev/null +++ b/src/MANYBODY/pair_vashishta_table.h @@ -0,0 +1,105 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(vashishta/table,PairVashishtaTable) + +#else + +#ifndef LMP_PAIR_VASHISHITA_TABLE_H +#define LMP_PAIR_VASHISHITA_TABLE_H + +#include "pair_vashishta.h" + +namespace LAMMPS_NS { + +class PairVashishtaTable : public PairVashishta { + public: + PairVashishtaTable(class LAMMPS *); + ~PairVashishtaTable(); + void compute(int, int); + void settings(int, char **); + double memory_usage(); + + private: + int ntable; + double deltaR2; + double oneOverDeltaR2; + double ***forceTable; // table of forces per element pair + double ***potentialTable; // table of potential energies + + int neigh3BodyMax; // max size of short neighborlist + int *neigh3BodyCount; // # of neighbors in short range + // 3 particle forces neighbor list + int **neigh3Body; // neighlist for short range 3 particle forces + + void twobody_table(const Param &, double, double &, int, double &); + void setup_params(); + void create_tables(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: Pair style Vashishta requires atom IDs + +This is a requirement to use the Vashishta potential. + +E: Pair style Vashishta requires newton pair on + +See the newton command. This is a restriction to use the Vashishta +potential. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +E: Cannot open Vashishta potential file %s + +The specified Vashishta potential file cannot be opened. Check that the path +and name are correct. + +E: Incorrect format in Vashishta potential file + +Incorrect number of words per line in the potential file. + +E: Illegal Vashishta parameter + +One or more of the coefficients defined in the potential file is +invalid. + +E: Potential file has duplicate entry + +The potential file has more than one entry for the same element. + +E: Potential file is missing an entry + +The potential file does not have a needed entry. + +*/