**From:** Sadegh Faramarzi Ganjabad (*safaramarziganjabad_at_mix.wvu.edu*)

**Date:** Mon Apr 23 2018 - 11:01:02 CDT

**Next message:**Brian Radak: "Re: Question about "useConstantArea" in puresure control model"**Previous message:**Vermaas, Joshua: "RE: Question about "useConstantArea" in puresure control model"**In reply to:**Brian Radak: "Re: How to interpret Thermodynamic Integration output"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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--
*

**Next message:**Brian Radak: "Re: Question about "useConstantArea" in puresure control model"**Previous message:**Vermaas, Joshua: "RE: Question about "useConstantArea" in puresure control model"**In reply to:**Brian Radak: "Re: How to interpret Thermodynamic Integration output"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

*
This archive was generated by hypermail 2.1.6
: Mon Dec 31 2018 - 23:21:01 CST
*