[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