Revamp of alchemical WCA decomposition 18/5018/2
authorradakb <brian.radak@gmail.com>
Wed, 13 Mar 2019 19:35:38 +0000 (15:35 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Mon, 1 Apr 2019 17:39:51 +0000 (12:39 -0500)
The old, undocumented, WCA strategy for LJ interactions is now
completely re-implemented and better integrated. All of the old
keywords are removed as well as all as the redundnant modifications
to electrostatics. The new code is compatible with TI (well, almost,
see below), IDWS, and all switching and LJcorrection schemes. There
re now two new WCA specific keywords:

alchWCA -- boolean for the WCA code path
alchRepLambdaEnd -- new lambda scheduling parameter

The WCA code separates the LJ terms into "repulsive" and
"dispersive" parts. For an exnihilated particle, the repulsion
is now turned on between 0 and alchRepLambdaEnd and the dispersion
is turned on between alchRepLambdaEnd and alchVdwLambdaEnd. By
construction, the LJcorrection and switching options only affect
the alchemical atoms when alchLambda >= alchRepLambdaEnd.

A major change from the old code is that electrostatic scaling is
completely decoupled from the LJ scaling.

Tests over a large set of molecules show that free energies agree
with the conventional soft-core routines, as expected.

TODO:

- There are now a LOT of if/else statements in the alchemical LJ
  routines. Many of these are global booleans and can probably be
  optimized (all combinations of alchWCA, switching, and
  vdwForceSwitching lead to globally different paths).

- Many of the generic inline routines in ComputeNonbondedAlch.h
  could probably be optimized.

- WCA has not historically been used with TI, but there is no
  theoretical limitation here. However, the additional decomposition
  implies TWO derivatives, although these are only both needed at
  alchLambda == alchRepLambdaEnd. This is not correctly implemented
  as is.

Change-Id: Ia3eef1633dc6dd21ffe0101428e2fb219d8037be

src/ComputeNonbondedAlch.h [new file with mode: 0644]
src/ComputeNonbondedBase.h
src/ComputeNonbondedBase2.h
src/ComputeNonbondedFEP.C
src/ComputeNonbondedTI.C
src/ComputeNonbondedUtil.C
src/ComputeNonbondedUtil.h
src/ComputePme.C
src/Controller.C
src/SimParameters.C
src/SimParameters.h

diff --git a/src/ComputeNonbondedAlch.h b/src/ComputeNonbondedAlch.h
new file mode 100644 (file)
index 0000000..db85505
--- /dev/null
@@ -0,0 +1,106 @@
+/**
+***  Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
+***  The Board of Trustees of the University of Illinois.
+***  All rights reserved.
+**/
+
+/*****************************************************************************
+ *  Inline Lennard-Jones energy/force/utility functions for alchemy          *
+ *****************************************************************************/
+
+/* Lennard-Jones force and energy.
+ *
+ * If using conventional soft-core, r2 should be passed as r2 + myVdwShift and
+ * the final result should be scaled by myVdwLambda.
+ */
+inline void vdw_forceandenergy(const BigReal A, const BigReal B,
+    const BigReal r2, BigReal *U, BigReal *F) {
+  const BigReal inv_r2 = 1. / r2;
+  const BigReal inv_r6 = inv_r2*inv_r2*inv_r2;
+  *U = inv_r6*(A*inv_r6 - B);
+  *F = 6*inv_r2*(2*(*U) + B*inv_r6);
+}
+
+// Same as above, energy only.
+inline void vdw_energy(const BigReal A, const BigReal B,
+    const BigReal r2, BigReal *U) {
+  const BigReal inv_r2 = 1. / r2;
+  const BigReal inv_r6 = inv_r2*inv_r2*inv_r2;
+  *U = inv_r6*(A*inv_r6 - B);
+}
+
+/* Multipliers for the vdW potential switching function and its derivative.
+ * This is correct even if switching is not active.
+ */
+inline void vdw_switch(const BigReal r2, const BigReal switchdist2,
+    const BigReal cutoff2, const BigReal switchfactor, BigReal* switchmul,
+    BigReal* switchmul2) {
+  if (r2 <= switchdist2) {
+    // "switching off" will always fall through here.
+    *switchmul = 1.0;
+    *switchmul2 = 0.0;
+  } else {
+    const BigReal cdiff2 = (cutoff2 - r2);
+    const BigReal sdiff2 = (r2 - switchdist2);
+    *switchmul = switchfactor*cdiff2*cdiff2*(cdiff2 + 3*sdiff2);
+    *switchmul2 = 12*switchfactor*cdiff2*sdiff2;
+  }
+}
+
+/* Force shifted Lennard-Jones force and energy.
+ *
+ * If using conventional soft-core, r2 and switchdist2 should be passed as
+ * r2 + myVdwShift and switchdist2 + myVdwShift and the final result should be
+ * scaled by myVdwLambda
+ */
+inline void vdw_fswitch_forceandenergy(const BigReal A, const BigReal B,
+    const BigReal r2, const BigReal switchdist2, const BigReal cutoff2,
+    BigReal* U, BigReal* F) {
+  // TODO: Some rearrangement could be done to use rsqrt() instead?
+  const BigReal inv_r2 = 1. / r2;
+  const BigReal inv_r6 = inv_r2*inv_r2*inv_r2;
+  const BigReal inv_r3 = sqrt(inv_r6);
+  const BigReal inv_cutoff6 = 1. / (cutoff2*cutoff2*cutoff2);
+  const BigReal inv_cutoff3 = sqrt(inv_cutoff6);
+  const BigReal switchdist6 = switchdist2*switchdist2*switchdist2;
+  const BigReal k_vdwa = A / (1 - switchdist6*inv_cutoff6);
+  const BigReal k_vdwb = B / (1 - sqrt(switchdist6*inv_cutoff6));
+  const BigReal tmpa = inv_r6 - inv_cutoff6;
+  const BigReal tmpb = inv_r3 - inv_cutoff3;
+  *U = k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb;
+  *F = 6*inv_r2*(2*k_vdwa*tmpa*inv_r6 - k_vdwb*tmpb*inv_r3);
+}
+
+// Same as above, energy only.
+inline void vdw_fswitch_energy(const BigReal A, const BigReal B,
+    const BigReal r2, const BigReal switchdist2, const BigReal cutoff2,
+    BigReal* U) {  
+  // TODO: Some rearrangement could be done to use rsqrt() instead?
+  const BigReal inv_r2 = 1. / r2;
+  const BigReal inv_r6 = inv_r2*inv_r2*inv_r2;
+  const BigReal inv_r3 = sqrt(inv_r6);
+  const BigReal inv_cutoff6 = 1. / (cutoff2*cutoff2*cutoff2);
+  const BigReal inv_cutoff3 = sqrt(inv_cutoff6);
+  const BigReal switchdist6 = switchdist2*switchdist2*switchdist2;
+  const BigReal k_vdwa = A / (1 - switchdist6*inv_cutoff6);
+  const BigReal k_vdwb = B / (1 - sqrt(switchdist6*inv_cutoff6));
+  const BigReal tmpa = inv_r6 - inv_cutoff6;
+  const BigReal tmpb = inv_r3 - inv_cutoff3;
+  *U = k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb;
+}
+
+/* Energy shifts for the vdW force/potential switching functions.
+ *
+ * If using conventional soft-core, switchdist2 should be passed as
+ * switchdist2 + myVdwShift.
+ */
+inline void vdw_fswitch_shift(const BigReal A, const BigReal B,
+    const BigReal switchdist2, const BigReal cutoff2, BigReal* dU) {
+  // TODO: Some rearrangement could be done to use rsqrt() instead.
+  const BigReal cutoff6 = cutoff2*cutoff2*cutoff2;
+  const BigReal switchdist6 = switchdist2*switchdist2*switchdist2;
+  const BigReal v_vdwa = -A / (cutoff6*switchdist6);
+  const BigReal v_vdwb = -B / sqrt(cutoff6*switchdist6);
+  *dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
+}
+
index 83244c6..7acd5f4 100644 (file)
@@ -294,7 +294,7 @@ void ComputeNonbondedUtil :: NAME
 
   FEP
   (
-  ENERGY( BigReal vdwEnergy_s = 0; BigReal vdwEnergy_s_Left = 0; )
+  ENERGY( BigReal vdwEnergy_s = 0; )
   SHORT
   (
   ENERGY( BigReal electEnergy_s = 0; )
@@ -380,16 +380,11 @@ void ComputeNonbondedUtil :: NAME
   ALCH(
     const BigReal switchdist2 = ComputeNonbondedUtil::switchOn2;
     const BigReal cutoff2 = ComputeNonbondedUtil::cutoff2;
-    const BigReal switchfactor = 1./((cutoff2 - switchdist2)*(cutoff2 - switchdist2)*(cutoff2 - switchdist2));
+    const BigReal diff = cutoff2 - switchdist2;
+    const BigReal switchfactor = 1./(diff*diff*diff);
     const BigReal alchVdwShiftCoeff = ComputeNonbondedUtil::alchVdwShiftCoeff;
+    const Bool alchWCAOn = ComputeNonbondedUtil::alchWCAOn;
     const Bool vdwForceSwitching = ComputeNonbondedUtil::vdwForceSwitching;
-    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;
-    const Bool Fep_Wham = ComputeNonbondedUtil::Fep_Wham;
-    const BigReal WCA_rcut1 = ComputeNonbondedUtil::WCA_rcut1;
-    const BigReal WCA_rcut2 = ComputeNonbondedUtil::WCA_rcut2;
-    const BigReal WCA_rcut3 = ComputeNonbondedUtil::WCA_rcut3;
 
     /* BKR - Enable dynamic lambda by passing the timestep to simParams.
        The SimParameters object now has all knowledge of the lambda schedule
@@ -402,53 +397,25 @@ void ComputeNonbondedUtil :: NAME
     BigReal lambdaUp = alchLambda; // needed in ComputeNonbondedBase2.h!
     BigReal elecLambdaUp = simParams->getElecLambda(lambdaUp);
     BigReal vdwLambdaUp = simParams->getVdwLambda(lambdaUp);
+    BigReal repLambdaUp = simParams->getRepLambda(lambdaUp);
     BigReal vdwShiftUp = alchVdwShiftCoeff*(1 - vdwLambdaUp);
     FEP(BigReal lambda2Up = simParams->getCurrentLambda2(params->step);)
     FEP(BigReal elecLambda2Up = simParams->getElecLambda(lambda2Up);)
     FEP(BigReal vdwLambda2Up = simParams->getVdwLambda(lambda2Up);)
+    FEP(BigReal repLambda2Up = simParams->getRepLambda(lambda2Up);)
     FEP(BigReal vdwShift2Up = alchVdwShiftCoeff*(1 - vdwLambda2Up);)
 
-    FEP( if( (Fep_Wham) && (Fep_WCA_repuOn) ) {
-       elecLambdaUp=0.0; 
-       vdwLambdaUp=ComputeNonbondedUtil::alchRepLambda;
-    })
-    FEP( if( (Fep_Wham) && (Fep_WCA_dispOn) ) {
-       elecLambdaUp=0.0; 
-       vdwLambdaUp=ComputeNonbondedUtil::alchDispLambda;
-    })
-    FEP( if( (Fep_Wham) && (Fep_ElecOn) ) {
-       elecLambdaUp=ComputeNonbondedUtil::alchElecLambda; 
-       vdwLambdaUp=1.0;
-       vdwLambda2Up=1.0;
-       vdwShiftUp = 0.0;
-       vdwShift2Up = 0.0;
-    })
-        
     /*lambda values 'down' are for atoms scaled down with lambda (partition 2)*/
     BigReal lambdaDown = 1 - alchLambda; // needed in ComputeNonbondedBase2.h!
     BigReal elecLambdaDown = simParams->getElecLambda(lambdaDown);
     BigReal vdwLambdaDown = simParams->getVdwLambda(lambdaDown);
+    BigReal repLambdaDown = simParams->getRepLambda(lambdaDown);
     BigReal vdwShiftDown = alchVdwShiftCoeff*(1 - vdwLambdaDown);
     FEP(BigReal lambda2Down = 1 - simParams->getCurrentLambda2(params->step);)
     FEP(BigReal elecLambda2Down = simParams->getElecLambda(lambda2Down);)
     FEP(BigReal vdwLambda2Down = simParams->getVdwLambda(lambda2Down);)
+    FEP(BigReal repLambda2Down = simParams->getRepLambda(lambda2Down);)
     FEP(BigReal vdwShift2Down = alchVdwShiftCoeff*(1 - vdwLambda2Down);)
-
-    FEP( if( (Fep_Wham) && (Fep_WCA_repuOn) ) {
-       elecLambdaDown=0.0; 
-       vdwLambdaDown = 1.0 - ComputeNonbondedUtil::alchRepLambda;
-    })
-    FEP( if( (Fep_Wham) && (Fep_WCA_dispOn) ) {
-       elecLambdaDown=0.0; 
-       vdwLambdaDown = 1.0 - ComputeNonbondedUtil::alchDispLambda;
-    })
-    FEP( if( (Fep_Wham) && (Fep_ElecOn) ) {
-       elecLambdaDown = 1.0 - ComputeNonbondedUtil::alchElecLambda; 
-       vdwLambdaDown = 1.0;
-       vdwLambda2Down = 1.0;
-       vdwShiftDown = 0.0;
-       vdwShift2Down = 0.0;
-    })
   
   // Thermodynamic Integration Notes: 
   // Separation of atom pairs into different pairlist according to lambda
@@ -2038,9 +2005,10 @@ void ComputeNonbondedUtil :: NAME
     BigReal myLambda; FEP(BigReal myLambda2;)
     BigReal myElecLambda;  FEP(BigReal myElecLambda2;)
     BigReal myVdwLambda; FEP(BigReal myVdwLambda2;)
+    BigReal myRepLambda; FEP(BigReal myRepLambda2;)
     BigReal myVdwShift; FEP(BigReal myVdwShift2;)
     BigReal alch_vdw_energy; BigReal alch_vdw_force;
-    FEP(BigReal alch_vdw_energy_2; BigReal alch_vdw_energy_2_Left;) TI(BigReal alch_vdw_dUdl;)
+    FEP(BigReal alch_vdw_energy_2; ) TI(BigReal alch_vdw_dUdl;)
     BigReal shiftedElec; BigReal shiftedElecForce;
     
     /********************************************************************/
@@ -2267,7 +2235,6 @@ PAIR(
     TI(reduction[vdwEnergyIndex_ti_1] += vdwEnergy_ti_1;) 
     TI(reduction[vdwEnergyIndex_ti_2] += vdwEnergy_ti_2;) 
     FEP( reduction[vdwEnergyIndex_s] += vdwEnergy_s; )
-    FEP( reduction[vdwEnergyIndex_s_Left] += vdwEnergy_s_Left; )
   SHORT
   (
     FEP( reduction[electEnergyIndex_s] += electEnergy_s; )
index a66d3c9..88fd2da 100644 (file)
@@ -18,6 +18,8 @@ ALCHPAIR(
   FEP(myElecLambda2 = ALCH1(elecLambda2Up) ALCH2(elecLambda2Down);)
   myVdwLambda =  ALCH1(vdwLambdaUp) ALCH2(vdwLambdaDown); 
   FEP(myVdwLambda2 = ALCH1(vdwLambda2Up) ALCH2(vdwLambda2Down);) 
+  myRepLambda =  ALCH1(repLambdaUp) ALCH2(repLambdaDown);
+  FEP(myRepLambda2 = ALCH1(repLambda2Up) ALCH2(repLambda2Down);)
   myVdwShift =  ALCH1(vdwShiftUp) ALCH2(vdwShiftDown); 
   FEP(myVdwShift2 =  ALCH1(vdwShift2Up) ALCH2(vdwShift2Down);) 
 )
@@ -225,10 +227,10 @@ MODIFIED(
         // Alchemical free energy calculation
         // Pairlists are separated so that lambda-coupled pairs are handled
         // independently from normal nonbonded (inside ALCHPAIR macro).
-        // The separation-shifted van der Waals potential and a shifted 
-        // electrostatics potential for decoupling are calculated explicitly.
-        // Would be faster with lookup tables but because only a small minority
-        // of nonbonded pairs are lambda-coupled the impact is minimal. 
+        // The separation-shifted van der Waals potential for decoupling is
+        // calculated explicitly. This would be faster with lookup tables but
+        // because only a small minority of nonbonded pairs are lambda-coupled
+        // the impact is minimal. 
         // Explicit calculation also makes things easier to modify.
         
         const BigReal r2 = r2list[k] - r2_delta;
@@ -236,14 +238,13 @@ MODIFIED(
         // These are now inline functions (in ComputeNonbondedFep.C) to 
         // tidy the code
 
-        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, &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, &alch_vdw_energy, &alch_vdw_force,
+        FEP(fep_vdw_forceandenergies(A, B, r2, myVdwShift, myVdwShift2,
+          switchdist2, cutoff2, switchfactor, vdwForceSwitching, myVdwLambda,
+          myVdwLambda2, alchWCAOn, myRepLambda, myRepLambda2, &alch_vdw_energy,
+          &alch_vdw_force, &alch_vdw_energy_2);)
+        TI(ti_vdw_force_energy_dUdl(A, B, r2, myVdwShift, switchdist2, cutoff2,
+          switchfactor, vdwForceSwitching, myVdwLambda, alchVdwShiftCoeff,
+          alchWCAOn, myRepLambda, &alch_vdw_energy, &alch_vdw_force,
           &alch_vdw_dUdl);)
       )
       
@@ -275,7 +276,6 @@ MODIFIED(
             //CkPrintf("Found vdw energy of %f\n", vdw_val);
             vdwEnergy += LAM(lambda_pair *) vdw_val;
             FEP( vdwEnergy_s += d_lambda_pair * vdw_val; )
-            FEP( vdwEnergy_s_Left += d_lambda_pair * vdw_val; )
           )
         }  else {
          //)
@@ -287,7 +287,6 @@ MODIFIED(
         vdwEnergy -= LAM(lambda_pair *) vdw_val;
 
         FEP(vdwEnergy_s -= vdw_val;)
-        FEP(vdwEnergy_s_Left -= vdw_val;)
       )
        //TABENERGY( } ) /* endif (tabtype >= 0) */
 #if (TABENERGY (1+) 0)
@@ -299,7 +298,6 @@ MODIFIED(
       ALCHPAIR(
         ENERGY(vdwEnergy   += alch_vdw_energy;)
         FEP(vdwEnergy_s += alch_vdw_energy_2;)
-        FEP(vdwEnergy_s_Left += alch_vdw_energy_2_Left;)
         TI(ALCH1(vdwEnergy_ti_1) ALCH2(vdwEnergy_ti_2) += alch_vdw_dUdl;)
       ) // ALCHPAIR
       
@@ -350,17 +348,7 @@ MODIFIED(
       ) //ENERGY
       ALCHPAIR(
         ENERGY(electEnergy   -= myElecLambda * fast_val;)
-        if( Fep_Wham ) {
-               if( Fep_ElecOn ) {
-            FEP(electEnergy_s -= (myElecLambda * fast_val + fast_val);)
-               }
-               else {  // repulsive or dispersive, no constribution from charges
-            FEP(electEnergy_s -= (myElecLambda * fast_val);)
-               }
-        }
-        else{ // the default (orginal) FEP
-          FEP(electEnergy_s -= myElecLambda2 * fast_val;)
-        } 
+        FEP(electEnergy_s -= myElecLambda2 * fast_val;)
         TI(
           NOENERGY(register BigReal fast_val = 
             ( ( diffa * fast_d * (1/6.)+ fast_c * (1/4.)) * diffa + fast_b *(1/2.)) * diffa + fast_a;)
@@ -639,17 +627,7 @@ MODIFIED(
           
       ALCHPAIR(
         ENERGY(fullElectEnergy   -= myElecLambda * slow_val;)
-        if( Fep_Wham ) {
-               if(Fep_ElecOn)  {
-                       FEP(fullElectEnergy_s -= (myElecLambda * slow_val + slow_val);)
-               }
-               else    {
-            FEP(fullElectEnergy_s -= myElecLambda * slow_val;)
-               }
-        }
-        else{ // orignal FEP
-          FEP(fullElectEnergy_s -= myElecLambda2 * slow_val;)
-        }
+        FEP(fullElectEnergy_s -= myElecLambda2 * slow_val;)
         TI(
           NOENERGY(register BigReal slow_val =
                ( ( diffa * slow_d *(1/6.)+ slow_c * (1/4.)) * diffa + slow_b *(1/2.)) * diffa + slow_a;)
index 24bac70..e85c025 100644 (file)
 */
 
 #include "ComputeNonbondedInl.h"
-
-// 3 inline functions to handle explicit calculation of separation-shifted
-// vdW for FEP and TI, and a shifted electrostatics potential for decoupling
+#include "ComputeNonbondedAlch.h"
 
 /* ******************************************** */
 /* vdW energy, force and lambda2 energy for FEP */
 /* ******************************************** */
 inline void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2, 
-  BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, 
-  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, 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)                            \
-                            *(cutoff2 - 3.*switchdist2 + 2.*r2) : 1.);
-  const BigReal switchmul2 = (r2 > switchdist2 ?                       \
-                             12.*switchfactor*(cutoff2 - r2)           \
-                             *(r2 - switchdist2) : 0.);
-  
-  *alch_vdw_energy_2_Left = 0.;
-
-  if(Fep_WCA_repuOn){
-    if(B != 0.0) {// Excluded when some atoms like H, drude, lone pair are involved
-      const BigReal Emin = B*B/(4.0*A);
-      const BigReal Rmin_SQ = powf(2.0*A/B, 1.f/3);
-      const BigReal r2_1 = r2 + (1.-WCA_rcut1)*(1.-WCA_rcut1)*Rmin_SQ;
-      const BigReal r2_2 = r2 + (1.-WCA_rcut2)*(1.-WCA_rcut2)*Rmin_SQ;
-      const BigReal r2_3 = r2 + (1.-WCA_rcut3)*(1.-WCA_rcut3)*Rmin_SQ;
-      const BigReal r6_1 = r2_1*r2_1*r2_1;
-      const BigReal r6_2 = r2_2*r2_2*r2_2;
-      const BigReal r6_3 = r2_3*r2_3*r2_3;
-      const BigReal WCA_rcut1_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut1) \
-                                                      *(1.0-WCA_rcut1)) ? \
-                                       A/(r6_1*r6_1) - B/r6_1 + Emin : 0.);
-      const BigReal WCA_rcut2_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut2) \
-                                                      *(1.0-WCA_rcut2)) ? \
-                                       A/(r6_2*r6_2) - B/r6_2 + Emin : 0.);
-      const BigReal WCA_rcut3_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut3) \
-                                                      *(1.0-WCA_rcut3)) ? \
-                                       A/(r6_3*r6_3) - B/r6_3 + Emin : 0.);
-      const BigReal WCA_rcut1_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut1) \
-                                                     *(1.0-WCA_rcut1)) ? \
-                                      (12.*(WCA_rcut1_energy)          \
-                                       + 6.*B/r6_1 - 12.0 * Emin )/r2_1: 0.);
-      const BigReal WCA_rcut2_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut2) \
-                                                     *(1.0-WCA_rcut2))? \
-                                      (12.*(WCA_rcut2_energy)          \
-                                       + 6.*B/r6_2 - 12.0 * Emin )/r2_2: 0.);
-      const BigReal WCA_rcut3_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut3) \
-                                                     *(1.0-WCA_rcut3)) ? \
-                                      (12.*(WCA_rcut3_energy)          \
-                                       + 6.*B/r6_3 - 12.0 * Emin )/r2_3: 0.);
-      // separation-shifted repulsion force and energy
-      *alch_vdw_energy = WCA_rcut2_energy; 
-      *alch_vdw_force = WCA_rcut2_force;
-      if(WCA_rcut1 < WCA_rcut2) {
-       *alch_vdw_energy_2_Left = *alch_vdw_energy + WCA_rcut2_energy - WCA_rcut1_energy; 
+    BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, 
+    BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching,
+    BigReal myVdwLambda, BigReal myVdwLambda2, Bool alchWCAOn,
+    BigReal myRepLambda, BigReal myRepLambda2, BigReal* alch_vdw_energy,
+    BigReal* alch_vdw_force, BigReal* alch_vdw_energy_2) {
+  BigReal U, F, U_2, dU, dU_2, switchmul, switchmul2;
+  /*
+   * The variable naming here is unavoidably awful. A number immediately after
+   * a variable implies a distance to that power. The suffix "_2" indicates an
+   * evaluation at alchLambda2. The prefix "inv" means the inverse (previously
+   * this also used underscores - it was awful).
+   */
+  if (alchWCAOn) {
+    // WCA-on, auxilliary, lambda-dependent cutoff based on Rmin
+    //
+    // Avoid divide by zero - corectly zeroes interaction below.
+    const BigReal Rmin2 = (B <= 0.0 ? 0.0 : powf(2.0*A/B, 1.f/3));
+    if (myRepLambda < 1.0) {
+      // modified repulsive soft-core
+      // All of the scaling is baked into the shift and cutoff values. There
+      // are no linear coupling terms out front.
+      const BigReal WCAshift = Rmin2*(1 - myRepLambda)*(1 - myRepLambda);
+      if (r2 <= Rmin2 - WCAshift) {
+        const BigReal epsilon = B*B/(4.0*A);
+        vdw_forceandenergy(A, B, r2 + WCAshift, &U, &F);
+        *alch_vdw_energy = U + epsilon;
+        *alch_vdw_force = F;
+      } else {
+        *alch_vdw_energy = 0.0;
+        *alch_vdw_force = 0.0;
       }
-      if(WCA_rcut2 < WCA_rcut3) {
-       *alch_vdw_energy_2 = *alch_vdw_energy + WCA_rcut3_energy - WCA_rcut2_energy; 
+    } else { // myRepLambda == 1.0
+      if (vdwForceSwitching) {
+        // force and potential switching
+        if (r2 <= Rmin2) {
+          // normal LJ, potential w/two shifts
+          const BigReal epsilon = B*B/(4.0*A);
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon + myVdwLambda*dU;
+          *alch_vdw_force = F;
+        } else if (r2 <= switchdist2) {
+          // normal LJ potential w/shift
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*(U + dU);
+          *alch_vdw_force = myVdwLambda*F;
+        } else { // r2 > switchdist
+          // normal LJ potential with linear coupling   
+          vdw_fswitch_forceandenergy(A, B, r2, switchdist2, cutoff2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*U;
+          *alch_vdw_force = myVdwLambda*F;
+        }
+      } else {
+        // potential switching - also correct for no switching
+        if (r2 <= Rmin2) {
+          // normal LJ potential w/shift
+          const BigReal epsilon = B*B/(4.0*A);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon;
+          *alch_vdw_force = F; 
+        } else { // r2 > Rmin2
+          // normal LJ potential with linear coupling
+          vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
+              &switchmul2);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*switchmul*U;
+          *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
+        }
       }
     }
-    else {
-      *alch_vdw_energy = 0.0;
-      *alch_vdw_force = 0.0;
-      *alch_vdw_energy_2_Left = 0.0;
-      *alch_vdw_energy_2 = 0.0;
-    }
-  }
-  else if(Fep_WCA_dispOn) {
-    // separation-shifted dispersion force and energy
-    if(B == 0.0) {     // some atoms like H, drude, lone pair are involved
-      *alch_vdw_energy = 0.0;
-      *alch_vdw_force = 0.0;
-      *alch_vdw_energy_2 = 0.0;
-    }
-    else {
-      const BigReal Emin = B*B/(4.0*A);
-      const BigReal Rmin_SQ = powf(2.0*A/B, 1.f/3);
-      const BigReal r2_1 = r2;
-      const BigReal r2_2 = r2;
-      const BigReal r6_1 = r2_1*r2_1*r2_1;
-      const BigReal r6_2 = r2_2*r2_2*r2_2;
-      *alch_vdw_energy = r2 > Rmin_SQ? \
-                         myVdwLambda*(A/(r6_1*r6_1) - B/r6_1): \
-                         A/(r6_1*r6_1) - B/r6_1 + (1.-myVdwLambda)* Emin;
-      *alch_vdw_force =  r2 > Rmin_SQ? \
-                         (12.*(*alch_vdw_energy) + 6.*myVdwLambda*B/r6_1)/r2_1 * switchmul \
-                         + (*alch_vdw_energy) * switchmul2:(12.*(*alch_vdw_energy) \
-                         + 6.*B/r6_1 - 12.0*(1.-myVdwLambda)* Emin )/r2_1;
-      *alch_vdw_energy *= switchmul; 
-                       if(!Fep_Wham){ 
-        *alch_vdw_energy_2 = r2 > Rmin_SQ? \
-                             myVdwLambda2*switchmul*(A/(r6_1*r6_1) - B/r6_1): \
-                             A/(r6_1*r6_1) - B/r6_1 + (1.-myVdwLambda2)* Emin;
-                       }
-                       else{
-                               *alch_vdw_energy_2 = r2 > Rmin_SQ? \
-                                                          switchmul*(A/(r6_1*r6_1) - B/r6_1): - Emin;
-        *alch_vdw_energy_2 += *alch_vdw_energy;
-                       }
+    if (myRepLambda2 < 1.0) {
+      // Same as above, but energy only. This is done separately bc the
+      // cutoffs are lambda dependent.
+      const BigReal WCAshift_2 = Rmin2*(1 - myRepLambda2)*(1 - myRepLambda2);
+      if (r2 <= Rmin2 - WCAshift_2) {
+        const BigReal epsilon = B*B/(4.0*A);
+        vdw_energy(A, B, r2 + WCAshift_2, &U_2);
+        *alch_vdw_energy_2 = U_2 + epsilon;
+      } else {
+        *alch_vdw_energy_2 = 0.0;
+      }
+    } else { // myRepLambda2 == 1.0
+      if (vdwForceSwitching) {
+        // force and potential switching
+        if (r2 <= Rmin2) {
+          // normal LJ, potential w/two shifts
+          const BigReal epsilon = B*B/(4.0*A);
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU_2);
+          vdw_energy(A, B, r2, &U_2);
+          *alch_vdw_energy_2 = \
+              U_2 + (1 - myVdwLambda2)*epsilon + myVdwLambda2*dU_2;
+        } else if (r2 <= switchdist2) {
+          // normal LJ potential w/shift
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU_2);
+          vdw_energy(A, B, r2, &U_2);
+          *alch_vdw_energy_2 = myVdwLambda2*(U_2 + dU_2);
+        } else { // r2 > switchdist
+          vdw_fswitch_energy(A, B, r2, switchdist2, cutoff2, &U_2);
+          *alch_vdw_energy_2 = myVdwLambda2*U_2;
+        }
+      } else {
+        // potential switching - also correct for no switching
+        if (r2 <= Rmin2) {
+          // normal LJ potential w/shift
+          const BigReal epsilon = B*B/(4.0*A);
+          vdw_energy(A, B, r2, &U_2);
+          *alch_vdw_energy_2 = U_2 + (1 - myVdwLambda2)*epsilon;
+        } else { // r2 > Rmin2
+          // normal LJ potential with linear coupling
+          vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
+              &switchmul2);
+          vdw_energy(A, B, r2, &U_2);
+          *alch_vdw_energy_2 = myVdwLambda2*switchmul*U_2;
+        }
+      }
     }
-  } 
-  else {
-    //myVdwShift already multplied by relevant (1-vdwLambda)
-    const BigReal r2_1 = 1./(r2 + myVdwShift);
-    const BigReal r2_2 = 1./(r2 + myVdwShift2);
-    const BigReal r6_1 = r2_1*r2_1*r2_1;
-    const BigReal r6_2 = r2_2*r2_2*r2_2;
-    // separation-shifted vdW force and energy
-    const BigReal U1 = A*r6_1*r6_1 - B*r6_1; // NB: unscaled, shorthand only!
-    const BigReal U2 = A*r6_2*r6_2 - B*r6_2;
-    *alch_vdw_energy = myVdwLambda*switchmul*U1;
-    *alch_vdw_energy_2 = myVdwLambda2*switchmul*U2;
-    *alch_vdw_force = (myVdwLambda*(switchmul*(12.*U1 + 6.*B*r6_1)*r2_1 \
-                                    + switchmul2*U1));
-    // BKR - separation-shifted vdW force switching and potential shifting
-    if(vdwForceSwitching){ // add potential shifts term
-      const BigReal cutoff6 = cutoff2*cutoff2*cutoff2;
-      const BigReal cutoff3 = cutoff2*sqrt(cutoff2);
-
-      const BigReal shifted_switchdist2 = switchdist2 + myVdwShift;
-      const BigReal shifted_switchdist = sqrt(shifted_switchdist2);
-      const BigReal shifted_switchdist6 = \
-          shifted_switchdist2*shifted_switchdist2*shifted_switchdist2;
-      const BigReal shifted_switchdist3 = shifted_switchdist2*shifted_switchdist;
-     
-      const BigReal shifted_switchdist2_2 = switchdist2 + myVdwShift2;
-      const BigReal shifted_switchdist_2 = sqrt(shifted_switchdist2_2);
-      const BigReal shifted_switchdist6_2 = \
-          shifted_switchdist2_2*shifted_switchdist2_2*shifted_switchdist2_2;
-      const BigReal shifted_switchdist3_2 = shifted_switchdist2_2*shifted_switchdist_2;
-
-      if(r2 > switchdist2) {
-       const BigReal k_vdwa = A*cutoff6 / (cutoff6 - shifted_switchdist6);
-       const BigReal k_vdwb = B*cutoff3 / (cutoff3 - shifted_switchdist3);
-        const BigReal r_1 = sqrt(r2_1);
-        const BigReal tmpa = r6_1 - (1./cutoff6);
-        const BigReal tmpb = r2_1*r_1 - (1./cutoff3);
-
-        const BigReal k_vdwa2 = A*cutoff6 / (cutoff6 - shifted_switchdist6_2);
-        const BigReal k_vdwb2 = B*cutoff3 / (cutoff3 - shifted_switchdist3_2);
-        const BigReal r_2 = sqrt(r2_2);
-        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);
-        *alch_vdw_energy_2 = (myVdwLambda2*(k_vdwa2*tmpa2*tmpa2 \
-                                            - 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{
-        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;
+  } else { // WCA-off
+    if (vdwForceSwitching) {
+      // force and potential switching
+      if (r2 <= switchdist2) {
+        // normal LJ potential w/shift
+        vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
+        vdw_fswitch_shift(A, B, switchdist2 + myVdwShift, cutoff2, &dU);
+        vdw_energy(A, B, r2 + myVdwShift2, &U_2);
+        vdw_fswitch_shift(A, B, switchdist2 + myVdwShift2, cutoff2, &dU_2);
+        *alch_vdw_energy = myVdwLambda*(U + dU);
+        *alch_vdw_energy_2 = myVdwLambda2*(U_2 + dU_2);
+        *alch_vdw_force = myVdwLambda*F;
+      } else { // r2 > switchdist2
+        vdw_fswitch_forceandenergy(A, B, r2 + myVdwShift, \
+            switchdist2 + myVdwShift, cutoff2, &U, &F);
+        vdw_fswitch_energy(A, B, r2 + myVdwShift2, switchdist2 + myVdwShift2, \
+            cutoff2, &U_2);
+        *alch_vdw_energy = myVdwLambda*U;
+        *alch_vdw_energy_2 = myVdwLambda2*U_2;
+        *alch_vdw_force = myVdwLambda*F;
       }
+    } else {
+      // potential switching - also correct for no switching
+      vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
+          &switchmul2);
+      vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
+      vdw_energy(A, B, r2 + myVdwShift2, &U_2); 
+      *alch_vdw_energy = myVdwLambda*switchmul*U;
+      *alch_vdw_energy_2 = myVdwLambda2*switchmul*U_2;
+      *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
     }
   }
 }
index 71b972e..c37f6dd 100644 (file)
 */
 
 #include "ComputeNonbondedInl.h"
-
-// 3 inline functions to handle explicit calculation of separation-shifted
-// vdW for FEP and TI, and a shifted electrostatics potential for decoupling
+#include "ComputeNonbondedAlch.h"
 
 /* ********************************** */
 /* vdW energy, force and dU/dl for TI */
 /* ********************************** */
 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, 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;
-    
-  // switching function (this is correct whether switching is active or not)
-  const BigReal switchmul = (r2 > switchdist2 ? switchfactor*(cutoff2 - r2) \
-                            *(cutoff2 - r2) \
-                            *(cutoff2 - 3.*switchdist2 + 2.*r2) : 1.);
-  const BigReal switchmul2 = (r2 > switchdist2 ?                      \
-                              12.*switchfactor*(cutoff2 - r2)       \
-                              *(r2 - switchdist2) : 0.);
-
-  // separation-shifted vdW force and energy
-  const BigReal U = A*r6_1*r6_1 - B*r6_1; // NB: unscaled! for shorthand only!
-  *alch_vdw_energy = myVdwLambda*switchmul*U;
-  *alch_vdw_force = (myVdwLambda*(switchmul*(12.*U + 6.*B*r6_1)*r2_1 \
-                                  + switchmul2*U));
-  *alch_vdw_dUdl = (switchmul*(U + myVdwLambda*alchVdwShiftCoeff \
-                                   *(6.*U + 3.*B*r6_1)*r2_1));
-
-  // BKR - separation-shifted vdW force switching and potential shifting
-  if(vdwForceSwitching){ // add potential shifts and additional dU/dl terms
-    const BigReal cutoff6 = cutoff2*cutoff2*cutoff2;
-    const BigReal cutoff3 = cutoff2*sqrt(cutoff2);
-    const BigReal shifted_switchdist2 = switchdist2 + myVdwShift;
-    const BigReal shifted_switchdist = sqrt(shifted_switchdist2);
-    const BigReal shifted_switchdist6 = \
-        shifted_switchdist2*shifted_switchdist2*shifted_switchdist2;
-    const BigReal shifted_switchdist3 = shifted_switchdist2*shifted_switchdist;
-
-    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;
-      *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));
-      *alch_vdw_dUdl = (Uh + alchVdwShiftCoeff*myVdwLambda*r2_1 \
-                        *(6.*k_vdwa*tmpa*r6_1                   \
-                          - 3.*k_vdwb*tmpb*r2_1*r_1));
-    }else{
-      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;
+    BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2,
+    BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda,
+    BigReal alchVdwShiftCoeff, Bool alchWCAOn, BigReal myRepLambda,
+    BigReal* alch_vdw_energy, BigReal* alch_vdw_force,
+    BigReal* alch_vdw_dUdl) {
+  BigReal U, F, dU, switchmul, switchmul2;
+  /*
+   * The variable naming here is unavoidably awful. A number immediately after
+   * a variable implies a distance to that power. The suffix "_2" indicates an
+   * evaluation at alchLambda2. The prefix "inv" means the inverse (previously
+   * this also used underscores - it was awful).
+   */
+  if (alchWCAOn) {
+    // WCA-on, auxilliary, lambda-dependent cutoff based on Rmin
+    //
+    // Avoid divide by zero - correctly zeroes interaction below.
+    const BigReal Rmin2 = (B <= 0.0 ? 0.0 : powf(2.0*A/B, 1.f/3));
+    if (myRepLambda < 1.0) {
+      // modified repulsive soft-core
+      // All of the scaling is baked into the shift and cutoff values. There
+      // are no linear coupling terms out front.
+      const BigReal WCAshift = Rmin2*(1 - myRepLambda)*(1 - myRepLambda);
+      if (r2 <= Rmin2 - WCAshift) {
+        const BigReal epsilon = B*B/(4.0*A);
+        vdw_forceandenergy(A, B, r2 + WCAshift, &U, &F);
+        *alch_vdw_energy = U + epsilon;
+        *alch_vdw_force = F;
+        *alch_vdw_dUdl = Rmin2*(1 - myRepLambda)*F;
+      } else {
+        *alch_vdw_energy = 0.0;
+        *alch_vdw_force = 0.0;
+        *alch_vdw_dUdl = 0.0;
+      }
+    } else { // myRepLambda == 1.0
+      if (vdwForceSwitching) {
+        // force and potential switching
+        if (r2 <= Rmin2) {
+          // normal LJ, potential w/two shifts 
+          const BigReal epsilon = B*B/(4.0*A); 
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon + myVdwLambda*dU;
+          *alch_vdw_force = F;
+          *alch_vdw_dUdl = dU - epsilon;
+        } else if (r2 <= switchdist2) {
+          // normal LJ potential w/shift
+          vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*(U + dU);
+          *alch_vdw_force = myVdwLambda*F;
+          *alch_vdw_dUdl = U + dU;
+        } else { // r2 > switchdist
+          // normal LJ potential with linear coupling
+          vdw_fswitch_forceandenergy(A, B, r2, switchdist2, cutoff2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*U;
+          *alch_vdw_force = myVdwLambda*F;
+          *alch_vdw_dUdl = U;
+        }
+      } else {
+        // potential switching - also correct for no switching
+        if (r2 <= Rmin2) {
+          // normal LJ potential w/shift
+          const BigReal epsilon = B*B/(4.0*A);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon;
+          *alch_vdw_force = F;
+          *alch_vdw_dUdl = -epsilon;
+        } else { // r2 > Rmin2
+          // normal LJ potential with linear coupling
+          vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
+              &switchmul2);
+          vdw_forceandenergy(A, B, r2, &U, &F);
+          *alch_vdw_energy = myVdwLambda*switchmul*U;
+          *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
+          *alch_vdw_dUdl = switchmul*U;
+        }
+      }
+    }
+  } else { // WCA-off
+    if (vdwForceSwitching) {
+      // force and potential switching
+      if (r2 <= switchdist2) {
+        // normal LJ potential w/shift
+        vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
+        vdw_fswitch_shift(A, B, switchdist2 + myVdwShift, cutoff2, &dU);
+        *alch_vdw_energy = myVdwLambda*(U + dU);
+        *alch_vdw_force = myVdwLambda*F;
+        *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F + dU;
+      } else { // r2 > switchdist2
+        vdw_fswitch_forceandenergy(A, B, r2 + myVdwShift, \
+            switchdist2 + myVdwShift, cutoff2, &U, &F);
+        *alch_vdw_energy = myVdwLambda*U;
+        *alch_vdw_force = myVdwLambda*F;
+        *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F;
+      }
+    } else {
+      // potential switching - also correct for no switching
+      vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
+          &switchmul2);
+      vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
+      *alch_vdw_energy = myVdwLambda*switchmul*U;
+      *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
+      *alch_vdw_dUdl = switchmul*(U + 0.5*myVdwLambda*alchVdwShiftCoeff*F);
     }
   }
 }
