Merge pull request #3590 from akohlmey/lepton-zbl

Add a custom zbl() function to lepton pair styles

doc/src/pair_lepton.rst

@@ -37,6 +37,8 @@ Examples
    pair_coeff  * *  "k*((r-r0)^2*step(r0-r)); k=200; r0=1.5" 2.0
    pair_coeff  1 2  "4.0*eps*((sig/r)^12 - (sig/r)^6);eps=1.0;sig=1.0" 1.12246204830937
    pair_coeff  2 2  "eps*(2.0*(sig/r)^9 - 3.0*(sig/r)^6);eps=1.0;sig=1.0"
+   pair_coeff  1 3  "zbl(13,6,r)"
+   pair_coeff  3 3  "(1.0-switch)*zbl(6,6,r)-switch*4.0*eps*((sig/r)^6);switch=0.5*(tanh(10.0*(r-sig))+1.0);eps=0.05;sig=3.20723"
 
    pair_style lepton/coul 2.5
    pair_coeff 1 1 "qi*qj/r" 4.0
@@ -103,6 +105,15 @@ For pair style *lepton/coul* only the "coul" value of the :doc:`special_bonds <s
 settings apply in case the interacting pair is also connected with a bond.
 The potential energy will *only* be added to the "ecoul" property.
 
+In addition to the functions listed below, both pair styles support in
+addition a custom "zbl(zi,zj,r)" function which computes the
+Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion for
+describing high-energy collisions between atoms.  For details of the
+function please see the documentation for :doc:`pair style zbl
+<pair_zbl>`. The arguments of the function are the atomic numbers of
+atom i (zi), atom j (zj) and the distance r.  Please see the examples
+above.
+
 ----------
 
 .. include:: lepton_expression.rst

src/LEPTON/lepton_utils.cpp

@@ -20,14 +20,73 @@
 #include "error.h"
 #include "input.h"
 #include "lammps.h"
+#include "pair_zbl_const.h"
 #include "variable.h"
 
 #include "fmt/args.h"
 
 #include <cctype>
+#include <cmath>
 #include <exception>
 #include <unordered_set>
 
+using namespace LAMMPS_NS;
+using namespace PairZBLConstants;
+
+namespace Lepton {
+class DerivativeException : public std::exception {
+  std::string message;
+
+ public:
+  // remove unused default constructor
+  DerivativeException() = delete;
+
+  explicit DerivativeException(int deg, const std::string &fn, const std::string &vn)
+  {
+    message = fmt::format("Order {} derivative of function {} in {} is not supported", deg, fn, vn);
+  }
+  const char *what() const noexcept override { return message.c_str(); }
+};
+}    // namespace Lepton
+
+double Lepton::ZBLFunction::evaluate(const double *args) const
+{
+  const double zi = args[0];
+  const double zj = args[1];
+  const double r = args[2];
+
+  const double rbya = r * (pow(zi, pzbl) + pow(zj, pzbl)) / (a0 * angstrom);
+  return zi * zj * qqr2e * qelectron * qelectron / r *
+      (c4 * exp(-d4 * rbya) + c3 * exp(-d3 * rbya) + c2 * exp(-d2 * rbya) + c1 * exp(-d1 * rbya));
+}
+
+double Lepton::ZBLFunction::evaluateDerivative(const double *args, const int *order) const
+{
+  if (order[0] > 0)
+    throw DerivativeException(order[0], "zbl()", "'zi'");
+  if (order[1] > 0)
+    throw DerivativeException(order[0], "zbl()", "'zj'");
+  if (order[2] > 1)
+    throw DerivativeException(order[0], "zbl()", "'r'");
+
+  if (order[2] == 1) {
+    const double zi = args[0];
+    const double zj = args[1];
+    const double r = args[2];
+
+    const double ainv = (pow(zi, pzbl) + pow(zj, pzbl)) / (a0 * angstrom);
+    const double e1 = exp(-d1 * ainv * r);
+    const double e2 = exp(-d2 * ainv * r);
+    const double e3 = exp(-d3 * ainv * r);
+    const double e4 = exp(-d4 * ainv * r);
+
+    const double sum1 = c1 * e1 + c2 * e2 + c3 * e3 + c4 * e4;
+    const double sum2 = ainv * (-c1 * d1 * e1 - c2 * d2 * e2 - c3 * d3 * e3 - c4 * d4 * e4);
+    return (zi * zj * qqr2e * qelectron * qelectron) * (sum2 - sum1 / r) / r;
+  }
+  return 0.0;
+}
+
 namespace LeptonUtils {
 class VariableException : public std::exception {
   std::string message;

src/LEPTON/lepton_utils.h

@@ -11,6 +11,8 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include "Lepton.h"
+
 #include <string>
 
 // forward declarations
@@ -19,6 +21,25 @@ namespace LAMMPS_NS {
 class LAMMPS;
 }
 
+// custom function zbl(zi,zj,r)
+
+namespace Lepton {
+class ZBLFunction : public CustomFunction {
+ public:
+  ZBLFunction(double _qqr2e, double _angstrom, double _qelectron) :
+      qqr2e(_qqr2e), angstrom(_angstrom), qelectron(_qelectron){};
+  ZBLFunction() = delete;
+
+  int getNumArguments() const override { return 3; }
+  CustomFunction *clone() const override { return new ZBLFunction(qqr2e, angstrom, qelectron); }
+  double evaluate(const double *) const override;
+  double evaluateDerivative(const double *, const int *) const override;
+
+ private:
+  double qqr2e, angstrom, qelectron;
+};
+}    // namespace Lepton
+
 // utility functions and classes
 
 namespace LeptonUtils {

src/LEPTON/pair_lepton.cpp

@@ -28,6 +28,7 @@
 #include "Lepton.h"
 #include "lepton_utils.h"
 #include <cmath>
+#include <map>
 
 using namespace LAMMPS_NS;
 
@@ -43,12 +44,15 @@ PairLepton::PairLepton(LAMMPS *lmp) :
   reinitflag = 0;
   cut_global = 0.0;
   centroidstressflag = CENTROID_SAME;
+
+  functions["zbl"] = new Lepton::ZBLFunction(force->qqr2e, force->angstrom, force->qelectron);
 }
 
 /* ---------------------------------------------------------------------- */
 
 PairLepton::~PairLepton()
 {
+  for (auto &f : functions) delete f.second;
   if (allocated) {
     memory->destroy(cut);
     memory->destroy(cutsq);
@@ -103,7 +107,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
   std::vector<Lepton::CompiledExpression> pairpot;
   try {
     for (const auto &expr : expressions) {
-      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
+      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
       pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
       pairforce.back().getVariableReference("r");
       if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
@@ -222,7 +226,7 @@ void PairLepton::coeff(int narg, char **arg)
   // check if the expression can be parsed and evaluated without error
   auto exp_one = LeptonUtils::condense(arg[2]);
   try {
-    auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp));
+    auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(exp_one, lmp), functions);
     auto pairforce = parsed.differentiate("r").createCompiledExpression();
     auto pairpot = parsed.createCompiledExpression();
     pairpot.getVariableReference("r") = 1.0;
@@ -265,7 +269,7 @@ double PairLepton::init_one(int i, int j)
   if (offset_flag) {
     try {
       auto expr = LeptonUtils::substitute(expressions[type2expression[i][j]], lmp);
-      auto pairpot = Lepton::Parser::parse(expr).createCompiledExpression();
+      auto pairpot = Lepton::Parser::parse(expr, functions).createCompiledExpression();
       pairpot.getVariableReference("r") = cut[i][j];
       offset[i][j] = pairpot.evaluate();
     } catch (std::exception &) {
@@ -420,7 +424,7 @@ double PairLepton::single(int /* i */, int /* j */, int itype, int jtype, double
                           double /* factor_coul */, double factor_lj, double &fforce)
 {
   auto expr = expressions[type2expression[itype][jtype]];
-  auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
+  auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
   auto pairpot = parsed.createCompiledExpression();
   auto pairforce = parsed.differentiate("r").createCompiledExpression();
 

src/LEPTON/pair_lepton.h

@@ -22,6 +22,12 @@ PairStyle(lepton,PairLepton);
 
 #include "pair.h"
 
+#include <map>
+
+namespace Lepton {
+class CustomFunction;
+}
+
 namespace LAMMPS_NS {
 
 class PairLepton : public Pair {
@@ -42,6 +48,8 @@ class PairLepton : public Pair {
 
  protected:
   std::vector<std::string> expressions;
+  std::map<std::string, Lepton::CustomFunction *> functions;
+
   double **cut;
   int **type2expression;
   double **offset;

src/LEPTON/pair_lepton_coul.cpp

@@ -82,7 +82,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLeptonCoul::eval()
   std::vector<std::pair<bool, bool>> have_q;
   try {
     for (const auto &expr : expressions) {
-      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
+      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");
@@ -243,7 +243,7 @@ double PairLeptonCoul::single(int i, int j, int itype, int jtype, double rsq, do
                               double /* factor_lj */, double &fforce)
 {
   auto expr = expressions[type2expression[itype][jtype]];
-  auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp));
+  auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, lmp), functions);
   auto pairpot = parsed.createCompiledExpression();
   auto pairforce = parsed.differentiate("r").createCompiledExpression();
 

src/OPENMP/pair_lepton_coul_omp.cpp

@@ -104,7 +104,7 @@ void PairLeptonCoulOMP::eval(int iifrom, int iito, ThrData *const thr)
   std::vector<std::pair<bool, bool>> have_q;
   try {
     for (const auto &expr : expressions) {
-      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
+      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
       pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
       if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
       pairforce.back().getVariableReference("r");

src/OPENMP/pair_lepton_omp.cpp

@@ -98,7 +98,7 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
   std::vector<Lepton::CompiledExpression> pairpot;
   try {
     for (const auto &expr : expressions) {
-      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp));
+      auto parsed = Lepton::Parser::parse(LeptonUtils::substitute(expr, Pointers::lmp), functions);
       pairforce.emplace_back(parsed.differentiate("r").createCompiledExpression());
       if (EFLAG) pairpot.emplace_back(parsed.createCompiledExpression());
     }

unittest/force-styles/tests/mol-pair-lepton_zbl.yaml

@@ -0,0 +1,102 @@
+---
+lammps_version: 22 Dec 2022
+date_generated: Sun Jan  8 01:13:49 2023
+epsilon: 5e-10
+skip_tests:
+prerequisites: ! |
+  atom full
+  pair lepton
+pre_commands: ! |
+  variable bond_factor     index 1.0
+  variable angle_factor    index 1.0
+  variable dihedral_factor index 1.0
+post_commands: ! ""
+input_file: in.fourmol
+pair_style: lepton 8.0
+pair_coeff: ! |
+  1 1  zbl(6,6,r)
+  1 2  zbl(6,7,r)
+  1 3  zbl(6,2,r)
+  1 4  zbl(6,8,r)
+  1 5  zbl(6,8,r)
+  2 2  zbl(7,7,r)
+  2 3  zbl(7,2,r)
+  2 4  zbl(7,8,r)
+  2 5  zbl(7,8,r)
+  3 3  zbl(2,2,r)
+  3 4  zbl(2,8,r)
+  3 5  zbl(2,8,r)
+  4 4  zbl(8,8,r)
+  4 5  zbl(8,8,r)
+  5 5  zbl(8,8,r)
+extract: ! ""
+natoms: 29
+init_vdwl: 7661.0923034839425
+init_coul: 0
+init_stress: ! |2-
+   1.0415723081438835e+04  7.6654640475740289e+03  1.0439568803573100e+04  2.5596340687548091e+03  3.0124073227376857e+03  3.1819159178174647e+02
+init_forces: ! |2
+    1 -2.2542878982886191e+02 -5.1056064396745825e+01  4.0347728160298351e+02
+    2  3.4577199265194707e+02  3.3685870288327618e+02 -3.5822722442776552e+02
+    3  6.7803530899610132e+02  3.5041867763592933e+01 -3.5973378440582275e+02
+    4 -8.1435522635968925e+02  2.2040390876529949e+02 -5.9600610062358885e+02
+    5 -2.1098162371264527e+02 -3.3211428613105602e+02  1.0436615327229824e+03
+    6 -8.2893212979873752e+01  1.2663096062116030e+02  4.4419304268422735e+02
+    7 -3.4934572486178489e+01 -1.2862622365409325e+02 -7.5679463440025643e+02
+    8 -1.4121290372363833e+02 -1.1033884290619251e+02 -3.6118050221036594e+02
+    9  9.6547187735733928e+01  2.0268265922256009e+02  6.9609964660581829e+02
+   10  8.2081753362417444e+01  2.1490541800870224e+02  6.7445631928792000e+02
+   11 -2.7501706654968336e+02 -6.4133882479619604e+02 -8.2554240753015847e+02
+   12 -2.3615126531081731e+02 -3.7888094263940093e+02  4.2996552265607966e+02
+   13  9.8935659849196156e+02 -3.9189784587066578e+02  1.6979303593289488e+00
+   14 -3.2831677005156325e+02  9.7498291303146885e+01 -1.0618689865846111e+03
+   15  2.3059516739411034e+01  1.0120552702611462e+03  3.8891358018891088e+02
+   16  3.2227110665984651e+02 -3.2696296221307608e+02 -3.9821264616502327e+02
+   17 -1.9854928611614176e+02  1.0231604600726870e+02  6.4707382129491782e+02
+   18 -1.1566833968062178e+02 -5.0500912216600113e+02  1.7898146720451105e+03
+   19 -1.0355744892944338e+03 -6.8474702051102486e+02 -1.1421521343728634e+03
+   20  1.1567664608091418e+03  1.2055463642684399e+03 -6.6488505245382839e+02
+   21 -5.1058201416833657e+02 -6.2261916138266497e+02  1.8034425299414422e+03
+   22 -1.0303537305440429e+03 -2.8540172556947783e+02 -1.4390471053430481e+03
+   23  1.5447065970289084e+03  9.0276100147850616e+02 -3.6857879996771908e+02
+   24  3.7449763655418099e+02 -1.6372799599084626e+03  9.3993003918438728e+02
+   25 -1.3650956291199398e+03  1.2266826749083317e+02 -1.1043020934470092e+03
+   26  9.9428586164565161e+02  1.5187270928735547e+03  1.7246869495941340e+02
+   27  2.8070997771034445e+02 -1.7944696249727249e+03  6.4639058889384091e+02
+   28 -1.4822994833123780e+03  5.2787642883033504e+02 -9.7293790834683261e+02
+   29  1.1993244048531997e+03  1.2647703273399616e+03  3.2788417785152956e+02
+run_vdwl: 7355.7708040367215
+run_coul: 0
+run_stress: ! |2-
+   1.0014026118319458e+04  7.4146110545302845e+03  1.0129433554690771e+04  2.4014185997980703e+03  2.8618718144358500e+03  3.1547123090533302e+02
+run_forces: ! |2
+    1 -2.1963909767996165e+02 -4.7230048342286665e+01  3.9447474330057327e+02
+    2  3.3823602012753361e+02  3.3060841467248804e+02 -3.4909191972075922e+02
+    3  6.4685420887736268e+02  3.4155483339278817e+01 -3.4174127433776528e+02
+    4 -7.8636648191524102e+02  2.1116614600799022e+02 -5.7402164047826784e+02
+    5 -2.0592655691757219e+02 -3.1948425705952934e+02  1.0031851850554129e+03
+    6 -8.1885702477873082e+01  1.2476420007196047e+02  4.3622385152526408e+02
+    7 -3.5214786990923372e+01 -1.2723182110324109e+02 -7.4725777446879317e+02
+    8 -1.3797341261617603e+02 -1.0667961036533134e+02 -3.4719859501013519e+02
+    9  9.3157498402761519e+01  1.9853026083255594e+02  6.7970067058747793e+02
+   10  7.4694873724181420e+01  1.9549844828167281e+02  6.3962599563799461e+02
+   11 -2.6505567174539726e+02 -6.1810691652881007e+02 -7.8981506764583025e+02
+   12 -2.2346480012776817e+02 -3.6198854997039064e+02  4.0604938064590891e+02
+   13  9.5587946001324235e+02 -3.7492976806829023e+02  4.0412877209292042e-01
+   14 -3.1069307401641237e+02  9.5118496292191480e+01 -1.0209459166022623e+03
+   15  2.3074811498468790e+01  9.7760720729835271e+02  3.7462072751151828e+02
+   16  3.2004057723494378e+02 -3.2580441129415715e+02 -3.9361639626975176e+02
+   17 -1.9645839108898824e+02  1.0114353674483996e+02  6.4131023697143405e+02
+   18 -9.0905016804868524e+01 -4.5649801918985156e+02  1.6740121681099276e+03
+   19 -9.7189809015843355e+02 -6.4245306822849329e+02 -1.0761208024504872e+03
+   20  1.0683309313831251e+03  1.1147511786417067e+03 -6.1516522889263013e+02
+   21 -4.7823440562179218e+02 -5.8004068421284194e+02  1.6903616731086115e+03
+   22 -9.6232775232513484e+02 -2.6660807378911699e+02 -1.3469354884213956e+03
+   23  1.4443362128008671e+03  8.4138206482572116e+02 -3.4761279951083094e+02
+   24  3.5961673762879963e+02 -1.5399088881652940e+03  8.8844237731370015e+02
+   25 -1.2856259100984250e+03  1.1769150558030140e+02 -1.0419813334649532e+03
+   26  9.2973443424366496e+02  1.4263851950698765e+03  1.6174119101406066e+02
+   27  2.6191415819111842e+02 -1.6858783749946833e+03  6.0207717487186846e+02
+   28 -1.3837052677973970e+03  4.9654194958937097e+02 -9.0783879641816793e+02
+   29  1.1195044942562959e+03  1.1874984040640106e+03  3.0711352926618406e+02
+...

unittest/utils/test_lepton.cpp

@@ -85,6 +85,47 @@ TEST_F(LeptonUtilsTest, substitute)
     }
 }
 
+// zbl() custom function
+
+TEST(LeptonCustomFunction, zbl)
+{
+    Lepton::ZBLFunction zbl(1.0, 1.0, 1.0);
+    std::map<std::string, Lepton::CustomFunction *> functions = {std::make_pair("zbl", &zbl)};
+    std::map<std::string, double> variables = {std::make_pair("zi", 6), std::make_pair("zj", 6),
+                                               std::make_pair("r", 2.0)};
+
+    auto parsed = Lepton::Parser::parse("zbl(zi, zj, r)", functions);
+    auto zbldzi = parsed.differentiate("zi");
+    auto zbldzj = parsed.differentiate("zj");
+    auto zbldr  = parsed.differentiate("r");
+    auto zbld2r = zbldr.differentiate("r");
+
+    double value = parsed.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, 0.065721538245489763);
+    value = zbldr.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, -0.15481915325334394);
+    variables["r"] = 1.0;
+    value          = parsed.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, 1.0701488641432269);
+    value = zbldr.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, -3.6376386525054412);
+    variables["zi"] = 13.0;
+    value           = parsed.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, 1.8430432789454971);
+    value = zbldr.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, -6.5373118484557642);
+    variables["zj"] = 13.0;
+    value           = parsed.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, 3.1965196467438446);
+    value = zbldr.evaluate(variables);
+    ASSERT_DOUBLE_EQ(value, -11.804490148948526);
+
+    // check for unsupported derivatives
+    ASSERT_ANY_THROW(value = zbldzi.evaluate(variables));
+    ASSERT_ANY_THROW(value = zbldzj.evaluate(variables));
+    ASSERT_ANY_THROW(value = zbld2r.evaluate(variables));
+}
+
 /**
  * This is a custom function equal to f(x,y) = 2*x*y.
  */