# *Help* with Thermodynamic Integration, please...

• Subject: *Help* with Thermodynamic Integration, please...
• Date: Wed, 25 Mar 1992 13:09:09 +1000

```To all practitioners of the black art of Thermodynamic Integration,
(and devotees of 'slow growth' and Free Energy Perturbation) ...
Greetings,
in an effort to generate some discussion about molecular dynamics
methodology for calculating free-energies (and so that I might
increase my scant knowledge), I offer the following data based on
a tutorial on 'relative free energy' in Biosym's Discover manual -
"Performing a Methanol to Ethane Free Energy Calculation".
I have performed the aforementioned gas-phase calculation, in which
one molecule of methanol is alchemically 'warped' into a molecule of
ethane. Biosym's cvff forcefield was used, and the free-energy
change was calculated by Thermodynamic Integration with 6 quadrature
points, dlambda=0.005.
The example in the manual used 10,000 x 1fs (=10ps) equilibration,
followed by 10,000 x 1fs (=10ps) dynamics at 300K and the free-energy
change is reported as 1.45 +/- 0.19 kcal/mol.
After installation of our new version of Discover, I performed this
calculation, giving a free-energy change of -0.98 +/- 0.15 kcal/mol.
I was worried that this difference was due to the change to the new
version of Discover (2.08b) but then realised that the initial velocities
are randomised with different seeds.  Even so, the difference between
these two numbers was intriguing, so I did a few more runs (identical
integration parameters) which are summarised below.
run#  stepsize/fs  equil.steps  dynamics steps      dA   +/-
0      1.0           10,000        10,000        +1.45  0.19 (tutorial)
1      1.0           10,000        10,000        -0.98  0.15
2      1.0           10,000        10,000        +0.67  0.22
.. "obviously" (?!) equilibrium not reached, so
3      1.0           20,000        10,000        +1.24  0.19
4      1.0           20,000        10,000        +1.51  0.18
.. Ah, much better, but let's try to get those std.devs smaller!
5      1.0           20,000        20,000        +0.53  0.17
6      1.0           20,000        20,000        -0.43  0.12
.. Maybe the steps are too big for such small molecules :-)
7      0.5           20,000        20,000        +0.93  0.17
8      0.5           20,000        20,000        +1.75  0.16
I have plotted Temp (and E_Total and E_Potential) vs time for several of
these runs (during both the equilibrium and dynamics phases, for each
of the lambda steps) and am worried by the wide temperature variation.
For example, run 1:
lambda    Temp.  +/-      E_Tot  +/-      E_Pot  +/-
0.05635   299.9  7.0      23.03  1.53     15.88  1.54  <- equm phase
300.0  7.7      22.08  0.03     14.92  0.18  <- dynamics phase
0.20000   299.9  7.4      22.96  0.82     15.81  0.84
299.6 23.7      24.28  0.41     17.14  0.69
0.44365   299.9 10.2      26.81  0.91     19.66  0.96
300.3 36.5      27.21  0.45     20.05  0.88
0.55635   299.9 10.0      25.84  1.28     18.69  1.29
300.0 38.3      26.02  0.21     18.87  0.92
0.75000   299.8 11.0      26.55  1.26     19.40  1.28
300.2 50.2      26.36  0.30     19.21  1.19
0.94365   299.9  8.6      25.37  0.61     18.22  0.63
300.4 49.3      25.44  0.22     18.28  1.16
At this stage I realised that I don't really understand what
is going on here and such powerful tools are dangerous in the hands
of incompetents.
A few questions:
1) What steps must one take to ensure that more than two successive
runs will give agreement within a few +/-'s?  I conclude that agreement
between runs 0,3&4 is fortuitous. Simply making the run longer doesn't
seem to be the answer (see runs 5&6 cf 3&4).
2) What about the temperature (and energy variations) - is this normal?
3) Is the number of quadrature points a problem? If so, how do
you choose the correct number?
4) What about d(lambda)=0.005? Too small, too big?  How do you know?
5) How do you make the +/-'s smaller?
My apologies for the length of the posting, but I think that many
of us would appreciate pointers to pitfalls in this area, by some
of you who are more expert.
* All* comments are very welcome (by e-mail if they are abusive!).
Thanks for listening,
Alan Arnold                          |  e-mail: apa[ AT