1) Added born_matrix method to bond_harmonic, angle_cosine_squared and dihedral_nharmonic

2) Changed compute_born_matrix so it can takes arguments pair bond angle dihedral to output only selected elements. Also moved flags to creator method and added error/warning messages

src/EXTRA-COMPUTE/compute_born_matrix.cpp

@@ -43,13 +43,26 @@ using namespace LAMMPS_NS;
 
 #define BIG 1000000000
 
-static int const albe[21][2] = {
+// This table is used to pick the 3d rij vector indices used to
+// compute the 6 indices long Voigt stress vector
+static int const sigma_albe[6][2] = {
+    {0, 0},    // s11
+    {1, 1},    // s22
+    {2, 2},    // s33
+    {1, 2},    // s44
+    {0, 2},    // s55
+    {0, 1},    // s66
+};
+
+// This table is used to pick the correct indices from the Voigt
+// stress vector to compute the Cij matrix (21 terms, see doc) contribution
+static int const C_albe[21][2] = {
     {0, 0},    // C11
     {1, 1},    // C22
     {2, 2},    // C33
-    {1, 2},    // C44
+    {3, 3},    // C44
-    {0, 2},    // C55
+    {4, 4},    // C55
-    {0, 1},    // C66
+    {5, 5},    // C66
     {0, 1},    // C12
     {0, 2},    // C13
     {0, 3},    // C14
@@ -67,6 +80,8 @@ static int const albe[21][2] = {
     {4, 5}     // C56
 };
 
+// This table is used to pick the 3d rij vector indices used to
+// compute the 21 indices long Cij matrix
 static int const albemunu[21][4] = {
     {0, 0, 0, 0},    // C11
     {1, 1, 1, 1},    // C22
@@ -103,6 +118,75 @@ ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : Comput
 
   // Error check
 
+  if (narg < 3) error->all(FLERR,"Illegal compute pressure command");
+
+  pairflag = 0;
+  bondflag = 0;
+  angleflag = 0;
+  dihedflag = 0;
+  impflag = 0;
+  kspaceflag = 0;
+  if (narg == 3) {
+    pairflag = 1;
+    bondflag = 1;
+    angleflag = 1;
+    dihedflag = 1;
+    impflag = 1;
+    kspaceflag = 1;
+  } else {
+    int iarg = 3;
+    while (iarg < narg) {
+      if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
+      else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
+      else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
+      else if (strcmp(arg[iarg],"dihedral") == 0) dihedflag = 1;
+      else if (strcmp(arg[iarg],"improper") == 0) impflag = 1;
+      else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
+      else error->all(FLERR,"Illegal compute pressure command");
+      ++iarg;
+    }
+  }
+
+  // Error check
+
+  if (pairflag) {
+    if (force->pair) {
+      if (force->pair->born_matrix_enable == 0) {
+        if (comm->me == 0) error->warning(FLERR, "Pair style does not support compute born/matrix");
+      }
+    }
+  }
+  if (bondflag) {
+    if (force->bond) {
+      if (force->bond->born_matrix_enable == 0) {
+        if (comm->me == 0) error->warning(FLERR, "Bond style does not support compute born/matrix");
+      }
+    }
+  }
+  if (angleflag) {
+    if (force->angle) {
+      if (force->angle->born_matrix_enable == 0) {
+        if (comm->me == 0) error->warning(FLERR, "Angle style does not support compute born/matrix");
+      }
+    }
+  }
+  if (dihedflag) {
+    if (force->dihedral) {
+     if (force->dihedral->born_matrix_enable == 0) {
+       if (comm->me == 0) error->warning(FLERR, "Dihedral style does not support compute born/matrix");
+     }
+    }
+  }
+  if (impflag) {
+    if (force->improper) {
+      if (force->improper->born_matrix_enable == 0) {
+        if (comm->me == 0) error->warning(FLERR, "Improper style does not support compute born/matrix");
+      }
+    }
+  }
+  if (kspaceflag) {
+    error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix");
+  }
   // Initialize some variables
 
   values_local = values_global = vector = nullptr;
@@ -138,42 +222,16 @@ void ComputeBornMatrix::init()
   // Timestep Value
   dt = update->dt;
 
-  pairflag = 0;
+  // IF everything works fine,
-  bondflag = 0;
+  // this is to be removed
-  angleflag = 0;
+  //
-  dihedflag = 0;
+  // if (force->pair) pairflag = 1;
-  impflag = 0;
+  // if (force->bond) bondflag = 1;
-  kspaceflag = 0;
+  // if (force->angle) angleflag = 1;
+  // if (force->dihedral) dihedflag = 1;
+  // if (force->improper) impflag = 1;
+  // if (force->kspace) kspaceflag = 1;
 
-  if (force->pair) pairflag = 1;
-  if (force->bond) bondflag = 1;
-  if (force->angle) angleflag = 1;
-  if (force->dihedral) dihedflag = 1;
-  if (force->improper) impflag = 1;
-  if (force->kspace) kspaceflag = 1;
-
-  // Error check
-
-  if (comm->me == 0) {
-    if (pairflag && (force->pair->born_matrix_enable == 0)) {
-      error->all(FLERR, "Pair style does not support compute born/matrix");
-    }
-    if (bondflag && (force->bond->born_matrix_enable == 0)) {
-      error->warning(FLERR, "Bond style does not support compute born/matrix");
-    }
-    if (angleflag && (force->angle->born_matrix_enable == 0)) {
-      error->warning(FLERR, "Angle style does not support compute born/matrix");
-    }
-    if (dihedflag && (force->dihedral->born_matrix_enable == 0)) {
-      error->warning(FLERR, "Dihedral style does not support compute born/matrix");
-    }
-    if (impflag && (force->improper->born_matrix_enable == 0)) {
-      error->warning(FLERR, "Improper style does not support compute born/matrix");
-    }
-    if (kspaceflag) {
-      error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix");
-    }
-  }
 
   // need an occasional half neighbor list
   int irequest = neighbor->request((void *) this);
@@ -204,7 +262,6 @@ void ComputeBornMatrix::compute_vector()
 
   // Compute Born contribution
   if (pairflag) compute_pairs();
-  // For now these functions are commented
   if (bondflag) compute_bonds();
   if (angleflag) compute_angles();
   if (dihedflag) compute_dihedrals();
@@ -331,11 +388,9 @@ void ComputeBornMatrix::compute_pairs()
    all atoms in interaction must be known to proc
    if bond is deleted or turned off (type <= 0)
    do not count or count contribution
-   COMMENTED FOR NOW
 ---------------------------------------------------------------------- */
 void ComputeBornMatrix::compute_bonds()
 {
-  /* ----------------------------------------------------------------------
   int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar;
   tagint tagprev;
   double dx,dy,dz,rsq;
@@ -426,7 +481,6 @@ void ComputeBornMatrix::compute_bonds()
     }
     m += 21;
   }
-------------------------------------------------------------------------- */
 }
 
 /* ----------------------------------------------------------------------
@@ -436,11 +490,9 @@ void ComputeBornMatrix::compute_bonds()
    all atoms in interaction must be known to proc
    if bond is deleted or turned off (type <= 0)
    do not count or count contribution
-   COMMENTED FOR NOW
 ---------------------------------------------------------------------- */
 void ComputeBornMatrix::compute_angles()
 {
-  /* ----------------------------------------------------------------------
   int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar;
   tagint tagprev;
   double delx1,dely1,delz1,delx2,dely2,delz2;
@@ -568,8 +620,8 @@ void ComputeBornMatrix::compute_angles()
         c = 0;
         d = 0;
         for (i = 0; i<6; i++) {
-          a = albe[i][0];
+          a = sigma_albe[i][0];
-          b = albe[i][1];
+          b = sigma_albe[i][1];
           dcos[i] = cost*(del1[a]*del2[b]+del1[b]*del2[a]*r1r2inv -
                           del1[a]*del1[b]*rsq1inv - del2[a]*del2[b]*rsq2inv);
         }
@@ -578,8 +630,8 @@ void ComputeBornMatrix::compute_angles()
           b = albemunu[i][1];
           c = albemunu[i][2];
           d = albemunu[i][3];
-          e = albe[i][0];
+          e = C_albe[i][0];
-          f = albe[i][1];
+          f = C_albe[i][1];
           d2lncos[i] = 2*(del1[a]*del1[b]*del1[c]*del1[d]*rsq1inv*rsq1inv +
                           del2[a]*del2[b]*del2[c]*del2[d]*rsq2inv*rsq2inv) -
                          (del1[a]*del2[b]+del1[b]*del2[a]) *
@@ -592,7 +644,6 @@ void ComputeBornMatrix::compute_angles()
     }
   m+=21;
   }
-------------------------------------------------------------------------- */
 }
 
 /* ----------------------------------------------------------------------
@@ -601,12 +652,10 @@ void ComputeBornMatrix::compute_angles()
    all atoms in interaction must be in group
    all atoms in interaction must be known to proc
    if flag is set, compute requested info about dihedral
-   COMMENTED FOR NOW
 ------------------------------------------------------------------------- */
 
 void ComputeBornMatrix::compute_dihedrals()
 {
-  /* ----------------------------------------------------------------------
   int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,dtype,ivar;
   tagint tagprev;
   double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
@@ -781,8 +830,8 @@ void ComputeBornMatrix::compute_dihedrals()
         e = 0;
         f = 0;
         for (i = 0; i<6; i++) {
-          a = albe[i][0];
+          a = sigma_albe[i][0];
-          b = albe[i][1];
+          b = sigma_albe[i][1];
           dmm[i] = 2*(b2sq*b1[a]*b1[b]+b1sq*b2[a]*b2[b] - b1b2*(b1[a]*b2[b]+b1[b]*b2[a]));
           dnn[i] = 2*(b3sq*b2[a]*b2[b]+b2sq*b3[a]*b3[b] - b2b3*(b2[a]*b3[b]+b2[b]*b3[a]));
           dmn[i] = b1b2*(b2[a]*b3[b]+b2[b]*b3[a]) + b2b3*(b1[a]*b2[b]+b1[b]*b2[a])
@@ -794,8 +843,8 @@ void ComputeBornMatrix::compute_dihedrals()
           b = albemunu[i][1];
           c = albemunu[i][2];
           d = albemunu[i][3];
-          e = albe[i][0];
+          e = C_albe[i][0];
-          f = albe[i][1];
+          f = C_albe[i][1];
           d2mm[i] = 4*(b1[a]*b1[b]*b2[c]*b2[d] + b1[c]*b1[d]*b2[a]*b2[b])
                   - 8*(b1[a]*b2[b]+b1[b]*b2[a])*(b1[c]*b2[d]+b1[d]*b2[c]);
           d2nn[i] = 4*(b2[a]*b2[b]*b3[c]*b3[d] + b2[c]*b2[d]*b3[a]*b3[b])
@@ -810,12 +859,11 @@ void ComputeBornMatrix::compute_dihedrals()
                   + ra2inv*ra2inv*dmm[e]*dmm[f]
                   - ra2inv*d2mm[i]
                   + rb2inv*rb2inv*dnn[e]*dnn[f]
-                  - rb2inv*d2nn[i] );
+                  - rb2inv*d2nn[i]);
           values_local[m+i] += dudih*d2cos[i] + du2dih*dcos[e]*dcos[f];
         }
       }
     }
   m+=21;
   }
-------------------------------------------------------------------------- */
 }

src/EXTRA-MOLECULE/dihedral_nharmonic.cpp

@@ -39,6 +39,7 @@ DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp)
 {
   writedata = 1;
   a = nullptr;
+  born_matrix_enable = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -336,3 +337,63 @@ void DihedralNHarmonic::write_data(FILE *fp)
   }
 
 }
+
+/* ----------------------------------------------------------------------*/
+
+void DihedralNHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4,
+                             double &dudih, double &du2dih) {
+  int i,type;
+  double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
+  double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv;
+  double c,s,kf;
+
+  int **dihedrallist = neighbor->dihedrallist;
+  double **x = atom->x;
+
+  type = dihedrallist[nd][4];
+
+  vb1x = x[i1][0] - x[i2][0];
+  vb1y = x[i1][1] - x[i2][1];
+  vb1z = x[i1][2] - x[i2][2];
+
+  vb2x = x[i3][0] - x[i2][0];
+  vb2y = x[i3][1] - x[i2][1];
+  vb2z = x[i3][2] - x[i2][2];
+
+  vb2xm = -vb2x;
+  vb2ym = -vb2y;
+  vb2zm = -vb2z;
+
+  vb3x = x[i4][0] - x[i3][0];
+  vb3y = x[i4][1] - x[i3][1];
+  vb3z = x[i4][2] - x[i3][2];
+
+  // c,s calculation
+
+  ax = vb1y*vb2zm - vb1z*vb2ym;
+  ay = vb1z*vb2xm - vb1x*vb2zm;
+  az = vb1x*vb2ym - vb1y*vb2xm;
+  bx = vb3y*vb2zm - vb3z*vb2ym;
+  by = vb3z*vb2xm - vb3x*vb2zm;
+  bz = vb3x*vb2ym - vb3y*vb2xm;
+
+  rasq = ax*ax + ay*ay + az*az;
+  rbsq = bx*bx + by*by + bz*bz;
+
+  ra2inv = rb2inv = 0.0;
+  if (rasq > 0) ra2inv = 1.0/rasq;
+  if (rbsq > 0) rb2inv = 1.0/rbsq;
+  rabinv = sqrt(ra2inv*rb2inv);
+
+  c = (ax*bx + ay*by + az*bz)*rabinv;
+
+  dudih = 0.0;
+  du2dih = 0.0;
+  for (i = 1; i < nterms[type]; i++) {
+    dudih += this->a[type][i]*i*pow(c,i-1);
+  }
+
+  for (i = 2; i < nterms[type]; i++) {
+    du2dih += this->a[type][i]*i*(i-1)*pow(c, i-2);
+  }
+}

