diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 5328d1163f..d0f83be3ba 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -137,6 +137,7 @@ OPT. * :doc:`lennard/mdf ` * :doc:`lepton (o) ` * :doc:`lepton/coul (o) ` + * :doc:`lepton/sphere ` * :doc:`line/lj ` * :doc:`lj/charmm/coul/charmm (giko) ` * :doc:`lj/charmm/coul/charmm/implicit (ko) ` diff --git a/doc/src/pair_lepton.rst b/doc/src/pair_lepton.rst index dbe3d74588..b83366d752 100644 --- a/doc/src/pair_lepton.rst +++ b/doc/src/pair_lepton.rst @@ -2,6 +2,7 @@ .. index:: pair_style lepton/omp .. index:: pair_style lepton/coul .. index:: pair_style lepton/coul/omp +.. index:: pair_style lepton/sphere pair_style lepton command ========================= @@ -15,7 +16,7 @@ Syntax pair_style style args -* style = *lepton* or *lepton/coul* +* style = *lepton* or *lepton/coul* or *lepton/sphere* * args = list of arguments for a particular style .. parsed-literal:: @@ -26,6 +27,8 @@ Syntax cutoff = global cutoff for the interactions (distance units) zero or more keywords may be appended keyword = *ewald* or *pppm* or *msm* or *dispersion* or *tip4p* + *lepton/sphere* args = cutoff + cutoff = global cutoff for the interactions (distance units) Examples """""""" @@ -48,19 +51,29 @@ Examples kspace_style pppm 1.0e-4 pair_coeff 1 1 "qi*qj/r*erfc(alpha*r); alpha=1.067" + pair_style lepton/sphere 2.5 + pair_coeff 1 * "k*((r-r0)^2*step(r0-r)); k=200; r0=radi+radj" + pair_coeff 2 2 "4.0*eps*((sig/r)^12 - (sig/r)^6); eps=1.0; sig=sqrt(radi*radj)" + Description """"""""""" .. versionadded:: 8Feb2023 -Pair styles *lepton* and *lepton/coul* compute pairwise interactions -between particles which depend solely on the distance and have a cutoff. -The potential function must be provided as an expression string using -"r" as the distance variable. With pair style *lepton/coul* one may -additionally reference the charges of the two atoms of the pair with -"qi" and "qj", respectively. Note that further constants in the -expression can be defined in the same string as additional expressions -separated by semi-colons as shown in the examples above. +.. versionchanged:: TBD + + added pair style lepton/sphere + +Pair styles *lepton*, *lepton/coul*, *lepton/sphere* compute pairwise +interactions between particles which depend on the distance and have a +cutoff. The potential function must be provided as an expression string +using "r" as the distance variable. With pair style *lepton/coul* one +may additionally reference the charges of the two atoms of the pair with +"qi" and "qj", respectively. With pair style *lepton/coul* one may +instead reference the radii of the two atoms of the pair with "radi" and +"radj", respectively. Note that further constants in the expression can +be defined in the same string as additional expressions separated by +semi-colons as shown in the examples above. The expression `"200.0*(r-1.5)^2"` represents a harmonic potential around the pairwise distance :math:`r_0` of 1.5 distance units and a @@ -76,6 +89,14 @@ The expression `"qi*qj/r"` represents a regular Coulombic potential with cutoff: U_{ij} = \frac{C q_i q_j}{\epsilon r} \qquad r < r_c +The expression `"200.0*(r-(radi+radj)^2"` represents a harmonic potential +that has the equilibrium distance chosen so that the radii of the two +atoms touch: + +.. math:: + + U_{ij} = K (r-(r_i+r_j))^2 + The `Lepton library `_, that the *lepton* pair style interfaces with, evaluates this expression string at run time to compute the pairwise energy. It also creates an analytical @@ -97,13 +118,20 @@ More on valid Lepton expressions below. The last coefficient is optional; it allows to set the cutoff for a pair of atom types to a different value than the global cutoff. -For pair style *lepton* only the "lj" value of the :doc:`special_bonds ` -settings apply in case the interacting pair is also connected with a bond. -The potential energy will *only* be added to the "evdwl" property. +For pair style *lepton* only the "lj" values of the :doc:`special_bonds +` settings apply in case the interacting pair is also +connected with a bond. The potential energy will *only* be added to the +"evdwl" property. -For pair style *lepton/coul* only the "coul" value of the :doc:`special_bonds ` -settings apply in case the interacting pair is also connected with a bond. -The potential energy will *only* be added to the "ecoul" property. +For pair style *lepton/coul* only the "coul" values of the +:doc:`special_bonds ` settings apply in case the +interacting pair is also connected with a bond. The potential energy +will *only* be added to the "ecoul" property. + +For pair style *lepton/sphere* only the "lj" values of the +:doc:`special_bonds ` settings apply in case the +interacting pair is also connected with a bond. The potential energy +will *only* be added to the "evdwl" property. In addition to the functions listed below, both pair styles support in addition a custom "zbl(zi,zj,r)" function which computes the @@ -127,12 +155,13 @@ above. Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" -Pair styles *lepton* and *lepton/coul* do not support mixing. Thus, -expressions for *all* I,J pairs must be specified explicitly. +Pair styles *lepton*, *lepton/coul*, and *lepton/sphere* do not support +mixing. Thus, expressions for *all* I,J pairs must be specified +explicitly. Only pair style *lepton* supports the :doc:`pair_modify shift ` option for shifting the energy of the pair interaction so that it is -0 at the cutoff, pair style *lepton/coul* does *not*. +0 at the cutoff, pair styles *lepton/coul* and *lepton/sphere* do *not*. The :doc:`pair_modify table ` options are not relevant for the these pair styles. @@ -141,7 +170,7 @@ These pair styles do not support the :doc:`pair_modify tail ` option for adding long-range tail corrections to energy and pressure. -These pair styles write its information to :doc:`binary restart files +These pair styles write their information to :doc:`binary restart files `, so pair_style and pair_coeff commands do not need to be specified in an input script that reads a restart file. @@ -158,6 +187,12 @@ These pair styles are part of the LEPTON package and only enabled if LAMMPS was built with this package. See the :doc:`Build package ` page for more info. +Pair style *lepton/coul* requires that atom atoms have a charge +property, e.g. via :doc:`atom_style charge `. + +Pair style *lepton/sphere* requires that atom atoms have a radius +property, e.g. via :doc:`atom_style sphere `. + Related commands """""""""""""""" diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index a0fb6eeefd..96488ee154 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -215,6 +215,7 @@ accelerated styles exist. * :doc:`lennard/mdf ` - LJ potential in A/B form with a taper function * :doc:`lepton ` - pair potential from evaluating a string * :doc:`lepton/coul ` - pair potential from evaluating a string with support for charges +* :doc:`lepton/sphere ` - pair potential from evaluating a string with support for radii * :doc:`line/lj ` - LJ potential between line segments * :doc:`list ` - potential between pairs of atoms explicitly listed in an input file * :doc:`lj/charmm/coul/charmm ` - CHARMM potential with cutoff Coulomb diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c2ecce84d2..996b7440f7 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -2973,11 +2973,13 @@ qx qy qz Rackers +radi radialscreened radialscreenedspin radialspin radian radians +radj Rafferty rahman Rahman diff --git a/src/LEPTON/pair_lepton_sphere.cpp b/src/LEPTON/pair_lepton_sphere.cpp new file mode 100644 index 0000000000..f29fd0b48a --- /dev/null +++ b/src/LEPTON/pair_lepton_sphere.cpp @@ -0,0 +1,231 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "pair_lepton_sphere.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "update.h" + +#include "Lepton.h" +#include "lepton_utils.h" +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +void PairLeptonSphere::compute(int eflag, int vflag) +{ + ev_init(eflag, vflag); + if (evflag) { + if (eflag) { + if (force->newton_pair) + eval<1, 1, 1>(); + else + eval<1, 1, 0>(); + } else { + if (force->newton_pair) + eval<1, 0, 1>(); + else + eval<1, 0, 0>(); + } + } else { + if (force->newton_pair) + eval<0, 0, 1>(); + else + eval<0, 0, 0>(); + } + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +template void PairLeptonSphere::eval() +{ + const double *const *const x = atom->x; + double *const *const f = atom->f; + const double *const radius = atom->radius; + const int *const type = atom->type; + const int nlocal = atom->nlocal; + const double *const special_lj = force->special_lj; + + const int inum = list->inum; + const int *const ilist = list->ilist; + const int *const numneigh = list->numneigh; + const int *const *const firstneigh = list->firstneigh; + double fxtmp, fytmp, fztmp; + + std::vector pairforce; + std::vector pairpot; + std::vector> have_rad; + try { + for (const auto &expr : expressions) { + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); + pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression()); + if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression()); + pairforce.back().getVariableReference("r"); + have_rad.emplace_back(true, true); + + // check if there are references to charges + try { + pairforce.back().getVariableReference("radi"); + } catch (std::exception &) { + have_rad.back().first = false; + } + try { + pairforce.back().getVariableReference("radj"); + } catch (std::exception &) { + have_rad.back().second = false; + } + } + } catch (std::exception &e) { + error->all(FLERR, e.what()); + } + + // loop over neighbors of my atoms + + for (int ii = 0; ii < inum; ii++) { + const int i = ilist[ii]; + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int *jlist = firstneigh[i]; + const int jnum = numneigh[i]; + fxtmp = fytmp = fztmp = 0.0; + + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + const double factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; + const int jtype = type[j]; + + const double delx = xtmp - x[j][0]; + const double dely = ytmp - x[j][1]; + const double delz = ztmp - x[j][2]; + const double rsq = delx * delx + dely * dely + delz * delz; + + if (rsq < cutsq[itype][jtype]) { + const double r = sqrt(rsq); + const int idx = type2expression[itype][jtype]; + pairforce[idx].getVariableReference("r") = r; + if (have_rad[idx].first) pairforce[idx].getVariableReference("radi") = radius[i]; + if (have_rad[idx].second) pairforce[idx].getVariableReference("radj") = radius[j]; + const double fpair = -pairforce[idx].evaluate() / r * factor_lj; + + fxtmp += delx * fpair; + fytmp += dely * fpair; + fztmp += delz * fpair; + if (NEWTON_PAIR || (j < nlocal)) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + } + + double evdwl = 0.0; + if (EFLAG) { + pairpot[idx].getVariableReference("r") = r; + if (have_rad[idx].first) pairpot[idx].getVariableReference("radi") = radius[i]; + if (have_rad[idx].second) pairpot[idx].getVariableReference("radj") = radius[j]; + evdwl = pairpot[idx].evaluate(); + evdwl *= factor_lj; + } + + if (EVFLAG) ev_tally(i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz); + } + } + f[i][0] += fxtmp; + f[i][1] += fytmp; + f[i][2] += fztmp; + } +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairLeptonSphere::settings(int narg, char **arg) +{ + if (narg != 1) + error->all(FLERR, "Incorrect number of arguments for pair_style lepton/sphere command"); + cut_global = utils::numeric(FLERR, arg[0], false, lmp); +} + +/* ---------------------------------------------------------------------- */ + +void PairLeptonSphere::init_style() +{ + if (!atom->radius_flag) + error->all(FLERR, "Pair style lepton/sphere requires atom attribute radius"); + if (offset_flag) error->all(FLERR, "Pair style lepton/sphere does not suport pair_modify shift"); + neighbor->add_request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairLeptonSphere::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global, sizeof(double), 1, fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairLeptonSphere::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); } + MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world); +} + +/* ---------------------------------------------------------------------- */ + +double PairLeptonSphere::single(int i, int j, int itype, int jtype, double rsq, + double /* factor_coul */, double factor_lj, double &fforce) +{ + auto expr = expressions[type2expression[itype][jtype]]; + auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions); + auto pairpot = parsed.createCompiledExpression(); + auto pairforce = parsed.differentiate("r").createCompiledExpression(); + + const double r = sqrt(rsq); + pairpot.getVariableReference("r") = r; + pairforce.getVariableReference("r") = r; + try { + pairpot.getVariableReference("radi") = atom->radius[i]; + pairforce.getVariableReference("radi") = atom->radius[i]; + } catch (std::exception &) { + /* ignore */ + } + try { + pairpot.getVariableReference("radj") = atom->radius[j]; + pairforce.getVariableReference("radj") = atom->radius[j]; + } catch (std::exception &) { + /* ignore */ + } + + fforce = -pairforce.evaluate() / r * factor_lj; + return pairpot.evaluate() * factor_lj; +} diff --git a/src/LEPTON/pair_lepton_sphere.h b/src/LEPTON/pair_lepton_sphere.h new file mode 100644 index 0000000000..ab586a309b --- /dev/null +++ b/src/LEPTON/pair_lepton_sphere.h @@ -0,0 +1,44 @@ +/* -*- 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 PAIR_CLASS +// clang-format off +PairStyle(lepton/sphere,PairLeptonSphere); +// clang-format on +#else + +#ifndef LMP_PAIR_LEPTON_SPHERE_H +#define LMP_PAIR_LEPTON_SPHERE_H + +#include "pair_lepton.h" + +namespace LAMMPS_NS { + +class PairLeptonSphere : public PairLepton { + public: + PairLeptonSphere(class LAMMPS *_lmp) : PairLepton(_lmp){}; + + void compute(int, int) override; + void settings(int, char **) override; + void init_style() override; + void write_restart_settings(FILE *) override; + void read_restart_settings(FILE *) override; + double single(int, int, int, int, double, double, double, double &) override; + + private: + template void eval(); +}; +} // namespace LAMMPS_NS +#endif +#endif