diff --git a/src/USER-OMP/pair_sw_omp.cpp b/src/USER-OMP/pair_sw_omp.cpp
index ba7eeacae1df8be1c050834610e83e4c05a10b5a..d6bd76a8972e87e5b26d0a4ef296046f7e2f0dbd 100644
--- a/src/USER-OMP/pair_sw_omp.cpp
+++ b/src/USER-OMP/pair_sw_omp.cpp
@@ -17,6 +17,7 @@
#include "atom.h"
#include "comm.h"
#include "force.h"
+#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
@@ -71,13 +72,13 @@ void PairSWOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG>
void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
{
- int i,j,k,ii,jj,kk,jnum,jnumm1;
+ int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
tagint itag,jtag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
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;
+ int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
evdwl = 0.0;
@@ -86,10 +87,13 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
+ const double cutshortsq = cutmax*cutmax;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
+ maxshort_thr = maxshort;
+ memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
double fxtmp,fytmp,fztmp;
@@ -109,12 +113,26 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
jlist = firstneigh[i];
jnum = numneigh[i];
+ int numshort = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
- jtag = tag[j];
+ delx = xtmp - x[j].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutshortsq) {
+ neighshort_thr[numshort++] = j;
+ if (numshort > maxshort_thr) {
+ maxshort_thr += maxshort_thr/2;
+ memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+ }
+ }
+
+ jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
@@ -126,12 +144,6 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
}
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;
@@ -148,11 +160,10 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
evdwl,0.0,fpair,delx,dely,delz,thr);
}
- jnumm1 = jnum - 1;
+ jnumm1 = numshort - 1;
for (jj = 0; jj < jnumm1; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
+ j = neighshort_thr[jj];
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
delr1[0] = x[j].x - xtmp;
@@ -164,9 +175,8 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0;
- for (kk = jj+1; kk < jnum; kk++) {
- k = jlist[kk];
- k &= NEIGHMASK;
+ for (kk = jj+1; kk < numshort; kk++) {
+ k = neighshort_thr[kk];
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
@@ -200,6 +210,7 @@ void PairSWOMP::eval(int iifrom, int iito, ThrData * const thr)
f[i].y += fytmp;
f[i].z += fztmp;
}
+ memory->destroy(neighshort_thr);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_tersoff_omp.cpp b/src/USER-OMP/pair_tersoff_omp.cpp
index 71d00108dfdec0d0d0750f27538c68f90bf728fa..f42f07830bce10e3140185d5fe1ca37249684767 100644
--- a/src/USER-OMP/pair_tersoff_omp.cpp
+++ b/src/USER-OMP/pair_tersoff_omp.cpp
@@ -17,6 +17,7 @@
#include "atom.h"
#include "comm.h"
#include "force.h"
+#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
@@ -73,14 +74,14 @@ void PairTersoffOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
{
- int i,j,k,ii,jj,kk,jnum;
+ int i,j,k,ii,jj,kk,jnum,maxshort_thr;
tagint itag,jtag;
int itype,jtype,ktype,iparam_ij,iparam_ijk;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fi[3],fj[3],fk[3];
double zeta_ij,prefactor;
- int *ilist,*jlist,*numneigh,**firstneigh;
+ int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
evdwl = 0.0;
@@ -89,10 +90,13 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
+ const double cutshortsq = cutmax*cutmax;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
+ maxshort_thr = maxshort;
+ memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
double fxtmp,fytmp,fztmp;
@@ -112,12 +116,26 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
jlist = firstneigh[i];
jnum = numneigh[i];
+ int numshort = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
- jtag = tag[j];
+ delx = xtmp - x[j].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutshortsq) {
+ neighshort_thr[numshort++] = j;
+ if (numshort > maxshort_thr) {
+ maxshort_thr += maxshort_thr/2;
+ memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+ }
+ }
+
+ jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
@@ -129,14 +147,8 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
}
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;
-
iparam_ij = elem2param[itype][jtype][jtype];
- if (rsq > params[iparam_ij].cutsq) continue;
+ if (rsq >= params[iparam_ij].cutsq) continue;
repulsive(¶ms[iparam_ij],rsq,fpair,EFLAG,evdwl);
@@ -155,9 +167,8 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
// skip immediately if I-J is not within cutoff
double fjxtmp,fjytmp,fjztmp;
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
+ for (jj = 0; jj < numshort; jj++) {
+ j = neighshort_thr[jj];
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype][jtype];
@@ -165,17 +176,16 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
delr1[1] = x[j].y - ytmp;
delr1[2] = x[j].z - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
- if (rsq1 > params[iparam_ij].cutsq) continue;
+ if (rsq1 >= params[iparam_ij].cutsq) continue;
// accumulate bondorder zeta for each i-j interaction via loop over k
fjxtmp = fjytmp = fjztmp = 0.0;
zeta_ij = 0.0;
- for (kk = 0; kk < jnum; kk++) {
+ for (kk = 0; kk < numshort; kk++) {
if (jj == kk) continue;
- k = jlist[kk];
- k &= NEIGHMASK;
+ k = neighshort_thr[kk];
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
@@ -183,7 +193,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
- if (rsq2 > params[iparam_ijk].cutsq) continue;
+ if (rsq2 >= params[iparam_ijk].cutsq) continue;
zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
}
@@ -204,10 +214,9 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
// attractive term via loop over k
- for (kk = 0; kk < jnum; kk++) {
+ for (kk = 0; kk < numshort; kk++) {
if (jj == kk) continue;
- k = jlist[kk];
- k &= NEIGHMASK;
+ k = neighshort_thr[kk];
ktype = map[type[k]];
iparam_ijk = elem2param[itype][jtype][ktype];
@@ -215,7 +224,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
delr2[1] = x[k].y - ytmp;
delr2[2] = x[k].z - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
- if (rsq2 > params[iparam_ijk].cutsq) continue;
+ if (rsq2 >= params[iparam_ijk].cutsq) continue;
attractive(¶ms[iparam_ijk],prefactor,
rsq1,rsq2,delr1,delr2,fi,fj,fk);
@@ -240,6 +249,7 @@ void PairTersoffOMP::eval(int iifrom, int iito, ThrData * const thr)
f[i].y += fytmp;
f[i].z += fztmp;
}
+ memory->destroy(neighshort_thr);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_vashishta_omp.cpp b/src/USER-OMP/pair_vashishta_omp.cpp
index 679eceab503669113a532401e355702968ad7bb3..4d05271bb29fab76b9a77d4f72de0e5d8041153f 100644
--- a/src/USER-OMP/pair_vashishta_omp.cpp
+++ b/src/USER-OMP/pair_vashishta_omp.cpp
@@ -17,6 +17,7 @@
#include "atom.h"
#include "comm.h"
#include "force.h"
+#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
@@ -71,13 +72,13 @@ void PairVashishtaOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG>
void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
{
- int i,j,k,ii,jj,kk,jnum,jnumm1;
+ int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
tagint itag,jtag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
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;
+ int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
evdwl = 0.0;
@@ -86,10 +87,13 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
+ const double cutshortsq = r0max*r0max;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
+ maxshort_thr = maxshort;
+ memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
double fxtmp,fytmp,fztmp;
@@ -109,12 +113,26 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
jlist = firstneigh[i];
jnum = numneigh[i];
+ int numshort = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
- jtag = tag[j];
+ delx = xtmp - x[j].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ if (rsq < cutshortsq) {
+ neighshort_thr[numshort++] = j;
+ if (numshort > maxshort_thr) {
+ maxshort_thr += maxshort_thr/2;
+ memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+ }
+ }
+
+ jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
@@ -126,12 +144,6 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
}
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;
@@ -148,11 +160,10 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
evdwl,0.0,fpair,delx,dely,delz,thr);
}
- jnumm1 = jnum - 1;
+ jnumm1 = numshort - 1;
for (jj = 0; jj < jnumm1; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
+ j = neighshort_thr[jj];
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
delr1[0] = x[j].x - xtmp;
@@ -164,9 +175,8 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0;
- for (kk = jj+1; kk < jnum; kk++) {
- k = jlist[kk];
- k &= NEIGHMASK;
+ for (kk = jj+1; kk < numshort; kk++) {
+ k = neighshort_thr[kk];
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
@@ -200,6 +210,7 @@ void PairVashishtaOMP::eval(int iifrom, int iito, ThrData * const thr)
f[i].y += fytmp;
f[i].z += fztmp;
}
+ memory->destroy(neighshort_thr);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_vashishta_table_omp.cpp b/src/USER-OMP/pair_vashishta_table_omp.cpp
index b5769438c4709789f7bebec4416dd6eebbde04b4..2ec076ad69f9e520cb1ca7b1fd1f1399328704df 100644
--- a/src/USER-OMP/pair_vashishta_table_omp.cpp
+++ b/src/USER-OMP/pair_vashishta_table_omp.cpp
@@ -41,20 +41,6 @@ void PairVashishtaTableOMP::compute(int eflag, int vflag)
ev_setup(eflag,vflag);
} else evflag = vflag_fdotr = 0;
- // reallocate 3-body neighbor list if necessary
- // NOTE: using 1000 is inefficient
- // could make this a LAMMPS paged neighbor list
-
- 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");
- }
-
const int nall = atom->nlocal + atom->nghost;
const int nthreads = comm->nthreads;
const int inum = list->inum;
@@ -86,13 +72,13 @@ void PairVashishtaTableOMP::compute(int eflag, int vflag)
template <int EVFLAG, int EFLAG>
void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
{
- int i,j,k,ii,jj,kk,jnum,jnumm1;
+ int i,j,k,ii,jj,kk,jnum,jnumm1,maxshort_thr;
tagint itag,jtag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
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;
+ int *ilist,*jlist,*numneigh,**firstneigh,*neighshort_thr;
evdwl = 0.0;
@@ -101,13 +87,16 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
const tagint * _noalias const tag = atom->tag;
const int * _noalias const type = atom->type;
const int nlocal = atom->nlocal;
+ const double cutshortsq = r0max*r0max;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
+ maxshort_thr = maxshort;
+ memory->create(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
double fxtmp,fytmp,fztmp;
-
+
// loop over full neighbor list of my atoms
for (ii = iifrom; ii < iito; ++ii) {
@@ -120,35 +109,30 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
ztmp = x[i].z;
fxtmp = fytmp = fztmp = 0.0;
- // reset the 3-body neighbor list
-
- neigh3BodyCount[i] = 0;
-
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
+ int numshort = 0;
for (jj = 0; jj < jnum; jj++) {
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 < cutshortsq) {
+ neighshort_thr[numshort++] = j;
+ if (numshort > maxshort_thr) {
+ maxshort_thr += maxshort_thr/2;
+ memory->grow(neighshort_thr,maxshort_thr,"pair_thr:neighshort_thr");
+ }
}
- if (rsq >= params[ijparam].cutsq) continue;
-
+ jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
@@ -159,6 +143,10 @@ void PairVashishtaTableOMP::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]];
+ ijparam = elem2param[itype][jtype][jtype];
+ if (rsq >= params[ijparam].cutsq) continue;
+
twobody_table(params[ijparam],rsq,fpair,EFLAG,evdwl);
fxtmp += delx*fpair;
@@ -172,13 +160,10 @@ void PairVashishtaTableOMP::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;
+ jnumm1 = numshort - 1;
for (jj = 0; jj < jnumm1; jj++) {
- j = jlist[jj];
- j &= NEIGHMASK;
+ j = neighshort_thr[jj];
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
delr1[0] = x[j].x - xtmp;
@@ -190,9 +175,8 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
double fjxtmp,fjytmp,fjztmp;
fjxtmp = fjytmp = fjztmp = 0.0;
- for (kk = jj+1; kk < jnum; kk++) {
- k = jlist[kk];
- k &= NEIGHMASK;
+ for (kk = jj+1; kk < numshort; kk++) {
+ k = neighshort_thr[kk];
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
@@ -226,6 +210,7 @@ void PairVashishtaTableOMP::eval(int iifrom, int iito, ThrData * const thr)
f[i].y += fytmp;
f[i].z += fztmp;
}
+ memory->destroy(neighshort_thr);
}
/* ---------------------------------------------------------------------- */