energy/virial tallying for UB bonds

src/AMOEBA/angle_amoeba.cpp

@@ -574,10 +574,7 @@ void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag)
     f[i2][2] -= delz*fbond;
   }
 
-  // NOTE: there is no force on i2 for UB
+  if (evflag) ev_tally2(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
-
-  //if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,ebond,f1,f3,
-  //                     delx1,dely1,delz1,delx2,dely2,delz2);
 }
 
 /* ---------------------------------------------------------------------- */

src/angle.cpp

@@ -179,6 +179,7 @@ void Angle::ev_setup(int eflag, int vflag, int alloc)
 /* ----------------------------------------------------------------------
    tally energy and virial into global and per-atom accumulators
    virial = r1F1 + r2F2 + r3F3 = (r1-r2) F1 + (r3-r2) F3 = del1*f1 + del2*f3
+   called by standard 3-body angles
 ------------------------------------------------------------------------- */
 
 void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
@@ -341,6 +342,7 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
 /* ----------------------------------------------------------------------
    tally energy and virial into global and per-atom accumulators
    virial = r1F1 + r2F2 + r3F3 + r4F4
+   called by AngleAmoeba for its 4-body angle term
 ------------------------------------------------------------------------- */
 
 void Angle::ev_tally4(int i, int j, int k, int m, int nlocal, int newton_bond,
@@ -437,6 +439,92 @@ void Angle::ev_tally4(int i, int j, int k, int m, int nlocal, int newton_bond,
   }
 }
 
+
+/* ----------------------------------------------------------------------
+   tally energy and virial into global and per-atom accumulators
+   called by AngleAmoeba for its 2-body Urey-Bradley H-H bond term
+   identical to Bond:ev_tally()
+------------------------------------------------------------------------- */
+
+void Angle::ev_tally2(int i, int j, int nlocal, int newton_bond,
+                      double ebond, double fbond,
+                      double delx, double dely, double delz)
+{
+  double ebondhalf,v[6];
+
+  if (eflag_either) {
+    if (eflag_global) {
+      if (newton_bond) energy += ebond;
+      else {
+        ebondhalf = 0.5*ebond;
+        if (i < nlocal) energy += ebondhalf;
+        if (j < nlocal) energy += ebondhalf;
+      }
+    }
+    if (eflag_atom) {
+      ebondhalf = 0.5*ebond;
+      if (newton_bond || i < nlocal) eatom[i] += ebondhalf;
+      if (newton_bond || j < nlocal) eatom[j] += ebondhalf;
+    }
+  }
+
+  if (vflag_either) {
+    v[0] = delx*delx*fbond;
+    v[1] = dely*dely*fbond;
+    v[2] = delz*delz*fbond;
+    v[3] = delx*dely*fbond;
+    v[4] = delx*delz*fbond;
+    v[5] = dely*delz*fbond;
+
+    if (vflag_global) {
+      if (newton_bond) {
+        virial[0] += v[0];
+        virial[1] += v[1];
+        virial[2] += v[2];
+        virial[3] += v[3];
+        virial[4] += v[4];
+        virial[5] += v[5];
+      } else {
+        if (i < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+        if (j < nlocal) {
+          virial[0] += 0.5*v[0];
+          virial[1] += 0.5*v[1];
+          virial[2] += 0.5*v[2];
+          virial[3] += 0.5*v[3];
+          virial[4] += 0.5*v[4];
+          virial[5] += 0.5*v[5];
+        }
+      }
+    }
+
+    if (vflag_atom) {
+      if (newton_bond || i < nlocal) {
+        vatom[i][0] += 0.5*v[0];
+        vatom[i][1] += 0.5*v[1];
+        vatom[i][2] += 0.5*v[2];
+        vatom[i][3] += 0.5*v[3];
+        vatom[i][4] += 0.5*v[4];
+        vatom[i][5] += 0.5*v[5];
+      }
+      if (newton_bond || j < nlocal) {
+        vatom[j][0] += 0.5*v[0];
+        vatom[j][1] += 0.5*v[1];
+        vatom[j][2] += 0.5*v[2];
+        vatom[j][3] += 0.5*v[3];
+        vatom[j][4] += 0.5*v[4];
+        vatom[j][5] += 0.5*v[5];
+      }
+    }
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 double Angle::memory_usage()

src/angle.h

@@ -78,6 +78,7 @@ class Angle : protected Pointers {
   void ev_tally(int, int, int, int, int, double, double *, double *, double, double, double, double,
                 double, double);
   void ev_tally4(int, int, int, int, int, int, double, double *, double *, double *, double *);
+  void ev_tally2(int, int, int, int, double, double, double, double, double);
 };
 
 }    // namespace LAMMPS_NS