use pairwise cutoff for trimming unless a custom cutoff if given

This commit is contained in:
Axel Kohlmeyer
2025-03-05 04:17:38 -05:00
parent b433c8cb77
commit 47cea8e1ba
5 changed files with 47 additions and 25 deletions

View File

@ -13,27 +13,29 @@
------------------------------------------------------------------------- */
#include "npair.h"
#include <cmath>
#include "neighbor.h"
#include "neigh_request.h"
#include "nbin.h"
#include "nstencil.h"
#include "atom.h"
#include "update.h"
#include "memory.h"
#include "error.h"
#include "memory.h"
#include "nbin.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "nstencil.h"
#include "update.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NPair::NPair(LAMMPS *lmp)
: Pointers(lmp), nb(nullptr), ns(nullptr), bins(nullptr), stencil(nullptr)
NPair::NPair(LAMMPS *lmp) : Pointers(lmp), nb(nullptr), ns(nullptr), bins(nullptr), stencil(nullptr)
{
last_build = -1;
mycutneighsq = nullptr;
molecular = atom->molecular;
copymode = 0;
cutoff_custom = 0.0;
execution_space = Host;
}
@ -50,7 +52,6 @@ NPair::~NPair()
void NPair::post_constructor(NeighRequest *nrq)
{
cutoff_custom = 0.0;
if (nrq->cut) cutoff_custom = nrq->cutoff;
}

View File

@ -50,7 +50,7 @@ NPairHalffull<NEWTON, TRI, TRIM>::NPairHalffull(LAMMPS *lmp) : NPair(lmp) {}
template<int NEWTON, int TRI, int TRIM>
void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
{
int i, j, ii, jj, n, jnum, joriginal;
int i, j, ii, jj, n, jnum, joriginal, itype, jtype;
int *neighptr, *jlist;
double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
@ -58,6 +58,7 @@ void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
double **x = atom->x;
int nlocal = atom->nlocal;
int *type = atom->type;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -85,6 +86,7 @@ void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
// loop over parent full list
i = ilist_full[ii];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
@ -95,6 +97,7 @@ void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
jtype = type[j];
if (NEWTON) {
if (j < nlocal) {
@ -121,7 +124,8 @@ void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cutneighsq[itype][jtype];
if (rsq > cutsq_trim) continue;
}
neighptr[n++] = joriginal;
@ -133,7 +137,8 @@ void NPairHalffull<NEWTON, TRI, TRIM>::build(NeighList *list)
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cutneighsq[itype][jtype];
if (rsq > cutsq_trim) continue;
}
neighptr[n++] = joriginal;

View File

@ -37,7 +37,7 @@ NPairSkipTemp<TRIM>::NPairSkipTemp(LAMMPS *lmp) : NPair(lmp) {}
template<int TRIM>
void NPairSkipTemp<TRIM>::build(NeighList *list)
{
int i, j, ii, jj, n, itype, jnum, joriginal;
int i, j, ii, jj, n, itype, jtype, jnum, joriginal;
int *neighptr, *jlist;
int *type = atom->type;
@ -94,7 +94,8 @@ void NPairSkipTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (!molskip && ijskip[itype][type[j]]) continue;
jtype = type[j];
if (!molskip && ijskip[itype][jtype]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
@ -103,7 +104,9 @@ void NPairSkipTemp<TRIM>::build(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cutneighsq[itype][jtype];
if (rsq > cutsq_trim) continue;
}
neighptr[n++] = joriginal;

View File

@ -36,7 +36,7 @@ NPairSkipRespaTemp<TRIM>::NPairSkipRespaTemp(LAMMPS *lmp) : NPair(lmp) {}
template<int TRIM>
void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
{
int i, j, ii, jj, n, itype, jnum, joriginal, n_inner, n_middle;
int i, j, ii, jj, n, itype, jtype, jnum, joriginal, n_inner, n_middle;
int *neighptr, *jlist, *neighptr_inner, *neighptr_middle;
int *type = atom->type;
@ -117,7 +117,8 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (!molskip && ijskip[itype][type[j]]) continue;
jtype = type[j];
if (!molskip && ijskip[itype][jtype]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
@ -126,7 +127,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cutneighsq[itype][jtype];
if (rsq > cutsq_trim) continue;
}
neighptr[n++] = joriginal;
@ -140,7 +143,8 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (!molskip && ijskip[itype][type[j]]) continue;
jtype = type[j];
if (!molskip && ijskip[itype][jtype]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
@ -149,7 +153,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cut_inner_sq;
if (rsq > cutsq_trim) continue;
}
neighptr_inner[n_inner++] = joriginal;
@ -164,7 +170,8 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
j = joriginal & NEIGHMASK;
if (!molskip && ijskip[itype][type[j]]) continue;
jtype = type[j];
if (!molskip && ijskip[itype][jtype]) continue;
if ((molskip == NeighRequest::INTRA) && (molecule[i] != molecule[j])) continue;
if ((molskip == NeighRequest::INTER) && (molecule[i] == molecule[j])) continue;
@ -173,7 +180,9 @@ void NPairSkipRespaTemp<TRIM>::build(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cut_middle_sq;
if (rsq > cutsq_trim) continue;
}
neighptr_middle[n_middle++] = joriginal;

View File

@ -34,12 +34,13 @@ void NPairTrim::build(NeighList *list)
double cutsq_custom = cutoff_custom * cutoff_custom;
int ii, jj, n, jnum, joriginal;
int ii, jj, n, jnum, joriginal, itype, jtype;
int *neighptr, *jlist;
double xtmp, ytmp, ztmp;
double delx, dely, delz, rsq;
double **x = atom->x;
int *type = atom->type;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
@ -64,6 +65,7 @@ void NPairTrim::build(NeighList *list)
neighptr = ipage->vget();
const int i = ilist_copy[ii];
itype = type[i];
ilist[ii] = i;
xtmp = x[i][0];
ytmp = x[i][1];
@ -77,13 +79,15 @@ void NPairTrim::build(NeighList *list)
for (jj = 0; jj < jnum; jj++) {
joriginal = jlist[jj];
const int j = joriginal & NEIGHMASK;
jtype = type[j];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
if (rsq > cutsq_custom) continue;
double cutsq_trim = (cutsq_custom > 0.0) ? cutsq_custom : cutneighsq[itype][jtype];
if (rsq > cutsq_trim) continue;
neighptr[n++] = joriginal;
}