index 3356c8b..bba7b5e 100644 (file)
@@ -113,16 +113,7 @@ BigReal         ComputeNonbondedUtil::c8;
 // fepb
 Bool      ComputeNonbondedUtil::alchFepOn;
 Bool      ComputeNonbondedUtil::alchThermIntOn;
-Bool      ComputeNonbondedUtil::Fep_WCA_repuOn;
-Bool      ComputeNonbondedUtil::Fep_WCA_dispOn;
-Bool      ComputeNonbondedUtil::Fep_ElecOn;
-Bool      ComputeNonbondedUtil::Fep_Wham; 
-BigReal   ComputeNonbondedUtil::WCA_rcut1;
-BigReal   ComputeNonbondedUtil::WCA_rcut2;
-BigReal   ComputeNonbondedUtil::WCA_rcut3;
-BigReal   ComputeNonbondedUtil::alchRepLambda;
-BigReal   ComputeNonbondedUtil::alchDispLambda;
-BigReal   ComputeNonbondedUtil::alchElecLambda;
+Bool      ComputeNonbondedUtil::alchWCAOn;
 BigReal   ComputeNonbondedUtil::alchVdwShiftCoeff;
 Bool      ComputeNonbondedUtil::vdwForceSwitching;
 Bool      ComputeNonbondedUtil::alchDecouple;
