[gmx-users] Very large fluctuations in dg/dl
Berk Hess
gmx3 at hotmail.com
Tue May 8 10:32:46 CEST 2007
>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
More information about the gromacs.org_gmx-users
mailing list