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

Sun May 4 23:46:02 CEST 2014

No solution here, just wanted to mention that I see this also with simulations of membrane proteins.
With 4.6.1, I get a similar behaviour: actual temperature reported by g_energy exceeds reference temperature.

Average temperatures reported are from 4 repeats of 9 us simulations after discarding 1 us as equilibration; stdev from 4 repeats. Reference temperature was 310 K.

CHARMM + CHARMM lipids + Tips3p:   311.83 +/- 0.01
Amber99SB-ILDN + S-LIPIDS + Tip3p: 311.59 +/- 0.01
OPLS + Berger lipids + Tip4p:      310.63 +/- 0.02

integrator = sd
dt = 0.002
constraints = all-bonds
constraint_algorithm =  lincs
lincs-iter =  1
lincs-order =  6
tc_grps             =  System
tau_t               =  1.0
ld_seed             =  -1
ref_t = 310

Chris.

________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Michael Shirts <mrshirts at gmail.com>
Sent: 04 May 2014 14:01
To: Discussion list for GROMACS users; Andrey Frolov
Subject: Re: [gmx-users] Langevin sd integrator: wrong average temperature

The error estimate is a rough guide, but is not always 100% accurate when
one is making.

I'd consider looking at:

http://pubs.acs.org/doi/abs/10.1021/ct300688p

With accompanying code here:

https://github.com/shirtsgroup/checkensemble

This might help ascertain a bit more how the ensemble that is being
generates deviates from the canonical ensemble.

The time step length does seemsworrisome.  Perhaps you should increase the
number of lincs iterations as well as using higher lincs order?

On Sun, May 4, 2014 at 10:40 AM, Andrey Frolov <andrey.i.frolov at mail.ru>wrote:

>
> 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
> >
> >==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.