From b12facbd5470d162b4ff779770bb25621f337d47 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 22 Feb 2023 18:46:31 -0500 Subject: [PATCH] add fix wall/table --- src/LEPTON/fix_wall_lepton.cpp | 3 +- src/fix_wall.cpp | 33 ++- src/fix_wall.h | 2 +- src/fix_wall_table.cpp | 467 +++++++++++++++++++++++++++++++++ src/fix_wall_table.h | 67 +++++ 5 files changed, 552 insertions(+), 20 deletions(-) create mode 100644 src/fix_wall_table.cpp create mode 100644 src/fix_wall_table.h diff --git a/src/LEPTON/fix_wall_lepton.cpp b/src/LEPTON/fix_wall_lepton.cpp index 8438fa8852..a81d3c4edb 100644 --- a/src/LEPTON/fix_wall_lepton.cpp +++ b/src/LEPTON/fix_wall_lepton.cpp @@ -93,8 +93,7 @@ void FixWallLepton::wall_particle(int m, int which, double coord) ; // do nothing } - double delta, fwall; - double vn; + double delta, fwall, vn; double **x = atom->x; double **f = atom->f; diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index cbea053dab..f7293ffd06 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -52,21 +52,20 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall fldflag = 0; int pbcflag = 0; - for (int i = 0; i < 6; i++) xstr[i] = estr[i] = sstr[i] = lstr[i] = tstr[i] = nullptr; + for (int i = 0; i < 6; i++) xstr[i] = estr[i] = sstr[i] = lstr[i] = fstr[i] = kstr[i] = nullptr; int iarg = 3; - tabfile = nullptr; - if (utils::strmatch(style, "^wall/table")) { - if (iarg + 1 > narg) error->all(FLERR, "Missing argument for fix {} command", style); - tabfile = arg[iarg]; - ++iarg; - } + if (utils::strmatch(style, "^wall/table")) iarg = 5; while (iarg < narg) { + int wantargs = 5; + if (utils::strmatch(style, "^wall/lepton")) wantargs = 4; + if (utils::strmatch(style, "^wall/morse")) wantargs = 6; + if ((strcmp(arg[iarg], "xlo") == 0) || (strcmp(arg[iarg], "xhi") == 0) || (strcmp(arg[iarg], "ylo") == 0) || (strcmp(arg[iarg], "yhi") == 0) || (strcmp(arg[iarg], "zlo") == 0) || (strcmp(arg[iarg], "zhi") == 0)) { - if (iarg + 4 > narg) error->all(FLERR, "Missing argument for fix {} command", style); + if (iarg + wantargs > narg) error->all(FLERR, "Missing argument for fix {} command", style); int newwall; if (strcmp(arg[iarg], "xlo") == 0) { @@ -108,13 +107,10 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall if (utils::strmatch(style, "^wall/lepton")) { lstr[nwall] = utils::strdup(arg[iarg + 2]); cutoff[nwall] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - nwall++; - iarg += 4; } else if (utils::strmatch(style, "^wall/table")) { - tstr[nwall] = utils::strdup(arg[iarg + 2]); - cutoff[nwall] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - nwall++; - iarg += 4; + fstr[nwall] = utils::strdup(arg[iarg + 2]); + kstr[nwall] = utils::strdup(arg[iarg + 3]); + cutoff[nwall] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); } else { if (iarg + 5 > narg) error->all(FLERR, "Missing argument for fix {} command", style); @@ -134,7 +130,9 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall alpha[nwall] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); astyle[nwall] = CONSTANT; } + // adjust so we can share the regular code path ++iarg; + --wantargs; } if (utils::strmatch(arg[iarg + 3], "^v_")) { @@ -145,9 +143,9 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), nwall sstyle[nwall] = CONSTANT; } cutoff[nwall] = utils::numeric(FLERR, arg[iarg + 4], false, lmp); - nwall++; - iarg += 5; } + nwall++; + iarg += wantargs; } else if (strcmp(arg[iarg], "units") == 0) { if (iarg + 2 > narg) error->all(FLERR, "Illegal fix {} command", style); if (strcmp(arg[iarg + 1], "box") == 0) @@ -248,7 +246,8 @@ FixWall::~FixWall() delete[] estr[m]; delete[] sstr[m]; delete[] lstr[m]; - delete[] tstr[m]; + delete[] fstr[m]; + delete[] kstr[m]; } } diff --git a/src/fix_wall.h b/src/fix_wall.h index 8fd0ef6910..12ceb17b49 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -50,7 +50,7 @@ class FixWall : public Fix { double xscale, yscale, zscale; int estyle[6], sstyle[6], astyle[6], wstyle[6]; int eindex[6], sindex[6]; - char *estr[6], *sstr[6], *astr[6], *lstr[6], *tstr[6], *tabfile; + char *estr[6], *sstr[6], *astr[6], *lstr[6], *fstr[6], *kstr[6]; int varflag; // 1 if any wall position,epsilon,sigma is a variable int eflag; // per-wall flag for energy summation int ilevel_respa; diff --git a/src/fix_wall_table.cpp b/src/fix_wall_table.cpp new file mode 100644 index 0000000000..c44864a7bd --- /dev/null +++ b/src/fix_wall_table.cpp @@ -0,0 +1,467 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_wall_table.h" +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "memory.h" +#include "table_file_reader.h" +#include "tokenizer.h" + +#include + +using namespace LAMMPS_NS; + +enum { NONE, LINEAR, SPLINE }; + +static constexpr double BIGNUM = 1.0e300; + +/* ---------------------------------------------------------------------- */ + +FixWallTable::FixWallTable(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg), tables(nullptr) +{ + dynamic_group_allow = 1; + tabstyle = NONE; + if (strcmp(arg[3], "linear") == 0) + tabstyle = LINEAR; + else if (strcmp(arg[3], "spline") == 0) + tabstyle = SPLINE; + else + error->all(FLERR, "Unknown table style {} in fix {}", arg[3], style); + + tablength = utils::inumeric(FLERR, arg[4], false, lmp); + if (tablength < 2) error->all(FLERR, "Illegal number of fix {} table entries", style); +} + +/* ---------------------------------------------------------------------- */ + +void FixWallTable::post_constructor() +{ + memory->sfree(tables); + tables = (Table *) memory->smalloc(nwall * sizeof(Table), "wall:tables"); + + for (int m = 0; m < nwall; ++m) { + Table &tb = tables[m]; + null_table(tb); + + if (comm->me == 0) read_table(tb, fstr[m], kstr[m]); + bcast_table(tb); + + // error check on table parameters + if (tb.ninput <= 1) error->all(FLERR, "Invalid fix {} table length: {}", style, tb.ninput); + tb.lo = tb.rfile[0]; + tb.hi = tb.rfile[tb.ninput - 1]; + if (tb.lo >= tb.hi) error->all(FLERR, "Fix {} table distance values do not increase", style); + if (cutoff[m] > tb.hi) + error->all(FLERR, "Fix {} wall cutoff {} is larger than table outer cutoff {}", style, + cutoff[m], tb.hi); + + // spline read-in data and compute r,e,f arrays within table + spline_table(tb); + compute_table(tb); + } +} + +/* ---------------------------------------------------------------------- */ + +FixWallTable::~FixWallTable() +{ + for (int m = 0; m < nwall; m++) free_table(tables[m]); + memory->sfree(tables); +} + +/* ---------------------------------------------------------------------- + compute the potential energy offset so it can be shifted to zero at the cutoff +------------------------------------------------------------------------- */ + +void FixWallTable::precompute(int m) +{ + double u, mdu; + + uf_lookup(m, cutoff[m], u, mdu); + offset[m] = u; +} + +/* ---------------------------------------------------------------------- + interaction of all particles in group with a wall + m = index of wall coeffs + which = xlo,xhi,ylo,yhi,zlo,zhi + error if any particle is on or behind wall +------------------------------------------------------------------------- */ + +void FixWallTable::wall_particle(int m, int which, double coord) +{ + double delta, fwall, vn, u, mdu; + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int dim = which / 2; + int side = which % 2; + if (side == 0) side = -1; + + int onflag = 0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (side < 0) + delta = x[i][dim] - coord; + else + delta = coord - x[i][dim]; + if (delta >= cutoff[m]) continue; + if (delta <= 0.0) { + onflag = 1; + continue; + } + + uf_lookup(m, delta, u, mdu); + fwall = side * mdu; + f[i][dim] -= fwall; + ewall[0] += u - offset[m]; + ewall[m + 1] += fwall; + + if (evflag) { + if (side < 0) + vn = -fwall * delta; + else + vn = fwall * delta; + v_tally(dim, i, vn); + } + } + + if (onflag) error->one(FLERR, "Particle on or inside fix {} surface", style); +} + +/* ---------------------------------------------------------------------- */ + +void FixWallTable::null_table(Table &tb) +{ + tb.rfile = tb.efile = tb.ffile = nullptr; + tb.e2file = tb.f2file = nullptr; + tb.r = tb.e = tb.de = nullptr; + tb.f = tb.df = tb.e2 = tb.f2 = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallTable::free_table(Table &tb) +{ + memory->destroy(tb.rfile); + memory->destroy(tb.efile); + memory->destroy(tb.ffile); + memory->destroy(tb.e2file); + memory->destroy(tb.f2file); + + memory->destroy(tb.r); + memory->destroy(tb.e); + memory->destroy(tb.de); + memory->destroy(tb.f); + memory->destroy(tb.df); + memory->destroy(tb.e2); + memory->destroy(tb.f2); +} + +/* ---------------------------------------------------------------------- */ + +void FixWallTable::read_table(Table &tb, const std::string &tabfile, const std::string &keyword) +{ + TableFileReader reader(lmp, tabfile, "fix wall/table"); + + char *line = reader.find_section_start(keyword); + if (!line) error->one(FLERR, "Did not find keyword {} in table file", keyword); + + // read args on 2nd line of section + // allocate table arrays for file values + + line = reader.next_line(); + param_extract(tb, line); + memory->create(tb.rfile, tb.ninput, "wall:rfile"); + memory->create(tb.efile, tb.ninput, "wall:efile"); + memory->create(tb.ffile, tb.ninput, "wall:ffile"); + + // read r,e,f table values from file + + reader.skip_line(); + for (int i = 0; i < tb.ninput; i++) { + line = reader.next_line(); + if (!line) + error->one(FLERR, "Data missing when parsing wall table '{}' line {} of {}.", keyword, i + 1, + tb.ninput); + try { + ValueTokenizer values(line); + values.next_int(); + tb.rfile[i] = values.next_double(); + tb.efile[i] = values.next_double(); + tb.ffile[i] = values.next_double(); + } catch (TokenizerException &e) { + error->one(FLERR, "Error parsing wall table '{}' line {} of {}. {}\nLine was: {}", keyword, + i + 1, tb.ninput, e.what(), line); + } + } + + // warn if force != dE/dr at any point that is not an inflection point + // check via secant approximation to dE/dr + // skip two end points since do not have surrounding secants + // inflection point is where curvature changes sign + + double r, e, f, rprev, rnext, eprev, enext, fleft, fright; + + int ferror = 0; + for (int i = 1; i < tb.ninput - 1; i++) { + r = tb.rfile[i]; + rprev = tb.rfile[i - 1]; + rnext = tb.rfile[i + 1]; + e = tb.efile[i]; + eprev = tb.efile[i - 1]; + enext = tb.efile[i + 1]; + f = tb.ffile[i]; + fleft = -(e - eprev) / (r - rprev); + fright = -(enext - e) / (rnext - r); + if (f < fleft && f < fright) ferror++; + if (f > fleft && f > fright) ferror++; + } + + if (ferror) + error->warning(FLERR, + "{} of {} force values in table are inconsistent with -dE/dr.\n" + "WARNING: Should only be flagged at inflection points", + ferror, tb.ninput); +} + +/* ---------------------------------------------------------------------- + build spline representation of e,f over entire range of read-in table + this function sets these values in e2file,f2file +------------------------------------------------------------------------- */ + +void FixWallTable::spline_table(Table &tb) +{ + memory->create(tb.e2file, tb.ninput, "wall:e2file"); + memory->create(tb.f2file, tb.ninput, "wall:f2file"); + + double ep0 = -tb.ffile[0]; + double epn = -tb.ffile[tb.ninput - 1]; + spline(tb.rfile, tb.efile, tb.ninput, ep0, epn, tb.e2file); + + if (tb.fpflag == 0) { + tb.fplo = (tb.ffile[1] - tb.ffile[0]) / (tb.rfile[1] - tb.rfile[0]); + tb.fphi = (tb.ffile[tb.ninput - 1] - tb.ffile[tb.ninput - 2]) / + (tb.rfile[tb.ninput - 1] - tb.rfile[tb.ninput - 2]); + } + + double fp0 = tb.fplo; + double fpn = tb.fphi; + spline(tb.rfile, tb.ffile, tb.ninput, fp0, fpn, tb.f2file); +} + +/* ---------------------------------------------------------------------- + compute r,e,f vectors from splined values +------------------------------------------------------------------------- */ + +void FixWallTable::compute_table(Table &tb) +{ + // delta = table spacing for N-1 bins + int tlm1 = tablength - 1; + + tb.delta = (tb.hi - tb.lo) / tlm1; + tb.invdelta = 1.0 / tb.delta; + tb.deltasq6 = tb.delta * tb.delta / 6.0; + + // N-1 evenly spaced bins in r from min to max + // r,e,f = value at lower edge of bin + // de,df values = delta values of e,f + // r,e,f are N in length so de,df arrays can compute difference + + memory->create(tb.r, tablength, "wall:r"); + memory->create(tb.e, tablength, "wall:e"); + memory->create(tb.de, tablength, "wall:de"); + memory->create(tb.f, tablength, "wall:f"); + memory->create(tb.df, tablength, "wall:df"); + memory->create(tb.e2, tablength, "wall:e2"); + memory->create(tb.f2, tablength, "wall:f2"); + + double a; + for (int i = 0; i < tablength; i++) { + a = tb.lo + i * tb.delta; + tb.r[i] = a; + tb.e[i] = splint(tb.rfile, tb.efile, tb.e2file, tb.ninput, a); + tb.f[i] = splint(tb.rfile, tb.ffile, tb.f2file, tb.ninput, a); + } + + for (int i = 0; i < tlm1; i++) { + tb.de[i] = tb.e[i + 1] - tb.e[i]; + tb.df[i] = tb.f[i + 1] - tb.f[i]; + } + // get final elements from linear extrapolation + tb.de[tlm1] = 2.0 * tb.de[tlm1 - 1] - tb.de[tlm1 - 2]; + tb.df[tlm1] = 2.0 * tb.df[tlm1 - 1] - tb.df[tlm1 - 2]; + + double ep0 = -tb.f[0]; + double epn = -tb.f[tlm1]; + spline(tb.r, tb.e, tablength, ep0, epn, tb.e2); + spline(tb.r, tb.f, tablength, tb.fplo, tb.fphi, tb.f2); +} + +/* ---------------------------------------------------------------------- + extract attributes from parameter line in table section + format of line: N value FP fplo fphi + N is required, other params are optional +------------------------------------------------------------------------- */ + +void FixWallTable::param_extract(Table &tb, char *line) +{ + tb.ninput = 0; + tb.fpflag = 0; + + try { + ValueTokenizer values(line); + + while (values.has_next()) { + std::string word = values.next_string(); + + if (word == "N") { + tb.ninput = values.next_int(); + } else if (word == "FP") { + tb.fpflag = 1; + tb.fplo = values.next_double(); + tb.fphi = values.next_double(); + } else { + error->one(FLERR, "Invalid keyword {} in fix {} table parameters", word, style); + } + } + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); + } + + if (tb.ninput == 0) error->one(FLERR, "Fix {} parameters did not set N", style); +} + +/* ---------------------------------------------------------------------- + broadcast read-in table info from proc 0 to other procs + this function communicates these values in Table: + ninput,rfile,efile,ffile,fpflag,fplo,fphi +------------------------------------------------------------------------- */ + +void FixWallTable::bcast_table(Table &tb) +{ + MPI_Bcast(&tb.ninput, 1, MPI_INT, 0, world); + + if (comm->me > 0) { + memory->create(tb.rfile, tb.ninput, "wall:rfile"); + memory->create(tb.efile, tb.ninput, "wall:efile"); + memory->create(tb.ffile, tb.ninput, "wall:ffile"); + } + + MPI_Bcast(tb.rfile, tb.ninput, MPI_DOUBLE, 0, world); + MPI_Bcast(tb.efile, tb.ninput, MPI_DOUBLE, 0, world); + MPI_Bcast(tb.ffile, tb.ninput, MPI_DOUBLE, 0, world); + + MPI_Bcast(&tb.fpflag, 1, MPI_INT, 0, world); + if (tb.fpflag) { + MPI_Bcast(&tb.fplo, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&tb.fphi, 1, MPI_DOUBLE, 0, world); + } +} + +/* ---------------------------------------------------------------------- + spline and splint routines modified from Numerical Recipes +------------------------------------------------------------------------- */ + +void FixWallTable::spline(double *x, double *y, int n, double yp1, double ypn, double *y2) +{ + int i, k; + double p, qn, sig, un; + auto u = new double[n]; + + if (yp1 > BIGNUM) + y2[0] = u[0] = 0.0; + else { + y2[0] = -0.5; + u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1); + } + for (i = 1; i < n - 1; i++) { + sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]); + p = sig * y2[i - 1] + 2.0; + y2[i] = (sig - 1.0) / p; + 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; + } + if (ypn > BIGNUM) + qn = un = 0.0; + else { + 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])); + } + y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0); + for (k = n - 2; k >= 0; k--) y2[k] = y2[k] * y2[k + 1] + u[k]; + + delete[] u; +} + +/* ---------------------------------------------------------------------- */ + +double FixWallTable::splint(double *xa, double *ya, double *y2a, int n, double x) +{ + int klo, khi, k; + double h, b, a, y; + + klo = 0; + khi = n - 1; + while (khi - klo > 1) { + k = (khi + klo) >> 1; + if (xa[k] > x) + khi = k; + else + klo = k; + } + h = xa[khi] - xa[klo]; + a = (xa[khi] - x) / h; + b = (x - xa[klo]) / h; + y = a * ya[klo] + b * ya[khi] + + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0; + return y; +} + +/* ---------------------------------------------------------------------- + calculate potential u and force f at distance x + ensure x is between wall min/max, exit with error if not +------------------------------------------------------------------------- */ + +void FixWallTable::uf_lookup(int m, double x, double &u, double &f) +{ + double fraction, a, b; + const Table &tb = tables[m]; + const int itable = static_cast((x - tb.lo) * tb.invdelta); + if (itable < 0) + error->one(FLERR, "Particle / wall distance < table inner cutoff: {:.8}", x); + else if (itable >= tablength) + error->one(FLERR, "Particle / wall distance > table outer cutoff: {:.8}", x); + + if (tabstyle == LINEAR) { + fraction = (x - tb.r[itable]) * tb.invdelta; + u = tb.e[itable] + fraction * tb.de[itable]; + f = tb.f[itable] + fraction * tb.df[itable]; + } else if (tabstyle == SPLINE) { + fraction = (x - tb.r[itable]) * tb.invdelta; + + b = (x - tb.r[itable]) * tb.invdelta; + a = 1.0 - b; + u = a * tb.e[itable] + b * tb.e[itable + 1] + + ((a * a * a - a) * tb.e2[itable] + (b * b * b - b) * tb.e2[itable + 1]) * tb.deltasq6; + f = a * tb.f[itable] + b * tb.f[itable + 1] + + ((a * a * a - a) * tb.f2[itable] + (b * b * b - b) * tb.f2[itable + 1]) * tb.deltasq6; + } +} diff --git a/src/fix_wall_table.h b/src/fix_wall_table.h new file mode 100644 index 0000000000..368c2e3cd1 --- /dev/null +++ b/src/fix_wall_table.h @@ -0,0 +1,67 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS +// clang-format off +FixStyle(wall/table,FixWallTable); +// clang-format on +#else + +#ifndef LMP_FIX_WALL_TABLE_H +#define LMP_FIX_WALL_TABLE_H + +#include "fix_wall.h" + +namespace LAMMPS_NS { + +class FixWallTable : public FixWall { + public: + FixWallTable(class LAMMPS *, int, char **); + ~FixWallTable() override; + + void post_constructor() override; + void precompute(int) override; + void wall_particle(int, int, double) override; + + protected: + double offset[6]; + + int tabstyle, tablength; + struct Table { + int ninput, fpflag; + double fplo, fphi; + double lo, hi; + double *rfile, *efile, *ffile; + double *e2file, *f2file; + double delta, invdelta, deltasq6; + double *r, *e, *de, *f, *df, *e2, *f2; + }; + Table *tables; + + void null_table(Table &); + void free_table(Table &); + void read_table(Table &, const std::string &, const std::string &); + void bcast_table(Table &); + void spline_table(Table &); + void compute_table(Table &); + + void param_extract(Table &, char *); + void spline(double *, double *, int, double, double, double *); + double splint(double *, double *, double *, int, double); + void uf_lookup(int, double, double &, double &); +}; + +} // namespace LAMMPS_NS + +#endif +#endif