Add single2 method to Pair that can compute and return the Hessian matrix

src/pair.cpp

@@ -301,7 +301,7 @@ void Pair::init_style()
    specific pair style can override this function
 ------------------------------------------------------------------------- */
 
-void Pair::init_list(int /*which*/, NeighList *ptr)
+void Pair::init_list(int which, NeighList *ptr)
 {
   list = ptr;
 }
@@ -1740,6 +1740,18 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
 
 /* ---------------------------------------------------------------------- */
 
+void Pair::pairTensor(double fforce, double dfac, double delr[3], double phiTensor[6]) {
+  int m = 0;
+  for (int k=0; k<3; k++) {
+    phiTensor[m] = fforce;
+    for (int l=k; l<3; l++) {
+      if (l>k) phiTensor[m] = 0;
+      phiTensor[m++] += delr[k]*delr[l] * dfac;
+    }
+  }
+}
+/* ---------------------------------------------------------------------- */
+
 double Pair::memory_usage()
 {
   double bytes = comm->nthreads*maxeatom * sizeof(double);

src/pair.h

@@ -145,6 +145,16 @@ class Pair : protected Pointers {
     return 0.0;
   }
 
+  void pairTensor(double fforce, double dfac, double delr[3], double phiTensor[6]);
+
+  virtual double single2(int, int, int, int,
+                        double, double[3], double, double,
+                        double& fforce, double d2u[6]) {
+    fforce = 0.0;
+    for (int i=0; i<6; i++) d2u[i] = 0;
+    return 0.0;
+  }
+
   virtual void settings(int, char **) = 0;
   virtual void coeff(int, char **) = 0;
 
@@ -207,10 +217,6 @@ class Pair : protected Pointers {
 
   typedef union {int i; float f;} union_int_float_t;
 
-  // Accessor for the user-intel package to determine virial calc for hybrid
-
-  inline int fdotr_is_set() const { return vflag_fdotr; }
-
  protected:
   int vflag_fdotr;
   int maxeatom,maxvatom;

src/pair_lj_smooth_linear.cpp

@@ -326,9 +326,9 @@ void PairLJSmoothLinear::read_restart_settings(FILE *fp)
 
 /* ---------------------------------------------------------------------- */
 
-double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype,
+double PairLJSmoothLinear::single(int i, int j, int itype, int jtype,
                                   double rsq,
-                                  double /*factor_coul*/, double factor_lj,
+                                  double factor_coul, double factor_lj,
                                   double &fforce)
 {
   double r2inv,r6inv,forcelj,philj,r,rinv;
@@ -347,3 +347,26 @@ double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype,
 
   return factor_lj*philj;
 }
+
+double PairLJSmoothLinear::single2(int i, int j, int itype, int jtype, double rsq,
+                         double delr[3], double factor_coul, double factor_lj,
+                         double &fforce, double d2u[6])
+{
+  double r2inv,r6inv,forcelj,philj,r,rinv;
+
+  r2inv = 1.0/rsq;
+  r6inv = r2inv*r2inv*r2inv;
+  rinv  = sqrt(r2inv);
+  r     = sqrt(rsq);
+  forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
+  forcelj = rinv*forcelj - dljcut[itype][jtype];
+  fforce = factor_lj*forcelj*rinv;
+
+  philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
+  philj = philj - ljcut[itype][jtype]
+                + (r-cut[itype][jtype])*dljcut[itype][jtype];
+
+  double d2r = factor_lj * r6inv * (13.0*lj1[itype][jtype]*r6inv - 7.0*lj2[itype][jtype])/rsq;
+  pairTensor(fforce, -(fforce + d2r) / rsq, delr, d2u);
+  return factor_lj*philj;
+}

src/pair_lj_smooth_linear.h

@@ -38,6 +38,7 @@ class PairLJSmoothLinear : public Pair {
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
   double single(int, int, int, int, double, double, double, double &);
+  double single2(int, int, int, int, double, double[3], double, double, double&, double[6]);
 
  protected:
   double cut_global;