diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 41a7e9afe0..7f77eaf81c 100755 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -54,6 +54,10 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) elem2param = NULL; map = NULL; + neigh3BodyMax = 0; + neigh3BodyCount = NULL; + neigh3Body = NULL; + useTable = true; tableSize = 65536; forceTable = NULL; @@ -78,8 +82,7 @@ PairVashishta::~PairVashishta() delete [] map; } } -#include -using namespace std; + void PairVashishta::modify_params(int narg, char **arg) { if(narg < 2 || narg > 3) // We accept table yes/no [tableSize] @@ -102,6 +105,16 @@ void PairVashishta::modify_params(int narg, char **arg) } } +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::createTable() { int ntypes = atom->ntypes+1; @@ -140,6 +153,7 @@ void PairVashishta::createTable() 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; @@ -174,6 +188,8 @@ 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]; @@ -184,6 +200,21 @@ 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) { @@ -194,16 +225,6 @@ void PairVashishta::compute(int eflag, int vflag) if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } - 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, useTable); f[i][0] += delx*fpair; @@ -217,6 +238,8 @@ 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++) { diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index a8b971b668..47a36e39a4 100755 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -53,6 +53,10 @@ class PairVashishta : public Pair { double ***forceTable; // Table containing the forces double ***potentialTable; // Table containing the potential energies + int neigh3BodyMax; // Max size of short neighborlist + int *neigh3BodyCount; // Number of neighbors in short range 3 particle forces neighbor list + int **neigh3Body; // Neighborlist for short range 3 particle forces + double cutmax; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements @@ -66,6 +70,7 @@ class PairVashishta : public Pair { void read_file(char *); void setup_params(); void createTable(); + void validateNeigh3Body(); void twobody(Param *, double, double &, int, double &, bool); void threebody(Param *, Param *, Param *, double, double, double *, double *, double *, double *, int, double &); diff --git a/src/USER-OMP/pair_vashishta_omp.cpp b/src/USER-OMP/pair_vashishta_omp.cpp index 99653803e7..e225c9aacb 100644 --- a/src/USER-OMP/pair_vashishta_omp.cpp +++ b/src/USER-OMP/pair_vashishta_omp.cpp @@ -36,6 +36,8 @@ PairVashishtaOMP::PairVashishtaOMP(LAMMPS *lmp) : void PairVashishtaOMP::compute(int eflag, int vflag) { + validateNeigh3Body(); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; @@ -105,6 +107,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr) ztmp = x[i].z; fxtmp = fytmp = fztmp = 0.0; + neigh3BodyCount[i] = 0; + // two-body interactions, skip half of them jlist = firstneigh[i]; @@ -114,6 +118,21 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr) j = jlist[jj]; j &= NEIGHMASK; jtag = tag[j]; + jtype = map[type[j]]; + + delx = xtmp - x[j].x; + dely = ytmp - x[j].y; + delz = ztmp - x[j].z; + 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; @@ -125,16 +144,6 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr) if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue; } - jtype = map[type[j]]; - - delx = xtmp - x[j].x; - dely = ytmp - x[j].y; - delz = ztmp - x[j].z; - 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, useTable); fxtmp += delx*fpair; @@ -148,6 +157,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr) evdwl,0.0,fpair,delx,dely,delz,thr); } + jlist = neigh3Body[i]; + jnum = neigh3BodyCount[i]; jnumm1 = jnum - 1; for (jj = 0; jj < jnumm1; jj++) {