[gmx-developers] Conserved energy drifts more with shorter time steps
Berk Hess
hess at kth.se
Wed Jun 3 14:53:18 CEST 2015
On 03/06/15 14:31 , David van der Spoel wrote:
> On 03/06/15 11:58, Berk Hess wrote:
>> 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.
>
> Problem is that we compute the new velocity in double (since dt is
> double), but then round it to float in vn, then we mulitply vn by dt
> again (a double operation). So if we remove the rounding in the
> intermediate step it may help. We're testing this now.
That completely negligible.
Many bits of vn*dt are lost after addition to x. The rounding that is
relevant is in de addition to x, that completely drown any rounding in
vn*dt.
But I now realize that this is could affect performance, we should use a
(local) real dt variable.
Berk
>
>
>>
>> Berk
>>>
>>>>
>>>> Berk
>>>>>
>>>>>
>>>>>>
>>>>>> Berk
>>>>>>
>>>>>
>>>>>
>>>>
>>>
>>>
>>
>
>
More information about the gromacs.org_gmx-developers
mailing list