[gmx-developers] Conserved energy drifts more with shorter time steps

[Berk Hess]

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.

Berk
>
>
>>
>> Berk
>>
>
>