Improve energy conservation of vdwForceSwitching/LJcorrection 69/4869/4
authorradakb <brian.radak@gmail.com>
Wed, 19 Dec 2018 16:20:23 +0000 (11:20 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Wed, 19 Dec 2018 21:32:46 +0000 (15:32 -0600)
The existing implementation for vdwForceSwitching + LJcorrection
introduces a discontinuity into the potential (but not the force)
in order to make the long-range correction more well-defined.
While this discontinuity is very small for normal cutoffs, it
can cause detectable issues at moderate lengths (~8 A).

This revision reverts to the old vdwForceSwitching behavior and
modifies the long-range correction to include a "core" correction
term based on the same approximation -- this is theoretically
dubious, but works quite well in practice.

On a practical note, this change greatly simplifies the code
because vdwForceSwitching behavior is now uniform. Supporting
multiple branches in the alchemical code will now be substantially
easier.

Change-Id: Ia53bc80f79e1976f2f0ca62f2f0d6988a2fe960a

src/ComputeNonbondedBase.h
src/ComputeNonbondedBase2.h
src/ComputeNonbondedFEP.C
src/ComputeNonbondedTI.C
src/ComputeNonbondedUtil.C
src/ComputeNonbondedUtil.h
src/Molecule.C

index 2e04e8f..83244c6 100644 (file)
@@ -383,7 +383,6 @@ void ComputeNonbondedUtil :: NAME
     const BigReal switchfactor = 1./((cutoff2 - switchdist2)*(cutoff2 - switchdist2)*(cutoff2 - switchdist2));
     const BigReal alchVdwShiftCoeff = ComputeNonbondedUtil::alchVdwShiftCoeff;
     const Bool vdwForceSwitching = ComputeNonbondedUtil::vdwForceSwitching;
-    const Bool LJcorrection = ComputeNonbondedUtil::LJcorrection;
     const Bool Fep_WCA_repuOn = ComputeNonbondedUtil::Fep_WCA_repuOn;
     const Bool Fep_WCA_dispOn = ComputeNonbondedUtil::Fep_WCA_dispOn;
     const Bool Fep_ElecOn = ComputeNonbondedUtil::Fep_ElecOn;
index 1416ae0..a66d3c9 100644 (file)
@@ -239,11 +239,11 @@ MODIFIED(
         FEP(fep_vdw_forceandenergies(A,B,r2,myVdwShift,myVdwShift2,switchdist2,
           cutoff2, myVdwLambda, myVdwLambda2, Fep_WCA_repuOn, Fep_WCA_dispOn,
           Fep_Wham, WCA_rcut1, WCA_rcut2, WCA_rcut3, switchfactor,
-          vdwForceSwitching, LJcorrection, &alch_vdw_energy, &alch_vdw_force, 
+          vdwForceSwitching, &alch_vdw_energy, &alch_vdw_force,
           &alch_vdw_energy_2, &alch_vdw_energy_2_Left);)
         TI(ti_vdw_force_energy_dUdl(A,B,r2,myVdwShift,switchdist2,
           cutoff2, myVdwLambda, alchVdwShiftCoeff, switchfactor,
-          vdwForceSwitching, LJcorrection, &alch_vdw_energy, &alch_vdw_force, 
+          vdwForceSwitching, &alch_vdw_energy, &alch_vdw_force,
           &alch_vdw_dUdl);)
       )
       
index 8f318fb..24bac70 100644 (file)
@@ -21,9 +21,8 @@ inline void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2,
   BigReal cutoff2, BigReal myVdwLambda, BigReal myVdwLambda2, 
   Bool Fep_WCA_repuOn, Bool Fep_WCA_dispOn, bool Fep_Wham, BigReal WCA_rcut1, 
   BigReal WCA_rcut2, BigReal WCA_rcut3, BigReal switchfactor, 
-  Bool vdwForceSwitching, Bool LJcorrection, BigReal* alch_vdw_energy,
-  BigReal* alch_vdw_force, BigReal* alch_vdw_energy_2,
-  BigReal* alch_vdw_energy_2_Left) {
+  Bool vdwForceSwitching, BigReal* alch_vdw_energy, BigReal* alch_vdw_force,
+  BigReal* alch_vdw_energy_2, BigReal* alch_vdw_energy_2_Left) {
   // switching function (this is correct whether switching is active or not)
   const BigReal switchmul = (r2 > switchdist2 ? switchfactor*(cutoff2 - r2) \
                             *(cutoff2 - r2)                            \
@@ -146,14 +145,6 @@ inline void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2,
           shifted_switchdist2_2*shifted_switchdist2_2*shifted_switchdist2_2;
       const BigReal shifted_switchdist3_2 = shifted_switchdist2_2*shifted_switchdist_2;
 
-      const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
-      const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
-      const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
-      const BigReal v_vdwa2 = -A / (cutoff6*shifted_switchdist6_2);
-      const BigReal v_vdwb2 = -B / (cutoff3*shifted_switchdist3_2);
-      const BigReal dU2 = v_vdwa2 - v_vdwb2; //deltaV2 from Steinbach & Brooks
-
       if(r2 > switchdist2) {
        const BigReal k_vdwa = A*cutoff6 / (cutoff6 - shifted_switchdist6);
        const BigReal k_vdwb = B*cutoff3 / (cutoff3 - shifted_switchdist3);
@@ -167,18 +158,22 @@ inline void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2,
         const BigReal tmpa2 = r6_2 - (1./cutoff6);
         const BigReal tmpb2 = r2_2*r_2 - (1./cutoff3);
 
-        *alch_vdw_energy = (myVdwLambda*(k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb
-                                         - (LJcorrection ? dU : 0.)));
+        *alch_vdw_energy = myVdwLambda*(k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb);
         *alch_vdw_energy_2 = (myVdwLambda2*(k_vdwa2*tmpa2*tmpa2 \
-                                            - k_vdwb2*tmpb2*tmpb2
-                                            - (LJcorrection ? dU2 : 0.)));
+                                            - k_vdwb2*tmpb2*tmpb2));
         *alch_vdw_force = (myVdwLambda*r2_1*(12.*k_vdwa*tmpa*r6_1
                                              - 6.*k_vdwb*tmpb*r2_1*r_1));
       }else{
-        if(!LJcorrection) {
-         *alch_vdw_energy += myVdwLambda*dU;
-         *alch_vdw_energy_2 += myVdwLambda2*dU2;
-        }
+        const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
+        const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
+        const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
+
+        const BigReal v_vdwa2 = -A / (cutoff6*shifted_switchdist6_2);
+        const BigReal v_vdwb2 = -B / (cutoff3*shifted_switchdist3_2);
+        const BigReal dU2 = v_vdwa2 - v_vdwb2; //deltaV2 from Steinbach & Brooks
+
+        *alch_vdw_energy += myVdwLambda*dU;
+        *alch_vdw_energy_2 += myVdwLambda2*dU2;
       }
     }
   }
