[gmx-developers] Conserved energy drifts more with shorter time steps
David van der Spoel
spoel at xray.bmc.uu.se
Wed Jun 3 11:30:03 CEST 2015
On 03/06/15 11:22, Berk Hess wrote:
> 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.
Looking at the update code
vn = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
v[n][d] = vn;
xprime[n][d] = x[n][d]+vn*dt;
Makeing vn double precision could be sufficient, I'll give it a go.
>
> Berk
>>
>>
>>>
>>> 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
More information about the gromacs.org_gmx-developers
mailing list