The behaviour of Ti estimator on amber output

[xiki-tempula]

In light of the recent PR #304, I think there is a question worth discussing, which is how we expect the TI estimator to behave when used on AMBER output.

For Gromacs, it is easy, the Gromacs will generate the dhdl for all lambdas including lambda=0 and lambda=1, so one just uses the trapezoidal rule. The implementation of the TI estimator is geared toward this kind of input and will give the free energy at lambda=0 and lambda=1 and retiring the result is also easy as one would have the value at lambda=1.

For amber (I don't use AMBER very often so it is likely that I will be wrong), it is more problematic as AMBER doesn't print out dhdl for lambda=0 (maybe for lambda=1 as well). In some sense, it is understandable as the dhdl is ill-defined at both end states. However, this would cause trouble to our Ti estimator, as the input lambda will be for example [0.21132, 0.78867] instead of [0,0, 0.21132, 0.78867, 1.0]. Thus, the output will also be defined at [0.21132, 0.78867] instead of [0,0, 0.21132, 0.78867, 1.0], which means that one cannot really retrieve the difference between lambda 0 and lambda 1.

I guess the Ti estimator could be changed such that if lambda 0 and lambda 1 is not defined, these two lambda will be automatically appended and the result being extrapolated.

I wonder what is the thoughts on this? @DrDomenicoMarson It strikes me that you are quite familiar with AMBER.

[xiki-tempula]

Ok, I guess amber does write out the dhdl at 0 and 1 like all other MD engines. No idea where I got that impression.

[hl2500]

@xiki-tempula There are two MD engines, sander and pmemd (pmemd is more commonly used), in Amber. It's my impression that sander did not allow sampling at 0.0 and 1.0. I couldn't find a reference at this moment, but it was mentioned in the tutorial, under Run TI and the Commonalities and Differences between Pmemd and Sander section.

A lambda value very close to 0.0 or 1.0 could be used, so it might not significantly influence the results from TI().fit(). As another option, gaussian quadrature could be applied and was recommended in Amber manual, page 487.

[xiki-tempula]

@hl2500 Thanks for the explanation.

[DrDomenicoMarson]

Hi, following on @hl2500 comment, I think nowadays nobody is using sander, as it is really slow (moreover now pmemd is available for free now).

Moreover, the more advanced functionalities regarding TI are implemented only in pmemd, and in this implementation is possible to have simulations at lambda 0 and 1, so there should be no issue at all!