index 89a8288..71b972e 100644 (file)
@@ -19,8 +19,8 @@
 inline void ti_vdw_force_energy_dUdl (BigReal A, BigReal B, BigReal r2, 
   BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2, 
   BigReal myVdwLambda, BigReal alchVdwShiftCoeff, BigReal switchfactor, 
-  Bool vdwForceSwitching, Bool LJcorrection, BigReal* alch_vdw_energy,
-  BigReal* alch_vdw_force, BigReal* alch_vdw_dUdl) {
+  Bool vdwForceSwitching, BigReal* alch_vdw_energy, BigReal* alch_vdw_force,
+  BigReal* alch_vdw_dUdl) {
   //myVdwShift already multplied by relevant (1-vdwLambda)
   const BigReal r2_1 = 1./(r2 + myVdwShift);
   const BigReal r6_1 = r2_1*r2_1*r2_1;
@@ -51,18 +51,13 @@ inline void ti_vdw_force_energy_dUdl (BigReal A, BigReal B, BigReal r2,
         shifted_switchdist2*shifted_switchdist2*shifted_switchdist2;
     const BigReal shifted_switchdist3 = shifted_switchdist2*shifted_switchdist;
 
-    const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
-    const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
-    const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
-
     if(r2 > switchdist2) {
       const BigReal k_vdwa = A*cutoff6 / (cutoff6 - shifted_switchdist6);
       const BigReal k_vdwb = B*cutoff3 / (cutoff3 - shifted_switchdist3);
       const BigReal tmpa = r6_1 - (1./cutoff6);
       const BigReal r_1 = sqrt(r2_1);
       const BigReal tmpb = r2_1*r_1 - (1./cutoff3);
-      const BigReal Uh = (k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb \
-                          - (LJcorrection ? dU : 0.));
+      const BigReal Uh = k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb;
       *alch_vdw_energy = myVdwLambda*Uh;
       *alch_vdw_force = (myVdwLambda*r2_1*(12.*k_vdwa*tmpa*r6_1
                                            - 6.*k_vdwb*tmpb*r2_1*r_1));
@@ -70,10 +65,12 @@ inline void ti_vdw_force_energy_dUdl (BigReal A, BigReal B, BigReal r2,
                         *(6.*k_vdwa*tmpa*r6_1                   \
                           - 3.*k_vdwb*tmpb*r2_1*r_1));
     }else{
-      if(!LJcorrection) {
-        *alch_vdw_energy += myVdwLambda*dU;
-        *alch_vdw_dUdl += dU;
-      }
+      const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
+      const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
+      const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
+
+      *alch_vdw_energy += myVdwLambda*dU;
+      *alch_vdw_dUdl += dU;
     }
   }
 }
