[gmx-users] Testing gromacs: Thermodynamics integration using tabulated potential
Fahimeh
fahimeh.baftizadeh at gmail.com
Tue May 27 22:48:01 CEST 2014
Sorry the solute is ethanol!
Fahimeh
On May 27, 2014, at 4:44 PM, Fahimeh Baftizadeh <fahimeh.baftizadeh at gmail.com> wrote:
> 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-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.
> > >
> >
> >
> >
> > --
> >
> > ---------------------------------------------------------------------------------------------------
> > "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