Re: How to interpret Thermodynamic Integration output

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Sun Feb 25 2018 - 17:23:39 CST

Hi Sadegh,

To add to what Brian said, the way I tend to parse the output is in a two-step process. In the first step, you are only looking to extract the non-averaged columns of the output file, which may be a simple call to awk or numpy.loadtxt, depending on your comfort level in bash or python. Taking the non-averaged columns let you have some idea of how noisy your averaging process needs to be, and lets you rigorously compute error estimates for each point. From there, you have the value for dG/dlambda along your lambda coordinate, which you can plot/integrate to get your final answer. I would, however, NOT recommend a simple trapezoidal rule in general for the integration. While it is normally ok for TI, if you ever decide to switch to FEP, I find that a spline interpolation deals better with the spiky end points I see on occasion, so learning how to integrate an interpolated curve is a useful skill to have (see https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.html#scipy.interpolate.B
Spline).

-Josh

On 02/25/2018 11:52 AM, Brian Radak wrote:
Hi Sadegh,

Yes, this is indeed a slightly undocumented update in the latest version (though not much less documented than the old version).

The derivatives are separated by alchemical groups (group 1 is scaled with lambda and group 2 is scaled with 1 - lambda, these are specified in the input PDB with 1 and -1).

Note also that everything in the alchOutFile is an average. The normal output is a local/block average over the last alchOutFreq steps and the AVG columns are cumulative averages. There's really no good reason to have alchOutFreq as small as you have it - it just hurts performance due to frequent I/O.

It is up to the user to decide what they want to do with the output. If you are using simple trapezoidal integration, you would have to scale each column according to the appropriate deltaLambda increment (which you have to compute yourself).

HTH,
BKR

P.S. Note also that the instantaneous TI derivatives are now also sent to stdout according to outputEnergies. You can safely suppress all output to alchOutFile by setting alchOutFreq 0 and just use the stdout if you prefer - I sincerely doubt you would be able to detect a difference in the analysis of any reasonable calculation.

On Fri, Feb 23, 2018 at 7:09 PM, Sadegh Faramarzi Ganjabad <safaramarziganjabad_at_mix.wvu.edu<mailto:safaramarziganjabad_at_mix.wvu.edu>> wrote:
Hello all,

I have done a TI calculation, but can't find any tutorial or document about how to interpret the output. Here is how my output file starts

#TITITLE: TS BOND1 AVGBOND1 ELECT1 AVGELECT1 VDW1 AVGVDW1 BOND2 AVGBOND2 ELECT2 AVGELECT2 VDW2 AVGVDW2
#NEW TI WINDOW: LAMBDA 0
#PARTITION 1 SCALING: BOND 0 VDW 0 ELEC 0
#PARTITION 2 SCALING: BOND 1 VDW 1 ELEC 1
#CONSTANT TEMPERATURE: 310 K
TI: 0 19.7158 19.7158 -31.3316 -31.3316 -6.5855 -6.5855 12.4099 12.4099 -40.2905 -40.2905 21.8477 21.8477
TI: 10 19.4744 19.4964 -31.2402 -31.2486 -6.5801 -6.5806 10.7752 10.9238 -40.0014 -40.0277 18.3626 18.6795
TI: 20 24.9347 22.0861 -29.0597 -30.2062 -6.5748 -6.5778 10.0471 10.5063 -38.9418 -39.5106 12.0912 15.5422
TI: 30 32.2653 25.3697 -25.5063 -28.6901 -6.5840 -6.5798 9.2265 10.0935 -38.3273 -39.1289 10.4395 13.8962
TI: 40 42.3413 29.5091 -20.9719 -26.8077 -6.6005 -6.5848 9.9239 10.0521 -38.7601 -39.0389 13.5623 13.8147
TI: 50 58.5311 35.1997 -15.9560 -24.6799 -6.6079 -6.5894 11.9345 10.4212 -39.2660 -39.0835 12.5075 13.5584
TI: 60 80.3605 42.6031 -12.8759 -22.7448 -6.6090 -6.5926 10.8168 10.4861 -39.5107 -39.1535 6.7353 12.4399
TI: 70 102.0040 50.9694 -14.6704 -21.6076 -6.6234 -6.5969 10.6856 10.5142 -39.9275 -39.2625 5.4385 11.4538
TI: 80 119.1639 59.3885 -20.0657 -21.4172 -6.6514 -6.6037 12.3398 10.7396 -39.4442 -39.2849 8.3762 11.0738
TI: 90 130.3328 67.1846 -23.9356 -21.6940 -6.6710 -6.6111 9.7465 10.6304 -37.6627 -39.1067 11.6755 11.1399
#100 STEPS OF EQUILIBRATION AT LAMBDA 0 COMPLETED
#STARTING COLLECTION OF ENSEMBLE AVERAGE
TI: 100 127.3925 127.3925 -23.4962 -23.4962 -6.6485 -6.6485 9.9449 9.9449 -37.7934 -37.7934 9.3023 9.3023
TI: 110 123.5767 123.9236 -23.1205 -23.1546 -6.6329 -6.6343 10.3050 10.2723 -39.3857 -39.2410 6.6564 6.8969

It seems like BOND1 and BOND2 etc. are for the initial and final values of lambda, respectively. If so, is AVGBOND1+AVGELECT1+AVGVDW1 the total TI energy for each timestep? then if I sum the total energies of all timesteps I will get the total TI for the total transition from lambda=0 to lambda=1?

Thanks in advance,
Sadegh

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2019 - 23:19:43 CST