index d84b830..3356c8b 100644 (file)
@@ -125,7 +125,6 @@ BigReal   ComputeNonbondedUtil::alchDispLambda;
 BigReal   ComputeNonbondedUtil::alchElecLambda;
 BigReal   ComputeNonbondedUtil::alchVdwShiftCoeff;
 Bool      ComputeNonbondedUtil::vdwForceSwitching;
-Bool      ComputeNonbondedUtil::LJcorrection;
 Bool      ComputeNonbondedUtil::alchDecouple;
 //fepe
 Bool      ComputeNonbondedUtil::lesOn;
@@ -304,7 +303,6 @@ void ComputeNonbondedUtil::select(void)
   Bool tabulatedEnergies = simParams->tabulatedEnergies;
   alchVdwShiftCoeff = simParams->alchVdwShiftCoeff;
   vdwForceSwitching = simParams->vdwForceSwitching;
-  LJcorrection = simParams->LJcorrection;
   WCA_rcut1 = simParams->alchFepWCArcut1;
   WCA_rcut2 = simParams->alchFepWCArcut2;
   WCA_rcut3 = simParams->alchFepWCArcut3;
@@ -867,14 +865,14 @@ void ComputeNonbondedUtil::select(void)
     // from Steinbach & Brooks, JCC 15, pgs 667-683, 1994, eqns 10-13
     if ( r2 > switchOn2 ) {
       BigReal tmpa = r_6 - cutoff_6;
-      vdwa_energy = k_vdwa * tmpa * tmpa - (LJcorrection ? v_vdwa : 0.0);
+      vdwa_energy = k_vdwa * tmpa * tmpa;
       BigReal tmpb = r_1 * r_2 - cutoff_3;
-      vdwb_energy = k_vdwb * tmpb * tmpb - (LJcorrection ? v_vdwb : 0.0);
+      vdwb_energy = k_vdwb * tmpb * tmpb;
       vdwa_gradient = -6.0 * k_vdwa * tmpa * r_2 * r_6;
       vdwb_gradient = -3.0 * k_vdwb * tmpb * r_2 * r_2 * r_1;
     } else {
-      vdwa_energy = r_12 + (LJcorrection ? 0.0 : v_vdwa);
-      vdwb_energy = r_6 + (LJcorrection ? 0.0 : v_vdwb);
+      vdwa_energy = r_12 + v_vdwa;
+      vdwb_energy = r_6 + v_vdwb;
       vdwa_gradient = -6.0 * r_2 * r_12;
       vdwb_gradient = -3.0 * r_2 * r_6;
     }
