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

Szilárd Páll pall.szilard at gmail.com
Tue May 6 12:19:10 CEST 2014 

Andrey,

The above mdp options contain hbonds constraining only which in case
of CO2 don't do anything. Have you really not used all-bonds
constraints in your runs?

Could you please either amend issue 1419 to contain the full and
up-to-date information or alternatively file a new issue? It would
also be helpful to have your latest input files that you obtained the
above results with, preferably mdp,gro, & top.

Cheers,
--
Szilárd


On Mon, May 5, 2014 at 5:32 PM, Szilárd Páll <pall.szilard at gmail.com> wrote:
> 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