Heterogeneous Parallelization and Acceleration of Molecular Dynamics Simulations in GROMACS

Szilárd Páll; Artem Zhmurov; Paul Bauer; Mark Abraham; Magnus Lundborg; Alan Gray; Berk Hess; Erik Lindahl

Source-linked AI summary

Heterogeneous Parallelization and Acceleration of Molecular Dynamics Simulations in GROMACS

Szilárd Páll, Artem Zhmurov, Paul Bauer, Mark Abraham, Magnus Lundborg, Alan Gray, Berk Hess, Erik Lindahl

arXiv:2006.09167v2Jun 16, 2020physics.comp-phcs.DCcs.DScs.PF

TL;DR

GROMACS needed algorithms that could exploit increasingly parallel CPUs and GPUs while managing computation and data-movement constraints. The paper presents heterogeneous acceleration across CPUs and GPUs, achieving at least doubled performance on the same hardware and excellent scaling to 16 GPUs in one benchmark.

Problem
Increasing hardware parallelism and widening computation–data-movement gaps require molecular-dynamics algorithms to exploit multiple levels of parallelism efficiently.
Method
The paper presents a heterogeneous GROMACS design combining cluster pair algorithms, CPU–GPU workload balancing, tolerance-based buffers, and multi-device acceleration.
Results
Performance at least doubled on the same hardware, while direct GPU communication enabled excellent scaling to 16 GPUs for the reaction-field STMV benchmark.
Takeaways & Limitations
The implementation provides substantial absolute performance and relative scaling across flexible combinations of CPU and GPU resources.
Takeaways & Limitations
The authors identify 3D FFT many-to-many communication in PME as limiting strong scaling and anticipate fast multipole methods for the largest parallel runs.

Abstract

from arXiv · show

The introduction of accelerator devices such as graphics processing units (GPUs) has had profound impact on molecular dynamics simulations and has enabled order-of-magnitude performance advances using commodity hardware. To fully reap these benefits, it has been necessary to reformulate some of the most fundamental algorithms, including the Verlet list, pair searching and cut-offs. Here, we present the heterogeneous parallelization and acceleration design of molecular dynamics implemented in the GROMACS codebase over the last decade. The setup involves a general cluster-based approach to pair lists and non-bonded pair interactions that utilizes both GPUs and CPU SIMD acceleration efficiently, including the ability to load-balance tasks between CPUs and GPUs. The algorithm work efficiency is tuned for each type of hardware, and to use accelerators more efficiently we introduce dual pair lists with rolling pruning updates. Combined with new direct GPU-GPU communication as well as GPU integration, this enables excellent performance from single GPU simulations through strong scaling across multiple GPUs and efficient multi-node parallelization.

I. INTRODUCTION

GROMACS’s heterogeneous acceleration redesign targets high absolute performance and broad feature support by combining CPUs, GPUs, and multiple accelerator devices. Its bottom-up approach uses hardware-tailored algorithms and parallel implementations to improve flexibility, portability, and performance from laptops to supercomputers.

Motivation: The design prioritizes cost-efficiency and throughput because producing trajectories at a high rate is imperative for sampling adequate timescales.The stated motivation is to make effective use of limited compute resources and improve research productivity.
Heterogeneous acceleration: GROMACS redesigned its algorithms from scratch to exploit CPUs and GPUs, including multiple devices, while supporting virtually all simulation features.The redesign aims to provide very high absolute performance regardless of available hardware.
Algorithm redesign: The cluster pair algorithm replaces particle-based Verlet lists and enables hardware-optimized clusters, accuracy-based cut-offs, load balancing, and on-the-fly tuning.CPU SIMD parallelization and multithreading support efficient GPU offloading of short-range non-bonded calculations.
Heterogeneous acceleration: The framework balances CPU and GPU utilization by offloading more work to increasingly powerful accelerators for dense heterogeneous machines.The approach is designed around the differing strengths of latency-optimized CPUs and throughput-optimized GPUs.
Scalability and portability: The latest release supports most MD-engine features across multiple accelerator architectures, with scaling across multiple accelerators and nodes.This extends the approach from single-device simulations to multi-device and multi-node execution.

II. COMPUTATIONAL CHALLENGES IN MD SIMULATIONS

MD simulations are dominated by per-step force calculations, although integration remains inherently iterative. As optimized engines reach sub-millisecond step rates, performance increasingly depends on strong scaling, fine-grained concurrency, and efficient data movement across evolving multicore and accelerator hardware.

