From chemistry-request@ccl.net Tue Mar 24 23:53:20 1992 Date: Wed, 25 Mar 1992 13:09:09 +1000 From: apa@ccadfa.cc.adfa.oz.AU (Alan P Arnold) Subject: *Help* with Thermodynamic Integration, please... To: chemistry@ccl.net Status: R 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@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!