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

[Berk Hess]

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?

Berk