From: Patrick Fuchs <[EMAIL PROTECTED]>
Reply-To: Discussion list for GROMACS users <gmx-users@gromacs.org>
To: gmx-users@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@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 [EMAIL PROTECTED]
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
