Optimized M-SHAKE Constraint Implementations for GPU-Accelerated Molecular Dynamics: Balancing Precision and...

In molecular dynamics (MD) simulations, enforcing holonomic constraints on specific interatomic distances is crucial to extend the integration time step while maintaining numerical stability. While double-precision arithmetic has traditionally been used to ensure numerical accuracy, there is a growing demand for single- and mixed-precision calculations as GPU-accelerated MD simulations become increasingly widespread. In this study, we present improvements to the constraint solvers through an analytical treatment for bonds involving a single hydrogen atom ( $X{H}_1$ groups). For more complex groups with multiple hydrogens ( $X{H}_2$ and $X{H}_3$ ), we introduce two precision-level implementation strategies: (1) a mixed-precision approach, where variables are generally described in single-precision while temporary coordinate updates are handled in double-precision and (2) a pairwise arithmetic approach, where temporary coordinate variables are represented as pairs of single-precision values to reduce rounding errors. These implementation strategies reduce velocity-update errors compared with schemes using single precision for all variables. Energy drift analysis in microcanonical simulations confirms that these strategies maintain accuracy comparable to full double-precision implementations. Performance benchmarks on GPUs with limited FP64 throughput (RTX 3080 and A6000 Ada) show that our proposed implementation strategies substantially reduce runtime compared with full double-precision calculations. Conversely, on architectures with balanced FP32 and FP64 performance, such as CPUs or H100 GPUs, full double-precision arithmetic remains the preferred choice for the constraint kernels evaluated in this study. Overall, these findings demonstrate that the proposed strategies provide an effective balance between numerical accuracy and computational performance on GPUs with limited FP64 throughput.