Merge pull request lammps#4048 from akohlmey/lepton_updates

Update LEPTON styles for more flexibility

doc/src/angle_lepton.rst

@@ -11,14 +11,24 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   angle_style lepton
+   angle_style style args
+
+* style = *lepton*
+* args = optional arguments
+
+.. parsed-literal::
+
+   args = *auto_offset* or *no_offset*
+     *auto_offset* = offset the potential energy so that the value at theta0 is 0.0 (default)
+     *no_offset* = do not offset the potential energy
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
    angle_style lepton
+   angle_style lepton no_offset
 
    angle_coeff  1  120.0  "k*theta^2; k=250.0"
    angle_coeff  2   90.0  "k2*theta^2 + k3*theta^3 + k4*theta^4; k2=300.0; k3=-100.0; k4=50.0"
@@ -41,6 +51,13 @@ angle coefficient.  For example `"200.0*theta^2"` represents a
 
    U_{angle,i} = K (\theta_i - \theta_0)^2 = K \theta^2 \qquad \theta = \theta_i - \theta_0
 
+.. versionchanged:: TBD
+
+By default the potential energy U is shifted so that the value U is 0.0
+for $theta = theta_0$.  This is equivalent to using the optional keyword
+*auto_offset*.  When using the keyword *no_offset* instead, the
+potential energy is not shifted.
+
 The `Lepton library <https://simtk.org/projects/lepton>`_, that the
 *lepton* angle style interfaces with, evaluates this expression string
 at run time to compute the pairwise energy.  It also creates an

doc/src/bond_lepton.rst

@@ -11,14 +11,24 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   bond_style lepton
+   bond_style style args
+
+* style = *lepton*
+* args = optional arguments
+
+.. parsed-literal::
+
+   args = *auto_offset* or *no_offset*
+     *auto_offset* = offset the potential energy so that the value at r0 is 0.0 (default)
+     *no_offset* = do not offset the potential energy
 
 Examples
 """"""""
 
 .. code-block:: LAMMPS
 
    bond_style lepton
+   bond_style lepton no_offset
 
    bond_coeff  1  1.5 "k*r^2; k=250.0"
    bond_coeff  2  1.1 "k2*r^2 + k3*r^3 + k4*r^4; k2=300.0; k3=-100.0; k4=50.0"
@@ -40,6 +50,13 @@ constant *K* of 200.0 energy units:
 
    U_{bond,i} = K (r_i - r_0)^2 = K r^2 \qquad r = r_i - r_0
 
+.. versionchanged:: TBD
+
+By default the potential energy U is shifted so that he value U is 0.0
+for $r = r_0$.  This is equivalent to using the optional keyword
+*auto_offset*.  When using the keyword *no_offset* instead, the
+potential energy is not shifted.
+
 The `Lepton library <https://simtk.org/projects/lepton>`_, that the
 *lepton* bond style interfaces with, evaluates this expression string at
 run time to compute the pairwise energy.  It also creates an analytical

doc/src/pair_lepton.rst

@@ -72,7 +72,7 @@ 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
+"qi" and "qj", respectively.  With pair style *lepton/sphere* one may
 instead reference the radii of the two atoms of the pair with "radi" and
 "radj", respectively; this is half of the diameter that can be set in
 :doc:`data files <read_data>` or the :doc:`set command <set>`.
@@ -166,8 +166,8 @@ 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 styles *lepton/coul* and *lepton/sphere* do *not*.
+option for shifting the potential energy of the pair interaction so that
+it is 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.

src/LEPTON/angle_lepton.cpp

