[gmx-users] Testing gromacs: Thermodynamics integration using tabulated potential

Fahimeh Baftizadeh fahimeh.baftizadeh at gmail.com
Tue May 27 22:44:46 CEST 2014


I will explain everything in details then:

Test1: Performing solvation free energy calculation using 20 lambda points
for decoupling Q. followed by another solvation free energy calculation
using 20 lambda points for decoupling VDW.
The mdp file of the fist calculation is:
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 2500000  ; 5 ns
nstcomm                  = 100
; Output control
nstxout                  = 500
nstvout                  = 500
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxtcout                = 0
xtc-precision            = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.0
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.0
; van der Waals
vdw-type                 = switch
rvdw-switch              = 0.8
rvdw                     = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
optimize_fft             = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 300
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 0.5
compressibility          = 4.5e-05ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda              = 0.00
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0.5
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = A  ; name of moleculetype to decouple
couple-lambda0           = vdw-q      ; only van der Waals interactions
couple-lambda1           = vdw     ; turn off everything, in this case only
vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = all-bonds  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes
; Highest order in the expansion of the constraint coupling matrix
;lincs-order              = 12
lincs-order              = 4

The result of g_bar for this part was:
point  0.000 -  0.050,   DG  5.51 +/-  0.04
point  0.050 -  0.100,   DG  5.09 +/-  0.04
point  0.100 -  0.150,   DG  4.67 +/-  0.02
point  0.150 -  0.200,   DG  4.19 +/-  0.03
point  0.200 -  0.250,   DG  3.76 +/-  0.02
point  0.250 -  0.300,   DG  3.43 +/-  0.04
point  0.300 -  0.350,   DG  3.05 +/-  0.03
point  0.350 -  0.400,   DG  2.65 +/-  0.02
point  0.400 -  0.450,   DG  2.32 +/-  0.01
point  0.450 -  0.500,   DG  2.03 +/-  0.01
point  0.500 -  0.550,   DG  1.76 +/-  0.01
point  0.550 -  0.600,   DG  1.50 +/-  0.02
point  0.600 -  0.650,   DG  1.25 +/-  0.01
point  0.650 -  0.700,   DG  1.06 +/-  0.01
point  0.700 -  0.750,   DG  0.88 +/-  0.01
point  0.750 -  0.800,   DG  0.69 +/-  0.01
point  0.800 -  0.850,   DG  0.52 +/-  0.01
point  0.850 -  0.900,   DG  0.36 +/-  0.01
point  0.900 -  0.950,   DG  0.23 +/-  0.01
point  0.950 -  1.000,   DG  0.10 +/-  0.01

total  0.000 -  1.000,   DG 45.05 +/-  0.09


For decoupling VDW, the charges in the top file was set to zero.  and the
mdp file was the same as above, except these two lines:

couple-lambda0           = vdw      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case
only vdW

The results of the g_bar of this part is:
point  0.000 -  0.050,   DG  2.17 +/-  0.01
point  0.050 -  0.100,   DG  2.25 +/-  0.01
point  0.100 -  0.150,   DG  2.40 +/-  0.01
point  0.150 -  0.200,   DG  2.53 +/-  0.00
point  0.200 -  0.250,   DG  2.60 +/-  0.00
point  0.250 -  0.300,   DG  2.64 +/-  0.00
point  0.300 -  0.350,   DG  2.62 +/-  0.01
point  0.350 -  0.400,   DG  2.56 +/-  0.01
point  0.400 -  0.450,   DG  2.46 +/-  0.00
point  0.450 -  0.500,   DG  2.34 +/-  0.01
point  0.500 -  0.550,   DG  2.17 +/-  0.01
point  0.550 -  0.600,   DG  1.91 +/-  0.01
point  0.600 -  0.650,   DG  1.57 +/-  0.01
point  0.650 -  0.700,   DG  1.10 +/-  0.02
point  0.700 -  0.750,   DG  0.48 +/-  0.03
point  0.750 -  0.800,   DG -0.39 +/-  0.02
point  0.800 -  0.850,   DG -1.07 +/-  0.01
point  0.850 -  0.900,   DG -1.11 +/-  0.03
point  0.900 -  0.950,   DG -0.37 +/-  0.02
point  0.950 -  1.000,   DG  0.88 +/-  0.01

total  0.000 -  1.000,   DG 29.73 +/-  0.06

So I concluded that the total DG=45.05+ 29.73 = 74.78

Test2: decoupling both VDW and Q at the same time. For this simulation, the
charges were the real charges and the mdp file had two modified lines:

