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 ?
-
$\begingroup$ 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? $\endgroup$– Lutz LehmannMay 16, 2021 at 10:24
-
$\begingroup$ 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. $\endgroup$– Daniel WiczewMay 16, 2021 at 10:45
-
$\begingroup$ Before you consider speed can you answer how many significant figures you actually need to do the science you want to do? $\endgroup$– Ian BushMay 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:
- the theoretical method truncation error,
- the floating point error from evaluating the ODE function and composing the method step and
- 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.
-
$\begingroup$ 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 $\endgroup$– njuffaMay 16, 2021 at 21:24 -
$\begingroup$ 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 ? $\endgroup$ May 18, 2021 at 7:57
-
1$\begingroup$ 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. $\endgroup$– njuffaMay 18, 2021 at 11:08