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

Szilárd Páll pall.szilard at gmail.com
Mon May 5 17:32:25 CEST 2014 

On Mon, May 5, 2014 at 10:45 AM, Mark Abraham <mark.j.abraham at gmail.com>wrote:

> Thanks for the reports. Szilárd and I have also observed similar problems
> in recent months, so I will try to spend some Quality Time with it this
> week.
>

Indeed, just to confirm, I have observed similar issues and got to a
reduced test case that showed (qualitatively) similar time-step dependence
and tau independence.

We have been discussing possible explanations, but for a definitive answer
that explains what exactly is happening here we'll need to take some more
time to think.

Cheers.
--
Szilárd


>
> Mark
> On May 4, 2014 11:48 PM, "Christopher Neale" <chris.neale at alum.utoronto.ca
> >
> wrote:
>
> > 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.
> > --
> > 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.
>

More information about the gromacs.org_gmx-users mailing list