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; + } + } } } }