add pair style lepton/sphere

doc/src/Commands_pair.rst

@@ -137,6 +137,7 @@ OPT.
    * :doc:`lennard/mdf <pair_mdf>`
    * :doc:`lepton (o) <pair_lepton>`
    * :doc:`lepton/coul (o) <pair_lepton>`
+   * :doc:`lepton/sphere <pair_lepton>`
    * :doc:`line/lj <pair_line_lj>`
    * :doc:`lj/charmm/coul/charmm (giko) <pair_charmm>`
    * :doc:`lj/charmm/coul/charmm/implicit (ko) <pair_charmm>`

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
+.. versionchanged:: TBD
-between particles which depend solely on the distance and have a cutoff.
+
-The potential function must be provided as an expression string using
+   added pair style lepton/sphere
-"r" as the distance variable.  With pair style *lepton/coul* one may
+
-additionally reference the charges of the two atoms of the pair with
+Pair styles *lepton*, *lepton/coul*, *lepton/sphere* compute pairwise
-"qi" and "qj", respectively.  Note that further constants in the
+interactions between particles which depend on the distance and have a
-expression can be defined in the same string as additional expressions
+cutoff.  The potential function must be provided as an expression string
-separated by semi-colons as shown in the examples above.
+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 <https://simtk.org/projects/lepton>`_, 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 <special_bonds>`
+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.
+<special_bonds>` settings apply in case the interacting pair is also
-The potential energy will *only* be added to the "evdwl" property.
+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 <special_bonds>`
+For pair style *lepton/coul* only the "coul" values of the
-settings apply in case the interacting pair is also connected with a bond.
+:doc:`special_bonds <special_bonds>` settings apply in case the
-The potential energy will *only* be added to the "ecoul" property.
+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 <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,
+Pair styles *lepton*, *lepton/coul*, and *lepton/sphere* do not support
-expressions for *all* I,J pairs must be specified explicitly.
+mixing.  Thus, expressions for *all* I,J pairs must be specified
+explicitly.
 
 Only pair style *lepton* supports the :doc:`pair_modify shift <pair_modify>`
 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 <pair_modify>` options are not relevant for
 the these pair styles.
@@ -141,7 +170,7 @@ These pair styles do not support the :doc:`pair_modify tail
 <pair_modify>` 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
 <restart>`, 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
 <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 <atom_style>`.
+
+Pair style *lepton/sphere* requires that atom atoms have a radius
+property, e.g. via :doc:`atom_style sphere <atom_style>`.
+
 Related commands
 """"""""""""""""
 

doc/src/pair_style.rst

@@ -215,6 +215,7 @@ accelerated styles exist.
 * :doc:`lennard/mdf <pair_mdf>` - LJ potential in A/B form with a taper function
 * :doc:`lepton <pair_lepton>` - pair potential from evaluating a string
 * :doc:`lepton/coul <pair_lepton>` - pair potential from evaluating a string with support for charges
+* :doc:`lepton/sphere <pair_lepton>` - pair potential from evaluating a string with support for radii
 * :doc:`line/lj <pair_line_lj>` - LJ potential between line segments
 * :doc:`list <pair_list>` - potential between pairs of atoms explicitly listed in an input file
 * :doc:`lj/charmm/coul/charmm <pair_charmm>` - CHARMM potential with cutoff Coulomb

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

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 <cmath>
+
+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 <int EVFLAG, int EFLAG, int NEWTON_PAIR> 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<Lepton::CompiledExpression> pairforce;
+  std::vector<Lepton::CompiledExpression> pairpot;
+  std::vector<std::pair<bool, bool>> 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;
+}

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 <int EVFLAG, int EFLAG, int NEWTON_PAIR> void eval();
+};
+}    // namespace LAMMPS_NS
+#endif
+#endif