[gmx-developers] Conserved energy drifts more with shorter time steps
Shirts, Michael R. (mrs5pt)
mrs5pt at eservices.virginia.edu
Wed Jun 3 14:07:09 CEST 2015
One thing to note is that the Lagranian for NH + any verlet scheme gives
an algorithm that is not sympletic. So the nice absence of long-term drift
one sees in NVE does not occur in the conserved quantity NH; there's no
shadow Hamiltonian that the the discretized system follows, so it's only
locally conserved.
Now, drift should at least get smaller as the timestep gets smaller up to
rounding errors. Because NH has some ringing, I would expect that the
drift would behave more consistently for larger systems.
FWIW, I always recall seeing rounding error in single precision around
0.25 when was testing integrators, though you're looking at it more
carefully here.
Best,
~~~~~~~~~~~~
Michael Shirts
Associate Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434) 243-1821
On 6/3/15, 6:36 AM, "Berk Hess" <hess at kth.se> wrote:
>On 03/06/15 12:14 , Berk Hess wrote:
>> On 03/06/15 10:36 , Berk Hess wrote:
>>> On 03/06/15 10:29 , David van der Spoel wrote:
>>>> On 03/06/15 10:02, Berk Hess wrote:
>>>>> On 03/06/15 9:57 , David van der Spoel wrote:
>>>>>> On 02/06/15 15:24, Berk Hess wrote:
>>>>>>> On 2015-06-02 15:20, Anton Feenstra wrote:
>>>>>>>> On 02-06-15 15:14, Berk Hess wrote:
>>>>>>>>> On 2015-06-02 15:08, David van der Spoel wrote:
>>>>>>>>>> On 02/06/15 14:30, Berk Hess wrote:
>>>>>>>>>>> On 2015-06-02 14:26, David van der Spoel wrote:
>>>>>>>>>>>> On 02/06/15 14:05, Berk Hess wrote:
>>>>>>>>>>>>> I assume David is looking at the conserved energy quantity.
>>>>>>>>>>>> Indeed, averages and drift over 900 ps, 1000 molecules:
>>>>>>>>>>>>
>>>>>>>>>>>> Energy Average Err.Est. RMSD Tot-Drift
>>>>>>>>>>>>
>>>>>>>>>>>>---------------------------------------------------------------
>>>>>>>>>>>>----------------
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> Total Energy 20460.7 2.5 778.213 12.4265 (kJ/mol)
>>>>>>>>>>>> Conserved En. -4384.7 380 778.283 2665.29 (kJ/mol)
>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> But I just noticed that Verlet - buffer - tolerance is not
>>>>>>>>>>>>> mentioned,
>>>>>>>>>>>>> so I assume it is default. The drift estimate is an
>>>>>>>>>>>>> overestimate
>>>>>>>>>>>>> which is tighter for shorter pair list lifetimes. I guess
>>>>>>>>>>>>> that is
>>>>>>>>>>>>> the
>>>>>>>>>>>>> explanation. Set it very low, so you are sure the energy
>>>>>>>>>>>>> drift is
>>>>>>>>>>>>> not
>>>>>>>>>>>>> affected by the (too small) buffer.
>>>>>>>>>>>>>
>>>>>>>>>>>> This may be the default.
>>>>>>>>>>>> verlet-buffer-tolerance = 0.005
>>>>>>>>>>>> should we set it to 1e-5 or something like that?
>>>>>>>>>>> It depends on how low drift you want to be able to measure. You
>>>>>>>>>>> need to
>>>>>>>>>>> set it lower than the drift you want to be able to measure.
>>>>>>>>>> Ok, will try. What about the time step dependence? Any clue?
>>>>>>>>> I explained that before, see 4 steps above.
>>>>>>>>>
>>>>>>>>> Berk
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Berk
>>>>>>>>>>>>
>>>>>>>>>>>>> Berk
>>>>>>>>>>>>>
>>>>>>>>>>>>> On Jun 2, 2015 1:45 PM, "Shirts, Michael R. (mrs5pt)"
>>>>>>>>>>>>> <mrs5pt at eservices.virginia.edu> wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> How do we run more NVE than with the settings below?
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> tcoupl = Nose-Hoover
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Would make it not NVE.
>>>>>>>>>>>> Turning tcoupl off would not keep the temperature constant of
>>>>>>>>>>>> course,
>>>>>>>>>>>> seeing that there is a drift in the conserved energy this
>>>>>>>>>>>>means
>>>>>>>>>>>> that
>>>>>>>>>>>> the temperature will change.
>>>>>>>>>>>>
>>>>>>>>>>>> However, the main question was why does "conserved-energy
>>>>>>>>>>>> conservation" get worse when the time step is reduced?
>>>>>>>>>>>> Integration
>>>>>>>>>>>> should presumably become more accurate, or not? Or does one
>>>>>>>>>>>> need
>>>>>>>>>>>> double precision for short time steps?
>>>>>>>>
>>>>>>>>
>>>>>>>> Yes, at least that is what I saw for my and Berk's 1999 JCC paper.
>>>>>>>> At single precision, drift would not go down below a certain time
>>>>>>>> step. I have no idea anymore which timestep precisely, and I'm
>>>>>>>>sure
>>>>>>>> this would be system specific.
>>>>>>>>
>>>>>>>> I also recall drift indeed going up at (extreme) short
>>>>>>>> timesteps. My
>>>>>>>> intuition for that at the time was that arithmetic errors would
>>>>>>>> accumulate faster for calculating many very small delta's.
>>>>>>> In that case it was due to the constraint algorithms. The
>>>>>>> velocity was
>>>>>>> corrected with the coordinate difference divided by the time
>>>>>>> step. Thus
>>>>>>> bit noise is amplified inversely proportional with the time step. I
>>>>>>> fixed this some years ago for LINCS and SHAKE, but SETTLE still
>>>>>>> has this
>>>>>>> property.
>>>>>>>
>>>>>>> David is not using constraints. So I think the issue here is, as I
>>>>>>> mailed before, that the pairlist buffer estimate is more less
>>>>>>> accurate,
>>>>>>> more of an overestimate, for longer pairlist lifetimes (larger time
>>>>>>> steps with fixed nslist), which leads to less drift.
>>>>>>>
>>>>>>> Berk
>>>>>> Unfortunately setting the verlet-buffer to 1e-6 (instead of default
>>>>>> 0.005) does not help. We get the following drift:
>>>>>>
>>>>>> Drift:
>>>>>> Dt Single(5e-3) Double(5-e3) Single(1e-6)
>>>>>> 1.0 240 166 191
>>>>>> 0.5 46 -12 68
>>>>>> 0.2 170 -7 137
>>>>>> 0.1 328 -6 353
>>>>>> 0.005 564 -4 573
>>>>>>
>>>>>> Seems this is still a precision issue, or not?
>>>>>>
>>>>> Yes, this still seems a precision issue.
>>>>> But from this I can't see how large the drift is. What is the energy
>>>>> drift per atom per ns?
>>>>> The drift could be so small that you are affected by bit noise in the
>>>>> coordinate update.
>>>> 14 atoms, 1000 molecules, 900 ps, so divided by 12600, that means
>>>> that for 0.05 fs the drift is ~0.04 kJ/mol/atom/ns
>>>>
>>>> But why would the 0.5 fs be 10 times more accurate than the 0.05 fs?
>>>> Would that be because one does many different additions instead of
>>>> one big one? But why would that not cancel out?
>>>> We see this for 6 different systems.
>>>>
>>> My guess is that, but that's just a guess:
>>> x_i+1 = x_i + dt*v_i+1/2
>>> results in incrementing few bits of x_i+1. Each increment gives a bit
>>> rounding error. 10x smaller time step gives 10x more bit rouding. But
>>> I think that should lead to a diffusive errors, so sqrt(10) times
>>> larger drift.
>>> I have a master student working on investigating exactly this, but
>>> progress is slow.
>>>
>>> Berk
>>>
>> Are you now running NVE, or did you check that NVE gives the same
>> drift as with NH?
>I also only noticed now that you are running NH with a tau_t of 0.5 ps,
>this is too short.
>I did a quick check on a water box with 1 fs time step and NH with
>tau_t=0.5 causes a relatively large energy drift.
>
>Berk
>
>
>--
>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