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

Adding triclinic patch (except for intel)

Browse Source
This commit is contained in:
jtclemm
2023-10-22 21:52:43 -06:00
parent b456beb62f 49f08e6231
commit 1571c96480
3033 changed files with 125592 additions and 129201 deletions
Download Patch File Download Diff File Expand all files Collapse all files

109
src/npair_bin.cpp
View File

@ -1,5 +1,4 @@
// clang-format off
/* -*- c++ -*- ----------------------------------------------------------
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
@ -13,14 +12,16 @@
------------------------------------------------------------------------- */
#include "npair_bin.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "neighbor.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"
using namespace LAMMPS_NS;
using namespace NeighConst;
@ -48,12 +49,13 @@ NPairBin<HALF, NEWTON, TRI, SIZE>::NPairBin(LAMMPS *lmp) : NPair(lmp) {}
template<int HALF, int NEWTON, int TRI, int SIZE>
void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
{
int i,j,jh,k,n,itype,jtype,ibin,bin_start,which,imol,iatom,moltemplate;
tagint tagprev;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
double radsum,cut,cutsq;
int i, j, jh, k, n, itype, jtype, ibin, bin_start, which, imol, iatom, moltemplate;
tagint itag, jtag, tagprev;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq, radsum,cut,cutsq;
int *neighptr;
const double delta = 0.01 * force->angstrom;
double **x = atom->x;
double *radius = atom->radius;
int *type = atom->type;
@ -68,8 +70,10 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
int *molindex = atom->molindex;
int *molatom = atom->molatom;
Molecule **onemols = atom->avec->onemols;
if (molecular == Atom::TEMPLATE) moltemplate = 1;
else moltemplate = 0;
if (molecular == Atom::TEMPLATE)
moltemplate = 1;
else
moltemplate = 0;
int history = list->history;
int mask_history = 1 << HISTBITS;
@ -86,6 +90,7 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
n = 0;
neighptr = ipage->vget();
itag = tag[i];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
@ -99,7 +104,7 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
ibin = atom2bin[i];
for (k = 0; k < nstencil; k++) {
bin_start = binhead[ibin+stencil[k]];
bin_start = binhead[ibin + stencil[k]];
if (stencil[k] == 0) {
if (HALF && NEWTON && (!TRI)) {
// Half neighbor list, newton on, orthonormal
@ -120,17 +125,26 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
// stores own/ghost pairs on both procs
if (j <= i) continue;
} else if (TRI) {
// Half neighbor list, newton on, triclinic
// 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
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;
// 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 (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;
}
}
}
} else {
@ -151,7 +165,7 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
}
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];
@ -170,33 +184,39 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
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))
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 (which > 0) neighptr[n++] = jh ^ (which << SBBITS);
} else 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;
}
} else {
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular != Atom::ATOMIC) {
if (!moltemplate)
which = find_special(special[i],nspecial[i],tag[j]);
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))
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 (which > 0) neighptr[n++] = j ^ (which << SBBITS);
} else 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;
}
}
}
@ -206,8 +226,7 @@ void NPairBin<HALF, NEWTON, TRI, SIZE>::build(NeighList *list)
firstneigh[i] = neighptr;
numneigh[i] = n;
ipage->vgot(n);
if (ipage->status())
error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
if (ipage->status()) error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
}
list->inum = inum;