@@ -44,6 +44,7 @@ AngleLepton::AngleLepton(LAMMPS *_lmp) :
 {
   writedata = 1;
   reinitflag = 0;
+  auto_offset = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -90,10 +91,21 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
 {
   std::vector<Lepton::CompiledExpression> angleforce;
   std::vector<Lepton::CompiledExpression> anglepot;
-  for (const auto &expr : expressions) {
-    auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
-    angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
-    if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
+  std::vector<bool> has_ref;
+  try {
+    for (const auto &expr : expressions) {
+      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
+      angleforce.emplace_back(parsed.differentiate("theta").createCompiledExpression());
+      has_ref.push_back(true);
+      try {
+        angleforce.back().getVariableReference("theta");
+      } catch (Lepton::Exception &) {
+        has_ref.back() = false;
+      }
+      if (EFLAG) anglepot.emplace_back(parsed.createCompiledExpression());
+    }
+  } catch (std::exception &e) {
+    error->all(FLERR, e.what());
   }
 
   const double *const *const x = atom->x;
@@ -142,8 +154,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
 
     const double dtheta = acos(c) - theta0[type];
     const int idx = type2expression[type];
-    angleforce[idx].getVariableReference("theta") = dtheta;
-
+    if (has_ref[idx]) angleforce[idx].getVariableReference("theta") = dtheta;
     const double a = -angleforce[idx].evaluate() * s;
     const double a11 = a * c / rsq1;
     const double a12 = -a / (r1 * r2);
@@ -179,7 +190,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void AngleLepton::eval()
 
     double eangle = 0.0;
     if (EFLAG) {
-      anglepot[idx].getVariableReference("theta") = dtheta;
+      try {
+        anglepot[idx].getVariableReference("theta") = dtheta;
+      } catch (Lepton::Exception &) {
+        ;    // ignore -> constant force
+      }
       eangle = anglepot[idx].evaluate() - offset[type];
     }
     if (EVFLAG)
@@ -202,6 +217,24 @@ void AngleLepton::allocate()
   for (int i = 1; i < np1; i++) setflag[i] = 0;
 }
 
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void AngleLepton::settings(int narg, char **arg)
+{
+  auto_offset = 1;
+  if (narg > 0) {
+    if (strcmp(arg[0],"auto_offset") == 0) {
+      auto_offset = 1;
+    } else if (strcmp(arg[0],"no_offset") == 0) {
+      auto_offset = 0;
+    } else {
+      error->all(FLERR, "Unknown angle style lepton setting {}", arg[0]);
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
@@ -224,9 +257,20 @@ void AngleLepton::coeff(int narg, char **arg)
     auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
     auto anglepot = parsed.createCompiledExpression();
     auto angleforce = parsed.differentiate("theta").createCompiledExpression();
-    anglepot.getVariableReference("theta") = 0.0;
-    angleforce.getVariableReference("theta") = 0.0;
-    offset_one = anglepot.evaluate();
+    try {
+      anglepot.getVariableReference("theta") = 0.0;
+    } catch (Lepton::Exception &) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Lepton potential expression {} does not depend on 'theta'", exp_one);
+    }
+    try {
+      angleforce.getVariableReference("theta") = 0.0;
+    } catch (Lepton::Exception &) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Force from Lepton expression {} does not depend on 'theta'",
+                       exp_one);
+    }
+    if (auto_offset) offset_one = anglepot.evaluate();
     angleforce.evaluate();
   } catch (std::exception &e) {
     error->all(FLERR, e.what());
@@ -284,6 +328,7 @@ void AngleLepton::write_restart(FILE *fp)
     fwrite(&n, sizeof(int), 1, fp);
     fwrite(exp.c_str(), sizeof(char), n, fp);
   }
+  fwrite(&auto_offset, sizeof(int), 1, fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -323,6 +368,9 @@ void AngleLepton::read_restart(FILE *fp)
     expressions.emplace_back(buf);
   }
 
+  if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
+  MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
+
   delete[] buf;
 }
 
@@ -363,7 +411,11 @@ double AngleLepton::single(int type, int i1, int i2, int i3)
   const auto &expr = expressions[type2expression[type]];
   auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
   auto anglepot = parsed.createCompiledExpression();
-  anglepot.getVariableReference("theta") = dtheta;
+  try {
+    anglepot.getVariableReference("theta") = dtheta;
+  } catch (Lepton::Exception &) {
+    ;    // ignore -> constant potential
+  }
   return anglepot.evaluate() - offset[type];
 }
 

src/LEPTON/angle_lepton.h

