[gmx-users] Langevin sd integrator: wrong average temperature
Mark Abraham
mark.j.abraham at gmail.com
Mon May 5 10:45:53 CEST 2014
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.
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.
>
More information about the gromacs.org_gmx-users
mailing list