Added quartic function for explicit pair interaction in pair_style list.

doc/src/pair_list.rst

@@ -69,6 +69,7 @@ Here is an example file:
    15 259 lj126     1.0 1.0      50.0
    15 603 morse    10.0 1.2 2.0  10.0 # and another comment
    18 470 harmonic 50.0 1.2       5.0
+   19 332 quartic   5.0 -1.2 1.2  5.0
 
 The style *lj126* computes pairwise interactions with the formula
 
@@ -106,6 +107,21 @@ and the coefficients:
 
 Note that the usual 1/2 factor is included in :math:`K`.
 
+The style *quartic* computes pairwise interactions with the formula
+
+.. math::
+
+   E = K (r - r_c)^2 (r - r_c -b_1) (r - r_c - b_2) \qquad r < r_c
+
+and the coefficients:
+
+* :math:`K` (energy units)
+* :math:`b_1` (distance units)
+* :math:`b_2` (distance units)
+* :math:`r_c` (distance units)
+
+Note that the cutoff :math:`r_c` should always be specified to ensure zero energy and smooth force at cutoff.
+
 ----------
 
 Mixing, shift, table, tail correction, restart, rRESPA info

src/MISC/pair_list.cpp

@@ -32,13 +32,14 @@
 
 using namespace LAMMPS_NS;
 
-enum { NONE = 0, HARM, MORSE, LJ126 };
+enum { NONE = 0, HARM, MORSE, LJ126, QUARTIC };
 
 static std::map<std::string, int> stylename = {
   { "none",     NONE  },
   { "harmonic", HARM  },
   { "morse",    MORSE },
-  { "lj126",    LJ126 }
+  { "lj126",    LJ126 },
+  { "quartic",  QUARTIC }
 };
 
 // fast power function for integer exponent > 0
@@ -157,8 +158,21 @@ void PairList::compute(int eflag, int vflag)
         if (eflag_either)
           epair = 4.0*par.param.lj126.epsilon*r6inv
             * (sig6*sig6*r6inv - sig6) - par.offset;
+       
+      } else if (par.style == QUARTIC) {
+
+        const double r = sqrt(rsq);
+        double dr = r - sqrt(par.cutsq);
+	double ra = dr - par.param.quartic.b1;
+	double rb = dr - par.param.quartic.b2;
+	double r2 = dr * dr;
+        fpair = -par.param.quartic.k / r * (r2 * (ra + rb) + 2.0 * dr * ra * rb);
+
+        if (eflag_either)
+          epair = par.param.quartic.k * r2 * ra * rb;
       }
 
+
       if (newton_pair || i < nlocal) {
         f[i].x += dx*fpair;
         f[i].y += dy*fpair;
@@ -221,10 +235,10 @@ void PairList::settings(int narg, char **arg)
 
   // read and parse potential file only on MPI rank 0.
   if (comm->me == 0) {
-    int nharm, nmorse, nlj126, nskipped;
+    int nharm, nmorse, nlj126, nquartic, nskipped;
     FILE *fp = utils::open_potential(arg[0],lmp,nullptr);
     TextFileReader reader(fp,"pair list coeffs");
-    npairs = nharm = nmorse = nlj126 = nskipped = 0;
+    npairs = nharm = nmorse = nlj126 = nquartic = nskipped = 0;
 
     try {
       char *line;
@@ -258,6 +272,14 @@ void PairList::settings(int narg, char **arg)
           ++nlj126;
           break;
 
+        case QUARTIC:
+          oneparam.param.quartic.k = values.next_double();
+          oneparam.param.quartic.b1 = values.next_double();
+          oneparam.param.quartic.b2 = values.next_double();
+          oneparam.param.quartic.rc = values.next_double();
+          ++nquartic;
+          break;
+
         case NONE: // fallthrough
           error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}",
                          utils::trim(line));
@@ -265,7 +287,7 @@ void PairList::settings(int narg, char **arg)
           break;
         }
         if (values.has_next())
-          oneparam.cutsq = values.next_double();
+          oneparam.cutsq = mypow(values.next_double(), 2);
         else
           oneparam.cutsq = cut_global*cut_global;
 
@@ -274,8 +296,8 @@ void PairList::settings(int narg, char **arg)
     } catch (std::exception &e) {
       error->one(FLERR,"Error reading pair list coeffs file: {}", e.what());
     }
-    utils::logmesg(lmp, "Read {} ({}/{}/{}) interacting pair lines from {}. "
-                   "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, arg[0], nskipped);
+    utils::logmesg(lmp, "Read {} ({}/{}/{}/{}) interacting pair lines from {}. "
+                   "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, nquartic, arg[0], nskipped);
 
     memory->create(params,npairs,"pair_list:params");
     memcpy(params, myparams.data(),npairs*sizeof(list_param));
@@ -339,6 +361,12 @@ void PairList::init_style()
         const double r6inv = par.cutsq*par.cutsq*par.cutsq;
         const double sig6  = mypow(par.param.lj126.sigma,6);
         par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6);
+      
+      } else if (par.style == QUARTIC) {
+        // the offset is always 0 at rc
+        par.offset = 0.0;
+	// correct cutsq
+	par.cutsq = mypow(par.param.quartic.rc, 2);
       }
     }
   }

src/MISC/pair_list.h

@@ -49,11 +49,16 @@ class PairList : public Pair {
   struct lj126_p {
     double epsilon, sigma;
   };
+  struct quartic_p {
+    double k, b1, b2, rc;
+  };
+
 
   union param_u {
     harm_p harm;
     morse_p morse;
     lj126_p lj126;
+    quartic_p quartic;
   };
 
   struct list_param {