*Help* with Thermodynamic Integration, please...



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
 ]ccadfa.cc.adfa.oz.auChem.
  Department,University College  |  voice : +61 6 268 8080
 Australian Defence Force Academy     |  fax   : +61 6 268 8002
 CANBERRA  ACT 2600 Australia         |
 PS my apologies to the administrator/s of this discussion group -\
 for sending a copy to chemistry-request rather than chemistry!