index 228578b..91a0238 100644 (file)
@@ -357,7 +357,6 @@ public:
   static Bool alchThermIntOn;
   static BigReal alchVdwShiftCoeff;
   static Bool vdwForceSwitching;
-  static Bool LJcorrection;
   static Bool Fep_WCA_repuOn;
   static Bool Fep_WCA_dispOn;
   static Bool Fep_ElecOn;
index c245acd..0adb8ef 100644 (file)
@@ -10057,26 +10057,49 @@ void Molecule::compute_LJcorrection() {
   BigReal rswitch5 = rswitch2*rswitch3;
   BigReal rswitch6 = rswitch3*rswitch3;
 
-  /* Here we tabulate the integrals over the untruncated region. This assumes:
-
-   1.) The energy and virial contribution can be well described by a mean field
-       approximation (i.e. a constant).
-
-   2.) The pair distribution function, g(r), is very close to unity on the
-       interval (i.e. g(r) = 1 for r > switchdist).
-
-   The mean field integrals are of the form:
-
-   4*N^2*PI*int r^2 U(r) dr : for the energy
-
-   (4/3)*N^2*PI*int r^3 dU(r)/dr dr : for the virial
-
-   NB: An extra factor of 1/2 comes from double counting the number of
-       interaction pairs (N*(N-1)/2, approximated as N^2).
-  */
+  /*
+   * Tabulate the integrals over the untruncated region. This assumes:
+   *
+   * 1.) The energy and virial contribution can be well described by a mean
+   *     field approximation (i.e. a constant).
+   *
+   * 2.) The radial distribution function, g(r), is very close to unity on the
+   *     interval (i.e. g(r) = 1 for r > rswitch).
+   *
+   * The mean field integrals are, for the energy (int_U_gofr):
+   *
+   * 2\pi \int_0^{\infty} dr r^2 (1 - S(r; r_s, r_c)) U(r)
+   *
+   * and for the virial (int_rF_gofr):
+   *
+   * -\frac{2\pi}{3) \int_0^{\infty}
+   *     dr r^3 \frac{d}{dr} [(1 - S(r; r_s, r_c)) U(r)]
+   *
+   * Here U(r) is the "mean" LJ-potential and S(r; r_s, r_c) is the switching
+   * function (parameterized by r_s = rswitch and r_c = rcut). These equations
+   * include all factors except NumLJSites^2 and the volume. Because
+   * NumLJSites^2 derives from the approximation N*(N-1)/2 ~= N^2/2 for the
+   * number of interaction pairs, there is an "extra" factor of 1/2.
+   */
   BigReal int_U_gofr_A, int_rF_gofr_A, int_U_gofr_B, int_rF_gofr_B;
   if (simParams->switchingActive) {
     if (!simParams->vdwForceSwitching) {
+      /*
+       * S(r; r_s, r_c)
+       * =
+       * \begin{cases}
+       *   1 & 0 \le r \le r_s
+       *   \\
+       *   \frac{
+       *     (r_c^2 - r^2)^2 (r_c^2 - 3r_s^2 + 2r^2)
+       *   }{
+       *     (r_c^2 - r_s^2)^3
+       *   } 
+       *   & r_s < r < r_c
+       *   \\
+       *   0 & r_c \le r < \infty
+       * \end{cases}
+       */
       BigReal rsum3 = (rcut + rswitch)*(rcut + rswitch)*(rcut + rswitch);
       int_U_gofr_A = int_rF_gofr_A = (16*PI*(3*rcut4 + 9*rcut3*rswitch
                                       + 11*rcut2*rswitch2 + 9*rcut*rswitch3
@@ -10085,28 +10108,93 @@ void Molecule::compute_LJcorrection() {
       int_U_gofr_B = int_rF_gofr_B = -16*PI / (3*rsum3);
     }
     else {
-      /* BKR - This assumes a definition of the vdwForceSwitching potential
-       that differs from the original by Steinbach and Brooks. Here the
-       potential shift is _subtracted_ from the switching region instead of
-       being _added_ to the region inside switchdist (that would require a
-       short-range correction with no obviously good approximation). The new
-       potential is discontinuous at the cutoff, since it is rigorously zero
-       beyond that point but finite and negative inside it. Nonetheless, the
-       force is continuous, so this is ok.
-
-       NB: Because the difference is just a constant, the virial correction is
-       the same, regardless which definition is chosen.
+      /* BKR - There are two choices for the force switching strategy:
+
+         1) apply a potential shift for 0 <= r <= rswitch (standard)
+         or
+         2) apply the opposite shift for rswitch < r < rcut
+
+         Option 2 was previously implemented, but introduces a potential
+         discontinuity at rcut and thus still causes energy conservation
+         issues. The energy correction for option 1, on the other hand,
+         requires the dubious approximation that g(r) ~= 1 for
+         0 <= r <= rswitch. However, this approximation only needs to hold in
+         so far as the integral out to rswitch is roughly the same -- this is
+         actually sufficiently close in practice. Both options lead to the same
+         virial correction.
+
+         From Steinbach and Brooks:
+
+         U_h(r; r_s, r_c)
+         =
+         \frac{C_n r_c^{n/2}}{r_c^{n/2} - r_s^{n/2}}
+           \left( \frac{1}{r^{n/2}} - \frac{1}{r_c^{n/2}} \right)^2
+
+         \Delta U(r_s, r_c) = -\frac{C_n}{(r_s r_c)^{n/2}}
       */
       BigReal lnr = log(rcut/rswitch);
-      int_rF_gofr_A = 16*PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
-      int_rF_gofr_B = -4*PI*lnr / (rcut3 - rswitch3);
+      /*
+       * Option 1 (shift below rswitch)
+       *
+       * S(r; r_s, r_c)
+       * =
+       * \begin{cases}
+       *   \frac{
+       *     U(r) + \Delta U(r_s, r_c)
+       *   }{
+       *     U(r)
+       *   }
+       *   & 0 \le r \le r_s
+       *   \\
+       *   \frac{
+       *     U_h(r; r_s, r_c)  
+       *   }{
+       *     U(r)
+       *   }
+       *   & r_s < r < r_c
+       *   \\
+       *   0 & r_c \le r < \infty
+       * \end{cases}
+       */
+      int_U_gofr_A = int_rF_gofr_A =\
+          16*PI / (9*rswitch3*rcut3*(rcut3 + rswitch3));
+      int_U_gofr_B = int_rF_gofr_B = -4*PI*lnr / (rcut3 - rswitch3);
+      /*
+       * Option 2 (shift above rswitch and below rcut)
+       *
+       * S(r; r_s, r_c)
+       * =
+       * \begin{cases}
+       *   1 & 0 \le r \le r_s
+       *   \\
+       *   \frac{
+       *     U_h(r; r_s, r_c) - \Delta U(r_s, r_c)
+       *   }{
+       *     U(r)
+       *   }
+       *   & r_s < r < r_c
+       *   \\
+       *   0 & r_c \le r < \infty
+       * \end{cases}
+       */ 
+/*
       int_U_gofr_A = (2*PI*(5*rswitch3 - 3*rcut3)
                       / (9*rcut3*rswitch6*(rcut3 + rswitch3)));
       int_U_gofr_B = (-2*PI*(rswitch3 - rcut3 + 6*rswitch3*lnr)
                       / (3*rswitch3*(rcut3 - rswitch3)));
+*/
     }
   }
   else {
+    /*
+     * S(r; r_s, r_c)
+     * =
+     * \begin{cases}
+     *   1 & 0 \le r \le r_c
+     *   \\
+     *   0 & r_c \le r < \infty
+     * \end{cases}
+     */
     int_rF_gofr_A = 8*PI / (9*rcut9);
     int_rF_gofr_B = -4*PI / (3*rcut3);
     int_U_gofr_A = 2*PI / (9*rcut9);