From bacfcd205093c626a30f43e157704118d001f945 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 7 Jul 2023 07:36:05 -0700 Subject: [PATCH 01/41] change triclinic logic to not depend on exact I,J atom coords --- src/npair_half_bin_newton_tri.cpp | 39 ++++++++++++++++++++++++++++--- src/npair_half_nsq_newton.cpp | 14 ++++++++++- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 88ef993a41..227a25c321 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -16,6 +16,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" @@ -36,10 +37,12 @@ NPairHalfBinNewtonTri::NPairHalfBinNewtonTri(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfBinNewtonTri::build(NeighList *list) { int i,j,k,n,itype,jtype,ibin,which,imol,iatom,moltemplate; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - + + double angstrom = force->angstrom; + double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -68,6 +71,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) n = 0; neighptr = ipage->vget(); + itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -87,6 +91,34 @@ void NPairHalfBinNewtonTri::build(NeighList *list) ibin = atom2bin[i]; for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + + 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) > angstrom) { + if (x[j][2] < ztmp) continue; + } else if (fabs(x[j][1]-ytmp) > angstrom) { + if (x[j][1] < ytmp) continue; + } else { + if (x[j][0] < xtmp) continue; + } + + /* + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + */ + } + } + + /* if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -95,7 +127,8 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] == xtmp && j <= i) continue; } } - + */ + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index e5f3138f0a..20ffbf6977 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -16,6 +16,7 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" +#include "force.h" #include "group.h" #include "molecule.h" #include "domain.h" @@ -41,6 +42,9 @@ void NPairHalfNsqNewton::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; + double angstrom = force->angstrom; + int triclinic = domain->triclinic; + double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -96,7 +100,15 @@ void NPairHalfNsqNewton::build(NeighList *list) if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; - } else { + } else if (triclinic) { + if (fabs(x[j][2]-ztmp) > angstrom) { + if (x[j][2] < ztmp) continue; + } else if (fabs(x[j][1]-ytmp) > angstrom) { + 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; From 129264aa148bb2b10ce0d6dbfda866443f58d749 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 7 Jul 2023 08:42:46 -0700 Subject: [PATCH 02/41] debugging --- src/npair_half_bin_newton_tri.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 227a25c321..00b5722673 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -92,6 +92,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + /* if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { @@ -108,17 +109,15 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] < xtmp) continue; } - /* if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } - */ } } - - /* + */ + if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -127,7 +126,6 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] == xtmp && j <= i) continue; } } - */ jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; From 42f2a4b5b0e8021b01b1655b2c555cd392539c18 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 7 Jul 2023 08:58:14 -0700 Subject: [PATCH 03/41] exclude self interactions and double counting of own/own --- src/npair_half_bin_newton_tri.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 00b5722673..4cd12b2e7c 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -92,7 +92,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { - /* + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { @@ -109,15 +109,15 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] < xtmp) continue; } - if (x[j][2] < ztmp) continue; - if (x[j][2] == ztmp) { - if (x[j][1] < ytmp) continue; - if (x[j][1] == ytmp && x[j][0] < xtmp) continue; - } + //if (x[j][2] < ztmp) continue; + //if (x[j][2] == ztmp) { + // if (x[j][1] < ytmp) continue; + // if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + // } } } - */ - + + /* if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -126,6 +126,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] == xtmp && j <= i) continue; } } + */ jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; @@ -136,6 +137,8 @@ void NPairHalfBinNewtonTri::build(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { + //printf("NEIGH i,j %d %d ijtag %d %d dist %g\n", + // i,j,tag[i],tag[j],sqrt(rsq)); if (molecular != Atom::ATOMIC) { if (!moltemplate) which = find_special(special[i],nspecial[i],tag[j]); From abadf9412afb695eabeb996d507e6533ec686d07 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 8 Jul 2023 10:26:34 -0700 Subject: [PATCH 04/41] check old results --- src/npair_half_bin_newton_tri.cpp | 7 ++++--- src/npair_half_nsq_newton.cpp | 4 +++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 4cd12b2e7c..61d67fed2a 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -92,6 +92,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { + /* if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; @@ -116,8 +117,9 @@ void NPairHalfBinNewtonTri::build(NeighList *list) // } } } - - /* + */ + + if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -126,7 +128,6 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] == xtmp && j <= i) continue; } } - */ jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index 20ffbf6977..995c9cfbed 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -100,7 +100,8 @@ void NPairHalfNsqNewton::build(NeighList *list) if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; - } else if (triclinic) { + /* + } else if (triclinic) { if (fabs(x[j][2]-ztmp) > angstrom) { if (x[j][2] < ztmp) continue; } else if (fabs(x[j][1]-ytmp) > angstrom) { @@ -108,6 +109,7 @@ void NPairHalfNsqNewton::build(NeighList *list) } else { if (x[j][0] < xtmp) continue; } + */ } else { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { From ce1a084a0efba244310fd4d8448b726210030f8d Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Sat, 8 Jul 2023 18:31:55 -0700 Subject: [PATCH 05/41] expand stencil for triclinic neighbor build --- src/npair_half_bin_newton_tri.cpp | 15 ++++++--------- src/npair_half_nsq_newton.cpp | 10 ++++------ src/npair_halffull_newton.cpp | 15 +++++++++++++++ src/nstencil_half_bin_3d_tri.cpp | 3 ++- 4 files changed, 27 insertions(+), 16 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 61d67fed2a..17504ac3af 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -41,7 +41,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - double angstrom = force->angstrom; + double delta = 0.01 * force->angstrom; double **x = atom->x; int *type = atom->type; @@ -92,7 +92,6 @@ void NPairHalfBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { - /* if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; @@ -101,15 +100,13 @@ void NPairHalfBinNewtonTri::build(NeighList *list) } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; } else { - - if (fabs(x[j][2]-ztmp) > angstrom) { + if (fabs(x[j][2]-ztmp) > delta) { if (x[j][2] < ztmp) continue; - } else if (fabs(x[j][1]-ytmp) > angstrom) { + } else if (fabs(x[j][1]-ytmp) > delta) { if (x[j][1] < ytmp) continue; } else { if (x[j][0] < xtmp) continue; } - //if (x[j][2] < ztmp) continue; //if (x[j][2] == ztmp) { // if (x[j][1] < ytmp) continue; @@ -117,9 +114,8 @@ void NPairHalfBinNewtonTri::build(NeighList *list) // } } } - */ - - + + /* if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -128,6 +124,7 @@ void NPairHalfBinNewtonTri::build(NeighList *list) if (x[j][0] == xtmp && j <= i) continue; } } + */ jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index 995c9cfbed..be06393f58 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -42,7 +42,7 @@ void NPairHalfNsqNewton::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - double angstrom = force->angstrom; + double delta = 0.01 * force->angstrom; int triclinic = domain->triclinic; double **x = atom->x; @@ -100,16 +100,14 @@ 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) > angstrom) { + } else if (triclinic) { + if (fabs(x[j][2]-ztmp) > delta) { if (x[j][2] < ztmp) continue; - } else if (fabs(x[j][1]-ytmp) > angstrom) { + } 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_halffull_newton.cpp b/src/npair_halffull_newton.cpp index 407a71e614..d1e894943c 100644 --- a/src/npair_halffull_newton.cpp +++ b/src/npair_halffull_newton.cpp @@ -14,7 +14,9 @@ #include "npair_halffull_newton.h" #include "atom.h" +#include "domain.h" #include "error.h" +#include "force.h" #include "my_page.h" #include "neigh_list.h" @@ -37,6 +39,9 @@ void NPairHalffullNewton::build(NeighList *list) int *neighptr, *jlist; double xtmp, ytmp, ztmp; + double delta = 0.01 * force->angstrom; + int triclinic = domain->triclinic; + double **x = atom->x; int nlocal = atom->nlocal; @@ -72,8 +77,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; + } } else { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -81,6 +95,7 @@ void NPairHalffullNewton::build(NeighList *list) if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } } + neighptr[n++] = joriginal; } diff --git a/src/nstencil_half_bin_3d_tri.cpp b/src/nstencil_half_bin_3d_tri.cpp index d066a24ee6..f94bfc5e63 100644 --- a/src/nstencil_half_bin_3d_tri.cpp +++ b/src/nstencil_half_bin_3d_tri.cpp @@ -29,7 +29,8 @@ void NStencilHalfBin3dTri::create() nstencil = 0; - for (k = 0; k <= sz; k++) + //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(i, j, k) < cutneighmaxsq) From e3349581c757ad49a187bc86dd429abb4420dd1f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 9 Jul 2023 15:14:36 -0400 Subject: [PATCH 06/41] fix whitespace and remove debug code --- src/npair_half_bin_newton_tri.cpp | 56 +++++++++++-------------------- src/npair_half_nsq_newton.cpp | 24 ++++++------- src/npair_halffull_newton.cpp | 22 ++++++------ src/nstencil_half_bin_3d_tri.cpp | 1 - 4 files changed, 42 insertions(+), 61 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 17504ac3af..2917055214 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; @@ -92,40 +92,24 @@ void NPairHalfBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { - 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 (x[j][2] < ztmp) continue; - //if (x[j][2] == ztmp) { - // if (x[j][1] < ytmp) continue; - // if (x[j][1] == ytmp && x[j][0] < xtmp) 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; + } } } - */ - + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; @@ -135,8 +119,6 @@ void NPairHalfBinNewtonTri::build(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { - //printf("NEIGH i,j %d %d ijtag %d %d dist %g\n", - // i,j,tag[i],tag[j],sqrt(rsq)); if (molecular != Atom::ATOMIC) { if (!moltemplate) which = find_special(special[i],nspecial[i],tag[j]); diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index be06393f58..295b7de18c 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -42,9 +42,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; @@ -100,15 +100,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 d1e894943c..0192ed5729 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; @@ -77,17 +77,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) { @@ -95,7 +95,7 @@ void NPairHalffullNewton::build(NeighList *list) if (x[j][1] == ytmp && x[j][0] < xtmp) continue; } } - + neighptr[n++] = joriginal; } diff --git a/src/nstencil_half_bin_3d_tri.cpp b/src/nstencil_half_bin_3d_tri.cpp index f94bfc5e63..5887e389fb 100644 --- a/src/nstencil_half_bin_3d_tri.cpp +++ b/src/nstencil_half_bin_3d_tri.cpp @@ -29,7 +29,6 @@ void NStencilHalfBin3dTri::create() nstencil = 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++) From 9db5d4523282972f30dd6cb172fddcbb12921c1d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 9 Jul 2023 15:29:18 -0400 Subject: [PATCH 07/41] port neighbor list build changes to corresponding OPENMP package files --- src/OPENMP/npair_half_bin_newton_tri_omp.cpp | 29 ++++++++++++++------ src/OPENMP/npair_half_nsq_newton_omp.cpp | 24 ++++++++++++---- src/OPENMP/npair_halffull_newton_omp.cpp | 13 +++++++++ src/npair_half_nsq_newton.cpp | 7 +++-- 4 files changed, 56 insertions(+), 17 deletions(-) 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_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index 295b7de18c..023ece69c9 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; From 07f42930ff711b09bf46f9641ae964617edf3cad Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 10 Jul 2023 12:53:02 -0700 Subject: [PATCH 08/41] clean up code and comments --- src/compute_property_local.cpp | 3 ++- src/npair_half_bin_newton_tri.cpp | 29 ++++++----------------------- src/npair_half_nsq_newton.cpp | 5 +++++ src/npair_halffull_newton.cpp | 5 +++++ src/nstencil_half_bin_2d_tri.cpp | 12 ++++++++++-- src/nstencil_half_bin_3d_tri.cpp | 8 +++++++- 6 files changed, 35 insertions(+), 27 deletions(-) diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index d0523a1bec..92036c4bd2 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -405,7 +405,8 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag) if (!(mask[j] & groupbit)) continue; // itag = jtag is possible for long cutoffs that include images of self - + // do not need triclinic logic here b/c neighbor list itself is correct + if (newton_pair == 0 && j >= nlocal) { jtag = tag[j]; if (itag > jtag) { diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 17504ac3af..71a15df59e 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -83,11 +83,12 @@ void NPairHalfBinNewtonTri::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 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]) { @@ -107,25 +108,9 @@ void NPairHalfBinNewtonTri::build(NeighList *list) } else { if (x[j][0] < xtmp) continue; } - //if (x[j][2] < ztmp) continue; - //if (x[j][2] == ztmp) { - // if (x[j][1] < ytmp) continue; - // if (x[j][1] == ytmp && x[j][0] < xtmp) 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; - } - } - */ - jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; @@ -135,8 +120,6 @@ void NPairHalfBinNewtonTri::build(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { - //printf("NEIGH i,j %d %d ijtag %d %d dist %g\n", - // i,j,tag[i],tag[j],sqrt(rsq)); if (molecular != Atom::ATOMIC) { if (!moltemplate) which = find_special(special[i],nspecial[i],tag[j]); diff --git a/src/npair_half_nsq_newton.cpp b/src/npair_half_nsq_newton.cpp index be06393f58..0174a78900 100644 --- a/src/npair_half_nsq_newton.cpp +++ b/src/npair_half_nsq_newton.cpp @@ -89,7 +89,12 @@ void NPairHalfNsqNewton::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 direct I/J 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; diff --git a/src/npair_halffull_newton.cpp b/src/npair_halffull_newton.cpp index d1e894943c..af15b27eac 100644 --- a/src/npair_halffull_newton.cpp +++ b/src/npair_halffull_newton.cpp @@ -70,6 +70,11 @@ void NPairHalffullNewton::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 direct I/J 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]; diff --git a/src/nstencil_half_bin_2d_tri.cpp b/src/nstencil_half_bin_2d_tri.cpp index 06831730fd..920918fe09 100644 --- a/src/nstencil_half_bin_2d_tri.cpp +++ b/src/nstencil_half_bin_2d_tri.cpp @@ -27,9 +27,17 @@ void NStencilHalfBin2dTri::create() { int i, j; + // for triclinic, need to use full stencil in all dims + // not a half stencil in y + // b/c transforming orthog -> lambda -> orthog for ghost atoms + // with an added PBC offset can shift both coords by epsilon + // thus for an I/J owned/ghost pair, the xy coords + // and bin assignments can be different on I proc vs J proc + nstencil = 0; - for (j = 0; j <= sy; j++) + for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance(i, j, 0) < cutneighmaxsq) stencil[nstencil++] = j * mbinx + i; + if (bin_distance(i, j, 0) < cutneighmaxsq) + stencil[nstencil++] = j * mbinx + i; } diff --git a/src/nstencil_half_bin_3d_tri.cpp b/src/nstencil_half_bin_3d_tri.cpp index f94bfc5e63..d146b92cd1 100644 --- a/src/nstencil_half_bin_3d_tri.cpp +++ b/src/nstencil_half_bin_3d_tri.cpp @@ -27,9 +27,15 @@ void NStencilHalfBin3dTri::create() { int i, j, k; + // for triclinic, need to use full stencil in all dims + // not a half stencil in z + // b/c transforming orthog -> lambda -> orthog for ghost atoms + // with an added PBC offset can shift all 3 coords by epsilon + // thus for an I/J owned/ghost pair, the xyz coords + // and bin assignments can be different on I proc vs J proc + nstencil = 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++) From fbbf44fb8e2948fea648f048f4ab1cb088cd93bc Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 10 Jul 2023 18:25:29 -0700 Subject: [PATCH 09/41] same changes to other NPair and NStencil methods --- src/npair_half_bin_newton_tri.cpp | 4 +- src/npair_half_multi_newton_tri.cpp | 119 +++++++++++-------- src/npair_half_multi_old_newton_tri.cpp | 44 ++++--- src/npair_half_respa_bin_newton_tri.cpp | 45 ++++--- src/npair_half_respa_nsq_newton.cpp | 24 +++- src/npair_half_respa_nsq_newton.h | 2 +- src/npair_half_size_bin_newton_tri.cpp | 41 +++++-- src/npair_half_size_multi_newton_tri.cpp | 107 ++++++++++------- src/npair_half_size_multi_old_newton_tri.cpp | 45 ++++--- src/npair_half_size_nsq_newton.cpp | 22 +++- src/npair_halffull_newton_trim.cpp | 19 +++ src/nstencil_half_multi_2d_tri.cpp | 2 +- src/nstencil_half_multi_3d_tri.cpp | 2 +- src/nstencil_half_multi_old_2d_tri.cpp | 2 +- src/nstencil_half_multi_old_3d_tri.cpp | 2 +- 15 files changed, 313 insertions(+), 167 deletions(-) diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index f6cbf1b1af..9c0688af68 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -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; diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index 9bebfe71e2..316acb5049 100644 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -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; diff --git a/src/npair_half_multi_old_newton_tri.cpp b/src/npair_half_multi_old_newton_tri.cpp index fbb9a8e504..9dcbcff9f4 100644 --- a/src/npair_half_multi_old_newton_tri.cpp +++ b/src/npair_half_multi_old_newton_tri.cpp @@ -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; diff --git a/src/npair_half_respa_bin_newton_tri.cpp b/src/npair_half_respa_bin_newton_tri.cpp index b2749bd7a7..05b839869a 100644 --- a/src/npair_half_respa_bin_newton_tri.cpp +++ b/src/npair_half_respa_bin_newton_tri.cpp @@ -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; diff --git a/src/npair_half_respa_nsq_newton.cpp b/src/npair_half_respa_nsq_newton.cpp index 77d6af141f..d0292eec92 100644 --- a/src/npair_half_respa_nsq_newton.cpp +++ b/src/npair_half_respa_nsq_newton.cpp @@ -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) { diff --git a/src/npair_half_respa_nsq_newton.h b/src/npair_half_respa_nsq_newton.h index e5233f5e9d..4a5ae23aef 100644 --- a/src/npair_half_respa_nsq_newton.h +++ b/src/npair_half_respa_nsq_newton.h @@ -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 diff --git a/src/npair_half_size_bin_newton_tri.cpp b/src/npair_half_size_bin_newton_tri.cpp index 47bb9d01e1..e6a236eecb 100644 --- a/src/npair_half_size_bin_newton_tri.cpp +++ b/src/npair_half_size_bin_newton_tri.cpp @@ -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; diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index 5d8a0f05ef..a363ae6e1e 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -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; diff --git a/src/npair_half_size_multi_old_newton_tri.cpp b/src/npair_half_size_multi_old_newton_tri.cpp index ea3f271956..974500d6b8 100644 --- a/src/npair_half_size_multi_old_newton_tri.cpp +++ b/src/npair_half_size_multi_old_newton_tri.cpp @@ -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; diff --git a/src/npair_half_size_nsq_newton.cpp b/src/npair_half_size_nsq_newton.cpp index 8b596e6968..abd2a4faff 100644 --- a/src/npair_half_size_nsq_newton.cpp +++ b/src/npair_half_size_nsq_newton.cpp @@ -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) { diff --git a/src/npair_halffull_newton_trim.cpp b/src/npair_halffull_newton_trim.cpp index b7bb72c990..7d420f88af 100644 --- a/src/npair_halffull_newton_trim.cpp +++ b/src/npair_halffull_newton_trim.cpp @@ -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) { diff --git a/src/nstencil_half_multi_2d_tri.cpp b/src/nstencil_half_multi_2d_tri.cpp index bf39c04099..85bbe94c86 100644 --- a/src/nstencil_half_multi_2d_tri.cpp +++ b/src/nstencil_half_multi_2d_tri.cpp @@ -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; diff --git a/src/nstencil_half_multi_3d_tri.cpp b/src/nstencil_half_multi_3d_tri.cpp index f2d4d051ad..9761e15854 100644 --- a/src/nstencil_half_multi_3d_tri.cpp +++ b/src/nstencil_half_multi_3d_tri.cpp @@ -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) diff --git a/src/nstencil_half_multi_old_2d_tri.cpp b/src/nstencil_half_multi_old_2d_tri.cpp index 1438aef843..0aeb65bebd 100644 --- a/src/nstencil_half_multi_old_2d_tri.cpp +++ b/src/nstencil_half_multi_old_2d_tri.cpp @@ -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) { diff --git a/src/nstencil_half_multi_old_3d_tri.cpp b/src/nstencil_half_multi_old_3d_tri.cpp index 836eee6039..3717b7836b 100644 --- a/src/nstencil_half_multi_old_3d_tri.cpp +++ b/src/nstencil_half_multi_old_3d_tri.cpp @@ -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); From 6bd965f0df72b23b1a4819f9382eefc642619e72 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 10 Jul 2023 22:35:36 -0400 Subject: [PATCH 10/41] fix whitespace (again) --- src/compute_property_local.cpp | 2 +- src/npair_half_bin_newton_tri.cpp | 36 +++---- src/npair_half_multi_newton_tri.cpp | 98 +++++++++--------- src/npair_half_multi_old_newton_tri.cpp | 32 +++--- src/npair_half_respa_bin_newton_tri.cpp | 36 +++---- src/npair_half_size_bin_newton_tri.cpp | 34 +++---- src/npair_half_size_multi_newton_tri.cpp | 100 +++++++++---------- src/npair_half_size_multi_old_newton_tri.cpp | 34 +++---- src/npair_halffull_newton_trim.cpp | 2 +- src/nstencil_half_bin_3d_tri.cpp | 2 +- 10 files changed, 188 insertions(+), 188 deletions(-) diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 92036c4bd2..87517a3e05 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -406,7 +406,7 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag) // itag = jtag is possible for long cutoffs that include images of self // do not need triclinic logic here b/c neighbor list itself is correct - + if (newton_pair == 0 && j >= nlocal) { jtag = tag[j]; if (itag > jtag) { diff --git a/src/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 9c0688af68..453d10096e 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -88,28 +88,28 @@ void NPairHalfBinNewtonTri::build(NeighList *list) // 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 (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 (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; diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index 316acb5049..1d75d6a3ef 100644 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -91,11 +91,11 @@ 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); @@ -109,60 +109,60 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) // 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]; 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; - } - } - } - } + 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), exclude half of interactions - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; + 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; + } + } + } + } - 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; - } - } + 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; + } + } } } diff --git a/src/npair_half_multi_old_newton_tri.cpp b/src/npair_half_multi_old_newton_tri.cpp index 9dcbcff9f4..72d46d042f 100644 --- a/src/npair_half_multi_old_newton_tri.cpp +++ b/src/npair_half_multi_old_newton_tri.cpp @@ -102,22 +102,22 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list) jtype = type[j]; if (cutsq[jtype] < distsq[k]) 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 (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; diff --git a/src/npair_half_respa_bin_newton_tri.cpp b/src/npair_half_respa_bin_newton_tri.cpp index 05b839869a..eac67b8bd5 100644 --- a/src/npair_half_respa_bin_newton_tri.cpp +++ b/src/npair_half_respa_bin_newton_tri.cpp @@ -118,24 +118,24 @@ void NPairHalfRespaBinNewtonTri::build(NeighList *list) ibin = atom2bin[i]; for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { - - 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 (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; diff --git a/src/npair_half_size_bin_newton_tri.cpp b/src/npair_half_size_bin_newton_tri.cpp index e6a236eecb..0d1a0a7329 100644 --- a/src/npair_half_size_bin_newton_tri.cpp +++ b/src/npair_half_size_bin_newton_tri.cpp @@ -101,23 +101,23 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) for (k = 0; k < nstencil; k++) { for (j = binhead[ibin+stencil[k]]; j >= 0; j = bins[j]) { - 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 (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; diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index a363ae6e1e..f597789dee 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -98,11 +98,11 @@ 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); @@ -119,61 +119,61 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) 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]) { + 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 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; - } - } - } + 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; + } + } } } diff --git a/src/npair_half_size_multi_old_newton_tri.cpp b/src/npair_half_size_multi_old_newton_tri.cpp index 974500d6b8..848a19aa39 100644 --- a/src/npair_half_size_multi_old_newton_tri.cpp +++ b/src/npair_half_size_multi_old_newton_tri.cpp @@ -108,23 +108,23 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list) jtype = type[j]; if (cutsq[jtype] < distsq[k]) 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 (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; diff --git a/src/npair_halffull_newton_trim.cpp b/src/npair_halffull_newton_trim.cpp index 7d420f88af..e758c04284 100644 --- a/src/npair_halffull_newton_trim.cpp +++ b/src/npair_halffull_newton_trim.cpp @@ -85,7 +85,7 @@ 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) { diff --git a/src/nstencil_half_bin_3d_tri.cpp b/src/nstencil_half_bin_3d_tri.cpp index d146b92cd1..72bef7fb76 100644 --- a/src/nstencil_half_bin_3d_tri.cpp +++ b/src/nstencil_half_bin_3d_tri.cpp @@ -33,7 +33,7 @@ void NStencilHalfBin3dTri::create() // with an added PBC offset can shift all 3 coords by epsilon // thus for an I/J owned/ghost pair, the xyz coords // and bin assignments can be different on I proc vs J proc - + nstencil = 0; for (k = -sz; k <= sz; k++) From 2eeea4332076c4a11193bd6c4a71d29ab649778c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 11 Jul 2023 01:09:44 -0400 Subject: [PATCH 11/41] port neighbor list changes to OPENMP package --- src/OPENMP/npair_half_bin_newton_tri_omp.cpp | 20 +-- .../npair_half_multi_newton_tri_omp.cpp | 118 ++++++++++-------- .../npair_half_multi_old_newton_tri_omp.cpp | 50 +++++--- src/OPENMP/npair_half_nsq_newton_omp.cpp | 7 +- .../npair_half_respa_bin_newton_tri_omp.cpp | 47 ++++--- .../npair_half_respa_nsq_newton_omp.cpp | 35 ++++-- .../npair_half_size_bin_newton_tri_omp.cpp | 40 +++--- .../npair_half_size_multi_newton_tri_omp.cpp | 43 ++++--- ...air_half_size_multi_old_newton_tri_omp.cpp | 42 ++++--- src/OPENMP/npair_half_size_nsq_newton_omp.cpp | 27 +++- src/npair_half_bin_newton_tri.cpp | 3 +- src/npair_half_multi_newton_tri.cpp | 4 +- src/npair_half_respa_bin_newton_tri.cpp | 3 +- src/npair_half_respa_nsq_newton.cpp | 7 +- src/npair_half_size_nsq_newton.cpp | 2 +- 15 files changed, 278 insertions(+), 170 deletions(-) diff --git a/src/OPENMP/npair_half_bin_newton_tri_omp.cpp b/src/OPENMP/npair_half_bin_newton_tri_omp.cpp index 3ad07acd56..47524474ed 100644 --- a/src/OPENMP/npair_half_bin_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_bin_newton_tri_omp.cpp @@ -12,17 +12,18 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "omp_compat.h" #include "npair_half_bin_newton_tri_omp.h" #include "npair_omp.h" -#include "neigh_list.h" +#include "omp_compat.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" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -54,8 +55,6 @@ void NPairHalfBinNewtonTriOmp::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - // loop over each atom, storing neighbors - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -93,10 +92,11 @@ void NPairHalfBinNewtonTriOmp::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++) { diff --git a/src/OPENMP/npair_half_multi_newton_tri_omp.cpp b/src/OPENMP/npair_half_multi_newton_tri_omp.cpp index a152d011a7..b18cba0261 100644 --- a/src/OPENMP/npair_half_multi_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_multi_newton_tri_omp.cpp @@ -12,17 +12,19 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ -#include "omp_compat.h" #include "npair_half_multi_newton_tri_omp.h" -#include "npair_omp.h" -#include "neighbor.h" -#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" +#include "neigh_list.h" +#include "neighbor.h" +#include "npair_omp.h" +#include "omp_compat.h" using namespace LAMMPS_NS; @@ -43,6 +45,7 @@ void NPairHalfMultiNewtonTriOmp::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) @@ -51,13 +54,11 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,k,n,itype,jtype,ibin,jbin,icollection,jcollection,which,ns,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; int js; - // loop over each atom, storing neighbors - int *collection = neighbor->collection; double **x = atom->x; int *type = atom->type; @@ -84,6 +85,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; itype = type[i]; icollection = collection[i]; xtmp = x[i][0]; @@ -98,65 +100,79 @@ void NPairHalfMultiNewtonTriOmp::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]) { + 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; + 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; + 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[i] = i; diff --git a/src/OPENMP/npair_half_multi_old_newton_tri_omp.cpp b/src/OPENMP/npair_half_multi_old_newton_tri_omp.cpp index e4895ff1a9..38f645abad 100644 --- a/src/OPENMP/npair_half_multi_old_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_multi_old_newton_tri_omp.cpp @@ -15,13 +15,15 @@ #include "omp_compat.h" #include "npair_half_multi_old_newton_tri_omp.h" #include "npair_omp.h" -#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" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -42,6 +44,7 @@ void NPairHalfMultiOldNewtonTriOmp::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) @@ -50,13 +53,11 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; - // loop over each atom, storing neighbors - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -82,6 +83,7 @@ void NPairHalfMultiOldNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -92,13 +94,12 @@ void NPairHalfMultiOldNewtonTriOmp::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]; @@ -109,12 +110,21 @@ void NPairHalfMultiOldNewtonTriOmp::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; + } } } @@ -129,7 +139,7 @@ void NPairHalfMultiOldNewtonTriOmp::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 726814c6f0..42cf63278a 100644 --- a/src/OPENMP/npair_half_nsq_newton_omp.cpp +++ b/src/OPENMP/npair_half_nsq_newton_omp.cpp @@ -58,8 +58,6 @@ void NPairHalfNsqNewtonOmp::build(NeighList *list) double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - // loop over each atom, storing neighbors - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -99,7 +97,12 @@ void NPairHalfNsqNewtonOmp::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; diff --git a/src/OPENMP/npair_half_respa_bin_newton_tri_omp.cpp b/src/OPENMP/npair_half_respa_bin_newton_tri_omp.cpp index c998f71290..78b3abdd66 100644 --- a/src/OPENMP/npair_half_respa_bin_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_respa_bin_newton_tri_omp.cpp @@ -15,13 +15,15 @@ #include "omp_compat.h" #include "npair_half_respa_bin_newton_tri_omp.h" #include "npair_omp.h" -#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" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -42,6 +44,7 @@ void NPairHalfRespaBinNewtonTriOmp::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; @@ -53,12 +56,10 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; - // loop over each atom, storing neighbors - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -111,6 +112,7 @@ void NPairHalfRespaBinNewtonTriOmp::build(NeighList *list) neighptr_middle = ipage_middle->vget(); } + itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -122,20 +124,31 @@ void NPairHalfRespaBinNewtonTriOmp::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; + } } } @@ -151,7 +164,7 @@ void NPairHalfRespaBinNewtonTriOmp::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_respa_nsq_newton_omp.cpp b/src/OPENMP/npair_half_respa_nsq_newton_omp.cpp index 6604861f74..a9745edc64 100644 --- a/src/OPENMP/npair_half_respa_nsq_newton_omp.cpp +++ b/src/OPENMP/npair_half_respa_nsq_newton_omp.cpp @@ -15,21 +15,22 @@ #include "omp_compat.h" #include "npair_half_respa_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; /* ---------------------------------------------------------------------- */ -NPairHalfRespaNsqNewtonOmp::NPairHalfRespaNsqNewtonOmp(LAMMPS *lmp) : - NPair(lmp) {} +NPairHalfRespaNsqNewtonOmp::NPairHalfRespaNsqNewtonOmp(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- multiple respa lists @@ -45,6 +46,8 @@ void NPairHalfRespaNsqNewtonOmp::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; @@ -55,13 +58,11 @@ void NPairHalfRespaNsqNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle,imol,iatom; - tagint tagprev; + int i,j,n,itype,jtype,n_inner,n_middle,imol,iatom; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; - // loop over each atom, storing neighbors - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; @@ -128,6 +129,12 @@ void NPairHalfRespaNsqNewtonOmp::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; @@ -138,6 +145,14 @@ void NPairHalfRespaNsqNewtonOmp::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) { @@ -159,7 +174,7 @@ void NPairHalfRespaNsqNewtonOmp::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_size_bin_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp index c320296442..7fcf07e9c8 100644 --- a/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp @@ -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" @@ -46,6 +47,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; const int mask_history = 1 << HISTBITS; + const double delta = 0.01 * force->angstrom; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -54,13 +56,11 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,jh,k,n,ibin,which,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; - // loop over each atom, storing neighbors - double **x = atom->x; double *radius = atom->radius; int *type = atom->type; @@ -87,6 +87,7 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -98,20 +99,31 @@ void NPairHalfSizeBinNewtonTriOmp::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; + } } } @@ -132,7 +144,7 @@ void NPairHalfSizeBinNewtonTriOmp::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_size_multi_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp index 9a0ead482b..916b7bfbc3 100644 --- a/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp @@ -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" @@ -48,6 +49,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; const int mask_history = 1 << HISTBITS; + const double delta = 0.01 * force->angstrom; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -55,15 +57,12 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js; int which,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; - int js; - - // loop over each atom, storing neighbors int *collection = neighbor->collection; double **x = atom->x; @@ -92,6 +91,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; itype = type[i]; icollection = collection[i]; xtmp = x[i][0]; @@ -107,12 +107,13 @@ void NPairHalfSizeMultiNewtonTriOmp::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); + if (icollection == jcollection) jbin = ibin; + else jbin = coord2bin(x[i], jcollection); // loop over all atoms in bins in stencil // stencil is empty if i larger than j @@ -130,14 +131,24 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) 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; + } } } } @@ -160,7 +171,7 @@ void NPairHalfSizeMultiNewtonTriOmp::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_size_multi_old_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp index c74b191f66..7faa210107 100644 --- a/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp @@ -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" @@ -32,7 +33,6 @@ NPairHalfSizeMultiOldNewtonTriOmp::NPairHalfSizeMultiOldNewtonTriOmp(LAMMPS *lmp NPair(lmp) {} /* ---------------------------------------------------------------------- - size particles binned neighbor list construction with Newton's 3rd law for triclinic each owned atom i checks its own bin and other bins in triclinic stencil multi-type stencil is itype dependent and is distance checked @@ -46,6 +46,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; const int mask_history = 1 << HISTBITS; + const double delta = 0.01 * force->angstrom; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -54,7 +55,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) NPAIR_OMP_SETUP(nlocal); int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom; - tagint tagprev; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -97,13 +98,12 @@ void NPairHalfSizeMultiOldNewtonTriOmp::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]; @@ -114,12 +114,22 @@ void NPairHalfSizeMultiOldNewtonTriOmp::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; + } } } @@ -140,7 +150,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::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_size_nsq_newton_omp.cpp b/src/OPENMP/npair_half_size_nsq_newton_omp.cpp index 35dc42ec5b..0628478c0b 100644 --- a/src/OPENMP/npair_half_size_nsq_newton_omp.cpp +++ b/src/OPENMP/npair_half_size_nsq_newton_omp.cpp @@ -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" @@ -30,13 +31,11 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) : - NPair(lmp) {} +NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles N^2 / 2 search for neighbor pairs with full Newton's 3rd law - shear history must be accounted for when a neighbor pair is added pair added to list if atoms i and j are both owned and i < j if j is ghost only me or other proc adds pair decision based on itag,jtag tests @@ -50,6 +49,8 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; const int mask_history = 1 << HISTBITS; + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; NPAIR_OMP_INIT; @@ -58,8 +59,8 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,jh,n,itag,jtag,which,imol,iatom; - tagint tagprev; + int i,j,jh,n,which,imol,iatom; + tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -104,6 +105,12 @@ void NPairHalfSizeNsqNewtonOmp::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; @@ -114,6 +121,14 @@ void NPairHalfSizeNsqNewtonOmp::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) { @@ -140,7 +155,7 @@ void NPairHalfSizeNsqNewtonOmp::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/npair_half_bin_newton_tri.cpp b/src/npair_half_bin_newton_tri.cpp index 453d10096e..d261363b0e 100644 --- a/src/npair_half_bin_newton_tri.cpp +++ b/src/npair_half_bin_newton_tri.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ #include "npair_half_bin_newton_tri.h" -#include "neigh_list.h" + #include "atom.h" #include "atom_vec.h" #include "domain.h" @@ -21,6 +21,7 @@ #include "force.h" #include "molecule.h" #include "my_page.h" +#include "neigh_list.h" using namespace LAMMPS_NS; diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index 1d75d6a3ef..2b753af499 100644 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -21,8 +21,8 @@ #include "force.h" #include "molecule.h" #include "my_page.h" -#include "neighbor.h" #include "neigh_list.h" +#include "neighbor.h" using namespace LAMMPS_NS; @@ -39,7 +39,7 @@ 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; + int i,j,k,n,itype,jtype,ibin,jbin,icollection,jcollection,which,ns,imol,iatom,moltemplate; tagint itag,jtag,tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; diff --git a/src/npair_half_respa_bin_newton_tri.cpp b/src/npair_half_respa_bin_newton_tri.cpp index eac67b8bd5..4cd4ead0fa 100644 --- a/src/npair_half_respa_bin_newton_tri.cpp +++ b/src/npair_half_respa_bin_newton_tri.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ #include "npair_half_respa_bin_newton_tri.h" -#include "neigh_list.h" + #include "atom.h" #include "atom_vec.h" #include "domain.h" @@ -21,6 +21,7 @@ #include "force.h" #include "molecule.h" #include "my_page.h" +#include "neigh_list.h" using namespace LAMMPS_NS; diff --git a/src/npair_half_respa_nsq_newton.cpp b/src/npair_half_respa_nsq_newton.cpp index d0292eec92..ae56d62fb5 100644 --- a/src/npair_half_respa_nsq_newton.cpp +++ b/src/npair_half_respa_nsq_newton.cpp @@ -13,15 +13,16 @@ ------------------------------------------------------------------------- */ #include "npair_half_respa_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 "my_page.h" -#include "error.h" +#include "neigh_list.h" using namespace LAMMPS_NS; @@ -132,7 +133,7 @@ void NPairHalfRespaNsqNewton::build(NeighList *list) if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; - } else if (triclinic) { + } else if (triclinic) { if (fabs(x[j][2]-ztmp) > delta) { if (x[j][2] < ztmp) continue; } else if (fabs(x[j][1]-ytmp) > delta) { diff --git a/src/npair_half_size_nsq_newton.cpp b/src/npair_half_size_nsq_newton.cpp index abd2a4faff..ce0c7f9562 100644 --- a/src/npair_half_size_nsq_newton.cpp +++ b/src/npair_half_size_nsq_newton.cpp @@ -113,7 +113,7 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; - } else if (triclinic) { + } else if (triclinic) { if (fabs(x[j][2]-ztmp) > delta) { if (x[j][2] < ztmp) continue; } else if (fabs(x[j][1]-ytmp) > delta) { From bb6e4d844088cea152ebd0cf5007a6ba5c2cbfe2 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 11 Jul 2023 14:30:14 -0700 Subject: [PATCH 12/41] add forgotten line to multi and mutli/old --- src/npair_half_multi_newton_tri.cpp | 1 + src/npair_half_multi_old_newton_tri.cpp | 1 + src/npair_half_size_multi_newton_tri.cpp | 1 + 3 files changed, 3 insertions(+) diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index 316acb5049..2595c3ce21 100644 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -121,6 +121,7 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]) { + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { diff --git a/src/npair_half_multi_old_newton_tri.cpp b/src/npair_half_multi_old_newton_tri.cpp index 9dcbcff9f4..8700db3b32 100644 --- a/src/npair_half_multi_old_newton_tri.cpp +++ b/src/npair_half_multi_old_newton_tri.cpp @@ -102,6 +102,7 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list) jtype = type[j]; if (cutsq[jtype] < distsq[k]) continue; + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index a363ae6e1e..4ff50870b6 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -126,6 +126,7 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]) { + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { From 2a7ac115d8563dab636f10ba516c461ba23f0694 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 11 Jul 2023 18:25:29 -0400 Subject: [PATCH 13/41] fix whitespace (one more time) --- src/npair_half_multi_newton_tri.cpp | 50 ++++++++++++------------ src/npair_half_multi_old_newton_tri.cpp | 34 ++++++++-------- src/npair_half_size_multi_newton_tri.cpp | 36 ++++++++--------- 3 files changed, 60 insertions(+), 60 deletions(-) diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index 1c95e73151..24300f6929 100644 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -114,31 +114,31 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) 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 <= 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; - } - } - } - } + 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 <= 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; diff --git a/src/npair_half_multi_old_newton_tri.cpp b/src/npair_half_multi_old_newton_tri.cpp index 8700db3b32..ce3149ebf5 100644 --- a/src/npair_half_multi_old_newton_tri.cpp +++ b/src/npair_half_multi_old_newton_tri.cpp @@ -102,23 +102,23 @@ void NPairHalfMultiOldNewtonTri::build(NeighList *list) jtype = type[j]; if (cutsq[jtype] < distsq[k]) 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 (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; diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index c521034406..aa0d8e3f42 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -125,24 +125,24 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) // if same size (same collection), exclude half of interactions if (cutcollectionsq[icollection][icollection] == - cutcollectionsq[jcollection][jcollection]) { - 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; - } - } - } + cutcollectionsq[jcollection][jcollection]) { + 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]; From 89fb236144902530091f279d012b90f907b93b50 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 11 Jul 2023 18:35:16 -0400 Subject: [PATCH 14/41] port bugfix for colloid test failure --- src/OPENMP/npair_half_multi_newton_tri_omp.cpp | 1 + src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/src/OPENMP/npair_half_multi_newton_tri_omp.cpp b/src/OPENMP/npair_half_multi_newton_tri_omp.cpp index b18cba0261..e26bea990f 100644 --- a/src/OPENMP/npair_half_multi_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_multi_newton_tri_omp.cpp @@ -130,6 +130,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]) { + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { diff --git a/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp index 916b7bfbc3..4765c918b7 100644 --- a/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp @@ -135,6 +135,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) if (cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]) { + if (j <= i) continue; if (j >= nlocal) { jtag = tag[j]; if (itag > jtag) { From 3fc809a1b9782430a4e41b8ca1d8a3c224762ef6 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 11 Jul 2023 18:45:45 -0700 Subject: [PATCH 15/41] add check for atom IDs when triclinic --- src/neighbor.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/neighbor.cpp b/src/neighbor.cpp index df1547e5eb..c6eea7e2f1 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -313,7 +313,10 @@ void Neighbor::init() triclinic = domain->triclinic; newton_pair = force->newton_pair; - // error check + // error checks + + if (triclinic && atom->tag_enable == 0) + error->all(FLERR, "Cannot build triclinic neighbor lists unless atoms have IDs"); if (delay > 0 && (delay % every) != 0) error->all(FLERR,"Neighbor delay must be 0 or multiple of every setting"); From a91b3dab963d0044885298955cc2edf8e6556ead Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 11 Jul 2023 18:50:23 -0700 Subject: [PATCH 16/41] doc atom ID requirement for triclinic --- doc/src/Howto_triclinic.rst | 3 ++- doc/src/atom_modify.rst | 5 +++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/doc/src/Howto_triclinic.rst b/doc/src/Howto_triclinic.rst index 0efadbcc8c..2983d013c6 100644 --- a/doc/src/Howto_triclinic.rst +++ b/doc/src/Howto_triclinic.rst @@ -12,7 +12,8 @@ is created, e.g. by the :doc:`create_box ` or :doc:`read_data ` or :doc:`read_restart ` commands. Additionally, LAMMPS defines box size parameters lx,ly,lz where lx = xhi-xlo, and similarly in the y and z dimensions. The 6 -parameters, as well as lx,ly,lz, can be output via the :doc:`thermo_style custom ` command. +parameters, as well as lx,ly,lz, can be output via the +:doc:`thermo_style custom ` command. LAMMPS also allows simulations to be performed in triclinic (non-orthogonal) simulation boxes shaped as a parallelepiped with diff --git a/doc/src/atom_modify.rst b/doc/src/atom_modify.rst index 1e5a3d49ff..21590e6680 100644 --- a/doc/src/atom_modify.rst +++ b/doc/src/atom_modify.rst @@ -65,6 +65,11 @@ switch. This is described on the :doc:`Build_settings ` doc page. If atom IDs are not used, they must be specified as 0 for all atoms, e.g. in a data or restart file. +.. note:: + + If a :doc:`triclinic simulation box ` is used, + atom IDs are required, due to how neighbor lists are built. + The *map* keyword determines how atoms with specific IDs are found when required. An example are the bond (angle, etc) methods which need to find the local index of an atom with a specific global ID From bc6fcdc61a550b645b25f12138a3bb4022153c9a Mon Sep 17 00:00:00 2001 From: "W. Michael Brown" Date: Fri, 4 Aug 2023 08:49:27 -0700 Subject: [PATCH 17/41] Applying triclinic neighbor fixes to intel package. --- src/INTEL/fix_intel.cpp | 2 + src/INTEL/npair_halffull_newton_intel.cpp | 66 ++++++++++---- .../npair_halffull_newton_trim_intel.cpp | 87 ++++++++++++++----- src/INTEL/npair_intel.cpp | 40 +++++++-- 4 files changed, 144 insertions(+), 51 deletions(-) diff --git a/src/INTEL/fix_intel.cpp b/src/INTEL/fix_intel.cpp index 4c46608677..8396904ffd 100644 --- a/src/INTEL/fix_intel.cpp +++ b/src/INTEL/fix_intel.cpp @@ -20,6 +20,7 @@ #include "fix_intel.h" #include "comm.h" +#include "domain.h" #include "error.h" #include "force.h" #include "neighbor.h" @@ -470,6 +471,7 @@ void FixIntel::pair_init_check(const bool cdmessage) int need_tag = 0; if (atom->molecular != Atom::ATOMIC || three_body_neighbor()) need_tag = 1; + if (domain->triclinic && force->newton_pair) need_tag = 1; // Clear buffers used for pair style char kmode[80]; diff --git a/src/INTEL/npair_halffull_newton_intel.cpp b/src/INTEL/npair_halffull_newton_intel.cpp index cd05d5f97a..adcf2527ab 100644 --- a/src/INTEL/npair_halffull_newton_intel.cpp +++ b/src/INTEL/npair_halffull_newton_intel.cpp @@ -20,7 +20,9 @@ #include "atom.h" #include "comm.h" +#include "domain.h" #include "error.h" +#include "force.h" #include "modify.h" #include "my_page.h" #include "neigh_list.h" @@ -56,6 +58,9 @@ void NPairHalffullNewtonIntel::build_t(NeighList *list, const int * _noalias const numneigh_full = list->listfull->numneigh; const int ** _noalias const firstneigh_full = (const int ** const)list->listfull->firstneigh; // NOLINT + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; + #if defined(_OPENMP) #pragma omp parallel #endif @@ -82,25 +87,50 @@ void NPairHalffullNewtonIntel::build_t(NeighList *list, const int * _noalias const jlist = firstneigh_full[i]; const int jnum = numneigh_full[i]; - #if defined(LMP_SIMD_COMPILER) - #pragma vector aligned - #pragma ivdep - #endif - for (int jj = 0; jj < jnum; jj++) { - const int joriginal = jlist[jj]; - const int j = joriginal & NEIGHMASK; - int addme = 1; - if (j < nlocal) { - if (i > j) addme = 0; - } else { - if (x[j].z < ztmp) addme = 0; - if (x[j].z == ztmp) { - if (x[j].y < ytmp) addme = 0; - if (x[j].y == ytmp && x[j].x < xtmp) addme = 0; + if (!triclinic) { + #if defined(LMP_SIMD_COMPILER) + #pragma vector aligned + #pragma ivdep + #endif + for (int jj = 0; jj < jnum; jj++) { + const int joriginal = jlist[jj]; + const int j = joriginal & NEIGHMASK; + int addme = 1; + if (j < nlocal) { + if (i > j) addme = 0; + } else { + if (x[j].z < ztmp) addme = 0; + if (x[j].z == ztmp) { + if (x[j].y < ytmp) addme = 0; + if (x[j].y == ytmp && x[j].x < xtmp) addme = 0; + } } + if (addme) + neighptr[n++] = joriginal; + } + } else { + #if defined(LMP_SIMD_COMPILER) + #pragma vector aligned + #pragma ivdep + #endif + for (int jj = 0; jj < jnum; jj++) { + const int joriginal = jlist[jj]; + const int j = joriginal & NEIGHMASK; + int addme = 1; + if (j < nlocal) { + if (i > j) addme = 0; + } else { + if (fabs(x[j].z-ztmp) > delta) { + if (x[j].z < ztmp) addme = 0; + } else if (fabs(x[j].y-ytmp) > delta) { + if (x[j].y < ytmp) addme = 0; + } else { + if (x[j].x < xtmp) addme = 0; + } + } + if (addme) + neighptr[n++] = joriginal; } - if (addme) - neighptr[n++] = joriginal; } ilist[ii] = i; @@ -203,7 +233,7 @@ void NPairHalffullNewtonIntel::build_t3(NeighList *list, int *numhalf) void NPairHalffullNewtonIntel::build(NeighList *list) { - if (_fix->three_body_neighbor() == 0) { + if (_fix->three_body_neighbor() == 0 || domain->triclinic) { if (_fix->precision() == FixIntel::PREC_MODE_MIXED) build_t(list, _fix->get_mixed_buffers()); else if (_fix->precision() == FixIntel::PREC_MODE_DOUBLE) diff --git a/src/INTEL/npair_halffull_newton_trim_intel.cpp b/src/INTEL/npair_halffull_newton_trim_intel.cpp index e38375f750..34b9b20e9c 100644 --- a/src/INTEL/npair_halffull_newton_trim_intel.cpp +++ b/src/INTEL/npair_halffull_newton_trim_intel.cpp @@ -20,7 +20,9 @@ #include "atom.h" #include "comm.h" +#include "domain.h" #include "error.h" +#include "force.h" #include "modify.h" #include "my_page.h" #include "neigh_list.h" @@ -57,6 +59,8 @@ void NPairHalffullNewtonTrimIntel::build_t(NeighList *list, const int ** _noalias const firstneigh_full = (const int ** const)list->listfull->firstneigh; // NOLINT const flt_t cutsq_custom = cutoff_custom * cutoff_custom; + const double delta = 0.01 * force->angstrom; + const int triclinic = domain->triclinic; #if defined(_OPENMP) #pragma omp parallel @@ -84,35 +88,70 @@ void NPairHalffullNewtonTrimIntel::build_t(NeighList *list, const int * _noalias const jlist = firstneigh_full[i]; const int jnum = numneigh_full[i]; - #if defined(LMP_SIMD_COMPILER) - #pragma vector aligned - #pragma ivdep - #endif - for (int jj = 0; jj < jnum; jj++) { - const int joriginal = jlist[jj]; - const int j = joriginal & NEIGHMASK; - int addme = 1; - if (j < nlocal) { - if (i > j) addme = 0; - } else { - if (x[j].z < ztmp) addme = 0; - if (x[j].z == ztmp) { - if (x[j].y < ytmp) addme = 0; - if (x[j].y == ytmp && x[j].x < xtmp) addme = 0; + if (!triclinic) { + #if defined(LMP_SIMD_COMPILER) + #pragma vector aligned + #pragma ivdep + #endif + for (int jj = 0; jj < jnum; jj++) { + const int joriginal = jlist[jj]; + const int j = joriginal & NEIGHMASK; + int addme = 1; + if (j < nlocal) { + if (i > j) addme = 0; + } else { + if (x[j].z < ztmp) addme = 0; + if (x[j].z == ztmp) { + if (x[j].y < ytmp) addme = 0; + if (x[j].y == ytmp && x[j].x < xtmp) addme = 0; + } } + + // trim to shorter cutoff + + const flt_t delx = xtmp - x[j].x; + const flt_t dely = ytmp - x[j].y; + const flt_t delz = ztmp - x[j].z; + const flt_t rsq = delx * delx + dely * dely + delz * delz; + + if (rsq > cutsq_custom) addme = 0; + + if (addme) + neighptr[n++] = joriginal; } + } else { + #if defined(LMP_SIMD_COMPILER) + #pragma vector aligned + #pragma ivdep + #endif + for (int jj = 0; jj < jnum; jj++) { + const int joriginal = jlist[jj]; + const int j = joriginal & NEIGHMASK; + int addme = 1; + if (j < nlocal) { + if (i > j) addme = 0; + } else { + if (fabs(x[j].z-ztmp) > delta) { + if (x[j].z < ztmp) addme = 0; + } else if (fabs(x[j].y-ytmp) > delta) { + if (x[j].y < ytmp) addme = 0; + } else { + if (x[j].x < xtmp) addme = 0; + } + } - // trim to shorter cutoff + // trim to shorter cutoff - const flt_t delx = xtmp - x[j].x; - const flt_t dely = ytmp - x[j].y; - const flt_t delz = ztmp - x[j].z; - const flt_t rsq = delx * delx + dely * dely + delz * delz; + const flt_t delx = xtmp - x[j].x; + const flt_t dely = ytmp - x[j].y; + const flt_t delz = ztmp - x[j].z; + const flt_t rsq = delx * delx + dely * dely + delz * delz; - if (rsq > cutsq_custom) addme = 0; + if (rsq > cutsq_custom) addme = 0; - if (addme) - neighptr[n++] = joriginal; + if (addme) + neighptr[n++] = joriginal; + } } ilist[ii] = i; @@ -235,7 +274,7 @@ void NPairHalffullNewtonTrimIntel::build_t3(NeighList *list, int *numhalf, void NPairHalffullNewtonTrimIntel::build(NeighList *list) { - if (_fix->three_body_neighbor() == 0) { + if (_fix->three_body_neighbor() == 0 || domain->triclinic) { if (_fix->precision() == FixIntel::PREC_MODE_MIXED) build_t(list, _fix->get_mixed_buffers()); else if (_fix->precision() == FixIntel::PREC_MODE_DOUBLE) diff --git a/src/INTEL/npair_intel.cpp b/src/INTEL/npair_intel.cpp index 600109d7ae..dcfb66e05f 100644 --- a/src/INTEL/npair_intel.cpp +++ b/src/INTEL/npair_intel.cpp @@ -204,6 +204,8 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, } const int special_bound = sb; + const double delta = 0.01 * force->angstrom; + #ifdef _LMP_INTEL_OFFLOAD const int * _noalias const binhead = this->binhead; const int * _noalias const bins = this->bins; @@ -229,7 +231,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, in(ncache_stride,maxnbors,nthreads,maxspecial,nstencil,e_nall,offload) \ in(offload_end,separate_buffers,astart,aend,nlocal,molecular) \ in(ntypes,xperiodic,yperiodic,zperiodic,xprd_half,yprd_half,zprd_half) \ - in(pack_width,special_bound) \ + in(pack_width,special_bound,delta) \ out(overflow:length(5) alloc_if(0) free_if(0)) \ out(timer_compute:length(1) alloc_if(0) free_if(0)) \ signal(tag) @@ -331,7 +333,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, const flt_t ztmp = x[i].z; const int itype = x[i].w; tagint itag; - if (THREE) itag = tag[i]; + if (THREE || (TRI && !FULL)) itag = tag[i]; const int ioffset = ntypes * itype; const int ibin = atombin[i]; @@ -365,7 +367,7 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, ty[u] = x[j].y; tz[u] = x[j].z; tjtype[u] = x[j].w; - if (THREE) ttag[u] = tag[j]; + if (THREE || (TRI && !FULL)) ttag[u] = tag[j]; } if (FULL == 0 && TRI != 1) { @@ -486,12 +488,32 @@ void NPairIntel::bin_newton(const int offload, NeighList *list, // Triclinic if (TRI) { - if (tz[u] < ztmp) addme = 0; - if (tz[u] == ztmp) { - if (ty[u] < ytmp) addme = 0; - if (ty[u] == ytmp) { - if (tx[u] < xtmp) addme = 0; - if (tx[u] == xtmp && j <= i) addme = 0; + if (FULL) { + if (tz[u] < ztmp) addme = 0; + if (tz[u] == ztmp) { + if (ty[u] < ytmp) addme = 0; + if (ty[u] == ytmp) { + if (tx[u] < xtmp) addme = 0; + if (tx[u] == xtmp && j <= i) addme = 0; + } + } + } else { + if (j <= i) addme = 0; + if (j >= nlocal) { + const tagint jtag = ttag[u]; + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) addme = 0; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) addme = 0; + } else { + if (fabs(tz[u]-ztmp) > delta) { + if (tz[u] < ztmp) addme = 0; + } else if (fabs(ty[u]-ytmp) > delta) { + if (ty[u] < ytmp) addme = 0; + } else { + if (tx[u] < xtmp) addme = 0; + } + } } } } From dbab5b69312a2d789f37973ce21873cc2e2757e4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 4 Aug 2023 22:24:48 -0400 Subject: [PATCH 18/41] possible workaround for unit test failure taken from: https://github.com/open-mpi/ompi/issues/9656 --- unittest/formats/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/unittest/formats/CMakeLists.txt b/unittest/formats/CMakeLists.txt index 93ea2f3b32..58c797b6e6 100644 --- a/unittest/formats/CMakeLists.txt +++ b/unittest/formats/CMakeLists.txt @@ -41,6 +41,8 @@ set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${ add_executable(test_file_operations test_file_operations.cpp) target_link_libraries(test_file_operations PRIVATE lammps GTest::GMock) add_test(NAME FileOperations COMMAND test_file_operations) +# try to mitigate possible OpenMPI bug +set_tests_properties(TextFileReader PROPERTIES ENVIRONMENT "OMPI_MCA_sharedfp=\"^sm\"") add_executable(test_dump_atom test_dump_atom.cpp) target_link_libraries(test_dump_atom PRIVATE lammps GTest::GMock) From 462a3935fea188aea0863f300a65fad0b8e35b1b Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Mon, 16 Oct 2023 14:07:42 -0600 Subject: [PATCH 19/41] Port bugfix to Kokkos --- src/KOKKOS/npair_halffull_kokkos.cpp | 17 +++ src/KOKKOS/npair_halffull_kokkos.h | 4 +- src/KOKKOS/npair_kokkos.cpp | 164 ++++++++++++++++++++------- src/KOKKOS/npair_kokkos.h | 6 +- 4 files changed, 143 insertions(+), 48 deletions(-) diff --git a/src/KOKKOS/npair_halffull_kokkos.cpp b/src/KOKKOS/npair_halffull_kokkos.cpp index ec17cec844..bc2549aa8d 100644 --- a/src/KOKKOS/npair_halffull_kokkos.cpp +++ b/src/KOKKOS/npair_halffull_kokkos.cpp @@ -18,6 +18,7 @@ #include "atom_masks.h" #include "atom_vec.h" #include "domain.h" +#include "force.h" #include "neigh_list_kokkos.h" #include @@ -66,6 +67,9 @@ void NPairHalffullKokkos::build(NeighList *list) d_numneigh = k_list->d_numneigh; d_neighbors = k_list->d_neighbors; + delta = 0.01 * force->angstrom; + triclinic = domain->triclinic; + // loop over parent full list copymode = 1; @@ -92,6 +96,11 @@ void NPairHalffullKokkos::operator()(TagNPairHalffullCom } // 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 const int jnum = d_numneigh_full(i); const AtomNeighbors neighbors_i = AtomNeighbors(&d_neighbors(i,0),d_numneigh(i), @@ -103,6 +112,14 @@ void NPairHalffullKokkos::operator()(TagNPairHalffullCom if (NEWTON) { 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/KOKKOS/npair_halffull_kokkos.h b/src/KOKKOS/npair_halffull_kokkos.h index c5a09f0b62..1249a9ce8a 100644 --- a/src/KOKKOS/npair_halffull_kokkos.h +++ b/src/KOKKOS/npair_halffull_kokkos.h @@ -257,8 +257,8 @@ class NPairHalffullKokkos : public NPair { void operator()(TagNPairHalffullCompute, const int&) const; private: - int nlocal; - double cutsq_custom; + int nlocal,triclinic; + double cutsq_custom,delta; typename AT::t_x_array_randomread x; diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index 06567cbeb6..8201ae028b 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -155,6 +155,8 @@ void NPairKokkos::build(NeighList *list_) list->grow(nall); + const double delta = 0.01 * force->angstrom; + NeighborKokkosExecute data(*list, k_cutneighsq.view(), @@ -176,7 +178,7 @@ void NPairKokkos::build(NeighList *list_) atomKK->molecular, nbinx,nbiny,nbinz,mbinx,mbiny,mbinz,mbinxlo,mbinylo,mbinzlo, bininvx,bininvy,bininvz, - exclude, nex_type, + delta, exclude, nex_type, k_ex1_type.view(), k_ex2_type.view(), k_ex_type.view(), @@ -239,7 +241,7 @@ void NPairKokkos::build(NeighList *list_) if (GHOST) { // assumes newton off - NPairKokkosBuildFunctorGhost f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor); + NPairKokkosBuildFunctorGhost f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor); // temporarily disable team policy for ghost due to known bug @@ -261,7 +263,7 @@ void NPairKokkos::build(NeighList *list_) //#endif } else { if (SIZE) { - NPairKokkosBuildFunctorSize f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor); + NPairKokkosBuildFunctorSize f(data,atoms_per_bin * 7 * sizeof(X_FLOAT) * factor); #ifdef LMP_KOKKOS_GPU if (ExecutionSpaceFromDevice::space == Device) { int team_size = atoms_per_bin*factor; @@ -279,7 +281,7 @@ void NPairKokkos::build(NeighList *list_) Kokkos::parallel_for(nall, f); #endif } else { - NPairKokkosBuildFunctor f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor); + NPairKokkosBuildFunctor f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor); #ifdef LMP_KOKKOS_GPU if (ExecutionSpaceFromDevice::space == Device) { int team_size = atoms_per_bin*factor; @@ -414,6 +416,7 @@ void NeighborKokkosExecute:: const X_FLOAT ytmp = x(i, 1); const X_FLOAT ztmp = x(i, 2); const int itype = type(i); + const tagint itag = tag(i); const int ibin = c_atom2bin(i); @@ -484,13 +487,29 @@ void NeighborKokkosExecute:: if (HalfNeigh && !Newton && j <= i) continue; if (!HalfNeigh && j == i) continue; + + // 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 + if (HalfNeigh && Newton && Tri) { - 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) { + const tagint 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; + } } } } @@ -568,8 +587,9 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic size_t sharedsize) const { auto* sharedmem = static_cast(dev.team_shmem().get_shmem(sharedsize)); - /* loop over atoms in i's bin, - */ + + // loop over atoms in i's bin + const int atoms_per_bin = c_bins.extent(1); const int BINS_PER_TEAM = dev.team_size()/atoms_per_bin <1?1:dev.team_size()/atoms_per_bin; const int TEAMS_PER_BIN = atoms_per_bin/dev.team_size()<1?1:atoms_per_bin/dev.team_size(); @@ -579,15 +599,14 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic if (ibin >= mbins) return; - X_FLOAT* other_x = sharedmem + 5*atoms_per_bin*MY_BIN; - int* other_id = (int*) &other_x[4 * atoms_per_bin]; + X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN; + int* other_id = (int*) &other_x[5 * atoms_per_bin]; int bincount_current = c_bincount[ibin]; for (int kk = 0; kk < TEAMS_PER_BIN; kk++) { const int MY_II = dev.team_rank()%atoms_per_bin+kk*dev.team_size(); const int i = MY_II < bincount_current ? c_bins(ibin, MY_II) : -1; - /* if necessary, goto next page and add pages */ int n = 0; @@ -608,6 +627,7 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; + other_x[MY_II + 4 * atoms_per_bin] = itag; } other_id[MY_II] = i; @@ -695,6 +715,7 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); + other_x[MY_II + 4 * atoms_per_bin] = tag(j); } other_id[MY_II] = j; @@ -708,13 +729,29 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic if (HalfNeigh && !Newton && j <= i) continue; if (!HalfNeigh && j == i) continue; + + // 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 + if (HalfNeigh && Newton && Tri) { - 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) { + const tagint jtag = other_x[m + 4 * atoms_per_bin]; + 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; + } } } } @@ -794,6 +831,7 @@ void NeighborKokkosExecute:: const X_FLOAT ytmp = x(i, 1); const X_FLOAT ztmp = x(i, 2); const int itype = type(i); + const tagint itag = tag(i); const typename ArrayTypes::t_int_1d_const_um stencil = d_stencil; @@ -905,6 +943,7 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team size_t sharedsize) const { auto* sharedmem = static_cast(dev.team_shmem().get_shmem(sharedsize)); + // loop over atoms in i's bin const int atoms_per_bin = c_bins.extent(1); @@ -916,8 +955,8 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team if (ibin >= mbins) return; - X_FLOAT* other_x = sharedmem + 5*atoms_per_bin*MY_BIN; - int* other_id = (int*) &other_x[4 * atoms_per_bin]; + X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN; + int* other_id = (int*) &other_x[5 * atoms_per_bin]; int bincount_current = c_bincount[ibin]; @@ -944,6 +983,7 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; + other_x[MY_II + 4 * atoms_per_bin] = itag; } other_id[MY_II] = i; #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) @@ -999,6 +1039,7 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); + other_x[MY_II + 4 * atoms_per_bin] = tag(j); } other_id[MY_II] = j; @@ -1084,6 +1125,7 @@ void NeighborKokkosExecute:: const X_FLOAT ztmp = x(i, 2); const X_FLOAT radi = radius(i); const int itype = type(i); + const tagint itag = tag(i); const int ibin = c_atom2bin(i); @@ -1167,13 +1209,29 @@ void NeighborKokkosExecute:: if (HalfNeigh && !Newton && j <= i) continue; if (!HalfNeigh && j == i) continue; + + // 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 + if (HalfNeigh && Newton && Tri) { - 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) { + const tagint 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; + } } } } @@ -1245,8 +1303,9 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP size_t sharedsize) const { auto* sharedmem = static_cast(dev.team_shmem().get_shmem(sharedsize)); - /* loop over atoms in i's bin, - */ + + // loop over atoms in i's bin + const int atoms_per_bin = c_bins.extent(1); const int BINS_PER_TEAM = dev.team_size()/atoms_per_bin <1?1:dev.team_size()/atoms_per_bin; const int TEAMS_PER_BIN = atoms_per_bin/dev.team_size()<1?1:atoms_per_bin/dev.team_size(); @@ -1256,15 +1315,14 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP if (ibin >= mbins) return; - X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN; - int* other_id = (int*) &other_x[5 * atoms_per_bin]; + X_FLOAT* other_x = sharedmem + 7*atoms_per_bin*MY_BIN; + int* other_id = (int*) &other_x[6 * atoms_per_bin]; int bincount_current = c_bincount[ibin]; for (int kk = 0; kk < TEAMS_PER_BIN; kk++) { const int MY_II = dev.team_rank()%atoms_per_bin+kk*dev.team_size(); const int i = MY_II < bincount_current ? c_bins(ibin, MY_II) : -1; - /* if necessary, goto next page and add pages */ int n = 0; @@ -1288,7 +1346,8 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; - other_x[MY_II + 4 * atoms_per_bin] = radi; + other_x[MY_II + 4 * atoms_per_bin] = itag; + other_x[MY_II + 5 * atoms_per_bin] = radi; } other_id[MY_II] = i; #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) @@ -1323,7 +1382,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP const X_FLOAT dely = ytmp - other_x[m + atoms_per_bin]; const X_FLOAT delz = ztmp - other_x[m + 2 * atoms_per_bin]; const X_FLOAT rsq = delx*delx + dely*dely + delz*delz; - const X_FLOAT radsum = radi + other_x[m + 4 * atoms_per_bin]; + const X_FLOAT radsum = radi + other_x[m + 5 * atoms_per_bin]; const X_FLOAT cutsq = (radsum + skin) * (radsum + skin); if (rsq <= cutsq) { @@ -1380,7 +1439,8 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); - other_x[MY_II + 4 * atoms_per_bin] = radius(j); + other_x[MY_II + 4 * atoms_per_bin] = tag(j); + other_x[MY_II + 5 * atoms_per_bin] = radius(j); } other_id[MY_II] = j; @@ -1394,13 +1454,29 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP if (HalfNeigh && !Newton && j <= i) continue; if (!HalfNeigh && j == i) continue; + + // 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 + if (HalfNeigh && Newton && Tri) { - 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) { + const tagint jtag = other_x[m + 4 * atoms_per_bin]; + 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; + } } } } @@ -1412,7 +1488,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP const X_FLOAT dely = ytmp - other_x[m + atoms_per_bin]; const X_FLOAT delz = ztmp - other_x[m + 2 * atoms_per_bin]; const X_FLOAT rsq = delx*delx + dely*dely + delz*delz; - const X_FLOAT radsum = radi + other_x[m + 4 * atoms_per_bin]; + const X_FLOAT radsum = radi + other_x[m + 5 * atoms_per_bin]; const X_FLOAT cutsq = (radsum + skin) * (radsum + skin); if (rsq <= cutsq) { diff --git a/src/KOKKOS/npair_kokkos.h b/src/KOKKOS/npair_kokkos.h index 4427012926..fe5484a771 100644 --- a/src/KOKKOS/npair_kokkos.h +++ b/src/KOKKOS/npair_kokkos.h @@ -189,6 +189,8 @@ class NeighborKokkosExecute public: NeighListKokkos neigh_list; + const double delta; + // data from Neighbor class const typename AT::t_xfloat_2d_randomread cutneighsq; @@ -282,7 +284,7 @@ class NeighborKokkosExecute const int & _mbinx,const int & _mbiny,const int & _mbinz, const int & _mbinxlo,const int & _mbinylo,const int & _mbinzlo, const X_FLOAT &_bininvx,const X_FLOAT &_bininvy,const X_FLOAT &_bininvz, - const int & _exclude,const int & _nex_type, + const double &_delta,const int & _exclude,const int & _nex_type, const typename AT::t_int_1d_const & _ex1_type, const typename AT::t_int_1d_const & _ex2_type, const typename AT::t_int_2d_const & _ex_type, @@ -301,7 +303,7 @@ class NeighborKokkosExecute const typename ArrayTypes::t_int_scalar _h_resize, const typename AT::t_int_scalar _new_maxneighs, const typename ArrayTypes::t_int_scalar _h_new_maxneighs): - neigh_list(_neigh_list), cutneighsq(_cutneighsq),exclude(_exclude), + neigh_list(_neigh_list), cutneighsq(_cutneighsq),delta(_delta),exclude(_exclude), nex_type(_nex_type),ex1_type(_ex1_type),ex2_type(_ex2_type), ex_type(_ex_type),nex_group(_nex_group), ex1_bit(_ex1_bit),ex2_bit(_ex2_bit), From 4ae0fc83124eff9e400ea43166fc21b4d1ac1ca1 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 16 Oct 2023 15:47:42 -0600 Subject: [PATCH 20/41] Fix GPU compile --- src/KOKKOS/npair_kokkos.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index 8201ae028b..83e60768cd 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -614,6 +614,7 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic X_FLOAT ytmp; X_FLOAT ztmp; int itype; + tagint itag; const int index = (i >= 0 && i < nlocal) ? i : 0; const AtomNeighbors neighbors_i = neigh_transpose ? neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index); @@ -623,6 +624,7 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic ytmp = x(i, 1); ztmp = x(i, 2); itype = type(i); + itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; @@ -970,6 +972,7 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team X_FLOAT ytmp; X_FLOAT ztmp; int itype; + tagint itag; const int index = (i >= 0 && i < nall) ? i : 0; const AtomNeighbors neighbors_i = neigh_transpose ? neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index); @@ -979,6 +982,7 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team ytmp = x(i, 1); ztmp = x(i, 2); itype = type(i); + itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; @@ -1331,6 +1335,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP X_FLOAT ztmp; X_FLOAT radi; int itype; + tagint itag; const int index = (i >= 0 && i < nlocal) ? i : 0; const AtomNeighbors neighbors_i = neigh_transpose ? neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index); @@ -1342,6 +1347,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP ztmp = x(i, 2); radi = radius(i); itype = type(i); + itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; From 3b4fff4164a4dc73433b11f9e7fac65c279a2fb2 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 17 Oct 2023 12:01:01 -0600 Subject: [PATCH 21/41] Need to sync tag, add template param --- src/KOKKOS/npair_halffull_kokkos.cpp | 41 +++-- src/KOKKOS/npair_halffull_kokkos.h | 249 +++++++++++++++++++-------- src/KOKKOS/npair_kokkos.cpp | 37 ++-- 3 files changed, 223 insertions(+), 104 deletions(-) diff --git a/src/KOKKOS/npair_halffull_kokkos.cpp b/src/KOKKOS/npair_halffull_kokkos.cpp index bc2549aa8d..ddd7362c4e 100644 --- a/src/KOKKOS/npair_halffull_kokkos.cpp +++ b/src/KOKKOS/npair_halffull_kokkos.cpp @@ -27,8 +27,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -template -NPairHalffullKokkos::NPairHalffullKokkos(LAMMPS *lmp) : NPair(lmp) { +template +NPairHalffullKokkos::NPairHalffullKokkos(LAMMPS *lmp) : NPair(lmp) { atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; } @@ -42,13 +42,19 @@ NPairHalffullKokkos::NPairHalffullKokkos(LAMMPS *lmp) : if ghost, also store neighbors of ghost atoms & set inum,gnum correctly ------------------------------------------------------------------------- */ -template -void NPairHalffullKokkos::build(NeighList *list) +template +void NPairHalffullKokkos::build(NeighList *list) { if (NEWTON || TRIM) { x = atomKK->k_x.view(); atomKK->sync(execution_space,X_MASK); } + + if (TRI) { + tag = atomKK->k_tag.view(); + atomKK->sync(execution_space,TAG_MASK); + } + nlocal = atom->nlocal; cutsq_custom = cutoff_custom*cutoff_custom; @@ -68,7 +74,6 @@ void NPairHalffullKokkos::build(NeighList *list) d_neighbors = k_list->d_neighbors; delta = 0.01 * force->angstrom; - triclinic = domain->triclinic; // loop over parent full list @@ -82,9 +87,9 @@ void NPairHalffullKokkos::build(NeighList *list) k_list->k_ilist.template modify(); } -template +template KOKKOS_INLINE_FUNCTION -void NPairHalffullKokkos::operator()(TagNPairHalffullCompute, const int &ii) const { +void NPairHalffullKokkos::operator()(TagNPairHalffullCompute, const int &ii) const { int n = 0; const int i = d_ilist_full(ii); @@ -112,7 +117,7 @@ void NPairHalffullKokkos::operator()(TagNPairHalffullCom if (NEWTON) { if (j < nlocal) { if (i > j) continue; - } else if (triclinic) { + } else if (TRI) { if (fabs(x(j,2)-ztmp) > delta) { if (x(j,2) < ztmp) continue; } else if (fabs(x(j,1)-ytmp) > delta) { @@ -158,14 +163,18 @@ void NPairHalffullKokkos::operator()(TagNPairHalffullCom } namespace LAMMPS_NS { -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; #ifdef LMP_KOKKOS_GPU -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; -template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; +template class NPairHalffullKokkos; #endif } diff --git a/src/KOKKOS/npair_halffull_kokkos.h b/src/KOKKOS/npair_halffull_kokkos.h index 1249a9ce8a..7e6c28aaba 100644 --- a/src/KOKKOS/npair_halffull_kokkos.h +++ b/src/KOKKOS/npair_halffull_kokkos.h @@ -16,53 +16,79 @@ // Trim off -// Newton +// Newton, no triclinic -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; NPairStyle(halffull/newton/kk/device, NPairKokkosHalffullNewtonDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE); + NP_ORTHO | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; NPairStyle(halffull/newton/kk/host, NPairKokkosHalffullNewtonHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_KOKKOS_HOST); + NP_ORTHO | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; NPairStyle(halffull/newton/skip/kk/device, NPairKokkosHalffullNewtonDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE); + NP_ORTHO | NP_SKIP | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; NPairStyle(halffull/newton/skip/kk/host, NPairKokkosHalffullNewtonHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_SKIP | NP_KOKKOS_HOST); + +// Newton, triclinic + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriDevice; +NPairStyle(halffull/newton/tri/kk/device, + NPairKokkosHalffullNewtonTriDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriHost; +NPairStyle(halffull/newton/tri/kk/host, + NPairKokkosHalffullNewtonTriHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_KOKKOS_HOST); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriDevice; +NPairStyle(halffull/newton/tri/skip/kk/device, + NPairKokkosHalffullNewtonTriDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriHost; +NPairStyle(halffull/newton/tri/skip/kk/host, + NPairKokkosHalffullNewtonTriHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_HOST); -// Newtoff +// Newtoff (can be triclinic but template param always set to 0) -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; NPairStyle(halffull/newtoff/kk/device, NPairKokkosHalffullNewtoffDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; NPairStyle(halffull/newtoff/kk/host, NPairKokkosHalffullNewtoffHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; NPairStyle(halffull/newtoff/skip/kk/device, NPairKokkosHalffullNewtoffDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; NPairStyle(halffull/newtoff/skip/kk/host, NPairKokkosHalffullNewtoffHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | @@ -70,166 +96,244 @@ NPairStyle(halffull/newtoff/skip/kk/host, //************ Ghost ************** -// Newton +// Newton, no triclinic -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonGhostDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; NPairStyle(halffull/newton/ghost/kk/device, - NPairKokkosHalffullNewtonGhostDevice, + NPairKokkosHalffullNewtonDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE); + NP_ORTHO | NP_GHOST | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; NPairStyle(halffull/newton/ghost/kk/host, NPairKokkosHalffullNewtonHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST); + NP_ORTHO | NP_GHOST | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonGhostDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; NPairStyle(halffull/newton/skip/ghost/kk/device, - NPairKokkosHalffullNewtonGhostDevice, + NPairKokkosHalffullNewtonDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE); + NP_ORTHO | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonHost; NPairStyle(halffull/newton/skip/ghost/kk/host, NPairKokkosHalffullNewtonHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST); + +// Newton, triclinic + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriDevice; +NPairStyle(halffull/newton/tri/ghost/kk/device, + NPairKokkosHalffullNewtonTriDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriHost; +NPairStyle(halffull/newton/tri/ghost/kk/host, + NPairKokkosHalffullNewtonTriHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriDevice; +NPairStyle(halffull/newton/tri/skip/ghost/kk/device, + NPairKokkosHalffullNewtonTriDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriHost; +NPairStyle(halffull/newton/tri/skip/ghost/kk/host, + NPairKokkosHalffullNewtonTriHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST); -// Newtoff +// Newtoff (can be triclinic but template param always set to 0) -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffGhostDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; NPairStyle(halffull/newtoff/ghost/kk/device, - NPairKokkosHalffullNewtoffGhostDevice, + NPairKokkosHalffullNewtoffDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; NPairStyle(halffull/newtoff/ghost/kk/host, NPairKokkosHalffullNewtoffHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffGhostDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffDevice; NPairStyle(halffull/newtoff/skip/ghost/kk/device, - NPairKokkosHalffullNewtoffGhostDevice, + NPairKokkosHalffullNewtoffDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffHost; NPairStyle(halffull/newtoff/skip/ghost/kk/host, NPairKokkosHalffullNewtoffHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_KOKKOS_HOST); - //************ Trim ************** -// Newton +// Newton, no triclinic -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; NPairStyle(halffull/newton/trim/kk/device, NPairKokkosHalffullNewtonTrimDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | - NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE); + NP_ORTHO | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; NPairStyle(halffull/newton/trim/kk/host, NPairKokkosHalffullNewtonTrimHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRIM | NP_KOKKOS_HOST); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; +NPairStyle(halffull/newton/trim/skip/kk/device, + NPairKokkosHalffullNewtonTrimDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; +NPairStyle(halffull/newton/trim/skip/kk/host, + NPairKokkosHalffullNewtonTrimHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); + +// Newton, triclinic + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimDevice; +NPairStyle(halffull/newton/tri/trim/kk/device, + NPairKokkosHalffullNewtonTriTrimDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimHost; +NPairStyle(halffull/newton/tri/trim/kk/host, + NPairKokkosHalffullNewtonTriTrimHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; -NPairStyle(halffull/newton/skip/trim/kk/device, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimDevice; +NPairStyle(halffull/newton/tri/trim/skip/kk/device, NPairKokkosHalffullNewtonTrimDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; -NPairStyle(halffull/newton/skip/trim/kk/host, - NPairKokkosHalffullNewtonTrimHost, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimHost; +NPairStyle(halffull/newton/tri/trim/skip/kk/host, + NPairKokkosHalffullNewtonTriTrimHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); -// Newtoff +// Newtoff (can be triclinic but template param always set to 0) -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; NPairStyle(halffull/newtoff/trim/kk/device, NPairKokkosHalffullNewtoffTrimDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; NPairStyle(halffull/newtoff/trim/kk/host, NPairKokkosHalffullNewtoffTrimHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_TRIM | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; -NPairStyle(halffull/newtoff/skip/trim/kk/device, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; +NPairStyle(halffull/newtoff/trim/skip/kk/device, NPairKokkosHalffullNewtoffTrimDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; -NPairStyle(halffull/newtoff/skip/trim/kk/host, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; +NPairStyle(halffull/newtoff/trim/skip/kk/host, NPairKokkosHalffullNewtoffTrimHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); //************ Ghost ************** -// Newton +// Newton, no triclinic -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonGhostTrimDevice; -NPairStyle(halffull/newton/ghost/trim/kk/device, - NPairKokkosHalffullNewtonGhostTrimDevice, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; +NPairStyle(halffull/newton/tri/trim/ghost/kk/device, + NPairKokkosHalffullNewtonTrimDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; +NPairStyle(halffull/newton/trim/ghost/kk/host, + NPairKokkosHalffullNewtonTrimHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimDevice; +NPairStyle(halffull/newton/trim/skip/ghost/kk/device, + NPairKokkosHalffullNewtonTrimDevice, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; +NPairStyle(halffull/newton/trim/skip/ghost/kk/host, + NPairKokkosHalffullNewtonTrimHost, + NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | + NP_ORTHO | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); + +// Newton, triclinic + +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimDevice; +NPairStyle(halffull/newton/tri/trim/ghost/kk/device, + NPairKokkosHalffullNewtonTriTrimDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; -NPairStyle(halffull/newton/ghost/trim/kk/host, - NPairKokkosHalffullNewtonTrimHost, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimHost; +NPairStyle(halffull/newton/tri/trim/ghost/kk/host, + NPairKokkosHalffullNewtonTriTrimHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonGhostTrimDevice; -NPairStyle(halffull/newton/skip/ghost/trim/kk/device, - NPairKokkosHalffullNewtonGhostTrimDevice, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimDevice; +NPairStyle(halffull/newton/tri/trim/skip/ghost/kk/device, + NPairKokkosHalffullNewtonTriTrimDevice, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTrimHost; -NPairStyle(halffull/newton/skip/ghost/trim/kk/host, - NPairKokkosHalffullNewtonTrimHost, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtonTriTrimHost; +NPairStyle(halffull/newton/tri/trim/skip/ghost/kk/host, + NPairKokkosHalffullNewtonTriTrimHost, NP_HALF_FULL | NP_NEWTON | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); -// Newtoff +// Newtoff (can be triclinic but template param always set to 0) -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffGhostTrimDevice; -NPairStyle(halffull/newtoff/ghost/trim/kk/device, - NPairKokkosHalffullNewtoffGhostTrimDevice, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; +NPairStyle(halffull/newtoff/trim/ghost/kk/device, + NPairKokkosHalffullNewtoffTrimDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; -NPairStyle(halffull/newtoff/ghost/trim/kk/host, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; +NPairStyle(halffull/newtoff/trim/ghost/kk/host, NPairKokkosHalffullNewtoffTrimHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_TRIM | NP_KOKKOS_HOST); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffGhostTrimDevice; -NPairStyle(halffull/newtoff/skip/ghost/trim/kk/device, - NPairKokkosHalffullNewtoffGhostTrimDevice, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimDevice; +NPairStyle(halffull/newtoff/trim/skip/ghost/kk/device, + NPairKokkosHalffullNewtoffTrimDevice, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_DEVICE); -typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; -NPairStyle(halffull/newtoff/skip/ghost/trim/kk/host, +typedef NPairHalffullKokkos NPairKokkosHalffullNewtoffTrimHost; +NPairStyle(halffull/newtoff/trim/skip/ghost/kk/host, NPairKokkosHalffullNewtoffTrimHost, NP_HALF_FULL | NP_NEWTOFF | NP_HALF | NP_NSQ | NP_BIN | NP_MULTI | NP_ORTHO | NP_TRI | NP_GHOST | NP_SKIP | NP_TRIM | NP_KOKKOS_HOST); + // clang-format on #else @@ -244,7 +348,7 @@ namespace LAMMPS_NS { struct TagNPairHalffullCompute{}; -template +template class NPairHalffullKokkos : public NPair { public: typedef DeviceType device_type; @@ -261,6 +365,7 @@ class NPairHalffullKokkos : public NPair { double cutsq_custom,delta; typename AT::t_x_array_randomread x; + typename AT::t_tagint_1d_randomread tag; typename AT::t_neighbors_2d_const d_neighbors_full; typename AT::t_int_1d_const d_ilist_full; diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index 83e60768cd..4e992fb2d7 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -219,6 +219,8 @@ void NPairKokkos::build(NeighList *list_) atomKK->sync(Device,X_MASK|RADIUS_MASK|TYPE_MASK); } + if (HALF && NEWTON && TRI) atomKK->sync(Device,TAG_MASK); + data.special_flag[0] = special_flag[0]; data.special_flag[1] = special_flag[1]; data.special_flag[2] = special_flag[2]; @@ -241,7 +243,7 @@ void NPairKokkos::build(NeighList *list_) if (GHOST) { // assumes newton off - NPairKokkosBuildFunctorGhost f(data,atoms_per_bin * 6 * sizeof(X_FLOAT) * factor); + NPairKokkosBuildFunctorGhost f(data,atoms_per_bin * 5 * sizeof(X_FLOAT) * factor); // temporarily disable team policy for ghost due to known bug @@ -416,7 +418,8 @@ void NeighborKokkosExecute:: const X_FLOAT ytmp = x(i, 1); const X_FLOAT ztmp = x(i, 2); const int itype = type(i); - const tagint itag = tag(i); + tagint itag; + if (HalfNeigh && Newton && Tri) itag = tag(i); const int ibin = c_atom2bin(i); @@ -624,12 +627,14 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic ytmp = x(i, 1); ztmp = x(i, 2); itype = type(i); - itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; - other_x[MY_II + 4 * atoms_per_bin] = itag; + if (HalfNeigh && Newton && Tri) { + itag = tag(i); + other_x[MY_II + 4 * atoms_per_bin] = itag; + } } other_id[MY_II] = i; @@ -717,7 +722,8 @@ void NeighborKokkosExecute::build_ItemGPU(typename Kokkos::TeamPolic other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); - other_x[MY_II + 4 * atoms_per_bin] = tag(j); + if (HalfNeigh && Newton && Tri) + other_x[MY_II + 4 * atoms_per_bin] = tag(j); } other_id[MY_II] = j; @@ -833,7 +839,6 @@ void NeighborKokkosExecute:: const X_FLOAT ytmp = x(i, 1); const X_FLOAT ztmp = x(i, 2); const int itype = type(i); - const tagint itag = tag(i); const typename ArrayTypes::t_int_1d_const_um stencil = d_stencil; @@ -957,8 +962,8 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team if (ibin >= mbins) return; - X_FLOAT* other_x = sharedmem + 6*atoms_per_bin*MY_BIN; - int* other_id = (int*) &other_x[5 * atoms_per_bin]; + X_FLOAT* other_x = sharedmem + 5*atoms_per_bin*MY_BIN; + int* other_id = (int*) &other_x[4 * atoms_per_bin]; int bincount_current = c_bincount[ibin]; @@ -972,7 +977,6 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team X_FLOAT ytmp; X_FLOAT ztmp; int itype; - tagint itag; const int index = (i >= 0 && i < nall) ? i : 0; const AtomNeighbors neighbors_i = neigh_transpose ? neigh_list.get_neighbors_transpose(index) : neigh_list.get_neighbors(index); @@ -982,12 +986,10 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team ytmp = x(i, 1); ztmp = x(i, 2); itype = type(i); - itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; - other_x[MY_II + 4 * atoms_per_bin] = itag; } other_id[MY_II] = i; #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) @@ -1043,7 +1045,6 @@ void NeighborKokkosExecute::build_ItemGhostGPU(typename Kokkos::Team other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); - other_x[MY_II + 4 * atoms_per_bin] = tag(j); } other_id[MY_II] = j; @@ -1129,7 +1130,8 @@ void NeighborKokkosExecute:: const X_FLOAT ztmp = x(i, 2); const X_FLOAT radi = radius(i); const int itype = type(i); - const tagint itag = tag(i); + tagint itag; + if (HalfNeigh && Newton && Tri) itag = tag(i); const int ibin = c_atom2bin(i); @@ -1347,12 +1349,14 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP ztmp = x(i, 2); radi = radius(i); itype = type(i); - itag = tag(i); other_x[MY_II] = xtmp; other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; - other_x[MY_II + 4 * atoms_per_bin] = itag; + if (HalfNeigh && Newton && Tri) { + itag = tag(i); + other_x[MY_II + 4 * atoms_per_bin] = itag; + } other_x[MY_II + 5 * atoms_per_bin] = radi; } other_id[MY_II] = i; @@ -1445,7 +1449,8 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); - other_x[MY_II + 4 * atoms_per_bin] = tag(j); + if (HalfNeigh && Newton && Tri) + other_x[MY_II + 4 * atoms_per_bin] = tag(j); other_x[MY_II + 5 * atoms_per_bin] = radius(j); } From 750957d585ede258dc025a9091d86cde28ce8c60 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 17 Oct 2023 12:27:35 -0600 Subject: [PATCH 22/41] Remove unused var --- src/KOKKOS/npair_halffull_kokkos.cpp | 5 ----- src/KOKKOS/npair_halffull_kokkos.h | 1 - 2 files changed, 6 deletions(-) diff --git a/src/KOKKOS/npair_halffull_kokkos.cpp b/src/KOKKOS/npair_halffull_kokkos.cpp index ddd7362c4e..c8c4d57fc9 100644 --- a/src/KOKKOS/npair_halffull_kokkos.cpp +++ b/src/KOKKOS/npair_halffull_kokkos.cpp @@ -50,11 +50,6 @@ void NPairHalffullKokkos::build(NeighList *list) atomKK->sync(execution_space,X_MASK); } - if (TRI) { - tag = atomKK->k_tag.view(); - atomKK->sync(execution_space,TAG_MASK); - } - nlocal = atom->nlocal; cutsq_custom = cutoff_custom*cutoff_custom; diff --git a/src/KOKKOS/npair_halffull_kokkos.h b/src/KOKKOS/npair_halffull_kokkos.h index 7e6c28aaba..98526c7fee 100644 --- a/src/KOKKOS/npair_halffull_kokkos.h +++ b/src/KOKKOS/npair_halffull_kokkos.h @@ -365,7 +365,6 @@ class NPairHalffullKokkos : public NPair { double cutsq_custom,delta; typename AT::t_x_array_randomread x; - typename AT::t_tagint_1d_randomread tag; typename AT::t_neighbors_2d_const d_neighbors_full; typename AT::t_int_1d_const d_ilist_full; From c051a4cf2df46de7437e7805e4b5598fded4ca27 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 17 Oct 2023 14:18:16 -0600 Subject: [PATCH 23/41] Fix perf regression --- src/KOKKOS/npair_kokkos.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index 4e992fb2d7..f677b3a1bf 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -1353,11 +1353,11 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + atoms_per_bin] = ytmp; other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; + other_x[MY_II + 4 * atoms_per_bin] = radi; if (HalfNeigh && Newton && Tri) { itag = tag(i); - other_x[MY_II + 4 * atoms_per_bin] = itag; + other_x[MY_II + 5 * atoms_per_bin] = itag; } - other_x[MY_II + 5 * atoms_per_bin] = radi; } other_id[MY_II] = i; #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) @@ -1392,7 +1392,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP const X_FLOAT dely = ytmp - other_x[m + atoms_per_bin]; const X_FLOAT delz = ztmp - other_x[m + 2 * atoms_per_bin]; const X_FLOAT rsq = delx*delx + dely*dely + delz*delz; - const X_FLOAT radsum = radi + other_x[m + 5 * atoms_per_bin]; + const X_FLOAT radsum = radi + other_x[m + 4 * atoms_per_bin]; const X_FLOAT cutsq = (radsum + skin) * (radsum + skin); if (rsq <= cutsq) { @@ -1449,9 +1449,9 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + atoms_per_bin] = x(j, 1); other_x[MY_II + 2 * atoms_per_bin] = x(j, 2); other_x[MY_II + 3 * atoms_per_bin] = type(j); + other_x[MY_II + 4 * atoms_per_bin] = radius(j); if (HalfNeigh && Newton && Tri) - other_x[MY_II + 4 * atoms_per_bin] = tag(j); - other_x[MY_II + 5 * atoms_per_bin] = radius(j); + other_x[MY_II + 5 * atoms_per_bin] = tag(j); } other_id[MY_II] = j; @@ -1475,7 +1475,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP if (HalfNeigh && Newton && Tri) { if (j <= i) continue; if (j >= nlocal) { - const tagint jtag = other_x[m + 4 * atoms_per_bin]; + const tagint jtag = other_x[m + 5 * atoms_per_bin]; if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { @@ -1499,7 +1499,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP const X_FLOAT dely = ytmp - other_x[m + atoms_per_bin]; const X_FLOAT delz = ztmp - other_x[m + 2 * atoms_per_bin]; const X_FLOAT rsq = delx*delx + dely*dely + delz*delz; - const X_FLOAT radsum = radi + other_x[m + 5 * atoms_per_bin]; + const X_FLOAT radsum = radi + other_x[m + 4 * atoms_per_bin]; const X_FLOAT cutsq = (radsum + skin) * (radsum + skin); if (rsq <= cutsq) { From 6b184e8079bb33310cd3bf441cff59f82a7305b9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 17 Oct 2023 19:39:31 -0400 Subject: [PATCH 24/41] copy-and-paste bugfix from @stanmoore1 --- src/balance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/balance.cpp b/src/balance.cpp index 3bd083e2b9..6f28081f13 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -473,7 +473,7 @@ void Balance::options(int iarg, int narg, char **arg, int sortflag_default) } iarg += 2+nopt; - } else if (strcmp(arg[iarg+1],"sort") == 0) { + } else if (strcmp(arg[iarg],"sort") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "balance sort", error); sortflag = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; From dbd5f93ed475333e61d5fa2914eb3345ffe27eba Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 17 Oct 2023 19:41:00 -0400 Subject: [PATCH 25/41] whitespace --- src/KOKKOS/npair_halffull_kokkos.h | 2 +- src/KOKKOS/npair_kokkos.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/npair_halffull_kokkos.h b/src/KOKKOS/npair_halffull_kokkos.h index 98526c7fee..3eee19b8c3 100644 --- a/src/KOKKOS/npair_halffull_kokkos.h +++ b/src/KOKKOS/npair_halffull_kokkos.h @@ -16,7 +16,7 @@ // Trim off -// Newton, no triclinic +// Newton, no triclinic typedef NPairHalffullKokkos NPairKokkosHalffullNewtonDevice; NPairStyle(halffull/newton/kk/device, diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index f677b3a1bf..45ec83e90e 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -1354,7 +1354,7 @@ void NeighborKokkosExecute::build_ItemSizeGPU(typename Kokkos::TeamP other_x[MY_II + 2 * atoms_per_bin] = ztmp; other_x[MY_II + 3 * atoms_per_bin] = itype; other_x[MY_II + 4 * atoms_per_bin] = radi; - if (HalfNeigh && Newton && Tri) { + if (HalfNeigh && Newton && Tri) { itag = tag(i); other_x[MY_II + 5 * atoms_per_bin] = itag; } From ea69d77b790f28632d6ae345eca6b6a143f010bf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 17 Oct 2023 20:02:01 -0400 Subject: [PATCH 26/41] fix issues with Qt library detection --- tools/lammps-gui/CMakeLists.txt | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tools/lammps-gui/CMakeLists.txt b/tools/lammps-gui/CMakeLists.txt index edfeeb1128..e83db05fdd 100644 --- a/tools/lammps-gui/CMakeLists.txt +++ b/tools/lammps-gui/CMakeLists.txt @@ -110,14 +110,15 @@ endif() # we require Qt 5 and at least version 5.12 at that. if(NOT LAMMPS_GUI_USE_QT5) - find_package(Qt6 6.2 COMPONENTS Widgets Charts) + find_package(Qt6 6.2 QUIET COMPONENTS Widgets Charts) endif() if(NOT Qt6_FOUND) find_package(Qt5 5.12 REQUIRED COMPONENTS Widgets Charts) - set(QT_VERSION_MAJOR "5") + set(QT_VERSION_MAJOR 5) else() - set(QT_VERSION_MAJOR "6") + set(QT_VERSION_MAJOR 6) endif() +message(STATUS "Using Qt version ${Qt${QT_VERSION_MAJOR}_VERSION} for LAMMPS GUI") set(PROJECT_SOURCES main.cpp @@ -188,7 +189,7 @@ else() endif() target_include_directories(lammps-gui PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) target_compile_definitions(lammps-gui PRIVATE LAMMPS_GUI_VERSION="${PROJECT_VERSION}") -target_link_libraries(lammps-gui PRIVATE Qt${QT_VERSION_MAJOR}::Widgets Qt${VERSION_MAJOR}::Charts) +target_link_libraries(lammps-gui PRIVATE Qt${QT_VERSION_MAJOR}::Widgets Qt${QT_VERSION_MAJOR}::Charts) if(BUILD_OMP) find_package(OpenMP COMPONENTS CXX REQUIRED) target_link_libraries(lammps-gui PRIVATE OpenMP::OpenMP_CXX) From 4c980eec91457763153b530a45095aa066cc96d1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 17 Oct 2023 20:52:13 -0400 Subject: [PATCH 27/41] correct table formatting --- doc/src/variable.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/variable.rst b/doc/src/variable.rst index f1a316da1f..3407b48156 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -1183,7 +1183,7 @@ table: +--------+------------+------------------------------------------+ | vector | c_ID | global vector | | vector | c_ID[I] | column of global array | ----------+------------+------------------------------------------+ ++--------+------------+------------------------------------------+ | atom | c_ID | per-atom vector | | atom | c_ID[I] | column of per-atom array | +--------+------------+------------------------------------------+ @@ -1247,7 +1247,7 @@ and atom-style variables are listed in the following table: +--------+------------+------------------------------------------+ | vector | f_ID | global vector | | vector | f_ID[I] | column of global array | ----------+------------+------------------------------------------+ ++--------+------------+------------------------------------------+ | atom | f_ID | per-atom vector | | atom | f_ID[I] | column of per-atom array | +--------+------------+------------------------------------------+ From 302e3be66978f81205d63a16b60d4eb5d9b5a34a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 18 Oct 2023 07:51:25 -0400 Subject: [PATCH 28/41] make sure itag is initialized --- src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp index 7faa210107..1c6d025fab 100644 --- a/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp @@ -87,6 +87,7 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) n = 0; neighptr = ipage.vget(); + itag = tag[i]; itype = type[i]; xtmp = x[i][0]; ytmp = x[i][1]; From 5cfd8b3c63d12051c988b1fe30c289d1d35a9bc5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 18 Oct 2023 07:56:50 -0400 Subject: [PATCH 29/41] silence coverity scan warning --- src/fix_property_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 3a53110839..9613523059 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -46,6 +46,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) : rmass_flag = 0; temperature_flag = 0; heatflow_flag = 0; + nmax_old = 0; nvalue = 0; values_peratom = 0; @@ -212,7 +213,6 @@ void FixPropertyAtom::post_constructor() { // perform initial allocation of atom-based array - nmax_old = 0; grow_arrays(atom->nmax); } From a5f61c5d4405a8af7c9c6962936a21c0622864b2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 18 Oct 2023 07:59:22 -0400 Subject: [PATCH 30/41] address bugs and issues flagged by static code analysis with coverity scan --- src/fix_press_langevin.cpp | 30 +++++++++++++++++------------- src/fix_press_langevin.h | 7 +++---- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/fix_press_langevin.cpp b/src/fix_press_langevin.cpp index d8b90d8b49..2f6e765cd5 100644 --- a/src/fix_press_langevin.cpp +++ b/src/fix_press_langevin.cpp @@ -46,7 +46,8 @@ enum { ISO, ANISO, TRICLINIC }; /* ---------------------------------------------------------------------- */ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), id_press(nullptr), pflag(0), random(nullptr), irregular(nullptr) + Fix(lmp, narg, arg), id_temp(nullptr), id_press(nullptr), temperature(nullptr), + pressure(nullptr), irregular(nullptr), random(nullptr) { if (narg < 5) utils::missing_cmd_args(FLERR, "fix press/langevin", error); @@ -62,6 +63,9 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : allremap = 1; pre_exchange_flag = 0; flipflag = 1; + seed = 111111; + pflag = 0; + kspace_flag = 0; p_ltime = 0.0; @@ -239,7 +243,7 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : t_start = utils::numeric(FLERR, arg[iarg + 1], false, lmp); t_stop = utils::numeric(FLERR, arg[iarg + 2], false, lmp); seed = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (seed <= 0.0) error->all(FLERR, "Fix press/langevin temp seed must be > 0"); + if (seed <= 0) error->all(FLERR, "Fix press/langevin temp seed must be > 0"); iarg += 4; } @@ -349,7 +353,7 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : // Kinetic contribution will be added by the fix style id_press = utils::strdup(std::string(id) + "_press"); - modify->add_compute(fmt::format("{} all pressure NULL virial", id_press)); + pressure = modify->add_compute(fmt::format("{} all pressure NULL virial", id_press)); pflag = 1; // p_fric is alpha coeff from GJF @@ -482,7 +486,7 @@ void FixPressLangevin::initial_integrate(int /* vflag */) if (delta != 0.0) delta /= update->endstep - update->beginstep; t_target = t_start + delta * (t_stop - t_start); - couple_beta(t_target); + couple_beta(); dt = update->dt; @@ -492,11 +496,12 @@ void FixPressLangevin::initial_integrate(int /* vflag */) displacement = dt * p_deriv[i] * gjfb[i]; displacement += 0.5 * dt * dt * f_piston[i] * gjfb[i] / p_mass[i]; displacement += 0.5 * dt * fran[i] * gjfb[i] / p_mass[i]; - dl = domain->boxhi[i] - domain->boxlo[i]; - if (i < 3) + if (i < 3) { + dl = domain->boxhi[i] - domain->boxlo[i]; dilation[i] = (dl + displacement) / dl; - else + } else { dilation[i] = displacement; + } } } } @@ -527,7 +532,7 @@ void FixPressLangevin::post_force(int /*vflag*/) } couple_pressure(); - couple_kinetic(t_target); + couple_kinetic(); for (int i = 0; i < 6; i++) { if (p_flag[i]) { @@ -594,10 +599,9 @@ void FixPressLangevin::couple_pressure() } /* ---------------------------------------------------------------------- */ -void FixPressLangevin::couple_kinetic(double t_target) +void FixPressLangevin::couple_kinetic() { double pk, volume; - nktv2p = force->nktv2p; // kinetic part @@ -607,7 +611,7 @@ void FixPressLangevin::couple_kinetic(double t_target) volume = domain->xprd * domain->yprd; pk = atom->natoms * force->boltz * t_target / volume; - pk *= nktv2p; + pk *= force->nktv2p; p_current[0] += pk; p_current[1] += pk; @@ -616,7 +620,7 @@ void FixPressLangevin::couple_kinetic(double t_target) /* ---------------------------------------------------------------------- */ -void FixPressLangevin::couple_beta(double t_target) +void FixPressLangevin::couple_beta() { double gamma[6]; int me = comm->me; @@ -812,7 +816,7 @@ int FixPressLangevin::modify_param(int narg, char **arg) id_press = utils::strdup(arg[1]); pressure = modify->get_compute_by_id(arg[1]); - if (pressure) error->all(FLERR, "Could not find fix_modify pressure compute ID: {}", arg[1]); + if (!pressure) error->all(FLERR, "Could not find fix_modify pressure compute ID: {}", arg[1]); if (pressure->pressflag == 0) error->all(FLERR, "Fix_modify pressure compute {} does not compute pressure", arg[1]); return 2; diff --git a/src/fix_press_langevin.h b/src/fix_press_langevin.h index fc9f39c434..868993b1f4 100644 --- a/src/fix_press_langevin.h +++ b/src/fix_press_langevin.h @@ -40,10 +40,9 @@ class FixPressLangevin : public Fix { int modify_param(int, char **) override; protected: - int dimension, which; + int dimension; int pstyle, pcouple, allremap; int p_flag[6]; // 1 if control P on this dim, 0 if not - double nktv2p; double t_start, t_stop, t_target; double p_fric[6], p_ltime; // Friction and Langevin charac. time double p_alpha[6]; @@ -68,8 +67,8 @@ class FixPressLangevin : public Fix { int seed; void couple_pressure(); - void couple_kinetic(double); - void couple_beta(double); + void couple_kinetic(); + void couple_beta(); void remap(); }; From 54ff01d86daa134e48feeef503da52b750cf9ffc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 18 Oct 2023 08:01:41 -0400 Subject: [PATCH 31/41] assign code owner to fix press/langevin --- .github/CODEOWNERS | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index f99a336dbb..9b316fbeb9 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -135,6 +135,7 @@ src/timer.* @akohlmey src/utils.* @akohlmey @rbberger src/verlet.* @sjplimp @stanmoore1 src/math_eigen_impl.h @jewettaij +src/fix_press_langevin.* @Bibobu # tools tools/coding_standard/* @akohlmey @rbberger From 45d2a91c6289466f3940571e5845aad527b912a8 Mon Sep 17 00:00:00 2001 From: Maria-Lesniewski Date: Wed, 18 Oct 2023 13:13:37 -0400 Subject: [PATCH 32/41] Barostat fix - see lammps PR 879 and 942 --- src/BOCS/fix_bocs.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/BOCS/fix_bocs.cpp b/src/BOCS/fix_bocs.cpp index d17884855a..52ca948657 100644 --- a/src/BOCS/fix_bocs.cpp +++ b/src/BOCS/fix_bocs.cpp @@ -1024,7 +1024,10 @@ void FixBocs::final_integrate() if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); - else pressure->compute_vector(); + else { + temperature->compute_vector(); + pressure->compute_vector(); + } couple(); pressure->addstep(update->ntimestep+1); } @@ -1961,7 +1964,7 @@ void FixBocs::nhc_press_integrate() int ich,i,pdof; double expfac,factor_etap,kecurrent; double kt = boltz * t_target; - + double lkt_press; // Update masses, to preserve initial freq, if flag set if (omega_mass_flag) { @@ -2006,7 +2009,8 @@ void FixBocs::nhc_press_integrate() } } - double lkt_press = pdof * kt; + if (pstyle == ISO) lkt_press = kt; + else lkt_press = pdof * kt; etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0]; double ncfac = 1.0/nc_pchain; From 27e0a7184954c9b5a2d9c5b3b68c700022c7c050 Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Thu, 19 Oct 2023 07:44:44 -0600 Subject: [PATCH 33/41] whitespace --- src/BOCS/fix_bocs.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/BOCS/fix_bocs.cpp b/src/BOCS/fix_bocs.cpp index 52ca948657..17bb1af002 100644 --- a/src/BOCS/fix_bocs.cpp +++ b/src/BOCS/fix_bocs.cpp @@ -1025,8 +1025,8 @@ void FixBocs::final_integrate() if (pstat_flag) { if (pstyle == ISO) pressure->compute_scalar(); else { - temperature->compute_vector(); - pressure->compute_vector(); + temperature->compute_vector(); + pressure->compute_vector(); } couple(); pressure->addstep(update->ntimestep+1); @@ -1965,6 +1965,7 @@ void FixBocs::nhc_press_integrate() double expfac,factor_etap,kecurrent; double kt = boltz * t_target; double lkt_press; + // Update masses, to preserve initial freq, if flag set if (omega_mass_flag) { From 0fc9b194dc36176e2eb2dd472cc50087c9e7e1da Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 19 Oct 2023 10:01:10 -0400 Subject: [PATCH 34/41] quote strings with special characters in keyword lists --- src/EXTRA-DUMP/dump_yaml.cpp | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-DUMP/dump_yaml.cpp b/src/EXTRA-DUMP/dump_yaml.cpp index 3ca5c59edf..6c21c24f77 100644 --- a/src/EXTRA-DUMP/dump_yaml.cpp +++ b/src/EXTRA-DUMP/dump_yaml.cpp @@ -24,6 +24,8 @@ using namespace LAMMPS_NS; +static constexpr char special_chars[] = "{}[],&:*#?|-<>=!%@\\"; + /* ---------------------------------------------------------------------- */ DumpYAML::DumpYAML(class LAMMPS *_lmp, int narg, char **args) : DumpCustom(_lmp, narg, args), thermo(false) @@ -67,7 +69,12 @@ void DumpYAML::write_header(bigint ndump) const auto &fields = th->get_fields(); thermo_data += "thermo:\n - keywords: [ "; - for (int i = 0; i < nfield; ++i) thermo_data += fmt::format("{}, ", keywords[i]); + for (int i = 0; i < nfield; ++i) { + if (keywords[i].find_first_of(special_chars) == std::string::npos) + thermo_data += fmt::format("{}, ", keywords[i]); + else + thermo_data += fmt::format("'{}', ", keywords[i]); + } thermo_data += "]\n - data: [ "; for (int i = 0; i < nfield; ++i) { @@ -107,7 +114,12 @@ void DumpYAML::write_header(bigint ndump) if (domain->triclinic) fmt::print(fp, " - [ {}, {}, {} ]\n", boxxy, boxxz, boxyz); fmt::print(fp, "keywords: [ "); - for (const auto &item : utils::split_words(columns)) fmt::print(fp, "{}, ", item); + for (const auto &item : utils::split_words(columns)) { + if (item.find_first_of(special_chars) == std::string::npos) + fmt::print(fp, "{}, ", item); + else + fmt::print(fp, "'{}', ", item); + } fputs(" ]\ndata:\n", fp); } else // reset so that the remainder of the output is not multi-proc filewriter = 0; From 62ea874d9df5727a97d90bf922291b0eebd103fc Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Thu, 19 Oct 2023 08:02:11 -0600 Subject: [PATCH 35/41] Add more error checks to Kokkos minimize --- src/KOKKOS/min_kokkos.cpp | 3 +++ src/min.cpp | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/KOKKOS/min_kokkos.cpp b/src/KOKKOS/min_kokkos.cpp index 4e1c3967ff..bbb9a0bd6e 100644 --- a/src/KOKKOS/min_kokkos.cpp +++ b/src/KOKKOS/min_kokkos.cpp @@ -59,6 +59,9 @@ void MinKokkos::init() { Min::init(); + if (!fix_minimize->kokkosable) + error->all(FLERR,"KOKKOS package requires fix minimize/kk"); + fix_minimize_kk = (FixMinimizeKokkos*) fix_minimize; } diff --git a/src/min.cpp b/src/min.cpp index 5a469a788b..acc7d17654 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -215,6 +215,9 @@ void Min::setup(int flag) } update->setupflag = 1; + if (lmp->kokkos) + error->all(FLERR,"KOKKOS package requires Kokkos-enabled min_style"); + // setup extra global dof due to fixes // cannot be done in init() b/c update init() is before modify init() From 2401cba84a5e9b07f50fe3a71a4cb6ec1048b025 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Thu, 19 Oct 2023 12:15:01 -0700 Subject: [PATCH 36/41] Fix bug in Kokkos SNAP on GPUs --- src/KOKKOS/pair_snap_kokkos_impl.h | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/KOKKOS/pair_snap_kokkos_impl.h b/src/KOKKOS/pair_snap_kokkos_impl.h index 91f432dbaf..68303f9b33 100644 --- a/src/KOKKOS/pair_snap_kokkos_impl.h +++ b/src/KOKKOS/pair_snap_kokkos_impl.h @@ -63,10 +63,6 @@ PairSNAPKokkos::PairSNAPKokkos(LAMMPS *lmp datamask_read = EMPTY_MASK; datamask_modify = EMPTY_MASK; - k_cutsq = tdual_fparams("PairSNAPKokkos::cutsq",atom->ntypes+1,atom->ntypes+1); - auto d_cutsq = k_cutsq.template view(); - rnd_cutsq = d_cutsq; - host_flag = (execution_space == Host); } @@ -101,6 +97,8 @@ void PairSNAPKokkos::init_style() if (force->newton_pair == 0) error->all(FLERR,"Pair style SNAP requires newton pair on"); + + // neighbor list request for KOKKOS neighflag = lmp->kokkos->neighflag; @@ -546,6 +544,9 @@ void PairSNAPKokkos::allocate() int n = atom->ntypes; MemKK::realloc_kokkos(d_map,"PairSNAPKokkos::map",n+1); + + MemKK::realloc_kokkos(k_cutsq,"PairSNAPKokkos::cutsq",n+1,n+1); + rnd_cutsq = k_cutsq.template view(); } From 1de15c38bf9375a77768c7b4530fc999a2ebd8a8 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Thu, 19 Oct 2023 12:16:26 -0700 Subject: [PATCH 37/41] whitespace --- src/KOKKOS/pair_snap_kokkos_impl.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/KOKKOS/pair_snap_kokkos_impl.h b/src/KOKKOS/pair_snap_kokkos_impl.h index 68303f9b33..7b9fda60db 100644 --- a/src/KOKKOS/pair_snap_kokkos_impl.h +++ b/src/KOKKOS/pair_snap_kokkos_impl.h @@ -97,8 +97,6 @@ void PairSNAPKokkos::init_style() if (force->newton_pair == 0) error->all(FLERR,"Pair style SNAP requires newton pair on"); - - // neighbor list request for KOKKOS neighflag = lmp->kokkos->neighflag; From 0f11a9dd70b8e839110b8c65554c3f0a91dc09ba Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 19 Oct 2023 23:21:18 -0400 Subject: [PATCH 38/41] require version newer than 2 Aug 2023 for LAMMPS GUI 1.5.9 --- tools/lammps-gui/lammpsgui.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index 11f2554b55..37b8aebf88 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -1432,12 +1432,12 @@ void LammpsGui::start_lammps() lammps.open(narg, args); lammpsstatus->show(); - // must have at least 2 August 2023 version of LAMMPS + // must have a version newer than the 2 August 2023 release of LAMMPS // TODO: must update this check before next feature release - if (lammps.version() < 20230802) { + if (lammps.version() <= 20230802) { QMessageBox::critical(this, "Incompatible LAMMPS Version", "LAMMPS-GUI version " LAMMPS_GUI_VERSION " requires\n" - "LAMMPS version 2 August 2023 or later"); + "a LAMMPS version more recent than 2 August 2023"); exit(1); } From 291defb45378eb759470347c97610335a26c9efb Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 20 Oct 2023 12:57:19 -0600 Subject: [PATCH 39/41] update refs on fix srd doc page --- doc/src/fix_srd.rst | 51 +++++++++++++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 18 deletions(-) diff --git a/doc/src/fix_srd.rst b/doc/src/fix_srd.rst index 8bfbcf2387..c15a566563 100644 --- a/doc/src/fix_srd.rst +++ b/doc/src/fix_srd.rst @@ -61,25 +61,30 @@ Description Treat a group of particles as stochastic rotation dynamics (SRD) particles that serve as a background solvent when interacting with big (colloidal) particles in groupbig-ID. The SRD formalism is described -in :ref:`(Hecht) `. The key idea behind using SRD particles as a -cheap coarse-grained solvent is that SRD particles do not interact -with each other, but only with the solute particles, which in LAMMPS -can be spheroids, ellipsoids, or line segments, or triangles, or rigid -bodies containing multiple spheroids or ellipsoids or line segments -or triangles. The collision and rotation properties of the model -imbue the SRD particles with fluid-like properties, including an -effective viscosity. Thus simulations with large solute particles can -be run more quickly, to measure solute properties like diffusivity -and viscosity in a background fluid. The usual LAMMPS fixes for such -simulations, such as :doc:`fix deform `, -:doc:`fix viscosity `, and :doc:`fix nvt/sllod `, -can be used in conjunction with the SRD model. +in :ref:`(Hecht) `. The same methodology is also called +multi-particle collision dynamics (MPCD) in the literature. -For more details on how the SRD model is implemented in LAMMPS, -:ref:`(Petersen) ` describes the implementation and usage of -pure SRD fluids. See the ``examples/srd`` directory for sample input -scripts using SRD particles for that and for mixture systems (solute -particles in an SRD fluid). +The key idea behind using SRD particles as a cheap coarse-grained +solvent is that SRD particles do not interact with each other, but +only with the solute particles, which in LAMMPS can be spheroids, +ellipsoids, or line segments, or triangles, or rigid bodies containing +multiple spheroids or ellipsoids or line segments or triangles. The +collision and rotation properties of the model imbue the SRD particles +with fluid-like properties, including an effective viscosity. Thus +simulations with large solute particles can be run more quickly, to +measure solute properties like diffusivity and viscosity in a +background fluid. The usual LAMMPS fixes for such simulations, such +as :doc:`fix deform `, :doc:`fix viscosity +`, and :doc:`fix nvt/sllod `, can be +used in conjunction with the SRD model. + +These 3 papers give more details on how the SRD model is implemented +in LAMMPS. :ref:`(Petersen) ` describes pure SRD fluid +systems. :ref:`(Bolintineanu1) ` describes models +where pure SRD fluids :ref:interact with boundary walls. +:ref:`(Bolintineanu2) :ref:` describes mixture models +where large colloidal particles are solvated by an SRD fluid. See the +``examples/srd`` :ref:directory for sample input scripts. This fix does two things: @@ -404,3 +409,13 @@ no, and rescale = yes. **(Petersen)** Petersen, Lechman, Plimpton, Grest, in' t Veld, Schunk, J Chem Phys, 132, 174106 (2010). + +.. _Bolintineanu1: + +**(Bolintineanu1)** +Bolintineanu, Lechman, Plimpton, Grest, Phys Rev E, 86, 066703 (2012). + +.. _Bolintineanu2: + +**(Bolintineanu2)** Bolintineanu, Grest, Lechman, Pierce, Plimpton, +Schunk, Comp Particle Mechanics, 1, 321-356 (2014). From 654a410b8c90840003c0dbaed27ae8672cccff27 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 20 Oct 2023 12:59:36 -0600 Subject: [PATCH 40/41] update refs on fix srd doc page --- doc/src/fix_srd.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_srd.rst b/doc/src/fix_srd.rst index c15a566563..59044a6e3b 100644 --- a/doc/src/fix_srd.rst +++ b/doc/src/fix_srd.rst @@ -82,8 +82,8 @@ These 3 papers give more details on how the SRD model is implemented in LAMMPS. :ref:`(Petersen) ` describes pure SRD fluid systems. :ref:`(Bolintineanu1) ` describes models where pure SRD fluids :ref:interact with boundary walls. -:ref:`(Bolintineanu2) :ref:` describes mixture models -where large colloidal particles are solvated by an SRD fluid. See the +:ref:`(Bolintineanu2) ` describes mixture models where +large colloidal particles are solvated by an SRD fluid. See the ``examples/srd`` :ref:directory for sample input scripts. This fix does two things: From fd05acec2fc0fc8756a2e8564a45a276f8f2e3c2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Oct 2023 20:54:05 -0400 Subject: [PATCH 41/41] fix spelling --- doc/src/compute.rst | 2 +- doc/src/compute_voronoi_atom.rst | 2 +- doc/src/fix.rst | 2 +- doc/src/fix_ave_histo.rst | 4 ++-- doc/src/thermo_style.rst | 2 +- doc/src/variable.rst | 4 ++-- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/doc/src/compute.rst b/doc/src/compute.rst index cc26d9acc9..6737203618 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -74,7 +74,7 @@ Global, per-atom, local, and per-grid quantities can also be of three for each atom, each local entity, or each grid cell. Note that a single compute can produce any combination of global, -per-atom, local, or per-grid values. Likewise it can prouduce any +per-atom, local, or per-grid values. Likewise it can produce any combination of scalar, vector, or array output for each style. The exception is that for per-atom, local, and per-grid output, either a vector or array can be produced, but not both. The doc page for each diff --git a/doc/src/compute_voronoi_atom.rst b/doc/src/compute_voronoi_atom.rst index 37e5386341..9607401ccd 100644 --- a/doc/src/compute_voronoi_atom.rst +++ b/doc/src/compute_voronoi_atom.rst @@ -232,4 +232,4 @@ Related commands Default """"""" -The default for the neighobrs keyword is no. +The default for the neighbors keyword is no. diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 3dd7e224e7..0889fe281f 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -106,7 +106,7 @@ Global, per-atom, local, and per-grid quantities can also be of three for each atom, each local entity, or each grid cell. Note that a single fix can produce any combination of global, -per-atom, local, or per-grid values. Likewise it can prouduce any +per-atom, local, or per-grid values. Likewise it can produce any combination of scalar, vector, or array output for each style. The exception is that for per-atom, local, and per-grid output, either a vector or array can be produced, but not both. The doc page for each diff --git a/doc/src/fix_ave_histo.rst b/doc/src/fix_ave_histo.rst index 31e5476f9e..9699e4238c 100644 --- a/doc/src/fix_ave_histo.rst +++ b/doc/src/fix_ave_histo.rst @@ -106,7 +106,7 @@ attributes are per-atom vector values. See the page for individual generate. Note that a compute or fix can produce multiple kinds of data (global, -per-atom, local). If LAMMPS cannot unambiguosly determine which kind +per-atom, local). If LAMMPS cannot unambiguously determine which kind of data to use, the optional *kind* keyword discussed below can force the desired disambiguation. @@ -263,7 +263,7 @@ keyword is set to *vector*, then all input values must be global or per-atom or local vectors, or columns of global or per-atom or local arrays. -The *kind* keyword only needs to be used if any of the specfied input +The *kind* keyword only needs to be used if any of the specified input computes or fixes produce more than one kind of output (global, per-atom, local). If not, LAMMPS will determine the kind of data all the inputs produce and verify it is all the same kind. If not, an diff --git a/doc/src/thermo_style.rst b/doc/src/thermo_style.rst index c3c607a479..89a2c0b740 100644 --- a/doc/src/thermo_style.rst +++ b/doc/src/thermo_style.rst @@ -442,7 +442,7 @@ equal-style and vector-style variables can be referenced; the latter requires a bracketed term to specify the Ith element of the vector calculated by the variable. However, an equal-style variable can use an atom-style variable in its formula indexed by the ID of an -individual atom. This is a way to output a speciic atom's per-atom +individual atom. This is a way to output a specific atom's per-atom coordinates or other per-atom properties in thermo output. See the :doc:`variable ` command for details. Note that variables of style *equal* and *vector* and *atom* define a formula which can diff --git a/doc/src/variable.rst b/doc/src/variable.rst index 3407b48156..92a78ee3c1 100644 --- a/doc/src/variable.rst +++ b/doc/src/variable.rst @@ -1167,7 +1167,7 @@ variables), or global vectors of values. The latter can also be a column of a global array. Atom-style variables can use scalar values (same as for equal-style -varaibles), or per-atom vectors of values. The latter can also be a +variables), or per-atom vectors of values. The latter can also be a column of a per-atom array. The various allowed compute references in the variable formulas for @@ -1232,7 +1232,7 @@ variables), or global vectors of values. The latter can also be a column of a global array. Atom-style variables can use scalar values (same as for equal-style -varaibles), or per-atom vectors of values. The latter can also be a +variables), or per-atom vectors of values. The latter can also be a column of a per-atom array. The allowed fix references in variable formulas for equal-, vector-,