# [gmx-developers] Conserved energy drifts more with shorter time steps

David van der Spoel spoel at xray.bmc.uu.se
Mon Jun 8 14:05:40 CEST 2015

```On 08/06/15 07:47, Berk Hess wrote:
> On 06/08/2015 02:09 AM, Shirts, Michael R. (mrs5pt) wrote:
>>> We need good energy conservation in order to compute heat capacities
>> >from fluctuations.
>>
>> But the important thing is not the energy conservation, but how closely
>> the simulation matches the Boltzmann ensemble.  For example, Langevin
>> dynamics does not conserve any energy, but gives the Boltzmann
>> distribution exactly.  Macroscopically, the heat capacity is defined in
>> terms of d<U>/dT, and the fluctuation formula can be derived from this
>> derivative definition if the energy distribution is Boltzmann
>> distribution
> I wanted to write the same thing.
> I was so focused on the (for me very interesting) topic of energy
> conservation that I only realized this after posting my last reply. A
> run with NH, v-rescale or SD with "moderately" good energy conservation
> should suffice.
We use NH but can not use any stochastic method since we analyze
velocity time correlations as well.

We will look into the ensemble checks though.
>
> Berk
>>
>> This is where I plug checkensemble,
>>
>> https://github.com/shirtsgroup/checkensemble
>>
>> Described here:
>>
>> http://dx.doi.org/10.1021/ct300688p
>>
>>
>> Which shows how one can perform tests to check whether a given set of
>> simulation conditions gives the Boltzmann distribution in the output
>> (also
>> works for NVT, NPT, etc).
>>
>> Also, for good estimation of the heat capacities:
>>
>> https://github.com/choderalab/pymbar/examples/heat-capacity
>>
>>
>> has an example of heat capacity calculated via MBAR in several different
>> ways - by fluctuation, by finite difference by perturbed expectation, and
>> by second derivative of the free energy.
>>
>> Best,
>> ~~~~~~~~~~~~
>> Michael Shirts
>> Associate Professor
>> Department of Chemical Engineering
>> University of Virginia
>> michael.shirts at virginia.edu
>> (434) 243-1821
>>
>>
>>
>> On 6/7/15, 5:16 PM, "David van der Spoel" <spoel at xray.bmc.uu.se> wrote:
>>
>>> On 07/06/15 23:06, Berk Hess wrote:
>>>> On 06/07/2015 11:00 PM, David van der Spoel wrote:
>>>>> On 07/06/15 21:51, Peter Ahlström wrote:
>>>>>> Hello,
>>>>>> I've been following the discussion for a while and came just to think
>>>>>> What about the PME - isn't that messing up the energy conservation?
>>>>>> I believe that that's also a non-conservative force -or am I
>>>>>> completely
>>>>>> wrong?
>>>>> It computes just the Coulomb energy (or Lennard Jones dispersion) in a
>>>>> complex way. We tried to turn off the Lennard Jones PME (using a
>>>>> cutoff) and that leads to poorer energy conservation. One can in fact
>>>>> control the error in the algorithm to be very small, so I think it is
>>>>> not a source of error. The known disadvantage of PME is of course that
>>>>> it enhances periodicity.
>>>> PME conserves energy perfectly (apart from the cut-off effect which you
>>>> control with ewald_rtol), so that's not an issue.
>>>>
>>>> Thermostats are tricky for energy conservation, since they mess with
>>>> the
>>>> velocities. NH with a short period is bad, since the fast oscillations
>>>> need a short time step for accuracy. Any thermostat with a long period
>>>> can cause issues, because the velocity scaling factor only differs from
>>>> 1 by the last few bits, thereby causing systematic rounding errors.
>>>> Using LF or VV does not influence these issues much.
>>>> If I need a thermostat and reasonably good energy conservation I use
>>>> v-rescale with a not too short period and a longer nsttcouple to avoid
>>>> systematic bit rounding.
>>>>
>>>> But why do you want "perfect" energy conservation with NVT? You usually
>>>> want that with NVE.
>>> We need good energy conservation in order to compute heat capacities
>> >from fluctuations. Of course one can do tricks like subtracting the
>>> drift or doing simulations at multiple temperatures, but that is
>>>
>>> It seems that double precision works much better than single. If your
>>> observation about the problem being x_n+1 = x_n + v_n dt is correct it
>>> means we may need to store x in double to be able to do this.
>>>
>>> Or would a unit change to pm and fs help? It seems too simple to be
>>> true.
>>>
>>> We need NVT because we want to compare our results to experimental data
>>> at 298.15K, not at 298.15 +/- 20K.
>>>
>>>> Berk
>>>>
>>>>>> Cheers,
>>>>>> Peter
>>>>>>
>>>>>>
>>>>>>
>>>>>> Dr Peter Ahlström
>>>>>> University of Borås
>>>>>> SE-501 90 Borås
>>>>>> Sweden
>>>>>> +46 33 435 4675
>>>>>>
>>>>>>   >>> David van der Spoel <spoel at xray.bmc.uu.se> 06/07/15 5:58 PM >>>
>>>>>> On 02/06/15 13:45, Shirts, Michael R. (mrs5pt) wrote:
>>>>>>   >> How do we run more NVE than with the settings below?
>>>>>>   >
>>>>>>   > tcoupl = Nose-Hoover
>>>>>>   >
>>>>>>   >
>>>>>>   > Would make it not NVE.
>>>>>>   >
>>>>>>   > Do you mean the conserved quantity isn't conserved? Energy
>>>>>> certainly
>>>>>>   > won't be conserved with Nose-Hoover.
>>>>>> Disturbing altogether.
>>>>>>
>>>>>> Is there any reason that md-vv would be better than the md integrator
>>>>>> for energy conservation?
>>>>>>
>>>>>>   >
>>>>>>   > Best,
>>>>>>   > ~~~~~~~~~~~~
>>>>>>   > Michael Shirts
>>>>>>   > Associate Professor
>>>>>>   > Department of Chemical Engineering
>>>>>>   > University of Virginia
>>>>>>   > michael.shirts at virginia.edu
>>>>>>   > (434) 243-1821
>>>>>>   >
>>>>>>   >
>>>>>>   >
>>>>>>   > On 6/2/15, 7:22 AM, "David van der Spoel" <spoel at xray.bmc.uu.se>
>>>>>> wrote:
>>>>>>   >
>>>>>>   >> On 02/06/15 12:47, Berk Hess wrote:
>>>>>>   >>> Hi,
>>>>>>   >>>
>>>>>>   >>> So this is not using SETTLE or LINCS?
>>>>>>   >>> With SETTLE this is a known issue.
>>>>>>   >> No this is flexible organic molecules at contant volume.
>>>>>>   >>
>>>>>>   >>>
>>>>>>   >>> Did you try running NVE?
>>>>>>   >> How do we run more NVE than with the settings below?
>>>>>>   >>
>>>>>>   >>
>>>>>>   >>>
>>>>>>   >>> Cheers,
>>>>>>   >>>
>>>>>>   >>> Berk
>>>>>>   >>>
>>>>>>   >>> On 2015-06-02 12:43, David van der Spoel wrote:
>>>>>>   >>>> Hi,
>>>>>>   >>>>
>>>>>>   >>>> this is in between a user and a developer query. We find that
>>>>>> for
>>>>>>   >>>> flexible liquid simulations using gromacs 5.0.4/single
>>>>>> precision the
>>>>>>   >>>> energy conservation gets WORSE with decreasing time step. A
>>>>>> plot
>>>>>>   >>>> showing this is in
>>>>>>   >>>> http://folding.bmc.uu.se/images/Econserved-timestep.xvg or
>>>>>>   >>>> http://folding.bmc.uu.se/images/Econserved-timestep.pdf
>>>>>>   >>>>
>>>>>>   >>>> MDP settings
>>>>>>   >>>> coulombtype = PME
>>>>>>   >>>> coulomb-modifier = Potential-shift
>>>>>>   >>>> rcoulomb-switch = 0
>>>>>>   >>>> rcoulomb = 1.1
>>>>>>   >>>> epsilon-r = 1
>>>>>>   >>>> epsilon-rf = inf
>>>>>>   >>>> vdw-type = PME
>>>>>>   >>>> vdw-modifier = Potential-shift
>>>>>>   >>>> rvdw-switch = 0
>>>>>>   >>>> rvdw = 1.1
>>>>>>   >>>> tcoupl = Nose-Hoover
>>>>>>   >>>> nsttcouple = 10
>>>>>>   >>>> nh-chain-length = 1
>>>>>>   >>>> pcoupl = No
>>>>>>   >>>> ref-t: 298.15
>>>>>>   >>>> tau-t: 0.5
>>>>>>   >>>>
>>>>>>   >>>> Any clues whether we are doing something wrong? Or is there a
>>>>>> bug?
>>>>>>   >>>>
>>>>>>   >>>> Cheers,
>>>>>>   >>>
>>>>>>   >>
>>>>>>   >>
>>>>>>   >> --
>>>>>>   >> David van der Spoel, Ph.D., Professor of Biology
>>>>>>   >> Dept. of Cell & Molec. Biol., Uppsala University.
>>>>>>   >> Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>>>>>>   >> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>>>>>   >> --
>>>>>>   >> Gromacs Developers mailing list
>>>>>>   >>
>>>>>>   >> * Please search the archive at
>>>>>>   >> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers
>>>>>>
>>>>>>   >> or send a mail to gmx-developers-request at gromacs.org.
>>>>>>   >
>>>>>>
>>>>>>
>>>>>> --
>>>>>> David van der Spoel, Ph.D., Professor of Biology
>>>>>> Dept. of Cell & Molec. Biol., Uppsala University.
>>>>>> Box 596, 75124 Uppsala, Sweden. Phone:    +46184714205.
>>>>>> spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>>>>> --
>>>>>> Gromacs Developers mailing list
>>>>>>
>>>>>> * Please search the archive at
>>>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers
>>>>>>
>>>>>> or send a mail to gmx-developers-request at gromacs.org.
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>
>>> --
>>> David van der Spoel, Ph.D., Professor of Biology
>>> Dept. of Cell & Molec. Biol., Uppsala University.
>>> Box 596, 75124 Uppsala, Sweden. Phone:    +46184714205.
>>> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
>>> --
>>> Gromacs Developers mailing list
>>>
>>> * Please search the archive at
>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_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-developers
>>> or send a mail to gmx-developers-request at gromacs.org.
>

--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
```