ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

same changes to other NPair and NStencil methods

Browse Source
This commit is contained in:
Steve Plimpton
2023-07-10 18:25:29 -07:00
parent 4ab95611db
commit fbbf44fb8e
15 changed files with 313 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/npair_half_bin_newton_tri.cpp
View File

@ -16,11 +16,11 @@
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
using namespace LAMMPS_NS;

119
src/npair_half_multi_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neighbor.h"
@ -39,11 +40,13 @@ NPairHalfMultiNewtonTri::NPairHalfMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalfMultiNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
int js;
const double delta = 0.01 * force->angstrom;
int *collection = neighbor->collection;
double **x = atom->x;
int *type = atom->type;
@ -72,6 +75,8 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -86,65 +91,79 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
// stencil is half if i same size as j
// stencil is full if i smaller than j
// if half: pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic:
// stencil is empty if i larger than j
// stencil is full if i smaller than j
// stencil is full if i same size as j
// for i smaller than j:
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
// if same size (same collection), use half stencil
if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
}
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = j;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = j;
else if (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else neighptr[n++] = j;
}
}
}
}
ilist[inum++] = i;

44
src/npair_half_multi_old_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,11 +39,13 @@ NPairHalfMultiOldNewtonTri::NPairHalfMultiOldNewtonTri(LAMMPS *lmp) : NPair(lmp)
void NPairHalfMultiOldNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*s;
double *cutsq,*distsq;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -71,6 +74,7 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -81,13 +85,12 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -98,14 +101,23 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

45
src/npair_half_respa_bin_newton_tri.cpp
View File

@ -16,10 +16,11 @@
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
using namespace LAMMPS_NS;
@ -38,10 +39,12 @@ NPairHalfRespaBinNewtonTri::NPairHalfRespaBinNewtonTri(LAMMPS *lmp) :
void NPairHalfRespaBinNewtonTri::build(NeighList *list)
{
int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -94,6 +97,7 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
neighptr_middle = ipage_middle->vget();
}
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -105,22 +109,33 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

24
src/npair_half_respa_nsq_newton.cpp
View File

@ -16,9 +16,10 @@
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "force.h"
#include "group.h"
#include "molecule.h"
#include "domain.h"
#include "my_page.h"
#include "error.h"
@ -38,12 +39,15 @@ NPairHalfRespaNsqNewton::NPairHalfRespaNsqNewton(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalfRespaNsqNewton::build(NeighList *list)
{
int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle,bitmask;
int i,j,n,itype,jtype,n_inner,n_middle,bitmask;
int imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *neighptr,*neighptr_inner,*neighptr_middle;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
@ -112,6 +116,12 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -122,6 +132,14 @@ void NPairHalfRespaNsqNewton::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

2
src/npair_half_respa_nsq_newton.h
View File

@ -15,7 +15,7 @@
// clang-format off
NPairStyle(half/respa/nsq/newton,
NPairHalfRespaNsqNewton,
NP_HALF | NP_RESPA | NP_NSQ | NP_NEWTON | NP_ORTHO);
NP_HALF | NP_RESPA | NP_NSQ | NP_NEWTON | NP_ORTHO | NP_TRI);
// clang-format on
#else

41
src/npair_half_size_bin_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -39,11 +40,13 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) :
void NPairHalfSizeBinNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -76,6 +79,7 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -87,22 +91,33 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list)
}
// loop over all atoms in bins in stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
if (exclude && exclusion(i,j,type[i],type[j],mask,molecule)) continue;

107
src/npair_half_size_multi_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neighbor.h"
@ -41,11 +42,13 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js;
int which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
const double delta = 0.01 * force->angstrom;
int *collection = neighbor->collection;
double **x = atom->x;
double *radius = atom->radius;
@ -78,6 +81,8 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
for (i = 0; i < nlocal; i++) {
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
icollection = collection[i];
xtmp = x[i][0];
@ -93,11 +98,13 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
ibin = atom2bin[i];
// loop through stencils for all collections
for (jcollection = 0; jcollection < ncollections; jcollection++) {
// if same collection use own bin
if (icollection == jcollection) jbin = ibin;
else jbin = coord2bin(x[i], jcollection);
else jbin = coord2bin(x[i], jcollection);
// loop over all atoms in bins in stencil
// stencil is empty if i larger than j
@ -108,56 +115,66 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
s = stencil_multi[icollection][jcollection];
ns = nstencil_multi[icollection][jcollection];
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
for (k = 0; k < ns; k++) {
js = binhead_multi[jcollection][jbin + s[k]];
for (j = js; j >= 0; j = bins[j]) {
// if same size (same collection), use half stencil
if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
// if same size (same collection), exclude half of interactions
if (cutcollectionsq[icollection][icollection] ==
cutcollectionsq[jcollection][jcollection]) {
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
jh = j;
if (history && rsq < radsum*radsum)
jh = jh ^ mask_history;
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = jh;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = jh;
else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
} else neighptr[n++] = jh;
}
}
}
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radi + radius[j];
cutdistsq = (radsum+skin) * (radsum+skin);
if (rsq <= cutdistsq) {
jh = j;
if (history && rsq < radsum*radsum)
jh = jh ^ mask_history;
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
else if (imol >= 0)
which = find_special(onemols[imol]->special[iatom],
onemols[imol]->nspecial[iatom],
tag[j]-tagprev);
else which = 0;
if (which == 0) neighptr[n++] = jh;
else if (domain->minimum_image_check(delx,dely,delz))
neighptr[n++] = jh;
else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
} else neighptr[n++] = jh;
}
}
}
}
ilist[inum++] = i;

