Sanku M wrote:
Hi,
I am trying to calculate the solvation free energy using thermodynamic
integration(TI) method . I am using gromacs 4.0.7 . But, I am having a
problem in getting accurate average temperature .
The following is what I did:
I found that when doing TI, grompp recommends using 'sd' ( stochastic
dynamics ) as integrator ( provides an warning that for decoupled
system, it is better to use sd in stead of md ).
Also I went through Justin Lemkul's tutorial which also recommends using
'sd' for TI method. So, I tried to use sd integrator.
The tutorial uses BAR for calculation free energy, not TI. Though the protocols
are similar, there is a difference, FYI.
I wanted to run NVT simulation at 300 K.
I also found from manual and also from Justin's website on FEP, sd
integrator implicitly controls temperature and so there is no need to
specify a thermostat. So, I did not specify any thermostat ( I wrote
tcoupl = No ). But, unfortunately, I found that after a long simulation
, the average temperature actually goes to 303 K in stead of 300 K( the
desired temperature). This happened for all Lambda values where the
average temperature turned out to be 303 K. The .mdp file is pasted at
the end of the email and I felt I am using reasonable cutoffs and PME as
electrostatics.
Here is the output from g_energy command to get the temperature.
Statistics over 4000001 steps [ 0.0000 thru 8000.0005 ps ], 1 data sets
All averages are exact over 4000001 steps
Energy Average RMSD Fluct. Drift
Tot-Drift
-------------------------------------------------------------------------------
Temperature 303.333 4.78084 4.78082 5.46937e-06
0.043755
Heat Capacity Cv: 12.4764 J/mol K (factor = 0.00024841)
I tried three further tests:
a) I removed tcoupl = No option . But, it is giving same 303K as average
temperature when using sd.
a) I tried to specify the Nose-Hoover thermostat along with sd. But
still, I found that when using sd, the average temperature is going to
303 K in stead of 300 K.
From the manual description of the sd integrator: "The parameter tcoupl is
ignored." This explains (a) and (b).
c) Finally, I tried to overlook grompp warning and went ahead and used
'md' along with Nose-Hoover thermostat. This time I found the right
average temperature 300 K is being achieved.B
But, I understand that for a decoupled system like here, I need to use a
method like stochastic dynamics.
But , I was wondering why it is reaching 303 K in stead of 300 K when
using sd but 300K is achieved when using md. I looked at the mailing
list where people had issues with sd and temperature control. But, I
could not find a good solution. So, any help on the right protocol will
be really appreciated.
What type of equilibration did you do prior to the data collection? If your
system isn't sampling the desired ensemble, then you shouldn't proceed.
-Justin
Here is my .mdp file.
; RUN CONTROL PARAMETERS
integrator = sd
; Start time and timestep in ps
tinit = 0.0
dt = 0.002
nsteps = 4000000
; For exact run continuation or redoing part of a run
; Part index is updated automatically on checkpointing (keeps files
separate)
simulation_part = 1
init_step = 0
; mode for center of mass motion removal
comm-mode = Linear
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =
; LANGEVIN DYNAMICS OPTIONS
; Friction coefficient (amu/ps) and random seed
bd-fric = 0
ld-seed = 1993
; OUTPUT CONTROL OPTIONS
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 10000
nstvout = 10000
nstfout = 0
; Output frequency for energies to log file and energy file
nstlog = 1000
nstenergy = 10
; Output frequency and precision for xtc file
nstxtcout = 1000
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps =
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 10
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = no
; nblist cut-off
rlist = 1.4
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = pme
rcoulomb-switch = 0
rcoulomb = 1.4
; Relative dielectric constant for the medium and the reaction field
epsilon_r = 1
epsilon_rf = 1
; Method for doing Van der Waals
epsilon_rf = 1
; Method for doing Van der Waals
vdw-type = shift
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Seperate tables between energy group pairs
energygrp_table =
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 6
ewald_rtol = 1e-6
ewald_geometry = 3d
epsilon_surface = 0
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = No
; Groups to couple separately
tc_grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 1.0
ref_t = 300
; Pressure coupling
Pcoupl = no
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0
compressibility = 4.5e-5
ref-p = 1
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = No
; Random seed for Andersen thermostat
andersen_seed = 124821
gen-vel = no
gen-temp = 300
gen-seed = -1
; OPTIONS FOR BONDS
constraints = hbonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
continuation = yes
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Free energy control stuff
free-energy = yes
init-lambda = .05
delta-lambda = 0
sc-alpha = 0.0
sc-power = 0
sc-sigma = 0.0
couple-moltype = Protein
couple-lambda0 = vdw
couple-lambda1 = vdw-q
couple-intramol = no
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
