handle pair_modify shift and enforce the bond lepton has zero energy at r0

src/LEPTON/bond_lepton.cpp

@@ -32,7 +32,8 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-BondLepton::BondLepton(LAMMPS *_lmp) : Bond(_lmp), r0(nullptr), type2expression(nullptr)
+BondLepton::BondLepton(LAMMPS *_lmp) :
+    Bond(_lmp), r0(nullptr), type2expression(nullptr), offset(nullptr)
 {
   writedata = 1;
   reinitflag = 0;
@@ -46,6 +47,7 @@ BondLepton::~BondLepton()
     memory->destroy(setflag);
     memory->destroy(r0);
     memory->destroy(type2expression);
+    memory->destroy(offset);
   }
 }
 
@@ -133,7 +135,7 @@ template <int EVFLAG, int EFLAG, int NEWTON_BOND> void BondLepton::eval()
     if (EFLAG) {
       double &r_pot = bondpot[idx].getVariableReference("r");
       r_pot = dr;
-      ebond = bondpot[idx].evaluate();
+      ebond = bondpot[idx].evaluate() - offset[type];
     }
     if (EVFLAG) ev_tally(i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz);
   }
@@ -148,6 +150,7 @@ void BondLepton::allocate()
 
   memory->create(r0, np1, "bond:r0");
   memory->create(type2expression, np1, "bond:type2expression");
+  memory->create(offset, np1, "bond:offset");
   memory->create(setflag, np1, "bond:setflag");
   for (int i = 1; i < np1; i++) setflag[i] = 0;
 }
@@ -169,14 +172,15 @@ void BondLepton::coeff(int narg, char **arg)
   // remove whitespace and quotes from expression string and then
   // check if the expression can be parsed and evaluated without error
   std::string exp_one = LMP_Lepton::condense(arg[2]);
