[gmx-users] Very large fluctuations in dg/dl

Patrick Fuchs Patrick.Fuchs at ebgm.jussieu.fr
Tue May 8 12:45:32 CEST 2007


Hi Berk and Tsjerk,
sounds reassuring since the last version of the GROMOS force field (G53a5 
and a6) has been parametrized to reproduce DGsol using reaction field.
But, you agree Berk that when you want to calculate the solvation free 
energy in the "GROMOS philosophy" (e.g. for reproducing or computing new 
parameters), you have to use RF even if energy conservation might not
be as good as with PME?
Thanks for your answers and your interest in the thread.
Cheers,

Patrick

On Tue, 8 May 2007, Berk Hess wrote:

>
>
>
>> From: "Tsjerk Wassenaar" <tsjerkw at gmail.com>
>> Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>> To: "Discussion list for GROMACS users" <gmx-users at gromacs.org>
>> Subject: Re: [gmx-users] Very large fluctuations in dg/dl
>> Date: Tue, 8 May 2007 06:39:50 +0200
>> 
>> Hi Berk,
>> 
>> Just a question relating to your reply (should be a new thread
>> maybe...). The Gromos force fields were parameterized (for use) with
>> cutoff/RF. PME is of course better for energy conservation, but the
>> inclusion of the interactions beyond the cut-off may give rise to
>> deviant behaviour from the force field. Could you enlighten me on this
>> issue?
>> 
>> Thanks,
>> 
>> Tsjerk
>
> Luckily it does not.
> And PME is not only better for energy conservation, but also captures
> long-range interaction, which for instance leads to more stable proteins.
> I have compared solvation free energies for G53A6 with reaction-field
> and PME and found no differences:
>
> B. Hess and N. van der Vegt
> Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic 
> Comparison of Biomolecular Force Fields and Water Models
> J. Phys. Chem. B 110, 17616 (2006)
> http://dx.doi.org/10.1021/jp0641029
>
> Berk.
>
>> 
>> On 5/7/07, Berk Hess <gmx3 at hotmail.com> wrote:
>>> 
>>> 
>>> 
>>> >From: Patrick Fuchs <Patrick.Fuchs at ebgm.jussieu.fr>
>>> >Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>>> >To: gmx-users at gromacs.org
>>> >Subject: [gmx-users] Very large fluctuations in dg/dl
>>> >Date: Mon, 07 May 2007 16:41:26 +0200
>>> >
>>> >Hi Gromacs users,
>>> >I have a few questions related to solvation free energy calculation via
>>> >thermodynamic integration.
>>> >I'm trying to reproduce some literature data (on e.g. methane,
>>> >methanol...) using the GROMOS G53a6 force field. I followed the tutorial
>>> >of David Mobley (thanks to him BTW), but I used the standard non bonded
>>> >options of the G53a6 force field (instead of OPLS). For each lambda
>>> >value I do a minimization, a 10 ps NVT followed by a 20 ps NPT
>>> >equilibration, and a 1 ns NVT production using the sd integrator. I used
>>> >21 lambda values (0.00, 0.05...1.00).
>>> >Here's my topology file:
>>> >----------------begining of methane.top------------------------
>>> >; topology for a methane molecule
>>> >
>>> >; include GROMOS53a6 force field
>>> >#include "ffG53a6.itp"
>>> >
>>> >;;;;;;; begin methane definition ;;;;;;;
>>> >[ moleculetype ]
>>> >; Name           nrexcl
>>> >METH             3
>>> >
>>> >[ atoms ]
>>> >;nr type resnr residue atom cgnr charge mass    typeB chargeB massB
>>> >  1  CH4  1     METH    C1   0    0.0000 16.0430 DUM   0.0000  16.04300
>>> >;;;;;; end methane definition ;;;;;;;;
>>> >
>>> >; include water topology
>>> >#ifdef FLEX_SPC
>>> >#include "flexspc.itp"
>>> >#else
>>> >#include "spc.itp"
>>> >#endif
>>> >
>>> >[ system ]
>>> >; name
>>> >1 methane molecule in water
>>> >
>>> >[ molecules ]
>>> >; name  number
>>> >METH    1
>>> >SOL     893
>>> >-----------------end of methane.top------------------------
>>> >
>>> >And here is my mdp file for lambda=0:
>>> >---------------begining of prod.mdp---------------------
>>> >title               = production NVT methane/water
>>> >cpp                 = /lib/cpp
>>> >; OPTIONS FOR BOND CONSTRAINTS
>>> >constraints         = all-bonds
>>> >; RUN CONTROL PARAMETERS
>>> >integrator          = sd
>>> >tinit               = 0
>>> >dt                  = 0.002
>>> >nsteps              = 500000 ; 1000 ps
>>> >; NUMBER OF STEPS FOR CENTER OF MASS MOTION REMOVAL
>>> >nstcomm             = 100
>>> >; OUTPUT CONTROL OPTIONS
>>> >nstxout                  = 500
>>> >nstvout                  = 500
>>> >nstfout                  = 0
>>> >nstlog                   = 500
>>> >nstenergy                = 100
>>> >nstxtcout                = 5000
>>> >xtc-precision            = 1000
>>> >; NON BONDED STUFF
>>> >ns_type             =  grid
>>> >nstlist                = 5
>>> >rlist                  = 0.8
>>> >coulombtype            = generalized-reaction-field
>>> >rcoulomb               = 1.4
>>> >rvdw                   = 1.4
>>> >epsilon_rf             = 54.0
>>> >;OPTIONS FOR TEMPERATURE COUPLING
>>> >tc_grps                  = system
>>> >tau_t                    = 0.1 ; inverse langevin friction cst
>>> >ref_t                    = 300
>>> >;OPTIONS FOR PRESSURE COUPLING
>>> >Pcoupl                   = no
>>> >tau_p                    = 0.5
>>> >compressibility          = 4.5e-5
>>> >ref_p                    = 1.0
>>> >; FREE ENERGY CONTROL STUFF
>>> >free_energy              = yes
>>> >init_lambda              = 0.00
>>> >delta_lambda             = 0
>>> >sc_alpha                 = 0.5
>>> >sc-power                 = 1.0
>>> >sc-sigma                 = 0.3
>>> >; VELOCITY GENERATION
>>> >gen_vel                  = yes
>>> >gen_temp                 = 300
>>> >gen_seed                 = -1
>>> >-----------------end of prod.mdp------------------------
>>> >I find a reasonable DeltaGsol value of 8.6 kJ/mol for methane (compared
>>> >to 8.7 in Geerke & van Gunsteren, ChemPhysChem 2006, 7, 671 – 678) but
>>> >I get really huge fluctuations in the values of dg/dl:
>>> >lambda=0.00: 5.0 +/- 10.8 (mean +/- standard deviation)
>>> >lambda=0.05: 4.3 +/- 11.2
>>> >...
>>> >lambda=1.00: -0.3 +/- 4.0
>>> >Furthermore, each of these mean value is very slow at converging (1 ns
>>> >seems a minimum for certain lambda values...).
>>> >I can't get reasonable fluctuations even if I sample more. In addition,
>>> >there are very frequent warnings in the log file such as:
>>> >----
>>> >Large VCM(group rest):      0.01363,      0.00818,      0.01147, ekin-cm:
>>> >3.09490e+00
>>> >----
>>> >Here are my questions:
>>> >1) Has someone an idea of what could be the cause of these [very] large
>>> >fluctuations? Does it come from my setup, or is this a normal behavior?
>>> >2) Are these 'Large VCM(group rest)' warnings related to the use of sd
>>> >integrator (when I switch to md integrator, I no longer get these 
>>> warnings)
>>> >?
>>> >Thanks for your answer,
>>> >
>>> >Patrick
>>> 
>>> Large fluctuations are very common, depending on your system
>>> and how your choose your intermediate states.
>>> 
>>> But sampling would improve a lot if you use a larger tau_t, for instance
>>> 1 ps, this will also solve your large vcm warnings.
>>> Note that with reaction-field your will have bad energy conservation.
>>> With a tau_t of 1 ps your temperature might go up somewhat.
>>> Using PME is always better.
>>> 
>>> Berk.
>>> 
>>> _________________________________________________________________
>>> Watch free concerts with Pink, Rod Stewart, Oasis and more. Visit MSN In
>>> Concert today. http://music.msn.com/presents?icid=ncmsnpresentstagline
>>> 
>>> _______________________________________________
>>> gmx-users mailing list    gmx-users at gromacs.org
>>> http://www.gromacs.org/mailman/listinfo/gmx-users
>>> Please search the archive at http://www.gromacs.org/search before posting!
>>> Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-users-request at gromacs.org.
>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>>> 
>> 
>> 
>> --
>> Tsjerk A. Wassenaar, Ph.D.
>> Junior UD (post-doc)
>> Biomolecular NMR, Bijvoet Center
>> Utrecht University
>> Padualaan 8
>> 3584 CH Utrecht
>> The Netherlands
>> P: +31-30-2539931
>> F: +31-30-2537623
>> _______________________________________________
>> gmx-users mailing list    gmx-users at gromacs.org
>> http://www.gromacs.org/mailman/listinfo/gmx-users
>> Please search the archive at http://www.gromacs.org/search before posting!
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
> _________________________________________________________________
> Live Search, for accurate results! http://www.live.nl
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface 
> or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>


More information about the gromacs.org_gmx-users mailing list