Computational structure: Per-step force calculation is the main computational cost, while numerical integration of Newton’s equations remains inherently iterative.The core simulation advances particle systems by integrating Newton’s equations of motion at every time step.
Computational structure: Total particle force combines non-bonded, bonded, restraint, and external-force terms that can be computed independently from particle coordinates.Non-bonded terms include Lennard-Jones and electrostatic interactions; bonded terms include bonds, angles, and torsions.
Scaling challenges: At or below 1 millisecond per step, optimized MD engines make strong scaling increasingly important for improving time-to-solution.Fixed, relatively small simulation sizes are commonly used to improve absolute performance, while force calculations historically dominated runtime.
Parallelization challenges: Multicore and accelerator architectures shift MD optimization toward fine-grained concurrency within force-calculation tasks rather than only coarse MPI and domain-decomposition parallelism.The shift is enabled by faster intra-node communication and is better suited to multicore and accelerator hardware.
Hardware and algorithmic challenges: Modern hardware performance increasingly depends on data layout, reduced data movement, communication-computation overlap, and recomputation instead of relying mainly on saving arithmetic operations.Processors, SIMD units, and GPUs provide more specialized compute capacity, but feeding compute units with data remains a challenge.

III. PARALLELIZATION OF MD IN GROMACS · A. The structure of the MD algorithm · B. Multi-level parallelism

GROMACS parallelizes molecular-dynamics force computation across independent tasks while explicitly targeting multiple hardware-parallelism levels. Its heterogeneous design coordinates SIMD/SIMT execution, multithreading, domain decomposition, and CPU–GPU offload to improve performance.

A. The structure of the MD algorithm: Force terms expose task parallelism within each molecular-dynamics step, with dependencies primarily on prior-step positions and, under domain decomposition, communicated coordinates.Force reduction for integration limits or coordinates the otherwise concurrent force calculations.
B. Multi-level parallelism: GROMACS targets multiple hardware-parallelism levels, matching computation and data movement to hierarchical memory and intra- and inter-node communication structures.The approach is motivated by the increasing importance of parallelism on petascale architectures.
B. Multi-level parallelism: CPUs provide fine-grained SIMD execution and multiple cores or hardware threads, while GPUs use SIMT groups executing threads in lockstep across execution units.GPU work distribution across execution units cannot be controlled directly in the same way as CPU work distribution.
B. Multi-level parallelism: Heterogeneous systems must manage CPU–GPU and GPU–GPU buses, separate global memory, network communication, and non-uniform data-movement costs.These communication and NUMA effects complicate coordination both within nodes and across nodes.
B. Multi-level parallelism: A redesign extended GROMACS parallelization to separately target each hardware level and introduced finer-grained algorithms and programming models.The redesign also made pair-interaction calculation flexible, accuracy-configurable, and capable of load balancing.
B. Multi-level parallelism: An accelerator-offload layer and multithreading schedule GPU tasks and data movement to enable concurrent CPU–GPU execution.The offload layer was added to extend intra-node parallelization across CPUs and GPUs.
B. Multi-level parallelism: OpenMP multithreading expanded from PME scaling to the MD engine, while larger MPI tasks running across multiple cores reduced scheduling bottlenecks for small GPU tasks.These changes assembled larger computation units for GPU execution instead of dispatching work from many MPI ranks per node.

IV. HETEROGENEOUS PARALLELIZATION · A. Offloading force computation

GROMACS developed heterogeneous force offloading around asynchronous CPU–GPU execution, explicit data movement, and progressively broader GPU task support. The design improves overlap and flexibility but can become limited by per-step transfers and CPU–GPU bandwidth imbalance.