couple-lambda0           = vdw-q      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case
only vdW

 The results of this simulation was:
point  0.000 -  0.050,   DG 15.17 +/-  0.16
point  0.050 -  0.100,   DG  6.94 +/-  0.04
point  0.100 -  0.150,   DG  4.91 +/-  0.03
point  0.150 -  0.200,   DG  4.28 +/-  0.02
point  0.200 -  0.250,   DG  4.03 +/-  0.01
point  0.250 -  0.300,   DG  3.88 +/-  0.01
point  0.300 -  0.350,   DG  3.78 +/-  0.02
point  0.350 -  0.400,   DG  3.67 +/-  0.02
point  0.400 -  0.450,   DG  3.53 +/-  0.02
point  0.450 -  0.500,   DG  3.40 +/-  0.01
point  0.500 -  0.550,   DG  3.20 +/-  0.02
point  0.550 -  0.600,   DG  2.94 +/-  0.02
point  0.600 -  0.650,   DG  2.62 +/-  0.02
point  0.650 -  0.700,   DG  2.15 +/-  0.02
point  0.700 -  0.750,   DG  1.60 +/-  0.05
point  0.750 -  0.800,   DG  1.02 +/-  0.06
point  0.800 -  0.850,   DG  0.62 +/-  0.03
point  0.850 -  0.900,   DG  1.05 +/-  0.02
point  0.900 -  0.950,   DG  2.18 +/-  0.02
point  0.950 -  1.000,   DG  3.53 +/-  0.02

total  0.000 -  1.000,   DG 74.51 +/-  0.26

Which is more or less like before!

Test3: I performed another simulation, decoupling Q and VDW at the same
time by 40 lambda points:

point  0.000 -  0.025,   DG  9.75 +/-  0.02
point  0.025 -  0.050,   DG  5.66 +/-  0.01
point  0.050 -  0.075,   DG  3.91 +/-  0.03
point  0.075 -  0.100,   DG  3.05 +/-  0.02
point  0.100 -  0.125,   DG  2.60 +/-  0.02
point  0.125 -  0.150,   DG  2.35 +/-  0.01
point  0.150 -  0.175,   DG  2.18 +/-  0.01
point  0.175 -  0.200,   DG  2.08 +/-  0.01
point  0.200 -  0.225,   DG  2.02 +/-  0.01
point  0.225 -  0.250,   DG  1.98 +/-  0.01
point  0.250 -  0.275,   DG  1.95 +/-  0.00
point  0.275 -  0.300,   DG  1.92 +/-  0.00
point  0.300 -  0.325,   DG  1.89 +/-  0.00
point  0.325 -  0.350,   DG  1.87 +/-  0.00
point  0.350 -  0.375,   DG  1.85 +/-  0.01
point  0.375 -  0.400,   DG  1.82 +/-  0.01
point  0.400 -  0.425,   DG  1.79 +/-  0.01
point  0.425 -  0.450,   DG  1.76 +/-  0.00
point  0.450 -  0.475,   DG  1.72 +/-  0.00
point  0.475 -  0.500,   DG  1.68 +/-  0.01
point  0.500 -  0.525,   DG  1.64 +/-  0.01
point  0.525 -  0.550,   DG  1.58 +/-  0.01
point  0.550 -  0.575,   DG  1.50 +/-  0.01
point  0.575 -  0.600,   DG  1.42 +/-  0.00
point  0.600 -  0.625,   DG  1.34 +/-  0.01
point  0.625 -  0.650,   DG  1.25 +/-  0.01
point  0.650 -  0.675,   DG  1.15 +/-  0.01
point  0.675 -  0.700,   DG  1.04 +/-  0.01
point  0.700 -  0.725,   DG  0.88 +/-  0.01
point  0.725 -  0.750,   DG  0.70 +/-  0.01
point  0.750 -  0.775,   DG  0.54 +/-  0.01
point  0.775 -  0.800,   DG  0.40 +/-  0.02
point  0.800 -  0.825,   DG  0.32 +/-  0.02
point  0.825 -  0.850,   DG  0.33 +/-  0.01
point  0.850 -  0.875,   DG  0.42 +/-  0.01
point  0.875 -  0.900,   DG  0.63 +/-  0.01
point  0.900 -  0.925,   DG  0.91 +/-  0.01
point  0.925 -  0.950,   DG  1.25 +/-  0.00
point  0.950 -  0.975,   DG  1.60 +/-  0.00
point  0.975 -  1.000,   DG  1.93 +/-  0.00

