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

Berk Hess hess at kth.se
Mon Jun 8 07:47:17 CEST 2015


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.

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
>>>>> about another problem:
>>>>> 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
>> sticking your head in the sand.
>>
>> 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.



More information about the gromacs.org_gmx-developers mailing list