[gmx-developers] Conserved energy drifts more with shorter time steps
Shirts, Michael R. (mrs5pt)
mrs5pt at eservices.virginia.edu
Mon Jun 8 02:09:22 CEST 2015
> 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
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