[gmx-users] Testing gromacs: Thermodynamics integration using tabulated potential
Justin Lemkul
jalemkul at vt.edu
Wed May 28 01:09:32 CEST 2014
On 5/27/14, 4:51 PM, Fahimeh wrote:
> Sorry again, solute is paracetamol and solvent is ethanol and I am using Cgenff force field.
>
In that case, I would consider your nonbonded settings a bit suspect. Those are
not the settings we normally use with CGenFF or CHARMM. I can't speak to the
issues with tables, but I would suggest you revisit these calculations using
normal parameters. Is your outcome "correct" as far as experimental data are
concerned?
-Justin
> 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.
>>
>
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list