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

Merge branch 'develop' into collected-small-changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-10-04 14:34:24 -04:00
parent bca8e6b85a dd072f7e08
commit 017ea1ace9
5 changed files with 321 additions and 138 deletions
Download Patch File Download Diff File Expand all files Collapse all files

282
src/MISC/pair_list.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -26,36 +25,42 @@
#include "memory.h"
#include "text_file_reader.h"
#include <cstring>
#include <cmath>
#include <cstring>
#include <exception>
#include <map>
using namespace LAMMPS_NS;
using MathSpecial::square;
enum { NONE = 0, HARM, MORSE, LJ126 };
enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC };
// clang-format off
static std::map<std::string, int> stylename = {
{ "none", NONE },
{ "harmonic", HARM },
{ "morse", MORSE },
{ "lj126", LJ126 }
{"none", NONE},
{"harmonic", HARM},
{"morse", MORSE},
{"lj126", LJ126},
{"quartic", QUARTIC}
};
// clang-format on
// fast power function for integer exponent > 0
static double mypow(double x, int n) {
static double mypow(double x, int n)
{
double yy;
if (x == 0.0) return 0.0;
for (yy = 1.0; n != 0; n >>= 1, x *=x)
for (yy = 1.0; n != 0; n >>= 1, x *= x)
if (n & 1) yy *= x;
return yy;
}
typedef struct { double x,y,z; } dbl3_t;
typedef struct {
double x, y, z;
} dbl3_t;
/* ---------------------------------------------------------------------- */
@ -86,19 +91,42 @@ PairList::~PairList()
void PairList::compute(int eflag, int vflag)
{
ev_init(eflag,vflag);
ev_init(eflag, vflag);
// get maximum allowed tag.
bigint maxtag_one, maxtag;
maxtag_one = maxtag = 0;
const int nlocal = atom->nlocal;
const int newton_pair = force->newton_pair;
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
dbl3_t * _noalias const f = (dbl3_t *) atom->f[0]; // NOLINT
const tagint *_noalias const tag = atom->tag;
for (int i = 0; i < nlocal; ++i) maxtag_one = MAX(maxtag_one, tag[i]);
MPI_Allreduce(&maxtag_one, &maxtag, 1, MPI_LMP_TAGINT, MPI_MAX, world);
double fpair,epair;
int i,j;
const int newton_pair = force->newton_pair;
const dbl3_t *_noalias const x = (dbl3_t *) atom->x[0];
dbl3_t *_noalias const f = (dbl3_t *) atom->f[0]; // NOLINT
double fpair, epair;
int i, j;
int pc = 0;
for (int n=0; n < npairs; ++n) {
for (int n = 0; n < npairs; ++n) {
const list_param &par = params[n];
// can only use valid tags or else atom->map() below will segfault.
if ((par.id1 < 1) || (par.id1 > maxtag)) {
if (check_flag)
error->all(FLERR, "Invalid pair list atom ID {}", par.id1);
else
continue;
}
if ((par.id2 < 1) || (par.id2 > maxtag)) {
if (check_flag)
error->all(FLERR, "Invalid pair list atom ID {}", par.id2);
else
continue;
}
i = atom->map(par.id1);
j = atom->map(par.id2);
@ -112,14 +140,14 @@ void PairList::compute(int eflag, int vflag)
// if id1 is a ghost, we skip if the sum of both ids is even.
// if id2 is a ghost, we skip if the sum of both ids is odd.
if (newton_pair) {
if ((i >= nlocal) && ((par.id1+par.id2) & 1) == 0) continue;
if ((j >= nlocal) && ((par.id1+par.id2) & 1) == 1) continue;
if ((i >= nlocal) && ((par.id1 + par.id2) & 1) == 0) continue;
if ((j >= nlocal) && ((par.id1 + par.id2) & 1) == 1) continue;
}
const double dx = x[i].x - x[j].x;
const double dy = x[i].y - x[j].y;
const double dz = x[i].z - x[j].z;
const double rsq = dx*dx + dy*dy + dz*dz;
const double rsq = dx * dx + dy * dy + dz * dz;
fpair = epair = 0.0;
if (check_flag) {
@ -128,61 +156,67 @@ void PairList::compute(int eflag, int vflag)
}
if (rsq < par.cutsq) {
const double r2inv = 1.0/rsq;
const double r2inv = 1.0 / rsq;
if (par.style == HARM) {
const double r = sqrt(rsq);
const double dr = par.param.harm.r0 - r;
fpair = 2.0*par.param.harm.k*dr/r;
fpair = 2.0 * par.param.harm.k * dr / r;
if (eflag_either)
epair = par.param.harm.k*dr*dr - par.offset;
if (eflag_either) epair = par.param.harm.k * dr * dr - par.offset;
} else if (par.style == MORSE) {
const double r = sqrt(rsq);
const double dr = par.param.morse.r0 - r;
const double dexp = exp(par.param.morse.alpha * dr);
fpair = 2.0*par.param.morse.d0*par.param.morse.alpha
* (dexp*dexp - dexp) / r;
const double dr = r - par.param.morse.r0;
const double dexp = exp(-par.param.morse.alpha * dr);
fpair = 2.0 * par.param.morse.d0 * par.param.morse.alpha * (dexp * dexp - dexp) / r;
if (eflag_either)
epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset;
if (eflag_either) epair = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp + 1.0) - par.offset;
} else if (par.style == LJ126) {
const double r6inv = r2inv*r2inv*r2inv;
const double sig6 = mypow(par.param.lj126.sigma,6);
fpair = 24.0*par.param.lj126.epsilon*r6inv
* (2.0*sig6*sig6*r6inv - sig6) * r2inv;
const double r6inv = r2inv * r2inv * r2inv;
const double sig6 = mypow(par.param.lj126.sigma, 6);
fpair = 24.0 * par.param.lj126.epsilon * r6inv * (2.0 * sig6 * sig6 * r6inv - sig6) * r2inv;
if (eflag_either)
epair = 4.0*par.param.lj126.epsilon*r6inv
* (sig6*sig6*r6inv - sig6) - par.offset;
epair = 4.0 * par.param.lj126.epsilon * r6inv * (sig6 * sig6 * r6inv - sig6) - par.offset;
} else if (par.style == QUARTIC) {
const double r = sqrt(rsq);
double dr = r - sqrt(par.cutsq);
double ra = dr - par.param.quartic.b1;
double rb = dr - par.param.quartic.b2;
double r2 = dr * dr;
fpair = -par.param.quartic.k / r * (r2 * (ra + rb) + 2.0 * dr * ra * rb);
if (eflag_either) epair = par.param.quartic.k * r2 * ra * rb;
}
if (newton_pair || i < nlocal) {
f[i].x += dx*fpair;
f[i].y += dy*fpair;
f[i].z += dz*fpair;
f[i].x += dx * fpair;
f[i].y += dy * fpair;
f[i].z += dz * fpair;
}
if (newton_pair || j < nlocal) {
f[j].x -= dx*fpair;
f[j].y -= dy*fpair;
f[j].z -= dz*fpair;
f[j].x -= dx * fpair;
f[j].y -= dy * fpair;
f[j].z -= dz * fpair;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,epair,0.0,fpair,dx,dy,dz);
if (evflag) ev_tally(i, j, nlocal, newton_pair, epair, 0.0, fpair, dx, dy, dz);
}
}
if (vflag_fdotr) virial_fdotr_compute();
if (check_flag) {
int tmp;
MPI_Allreduce(&pc,&tmp,1,MPI_INT,MPI_SUM,world);
if (tmp != 2*npairs)
error->all(FLERR,"Not all pairs processed in pair_style list");
MPI_Allreduce(&pc, &tmp, 1, MPI_INT, MPI_SUM, world);
if (tmp != 2 * npairs)
error->all(FLERR, "Not all pairs processed in pair_style list: {} vs {}", tmp, 2 * npairs);
}
}
@ -193,14 +227,13 @@ void PairList::compute(int eflag, int vflag)
void PairList::allocate()
{
allocated = 1;
int n = atom->ntypes;
int np1 = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutsq, np1, np1, "pair:cutsq");
}
/* ----------------------------------------------------------------------
@ -209,13 +242,19 @@ void PairList::allocate()
void PairList::settings(int narg, char **arg)
{
if (narg < 2)
error->all(FLERR,"Illegal pair_style command");
if (narg < 2) utils::missing_cmd_args(FLERR, "pair_style list", error);
cut_global = utils::numeric(FLERR,arg[1],false,lmp);
if (narg > 2) {
if (strcmp(arg[2],"nocheck") == 0) check_flag = 0;
if (strcmp(arg[2],"check") == 0) check_flag = 1;
cut_global = utils::numeric(FLERR, arg[1], false, lmp);
int iarg = 2;
while (iarg < narg) {
if (strcmp(arg[iarg], "nocheck") == 0) {
check_flag = 0;
++iarg;
} else if (strcmp(arg[2], "check") == 0) {
check_flag = 1;
++iarg;
} else
error->all(FLERR, "Unknown pair style list keyword: {}", arg[iarg]);
}
std::vector<int> mystyles;
@ -223,13 +262,15 @@ void PairList::settings(int narg, char **arg)
// read and parse potential file only on MPI rank 0.
if (comm->me == 0) {
int nharm, nmorse, nlj126, nskipped;
FILE *fp = utils::open_potential(arg[0],lmp,nullptr);
TextFileReader reader(fp,"pair list coeffs");
npairs = nharm = nmorse = nlj126 = nskipped = 0;
int nharm, nmorse, nlj126, nquartic, nskipped;
FILE *fp = utils::open_potential(arg[0], lmp, nullptr);
if (!fp)
error->one(FLERR, "Error opening pair list coeffs file {}: {}", arg[0], utils::getsyserror());
TextFileReader reader(fp, "pair list coeffs");
npairs = nharm = nmorse = nlj126 = nquartic = nskipped = 0;
char *line;
try {
char *line;
while ((line = reader.next_line())) {
ValueTokenizer values(line);
list_param oneparam;
@ -241,75 +282,86 @@ void PairList::settings(int narg, char **arg)
switch (oneparam.style) {
case HARM:
oneparam.param.harm.k = values.next_double();
oneparam.param.harm.r0 = values.next_double();
++nharm;
break;
case HARM:
oneparam.param.harm.k = values.next_double();
oneparam.param.harm.r0 = values.next_double();
++nharm;
break;
case MORSE:
oneparam.param.morse.d0 = values.next_double();
oneparam.param.morse.alpha = values.next_double();
oneparam.param.morse.r0 = values.next_double();
++nmorse;
break;
case MORSE:
oneparam.param.morse.d0 = values.next_double();
oneparam.param.morse.alpha = values.next_double();
oneparam.param.morse.r0 = values.next_double();
++nmorse;
break;
case LJ126:
oneparam.param.lj126.epsilon = values.next_double();
oneparam.param.lj126.sigma = values.next_double();
++nlj126;
break;
case LJ126:
oneparam.param.lj126.epsilon = values.next_double();
oneparam.param.lj126.sigma = values.next_double();
++nlj126;
break;
case NONE: // fallthrough
error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}",
utils::trim(line));
++nskipped;
break;
case QUARTIC:
oneparam.param.quartic.k = values.next_double();
oneparam.param.quartic.b1 = values.next_double();
oneparam.param.quartic.b2 = values.next_double();
if (!values.has_next())
throw FileReaderException("Must specify individual cutoff for quartic interaction");
++nquartic;
break;
case NONE: // fallthrough
error->warning(FLERR, "Skipping unrecognized pair list potential entry: {}",
utils::trim(line));
++nskipped;
break;
}
if (values.has_next())
oneparam.cutsq = square(values.next_double());
else
oneparam.cutsq = cut_global*cut_global;
oneparam.cutsq = cut_global * cut_global;
myparams.push_back(oneparam);
}
} catch (std::exception &e) {
error->one(FLERR,"Error reading pair list coeffs file: {}", e.what());
error->one(FLERR, "Error reading pair list coeffs file: {}\n{}", e.what(), line);
}
utils::logmesg(lmp, "Read {} ({}/{}/{}) interacting pair lines from {}. "
"{} skipped entries.\n", npairs, nharm, nmorse, nlj126, arg[0], nskipped);
utils::logmesg(lmp,
"Read {} ({}/{}/{}/{}) interacting pair lines from {}. "
"{} skipped entries.\n",
npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped);
memory->create(params,npairs,"pair_list:params");
memcpy(params, myparams.data(),npairs*sizeof(list_param));
memory->create(params, npairs, "pair_list:params");
memcpy(params, myparams.data(), npairs * sizeof(list_param));
fclose(fp);
}
MPI_Bcast(&npairs, 1, MPI_INT, 0, world);
if (comm->me != 0) memory->create(params,npairs,"pair_list:params");
MPI_Bcast(params, npairs*sizeof(list_param), MPI_BYTE, 0, world);
if (comm->me != 0) memory->create(params, npairs, "pair_list:params");
MPI_Bcast(params, npairs * sizeof(list_param), MPI_BYTE, 0, world);
}
/* ----------------------------------------------------------------------
there are no coeffs to be set, but we need to update setflag and pretend
there are no coeffs to be set, but we need to update setflag and pretend there are
------------------------------------------------------------------------- */
void PairList::coeff(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 2) utils::missing_cmd_args(FLERR, "pair_coeff list", error);
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -318,29 +370,31 @@ void PairList::coeff(int narg, char **arg)
void PairList::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style list requires atom IDs");
if (atom->tag_enable == 0) error->all(FLERR, "Pair style list requires atom IDs");
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Pair style list requires an atom map");
if (atom->map_style == Atom::MAP_NONE) error->all(FLERR, "Pair style list requires an atom map");
if (offset_flag) {
for (int n=0; n < npairs; ++n) {
for (int n = 0; n < npairs; ++n) {
list_param &par = params[n];
if (par.style == HARM) {
const double dr = sqrt(par.cutsq) - par.param.harm.r0;
par.offset = par.param.harm.k*dr*dr;
par.offset = par.param.harm.k * dr * dr;
} else if (par.style == MORSE) {
const double dr = par.param.morse.r0 - sqrt(par.cutsq);
const double dexp = exp(par.param.morse.alpha * dr);
par.offset = par.param.morse.d0 * (dexp*dexp - 2.0*dexp);
par.offset = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp - 1.0);
} else if (par.style == LJ126) {
const double r6inv = par.cutsq*par.cutsq*par.cutsq;
const double sig6 = mypow(par.param.lj126.sigma,6);
par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6);
const double r6inv = par.cutsq * par.cutsq * par.cutsq;
const double sig6 = mypow(par.param.lj126.sigma, 6);
par.offset = 4.0 * par.param.lj126.epsilon * r6inv * (sig6 * sig6 * r6inv - sig6);
} else if (par.style == QUARTIC) {
// the offset is always 0 at rc
par.offset = 0.0;
}
}
}
@ -362,10 +416,10 @@ double PairList::init_one(int, int)
double PairList::memory_usage()
{
double bytes = (double)npairs * sizeof(int);
bytes += (double)npairs * sizeof(list_param);
const int n = atom->ntypes+1;
bytes += (double)n*(n*sizeof(int) + sizeof(int *));
bytes += (double)n*(n*sizeof(double) + sizeof(double *));
double bytes = (double) npairs * sizeof(int);
bytes += (double) npairs * sizeof(list_param);
const int n = atom->ntypes + 1;
bytes += (double) n * (n * sizeof(int) + sizeof(int *));
bytes += (double) n * (n * sizeof(double) + sizeof(double *));
return bytes;
}

5
src/MISC/pair_list.h
View File

@ -49,11 +49,16 @@ class PairList : public Pair {
struct lj126_p {
double epsilon, sigma;
};
struct quartic_p {
double k, b1, b2;
};
union param_u {
harm_p harm;
morse_p morse;
lj126_p lj126;
quartic_p quartic;
};
struct list_param {