Kevin Schäfers is a PhD student at the University of
Wuppertal, in the Institute for Modelling, Analysis and
Computational Mathematics (IMACM).
In a long-standing collaboration at the University of Wuppertal, researchers in theoretical physics (Francesco Knechtli’s group), applied computer science (Andreas Frommer’s group), and numerical analysis (Michael Günther’s group) have been advancing the field of lattice Quantum Chromodynamics (QCD) simulations. This research, conducted as part of the FOR5269 research unit in collaboration with Stefan Schaefer (DESY Zeuthen) and Mike Peardon (Trinity College Dublin), focuses on developing “Future methods for studying confined gluons in QCD“.
This post highlights one aspect of the project: „Monte Carlo sampling of gluons by molecular dynamics“. A key focus is on designing geometric numerical integration schemes to handle the molecular dynamics step in the Hybrid Monte Carlo (HMC) algorithm [1] with greater efficiency and stability.
The HMC algorithm
The HMC algorithm is a computational method for sampling from complex probability distributions in high-dimensional spaces, crucial for physics simulations. By blending Monte Carlo sampling with molecular dynamics trajectories, HMC efficiently proposes new states with a high likelihood of acceptance, improving both accuracy and computational efficiency. This makes it especially valuable in lattice QCD, where accurate sampling in high dimensions is essential. Within HMC, the molecular dynamics step requires solving Hamiltonian equations of motion for a separable Hamiltonian H(p,q) = T(p) + V(q) with quadratic kinetic energy T(p). To satisfy the crucial detailed balance condition, the integration scheme must be volume-preserving and time-reversible. Additionally, as the phase space of the Hamiltonian system is the cotangent bundle over a base space that is the special unitary group SU(3), numerical approximations must remain on this manifold.
Splitting methods
A promising framework for geometric numerical integration in this context is splitting methods [2-3] Given the separability of the Hamiltonian, these methods can split the system into its kinetic and its potential part, each of which can be integrated exactly. Composing the flows of the respective subsystems yields a numerical approximation to the overall system that is symplectic (and thus volume-preserving) and, with a self-adjoint composition, time-reversible. These properties make splitting methods well-suited for the HMC’s molecular dynamics step.
For systems with quadratic kinetic energy, force-gradient integrators [4] can provide even greater efficiency. These integrators incorporate a certain commutator into the computational process that only depends on the coordinates q. Force-gradient integrators can be regarded as splitting methods, applied to a modified Hamiltonian with modified potential. Consequently, it does not compromise the geometric properties of the method. However, calculating the required Hessian of the potential can be computationally costly or unavailable in some software packages. In such cases, Hessian-free force-gradient integrators [5] approximate this step by using two force evaluations, eliminating the need for direct Hessian computation and simplifying integration into existing software. Due to the approximation, additional error terms occur and the integrators are no longer symplectic but still volume-preserving and hence still satisfy the detailed balance condition.
Challenges in developing efficient integrators
The main goal in lattice QCD simulations is achieving a high acceptance rate in HMC with minimal computational cost. This requires tackling several challenges:
• High-order integration: As lattice volumes grow, higher-order integrators become more efficient. For instance, on a comparatively smaller lattice (e.g., 48 x 243), fourth-order integration schemes are usually the most efficient choice.
• Numerical stability: To achieve a target acceptance rate (typically around 90%), the Hamiltonian equations only need a coarse approximation, but the proposed integrators are only conditionally stable, which can lead to instabilities at the desired acceptance rates.
• Multiple time stepping: The potential energy comprises a gauge and a fermion part, differing in magnitude and computational cost. As the fermion part causes the main part of the computational costs in the molecular dynamics step, the use of multiple time stepping techniques [6] employs different step sizes and saves costly evaluations of the fermion part while employing a sufficiently small step size to the gauge part.
• Further splitting of the fermion part: Techniques like Hasenbusch mass preconditioning [7] split the fermion part further, each partition evaluated at different step sizes to avoid instabilities and improve efficiency.
Ultimately, we aim to design high-order, Hessian-free force-gradient integrators that
• maximize accuracy at a given number of force evaluations per time step;
• expand the stability domain; and
• use multiple time stepping techniques to minimize evaluations of costly potential terms.
These goals are part of an overarching effort to understand how accuracy, stability, and multiple time stepping techniques interact, guiding the tuning process to derive efficient integrators for specific lattice QCD ensembles.
Conclusion
Hessian-free force-gradient integrators show great promise for efficiently solving equations of motion in the molecular dynamics step of the HMC algorithm. These integrators are straightforward to implement in existing software packages, as shown in openQCD (link to GitHub). While recent studies have examined their accuracy on a single time scale, ongoing research will address stability properties and extend these integrators to multiple time stepping integrators, contributing to the design of optimal integrators for long-term, large-scale lattice QCD simulations.
References
- [1] S. Duane, A.D. Kennedy, B.J. Pendleton, D. Roweth, Hybrid Monte Carlo, Physics Letters B 195.2 (1987), pp. 216-222.
- [2] R.I. McLachlan, D.R.W. Quispel, Splitting Methods, Acta Numerica 11 (2002), pp. 341-434.
- [3] S. Blanes, F. Casas, A. Murua, Splitting Methods for Differential Equations, arXiv preprint arXiv:2401.01722 (2024).
- [4] I.P. Omelyan, I.M. Mryglod, R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Computer Physics Communications 151.3 (2003), pp. 272-314.
- [5] K. Schäfers, J. Finkenrath, M. Günther, F. Knechtli, Hessian-free force-gradient integrators, arXiv preprint arXiv:2403.10370 (2024).
- [6] C. Urbach, K. Jansen, A. Shindler, U. Wenger, HMC algorithm with multiple time scale integration and mass preconditioning, Computer Physics Communications 174.2 (2006), pp. 87-98.
- [7] M. Hasenbusch, Speeding up the hybrid Monte Carlo algorithm for dynamical fermions, Physics Letters B 519.1-2 (2001), pp. 177-182.