total  0.000 -  1.000,   DG 74.68 +/-  0.06

They are all giving me same results!
What do you think?  am I missing something here?

After all, I would like to receive the same results with tabulated
potential. I can provide the details of that simulations as well. but lets
first converge on this!

Many thanks for your comments
Fahimeh






On Tue, May 27, 2014 at 3:46 PM, Mark Abraham <mark.j.abraham at gmail.com>wrote:

> On Tue, May 27, 2014 at 7:06 PM, Fahimeh Baftizadeh <
> fahimeh.baftizadeh at googlemail.com> wrote:
>
> > Dear Mark,
> >
> > Thanks for your reply.
> > Before this test, I have performed several thermodynamics integration
> > playing with the number of lambda points and also decoupling first
> coulomb
> > and then VDW or decoupling both at the same time. All my tests gave the
> > same results as long as I am not using tabulated form of potential. In
> all
> > of them, I obtained 74kj/mol as solvation free energy!
> >
>
> That does not sound healthy - what did g_bar report for the errors in each
> case?
>
> Now, my question is that, why using tabulated potential makes that much of
> > a difference?
> >
>
> Not by design, but the possibility of a bug is non-negligible. I would pay
> particular attention to the use and implelementation of dispersion
> correction.
>
>
> > Do you mean that decoupling both Q and VDW can be a problem while using
> > tabulated potential?
>
>
> No, decoupling together is a problem in general. But I am not the expert
> here! Anybody who can really help will definitely want to see your .mdp
> files and know what your solute is.
>
>
> > In any case, I will lunch a simulation decoupling them
> > separately and will update you.
> >
>
> That will help establish where the problem lies.
>
> Mark
>
> I really appreciate your help
> > Fahimeh
> >
> >
> >
> > On Tue, May 27, 2014 at 12:20 PM, Mark Abraham <mark.j.abraham at gmail.com
> > >wrote:
> >
> > > Hi,
> > >
> > > Decoupling both q and vdw at the same time is not a good idea. You
> should
> > > observe that at nearly-disappeared lambda, the error bars grow
> > > dramatically.
> > >
> > > Mark
> > >
> > >
> > > On Tue, May 27, 2014 at 5:38 PM, Fahimeh
> > > <fahimeh.baftizadeh at googlemail.com>wrote:
> > >
> > > > Hi everyone,
> > > >
> > > > I am computing solvation free energy of molecule A in solvent B. I
> used
> > > > two different procedure and I thought that the results should be
> > > > comparable, however they are very different.
> > > >
> > > > Procedure 1:
> > > > Thermodynamics integration using 20 lambda points. VDWtype was switch
> > and
> > > > I used PME coulombtype.  I turned off VDW and coulomb interactions
> both
> > > at
> > > > the same time. The free energy of solvation was  74kj/mol
> > > >
> > > > Procedure 2:
> > > > Thermodynamics integration using 20 lambda points. I introduce
> > tabulated
> > > > potential for the interaction between molecule A and molecule B.
> > > therefore
> > > > the vdwtype and columbtype was set to User. I provided two files,
> > called
> > > > table.xvg and table_A_B.xvg which are identical for now as a test.
> > These
> > > > files contain, r, 1/r, 1/r**2, -1/r**6, -6/r**7, 1/r**12, 12/r**13.
> >  The
> > > > other parameters for VDW interactions are untouched. The result of
> this
> > > > simulation however is very different ~ 50kj/mol ...
> > > >
> > > > I am expecting to receive not identical but something comparable ...
> > > could
> > > > you please help me to understand this?
> > > >
> > > > Many thanks in advance
> > > >
> > > > Fahimeh
> > > > --
> > > > Gromacs Users mailing list
> > > >
> > > > * Please search the archive at
> > > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > > > posting!
> > > >
> > > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> > > >
> > > > * For (un)subscribe requests visit
> > > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-usersor
> > > > send a mail to gmx-users-request at gromacs.org.
> > > >
> > > --
> > > Gromacs Users mailing list
> > >
> > > * Please search the archive at
> > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > > posting!
> > >
> > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> > >
> > > * For (un)subscribe requests visit
> > > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > > send a mail to gmx-users-request at gromacs.org.
> > >
> >
> >
> >
> > --
> >
> >
> ---------------------------------------------------------------------------------------------------
> > "If you torture the data long enough, they will confess to anything!"
> John
> > Tukey
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>


More information about the gromacs.org_gmx-users mailing list