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

[Berk Hess]

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