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

enable and process pair style table with clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-24 22:20:58 -04:00
parent d79c42ac41
commit 7ddfa382dc
1 changed files with 307 additions and 293 deletions
Download Patch File Download Diff File Expand all files Collapse all files

172
src/pair_table.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -18,20 +17,17 @@
#include "pair_table.h" #include "pair_table.h"
#include <cmath>
#include <cstring>
#include "atom.h" #include "atom.h"
#include "force.h"
#include "comm.h" #include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h" #include "error.h"
#include "force.h"
#include "tokenizer.h" #include "memory.h"
#include "neigh_list.h"
#include "table_file_reader.h" #include "table_file_reader.h"
#include "tokenizer.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -117,27 +113,27 @@ void PairTable::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
tb = &tables[tabindex[itype][jtype]]; tb = &tables[tabindex[itype][jtype]];
if (rsq < tb->innersq) if (rsq < tb->innersq)
error->one(FLERR,"Pair distance < table inner cutoff: " error->one(FLERR, "Pair distance < table inner cutoff: ijtype {} {} dist {}", itype,
"ijtype {} {} dist {}",itype,jtype,sqrt(rsq)); jtype, sqrt(rsq));
if (tabstyle == LOOKUP) { if (tabstyle == LOOKUP) {
itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta);
if (itable >= tlm1) if (itable >= tlm1)
error->one(FLERR,"Pair distance > table outer cutoff: " error->one(FLERR, "Pair distance > table outer cutoff: ijtype {} {} dist {}", itype,
"ijtype {} {} dist {}",itype,jtype,sqrt(rsq)); jtype, sqrt(rsq));
fpair = factor_lj * tb->f[itable]; fpair = factor_lj * tb->f[itable];
} else if (tabstyle == LINEAR) { } else if (tabstyle == LINEAR) {
itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta);
if (itable >= tlm1) if (itable >= tlm1)
error->one(FLERR,"Pair distance > table outer cutoff: " error->one(FLERR, "Pair distance > table outer cutoff: ijtype {} {} dist {}", itype,
"ijtype {} {} dist {}",itype,jtype,sqrt(rsq)); jtype, sqrt(rsq));
fraction = (rsq - tb->rsq[itable]) * tb->invdelta; fraction = (rsq - tb->rsq[itable]) * tb->invdelta;
value = tb->f[itable] + fraction * tb->df[itable]; value = tb->f[itable] + fraction * tb->df[itable];
fpair = factor_lj * value; fpair = factor_lj * value;
} else if (tabstyle == SPLINE) { } else if (tabstyle == SPLINE) {
itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta); itable = static_cast<int>((rsq - tb->innersq) * tb->invdelta);
if (itable >= tlm1) if (itable >= tlm1)
error->one(FLERR,"Pair distance > table outer cutoff: " error->one(FLERR, "Pair distance > table outer cutoff: ijtype {} {} dist {}", itype,
"ijtype {} {} dist {}",itype,jtype,sqrt(rsq)); jtype, sqrt(rsq));
b = (rsq - tb->rsq[itable]) * tb->invdelta; b = (rsq - tb->rsq[itable]) * tb->invdelta;
a = 1.0 - b; a = 1.0 - b;
value = a * tb->f[itable] + b * tb->f[itable + 1] + value = a * tb->f[itable] + b * tb->f[itable + 1] +
@ -174,8 +170,7 @@ void PairTable::compute(int eflag, int vflag)
evdwl *= factor_lj; evdwl *= factor_lj;
} }
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
evdwl,0.0,fpair,delx,dely,delz);
} }
} }
} }
@ -211,11 +206,16 @@ void PairTable::settings(int narg, char **arg)
// new settings // new settings
if (strcmp(arg[0],"lookup") == 0) tabstyle = LOOKUP; if (strcmp(arg[0], "lookup") == 0)
else if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR; tabstyle = LOOKUP;
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE; else if (strcmp(arg[0], "linear") == 0)
else if (strcmp(arg[0],"bitmap") == 0) tabstyle = BITMAP; tabstyle = LINEAR;
else error->all(FLERR,"Unknown table style in pair_style command: {}", arg[0]); else if (strcmp(arg[0], "spline") == 0)
tabstyle = SPLINE;
else if (strcmp(arg[0], "bitmap") == 0)
tabstyle = BITMAP;
else
error->all(FLERR, "Unknown table style in pair_style command: {}", arg[0]);
tablength = utils::inumeric(FLERR, arg[1], false, lmp); tablength = utils::inumeric(FLERR, arg[1], false, lmp);
if (tablength < 2) error->all(FLERR, "Illegal number of pair table entries"); if (tablength < 2) error->all(FLERR, "Illegal number of pair table entries");
@ -225,12 +225,18 @@ void PairTable::settings(int narg, char **arg)
int iarg = 2; int iarg = 2;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"ewald") == 0) ewaldflag = 1; if (strcmp(arg[iarg], "ewald") == 0)
else if (strcmp(arg[iarg],"pppm") == 0) pppmflag = 1; ewaldflag = 1;
else if (strcmp(arg[iarg],"msm") == 0) msmflag = 1; else if (strcmp(arg[iarg], "pppm") == 0)
else if (strcmp(arg[iarg],"dispersion") == 0) dispersionflag = 1; pppmflag = 1;
else if (strcmp(arg[iarg],"tip4p") == 0) tip4pflag = 1; else if (strcmp(arg[iarg], "msm") == 0)
else error->all(FLERR,"Illegal pair_style command"); msmflag = 1;
else if (strcmp(arg[iarg], "dispersion") == 0)
dispersionflag = 1;
else if (strcmp(arg[iarg], "tip4p") == 0)
tip4pflag = 1;
else
error->all(FLERR, "Illegal pair_style command");
iarg++; iarg++;
} }
@ -265,8 +271,7 @@ void PairTable::coeff(int narg, char **arg)
int me; int me;
MPI_Comm_rank(world, &me); MPI_Comm_rank(world, &me);
tables = (Table *) tables = (Table *) memory->srealloc(tables, (ntables + 1) * sizeof(Table), "pair:tables");
memory->srealloc(tables,(ntables+1)*sizeof(Table),"pair:tables");
Table *tb = &tables[ntables]; Table *tb = &tables[ntables];
null_table(tb); null_table(tb);
if (me == 0) read_table(tb, arg[2], arg[3]); if (me == 0) read_table(tb, arg[2], arg[3]);
@ -274,9 +279,12 @@ void PairTable::coeff(int narg, char **arg)
// set table cutoff // set table cutoff
if (narg == 5) tb->cut = utils::numeric(FLERR,arg[4],false,lmp); if (narg == 5)
else if (tb->rflag) tb->cut = tb->rhi; tb->cut = utils::numeric(FLERR, arg[4], false, lmp);
else tb->cut = tb->rfile[tb->ninput-1]; else if (tb->rflag)
tb->cut = tb->rhi;
else
tb->cut = tb->rfile[tb->ninput - 1];
// error check on table parameters // error check on table parameters
// insure cutoff is within table // insure cutoff is within table
@ -300,10 +308,10 @@ void PairTable::coeff(int narg, char **arg)
// for tabstyle SPLINE, always need to build spline tables // for tabstyle SPLINE, always need to build spline tables
tb->match = 0; tb->match = 0;
if (tabstyle == LINEAR && tb->ninput == tablength && if (tabstyle == LINEAR && tb->ninput == tablength && tb->rflag == RSQ && tb->rhi == tb->cut)
tb->rflag == RSQ && tb->rhi == tb->cut) tb->match = 1; tb->match = 1;
if (tabstyle == BITMAP && tb->ninput == 1 << tablength && if (tabstyle == BITMAP && tb->ninput == 1 << tablength && tb->rflag == BMP && tb->rhi == tb->cut)
tb->rflag == BMP && tb->rhi == tb->cut) tb->match = 1; tb->match = 1;
if (tb->rflag == BMP && tb->match == 0) if (tb->rflag == BMP && tb->match == 0)
error->all(FLERR, "Bitmapped table in file does not match requested table"); error->all(FLERR, "Bitmapped table in file does not match requested table");
@ -354,13 +362,10 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
// transparently convert units for supported conversions // transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert(); int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY, double conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
unit_convert);
char *line = reader.find_section_start(keyword); char *line = reader.find_section_start(keyword);
if (!line) { if (!line) { error->one(FLERR, "Did not find keyword in table file"); }
error->one(FLERR,"Did not find keyword in table file");
}
// read args on 2nd line of section // read args on 2nd line of section
// allocate table arrays for file values // allocate table arrays for file values
@ -410,8 +415,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
if (tb->rflag == RLINEAR) if (tb->rflag == RLINEAR)
rnew = tb->rlo + (tb->rhi - tb->rlo) * i / (tb->ninput - 1); rnew = tb->rlo + (tb->rhi - tb->rlo) * i / (tb->ninput - 1);
else if (tb->rflag == RSQ) { else if (tb->rflag == RSQ) {
rnew = tb->rlo*tb->rlo + rnew = tb->rlo * tb->rlo + (tb->rhi * tb->rhi - tb->rlo * tb->rlo) * i / (tb->ninput - 1);
(tb->rhi*tb->rhi - tb->rlo*tb->rlo)*i/(tb->ninput-1);
rnew = sqrt(rnew); rnew = sqrt(rnew);
} else if (tb->rflag == BMP) { } else if (tb->rflag == BMP) {
rsq_lookup.i = i << nshiftbits; rsq_lookup.i = i << nshiftbits;
@ -458,21 +462,24 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
} }
if (ferror) if (ferror)
error->warning(FLERR,"{} of {} force values in table {} are inconsistent " error->warning(FLERR,
"with -dE/dr.\nWARNING: Should only be flagged at " "{} of {} force values in table {} are inconsistent with -dE/dr.\n"
"inflection points",ferror,tb->ninput,keyword); "WARNING: Should only be flagged at inflection points",
ferror, tb->ninput, keyword);
// warn if re-computed distance values differ from file values // warn if re-computed distance values differ from file values
if (rerror) if (rerror)
error->warning(FLERR,"{} of {} distance values in table {} with relative " error->warning(FLERR,
"error\nWARNING: over {} to re-computed values", "{} of {} distance values in table {} with relative error\n"
"WARNING: over {} to re-computed values",
rerror, tb->ninput, EPSILONR, keyword); rerror, tb->ninput, EPSILONR, keyword);
// warn if data was read incompletely, e.g. columns were missing // warn if data was read incompletely, e.g. columns were missing
if (cerror) if (cerror)
error->warning(FLERR,"{} of {} lines in table {} were incomplete\n" error->warning(FLERR,
"{} of {} lines in table {} were incomplete\n"
"WARNING: or could not be parsed completely", "WARNING: or could not be parsed completely",
cerror, tb->ninput, keyword); cerror, tb->ninput, keyword);
} }
@ -556,9 +563,12 @@ void PairTable::param_extract(Table *tb, char *line)
if (word == "N") { if (word == "N") {
tb->ninput = values.next_int(); tb->ninput = values.next_int();
} else if ((word == "R") || (word == "RSQ") || (word == "BITMAP")) { } else if ((word == "R") || (word == "RSQ") || (word == "BITMAP")) {
if (word == "R") tb->rflag = RLINEAR; if (word == "R")
else if (word == "RSQ") tb->rflag = RSQ; tb->rflag = RLINEAR;
else if (word == "BITMAP") tb->rflag = BMP; else if (word == "RSQ")
tb->rflag = RSQ;
else if (word == "BITMAP")
tb->rflag = BMP;
tb->rlo = values.next_double(); tb->rlo = values.next_double();
tb->rhi = values.next_double(); tb->rhi = values.next_double();
} else if (word == "FPRIME") { } else if (word == "FPRIME") {
@ -589,8 +599,10 @@ void PairTable::compute_table(Table *tb)
// delta = table spacing in rsq for N-1 bins // delta = table spacing in rsq for N-1 bins
double inner; double inner;
if (tb->rflag) inner = tb->rlo; if (tb->rflag)
else inner = tb->rfile[0]; inner = tb->rlo;
else
inner = tb->rfile[0];
tb->innersq = inner * inner; tb->innersq = inner * inner;
tb->delta = (tb->cut * tb->cut - tb->innersq) / tlm1; tb->delta = (tb->cut * tb->cut - tb->innersq) / tlm1;
tb->invdelta = 1.0 / tb->delta; tb->invdelta = 1.0 / tb->delta;
@ -696,23 +708,24 @@ void PairTable::compute_table(Table *tb)
double fp0, fpn; double fp0, fpn;
double secant_factor = 0.1; double secant_factor = 0.1;
if (tb->fpflag) fp0 = (tb->fplo/sqrt(tb->innersq) - tb->f[0]/tb->innersq) / if (tb->fpflag)
(2.0 * sqrt(tb->innersq)); fp0 = (tb->fplo / sqrt(tb->innersq) - tb->f[0] / tb->innersq) / (2.0 * sqrt(tb->innersq));
else { else {
double rsq1 = tb->innersq; double rsq1 = tb->innersq;
double rsq2 = rsq1 + secant_factor * tb->delta; double rsq2 = rsq1 + secant_factor * tb->delta;
fp0 = (splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,sqrt(rsq2)) / fp0 = (splint(tb->rfile, tb->ffile, tb->f2file, tb->ninput, sqrt(rsq2)) / sqrt(rsq2) -
sqrt(rsq2) - tb->f[0] / sqrt(rsq1)) / (secant_factor*tb->delta); tb->f[0] / sqrt(rsq1)) /
(secant_factor * tb->delta);
} }
if (tb->fpflag && tb->cut == tb->rfile[tb->ninput-1]) fpn = if (tb->fpflag && tb->cut == tb->rfile[tb->ninput - 1])
(tb->fphi/tb->cut - tb->f[tlm1]/(tb->cut*tb->cut)) / (2.0 * tb->cut); fpn = (tb->fphi / tb->cut - tb->f[tlm1] / (tb->cut * tb->cut)) / (2.0 * tb->cut);
else { else {
double rsq2 = tb->cut * tb->cut; double rsq2 = tb->cut * tb->cut;
double rsq1 = rsq2 - secant_factor * tb->delta; double rsq1 = rsq2 - secant_factor * tb->delta;
fpn = (tb->f[tlm1] / sqrt(rsq2) - fpn = (tb->f[tlm1] / sqrt(rsq2) -
splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,sqrt(rsq1)) / splint(tb->rfile, tb->ffile, tb->f2file, tb->ninput, sqrt(rsq1)) / sqrt(rsq1)) /
sqrt(rsq1)) / (secant_factor*tb->delta); (secant_factor * tb->delta);
} }
for (int i = 0; i < tablength; i++) tb->f[i] /= sqrt(tb->rsq[i]); for (int i = 0; i < tablength; i++) tb->f[i] /= sqrt(tb->rsq[i]);
@ -856,14 +869,14 @@ void PairTable::free_table(Table *tb)
spline and splint routines modified from Numerical Recipes spline and splint routines modified from Numerical Recipes
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairTable::spline(double *x, double *y, int n, void PairTable::spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
double yp1, double ypn, double *y2)
{ {
int i, k; int i, k;
double p, qn, sig, un; double p, qn, sig, un;
double *u = new double[n]; double *u = new double[n];
if (yp1 > 0.99e30) y2[0] = u[0] = 0.0; if (yp1 > 0.99e30)
y2[0] = u[0] = 0.0;
else { else {
y2[0] = -0.5; y2[0] = -0.5;
u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1); u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
@ -875,7 +888,8 @@ void PairTable::spline(double *x, double *y, int n,
u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]); u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p; u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
} }
if (ypn > 0.99e30) qn = un = 0.0; if (ypn > 0.99e30)
qn = un = 0.0;
else { else {
qn = 0.5; qn = 0.5;
un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])); un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
@ -897,8 +911,10 @@ double PairTable::splint(double *xa, double *ya, double *y2a, int n, double x)
khi = n - 1; khi = n - 1;
while (khi - klo > 1) { while (khi - klo > 1) {
k = (khi + klo) >> 1; k = (khi + klo) >> 1;
if (xa[k] > x) khi = k; if (xa[k] > x)
else klo = k; khi = k;
else
klo = k;
} }
h = xa[khi] - xa[klo]; h = xa[khi] - xa[klo];
a = (xa[khi] - x) / h; a = (xa[khi] - x) / h;
@ -969,8 +985,7 @@ void PairTable::read_restart_settings(FILE *fp)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double PairTable::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, double PairTable::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double /*factor_coul*/, double factor_lj, double &fforce)
double &fforce)
{ {
int itable; int itable;
double fraction, value, a, b, phi; double fraction, value, a, b, phi;
@ -995,8 +1010,7 @@ double PairTable::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
b = (rsq - tb->rsq[itable]) * tb->invdelta; b = (rsq - tb->rsq[itable]) * tb->invdelta;
a = 1.0 - b; a = 1.0 - b;
value = a * tb->f[itable] + b * tb->f[itable + 1] + value = a * tb->f[itable] + b * tb->f[itable + 1] +
((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * ((a * a * a - a) * tb->f2[itable] + (b * b * b - b) * tb->f2[itable + 1]) * tb->deltasq6;
tb->deltasq6;
fforce = factor_lj * value; fforce = factor_lj * value;
} else { } else {
union_int_float_t rsq_lookup; union_int_float_t rsq_lookup;
@ -1036,9 +1050,9 @@ void *PairTable::extract(const char *str, int &dim)
double cut_coul = tables[0].cut; double cut_coul = tables[0].cut;
for (int m = 1; m < ntables; m++) for (int m = 1; m < ntables; m++)
if (tables[m].cut != cut_coul) if (tables[m].cut != cut_coul)
error->all(FLERR, error->all(FLERR, "Pair table cutoffs must all be equal to use with KSpace");
"Pair table cutoffs must all be equal to use with KSpace");
dim = 0; dim = 0;
return &tables[0].cut; return &tables[0].cut;
} else return nullptr; } else
return nullptr;
} }