Optimize alchemical scaling schedule for LJcorrection 61/4961/3
authorradakb <brian.radak@gmail.com>
Tue, 19 Feb 2019 15:27:48 +0000 (10:27 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Fri, 29 Mar 2019 17:57:31 +0000 (12:57 -0500)
The long-range Lennard-Jones correction is included and scaled
w/VDW energy during alchemical simulations. This change revises
the coupling schedule when alchVdwLambdaEnd != 0.0 or 1.0.

The old behavior was to compute the full correction at both
alchemical endpoints and then completely decouple one before
recoupling the other - this is extremely dramatic in the rare, but
reasonable, instance that the coupling schedule is not symmetric.
The new behavior always retains the contribution from non-alchemical
atoms and only couples/decouples the "correction" for the alchemical
atoms.

The TI derivatives are now also "cleaned up" so that there are no
spurious components coming from the opposing alchemical endpoint.
For example, previously, there were non-zero contributions in VDW2
that exactly cancelled part of VDW1, even if alchemical group 2 is
empty. The result is the same, but the new code is more intuitive.

Finally, this enables proper all-in-one decomposition of free energies
into VDW and ELECT components for arbitrary values of alchVdwLambdaEnd.
The old schedule could inadvertently contaminate the ELECT portion
with the LJcorrection unless alchVdwLambdaEnd was properly chosen.

Change-Id: I3ae5903434e68fb134fe24952bb35a175f9d0621

src/Controller.C
src/Molecule.C
src/Molecule.h

index af0a937..f2b4726 100644 (file)
@@ -3052,15 +3052,15 @@ void Controller::printEnergies(int step, int minimize)
 #else
       // Apply tail correction to energy.
       BigReal alchLambda = simParameters->getCurrentLambda(step);
-      ljEnergy += molecule->getEnergyTailCorr(alchLambda) / volume;
+      ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
 
       if (simParameters->alchOn) {
         if (simParameters->alchFepOn) {
           BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
-          ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2) / volume;
+          ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
         } else if (simParameters->alchThermIntOn) {
-          ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0) / volume;
-          ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0) / volume;
+          ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
+          ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
         }
       }
 #endif
