Re: How to interpret Thermodynamic Integration output

From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Mon Apr 23 2018 - 11:01:02 CDT

Brian,

I am pinning the two groups of atoms together with a harmonic restraint.

Thanks,
Sadegh Faramarzi,
Research Assistant
West Virginia University, Department of Chemistry
Email:safaramarziganjabad_at_mix.wvu.edu

On Thu, Apr 19, 2018 at 10:15 AM, Brian Radak <brian.radak_at_gmail.com> wrote:

> It's hard to say. One is inclined to believe that the force along the
> coupling coordinate (lambda) is smooth and continuous, but I can't recall a
> proof that says this must be true. It almost certainly depends on the path.
>
> Running simulations at lambda = 0 is a fraught proposition. It has nothing
> to do with TI or FEP per se, but rather the fact that you are asking MD,
> which samples phase space according to the force, to sample without forces
> (since the ligand is completely decoupled). This will almost assuredly
> yield nonsense unless you provide additional restraints - I assume you are
> doing this?
>
>
> On Wed, Apr 18, 2018, 3:58 PM Sadegh Faramarzi Ganjabad <
> safaramarziganjabad_at_mix.wvu.edu> wrote:
>
>> Brian,
>>
>> Thanks for explaining. I ran the extra lambda values between lambda = 0
>> and 0.125. Now the TI values drop more smoothly, but they are generally
>> lower than the original TI. I assume that is because the original
>> simulation wasn't long enough (?).
>>
>> Sadegh Faramarzi,
>> Research Assistant
>> West Virginia University, Department of Chemistry
>> Email:safaramarziganjabad_at_mix.wvu.edu
>>
>> On Fri, Mar 23, 2018 at 11:30 AM, Brian Radak <brian.radak_at_gmail.com>
>> wrote:
>>
>>>
>>>
>>> On Thu, Mar 22, 2018 at 3:13 PM, Sadegh Faramarzi Ganjabad <
>>> safaramarziganjabad_at_mix.wvu.edu> wrote:
>>>
>>>> Brian,
>>>>
>>>> I read that paper, although they discuss the theory rather than the
>>>> procedure to calculate TI energies from the output data. Anyways, I simply
>>>> plotted 'AVGBOND1 + AVGELECT1 + AVGVDW1 - (AVGBOND2 + AVGELECT2 + AVGVDW2)'
>>>> vs. lambda values. The plot looks smooth for each lambda values, except for
>>>> some noises when moving from one lambda to the next one. Here are the
>>>> averages and errors of the cumulative averages for all 16 lambda values
>>>>
>>>> 76.9±2.0 , 24.0±1.5 , 20.6± 1.7, 20.4± 1.8 , 22.3± 1.8 ,
>>>> 21.3± 2.6 , 21.5± 0.9 , 19.3± 1.0 , 16.6± 0.9 , 8.2± 0.9 ,
>>>> -6.0± 1.5 , -21.2± 2.3 , -15.5± 2.6 , -32.2± 1.7 , -21.1±
>>>> 1.8 , -32.5± 1.6
>>>>
>>>> As seen from data, from lambda = 0.0625 to 0.125 there is a sharp
>>>> decrease in TI values (from 76.9 to 24.0). I was wondering
>>>>
>>>> 1. What would be possible reasons for the discontinuity near lambda=
>>>> [0-0.125]
>>>>
>>>> There is generally a large derivative when an atom is almost entirely
>>> decoupled, especially for Lennard-Jones interactions. Using a non-zero
>>> alchVdwShiftCoeff (the default is 5.0 A**2) should help mitigate this.
>>> Otherwise there is no way around this.
>>>
>>> Really the problem here is sampling and has nothing to do with TI. The
>>> derivative has some dependence on lambda, but for very small values lambda
>>> has very little effect on the dynamics (at 0 it's just free diffusion + any
>>> bonded terms you might still have). This is a general problem in alchemical
>>> simulations and there is no general solution as of yet.
>>>
>>>
>>>> 2. How can I resolve it? if I am thinking of going back and running
>>>> some extra lambda values in between 0 and 0.125.
>>>>
>>>> This is recommended.
>>>
>>>
>>>> 3. As seen from the data above, for lambda > 0.5 I get negative values
>>>> for TI. So I can't calculate the logarithm for those values like this plot
>>>> on the manual
>>>>
>>>> http://www.ks.uiuc.edu/Research/namd/2.10/ug/node64.html
>>>>
>>>> should I consider the absolute values of TI?
>>>>
>>>> I have no idea why they chose a log plot other than for easier
>>> visualization. I would just use a logscale ("set logscale y" if you are
>>> using gnuplot).
>>>
>>>
>>>> 4. Regardless of the integration method, the area below the plot
>>>> obtained from the data above will give me DeltaG of the whole process. Is
>>>> that right?
>>>>
>>>> By the Fundamental Theorem of Calculus, the integral of the mean
>>> derivative over an interval will give you the free energy change on that
>>> interval (presumably 0 to 1).
>>>
>>>
>>>> Thanks,
>>>> Sadegh Faramarzi,
>>>> Research Assistant
>>>> West Virginia University, Department of Chemistry
>>>> Email:safaramarziganjabad_at_mix.wvu.edu
>>>>
>>>> On Tue, Feb 27, 2018 at 2:00 PM, Brian Radak <brian.radak_at_gmail.com>
>>>> wrote:
>>>>
>>>>> Unfortunately most of what you conveyed is not at all correct - I
>>>>> strongly recommend you look at the extensive literature on thermodynamic
>>>>> integration. The most recent systematic work I can think of is the
>>>>> comparison paper by Paliwal and Shirts in JCTC, but several other more
>>>>> "seminal" papers exist.
>>>>>
>>>>> You only need 6 numbers from each simulation - the columns beginning
>>>>> with AVG. These are simply the cumulative average over the course of the
>>>>> simulation and the thing that you literally want to integrate. You can also
>>>>> average the block average components or just average the instantaneous
>>>>> values. I doubt that any of the noise between those methods is anywhere
>>>>> near the overall standard error unless your simulation is quite short--0000000000006eb9be056a86271a--

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