derived class of sw

src/MANYBODY/pair_sw.cpp

@@ -28,19 +28,17 @@
 #include "neighbor.h"
 #include "potential_file_reader.h"
 #include "tokenizer.h"
-#include "math_const.h"
 
 #include <cmath>
 #include <cstring>
 
 using namespace LAMMPS_NS;
-using namespace MathConst;
 
 #define DELTA 4
 
 /* ---------------------------------------------------------------------- */
 
-PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
+PairSW::PairSW(LAMMPS *lmp) : Pair(lmp), variant(SW)
 {
   single_enable = 0;
   restartinfo = 0;
@@ -233,24 +231,9 @@ void PairSW::allocate()
    global settings
 ------------------------------------------------------------------------- */
 
-void PairSW::settings(int narg, char **arg)
+void PairSW::settings(int narg, char **/*arg*/)
 {
-  //if (narg != 0) error->all(FLERR,"Illegal pair_style command");
-  // default values
-
-  modify_flag = 0;
-
-  // process optional keywords
-
-  int iarg = 0;
-
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"modify") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style command");
-      modify_flag = utils::numeric(FLERR,arg[iarg+1],false,lmp);
-      iarg += 2;
-    } else error->all(FLERR,"Illegal pair_style command");
-  }
+  if (narg != 0) error->all(FLERR,"Illegal pair_style command");
 }
 
 /* ----------------------------------------------------------------------
@@ -309,8 +292,22 @@ void PairSW::read_file(char *file)
   // open file on proc 0
 
   if (comm->me == 0) {
-    PotentialFileReader reader(lmp, file, "sw", unit_convert_flag);
-    char *line;
+    std::string potential_name;
+    switch (variant) {
+    case SW:
+      potential_name = "sw";
+      break;
+
+    case SW_MOD:
+      potential_name = "sw/mod";
+      break;
+
+    default:
+      error->one(FLERR,"Unknown SW style variant {}",variant);
+    }
+
+    PotentialFileReader reader(lmp, file, potential_name, unit_convert_flag);
+    char * line;
 
     // transparently convert units for supported conversions
 
@@ -345,7 +342,8 @@ void PairSW::read_file(char *file)
 
         if (nparams == maxparam) {
           maxparam += DELTA;
-          params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),"pair:params");
+          params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+                                              "pair:params");
 
           // make certain all addional allocated storage is initialized
           // to avoid false positives when checking with valgrind
@@ -502,7 +500,7 @@ void PairSW::threebody(Param *paramij, Param *paramik, Param *paramijk,
   double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
   double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
   double rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2;
-  double facang,facang12,csfacang,csfac1,csfac2,factor;
+  double facang,facang12,csfacang,csfac1,csfac2;
 
   r1 = sqrt(rsq1);
   rinvsq1 = 1.0/rsq1;
@@ -521,14 +519,6 @@ void PairSW::threebody(Param *paramij, Param *paramik, Param *paramijk,
   rinv12 = 1.0/(r1*r2);
   cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
   delcs = cs - paramijk->costheta;
-  // Modification
-  if(modify_flag) {
-    if(fabs(delcs) >= 0.35) delcs = 0.0;
-    else if(fabs(delcs) < 0.35 && fabs(delcs) > 0.25) {
-      factor = 0.5 + 0.5*cos(MY_PI*(delcs-0.25)/(0.35-0.25));
-      delcs *= factor;
-    }
-  }
   delcssq = delcs*delcs;
 
   facexp = expgsrainv1*expgsrainv2;

src/MANYBODY/pair_sw.h

@@ -36,6 +36,8 @@ class PairSW : public Pair {
 
   static constexpr int NPARAMS_PER_LINE = 14;
 
+  enum { SW, SW_MOD };    // for telling class variants apart in shared code
+
   struct Param {
     double epsilon, sigma;
     double littlea, lambda, gamma, costheta;
@@ -49,13 +51,12 @@ class PairSW : public Pair {
   };
 
  protected:
+  int variant;
   double cutmax;      // max cutoff for all elements
   Param *params;      // parameter set for an I-J-K interaction
   int maxshort;       // size of short neighbor list array
   int *neighshort;    // short neighbor list array
 
-  int modify_flag;    // flag to turn on/off the modification
-
   virtual void allocate();
   void read_file(char *);
   virtual void setup_params();

src/MANYBODY/pair_sw_mod.cpp

@@ -0,0 +1,96 @@
+// clang-format off
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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 authors: Jin-Wu Jiang (SHU) and Wengen Ouyang (WHU)
+------------------------------------------------------------------------- */
+
+#include "pair_sw_mod.h"
+
+#include "math_const.h"
+
+#include <cmath>
+#include <cstring>
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+/* ---------------------------------------------------------------------- */
+
+PairSWMOD::PairSWMOD(LAMMPS *lmp) : PairSW(lmp) {
+  variant = SW_MOD;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSWMOD::threebody(Param *paramij, Param *paramik, Param *paramijk,
+                       double rsq1, double rsq2,
+                       double *delr1, double *delr2,
+                       double *fj, double *fk, int eflag, double &eng)
+{
+  double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
+  double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
+  double rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2;
+  double facang,facang12,csfacang,csfac1,csfac2,factor;
+
+  r1 = sqrt(rsq1);
+  rinvsq1 = 1.0/rsq1;
+  rainv1 = 1.0/(r1 - paramij->cut);
+  gsrainv1 = paramij->sigma_gamma * rainv1;
+  gsrainvsq1 = gsrainv1*rainv1/r1;
+  expgsrainv1 = exp(gsrainv1);
+
+  r2 = sqrt(rsq2);
+  rinvsq2 = 1.0/rsq2;
+  rainv2 = 1.0/(r2 - paramik->cut);
+  gsrainv2 = paramik->sigma_gamma * rainv2;
+  gsrainvsq2 = gsrainv2*rainv2/r2;
+  expgsrainv2 = exp(gsrainv2);
+
+  rinv12 = 1.0/(r1*r2);
+  cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
+  delcs = cs - paramijk->costheta;
+  // Modification to delcs
+  if(fabs(delcs) >= 0.35) delcs = 0.0;
+  else if(fabs(delcs) < 0.35 && fabs(delcs) > 0.25) {
+    factor = 0.5 + 0.5*cos(MY_PI*(delcs-0.25)/(0.35-0.25));
+    delcs *= factor;
+  }
+  delcssq = delcs*delcs;
+
+  facexp = expgsrainv1*expgsrainv2;
+
+  // facrad = sqrt(paramij->lambda_epsilon*paramik->lambda_epsilon) *
+  //          facexp*delcssq;
+
+  facrad = paramijk->lambda_epsilon * facexp*delcssq;
+  frad1 = facrad*gsrainvsq1;
+  frad2 = facrad*gsrainvsq2;
+  facang = paramijk->lambda_epsilon2 * facexp*delcs;
+  facang12 = rinv12*facang;
+  csfacang = cs*facang;
+  csfac1 = rinvsq1*csfacang;
+
+  fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12;
+  fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12;
+  fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12;
+
+  csfac2 = rinvsq2*csfacang;
+
+  fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12;
+  fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12;
+  fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12;
+
+  if (eflag) eng = facrad;
+}

src/MANYBODY/pair_sw_mod.h

@@ -0,0 +1,90 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   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(sw/mod,PairSWMOD);
+// clang-format on
+#else
+
+#ifndef LMP_PAIR_SW_MOD_H
+#define LMP_PAIR_SW_MOD_H
+
+#include "pair_sw.h"
+
+namespace LAMMPS_NS {
+
+class PairSWMOD : public PairSW {
+ public:
+  PairSWMOD(class LAMMPS *);
+  ~PairSWMOD() {}
+
+ protected:
+  void threebody(Param *, Param *, Param *, double, double, double *, double *, double *, double *,
+                 int, double &);
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair style Stillinger-Weber requires atom IDs
+
+This is a requirement to use the SW potential.
+
+E: Pair style Stillinger-Weber requires newton pair on
+
+See the newton command.  This is a restriction to use the SW
+potential.
+
+E: All pair coeffs are not set
+
+All pair coefficients must be set in the data file or by the
+pair_coeff command before running a simulation.
+
+E: Cannot open Stillinger-Weber potential file %s
+
+The specified SW potential file cannot be opened.  Check that the path
+and name are correct.
+
+E: Incorrect format in Stillinger-Weber potential file
+
+Incorrect number of words per line in the potential file.
+
+E: Illegal Stillinger-Weber parameter
+
+One or more of the coefficients defined in the potential file is
+invalid.
+
+E: Potential file has duplicate entry
+
+The potential file has more than one entry for the same element.
+
+E: Potential file is missing an entry
+
+The potential file does not have a needed entry.
+
+*/