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

[Berk Hess]

On 03/06/15 12:14 , Berk Hess 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.
>>
>> Berk
>>
> Are you now running NVE, or did you check that NVE gives the same 
> drift as with NH?
I also only noticed now that you are running NH with a tau_t of 0.5 ps, 
this is too short.
I did a quick check on a water box with 1 fs time step and NH with 
tau_t=0.5 causes a relatively large energy drift.

Berk