IV. HETEROGENEOUS PARALLELIZATION: Asynchronous offloading uses CUDA or OpenCL to control CPU–GPU data movement and schedule concurrent task execution and synchronization.The design moves data early, keeps it on the accelerator, overlaps transfers with computation, and optimizes task scheduling.
A. Offloading force computation: GROMACS initially offloaded non-bonded pair interactions to GPUs while evaluating PME and bonded interactions on CPUs.This approach supports domain decomposition and general CPU execution of special algorithms, but requires CPU resources.
A. Offloading force computation: Domain decomposition splits non-bonded GPU work into local-only and non-local interaction kernels to handle halo-exchanged particles.Stream-priority support was introduced through co-design with NVIDIA and exposed in CUDA.
A. Offloading force computation: PME long-range electrostatics was added to CUDA and OpenCL offloading, using a separate stream to overlap PME with short-range interactions.The main challenge is handling the two required 3D FFTs, with small simulation grids motivating further considerations.
A. Offloading force computation: Bonded interactions became the last force task to be offloaded, sharing the short-range non-bonded stream without domain decomposition.The passage describes this execution flow as part of the incremental GPU offload design.
A. Offloading force computation: Every offloaded force step transfers coordinates to and forces from the GPU before force reduction, followed by CPU integration.These transfers can make accelerator-heavy systems CPU-limited.
A. Offloading force computation: GPU compute-to-PCIe bandwidth imbalance and limited high-performance interconnect availability disadvantage offload-based execution.High-performance interconnects are typically reserved for GPU–GPU communication, constraining CPU–GPU transfer efficiency.

B. Offloading complete MD iterations · V. ALGORITHM DETAILS · A. The cluster pair algorithm

The design accelerates complete molecular-dynamics iterations on GPUs while using CPUs for infrequent supporting tasks and direct GPU communication to reduce transfers and synchronization. Its cluster pair algorithm regularizes neighbor interactions, improves data reuse, and tunes execution for GPU and SIMD architectures.

B. Offloading complete MD iterations: Complete MD iterations, including integration, can execute on accelerators, removing the CPU from the critical path and reducing synchronization events.The CPU remains available for infrequent pair search and domain decomposition, as well as special algorithms during GPU execution.
B. Offloading complete MD iterations: Communication can route directly between GPUs, including through dedicated high-bandwidth/low-latency interconnects, instead of staging through CPU memory.With CUDA-aware MPI, communication operates directly in GPU memory spaces.
A. The cluster pair algorithm: Traditional Verlet-list formulations expose irregular parallelism and varying list lengths, motivating fixed-size work units better suited to SIMD-like architectures.Padding and binning reduce imbalance, but fixed-size interaction groups have become the prevailing reformulation.
A. The cluster pair algorithm: The cluster pair algorithm increases arithmetic-to-memory ratio and data reuse by making fixed-size cluster pairs, rather than individual particles, the computational unit.Its goal is to load j-particle data efficiently and reuse it across multiple i-particles, analogous to blocking in matrix-matrix multiplication.
V. ALGORITHM DETAILS: GPU execution computes all particle pairs within an i–j cluster pair, avoiding the low-efficiency zero interactions that would result from matching cluster size to wide GPU execution widths.Matching the cluster size to GPU execution width would raise raw FLOP-rate but reduce efficiency.
V. ALGORITHM DETAILS: The GPU implementation swaps loop order to traverse joint j-cluster lists outside and eight i-clusters inside, reducing interaction-mask storage and enabling register accumulation.The interaction mask needs 8 bits per j-cluster, while i-particle partial forces require only 8 × 8 forces of memory.
A. The cluster pair algorithm: Approximately four of eight i-clusters typically interact with j-particles, so about half of inner-loop checks skip, costing an estimated 8–12%.Self-exclusions and force-field exclusions are enforced within the interaction kernels and cluster-pair lists.
V. ALGORITHM DETAILS: The cluster algorithm provides multiprocessor-level parallelism through super-cluster work assigned to thread blocks and uses separate kernels for electrostatics, Lennard-Jones interactions, cut-offs, and energy requirements.It is implemented in CUDA and OpenCL and tuned for multiple GPU architectures.

B. Algorithmic work efficiency and pair-list buffers

The cluster algorithm improves SIMD-oriented execution efficiency by accepting some zero-valued interaction work and exploiting particles outside the buffered cut-off as an implicit buffer. Accuracy is controlled through user-requested tolerance settings, which can reduce explicit buffering and support consistent optimization across systems.

Work efficiency: Parallel work efficiency measures the fraction of non-zero interactions evaluated, and remains ≤1 even for standard Verlet lists because finite buffers include non-interacting particles.The cluster algorithm augments this trade-off with additional work associated with clustered pair lists.
Pair-list buffers: Particles outside the buffered cut-off provide an implicit buffer, allowing shorter explicit buffers while maintaining algorithmic accuracy.This improves computational efficiency across modern architectures by offsetting the cluster algorithm’s work-efficiency trade-off.
Pair-list buffers: 0.218nm is required by the 1 × 1 setup to reach the same error tolerance achieved by the 8 × 4 cluster setup with a 0.105nm buffer.The comparison demonstrates the cluster setup’s ability to maintain the same error tolerance with a shorter explicit buffer.
Accuracy control: Tolerance-based settings are preferred over heuristic rules because the worst algorithmic component can dominate final simulation accuracy.The acceptable error varies between problems, so improving accuracy in only selected parts provides little benefit.
Accuracy control: The default tolerance is 0.005kJ/mol/ps per atom, while the observed energy drift is 0.0001kJ/mol/ps per atom for the stated water or solvated-biomolecule conditions.For PME electrostatics and 20fs pair-list update intervals, the implicit buffer is sufficient and no additional explicit buffer is required.