src/EXTRA-MOLECULE/dihedral_nharmonic.h

@@ -33,6 +33,7 @@ class DihedralNHarmonic : public Dihedral {
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_data(FILE *);
+  void born_matrix(int /*dtype*/, int, int, int, int, double&, double&);
 
  protected:
   int *nterms;

src/MOLECULE/angle_cosine_squared.cpp

@@ -40,6 +40,7 @@ AngleCosineSquared::AngleCosineSquared(LAMMPS *lmp) : Angle(lmp)
 {
   k = nullptr;
   theta0 = nullptr;
+  born_matrix_enable = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -263,3 +264,32 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3)
   double tk = k[type] * dcostheta;
   return tk*dcostheta;
 }
+
+/* ---------------------------------------------------------------------- */
+
+void AngleCosineSquared::born_matrix(int type, int i1, int i2, int i3, double& du, double& du2)
+{
+  double **x = atom->x;
+
+  double delx1 = x[i1][0] - x[i2][0];
+  double dely1 = x[i1][1] - x[i2][1];
+  double delz1 = x[i1][2] - x[i2][2];
+  domain->minimum_image(delx1,dely1,delz1);
+  double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
+
+  double delx2 = x[i3][0] - x[i2][0];
+  double dely2 = x[i3][1] - x[i2][1];
+  double delz2 = x[i3][2] - x[i2][2];
+  domain->minimum_image(delx2,dely2,delz2);
+  double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
+
+  double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
+  c /= r1*r2;
+  if (c > 1.0) c = 1.0;
+  if (c < -1.0) c = -1.0;
+
+  double dcostheta = c - cos(theta0[type]);
+  double tk = k[type] * dcostheta;
+  du2 = 2*k[type];
+  du = du2*dcostheta;
+}

