Lambda-ABF implementation for the NAMD interface

colvartools/abmd.tcl

@@ -53,7 +53,7 @@ proc calc_z_max_gradient { args } { return 0 }
 proc calc_z_min_gradient { args } { return 0 }
 
 proc setup_ABMD { colvar force_k z_stop {direction up} } {
-  # cv config "scriptedColvarForces on"
+  cv config "scriptedColvarForces on"
 
   namespace eval ::ABMD {}
   set ::ABMD::cvname $colvar

doc/colvars-code-refs.bib

@@ -230,6 +230,19 @@ @article{Humphrey1996
   url = {https://doi.org/10.1016/0263-7855(96)00018-5}
 }
 
+% alchLambda colvar component
+% alchFLambda colvar component
+% Tinker-HP interface
+@misc{Lagardere2023,
+      title={Lambda-ABF: Simplified, Accurate and Cost-effective Alchemical Free Energy Computations},
+      author={Louis Lagard\`ere and Lise Maurin and Olivier Adjoua and Krystel El Hage and Pierre Monmarch\'e and Jean-Philip Piquemal and J\'er\^ome H\'enin},
+      year={2023},
+      eprint={2307.08006},
+      archivePrefix={arXiv},
+      primaryClass={physics.chem-ph},
+      url = {https://arxiv.org/abs/2307.08006}
+}
+
 % eABF implementation
 % CZAR eABF estimator
 @article{Lesage2017,

doc/colvars-refman-main.tex

@@ -3007,34 +3007,82 @@
 \cvsubsubsec{\texttt{alchLambda}: alchemical lambda parameter.}{sec:cvc_alchlambda}
 \labelkey{colvar|alchLambda}
 
-The \texttt{alchLambda~\{\}} block defines a component returning the alchemical lambda parameter
-of the back-end simulation in \MDENGINE{}, provided that it is enabled using the appropriate options
-of \MDENGINE{}. This coordinate is obtained from the back-end at the beginning of a simulation, and
-synchronized from Colvars to \MDENGINE{} at every time step. This enables lambda-dynamics simulations.
-The \texttt{alchLambda} block itself takes no parameters, it should be left empty.
-In contrast, the \texttt{colvar} block containing it should define the relevant extended-system parameters
-to enable lambda dynamics, primarily \refkey{extendedMass}{colvar|extendedMass}:
+The \texttt{alchLambda} block defines a component returning the alchemical lambda parameter
+of the back-end simulation in \MDENGINE{}, provided that it is enabled within that software.
+\cvnamdonly{Specifically, the flag \texttt{alchOn} should be enabled, and \texttt{alchType}
+should be set to \texttt{TI}.}
+The current lambda value is obtained from the back-end at the beginning of a simulation, and
+synchronized from Colvars to \MDENGINE{} at every following time step.
+This enables lambda-dynamics simulations as documented in detail in \cite{Lagardere2023}.
+The \texttt{alchLambda} keyword itself takes no parameters.
+
+To run lambda-dynamics, the enclosing \texttt{colvar} block should contain extended-system
+parameters (\ref{sec:colvar_extended}) as in the following example:
 
 \begin{cvexampleinput}
 \-colvar~\{\\
 \-~~name~lambda\\
 \-~~extendedLagrangian~on\\
-\-~~extendedMass~400\\
-\-\\
-\-~~alchLambda~\{~\#~Keep~the~line~break\\
-\-~~\}\\
+\-~~extendedMass~150000\\
+\-~~extendedLangevinDamping~1000\\
+\\
+\-~~lowerBoundary~0\\
+\-~~upperBoundary~1\\
+\-~~reflectingLowerBoundary\\
+\-~~reflectingUpperBoundary\\
+\\
+\-~~alchLambda\\
 \-\}
 \end{cvexampleinput}
 
+These parameter values have been found  to give stable dynamics for a wide range of applications\cite{Lagardere2023}.
+Reflecting boundary conditions are recommended for lambda-dynamics to ensure that the alchemical
+variable remains precisely within its physical range.
+
+\cvnamdonly{
+The NAMD configuration file should also contain the parameters for setting up the alchemical simulation, for example:
+\begin{cvexampleinput}
+computeEnergies~~~~~~~~~1\\
+\\
+alch~~~~~~~~~~~~~~~~~~~~on\\
+alchType~~~~~~~~~~~~~~~~TI\\
+alchFile~~~~~~~~~~~~~~~~\$\{alchFile\}.pdb\\
+alchCol~~~~~~~~~~~~~~~~~B\\
+alchOutFile~~~~~~~~~~~~~\$\{outName\}.fepout\\
+alchOutFreq~~~~~~~~~~~~~100\\
+alchVdwLambdaEnd~~~~~~~~0.5\\
+alchElecLambdaStart~~~~~0.5\\
+alchVdWShiftCoeff~~~~~~~6.0\\
+alchlambda~~~~~~~~~~~~~~0\\
+\end{cvexampleinput}
+
+Here, the value of 1 for \texttt{computeEnergies} and TI for \texttt{alchType} are mandatory.
+All other parameters can be chosen freely as documented in the NAMD User's Guide.
+
+\paragraph{Special considerations for running lambda-dynamics in NAMD/Colvars.}
+As of NAMD 3.0, lambda-dynamics is not compatible with multiple timestepping. We recommend applying
+hydrogen mass repartitioning (HMR) to the topology and increasing the base timestep to 2~fs
+to limit the performance impact of disabling multiple timestepping.
+When running in parallel on CPUs, synchronization issues prevent lambda-dynamics trajectories from being
+fully reproducible.
+This can be circumvented by running NAMD in GPU-resident mode.
+} % end of \cvnamdonly
+
+
+
 \cvsubsubsec{\texttt{alchFlambda}: Force on the alchemical lambda parameter.}{sec:cvc_alchFlambda}
 \labelkey{colvar|alchFlambda}
 
 The \texttt{alchFlambda~\{\}} block defines a component returning the force exterted on the alchemical
 lambda parameter of the back-end simulation in \MDENGINE{}, provided that it is enabled using the appropriate options
 of \MDENGINE{}. This coordinate is obtained from the back-end at each time step of a simulation.
-The \texttt{alchfLambda} block itself takes no parameters, it should be left empty.
+The \texttt{alchFlambda} block itself takes no parameters, it should be left empty.
+
+
 } % end of \cvalchlambdaonly
 
+
+
 \cvsubsec{Raw data: building blocks for custom functions}{sec:cvc_raw}
 
 \cvsubsubsec{\texttt{cartesian}: vector of atomic Cartesian coordinates.}{sec:cvc_cartesian}
@@ -4473,7 +4521,7 @@
     positive decimal}{%
     This parameter specifies the fictitious mass in the case of an alchemical
     variable used in lambda-dynamics (\refkey{alchLambda}{colvar|alchLambda}).
-    Note that the dimension is not that of a usual mass. Its unit is the energy unit
+    Note that the dimension is not that of a usual mass: its unit is the energy unit
     (\ref{sec:units}) times fs$^2$.
     }
 }

doc/colvars-refman-namd.tex

@@ -32,7 +32,7 @@
 \newcommand{\cvnamebasedonly}[1]{#1}
 
 % only in programs were the alchemical coordinate (alchLambda) is available
-\newcommand{\cvalchlambdaonly}[1]{}
+\newcommand{\cvalchlambdaonly}[1]{#1}
 
 % File output prefixes
 \newcommand{\outputName}{\emph{outputName}}

lammps/src/COLVARS/colvarproxy_lammps.h

@@ -64,7 +64,8 @@ class colvarproxy_lammps : public colvarproxy {
  public:
 
   bool total_forces_enabled() const override { return total_force_requested; };
-  bool total_forces_same_step() const override { return true; };
+  // Total forces are saved at end of step, only processed at the next step
+  bool total_forces_same_step() const override { return false; };
   bool want_exit() const { return do_exit; };
 
   // perform colvars computation. returns biasing energy

lammps/tests/library/000_distance-grid_abf/AutoDiff/test.colvars.out

@@ -1,10 +1,12 @@
-colvars: Creating proxy instance
 colvars: ----------------------------------------------------------------------
 colvars: Please cite Fiorin et al, Mol Phys 2013:
-colvars:  https://dx.doi.org/10.1080/00268976.2013.813594
-colvars: in any publication based on this calculation.
+colvars:   https://doi.org/10.1080/00268976.2013.813594
+colvars: as well as all other papers listed below for individual features used.
 colvars: This version was built with the C++11 standard or higher.
-colvars: Using LAMMPS interface, version "2020-04-07".
+colvars: Summary of compile-time features available in this build:
+colvars:   - SMP parallelism: enabled (num. threads = 2)
+colvars:   - Lepton custom functions: available
+colvars:   - Tcl interpreter: not available
 colvars: ----------------------------------------------------------------------
 colvars: Reading new configuration from file "test.in":
 colvars: # units = "" [default]
@@ -56,27 +58,27 @@ colvars:     # forceNoPBC = off [default]
 colvars:     # scalable = on [default]
 colvars:       Initializing atom group "group1".
 colvars:       # name = "" [default]
-colvars:       # centerReference = off [default]
-colvars:       # rotateReference = off [default]
+colvars:       # centerToOrigin = off [default]
+colvars:       # centerToReference = off [default]
+colvars:       # rotateToReference = off [default]
 colvars:       # atomsOfGroup = "" [default]
 colvars:       # indexGroup = "group1"
 colvars:       # psfSegID =  [default]
 colvars:       # atomsFile = "" [default]
 colvars:       # dummyAtom = ( 0 , 0 , 0 ) [default]
-colvars:       # enableForces = on [default]
 colvars:       # enableFitGradients = on [default]
 colvars:       # printAtomIDs = off [default]
 colvars:       Atom group "group1" defined with 4 atoms requested.
 colvars:       Initializing atom group "group2".
 colvars:       # name = "" [default]
-colvars:       # centerReference = off [default]
-colvars:       # rotateReference = off [default]
+colvars:       # centerToOrigin = off [default]
+colvars:       # centerToReference = off [default]
+colvars:       # rotateToReference = off [default]
 colvars:       # atomsOfGroup = "" [default]
 colvars:       # indexGroup = "group2"
 colvars:       # psfSegID =  [default]
 colvars:       # atomsFile = "" [default]
 colvars:       # dummyAtom = ( 0 , 0 , 0 ) [default]
-colvars:       # enableForces = on [default]
 colvars:       # enableFitGradients = on [default]
 colvars:       # printAtomIDs = off [default]
 colvars:       Atom group "group2" defined with 4 atoms requested.
@@ -104,19 +106,20 @@ colvars: ----------------------------------------------------------------------
 colvars:   Initializing a new "abf" instance.
 colvars:   # name = "abf1" [default]
 colvars:   # colvars = { one }
-colvars:   # zeroStepData = off [default]
+colvars:   # stepZeroData = off [default]
 colvars:   # outputEnergy = off [default]
 colvars:   # outputFreq = 10 [default]
 colvars:   # timeStepFactor = 1 [default]
 colvars:   WARNING: ABF should not be run without a thermostat or at 0 Kelvin!
 colvars:   # applyBias = on [default]
-colvars:   # updateBias = on [default]
 colvars:   # hideJacobian = off [default]
 colvars:   Jacobian (geometric) forces will be included in reported free energy gradients.
 colvars:   # fullSamples = 10
+colvars:   # minSamples = 5 [default]
 colvars:   # inputPrefix =  [default]
 colvars:   # historyFreq = 0 [default]
 colvars:   # shared = off [default]
+colvars:   # updateBias = on [default]
 colvars:   # maxForce =  [default]
 colvars:   # integrate = on [default]
 colvars:   Finished ABF setup.
@@ -125,16 +128,14 @@ colvars: Collective variables biases initialized, 1 in total.
 colvars: ----------------------------------------------------------------------
 colvars: Collective variables module (re)initialized.
 colvars: ----------------------------------------------------------------------
-colvars: The restart output state file will be "rest.colvars.state".
-colvars: The final output state file will be "test.colvars.state".
-colvars: Opening trajectory file "test.colvars.traj".
+colvars: Current simulation parameters: initial step = 0, integration timestep = 1
+colvars: Updating atomic parameters (masses, charges, etc).
 colvars: Re-initialized atom group for variable "one":0/0. 4 atoms: total mass = 54.028, total charge = -0.72.
 colvars: Re-initialized atom group for variable "one":0/1. 4 atoms: total mass = 54.028, total charge = -0.4.
-colvars:   Prepared sample and gradient buffers at step 0.
+colvars: The final output state file will be "test.colvars.state".
 colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
 colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
-colvars: Saving collective variables state to "rest.colvars.state".
+colvars: Saving collective variables state to "test.colvars.state".
 colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj".
-colvars: Saving collective variables state to "rest.colvars.state".
 colvars: Saving collective variables state to "test.colvars.state".
-colvars: Resetting the Collective Variables module.
+colvars: Saving collective variables state to "test.colvars.state".

lammps/tests/library/000_distance-grid_abf/AutoDiff/test.colvars.state.stripped

@@ -5,13 +5,13 @@ configuration {
 
 colvar {
   name one
-  x  3.20036069218762e+00
+  x       3.2003843562355
 }
 
 abf {
   configuration {
-    step 20
-    name abf1
+step 20
+name abf1
   }
 
 samples
@@ -20,8 +20,8 @@ samples
  0 0 0 0
 
 gradient
- 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 -1.12776581376024e-01 0.00000000000000e+00
- 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
- 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00
+ 0 0 0 0 0 0 -0.10385485654788 0
+ 0 0 0 0 0 0 0 0
+ 0 0 0 0
 }
 

lammps/tests/library/000_distance-grid_abf/AutoDiff/test.colvars.traj

@@ -5,18 +5,18 @@
            3    3.20405779931486e+00  0.00000000000000e+00 
            4    3.20488715757478e+00  0.00000000000000e+00 
            5    3.20630753860023e+00  0.00000000000000e+00 
-           6    3.20822521484009e+00 -6.61626791923408e+00 
-           7    3.21041717634450e+00 -1.27156541500808e+01 
-           8    3.21265147474750e+00 -1.78199101059935e+01 
-           9    3.21470018041825e+00 -2.15984911523814e+01 
-          10    3.21635524474136e+00 -2.38649655797079e+01 
-          11    3.21744175542164e+00 -2.04771085292232e+01 
-          12    3.21789227153408e+00 -1.70074606132688e+01 
-          13    3.21766748124389e+00 -1.36104049871579e+01 
-          14    3.21676135390922e+00 -1.04308363259766e+01 
-          15    3.21520071953742e+00 -7.58252414070870e+00 
-          16    3.21304286766655e+00 -5.14775867496289e+00 
-          17    3.21037093904356e+00 -3.17629358378122e+00 
-          18    3.20728794724505e+00 -1.68764750172374e+00 
-          19    3.20391020908566e+00 -6.75514681148352e-01 
-          20    3.20036069218762e+00 -1.12776581376024e-01 
+           6    3.20822521484009e+00 -6.61157190854750e+00 
+           7    3.21041724907750e+00 -1.27050301328157e+01 
+           8    3.21265178240395e+00 -1.78028628544160e+01 
+           9    3.21470098068744e+00 -2.15753129933168e+01 
+          10    3.21635688320852e+00 -2.38366465889628e+01 
+          11    3.21744464629897e+00 -2.04504671342194e+01 
+          12    3.21789678756116e+00 -1.69832533952049e+01 
+          13    3.21767393972780e+00 -1.35890910889714e+01 
+          14    3.21677000635471e+00 -1.04125653142527e+01 
+          15    3.21521174831218e+00 -7.56715843574337e+00 
+          16    3.21305638863422e+00 -5.13491652636543e+00 
+          17    3.21038700938995e+00 -3.16541313007834e+00 
+          18    3.20730657726304e+00 -1.67806747099744e+00 
+          19    3.20393137633088e+00 -6.66563054214209e-01 
+          20    3.20038435623547e+00 -1.03854856547883e-01 

lammps/tests/library/000_distance-grid_abf/AutoDiff/test.grad

@@ -7,7 +7,7 @@
   1.75000000000000e+00   0.00000000000000e+00
   2.25000000000000e+00   0.00000000000000e+00
   2.75000000000000e+00   0.00000000000000e+00
-  3.25000000000000e+00  -1.12776581375979e-01
+  3.25000000000000e+00  -1.03854856547883e-01
   3.75000000000000e+00   0.00000000000000e+00
   4.25000000000000e+00   0.00000000000000e+00
   4.75000000000000e+00   0.00000000000000e+00

lammps/tests/library/000_distance-grid_abf/AutoDiff/test.pmf

@@ -1,13 +1,13 @@
 # 1
 # -2.50000000000000e-01  5.00000000000000e-01         21  0
 
-  0.00000000000000e+00   5.63882906879895e-02
-  5.00000000000000e-01   5.63882906879895e-02
-  1.00000000000000e+00   5.63882906879895e-02
-  1.50000000000000e+00   5.63882906879895e-02
-  2.00000000000000e+00   5.63882906879895e-02
-  2.50000000000000e+00   5.63882906879895e-02
-  3.00000000000000e+00   5.63882906879895e-02
+  0.00000000000000e+00   5.19274282739416e-02
+  5.00000000000000e-01   5.19274282739416e-02
+  1.00000000000000e+00   5.19274282739416e-02
+  1.50000000000000e+00   5.19274282739416e-02
+  2.00000000000000e+00   5.19274282739416e-02
+  2.50000000000000e+00   5.19274282739416e-02
+  3.00000000000000e+00   5.19274282739416e-02
   3.50000000000000e+00   0.00000000000000e+00
   4.00000000000000e+00   0.00000000000000e+00
   4.50000000000000e+00   0.00000000000000e+00