45
src/npair_half_size_multi_old_newton_tri.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,12 +39,14 @@ NPairHalfSizeMultiOldNewtonTri::NPairHalfSizeMultiOldNewtonTri(LAMMPS *lmp) : NP
void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
{
int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate;
tagint tagprev;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutdistsq;
int *neighptr,*s;
double *cutsq,*distsq;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -76,6 +79,7 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -87,13 +91,12 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
tagprev = tag[i] - iatom - 1;
}
// loop over all atoms in bins, including self, in stencil
// skip if i,j neighbor cutoff is less than bin distance
// bins below self are excluded from stencil
// pairs for atoms j "below" i are excluded
// below = lower z or (equal z and lower y) or (equal zy and lower x)
// (equal zyx and j <= i)
// latter excludes self-self interaction but allows superposed atoms
// loop over all atoms in bins in stencil
// for triclinic, bin stencil is full in all 3 dims
// must use itag/jtag to eliminate half the I/J interactions
// cannot use I/J exact coord comparision
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
ibin = atom2bin[i];
s = stencil_multi_old[itype];
@ -104,14 +107,24 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list)
for (j = binhead[ibin+s[k]]; j >= 0; j = bins[j]) {
jtype = type[j];
if (cutsq[jtype] < distsq[k]) continue;
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {
if (x[j][1] < ytmp) continue;
if (x[j][1] == ytmp) {
if (x[j][0] < xtmp) continue;
if (x[j][0] == xtmp && j <= i) continue;
}
}
if (j <= i) continue;
if (j >= nlocal) {
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
}
}
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;

22
src/npair_half_size_nsq_newton.cpp
View File

@ -18,6 +18,7 @@
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "molecule.h"
#include "group.h"
#include "my_page.h"
@ -39,12 +40,15 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {}
void NPairHalfSizeNsqNewton::build(NeighList *list)
{
int i,j,jh,n,itag,jtag,bitmask,which,imol,iatom,moltemplate;
tagint tagprev;
int i,j,jh,n,bitmask,which,imol,iatom,moltemplate;
tagint itag,jtag,tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radi,radsum,cutsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
double *radius = atom->radius;
tagint *tag = atom->tag;
@ -93,6 +97,12 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
}
// loop over remaining atoms, owned and ghost
// use itag/jtap comparision to eliminate half the interactions
// itag = jtag is possible for long cutoffs that include images of self
// for triclinic, must use delta to eliminate half the I/J interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
for (j = i+1; j < nall; j++) {
if (includegroup && !(mask[j] & bitmask)) continue;
@ -103,6 +113,14 @@ void NPairHalfSizeNsqNewton::build(NeighList *list)
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

19
src/npair_halffull_newton_trim.cpp
View File

@ -14,7 +14,9 @@
#include "npair_halffull_newton_trim.h"
#include "atom.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "my_page.h"
#include "neigh_list.h"
@ -38,6 +40,9 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
double xtmp, ytmp, ztmp;
double delx, dely, delz, rsq;
const double delta = 0.01 * force->angstrom;
const int triclinic = domain->triclinic;
double **x = atom->x;
int nlocal = atom->nlocal;
@ -68,6 +73,11 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
ztmp = x[i][2];
// loop over full neighbor list
// use i < j < nlocal to eliminate half the local/local interactions
// for triclinic, must use delta to eliminate half the local/ghost interactions
// cannot use I/J exact coord comparision as for orthog
// b/c transforming orthog -> lambda -> orthog for ghost atoms
// with an added PBC offset can shift all 3 coords by epsilon
jlist = firstneigh_full[i];
jnum = numneigh_full[i];
@ -75,8 +85,17 @@ void NPairHalffullNewtonTrim::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (j < nlocal) {
if (i > j) continue;
} else if (triclinic) {
if (fabs(x[j][2]-ztmp) > delta) {
if (x[j][2] < ztmp) continue;
} else if (fabs(x[j][1]-ytmp) > delta) {
if (x[j][1] < ytmp) continue;
} else {
if (x[j][0] < xtmp) continue;
}
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp) {

2
src/nstencil_half_multi_2d_tri.cpp
View File

@ -80,7 +80,7 @@ void NStencilHalfMulti2dTri::create()
cutsq = cutcollectionsq[icollection][jcollection];
if (flag_half_multi[icollection][jcollection]) {
for (j = 0; j <= sy; j++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i, j, 0, bin_collection) < cutsq)
stencil_multi[icollection][jcollection][ns++] = j * mbinx + i;

2
src/nstencil_half_multi_3d_tri.cpp
View File

@ -81,7 +81,7 @@ void NStencilHalfMulti3dTri::create()
cutsq = cutcollectionsq[icollection][jcollection];
if (flag_half_multi[icollection][jcollection]) {
for (k = 0; k <= sz; k++)
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++)
if (bin_distance_multi(i, j, k, bin_collection) < cutsq)

2
src/nstencil_half_multi_old_2d_tri.cpp
View File

@ -37,7 +37,7 @@ void NStencilHalfMultiOld2dTri::create()
s = stencil_multi_old[itype];
distsq = distsq_multi_old[itype];
n = 0;
for (j = 0; j <= sy; j++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++) {
rsq = bin_distance(i, j, 0);
if (rsq < typesq) {

2
src/nstencil_half_multi_old_3d_tri.cpp
View File

@ -37,7 +37,7 @@ void NStencilHalfMultiOld3dTri::create()
s = stencil_multi_old[itype];
distsq = distsq_multi_old[itype];
n = 0;
for (k = 0; k <= sz; k++)
for (k = -sz; k <= sz; k++)
for (j = -sy; j <= sy; j++)
for (i = -sx; i <= sx; i++) {
rsq = bin_distance(i, j, k);