src/MOLECULE/angle_cosine_squared.h

@@ -35,6 +35,7 @@ class AngleCosineSquared : public Angle {
   void read_restart(FILE *);
   void write_data(FILE *);
   virtual double single(int, int, int, int);
+  void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2);
 
  protected:
   double *k, *theta0;

src/MOLECULE/bond_harmonic.cpp

@@ -31,6 +31,7 @@ using namespace LAMMPS_NS;
 BondHarmonic::BondHarmonic(LAMMPS *lmp) : Bond(lmp)
 {
   reinitflag = 1;
+  born_matrix_enable = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -201,6 +202,20 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/,
   return rk*dr;
 }
 
+
+/* ---------------------------------------------------------------------- */
+
+void BondHarmonic::born_matrix(int type, double rsq, int /*i*/, int /*j*/,
+                            double &du, double& du2)
+{
+  double r = sqrt(rsq);
+  double dr = r - r0[type];
+  du2 = 0.0;
+  du = 0.0;
+  du2 = 2*k[type];
+  if (r > 0.0) du = du2*dr;
+}
+
 /* ----------------------------------------------------------------------
     Return ptr to internal members upon request.
 ------------------------------------------------------------------------ */

src/MOLECULE/bond_harmonic.h

@@ -35,6 +35,7 @@ class BondHarmonic : public Bond {
   virtual void read_restart(FILE *);
   void write_data(FILE *);
   double single(int, double, int, int, double &);
+  void born_matrix(int, double, int, int, double &, double &);
   virtual void *extract(const char *, int &);
 
  protected: