Initial version of IDWS 67/4267/9
authorJérôme Hénin <heninj@gmail.com>
Fri, 15 Jun 2018 13:55:29 +0000 (15:55 +0200)
committerJérôme Hénin <heninj@gmail.com>
Wed, 27 Jun 2018 12:50:48 +0000 (14:50 +0200)
Interleaved Double-Wide Sampling, based on an idea from Grace Brannigan.
Switches between lambda2 values, eg. every fullElectFrequency steps.
On-the-fly FEP average is calculated using only forward samples.
Use lib/alch/deinterleave_idws.py to get complete forward and backward data.
Also adds the fep.tcl driver script to lib/alch.

Change-Id: I2faffe5acffb0c8bad18573074c082664f79ced4

19 files changed:
lib/alch/README.txt [new file with mode: 0644]
lib/alch/deinterleave_idws.py [new file with mode: 0755]
lib/alch/fep.tcl [new file with mode: 0644]
src/ComputeAngles.C
src/ComputeAniso.C
src/ComputeBonds.C
src/ComputeCrossterms.C
src/ComputeDihedrals.C
src/ComputeImpropers.C
src/ComputeNonbondedBase.h
src/ComputeNonbondedUtil.C
src/ComputeNonbondedUtil.h
src/ComputePme.C
src/ComputeThole.C
src/Controller.C
src/Controller.h
src/SimParameters.C
src/SimParameters.h
ug/ug_alchemy.tex

diff --git a/lib/alch/README.txt b/lib/alch/README.txt
new file mode 100644 (file)
index 0000000..efc4b04
--- /dev/null
@@ -0,0 +1,23 @@
+####################################################
+#
+# Scripts for Alchemical transformations in NAMD
+#
+####################################################
+
+
+This directory contains scripts that are helpful for either running alchemical
+transformations, or analyzing their results.
+
+
+# fep.tcl
+Driver script to be sourced in a NAMD config file. Provides a concise syntax
+to run alchemical transformations in sequential windows.
+
+
+# deinterleave_idws.py
+A tool for post-processing fepout files obtained with NAMD's
+Interleaved Double-Wide Sampling (option alchLambdaIDWS)
+
+deinterleave_idws.py depends on SciPy.interpolate; it takes a fepout file containing
+interleaved backward and forward energy differences, and splits the data it into distinct
+fepout files that can be processed by other tools e.g. the ParseFEP plugin of VMD.
diff --git a/lib/alch/deinterleave_idws.py b/lib/alch/deinterleave_idws.py
new file mode 100755 (executable)
index 0000000..24817b0
--- /dev/null
@@ -0,0 +1,167 @@
+#!/usr/bin/python3
+
+# Input: fepout file from an interleaved double-wide sampling run from 0 to 1
+# (with the last window from 1 to the previous lambda value)
+# Output:
+# - <base>_fwd.fepout (forward deltaE only)
+# - <base>_bwd.fepout (backward deltaE only)
+# - <base>_bwd_r.fepout (backward data, with order of windows reversed to simulate a
+#   reverse sequential transformation (0 -> 1 -> 0)
+# Process fwd.fepout and bwd_r.fepout with ParseFEP for SOS and BAR estimators
+# Missing samples of deltaE are interpolated
+#
+# Note:
+# The first window samples forward only, the last window backward only
+#
+# Jérôme Hénin <henin@ibpc.fr> (2018)
+
+import sys
+import os.path
+import argparse
+import re
+import math
+from scipy.interpolate import interp1d
+from scipy.interpolate import InterpolatedUnivariateSpline
+
+
+def interpolated_window(steps_window, deltaE_window, lambda1, lambda2):
+  equil = True
+  acc_exp = 0
+  acc_N = 0
+  buf = ''
+  # Create interpolation function for deltaE data
+  func_dE = InterpolatedUnivariateSpline(steps_window, deltaE_window)
+  for t in range(len(steps)):
+    s = steps[t]
+    if equil and s >= equilSteps:
+      buf += '#' + str(s) + ' STEPS OF EQUILIBRATION AT LAMBDA ' + str(0) + ' COMPLETED\n#STARTING COLLECTION OF ENSEMBLE AVERAGE\n'
+      equil = False
+      acc_exp = 0
+      acc_N = 0
+    dE = float(func_dE(s))
+    acc_exp += math.exp(-beta * dE)
+    acc_N += 1
+    buf += 'FepEnergy: {:6d}  {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f} {:14.4f}\n'.format(s, 0, 0, 0, 0, dE, 0., T[t], -RT*math.log(acc_exp / acc_N))
+  return buf, -RT*math.log(acc_exp / acc_N)
+
+
+parser = argparse.ArgumentParser(description='Process fepout data from Interleaved Double-Wide Sampling')
+parser.add_argument('filename', help='fepout file name with lambda from 0 to 1')
+parser.add_argument('-T', type=float, default=300., help='temperature for FEP estimate')
+parser.add_argument('-r', type=bool, default=True, help='reverse order of backward-sampling windows')
+args = parser.parse_args()
+
+basename = os.path.splitext(args.filename)[0]
+outf = open(basename + '_fwd.fepout', 'w')
+
+# Only interpolate deltaE data
+quick = True
+
+# Buffer to hold the contents of the backwards fepout file, by window. Will be written in reverse.
+outb_buffer = []
+
+first_window = True
+last_window = False
+
+l_l1 = []
+l_l2 = []
+l_l3 = []
+
+RT = 1.98720650096e-3 * args.T # in kcal/mol
+beta = 1. / RT
+
+with open(args.filename, 'r') as fepout:
+  for line in fepout:
+    fields = re.split(' +', line.strip())
+    if line[:16] == '#NEW FEP WINDOW:':
+
+      # data will contain all dE data for the window
+      # make two copies, forward and backward
+      deltaE_f = []
+      steps_f = []
+      deltaE_b = []
+      steps_b = []
+      # Common data: step number, temperature
+      steps = []
+      T = []
+
+      lambda1 = float(fields[6])
+      lambda2 = float(fields[8])
+      if len(fields) == 11:
+        lambdaIDWS = float(fields[10])
+      else:
+        lambdaIDWS = -1.
+      print("Lambda window: " + str(lambda1) + " to " + str(lambda2) + (" IDWS: " + str(lambdaIDWS) if lambdaIDWS >= 0 else ""))
+
+      if first_window:
+        if lambda1 == 1:
+          start_from_one = True
+        elif lambda1 == 0:
+          start_from_one = False
+        else:
+          sys.exit('Lambda should start at zero or one, found ' + str(lambda1))
+
+      # Skip the last window for the "forward" output, as it is backward from 1
+      if ((not start_from_one and lambda1 == 1) or (start_from_one and lambda1 == 0)):
+        last_window = True
+      if (not last_window):
+        outf.write("#NEW FEP WINDOW: LAMBDA SET TO {} LAMBDA2 {}\n".format(lambda1, lambda2))
+      else:
+        # special case, last window is backward
+        lambdaIDWS = lambda2
+        lambda2 = -1
+
+      if not first_window:
+        # create new window in the backward output
+        # for every window but the first one, which is forward-only
+        outb_buffer.append("#NEW FEP WINDOW: LAMBDA SET TO {} LAMBDA2 {}\n".format(lambda1, lambdaIDWS))
+
+      # Done processing the header of the first window
+      first_window = False
+
+    if line.strip() == '#STARTING COLLECTION OF ENSEMBLE AVERAGE':
+      equilSteps = steps[-1] + 1 # actually a lower bound of equilSteps if fepOutFreq > 1
+
+    if line[:10] == 'FepEnergy:' and not last_window:
+      steps_f.append(int(fields[1]))
+      deltaE_f.append(float(fields[6]))
+
+    if line[:10] == 'FepE_back:' or (last_window and line[:10] == 'FepEnergy:'):
+      steps_b.append(int(fields[1]))
+      deltaE_b.append(float(fields[6]))
+
+    # Collect all timestep numbers (both fwd and bwd) for interpolation
+    if line[:4] == 'FepE':
+      steps.append(int(fields[1]))
+      T.append(float(fields[8]))
+
+    if line[:19] == '#Free energy change':
+      # End of window
+      # Create interpolated function from data and write new data for the whole window
+
+      if not last_window:
+        buf, dG = interpolated_window(steps_f, deltaE_f, lambda1, lambda2)
+        outf.write(buf)
+        outf.write(line)
+
+      if len(steps_b) > 0: # only if we do have backward data in this window
+        buf, dG = interpolated_window(steps_b, deltaE_b, lambda1, lambdaIDWS)
+        outb_buffer[-1] += buf
+        outb_buffer[-1] += '#Free energy change for lambda window [ {} {} ] is {} ; net change until now is 0.\n'.format(lambda1, lambda2, dG)
+
+
+outf.close()
+
+if args.r:
+  # Write backward windows in reverse order
+  outb_reversed = open(basename + '_bwd.fepout', 'w')
+  for string in reversed(outb_buffer):
+    outb_reversed.write(string)
+  outb_reversed.close()
+else:
+  # Write output in same order as run
+  outb = open(basename + '_bwd.fepout', 'w')
+  for string in outb_buffer:
+    outb.write(string)
+  outb.close()
+
diff --git a/lib/alch/fep.tcl b/lib/alch/fep.tcl
new file mode 100644 (file)
index 0000000..a84b7e7
--- /dev/null
@@ -0,0 +1,270 @@
+##############################################################
+# FEP SCRIPT
+# Jerome Henin <henin@ibpc.fr>
+#
+# Changes:
+# 2018-04-18: added interleaved double-wide sampling (IDWS)
+# 2010-04-24: added runFEPmin
+# 2009-11-17: changed for NAMD 2.7 keywords
+# 2008-06-25: added TI routines
+# 2007-11-01: fixed runFEP to handle backwards transformations
+#             (i.e. dLambda < 0)
+##############################################################
+
+##############################################################
+# Example NAMD input:
+#
+# source fep.tcl
+#
+# alch                  on
+# alchFile              system.fep
+# alchCol               B
+# alchOutFreq           10
+# alchOutFile           system.fepout
+# alchEquilSteps        500
+#
+# set nSteps      5000
+# set init {0 0.05 0.1}
+# set end {0.9 0.95 1.0}
+#
+# runFEPlist $init $nSteps
+# runFEP 0.1 0.9 0.1 $nSteps
+# runFEPlist $end $nSteps
+#
+## Alternately, in one step:
+#
+# runFEPlist [concat $init [FEPlist 0.1 0.9 0.1] $end] $nSteps
+#
+##############################################################
+
+##############################################################
+# Special usage for Interleaved Double-Wide sampling
+#
+# Example of a piecewise calculation with restarts
+# and a nonlinear lambda schedule
+#
+## Run individual points 0, 0.05 then the series from 0.1 to 0.4
+#
+# runFEPlist [concat {0. 0.05} [FEPlist 0.1 0.4 0.1]] $numSteps true
+#
+## Continue series from 0.5 to 0.9, sampling backward dE from 0.4
+#
+# runFEPlist [FEPlist 0.5 0.9 0.1] $numSteps true 0.4
+#
+## Add two values 0.95 and 1, sampling backward dE from 0.9
+## (automatically adds final backward window from 1. to 0.95)
+#
+# runFEPlist {0.95 1.} $numSteps true 0.9
+#
+##############################################################
+
+##############################################################
+# proc runFEPlist { lambdaList nSteps {IDWS} {prev_lambda} }
+#
+# Run n FEP windows joining (n + 1) lambda-points
+# Provide prev_lambda value if continuing a sequential
+# transformation and using IDWS
+##############################################################
+
+proc runFEPlist { lambdaList nSteps { IDWS false } { prev_lambda -1 } } {
+    set epsilon 1e-15
+
+    # Keep track of window number
+    global win
+    if {![info exists win]} {
+      set win 1
+    }
+
+    set l1 [lindex $lambdaList 0]
+    foreach l2 [lrange $lambdaList 1 end] {
+      print [format "Running FEP window %3s: Lambda1 %-6s Lambda2 %-6s \[dLambda %-6s\]"\
+        $win $l1 $l2 [expr $l2 - $l1]]
+      firsttimestep 0
+      alchLambda       $l1
+      alchLambda2      $l2
+      
+      if { $IDWS && ($prev_lambda >= 0.) } {
+        alchLambdaIDWS $prev_lambda
+      }
+      run $nSteps
+
+      # Keep track of previous value to set is as target for backward calculation in IDWS
+      set prev_lambda $l1
+      set l1 $l2
+      incr win
+    }
+
+    if { $IDWS && ($l1 > [expr {1. - $epsilon}] || $l1 < $epsilon)} {
+      # If the list ends at 1 or zero, we add a final window, which is backward from the end point
+      # to complete double-wide sampling
+      # this will be look like "forward" sampling in the fepout file ("FepEnergy:" keyword)
+      print [format "Running FEP window %3s: Lambda1 %-6s Lambda2 %-6s \[dLambda %-6s\]"\
+        $win $l1 $l2 [expr $l2 - $l1]]
+      firsttimestep 0
+      alchLambda       $l1
+      alchLambda2      $prev_lambda
+      alchLambdaIDWS   -1
+      run $nSteps
+    }
+}
+
+
+##############################################################
+# proc runFEP { start stop dLambda nSteps {IDWS} }
+#
+# run FEP windows of width dLambda between values start and stop
+##############################################################
+
+proc runFEP { start stop dLambda nSteps { IDWS false } } {
+
+    runFEPlist [FEPlist $start $stop $dLambda] $nSteps $IDWS
+}
+
+
+##############################################################
+# proc FEPlist { start stop dLambda nSteps }
+#
+# Create list of FEP windows
+##############################################################
+
+proc FEPlist { start stop dLambda } {
+    set epsilon 1e-15
+
+    if { ($stop < $start) && ($dLambda > 0) } {
+      set dLambda [expr {-$dLambda}]
+    }
+
+    if { $start == $stop } {
+      set ll [list $start $start]
+    } else {
+      set ll [list $start]
+      set l2 [increment $start $dLambda]
+
+      if { $dLambda > 0} {
+        # A small workaround for numerical rounding errors
+        while { [expr {$l2 <= ($stop + $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      } else {
+        while { [expr {$l2 >= ($stop - $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      }
+    }
+
+    return $ll
+}
+
+
+##############################################################
+##############################################################
+
+proc runFEPmin { start stop dLambda nSteps nMinSteps temp} {
+    set epsilon 1e-15
+
+    if { ($stop < $start) && ($dLambda > 0) } {
+      set dLambda [expr {-$dLambda}]
+    }
+
+    if { $start == $stop } {
+      set ll [list $start $start]
+    } else {
+      set ll [list $start]
+      set l2 [increment $start $dLambda]
+
+      if { $dLambda > 0} {
+        # A small workaround for numerical rounding errors
+        while { [expr {$l2 <= ($stop + $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      } else {
+        while { [expr {$l2 >= ($stop - $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      }
+    }
+
+    if { $nMinSteps > 0 } { 
+      alchLambda       $start
+      alchLambda2      $start
+      minimize $nMinSteps
+      reinitvels $temp
+    }
+
+    runFEPlist $ll $nSteps
+}
+
+##############################################################
+##############################################################
+
+proc runTIlist { lambdaList nSteps } {
+    # Keep track of window number
+    global win
+    if {![info exists win]} {
+           set win 1
+    }
+
+    foreach l $lambdaList {
+           print [format "Running TI window %3s: Lambda %-6s " $win $l ]
+           firsttimestep 0
+           alchLambda       $l
+           run $nSteps
+           incr win
+    }
+}
+
+
+##############################################################
+##############################################################
+
+proc runTI { start stop dLambda nSteps } {
+    set epsilon 1e-15
+
+    if { ($stop < $start) && ($dLambda > 0) } {
+      set dLambda [expr {-$dLambda}]
+    }
+
+    if { $start == $stop } {
+      set ll [list $start $start]
+    } else {
+      set ll [list $start]
+      set l2 [increment $start $dLambda]
+
+      if { $dLambda > 0} {
+        # A small workaround for numerical rounding errors
+        while { [expr {$l2 <= ($stop + $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      } else {
+        while { [expr {$l2 >= ($stop - $epsilon) } ] } {
+          lappend ll $l2
+          set l2 [increment $l2 $dLambda]
+        }
+      }
+    }
+
+    runTIlist $ll $nSteps
+}
+
+##############################################################
+# Increment lambda and try to correct truncation errors around
+# 0 and 1
+##############################################################
+
+proc increment { lambda dLambda } {
+    set epsilon 1e-15
+    set new [expr { $lambda + $dLambda }]
+
+    if { [expr $new > - $epsilon && $new < $epsilon] } {
+      return 0.0
+    }
+    if { [expr ($new - 1) > - $epsilon && ($new - 1) < $epsilon] } {
+      return 1.0
+    }
+    return $new
+}
index 1499ce1..93730d8 100644 (file)
@@ -52,7 +52,7 @@ void AngleElem::computeForce(AngleElem *tuples, int ntuple, BigReal *reduction,
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index 4c6db52..a38c702 100644 (file)
@@ -47,7 +47,7 @@ void AnisoElem::computeForce(AnisoElem *tuples, int ntuple, BigReal *reduction,
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index bc3c2e7..8692634 100644 (file)
@@ -46,7 +46,7 @@ void BondElem::computeForce(BondElem *tuples, int ntuple, BigReal *reduction,
  //fepb BKR
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index f4c99f2..6a46da6 100644 (file)
@@ -58,7 +58,7 @@ void CrosstermElem::computeForce(CrosstermElem *tuples, int ntuple, BigReal *red
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index 52612c1..3aa776a 100644 (file)
@@ -46,7 +46,7 @@ void DihedralElem::computeForce(DihedralElem *tuples, int ntuple, BigReal *reduc
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index 37a3a02..15db480 100644 (file)
@@ -46,7 +46,7 @@ void ImproperElem::computeForce(ImproperElem *tuples, int ntuple, BigReal *reduc
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
index a6c316f..2e04e8f 100644 (file)
@@ -404,7 +404,7 @@ void ComputeNonbondedUtil :: NAME
     BigReal elecLambdaUp = simParams->getElecLambda(lambdaUp);
     BigReal vdwLambdaUp = simParams->getVdwLambda(lambdaUp);
     BigReal vdwShiftUp = alchVdwShiftCoeff*(1 - vdwLambdaUp);
-    FEP(BigReal lambda2Up = ComputeNonbondedUtil::alchLambda2;)
+    FEP(BigReal lambda2Up = simParams->getCurrentLambda2(params->step);)
     FEP(BigReal elecLambda2Up = simParams->getElecLambda(lambda2Up);)
     FEP(BigReal vdwLambda2Up = simParams->getVdwLambda(lambda2Up);)
     FEP(BigReal vdwShift2Up = alchVdwShiftCoeff*(1 - vdwLambda2Up);)
@@ -430,7 +430,7 @@ void ComputeNonbondedUtil :: NAME
     BigReal elecLambdaDown = simParams->getElecLambda(lambdaDown);
     BigReal vdwLambdaDown = simParams->getVdwLambda(lambdaDown);
     BigReal vdwShiftDown = alchVdwShiftCoeff*(1 - vdwLambdaDown);
-    FEP(BigReal lambda2Down = 1 - ComputeNonbondedUtil::alchLambda2;)
+    FEP(BigReal lambda2Down = 1 - simParams->getCurrentLambda2(params->step);)
     FEP(BigReal elecLambda2Down = simParams->getElecLambda(lambda2Down);)
     FEP(BigReal vdwLambda2Down = simParams->getVdwLambda(lambda2Down);)
     FEP(BigReal vdwShift2Down = alchVdwShiftCoeff*(1 - vdwLambda2Down);)
index efb8479..9907347 100644 (file)
@@ -120,7 +120,6 @@ Bool      ComputeNonbondedUtil::Fep_Wham;
 BigReal   ComputeNonbondedUtil::WCA_rcut1;
 BigReal   ComputeNonbondedUtil::WCA_rcut2;
 BigReal   ComputeNonbondedUtil::WCA_rcut3;
-BigReal   ComputeNonbondedUtil::alchLambda2;
 BigReal   ComputeNonbondedUtil::alchRepLambda;
 BigReal   ComputeNonbondedUtil::alchDispLambda;
 BigReal   ComputeNonbondedUtil::alchElecLambda;
@@ -300,7 +299,6 @@ void ComputeNonbondedUtil::select(void)
   Fep_ElecOn = simParams->alchFepElecOn;
   Fep_Wham = simParams->alchFepWhamOn;
   alchThermIntOn = simParams->alchThermIntOn;
-  alchLambda2 = 0;
   lesOn = simParams->lesOn;
   lesScaling = lesFactor = 0;
   Bool tabulatedEnergies = simParams->tabulatedEnergies;
@@ -350,7 +348,6 @@ void ComputeNonbondedUtil::select(void)
 #ifdef NAMD_CUDA
     NAMD_die("Alchemical free-energy perturbation is not supported in CUDA version");
 #endif
-    alchLambda2 = simParams->alchLambda2;
     ComputeNonbondedUtil::calcPair = calc_pair_energy_fep;
     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_fep;
     ComputeNonbondedUtil::calcSelf = calc_self_energy_fep;
@@ -371,7 +368,6 @@ void ComputeNonbondedUtil::select(void)
 #ifdef NAMD_CUDA
     NAMD_die("Alchemical thermodynamic integration is not supported in CUDA version");
 #endif
-    alchLambda2 = simParams->alchLambda2;
     ComputeNonbondedUtil::calcPair = calc_pair_ti;
     ComputeNonbondedUtil::calcPairEnergy = calc_pair_energy_ti;
     ComputeNonbondedUtil::calcSelf = calc_self_ti;
index 4e296e8..228578b 100644 (file)
@@ -355,7 +355,6 @@ public:
 //sd-db
   static Bool alchFepOn;
   static Bool alchThermIntOn;
-  static BigReal alchLambda2;
   static BigReal alchVdwShiftCoeff;
   static Bool vdwForceSwitching;
   static Bool LJcorrection;
index bd24465..375431b 100644 (file)
@@ -495,6 +495,7 @@ private:
   int qsize, fsize, bsize;
   int offload;
   BigReal alchLambda;  // set on each step in ComputePme::ungridForces()
+  BigReal alchLambda2; // set on each step in ComputePme::ungridForces()
 
   float **q_arr;
   // q_list and q_count not used for offload
@@ -898,6 +899,7 @@ void ComputePmeMgr::initialize(CkQdMsg *msg) {
 #endif
 
   alchLambda = -1.;  // illegal value to catch if not updated
+  alchLambda2 = -1.;
   useBarrier = simParams->PMEBarrier;
 
   if ( numGrids != 1 || simParams->PMEPencils == 0 ) usePencils = 0;
@@ -4052,6 +4054,8 @@ void ComputePme::ungridForces() {
         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);
         }
@@ -4111,6 +4115,8 @@ void ComputePme::ungridForces() {
          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;
          }
@@ -4256,8 +4262,8 @@ void ComputePmeMgr::submitReductions() {
           }
         }
         else {
-          elecLambda2Up = simParams->getElecLambda(simParams->alchLambda2);
-          elecLambda2Down = simParams->getElecLambda(1.-simParams->alchLambda2);
+          elecLambda2Up = simParams->getElecLambda(alchLambda2);
+          elecLambda2Down = simParams->getElecLambda(1.-alchLambda2);
         }
         
         if ( g == 0 ) scale2 = elecLambda2Up;
index 2424d74..3ecf793 100644 (file)
@@ -47,7 +47,7 @@ void TholeElem::computeForce(TholeElem *tuples, int ntuple, BigReal *reduction,
  SimParameters *const simParams = Node::Object()->simParameters;
  const int step = tuples[0].p[0]->p->flags.step;
  const BigReal alchLambda = simParams->getCurrentLambda(step);
- const BigReal alchLambda2 = simParams->alchLambda2;
+ const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda2);
index 5d045c9..89ca427 100644 (file)
@@ -1276,11 +1276,15 @@ void Controller::printFepMessage(int step)
       && !simParams->alchLambdaFreq) {
     const BigReal alchLambda = simParams->alchLambda;
     const BigReal alchLambda2 = simParams->alchLambda2;
+    const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
     const BigReal alchTemp = simParams->alchTemp;
     const int alchEquilSteps = simParams->alchEquilSteps;
     iout << "FEP: RESETTING FOR NEW FEP WINDOW "
-         << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2
-         << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
+         << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
+    if ( alchLambdaIDWS >= 0. ) {
+      iout << " LAMBDA_IDWS " << alchLambdaIDWS;
+    }
+    iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
          << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
          << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp 
          << " K FOR FEP CALCULATION\n" << endi;
@@ -3586,7 +3590,6 @@ void Controller::outputFepEnergy(int step) {
   const int stepInRun = step - simParams->firstTimestep;
   const int alchEquilSteps = simParams->alchEquilSteps;
   const BigReal alchLambda = simParams->alchLambda;
-  const BigReal alchLambda2 = simParams->alchLambda2;
   const bool alchEnsembleAvg = simParams->alchEnsembleAvg; 
   const bool FepWhamOn = simParams->alchFepWhamOn;
 
@@ -3599,7 +3602,8 @@ void Controller::outputFepEnergy(int step) {
        electEnergy - electEnergySlow - ljEnergy;
   BigReal RT = BOLTZMANN * simParams->alchTemp;
 
-  if (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;
@@ -3628,10 +3632,14 @@ void Controller::outputFepEnergy(int step) {
         }
       }
       if(!step){
-        if(!FepWhamOn){ 
+        if(!FepWhamOn){
           fepFile << "#NEW FEP WINDOW: "
                   << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " 
-                  << alchLambda2 << std::endl;
+                  << simParams->alchLambda2;
+          if ( simParams->alchLambdaIDWS >= 0. ) {
+            fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
+          }
+          fepFile << std::endl;
         }
       }
     }
@@ -3647,7 +3655,7 @@ void Controller::outputFepEnergy(int step) {
     if (alchEnsembleAvg && (step == simParams->N)) {
       fepSum = fepSum + dG;
       fepFile << "#Free energy change for lambda window [ " << alchLambda
-              << " " << alchLambda2 << " ] is " << dG 
+              << " " << simParams->alchLambda2 << " ] is " << dG 
               << " ; net change until now is " << fepSum << std::endl;
       fepFile.flush();
     }
@@ -3742,7 +3750,7 @@ void Controller::outputTiEnergy(int step) {
 
       if (alchLambdaFreq > 0) {
         tiFile << "#ALCHEMICAL SWITCHING ACTIVE " 
-               << simParams->alchLambda << " --> " << simParams->alchLambda2
+               << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
                << "\n#LAMBDA SCHEDULE: " 
                << "dL: " << simParams->getLambdaDelta() 
                << " Freq: " << alchLambdaFreq;
@@ -3834,7 +3842,13 @@ void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
 
   if(stepInRun){
     if(!FepWhamOn){
-      fepFile << FEPTITLE(step);
+      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);
index 5e3d325..a63a88d 100644 (file)
@@ -337,6 +337,13 @@ static char *FEPTITLE(int X)
   return tmp_string;
 }
 
+static char *FEPTITLE_BACK(int X)
+{
+  static char tmp_string[21];
+  sprintf(tmp_string, "FepE_back: %6d ",X);
+  return tmp_string;
+}
+
 static char *FEPTITLE2(int X)
 {
   static char tmp_string[21];
index 27b6a2e..5c42860 100644 (file)
@@ -359,8 +359,26 @@ void SimParameters::scriptSet(const char *param, const char *value) {
     return;
   }
 
+  if ( ! strncasecmp(param,"alchLambdaIDWS",MAX_SCRIPT_PARAM_SIZE) ) {
+    alchLambdaIDWS = atof(value);
+    if ( alchLambdaIDWS > 1. ) {
+      NAMD_die("alchLambdaIDWS should be either in the range [0.0, 1.0], or negative (disabled).\n");
+    }
+    // Switch lambda2 every other cycle of fullElectFrequency steps
+    // or every other nonbondedFrequency steps if undefined
+    // or every alchOutFreq steps if larger (no need to switch faster that we output)
+    alchIDWSfreq = fullElectFrequency > 0 ? fullElectFrequency : nonbondedFrequency;
+    if ( alchOutFreq > alchIDWSfreq )
+      alchIDWSfreq = alchOutFreq;
+    ComputeNonbondedUtil::select();
+    return;
+  }
+
   if ( ! strncasecmp(param,"alchLambdaFreq",MAX_SCRIPT_PARAM_SIZE) ) {
     alchLambdaFreq = atoi(value);
+    if ( alchLambdaIDWS >= 0 ) {
+      NAMD_die("alchLambdaIDWS and alchLambdaFreq are not compatible.\n");
+    }
     ComputeNonbondedUtil::select();
     return;
   }
@@ -1119,6 +1137,8 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
      PARSE_STRING);
    opts.optional("alch", "alchLambda2", "Alchemical coupling comparison value",
      &alchLambda2, -1);
+   opts.optional("alch", "alchLambdaIDWS", "Alchemical coupling comparison value for interleaved double-wide sampling",
+     &alchLambdaIDWS, -1);
    opts.optional("alch", "alchLambdaFreq",
      "Frequency of increasing coupling parameter value", &alchLambdaFreq, 0);
    opts.range("alchLambdaFreq", NOT_NEGATIVE);
@@ -3598,7 +3618,19 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
          strcat(alchOutFile, ".ti"); 
        } 
      }
-   } 
+     if (alchLambdaIDWS >= 0.) {
+       if ( alchLambdaIDWS > 1. ) {
+         NAMD_die("alchLambdaIDWS should be either in the range [0.0, 1.0], or negative (disabled).\n");
+      }
+       // Switch lambda2 every other cycle of fullElectFrequency steps
+       // or every other nonbondedFrequency steps if undefined
+       // or every alchOutFreq steps if larger (no need to switch faster that we output)
+       alchIDWSfreq = fullElectFrequency > 0 ? fullElectFrequency : nonbondedFrequency;
+       if ( alchOutFreq > alchIDWSfreq )
+         alchIDWSfreq = alchOutFreq;
+     }
+   }
+        
 //fepe
 
    if ( alchOn && alchFepOn && alchThermIntOn )
@@ -5059,6 +5091,10 @@ if ( openatomOn )
           << alchLambda << "\n";
      iout << iINFO << "FEP COMPARISON LAMBDA VALUE  "
           << alchLambda2 << "\n";
+     if (alchLambdaIDWS >= 0.) {
+        iout << iINFO << "FEP ALTERNATE COMPARISON LAMBDA VALUE  "
+          << alchLambdaIDWS << "\n";
+     }
      if (alchLambdaFreq > 0) {
        iout << iINFO << "FEP CURRENT LAMBDA VALUE SET TO INCREASE IN EVERY  "
             << alchLambdaFreq << " STEPS\n";
@@ -7087,6 +7123,19 @@ void SimParameters::receive_SimParameters(MIStream *msg)
 }
 /*      END OF FUNCTION receive_SimParameters  */
 
+
+//fepb IDWS
+BigReal SimParameters::getCurrentLambda2(const int step) {
+  if ( alchLambdaIDWS >= 0. ) {
+    const BigReal lambda2 = ( (step / alchIDWSfreq) % 2 == 1 ) ? alchLambda2 : alchLambdaIDWS;
+    return lambda2;
+  } else {
+    return alchLambda2;
+  }
+}
+//fepe IDWS
+
+
 //fepb BKR
 BigReal SimParameters::getCurrentLambda(const int step) {
   /*Get lambda at the current step. 
index 1a6b611..7c4cab0 100644 (file)
@@ -385,9 +385,12 @@ public:
   int alchMethod;           //  Which alchemical method to use? fep or ti
   BigReal alchLambda;       //  lambda for dynamics
   BigReal alchLambda2;      //  lambda for comparison
+  BigReal alchLambdaIDWS;   //  alternate lambda for interleaved double-wide sampling
+  int alchIDWSfreq;         //  freq with which lambda2 changes to lambdaIDWS
   int alchLambdaFreq;       //  freq. (in steps) with which lambda changes 
                             //  from alchLambda to alchLambda2
   BigReal getCurrentLambda(const int); // getter for changing lambda
+  BigReal getCurrentLambda2(const int); // getter for alternating lambda2 in IDWS
   BigReal getLambdaDelta(void); // getter for lambda increment
   BigReal alchRepLambda;    //  lambda for WCA repulsive interaction
   BigReal alchDispLambda;   //  lambda for WCA dispersion interaction
index 644223c..0abbf2a 100644 (file)
@@ -355,6 +355,16 @@ total free energy difference may be determined. }
 %\item
 %\NAMDCONF{alchLambdaMinusDelta}{Backward projected value of the coupling parameter} {positive decimal between 0.0 and 1.0} {The {\tt alchLambdaMinusDelta} value corresponds to the coupling parameter to be used for sampling in the previous window, assuming that a double-wide sampling simulation is being performed, whereby both the forward and the backward transformations are computed explicitly.  The free energy difference between {\tt alchLambdaMinusDelta} and {\tt alchLambda} is calculated.  Through simulations at progressive values of {\tt alchLambda} and {\tt alchLambdaMinusDelta} the total free energy difference may be determined.}
 
+\item
+\NAMDCONF{alchLambdaIDWS}{Backward value of the coupling parameter for Interleaved Double-Wide Sampling}
+{decimal smaller than 1.0}
+{Setting this parameter between 0 and 1 activates Interleaved Double-Wide Sampling (IDWS), whereby the target lambda value alternates between {\tt alchLambda2} and {\tt alchLambdaIDWS}.
+The switch occurs every {\tt fullElectFrequency} steps if defined, or {\tt nonbondedFrequency} otherwise.
+Setting this parameter to a negative value (including between run statements) disables IDWS.
+When IDWS is active, the alchemy output file contains {\tt FepEnergy} line headers for the forward energy differences, and {\tt FepE\_Back} for backward energy differences.
+FEP free energy estimates are based on forward data only.
+This file may be postprocessed with deinterleave\_idws.py to produce separate files for the forward and backwards samples.}
+
 \item
 \NAMDCONFWDEF{alchEquilSteps}{Number of equilibration steps in a window,
 prior to data collection}