Skip to content

Commit

Permalink
Merge pull request lammps#3846 from lammps/triclinic-neighbor-bug
Browse files Browse the repository at this point in the history 
Fix occasional bug in neighbor list build for triclinic geometries due to round off
  • Loading branch information
@akohlmey
akohlmey authored Oct 17, 2023
2 parents 0fe6218 + c051a4c commit f080133
Show file tree
Hide file tree
Showing 43 changed files with 1,203 additions and 519 deletions.
3 changes: 2 additions & 1 deletion doc/src/Howto_triclinic.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,8 @@ is created, e.g. by the :doc:`create_box <create_box>` or
:doc:`read_data <read_data>` or :doc:`read_restart <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 <thermo_style>` command.
parameters, as well as lx,ly,lz, can be output via the
:doc:`thermo_style custom <thermo_style>` command.

LAMMPS also allows simulations to be performed in triclinic
(non-orthogonal) simulation boxes shaped as a parallelepiped with
Expand Down
5 changes: 5 additions & 0 deletions doc/src/atom_modify.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,11 @@ switch. This is described on the :doc:`Build_settings <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 <Howto_triclinic>` 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
Expand Down
2 changes: 2 additions & 0 deletions src/INTEL/fix_intel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "fix_intel.h"

#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "neighbor.h"
Expand Down Expand Up @@ -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];
Expand Down
66 changes: 48 additions & 18 deletions src/INTEL/npair_halffull_newton_intel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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)
Expand Down
87 changes: 63 additions & 24 deletions src/INTEL/npair_halffull_newton_trim_intel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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)
Expand Down
40 changes: 31 additions & 9 deletions src/INTEL/npair_intel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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)
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
}
}
}
}
}
Expand Down
Loading

0 comments on commit f080133

Please sign in to comment.