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

Fahimeh fahimeh.baftizadeh at gmail.com
Tue May 27 22:52:04 CEST 2014


Sorry again, solute is paracetamol and solvent is ethanol and I am using Cgenff force field.

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