C. Non-bonded pair interaction kernel throughput

Non-bonded pair-kernel throughput depends on particle density and hardware-specific saturation behavior. CPUs reach peak raw throughput at lower pair counts, whereas GPUs require more work and lose efficiency on very small systems.

Density dependence: Throughput varies with particle density across coarse-grained, all-atom, and liquid-crystal or gas simulations.Performance was investigated using a 157k-particle Lennard-Jones system to reduce input-size effects.
Kernel saturation: ∼150 pairs per particle reaches peak raw throughput on the CPU, while the GPU saturates only at ∼1000 pairs.The comparison uses CUDA GPU and AVX512 CPU kernels with same-size clusters and identical work-efficiency.
Kernel saturation: Effective pair throughput increases monotonically as pair counts grow because more pairs become non-zero interactions.The raw-throughput comparison includes zero interactions, whereas effective throughput counts only non-zero interactions.
System-size effects: GPU kernels achieve significantly higher throughput from moderate input sizes, while CPU kernels have less overhead for small systems.The system-size comparison uses SPC/E water with a 1 nm cut-off, Ewald electrostatics, 40 step search frequency, and default tolerances.
System-size effects: Sub-10k atom systems lack enough parallelism to execute in a balanced manner on the largest GPUs.Kernel invocation, fixed-cost operations, and pair-list balancing create small-system overhead; balancing is limited when there is insufficient work.

D. The pair list generation algorithm

GROMACS generates pair lists with a fixed lifetime and an accuracy-based buffer, allowing pair-search frequency to be chosen freely. The CPU-based, SIMD-optimized search builds hardware-adapted hierarchical lists whose pruning is reused across steps and reduces list size by 50–75%.

Pair-list lifetime and buffering: GROMACS uses a fixed pair-list lifetime and an accuracy-based buffer estimate, allowing pair-search frequency to be chosen freely.The fixed lifetime avoids displacement-based heuristic updates, while frequency affects buffer size rather than requiring careful tuning.
CPU pair search: Pair searching remains on the CPU to support parallelization, portability, and maintainability while reducing GPU idle time through algorithmic improvements.The implementation uses SIMD optimization and builds the hierarchical GPU list using cluster bounding-box distances.
Hardware-adapted pruning: 50–75%: pruning reduces the pair-list size by this amount, depending on the cut-off and buffer length.The search splits j-clusters for hardware execution width, performs initial pruning when the list is first processed, and reuses the pruned list on consecutive MD steps.

E. Dual pair list with dynamic pruning

The section addresses inefficient domain decomposition and pair-list generation by examining the trade-offs of increasing search frequency and introducing an extension to the cluster algorithm.

E. Dual pair list with dynamic pruning: Irregular data access makes domain decomposition and pair-list generation improve more slowly than compute kernels.Their performance has not improved at the same rate as compute kernels.
E. Dual pair list with dynamic pruning: Increasing search frequency trades pair-list costs for additional pair-interaction work but introduces significant drawbacks.The approach becomes costly as the buffer increases and has a narrow optimum.
E. Dual pair list with dynamic pruning: As the buffer increases, overhead becomes large, while the trade-off remains sensitive to inputs and runtime conditions.The passage reports a small optimal window for this trade-off.
E. Dual pair list with dynamic pruning: The authors developed an extension to the cluster algorithm to address these pair-list generation limitations.The supplied passage truncates immediately after identifying this extension.

F. Multi-level load balancing

GROMACS addresses load imbalance across data, task, GPU, PME, and domain-decomposition parallelism through indirect scheduling, overlapping execution, online PME rebalancing, and domain resizing. Dual pair lists with pruning reduce the cost and sensitivity of shifting work toward direct-space interactions, while PME balancing can improve performance by up to 25% on balanced CPU–GPU systems.

