From: Sadegh Faramarzi Ganjabad (safaramarziganjabad_at_mix.wvu.edu)
Date: Thu Mar 22 2018 - 14:13:03 CDT
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]
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.
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?
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?
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.
>
> If you are using a staged integration scheme (i.e. alchBondLambdaEnd
> and/or alchVdwLambdaEnd != 1 and/or alchElecLambdaStart != 0 - NB This is
> NOT the default) then you have to consider multiple integrals over multiple
> intervals (because the integral formally vanishes even though the integrand
> does not).
>
> If you are using pure linear coupling (NOT the default), then the
> integrand is simply AVGBOND1 + AVGELECT1 + AVGVDW1 - (AVGBOND2 + AVGELECT2
> + AVGVDW2).
>
> If your lambda values are evenly spaced such that you have a fixed
> "deltaLambda" AND you have collected data at your endpoints (i.e. 0 and 1),
> then you can perform trapezoidal integration by taking a weighted sum of
> the integrand (the first and last value get a weight of 1/2 and the rest
> get a weight of 1) and then multiplying through by deltaLambda.
>
> I would be VERY careful with any more complicated integration scheme - the
> efficacy is extremely dependent on the shape of the integrand. Spline
> interpolation can be especially abysmal if you have very few points and may
> not be any better than trapezoidal integration IMO. Furthermore, there is
> no way to assess the systematic error of any integration method by doing
> anything other than adding more lambda values, but this also adds more
> statistical error.
>
> HTH,
> BKR
>
>
> On Tue, Feb 27, 2018 at 1:26 PM, Sadegh Faramarzi Ganjabad <
> safaramarziganjabad_at_mix.wvu.edu> wrote:
>
>> Brian and Josh,
>>
>> Thanks for your replies. As a rough analysis, I plotted the
>> 'non-averaged' columns vs. timesteps. The plots look noisy for some energy
>> components. Also, I multiplied each entry of the columns by 10 (since I
>> sampled every 10 timesteps), and summed them up. I believe that is the same
>> as integration; correct me if I'm wrong. Then I summed the sums of
>> BOND+ELECT+VDW columns for each group (I guess that is the total DeltaG)--0000000000003405470568051b53--
This archive was generated by hypermail 2.1.6 : Tue Dec 31 2019 - 23:19:46 CST