[gmx-developers] Conserved energy drifts more with shorter time steps
Berk Hess
hess at kth.se
Mon Jun 8 14:12:53 CEST 2015
On 08/06/15 14:05 , David van der Spoel wrote:
> 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.
But v-rescale has less effect on velocities than NH.
Berk
>
> 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
>>>>>>> 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