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

[David van der Spoel]

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.


>
> Berk
>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se