+  double offset_one = 0.0;
   try {
     auto parsed = LMP_Lepton::Parser::parse(exp_one);
     auto bondpot = parsed.createCompiledExpression();
     auto bondforce = parsed.differentiate("r").createCompiledExpression();
     double &r_pot = bondpot.getVariableReference("r");
     double &r_for = bondforce.getVariableReference("r");
-    r_for = r_pot = 1.0;
-    bondpot.evaluate();
+    r_for = r_pot = r0_one;
+    offset_one = bondpot.evaluate();
     bondforce.evaluate();
   } catch (std::exception &e) {
     error->all(FLERR, e.what());
@@ -195,6 +199,7 @@ void BondLepton::coeff(int narg, char **arg)
   for (int i = ilo; i <= ihi; i++) {
     r0[i] = r0_one;
     type2expression[i] = idx;
+    offset[i] = offset_one;
     setflag[i] = 1;
     count++;
   }
@@ -302,7 +307,7 @@ double BondLepton::single(int type, double rsq, int /*i*/, int /*j*/, double &ff
   double ebond = 0.0;
   if (r > 0.0) {
     fforce = -bondforce.evaluate() / r;
-    ebond = bondpot.evaluate();
+    ebond = bondpot.evaluate() - offset[type];
   }
   return ebond;
 }

src/LEPTON/bond_lepton.h

@@ -41,6 +41,7 @@ class BondLepton : public Bond {
   std::vector<std::string> expressions;
   double *r0;
   int *type2expression;
+  double *offset;
 
   virtual void allocate();
 

src/LEPTON/pair_lepton.cpp

@@ -33,7 +33,8 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairLepton::PairLepton(LAMMPS *lmp) : Pair(lmp), cut(nullptr), type2expression(nullptr)
+PairLepton::PairLepton(LAMMPS *lmp) :
+    Pair(lmp), cut(nullptr), type2expression(nullptr), offset(nullptr)
 {
   respa_enable = 0;
   single_enable = 1;
@@ -53,6 +54,7 @@ PairLepton::~PairLepton()
     memory->destroy(cutsq);
     memory->destroy(setflag);
     memory->destroy(type2expression);
+    memory->destroy(offset);
   }
 }
 
@@ -148,7 +150,8 @@ template <int EVFLAG, int EFLAG, int NEWTON_PAIR> void PairLepton::eval()
         if (EFLAG) {
           double &r_pot = pairpot[idx].getVariableReference("r");
           r_pot = r;
-          evdwl = factor_lj * pairpot[idx].evaluate();
+          evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
+          evdwl *= factor_lj;
         }
 
         if (EVFLAG) ev_tally(i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz);
@@ -176,6 +179,7 @@ void PairLepton::allocate()
   memory->create(cut, np1, np1, "pair:cut");
   memory->create(cutsq, np1, np1, "pair:cutsq");
   memory->create(type2expression, np1, np1, "pair:type2expression");
+  memory->create(offset, np1, np1, "pair:offset");
 }
 
 /* ----------------------------------------------------------------------
@@ -249,8 +253,18 @@ double PairLepton::init_one(int i, int j)
 {
   if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
 
+  if (offset_flag) {
+    auto parsed = LMP_Lepton::Parser::parse(expressions[type2expression[i][j]]);
+    auto pairpot = parsed.createCompiledExpression();
+    double &r_pot = pairpot.getVariableReference("r");
+    r_pot = 1.0;
+    offset[i][j] = pairpot.evaluate();
+  } else
+    offset[i][j] = 0.0;
+
   cut[j][i] = cut[i][j];
   type2expression[j][i] = type2expression[i][j];
+  offset[j][i] = offset[i][j];
 
   return cut[i][j];
 }
@@ -341,6 +355,7 @@ void PairLepton::read_restart(FILE *fp)
 void PairLepton::write_restart_settings(FILE *fp)
 {
   fwrite(&cut_global, sizeof(double), 1, fp);
+  fwrite(&offset_flag, sizeof(int), 1, fp);
 }
 
 /* ----------------------------------------------------------------------
@@ -349,8 +364,12 @@ void PairLepton::write_restart_settings(FILE *fp)
 
 void PairLepton::read_restart_settings(FILE *fp)
 {
-  if (comm->me == 0) { utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error); }
+  if (comm->me == 0) {
+    utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
+    utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
+  }
   MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
+  MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
 }
 
 /* ----------------------------------------------------------------------
@@ -389,5 +408,5 @@ double PairLepton::single(int /* i */, int /* j */, int itype, int jtype, double
 
   r_pot = r_for = r;
   fforce = -pairforce.evaluate() / r * factor_lj;
-  return pairpot.evaluate() * factor_lj;
+  return (pairpot.evaluate() - offset[itype][jtype]) * factor_lj;
 }

src/LEPTON/pair_lepton.h

@@ -54,6 +54,7 @@ class PairLepton : public Pair {
   std::vector<std::string> expressions;
   double **cut;
   int **type2expression;
+  double **offset;
   double cut_global;
 
   virtual void allocate();

src/OPENMP/bond_lepton_omp.cpp

@@ -140,7 +140,7 @@ void BondLeptonOMP::eval(int nfrom, int nto, ThrData *const thr)
     if (EFLAG) {
       double &r_pot = bondpot[idx].getVariableReference("r");
       r_pot = dr;
-      ebond = bondpot[idx].evaluate();
+      ebond = bondpot[idx].evaluate() - offset[type];
     }
     if (EVFLAG)
       ev_tally_thr(this, i1, i2, nlocal, NEWTON_BOND, ebond, fbond, delx, dely, delz, thr);

src/OPENMP/pair_lepton_omp.cpp

@@ -144,7 +144,8 @@ void PairLeptonOMP::eval(int iifrom, int iito, ThrData *const thr)
         if (EFLAG) {
           double &r_pot = pairpot[idx].getVariableReference("r");
           r_pot = r;
-          evdwl = factor_lj * pairpot[idx].evaluate();
+          evdwl = pairpot[idx].evaluate() - offset[itype][jtype];
+          evdwl *= factor_lj;
         }
 
         if (EVFLAG)