F. Multi-level load balancing: Irregular pair lists, particle densities, bonded interactions, and heterogeneous CPU–GPU tasks create load imbalance across multiple parallelism levels.Examples include membrane-protein systems with united-atom lipids and uneven bond distributions between solvent and protein.
F. Multi-level load balancing: Sorting pair lists and heuristically splitting them increase GPU parallelism when irregular work or too few cluster lists cause inefficient scheduling.Sorting prevents long lists from being scheduled late, while list splitting increases available parallelism for the target hardware.
F. Multi-level load balancing: Asynchronous multiple streams and event-based dependencies maximize opportunities to overlap GPU kernels and mitigate kernel-tail effects.The approach requires enough concurrent work to fill otherwise idle GPU cores.
F. Multi-level load balancing: PME load balancing shifts work between reciprocal-space CPU tasks and direct-space GPU tasks by scaling the cut-off with PME grid spacing.The online load balancer automates this redistribution, including configurations with dedicated PME ranks.
F. Multi-level load balancing: Up to 25% performance gain is possible from PME load balancing on balanced CPU–GPU setups, despite increased communication and higher real-space computational cost.A cut-off scaling threshold avoids increasing GPU load in heavily CPU-bound runs.
F. Multi-level load balancing: Dual pair lists with pruning make load balancing less sensitive and enable hundreds of steps between updates, reducing CPU load and maximizing GPU utilization.This avoids most drawbacks of shifting work to direct space because pruning is substantially cheaper than evaluating interactions.

VI. PERFORMANCE BENCHMARKS

The benchmarks show that GROMACS’s algorithmic and software advances substantially improve performance on identical hardware. Performance depends on offload configuration, workload concurrency, and direct GPU communication for multi-GPU scaling.

Hardware evolution: At least doubled performance on identical hardware, with latest-generation consumer GPUs reaching almost fourfold improvement over five years.These gains apply even to older cost-efficient GPUs and are shown for the 142k-atom GluCl benchmark.
Single-GPU offload: Offloading non-bonded interactions significantly increases single-GPU performance, while also offloading PME removes CPU oversubscription and can make the GPU the bottleneck.Offloading bonded forces can instead regress performance, particularly for the membrane-protein workload with many torsions.
Ensemble throughput: Running many uncoupled or loosely coupled trajectories on one GPU can increase cumulative throughput by overlapping one trajectory’s transfers with another’s computation.For older GPUs, the difference between worst- and best-performing offload cases is only about 25%.
Multi-GPU scaling: Direct GPU communication enables efficient multi-GPU scaling with outstanding absolute performance across up to four nodes equipped with four NVIDIA Tesla V100 GPUs per node.The STMV benchmark contains 1M atoms and uses 2 fs steps, with NVLink within nodes and MPI over inter-node connections.

VII. DISCUSSION

The discussion presents GROMACS as a major step toward high-performance, flexibly scaled molecular dynamics on heterogeneous CPU–GPU systems. It also emphasizes remaining time-scale, tuning, parallelization, and large-system scaling challenges that require further algorithmic and hardware advances.

Current capabilities: Microsecond-scale simulations are routine on commodity hardware, but bridging hundreds of microseconds to millisecond single trajectories still requires special-purpose hardware.General-purpose codes retain advantages in flexibility, adaptability, and portability.
Exascale direction: Efficient heterogeneous parallelization and large-scale ensemble algorithms are expected to characterize molecular dynamics in the exascale era.The discussion links these advances with new algorithms and continued performance improvements in core MD codes.
Usability and future direction: GROMACS flexibility creates developer and user challenges because algorithmic and parallelization settings must adapt to increasingly diverse hardware.The approach uses heuristics-based defaults while moving toward tolerance-based settings.
Future challenges: Reaching < 100µs iterations for larger systems will require improved GPU scheduling, CPU tasking, communication overlap, and long-range interaction algorithms.GROMACS already achieves such iteration rates for small systems such as the ∼5000-atom villin headpiece.
Future challenges: PME is expected to remain the low-scale long-range method, while fast multipole methods may become preferable at the largest parallel scales because 3D FFT communication limits strong scaling.The discussion identifies 3D FFT many-to-many communication as the relevant strong-scaling limitation.
Current capabilities: GROMACS provides major absolute-performance and relative-scaling advances by supporting almost arbitrary-balance combinations of CPU and GPU functional units.The implementation is intended to support applications ranging from consumer GPU hardware to large HPC resources.