index 0adb8ef..01c4867 100644 (file)
@@ -9926,7 +9926,7 @@ void Molecule::compute_LJcorrection() {
     alternative would be to count "fractional interactions," but that makes
     TI derivatives a bit harder and for no obvious gain.
   */
-  if (simParams->alchOn && simParams->alchVdwLambdaEnd) {
+  if (simParams->alchOn && simParams->alchVdwLambdaEnd > 0.0) {
     BigReal sumOfAs1 = sumOfAs;
     BigReal sumOfAs2 = sumOfAs;
     BigReal sumOfBs1 = sumOfBs;
@@ -9981,7 +9981,10 @@ void Molecule::compute_LJcorrection() {
           B = 4.0*sigma_ij*epsilon_ij;
         }
         if (!A && !B) continue; // don't count zeroed interactions
-
+        // remove all alchemical interactions from group 0
+        sumOfAs -= 2*A;
+        sumOfBs -= 2*B;
+        count -= 2;
         if ( alchFlagSum > 0 ){ // in group 1, remove from group 2
           sumOfAs2 -= 2*A;
           sumOfBs2 -= 2*B;
@@ -10006,17 +10009,21 @@ void Molecule::compute_LJcorrection() {
       if ( alchFlagi == 1 || alchFlagi == -1 ) alch_counter++;
       if ( alch_counter == (numFepInitial + numFepFinal) ) break;
     }
-    if ( count1 ) { // Avoid divide by zero for empty groups.
+    LJAvgA = sumOfAs / count;
+    LJAvgB = sumOfBs / count;
+    if ( count1 ) {
       LJAvgA1 = sumOfAs1 / count1;
       LJAvgB1 = sumOfBs1 / count1;
-    } else { // Avoid trailing decimals due to finite precision.
-      LJAvgA1 = LJAvgB1 = 0.0;
+    } else { // This endpoint is only non-alchemical atoms.
+      LJAvgA1 = LJAvgA;
+      LJAvgB1 = LJAvgB;
     }
     if ( count2 ) {
       LJAvgA2 = sumOfAs2 / count2;
       LJAvgB2 = sumOfBs2 / count2;
-    } else {
-      LJAvgA2 = LJAvgB2 = 0.0;
+    } else { // This endpoint is only non-alchemical atoms.
+      LJAvgA2 = LJAvgA;
+      LJAvgB2 = LJAvgB;
     }
     if ( ! CkMyPe() ) {
       iout << iINFO << "LONG-RANGE LJ: APPLYING ANALYTICAL CORRECTIONS TO "
@@ -10026,6 +10033,7 @@ void Molecule::compute_LJcorrection() {
       iout << iINFO << "LONG-RANGE LJ: AVERAGE A1 AND B1 COEFFICIENTS "
            << LJAvgA1 << " AND " << LJAvgB1 << "\n" << endi;
     }
+    numLJsites = (numLJsites1 + numLJsites2 - numLJsites);
     LJAvgA1 *= BigReal(numLJsites1)*numLJsites1; 
     LJAvgB1 *= BigReal(numLJsites1)*numLJsites1;
     LJAvgA2 *= BigReal(numLJsites2)*numLJsites2;
@@ -10200,23 +10208,31 @@ void Molecule::compute_LJcorrection() {
     int_U_gofr_A = 2*PI / (9*rcut9);
     int_U_gofr_B = -2*PI / (3*rcut3);
   }
-  // For alchOn, these represent all atoms, with no sorting by group.
+  // For alchOn, these are the averages for all non-alchemical atoms.
   tail_corr_virial = int_rF_gofr_A*LJAvgA + int_rF_gofr_B*LJAvgB;
   tail_corr_ener = int_U_gofr_A*LJAvgA + int_U_gofr_B*LJAvgB;
 
-  tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1;
-  tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1;
-  tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2;
-  tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2;
+  tail_corr_dUdl_1 = int_U_gofr_A*LJAvgA1 + int_U_gofr_B*LJAvgB1 -
+      tail_corr_ener;
+  tail_corr_virial_1 = int_rF_gofr_A*LJAvgA1 + int_rF_gofr_B*LJAvgB1 -
+      tail_corr_virial;
+  tail_corr_dUdl_2 = int_U_gofr_A*LJAvgA2 + int_U_gofr_B*LJAvgB2 -
+      tail_corr_ener;
+  tail_corr_virial_2 = int_rF_gofr_A*LJAvgA2 + int_rF_gofr_B*LJAvgB2 -
+      tail_corr_virial;
 }
 
 // Convenience function to simplify lambda scaling.
-BigReal Molecule::getEnergyTailCorr(const BigReal alchLambda){
+// Setting doti TRUE and alchLambda = 0 or 1 will return the thermodynamic
+// derivative at the corresponding endpoint.
+BigReal Molecule::getEnergyTailCorr(const BigReal alchLambda, const int doti){
   const BigReal scl = simParams->nonbondedScaling;
-  if (simParams->alchOn && simParams->alchVdwLambdaEnd) {
+  if (simParams->alchOn && simParams->alchVdwLambdaEnd > 0.0) {
     const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
     const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
-    return scl*(vdw_lambda_1*tail_corr_dUdl_1 + vdw_lambda_2*tail_corr_dUdl_2);
+    const BigReal corr = (doti ? 0.0 : tail_corr_ener);
+    return scl*(corr + vdw_lambda_1*tail_corr_dUdl_1 +
+                       vdw_lambda_2*tail_corr_dUdl_2);
   }
   else {
     return scl*tail_corr_ener;
@@ -10226,11 +10242,11 @@ BigReal Molecule::getEnergyTailCorr(const BigReal alchLambda){
 // Convenience function to simplify lambda scaling.
 BigReal Molecule::getVirialTailCorr(const BigReal alchLambda){
   const BigReal scl = simParams->nonbondedScaling;
-  if (simParams->alchOn && simParams->alchVdwLambdaEnd) {
+  if (simParams->alchOn && simParams->alchVdwLambdaEnd > 0.0) {
     const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
     const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
-    return scl*(vdw_lambda_1*tail_corr_virial_1 +
-                vdw_lambda_2*tail_corr_virial_2);
+    return scl*(tail_corr_virial + vdw_lambda_1*tail_corr_virial_1 +
+                                   vdw_lambda_2*tail_corr_virial_2);
   }
   else {
     return scl*tail_corr_virial;
index 638ba33..02d636e 100644 (file)
@@ -470,7 +470,7 @@ public:
   BigReal tail_corr_dUdl_1, tail_corr_dUdl_2;
   BigReal tail_corr_virial_1, tail_corr_virial_2;
   void compute_LJcorrection();
-  BigReal getEnergyTailCorr(const BigReal);
+  BigReal getEnergyTailCorr(const BigReal, const int);
   BigReal getVirialTailCorr(const BigReal);
 
   int const * getLcpoParamType() {