@@ -202,7 +193,6 @@ void ComputeNonbondedUtil::submitReductionData(BigReal *data, SubmitReduction *r
   reduction->item(REDUCTION_ELECT_ENERGY_F) += data[electEnergyIndex_s];
   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += data[fullElectEnergyIndex_s];
   reduction->item(REDUCTION_LJ_ENERGY_F) += data[vdwEnergyIndex_s];
-  reduction->item(REDUCTION_LJ_ENERGY_F_LEFT) += data[vdwEnergyIndex_s_Left];
 
   reduction->item(REDUCTION_ELECT_ENERGY_TI_1) += data[electEnergyIndex_ti_1];
   reduction->item(REDUCTION_ELECT_ENERGY_SLOW_TI_1) += data[fullElectEnergyIndex_ti_1];
@@ -293,25 +283,16 @@ void ComputeNonbondedUtil::select(void)
 
 //fepb
   alchFepOn = simParams->alchFepOn;
-  Fep_WCA_repuOn = simParams->alchFepWCARepuOn;
-  Fep_WCA_dispOn = simParams->alchFepWCADispOn;
-  Fep_ElecOn = simParams->alchFepElecOn;
-  Fep_Wham = simParams->alchFepWhamOn;
   alchThermIntOn = simParams->alchThermIntOn;
+  alchWCAOn = simParams->alchWCAOn;
+  alchDecouple = simParams->alchDecouple;
+
   lesOn = simParams->lesOn;
   lesScaling = lesFactor = 0;
+
   Bool tabulatedEnergies = simParams->tabulatedEnergies;
   alchVdwShiftCoeff = simParams->alchVdwShiftCoeff;
   vdwForceSwitching = simParams->vdwForceSwitching;
-  WCA_rcut1 = simParams->alchFepWCArcut1;
-  WCA_rcut2 = simParams->alchFepWCArcut2;
-  WCA_rcut3 = simParams->alchFepWCArcut3;
-
-  alchRepLambda = simParams->alchRepLambda;
-  alchDispLambda = simParams->alchDispLambda;
-  alchElecLambda = simParams->alchElecLambda;
-
-  alchDecouple = simParams->alchDecouple;
 
   delete [] lambda_table;
   lambda_table = 0;
index 91a0238..f118300 100644 (file)
@@ -265,7 +265,7 @@ public:
         electEnergyIndex, fullElectEnergyIndex, vdwEnergyIndex,
         goNativeEnergyIndex, goNonnativeEnergyIndex, groLJEnergyIndex, groGaussEnergyIndex,
 //sd-db
-        electEnergyIndex_s, fullElectEnergyIndex_s, vdwEnergyIndex_s, vdwEnergyIndex_s_Left,
+        electEnergyIndex_s, fullElectEnergyIndex_s, vdwEnergyIndex_s,
         electEnergyIndex_ti_1, fullElectEnergyIndex_ti_1, vdwEnergyIndex_ti_1,
         electEnergyIndex_ti_2, fullElectEnergyIndex_ti_2, vdwEnergyIndex_ti_2,
 //sd-de
@@ -355,19 +355,10 @@ public:
 //sd-db
   static Bool alchFepOn;
   static Bool alchThermIntOn;
+  static Bool alchWCAOn;
   static BigReal alchVdwShiftCoeff;
   static Bool vdwForceSwitching;
-  static Bool Fep_WCA_repuOn;
-  static Bool Fep_WCA_dispOn;
-  static Bool Fep_ElecOn;
-  static Bool Fep_Wham;
-  static BigReal WCA_rcut1;
-  static BigReal WCA_rcut2; 
-  static BigReal WCA_rcut3; 
   static Bool alchDecouple;
-  static BigReal alchRepLambda;
-  static BigReal alchDispLambda;
-  static BigReal alchElecLambda;
 //sd-de
   static Bool lesOn;
   static int lesFactor;
index 8e30dea..5e8b37f 100644 (file)
@@ -4057,24 +4057,12 @@ void ComputePme::ungridForces() {
       if (alchOn) {
         float scale = 1.;
         BigReal elecLambdaUp, elecLambdaDown;
-        if ( simParams->alchFepWhamOn ) {
-         if ( simParams->alchFepElecOn ) {
-            elecLambdaUp = simParams->alchElecLambda;
-            elecLambdaDown = 1.0 - simParams->alchElecLambda;
-         }
-         else {
-            elecLambdaUp = 0.0;
-            elecLambdaDown = 1.0;
-         }
-        }
-        else {
-          BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
-          myMgr->alchLambda = alchLambda;
-          BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
-          myMgr->alchLambda2 = alchLambda2;
-         elecLambdaUp = simParams->getElecLambda(alchLambda);
-         elecLambdaDown = simParams->getElecLambda(1. - alchLambda);
-        }
+        BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
+        myMgr->alchLambda = alchLambda;
+        BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
+        myMgr->alchLambda2 = alchLambda2;
+        elecLambdaUp = simParams->getElecLambda(alchLambda);
+        elecLambdaDown = simParams->getElecLambda(1. - alchLambda);
        
         if ( g == 0 ) scale = elecLambdaUp;
         else if ( g == 1 ) scale = elecLambdaDown;
@@ -4117,25 +4105,13 @@ void ComputePme::ungridForces() {
         }
       } else if ( lesOn ) {
         float scale = 1.;
-        if ( alchFepOn ) {
-         if(simParams->alchFepWhamOn) {
-           if(simParams->alchFepElecOn) {
-             if ( g == 0 ) scale = simParams->alchElecLambda;
-             else if ( g == 1 ) scale = 1. - simParams->alchElecLambda;
-           }
-           else {
-             if ( g == 0 ) scale = 0.0;
-             else if ( g == 1 ) scale = 1.0;
-           }
-         }
-         else {
-            BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
-            myMgr->alchLambda = alchLambda;
-            BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
-            myMgr->alchLambda2 = alchLambda2;
-            if ( g == 0 ) scale = alchLambda;
-            else if ( g == 1 ) scale = 1. - alchLambda;
-         }
+        if ( alchFepOn ) { 
+          BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
+          myMgr->alchLambda = alchLambda;
+          BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
+          myMgr->alchLambda2 = alchLambda2;
+          if ( g == 0 ) scale = alchLambda;
+          else if ( g == 1 ) scale = 1. - alchLambda;
         } else if ( lesOn ) {
           scale = 1.0 / (float)lesFactor;
         }
@@ -4218,24 +4194,12 @@ void ComputePmeMgr::submitReductions() {
       float scale = 1.;
       if (alchOn) {
         BigReal elecLambdaUp, elecLambdaDown;
-        if( simParams->alchFepWhamOn ) {
-          if( simParams->alchFepElecOn ) {
-            elecLambdaUp = simParams->alchElecLambda;
-            elecLambdaDown = 1.0 - simParams->alchElecLambda;
-          }
-          else {
-            elecLambdaUp = 0.0;
-            elecLambdaDown = 1.0;
-          }
-        }
-        else {
-          // alchLambda set on each step in ComputePme::ungridForces()
-          if ( alchLambda < 0 || alchLambda > 1 ) {
-            NAMD_bug("ComputePmeMgr::submitReductions alchLambda out of range");
-          }
-          elecLambdaUp = simParams->getElecLambda(alchLambda);
-          elecLambdaDown = simParams->getElecLambda(1-alchLambda);
+        // alchLambda set on each step in ComputePme::ungridForces()
+        if ( alchLambda < 0 || alchLambda > 1 ) {
+          NAMD_bug("ComputePmeMgr::submitReductions alchLambda out of range");
         }
+        elecLambdaUp = simParams->getElecLambda(alchLambda);
+        elecLambdaDown = simParams->getElecLambda(1-alchLambda);
         if ( g == 0 ) scale = elecLambdaUp;
         else if ( g == 1 ) scale = elecLambdaDown;
         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
@@ -4267,21 +4231,8 @@ void ComputePmeMgr::submitReductions() {
 
       if (alchFepOn) {
        BigReal elecLambda2Up=0.0, elecLambda2Down=0.0;
-        if(simParams->alchFepWhamOn) {
-          if(simParams->alchFepElecOn) {
-            elecLambda2Up = simParams->alchElecLambda;
-            elecLambda2Down =  1.0 - simParams->alchElecLambda;
-          }
-          else {
-            elecLambda2Up = 0.0;
-            elecLambda2Down =  1.0;
-          }
-        }
-        else {
-          elecLambda2Up = simParams->getElecLambda(alchLambda2);
-          elecLambda2Down = simParams->getElecLambda(1.-alchLambda2);
-        }
-        
+        elecLambda2Up = simParams->getElecLambda(alchLambda2);
+        elecLambda2Down = simParams->getElecLambda(1.-alchLambda2);        
         if ( g == 0 ) scale2 = elecLambda2Up;
         else if ( g == 1 ) scale2 = elecLambda2Down;
         else if ( g == 2 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
@@ -4289,12 +4240,6 @@ void ComputePmeMgr::submitReductions() {
         else if (alchDecouple && g == 3 ) scale2 = 1 - elecLambda2Down;
         else if (alchDecouple && g == 4 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
       }
-      if(simParams->alchFepWhamOn && simParams->alchFepElecOn) {       // FEP with wham post-process
-       if( g==0 )      scale2 = scale + 1.0;
-       else if( g==1 ) scale2 = scale - 1.0;
-       else if( g==2 ) scale2 = scale - 1.0;
-       else if( g==3 ) scale2 = scale + 1.0;
-      }
       reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += evir[g][0] * scale2;
       
       if (alchThermIntOn) {
index f2b4726..080bc75 100644 (file)
@@ -3022,7 +3022,6 @@ void Controller::printEnergies(int step, int minimize)
       bondedEnergyDiff_f = reduction->item(REDUCTION_BONDED_ENERGY_F);
       electEnergy_f = reduction->item(REDUCTION_ELECT_ENERGY_F);
       ljEnergy_f = reduction->item(REDUCTION_LJ_ENERGY_F);
-      ljEnergy_f_left = reduction->item(REDUCTION_LJ_ENERGY_F_LEFT);
 
       bondedEnergy_ti_1 = reduction->item(REDUCTION_BONDED_ENERGY_TI_1);
       electEnergy_ti_1 = reduction->item(REDUCTION_ELECT_ENERGY_TI_1);
@@ -3613,53 +3612,50 @@ void Controller::enqueueCollections(int timestep)
 
 //Modifications for alchemical fep
 void Controller::outputFepEnergy(int step) {
- if (simParams->alchOn && simParams->alchFepOn) {
-  const int stepInRun = step - simParams->firstTimestep;
-  const int alchEquilSteps = simParams->alchEquilSteps;
-  const BigReal alchLambda = simParams->alchLambda;
-  const bool alchEnsembleAvg = simParams->alchEnsembleAvg; 
-  const bool FepWhamOn = simParams->alchFepWhamOn;
-
-  if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
-    FepNo = 0;
-    exp_dE_ByRT = 0.0;
-    net_dE = 0.0;
-  }
-  dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
-       electEnergy - electEnergySlow - ljEnergy;
-  BigReal RT = BOLTZMANN * simParams->alchTemp;
+  if (simParams->alchOn && simParams->alchFepOn) {
+    const int stepInRun = step - simParams->firstTimestep;
+    const int alchEquilSteps = simParams->alchEquilSteps;
+    const BigReal alchLambda = simParams->alchLambda;
+    const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
 
-  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. || (step / simParams->alchIDWSFreq) % 2 == 1 )){
-    // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
-    FepNo++;
-    exp_dE_ByRT += exp(-dE/RT);
-    net_dE += dE;
-  }
+    if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
+      FepNo = 0;
+      exp_dE_ByRT = 0.0;
+      net_dE = 0.0;
+    }
+    dE = bondedEnergyDiff_f + electEnergy_f + electEnergySlow_f + ljEnergy_f -
+         electEnergy - electEnergySlow - ljEnergy;
+    BigReal RT = BOLTZMANN * simParams->alchTemp;
+
+    if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
+        (step / simParams->alchIDWSFreq) % 2 == 1 )) {
+      // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
+      FepNo++;
+      exp_dE_ByRT += exp(-dE/RT);
+      net_dE += dE;
+    }
 
-  if (simParams->alchOutFreq) { 
-    if (stepInRun == 0) {
-      if (!fepFile.is_open()) {
-        NAMD_backup_file(simParams->alchOutFile);
-        fepFile.open(simParams->alchOutFile);
-        iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
-        if(alchEnsembleAvg){
-          fepSum = 0.0;
-          fepFile << "#            STEP                 Elec                            "
-                  << "vdW                    dE           dE_avg         Temp             dG\n"
-                  << "#                           l             l+dl      "
-                  << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
-        }
-        else{
-          if(!FepWhamOn){ 
+    if (simParams->alchOutFreq) {
+      if (stepInRun == 0) {
+        if (!fepFile.is_open()) {
+          NAMD_backup_file(simParams->alchOutFile);
+          fepFile.open(simParams->alchOutFile);
+          iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
+          if(alchEnsembleAvg){
+            fepSum = 0.0;
+            fepFile << "#            STEP                 Elec                            "
+                    << "vdW                    dE           dE_avg         Temp             dG\n"
+                    << "#                           l             l+dl      "
+                    << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
+          }
+          else{
             fepFile << "#            STEP                 Elec                            "
                     << "vdW                    dE         Temp\n"
                     << "#                           l             l+dl      "
                     << "       l            l+dl         E(l+dl)-E(l)" << std::endl;
-          } 
+          }
         }
-      }
-      if(!step){
-        if(!FepWhamOn){
+        if(!step){
           fepFile << "#NEW FEP WINDOW: "
                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
                   << simParams->alchLambda2;
@@ -3669,25 +3665,24 @@ void Controller::outputFepEnergy(int step) {
           fepFile << std::endl;
         }
       }
-    }
-    if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
-      fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
-              << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
-              << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
-    }
-    if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
-      writeFepEnergyData(step, fepFile);
-      fepFile.flush();
-    }
-    if (alchEnsembleAvg && (step == simParams->N)) {
-      fepSum = fepSum + dG;
-      fepFile << "#Free energy change for lambda window [ " << alchLambda
-              << " " << simParams->alchLambda2 << " ] is " << dG 
-              << " ; net change until now is " << fepSum << std::endl;
-      fepFile.flush();
+      if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
+        fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
+                << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
+                << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
+      }
+      if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
+        writeFepEnergyData(step, fepFile);
+        fepFile.flush();
+      }
+      if (alchEnsembleAvg && (step == simParams->N)) {
+        fepSum = fepSum + dG;
+        fepFile << "#Free energy change for lambda window [ " << alchLambda
+                << " " << simParams->alchLambda2 << " ] is " << dG
+                << " ; net change until now is " << fepSum << std::endl;
+        fepFile.flush();
+      }
     }
   }
- }
 }
 
 void Controller::outputTiEnergy(int step) {
@@ -3851,81 +3846,36 @@ BigReal Controller::computeAlchWork(const int step) {
 void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
   BigReal eeng = electEnergy + electEnergySlow;
   BigReal eeng_f = electEnergy_f + electEnergySlow_f;
-  BigReal dE_Left = eeng_f + ljEnergy_f_left - eeng - ljEnergy;
   BigReal RT = BOLTZMANN * simParams->alchTemp;
 
   const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
   const int stepInRun = step - simParams->firstTimestep;
-  const bool FepWhamOn = simParams->alchFepWhamOn;
-  const bool WCARepuOn = simParams->alchFepWCARepuOn;
-  const BigReal WCArcut1 = simParams->alchFepWCArcut1;
-  const BigReal WCArcut2 = simParams->alchFepWCArcut2;
-  const BigReal WCArcut3 = simParams->alchFepWCArcut3;
   const BigReal alchLambda = simParams->alchLambda;
-       
-  const BigReal alchRepLambda = simParams->alchRepLambda;
-  const BigReal alchDispLambda = simParams->alchDispLambda;
-  const BigReal alchElecLambda = simParams->alchElecLambda;
 
   if(stepInRun){
-    if(!FepWhamOn){
-      if ( simParams->alchLambdaIDWS >= 0. && (step / simParams->alchIDWSFreq) % 2 == 0 ) {
-        // IDWS is active and we are on a "backward" timestep
-        fepFile << FEPTITLE_BACK(step);
-      } else {
-        // "forward" timestep
-        fepFile << FEPTITLE(step);
-      }
-      fepFile << FORMAT(eeng);
-      fepFile << FORMAT(eeng_f);
-      fepFile << FORMAT(ljEnergy);
-      fepFile << FORMAT(ljEnergy_f);
-    }
-    else{ // FepWhamOn = ON
-      if(WCARepuOn){
-        if(WCArcut1<WCArcut2) { // [s1,s2]
-          fepFile << "FEP_WCA_REP  ";
-          fepFile << FORMAT(WCArcut1);
-          fepFile << FORMAT(WCArcut2);
-          fepFile << FORMAT(1.0);
-          fepFile << FORMAT(dE_Left);
-        }
-        if(WCArcut2<WCArcut3) { // [s2,s3]
-          if(WCArcut1<WCArcut2) fepFile << " BREAK ";
-          fepFile << "FEP_WCA_REP  ";
-          fepFile << FORMAT(WCArcut2);
-          fepFile << FORMAT(WCArcut3);
-          fepFile << FORMAT(0.0);
-          fepFile << FORMAT(dE);
-        }
-        fepFile << std::endl;
-      }
-      else if(simParams->alchFepWCADispOn) {
-        fepFile << "FEP_WCA_DISP ";
-        fepFile << FORMAT(alchDispLambda);
-      }
-      else if(simParams->alchFepElecOn)        {
-        fepFile << "FEP_ELEC     ";
-        fepFile << FORMAT(alchElecLambda);
-      }
-    }
-    if( ! WCARepuOn ) {
-      fepFile << FORMAT(dE);
+    if ( simParams->alchLambdaIDWS >= 0. &&
+        (step / simParams->alchIDWSFreq) % 2 == 0 ) {
+      // IDWS is active and we are on a "backward" timestep
+      fepFile << FEPTITLE_BACK(step);
+    } else {
+      // "forward" timestep
+      fepFile << FEPTITLE(step);
     }
+    fepFile << FORMAT(eeng);
+    fepFile << FORMAT(eeng_f);
+    fepFile << FORMAT(ljEnergy);
+    fepFile << FORMAT(ljEnergy_f);
+    fepFile << FORMAT(dE);
     if(alchEnsembleAvg){
-      BigReal dE_avg = net_dE/FepNo;
+      BigReal dE_avg = net_dE / FepNo;
       fepFile << FORMAT(dE_avg);
     }
-    if(!FepWhamOn){
-      fepFile << FORMAT(temperature);
-    }
+    fepFile << FORMAT(temperature);
     if(alchEnsembleAvg){
-      dG = -(RT * log(exp_dE_ByRT/FepNo));
+      dG = -(RT * log(exp_dE_ByRT / FepNo));
       fepFile << FORMAT(dG);
-    } 
-    if( ! WCARepuOn ) {
-      fepFile << std::endl;
     }
+    fepFile << std::endl;
   }
 }
 
index 1b3360f..126dcae 100644 (file)
@@ -296,52 +296,6 @@ void SimParameters::scriptSet(const char *param, const char *value) {
   }
   SCRIPT_PARSE_INT("alchEquilSteps",alchEquilSteps)
 
-  if ( ! strncasecmp(param,"alchRepLambda",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchRepLambda = atof(value);
-    alchFepWCARepuOn = true;
-    alchFepWCADispOn = false;
-    alchFepElecOn    = false;
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
-  if ( ! strncasecmp(param,"alchDispLambda",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchDispLambda = atof(value);
-    alchFepWCARepuOn = false;
-    alchFepWCADispOn = true;
-    alchFepElecOn    = false;
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
-  if ( ! strncasecmp(param,"alchElecLambda",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchElecLambda = atof(value);
-    alchFepWCARepuOn = false;
-    alchFepWCADispOn = false;
-    alchFepElecOn    = true;
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
-
-  if ( ! strncasecmp(param,"alchFepWCArcut1",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchFepWCArcut1 = atof(value);
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
-  if ( ! strncasecmp(param,"alchFepWCArcut2",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchFepWCArcut2 = atof(value);
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
-  if ( ! strncasecmp(param,"alchFepWCArcut3",MAX_SCRIPT_PARAM_SIZE) ) {
-    alchFepWCArcut3 = atof(value);
-    ComputeNonbondedUtil::select();
-    return;
-  }
-
   if ( ! strncasecmp(param,"alchLambda",MAX_SCRIPT_PARAM_SIZE) ) {
     alchLambda = atof(value);
     if ( alchLambda < 0.0 || 1.0 < alchLambda ) {
@@ -1111,15 +1065,23 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
      "the altered alchemical vDW interactions", &alchVdwShiftCoeff, 5.);
    opts.range("alchVdwShiftCoeff", NOT_NEGATIVE);
 
+   opts.optionalB("alch", "alchWCA", "Is WCA decomposition being performed?",
+     &alchWCAOn, FALSE);
+
    // scheduling options for different interaction types
    opts.optional("alch", "alchElecLambdaStart", "Lambda at which electrostatic"
       "scaling of exnihilated particles begins", &alchElecLambdaStart, 0.5);
    opts.range("alchElecLambdaStart", NOT_NEGATIVE);
 
    opts.optional("alch", "alchVdwLambdaEnd", "Lambda at which vdW"
-      "scaling of exnihilated particles begins", &alchVdwLambdaEnd, 1.0);
+      "scaling of exnihilated particles ends", &alchVdwLambdaEnd, 1.0);
    opts.range("alchVdwLambdaEnd", NOT_NEGATIVE);
 
+   opts.optional("alch", "alchRepLambdaEnd", "Lambda at which repulsive vdW"
+      "scaling of exnihilated particles ends and attractive vdW scaling"
+      "begins", &alchRepLambdaEnd, 0.5);
+   opts.range("alchRepLambdaEnd", NOT_NEGATIVE);
+
    opts.optional("alch", "alchBondLambdaEnd", "Lambda at which bonded"
       "scaling of exnihilated particles begins", &alchBondLambdaEnd, 0.0);
    opts.range("alchBondLambdaEnd", NOT_NEGATIVE);
@@ -1145,35 +1107,8 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
      "data collection in the alchemical window", &alchEquilSteps, 0);
    opts.range("alchEquilSteps", NOT_NEGATIVE);
 
-   // WCA decomposition options
-   opts.optionalB("alch", "alchFepWCARepuOn",
-     "WCA decomposition repu interaction in use?", &alchFepWCARepuOn, FALSE);
-   opts.optionalB("alch", "alchFepWCADispOn",
-     "WCA decomposition disp interaction in use?", &alchFepWCADispOn, FALSE);
    opts.optionalB("alch", "alchEnsembleAvg", "Ensemble Average in use?",
      &alchEnsembleAvg, TRUE);
-   opts.optionalB("alch", "alchFepWhamOn",
-     "Energy output for Wham postprocessing in use?", &alchFepWhamOn, FALSE);
-   opts.optional("alch", "alchFepWCArcut1",
-     "WCA repulsion Coeff1 used for generating the altered alchemical vDW "
-     "interactions", &alchFepWCArcut1, 0.0);
-   opts.range("alchFepWCArcut1", NOT_NEGATIVE);
-   opts.optional("alch", "alchFepWCArcut2", "WCA repulsion Coeff2 used for "
-     "generating the altered alchemical vDW interactions", &alchFepWCArcut2,
-     1.0);
-   opts.range("alchFepWCArcut2", NOT_NEGATIVE);
-   opts.optional("alch", "alchFepWCArcut3",
-     "WCA repulsion Coeff3 used for generating the altered alchemical vDW "
-     "interactions", &alchFepWCArcut3, 1.0);
-   opts.range("alchFepWCArcut3", NOT_NEGATIVE);
-   // These default to invalid lambda values.
-   opts.optional("alch", "alchRepLambda", "Lambda of WCA repulsion"
-     "Coupling parameter value for WCA repulsion", &alchRepLambda, -1.0);
-   opts.optional("alch", "alchDispLambda", "Lambda of WCA dispersion"
-     "Coupling parameter value for WCA dispersion", &alchDispLambda, -1.0);
-   opts.optional("alch", "alchElecLambda", "Lambda of electrostatic "
-     "perturbation Coupling parameter value for electrostatic perturbation",
-     &alchElecLambda, -1.0);
 //fepe
 
    opts.optionalB("main", "les", "Is locally enhanced sampling enabled?",
@@ -3519,10 +3454,6 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
    alchOnAtStartup = alchOn;
 
    if (alchOn) {
-     if (vdwForceSwitching && (alchFepWCARepuOn || alchFepWCADispOn)) {
-       iout << iWARN << "vdwForceSwitching not implemented for alchemical "
-         "interactions when WCA decomposition is on!\n" << endi;
-     }
      if (martiniSwitching) {
        iout << iWARN << "Martini switching disabled for alchemical "
          "interactions.\n" << endi;
@@ -3565,6 +3496,18 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
 
      if (alchElecLambdaStart > 1.0) 
        NAMD_die("alchElecLambdaStart should be in the range [0.0, 1.0]\n");
+
+     if (alchWCAOn) {
+       if (alchRepLambdaEnd > 1.0)
+         NAMD_die("alchRepLambdaEnd should be in the range [0.0, 1.0]\n");
+       if (alchVdwLambdaEnd < alchRepLambdaEnd)
+         NAMD_die("alchVdwLambdaEnd should be greater than alchRepLambdaEnd\n");
+       if (alchVdwShiftCoeff > 0.0) {
+         iout << iWARN << "alchVdwShiftCoeff is non-zero but not used when WCA"
+              << " is active. Setting it to zero now.\n" << endi;
+         alchVdwShiftCoeff = 0.0;
+       }
+     }
      
      if (alchFepOn) {
        if (alchLambda2 < 0.0 || alchLambda2 > 1.0)
@@ -3577,57 +3520,9 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
          strcat(alchOutFile, ".fep");
        }
 
-       if (!alchFepWhamOn && 
-           ((!opts.defined("alchLambda")) || (!opts.defined("alchLambda2")))) {
+       if (!opts.defined("alchLambda") || !opts.defined("alchLambda2")) {
          NAMD_die("alchFepOn is on, but alchLambda or alchLambda2 is not set.");
        }
-       
-       if(alchRepLambda > 1.0)
-         NAMD_die("alchRepLambda should be in the range [0.0, 1.0].");
-       else if(alchRepLambda >= 0.0)
-         alchFepWCARepuOn = true;
-       else
-         alchFepWCARepuOn = false;
-       if(alchDispLambda > 1.0)
-         NAMD_die("alchDispLambda should be in the range [0.0, 1.0].");
-       else if(alchDispLambda >= 0.0)
-         alchFepWCADispOn = true;
-       else
-         alchFepWCADispOn = false;
-       if(alchElecLambda > 1.0)
-         NAMD_die("alchElecLambda should be in the range [0.0, 1.0].");
-       else if(alchElecLambda >= 0.0)
-         alchFepElecOn = true;
-       else
-         alchFepElecOn = false;
-       
-       if ((alchFepWCARepuOn || alchFepWCADispOn || alchFepElecOn) && 
-           !alchFepWhamOn)
-         NAMD_die("alchFepWhamOn has to be on if one of alchFepWCARepuOn/alchFepWCADispOn/alchFepElecOn is set.");
-       if (alchFepWCARepuOn && alchFepWCADispOn)
-          NAMD_die("With WCA decomposition, repulsion and dispersion can NOT be in the same FEP stage");
-       if (alchFepWCARepuOn && alchFepElecOn)
-          NAMD_die("With WCA decomposition, repulsion and electrostatic perturbation can NOT be in the same FEP stage");
-       if (alchFepWCADispOn && alchFepElecOn)
-          NAMD_die("With WCA decomposition, dispersion and electrostatic perturbation can NOT be in the same FEP stage");
-       if (alchFepWCARepuOn && 
-           (!opts.defined("alchFepWCArcut1") || 
-            !opts.defined("alchFepWCArcut2") ||
-            !opts.defined("alchFepWCArcut3") ))
-          NAMD_die("When using WCA repulsion,  alchFepWCArcut1, alchFepWCArcut2, and alchFepWCArcut3 must be defined!");
-       if (alchFepWCARepuOn && 
-           ((alchFepWCArcut1 > alchFepWCArcut2) || 
-            (alchFepWCArcut2 > alchFepWCArcut3) ))
-           NAMD_die("When using WCA repulsion,  alchFepWCArcut2 must be larger than alchFEPWCArcut1, alchFepWCArcut3 must be larger than alchFEPWCArcut2!");
-       if (alchFepWhamOn && (alchRepLambda < 0.0) && (alchDispLambda < 0.0) && 
-           (alchElecLambda < 0.0) )    
-                  NAMD_die("One of alchRepLambda, alchDispLambda and alchElecLambda should be set up when alchFepWhamOn is true!");
-       if (alchFepWhamOn && (!alchFepElecOn)) {
-                alchElecLambda = 0.0;
-                ComputeNonbondedUtil::alchElecLambda = alchElecLambda;
-       }
      }
      else if (alchThermIntOn) {
        // alchLambda2 is only needed for nonequilibrium switching
@@ -3649,6 +3544,12 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
      NAMD_die("Sorry, combined LES with FEP or TI is not implemented.\n");
    if ( alchOn && alchThermIntOn && lesOn )
      NAMD_die("Sorry, combined LES and TI is not implemented.\n");
+   if ( alchWCAOn && !alchOn ) {
+     iout << iWARN << "Alchemical WCA decomposition was requested but \
+       alchemical free energy calculation is not active. Setting \
+       alchWCA to off.\n" << endi;
+     alchWCAOn = FALSE;
+   }
    if ( alchDecouple && !alchOn ) {
          iout << iWARN << "Alchemical decoupling was requested but \
            alchemical free energy calculation is not active. Setting \
@@ -5122,38 +5023,45 @@ if ( openatomOn )
        iout << iINFO << "FEP INTRA-ALCHEMICAL BONDED INTERACTIONS WILL BE "
             << "RETAINED\n";
      }
-     iout << iINFO << "FEP VDW SHIFTING COEFFICIENT "
-          << alchVdwShiftCoeff << "\n";
+     if (alchWCAOn) {
+       iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN (WCA) VDW DECOUPLING "
+            << "ACTIVE\n";
+     } else {
+       iout << iINFO << "FEP VDW SHIFTING COEFFICIENT "
+            << alchVdwShiftCoeff << "\n";
+     }
      iout << iINFO << "FEP ELEC. ACTIVE FOR ANNIHILATED "
           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
           << (1 - alchElecLambdaStart) << "\n";
      iout << iINFO << "FEP ELEC. ACTIVE FOR EXNIHILATED "
           << "PARTICLES BETWEEN LAMBDA = "
           << alchElecLambdaStart << " AND LAMBDA = 1\n";
-     iout << iINFO << "FEP VDW ACTIVE FOR ANNIHILATED "
-          << "PARTICLES BETWEEN LAMBDA = "
-          << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
-     iout << iINFO << "FEP VDW ACTIVE FOR EXNIHILATED "
-          << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
-          << alchVdwLambdaEnd << "\n";
+     if (alchWCAOn) {
+       iout << iINFO << "FEP VDW-REPU. ACTIVE FOR ANNIHILATED PARTICLES "
+            << "BETWEEN LAMBDA = " << (1 - alchRepLambdaEnd) << " AND LAMBDA "
+            << "= 1\n";
+       iout << iINFO << "FEP VDW-REPU. ACTIVE FOR EXNIHILATED PARTICLES "
+            << "BETWEEN LAMBDA = 0 AND LAMBDA " << alchRepLambdaEnd << "\n";
+       iout << iINFO << "FEP VDW-ATTR. ACTIVE FOR ANNIHILATED PARTICLES "
+            << "BETWEEN LAMBDA = " << (1 - alchVdwLambdaEnd) << " AND LAMBDA = "
+            << (1 - alchRepLambdaEnd) << "\n";
+       iout << iINFO << "FEP VDW-ATTR. ACTIVE FOR EXNIHILATED PARTICLES "
+            << "BETWEEN LAMBDA = " << alchRepLambdaEnd << " AND LAMBDA = "
+            << alchVdwLambdaEnd << "\n";
+     } else {
+       iout << iINFO << "FEP VDW ACTIVE FOR ANNIHILATED "
+            << "PARTICLES BETWEEN LAMBDA = "
+            << (1 - alchVdwLambdaEnd) << " AND LAMBDA = 1\n";
+       iout << iINFO << "FEP VDW ACTIVE FOR EXNIHILATED "
+            << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
+            << alchVdwLambdaEnd << "\n";
+     }
      iout << iINFO << "FEP BOND ACTIVE FOR ANNIHILATED "
           << "PARTICLES BETWEEN LAMBDA = "
           << (1 - alchBondLambdaEnd) << " AND LAMBDA = 1\n";
      iout << iINFO << "FEP BOND ACTIVE FOR EXNIHILATED "
           << "PARTICLES BETWEEN LAMBDA = 0 AND LAMBDA = "
           << alchBondLambdaEnd << "\n";
-
-     if (alchFepWCADispOn)
-     {
-       iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN DECOMPOSITION (DISPERSION) ON\n";
-     }
-     if (alchFepWCARepuOn)
-     {
-       iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN DECOMPOSITION (REPULSION) ON\n";
-       iout << iINFO << "FEP WEEKS-CHANDLER-ANDERSEN RCUT1 = " 
-            << alchFepWCArcut1 << " , RCUT2 = " 
-            << alchFepWCArcut2  << " AND RCUT3 = " << alchFepWCArcut3 << "\n";
-     }
    }
 //fepe
 
@@ -7219,9 +7127,35 @@ BigReal SimParameters::getElecLambda(const BigReal lambda) {
           : (lambda - alchElecLambdaStart) / (1. - alchElecLambdaStart));
 }
 
+/*
+ * Modifications for WCA decomposition of van der Waal interactions.
+ *
+ * WCA requires that repulsive and attractive components of the vdW 
+ * forces be treated separately. To keep the code clean, the same scaling
+ * function is always used and simply has its behavior modified. However,
+ * the new repluslive scaling only ever gets used when alchWCAOn.
+ */
 BigReal SimParameters::getVdwLambda(const BigReal lambda) {
   // Convenience function for staggered lambda scaling
-  return (lambda >= alchVdwLambdaEnd ? 1. : lambda / alchVdwLambdaEnd);
+  if ( alchWCAOn ) {
+    // Read this with the alias alchRepLambdaEnd --> alchAttLambdaStart.
+    // The second condition is needed when attractive interactions are inactive
+    // for the whole range, otherwise lambda = 0/1 are incorrect.
+    if ( lambda < alchRepLambdaEnd || alchRepLambdaEnd == 1.0 ) {
+      return 0.0;
+    } else if ( lambda >= alchVdwLambdaEnd ) {
+      return 1.0;
+    } else {
+      return (lambda - alchRepLambdaEnd) / (alchVdwLambdaEnd - alchRepLambdaEnd);
+    }
+  } else {
+    return (lambda >= alchVdwLambdaEnd ? 1. : lambda / alchVdwLambdaEnd);
+  }
+}
+
+BigReal SimParameters::getRepLambda(const BigReal lambda) {
+  // Convenience function for staggered lambda scaling
+  return (lambda >= alchRepLambdaEnd ? 1. : lambda / alchRepLambdaEnd);
 }
 
 BigReal SimParameters::getBondLambda(const BigReal lambda) {
index 35f11f7..94b8df0 100644 (file)
@@ -378,10 +378,7 @@ public:
   Bool alchOn;              //  Doing alchemical simulation?
   Bool alchFepOn;           //  Doing alchemical simulation?
   Bool alchThermIntOn;      //  Doing thermodynamic integration?
-  Bool alchFepWCARepuOn;    //  Doing WCA decomposition repulsion interaction?
-  Bool alchFepWCADispOn;    //  Doing WCA decomposition dispersion interaction?
-  Bool alchFepElecOn;       //  Doing electrostatic interaction perturbation?
-  Bool alchFepWhamOn;       //  Doing Wham postprocessing for FEP?
+  Bool alchWCAOn;           //  Using WCA decomposition for vdWs?
   int alchMethod;           //  Which alchemical method to use? fep or ti
   BigReal alchLambda;       //  lambda for dynamics
   BigReal alchLambda2;      //  lambda for comparison
@@ -393,12 +390,6 @@ public:
   BigReal getCurrentLambda2(const int); // getter for alternating lambda2 in IDWS
   int setupIDWS();          //  activates IDWS and sets alchIDWSFreq
   BigReal getLambdaDelta(void); // getter for lambda increment
-  BigReal alchRepLambda;    //  lambda for WCA repulsive interaction
-  BigReal alchDispLambda;   //  lambda for WCA dispersion interaction
-  BigReal alchElecLambda;   //  lambda for electrostatic perturbation
-  BigReal alchFepWCArcut1;  //  rcut1 of WCA decompistion repulsion
-  BigReal alchFepWCArcut2;  //  rcut2 of WCA decomposition repulsion
-  BigReal alchFepWCArcut3;  //  rcut3 of WCA decomposition repulsion
   BigReal alchTemp;         //  temperature for alchemical calculation
   int alchOutFreq;          //  freq. of alchemical output
   Bool alchEnsembleAvg;      //if do ensemble average for the net free energy difference 
@@ -417,6 +408,12 @@ public:
                              //  For annihilated particles the endpoint is 
                              //  (1-alchVdwLambdaEnd)
   BigReal getVdwLambda(const BigReal); // return max[1,x/vdwEnd]
+  BigReal alchRepLambdaEnd;  //  lambda value for endpoint of repulsive vdW
+                             //  interactions of exnihilated particles.
+                             //  For annihilated particles the endpoint is 
+                             //  (1-alchRepLambdaEnd). This also implies the
+                             //  START for attractive vdW interactions.
+  BigReal getRepLambda(const BigReal); // return max[1,x/repEnd]
   BigReal alchBondLambdaEnd; //  lambda value for endpoint of bonded 
                              //  interactions involving exnihilated particles.
                              //  For annihilated particles the endpoint is