# How error scales with numerical precision in molecular dynamics?

In terms of time-step, numerical error in molecular dynamics scales with square, i.e. $$error \approx dt^2$$. But how it look for numerical precision ? E.g. how much bigger will be numerical error when I use half precision, instead single and single instead double ?

• The method is Verlet? Which variant, velocity, leapfrog, something else? Did you de-singularize the potential? Is the step size small enough to have a sufficient sampling density during a collision or fly-by? Why would you prefer single or half precision? Is double precision readily available? May 16, 2021 at 10:24
• Which variant, velocity, leapfrog, something else - Leapfrog with BAOAB integration scheme. "Why would you prefer single or half precision" - speed. Nevertheless, more general answer would be preferable. May 16, 2021 at 10:45
• Before you consider speed can you answer how many significant figures you actually need to do the science you want to do? May 16, 2021 at 12:02

## 1 Answer

In general in a correctly implemented fixed-step ODE solver method you have 3 sources for numerical errors:

1. the theoretical method truncation error,
2. the floating point error from evaluating the ODE function and composing the method step and
3. the error from accumulating the single updates of size $$O(h)$$ to the integration result of size $$O(1)$$.

So the total error in each step is something like $$C_1h^{p+1}+C_2\mu_{eval}h+C_2\mu_{acc}$$ where $$\mu_{eval}$$ is the machine constant for the number type used for evaluation and $$\mu_{acc}$$ the one for the number type used in the accumulation of the steps. The global error, in a first approach, results from multiplying with $$N=T/h$$, giving $$C_1h^pT+C_2\mu_{eval}T+C_2\frac{\mu_{acc}T}{h}$$ So ideally (in an academic situation where variables and derivatives have a scale close to $$1$$) you would like to have $$h^p\ge\mu_{eval}\ge\frac{\mu_{acc}}{h}$$. For $$p=2$$ using single precision in evaluation, $$\mu_{eval}=10^{-8}$$, and double precision for the accumulation, $$\mu_{acc}=10^{-15}$$, this would give errors according to the second order of the method down to $$h=10^{-4}$$. Another variant to achieve increased precision in the accumulation without a different data type is Kahan or compensated summation.

• An alternative to IEEE-754 double-precision (binary64) accumulation seems to be Q24.40 fixed-point accumulation as used by the GPU version of AMBER: Scott Le Grand, Andreas W. Götz, Ross C. Walker, "SPFP: Speed without compromise—A mixed precision model for GPU accelerated molecular dynamics simulations", Computer Physics Communications 184 (2013) 374–380 May 16, 2021 at 21:24
• Thanks for your answer Lutz, especially that you considered accumulation and evaluation in different precision. @njuffa Isn't AMBER using double for evaluation and single for accumulation ? May 18, 2021 at 7:57
• GPU version started out with the SPDP model, see here: "In addition to the traditional full double precision implementation used in the PMEMD CPU code the Amber 18 midpoint also introduces a mixed precision model, first pioneered with the Amber GPU implementation, termed SPDP. This precision model uses single precision for each particle-particle interaction but sums the resulting forces into double precision accumulators." Because DP is slow on consumer GPUs, an alternative SPFP model was developed replacing DP with 64-bit fixed-point (FP) formats. May 18, 2021 at 11:08