[gmx-developers] multicomponent lambda free energy calculations.

[Michael Shirts]

I've put together a reworking of the foreign lambda calculation in the
CVS to make it easier to use best practices for free energy
calculations.  This code is caught up to CVS master (as of today).
This can be downloaded for testing with the command:

git clone git://git.gromacs.org/~mrshirts/gromacs_mrs.git

What is described below sounds userful, You should 1) try to see if
you can rewrite your topologies in a way that can use fewer separate
topologies with this new code and 2) see if you can set up tests that
-verify- that what is being done is actually correct, comparing
between the old and new codes.  I've checked for several test cases,
but there may always be some new test case that I missed.

Key points:

1.  Previous functionality is preserved (including couple keywords),
with the exception that instead of using the key word
'foreign_lambda', the key word 'fep-lambda' is used instead.  A bug in
the free energy calculations code 1,4 interactions without LJ terms
was found -- this has also been fixed in the main branch as well.

2. Lambda can now be expressed as a vector of the different
interaction terms; currently 'coul-lambda', 'vdw-lambda',
'bonded-lambda', 'restraint-lambda' and 'mass-lambda' are possible,
but finer/different divisions are possible in the future as needed.

This make it possible to specify a pathway of first turning off the
coulomb term while LJ is on, then turning off the LJ, all with the
same topology, for example, first turning off the charge, then the
Lennard-Jones, while restraints are being switched off.  This can be
cobbled together with messy topology switching, but this should make
it easier.  Additionally, having all states accessible in the same
simulation is required for good sampling in Hamilton replica exchange
or expanded ensemble methods (the last of which may not get
implemented for a bit, since classes start soon).

3. Softcore interactions are only applied to the LJ term, not the
coulomb term, unless specified by 'sc-coul = yes'. I can't think of
any reason one would -want- to do this, but it was the previous
functionality, so I wanted to continue to support it.

The number of components specified will be equal to the number of
components that are entered in "X-lambda" lines.

The free energy portion of the mdp file would look like:
free-energy              = yes
nstdhdl                  = 20
sc-power                 = 1
sc-alpha                 = 0.5
sc-coul                  = no
init_fep_state           = 4
coul-lambda              = 0.0 0.5 1.0 1.0 1.0 1.0 1.0
vdw-lambda               = 0.0 0.0 0.0 0.2 0.4 0.6 1.0

If the A state is fully coupled, and the B state has charged and LJ
terms turned off, then here, the coulomb goes from state A to state B
from states 1 to 3, and then while the coulomb state is B (usually,
turned off), then the LJ terms go to zero from stated 3 to 7.

The current state can be specified either by init_lambda (as was
previously done), setting all of the lambda values to that value, or
by setting init_fep_state, which will set the lambda value to the ith
column - i.e., above, init_fep_state = 4 sets coul-lambda=1.0 and
vdw-lambda=0.2

The default for each component are set by fep-lambda -- if it is not
specified, everything not explicitly specified is left in the A state;
if specified, then everything that is not explicitly set will be the
same as fep-lambda.

Examples:

Example A:
fep-lambda      = 0.0 0.0 0.0
coul-lambda     = 0.0 0.8 1.0
vdw-lambda      = 0.0 0.2 1.0

Coul and vdw terms change, all other terms remain in A state; the same
as if fep-lambda was omitted.

Example B:
fep-lambda        = 0.0 0.0 0.0
coul-lambda       = 0.0 0.8 1.0
vdw-lambda        = 0.0 0.2 0.4
Restraint-lambda  = 0.0 0.4 1.0
Bonded-lambda     = 0.0 0.1 1.0
Mass-lambda       = 0.0 0.0 1.0

Fep-lambda does nothing, since all defined components are specified.

Example C:
fep-lambda        = 0.0 0.2 0.8
coul-lambda       = 0.0 0.8 1.0

Vdw, restraint, bonded, mass all have 0.0 in the first state, 0.2 in
the second state, and 0.8 in the last state, copied from fep-lambda.

This scheme should be sufficient to handle almost all possible
transformations, but can be extended fairly easily as needed.

4. Simplified, non-time averaged distance restraints that work for
free energy calculations; using the same current format for dihedral
restraints (type and label are ignored).  This format may be altered
when merged, however.

5. Free energy enabled dihedral restraints are now possible; simply
specify the B state potential energy parameters.

6. Output format:

@ s0 legend " dH/d\8l\4(coul-lambda)"
@ s1 legend " dH/d\8l\4(vdw-lambda)"
@ s2 legend "\8D\4H \8l\4 (0,0)"
@ s3 legend "\8D\4H \8l\4 (0.5,0)"
@ s4 legend "\8D\4H \8l\4 (1,0)"
@ s5 legend "\8D\4H \8l\4 (1,0.2)"
@ s6 legend "\8D\4H \8l\4 (1,0.6)"
@ s7 legend "\8D\4H \8l\4 (1,1)"
@ s8 legend "pV (kJ/mol)"
0.0000  79.8707 -168.095 0 39.9354 79.8707 65.2871 81.172 112.139 0.55051
0.0400  63.297 -119.502 0 31.6485 63.297 55.0677 75.0248 106.971 0.550642
0.0800  79.9611 -145.203 0 39.9806 79.9611 69.0604 88.2085 120.111 0.549732
0.1200  72.0133 -108.395 0 36.0067 72.0133 65.8608 87.447 119.368 0.550139

The first column is the time in ps.  The next n columns are the value
of the dhdl components (however many were specified in the mdp).

The next N columns are the differences between the current value of
lambda ((0,0), as can be seen by the column that is zero) and the
other values.  These values of Delta E are needed to for the Bennett's
acceptance ratio and FEP free energy ratio methods.  The last column
is the pV energy term that must be added to the energy differences
when constant pressure is used, to convert between Helmholtz and Gibbs
free energy.

7. g_bar is not yet supported -- however, it should be simple, as it
only should require reparsing the dhdl file to handle more than one
dvdl component.  I have a pymbar script to handle the output to
compute free energies, and will assist anyone who wants to try out
this code analyze their simulations with pymbar.

8. I'm working on tutorials for this (and more documentation), so
definitely contact me for any help in the meantime.

~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821