diff --git a/src/OPENMP/npair_half_bin_newton_tri_omp.cpp b/src/OPENMP/npair_half_bin_newton_tri_omp.cpp index e754456ef1..3ad07acd56 100644 --- a/src/OPENMP/npair_half_bin_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_bin_newton_tri_omp.cpp @@ -18,6 +18,7 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" +#include "force.h" #include "molecule.h" #include "domain.h" #include "my_page.h" @@ -40,6 +41,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int molecular = atom->molecular; const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; + const double delta = 0.01 * force->angstrom; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -48,7 +50,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,k,n,itype,jtype,ibin,which,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; @@ -79,6 +81,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -98,12 +101,22 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) 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; + } } } @@ -119,7 +132,7 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) if (molecular != Atom::ATOMIC) { if (!moltemplate) which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >=0) + else if (imol >= 0) which = find_special(onemols[imol]->special[iatom], onemols[imol]->nspecial[iatom], tag[j]-tagprev); diff --git a/src/OPENMP/npair_half_nsq_newton_omp.cpp b/src/OPENMP/npair_half_nsq_newton_omp.cpp index cb08cb7f7a..726814c6f0 100644 --- a/src/OPENMP/npair_half_nsq_newton_omp.cpp +++ b/src/OPENMP/npair_half_nsq_newton_omp.cpp @@ -15,14 +15,16 @@ #include "omp_compat.h" #include "npair_half_nsq_newton_omp.h" #include "npair_omp.h" -#include "neigh_list.h" + #include "atom.h" #include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" #include "group.h" #include "molecule.h" -#include "domain.h" #include "my_page.h" -#include "error.h" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -42,6 +44,8 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list) const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; const int molecular = atom->molecular; const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -49,8 +53,8 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,itag,jtag,which,imol,iatom; - tagint tagprev; + int i,j,n,itype,jtype,which,imol,iatom; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; @@ -106,6 +110,14 @@ void NPairHalfNsqNewtonOmp::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) { @@ -127,7 +139,7 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list) if (molecular != Atom::ATOMIC) { if (!moltemplate) which = find_special(special[i],nspecial[i],tag[j]); - else if (imol >=0) + else if (imol >= 0) which = find_special(onemols[imol]->special[iatom], onemols[imol]->nspecial[iatom], tag[j]-tagprev); diff --git a/src/OPENMP/npair_halffull_newton_omp.cpp b/src/OPENMP/npair_halffull_newton_omp.cpp index abd5f7eacb..e833ab3095 100644 --- a/src/OPENMP/npair_halffull_newton_omp.cpp +++ b/src/OPENMP/npair_halffull_newton_omp.cpp @@ -15,7 +15,9 @@ #include "npair_halffull_newton_omp.h" #include "atom.h" +#include "domain.h" #include "error.h" +#include "force.h" #include "my_page.h" #include "neigh_list.h" #include "npair_omp.h" @@ -38,6 +40,8 @@ NPairHalffullNewtonOmp::NPairHalffullNewtonOmp(LAMMPS *lmp) : NPair(lmp) {} void NPairHalffullNewtonOmp::build(NeighList *list) { const int inum_full = list->listfull->inum; + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -83,8 +87,17 @@ void NPairHalffullNewtonOmp::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) { diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 71a15df59e..f6cbf1b1af 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -40,9 +40,9 @@ void NPairHalfBinNewtonTri::build(NeighList *list) tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - - double delta = 0.01 * force->angstrom; - + + const double delta = 0.01 * force->angstrom; + double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -85,7 +85,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) // 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 coord comparision + // 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 diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index 0174a78900..4d5afbdd3e 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -13,15 +13,16 @@ ------------------------------------------------------------------------- */ #include "npair_half_nsq_newton.h" -#include "neigh_list.h" + #include "atom.h" #include "atom_vec.h" +#include "domain.h" +#include "error.h" #include "force.h" #include "group.h" #include "molecule.h" -#include "domain.h" #include "my_page.h" -#include "error.h" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -42,9 +43,9 @@ void NPairHalfNsqNewton::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - double delta = 0.01 * force->angstrom; - int triclinic = domain->triclinic; - + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; + double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -92,7 +93,7 @@ void NPairHalfNsqNewton::build(NeighList *list) // 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 direct I/J coord comparision as for orthog + // 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 @@ -105,15 +106,15 @@ void NPairHalfNsqNewton::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 { + } 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) { if (x[j][1] < ytmp) continue; diff --git a/src/npair_halffull_newton.cpp b/src/npair_halffull_newton.cpp index af15b27eac..12320c46f3 100644 --- a/src/npair_halffull_newton.cpp +++ b/src/npair_halffull_newton.cpp @@ -39,8 +39,8 @@ void NPairHalffullNewton::build(NeighList *list) int *neighptr, *jlist; double xtmp, ytmp, ztmp; - double delta = 0.01 * force->angstrom; - int triclinic = domain->triclinic; + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; double **x = atom->x; int nlocal = atom->nlocal; @@ -72,7 +72,7 @@ void NPairHalffullNewton::build(NeighList *list) // 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 direct I/J coord comparision as for orthog + // 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 @@ -82,17 +82,17 @@ void NPairHalffullNewton::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; - } + 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) { @@ -100,7 +100,7 @@ void NPairHalffullNewton::build(NeighList *list) if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } } - + neighptr[n++] = joriginal; }