@@ -29,6 +29,7 @@ class AngleLepton : public Angle {
   AngleLepton(class LAMMPS *);
   ~AngleLepton() override;
   void compute(int, int) override;
+  void settings(int, char **) override;
   void coeff(int, char **) override;
   double equilibrium_angle(int) override;
   void write_restart(FILE *) override;
@@ -42,6 +43,7 @@ class AngleLepton : public Angle {
   double *theta0;
   int *type2expression;
   double *offset;
+  int auto_offset;
 
   virtual void allocate();
 

src/LEPTON/bond_lepton.cpp

@@ -37,6 +37,7 @@ BondLepton::BondLepton(LAMMPS *_lmp) :
 {
   writedata = 1;
   reinitflag = 0;
+  auto_offset = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -82,10 +83,17 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
 {
   std::vector<Lepton::CompiledExpression> bondforce;
   std::vector<Lepton::CompiledExpression> bondpot;
+  std::vector<bool> has_ref;
   try {
     for (const auto &expr : expressions) {
       auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
       bondforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
+      has_ref.push_back(true);
+      try {
+        bondforce.back().getVariableReference("r");
+      } catch (Lepton::Exception &) {
+        has_ref.back() = false;
+      }
       if (EFLAG) bondpot.emplace_back(parsed.createCompiledExpression());
     }
   } catch (std::exception &e) {
@@ -116,7 +124,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
 
     double fbond = 0.0;
     if (r > 0.0) {
-      bondforce[idx].getVariableReference("r") = dr;
+      if (has_ref[idx]) bondforce[idx].getVariableReference("r") = dr;
       fbond = -bondforce[idx].evaluate() / r;
     }
 
@@ -136,7 +144,11 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
 
     double ebond = 0.0;
     if (EFLAG) {
-      bondpot[idx].getVariableReference("r") = dr;
+      try {
+        bondpot[idx].getVariableReference("r") = dr;
+      } catch (Lepton::Exception &) {
+        ;    // ignore -> constant potential
+      }
       ebond = bondpot[idx].evaluate() - offset[type];
     }
     if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz);
@@ -157,6 +169,24 @@ void BondLepton::allocate()
   for (int i = 1; i < np1; i++) setflag[i] = 0;
 }
 
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void BondLepton::settings(int narg, char **arg)
+{
+  auto_offset = 1;
+  if (narg > 0) {
+    if (strcmp(arg[0],"auto_offset") == 0) {
+      auto_offset = 1;
+    } else if (strcmp(arg[0],"no_offset") == 0) {
+      auto_offset = 0;
+    } else {
+      error->all(FLERR, "Unknown bond style lepton setting {}", arg[0]);
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    set coeffs for one or more types
 ------------------------------------------------------------------------- */
@@ -179,9 +209,19 @@ void BondLepton::coeff(int narg, char **arg)
     auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
     auto bondpot = parsed.createCompiledExpression();
     auto bondforce = parsed.differentiate("r").createCompiledExpression();
-    bondpot.getVariableReference("r") = 0.0;
-    bondforce.getVariableReference("r") = 0.0;
-    offset_one = bondpot.evaluate();
+    try {
+      bondpot.getVariableReference("r") = 0.0;
+    } catch (Lepton::Exception &e) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Lepton potential expression {} does not depend on 'r'", exp_one);
+    }
+    try {
+      bondforce.getVariableReference("r") = 0.0;
+    } catch (Lepton::Exception &e) {
+      if (comm->me == 0)
+        error->warning(FLERR, "Force from Lepton expression {} does not depend on 'r'", exp_one);
+    }
+    if (auto_offset) offset_one = bondpot.evaluate();
     bondforce.evaluate();
   } catch (std::exception &e) {
     error->all(FLERR, e.what());
@@ -239,6 +279,7 @@ void BondLepton::write_restart(FILE *fp)
     fwrite(&n, sizeof(int), 1, fp);
     fwrite(exp.c_str(), sizeof(char), n, fp);
   }
+  fwrite(&auto_offset, sizeof(int), 1, fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -278,6 +319,9 @@ void BondLepton::read_restart(FILE *fp)
     expressions.emplace_back(buf);
   }
 
+  if (comm->me == 0) utils::sfread(FLERR, &auto_offset, sizeof(int), 1, fp, nullptr, error);
+  MPI_Bcast(&auto_offset, 1, MPI_INT, 0, world);
+
   delete[] buf;
 }
 
@@ -302,8 +346,12 @@ double BondLepton::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
   auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
   auto bondpot = parsed.createCompiledExpression();
   auto bondforce = parsed.differentiate("r").createCompiledExpression();
-  bondforce.getVariableReference("r") = dr;
-  bondpot.getVariableReference("r") = dr;
+  try {
+    bondpot.getVariableReference("r") = dr;
+    bondforce.getVariableReference("r") = dr;
+  } catch (Lepton::Exception &) {
+    ;    // ignore -> constant potential or force
+  }
 
   // force and energy
 

src/LEPTON/bond_lepton.h

@@ -29,6 +29,7 @@ class BondLepton : public Bond {
   BondLepton(class LAMMPS *);
   ~BondLepton() override;
   void compute(int, int) override;
+  void settings(int, char **) override;
   void coeff(int, char **) override;
   double equilibrium_distance(int) override;
   void write_restart(FILE *) override;
@@ -42,6 +43,7 @@ class BondLepton : public Bond {
   double *r0;
   int *type2expression;
   double *offset;
+  int auto_offset;
 
   virtual void allocate();