[gmx-users] Langevin sd integrator: wrong average temperature

Andrey Frolov andrey.i.frolov at mail.ru
Sun May 4 16:41:40 CEST 2014


Dear Justin,

Thank you very much for your comments.

Sun, 04 May 2014 09:26:55 -0400 от Justin Lemkul <jalemkul at vt.edu>:
>
>Is the effect dependent on tau_t at all?  How did you decide on a value of 2.0?  For larger "tau_t" the temperature is same: 

integrator,tcoupl,tau_t[ps],dt[ps]        Average   Err.Est.  RMSD  Tot-Drift
-------------------------------------------------------------------------------
sd,no,2,dt=0.002             Temperature    311.5      0.1      6.1      0.6      (K)
sd,no,5,dt=0.002             Temperature    311.9      0.2      6.3     -0.1      (K)

>
>  Is this the correct inverse friction coefficient for supercritical CO2?  Well, I set "tau_t" according to Gromacs manual suggestion: " When [sd] used as a thermostat, an appropriate value for  tau-t  is 2 ps, since this results in a friction that is lower than the internal friction of water, while it is high enough to remove excess heat (unless  cut-off  or  reaction-field  electrostatics is used). NOTE: temperature deviations decay twice as fast as with a Berendsen thermostat with the same  tau-t"

In this case, I aim to use Langevin dynamics as a thermostat, and therefore, to my understanding, the friction coefficient should NOT necessarily correspond to the experimental value. The realistic friction should result itself from the interaction between molecules in the simulation cell. Also, I aim to perform free energy calculations and therefore the dynamics of molecules is not of the highest importance. The important thing is that the sampling of the ensemble should be accurate. As far as I can see, the change of tau_t can really change the dynamics of the molecules, temperature and other properties relaxation times but can not influence much the accuracy of the sampling. Therefore, the tau_t is set arbitrary according to the manual suggestions for water. 

>I 
>have a system I'm working with (RNA in water) and with SD the temperature works 
>out perfectly fine, so likely there is something specific to your system that 
>needs to be considered.
>
>Posting a full .mdp file would also be helpful.  Though you have listed most of 
>the important settings, seeing everything is often useful for diagnostic purposes.
; RUN CONTROL PARAMETERS = 
integrator = sd ; changed by FillMDP.sh
; start time and timestep in ps = 
tinit = 0
dt = 0.002
; 6 ns - enought to sample small molecules (without slow internal degrees of freedom) [Shirts et al.] 
nsteps = 2500000 ; changed by FillMDP.sh
; mode for center of mass motion removal =
; We remove center of mass motion. In periodic boundary conditions, the center of mass motion is spurious; the periodic system is the same in all translational directions.
comm-mode = Linear
; number of steps for center of mass motion removal = 
nstcomm = 1
; OUTPUT CONTROL OPTIONS = 
nstlog = 1000
nstenergy = 100 ; changed by FillMDP.sh
nstxout = 100000
nstvout = 100000
nstfout = 0
nstxtcout = 10000 ; changed by FillMDP.sh
xtc_precision = 100000 ; changed by FillMDP.sh
xtc_grps =
; NEIGHBORSEARCHING PARAMETERS = 
; nblist update frequency = 
nstlist = 10
; ns algorithm (simple or grid) = 
ns_type = grid
; Periodic boundary conditions: xyz or no = 
pbc = xyz
; Neighbor list should be at least 2 A greater than the either rcut or rvdw
; nblist cut-off = 
rlist = 1.4

; OPTIONS FOR ELECTROSTATICS AND VDW: These parameters were all optimized for fast and accurate small molecule calculations.
; See Shirts and Paliwal (2011)
; Method for doing electrostatics = 
coulombtype = PME-Switch
rcoulomb-switch = 1.19999999
rcoulomb = 1.2
; Method for doing Van der Waals = 
vdw-type = Switch
; cut-off lengths = 
rvdw-switch = 1.0
rvdw = 1.2
; Spacing for the PME/PPPM FFT grid = 
fourierspacing = 0.1
; EWALD/PME/PPPM parameters = 
pme_order = 6
ewald_rtol = 1e-06
ewald_geometry = 3d
epsilon_surface = 0
; Apply long range dispersion corrections for Energy and Pressure = 
DispCorr = EnerPres

; Temperature coupling
; Groups to couple separately = 
tcoupl = no ; changed by FillMDP.sh
tc-grps = System
; Time constant (ps) and reference temperature (K) = 
; for equilibration - strong thermostat coupling
tau_t = 2 ; changed by FillMDP.sh
ref_t = 308.15 ; changed by FillMDP.sh
; Pressure coupling = 
;;; use Berendsen ONLY for equilibration run
pcoupl = no 
; Time constant (ps), compressibility (1/bar) and reference P (bar) = 
; for equilibration - strong barostat coupling
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 250.0 ; changed by FillMDP.sh
; We don't strictly need these, because it already has velocities
; that are at the right temperature. But including this is safer.
----------
gen_vel = yes
gen_temp = 308.15 ; changed by FillMDP.sh
gen_seed = 12 ; make sure you set the seed to be able to reproduce the simulation

;;; This contraints section taken from ethanol solvation example on alchemistry.org
;; constrain the hydrogen bonds, allowing longer timesteps.
;; Better to choose a higher lincs order just to be sure that 
;; the constraints are obeyed to high precision; it's not that expensive.
constraints = hbonds
;; Type of constraint algorithm = 
constraint-algorithm = lincs
;; Highest order in the expansion of the constraint coupling matrix = 
lincs-order = 12
 
Andrey



>
>-Justin
>
>-- 
>==================================================
>
>Justin A. Lemkul, Ph.D.
>Ruth L. Kirschstein NRSA Postdoctoral Fellow
>
>Department of Pharmaceutical Sciences
>School of Pharmacy
>Health Sciences Facility II, Room 601
>University of Maryland, Baltimore
>20 Penn St.
>Baltimore, MD 21201
>
>jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>http://mackerell.umaryland.edu/~jalemkul
>
>==================================================



More information about the gromacs.org_gmx-users mailing list