[gmx-developers] Conserved energy drifts more with shorter time steps
Berk Hess
hess at kth.se
Wed Jun 3 11:58:42 CEST 2015
On 03/06/15 11:30 , David van der Spoel wrote:
> On 03/06/15 11:22, Berk Hess wrote:
>> On 03/06/15 11:14 , David van der Spoel 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.
>>> OK, good to know. Did he/she try to just make the update algorithm
>>> double?
>> That's one aspect he should look into.
>>>
>>> By the way, if you grep for invdt in mdlib/* there are quite some real
>>> variables invdt (1/dt) in the constraint code, but that is another
>>> issue.
>> That's not an issue. We don't need the time step accurate to more than
>> single precision. The issue is that we are adding increments that only
>> affects the last few bits of the value your incrementing, so we could
>> use an invdt with even fewer bits than float without affecting the
>> accuracy.
>
> Looking at the update code
> vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
> v[n][d] = vn;
> xprime[n][d] = x[n][d]+vn*dt;
> Makeing vn double precision could be sufficient, I'll give it a go.
How would that help?
The problem is that the increment vn*dt is very small compared to x. You
need more bits in x, not in vn or dt.
Berk
>
>>
>> Berk
>>>
>>>
>>>>
>>>> Berk
>>>>
>>>
>>>
>>
>
>
More information about the gromacs.org_gmx-developers
mailing list