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

[Berk Hess]

>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