MATERIALS SCIENCE AND CHEMISTRY

Benchmarking and Optimization of Path MD in GROMACS and Application to Ligand-Protein Unbinding Processes

Principal Investigator:
Dr. Davide Mandelli

Affiliation:
Computational Biomedicine INM-9, Forschungszentrum Jülich, Germany

Local Project ID:
pathmd-gpu

HPC Platform used:
JUWELS BOOSTER at JSC

Date published:

Teaser Abstract

Predicting how long a drug stays bound to its protein target — the drug residence time — is one of the key challenges in computational pharmacology. Standard molecular dynamics (MD) simulations are inherently serial in time and therefore cannot exploit modern exascale supercomputers to their full potential. This project implemented, optimized, and benchmarked the Metadynamics of Paths (MoP) algorithm — a highly parallelizable path-sampling method — within the widely-used GROMACS biomolecular simulation package. Executed on the pre-exascale JUWELS Booster at the Jülich Supercomputing Centre, the implementation demonstrated unprecedented weak-scaling parallel efficiency of over 70% using up to 3,500 GPUs, paving the way for exascale-ready predictions of drug residence times in realistic pharmaceutical systems.

Introduction

Modern drug design increasingly recognizes that not only how strongly a drug binds to its target, but also how long it remains bound — the drug residence time (RT) — is a critical determinant of pharmacological efficacy. Drugs with longer RTs often show enhanced selectivity and prolonged action, while those with short RTs enable rapid dose adjustment and reduced side effects. Accurately computing RTs from atomic-scale simulations would therefore provide an invaluable tool for rational drug design.

The principal challenge is timescale: most biologically relevant unbinding events occur on timescales of milliseconds to seconds, orders of magnitude beyond what conventional MD simulations can access. This bottleneck is rooted in the intrinsically sequential nature of time integration: at each simulation step, forces must be computed and atomic positions updated before the next step can begin. This seriality means that — unlike force computation, which parallelizes efficiently across many processors — the propagation of time cannot benefit from massively parallel architectures such as the exascale JUWELS Booster supercomputer.

The path molecular dynamics (path-MD) algorithm [1] offers a radically different strategy: instead of propagating a single trajectory forward in time, it maps a discretized trajectory of N steps onto the statistical mechanics of a fictitious polymer made of N replicas of the system (Figure 1). Each replica — or bead — can be assigned to independent computational resources, converting the timescale problem into a size-scale problem that is inherently and infinitely parallelizable. When combined with the free-energy method metadynamics [2], this yields the Metadynamics of Paths (MoP) algorithm, which efficiently samples reactive trajectories and enables direct computation of kinetic rates and drug residence times.

Prior proof-of-concept applications of path-MD had been limited to small systems using an implementation in the LAMMPS code [1]. The present project aimed to port the algorithm to GROMACS — a high-performance, GPU-accelerated code widely used in the biomolecular simulation community — and to demonstrate its scalability on the pre-exascale JUWELS Booster at the Jülich Supercomputing Centre.

Figure 1: The path-MD algorithm maps a discretized molecular trajectory of N steps into a fictitious elastic polymer of N beads (replicas), each assigned to an independent processing element (PE). The beads interact via springs whose equilibrium lengths encode the physical forces, and the whole polymer is evolved by standard molecular dynamics. This MPMD parallelization strategy is infinitely scalable in principle. (From Ref. [3].)

Methods

The Path-MD Algorithm

Consider a molecular system of Q atoms governed by the Smoluchowski (overdamped Langevin) equation of motion. The probability of observing a discretized stochastic trajectory R1 → R2 → ... → RN of N steps follows a Boltzmann distribution in the enlarged phase space of all N configurations, with an effective potential Veff that is isomorphic to the potential energy of a harmonic polymer. This reformulation maps the dynamical sampling problem onto a static polymer sampling problem solvable by standard finite-temperature MD.

The effective forces on each polymer bead depend on the forces on neighboring beads, requiring point-to-point communication between adjacent replicas only. This minimal communication footprint is the key reason for the algorithm's outstanding parallel scalability: each bead constitutes an independent work unit of identical cost, making path-MD a textbook weak-scaling application.

To sample the transition path ensemble — the set of all trajectories that connect the bound and unbound states of a drug-protein complex — the MoP algorithm applies metadynamics using a generalized polymer end-to-end distance Se2e = s(RN) − s(R1) as the collective variable, where s(R) distinguishes the bound and unbound states. By enhancing fluctuations of Se2e, the metadynamics bias increases the probability of sampling reactive trajectories. Kinetic rates are then recovered by standard reweighting techniques [4].

GROMACS Implementation

The path-MD algorithm was implemented in GROMACS 2024 using a Multiple Program Multiple Data (MPMD) parallelization strategy. Each polymer bead is assigned to a separate GROMACS instance running on dedicated hardware resources. Within each instance, the computation of the effective forces on a bead requires:

  • Three calls to GROMACS's internal do_force routine to compute physical forces on the bead and on two displaced configurations needed for the finite-difference Hessian approximation.
  • Seven newly implemented path-MD-specific kernel functions performing displacement computation, force assembly, spring energy accumulation, and special-case handling for the first and last beads.

All kernels are optimized for OpenMP-based multi-threading and vectorization. Inter-bead communication for the spring-like terms uses non-blocking MPI send/receive operations to minimize synchronization overhead. The full GROMACS parallelization toolkit — domain decomposition and GPU force offloading — is available for each bead's force computation.

To enable MoP simulations, we have also developed a new patch for the PLUMED library [5] compatible with GROMACS 2024, together with an optimized CUSTOM collective variable that restricts inter-replica communication to only the two beads involved in the end-to-end distance definition, avoiding costly all-to-all communication.

Benchmark System

All benchmarks were performed on a pharmacologically relevant model system: the human adenosine type 2A receptor (hA2AR), a G-protein-coupled receptor (GPCR) of interest as a drug target for Parkinson's disease, in complex with its high-affinity antagonist ZM241385, embedded in a realistic lipid bilayer and solvated in explicit water — a total of approximately 150,000 atoms (Figure 2). For the path-MD benchmarks, a polymer of 512 beads was used, yielding an effective system of ~76.8 million atoms. The optimal hardware configuration of 1 bead per GPU (4 GPUs per node) was identified from preliminary benchmarks and used throughout.

Figure 2: The benchmark system: human adenosine receptor type 2A (green ribbon) in complex with its ZM241385 inhibitor (blue), embedded in a lipid bilayer and solvated in explicit water. The simulation box contains ~150,000 atoms. (From Ref. [3].)

Findings and Impact

Validation

Prior to large-scale benchmarks, the implementation was validated by verifying energy conservation in NVE (microcanonical) simulations and temperature stability in NVT (canonical) simulations, using several test systems of increasing complexity: a Lennard-Jones cluster, water boxes, and the T4 Lysozyme/benzene complex. For the hA2AR benchmark system, a total energy drift of only ~0.05% in double precision was observed in NVE runs — comparable to standard MD — confirming the correctness of the implementation. NVT runs demonstrated stable temperature maintenance at the target value of 310 K.

Weak Scaling Performance

The weak scaling performance of the MoP implementation was tested on the JUWELS Booster supercomputer in two configurations: unbiased path-MD and metadynamics simulations with the PLUMED interface. For the unbiased case, weak scaling parallel efficiency remained above 72% across all tested polymer sizes, up to the maximum of 3,500 GPUs (875 nodes) — the largest test corresponding to approximately 94% of the entire JUWELS Booster machine.

When metadynamics is activated via PLUMED, performance is slightly reduced due to additional inter-replica communication for the collective variable evaluation, but remains consistently above 70% efficiency throughout. These results, shown in Figure 3, represent a major advance in scaling over all existing MD-based path-sampling methods, which generally scale over at most tens to hundreds of cores. The path-MD algorithm can, in principle, simulate trajectories of arbitrary length by increasing the number of beads proportionally to the computational resources available — a property unique among MD-based enhanced sampling methods.

Figure 3: Weak scaling parallel efficiency of unbiased path-MD (left) and metadynamics of paths (right) on the JUWELS Booster. Efficiency remains above 70–72% up to 3,500 GPUs (875 nodes, ~94% of the machine), using the human A2A receptor system (~150,000 atoms per bead, 30,000 steps). (From Ref. [3].)

Scientific Impact and Outlook

The results of this project establish path-MD in GROMACS as a production-ready exascale algorithm for studying rare biomolecular events. The key contributions are:

  • A fully validated, GPU-accelerated implementation of path-MD in GROMACS 2024, with domain decomposition and multi-GPU offloading support.
  • A new PLUMED patch enabling MoP simulations with optimized inter-replica communication.
  • Demonstrated weak scaling efficiency >70% up to 94% of the JUWELS Booster — the largest weak scaling test ever reported for an MD-based path-sampling method.
  • A clear route toward computing drug residence times for pharmaceutically relevant systems, using the human adenosine receptor as a benchmark.

The implementation will enable future applications of MoP to the computation of unbinding kinetics of drug-receptor complexes — including the hA2AR system already benchmarked — opening a realistic path toward exascale-based prediction of drug residence times with atomic resolution. A complementary priority is the development of optimal collective variables for use in trajectory space. While the end-to-end distance Se2e is a robust starting choice, machine-learning-based CVs derived from the transition path ensemble are expected to dramatically accelerate convergence [6]. The full GROMACS implementation now makes it feasible to apply these data-driven approaches to complex biological systems for the first time. More broadly, the path-MD algorithm represents a paradigm shift in the use of exascale resources for molecular simulation: rather than simply computing larger systems, it allows computing longer processes — directly addressing the timescale bottleneck that has limited biomolecular simulation for decades. The JUWELS Booster and future exascale systems like JUPITER are ideally suited to this new class of applications.

References

[1] Mandelli, D.; Hirshberg, B.; Parrinello, M. Metadynamics of Paths. Phys. Rev. Lett. 2020, 125, 026001. link.aps.org/doi/10.1103/PhysRevLett.125.026001

[2] Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. 2002, 99, 12562. www.pnas.org/cgi/doi/10.1073/pnas.202427399

[3] Malapally, N.; Devodier, M.; Suarez, E.; Lippert, T.; Rossetti, G.; Carloni, P.; Mandelli, D. Molecular Simulations-based Predictions of Drug's Residence Times in the Exascale Era. In Computational Modeling of Biomolecular Interactions: Methods and Applications, Wiley, 2025 (in press).

[4] Malapally, N.; Devodier, M.; Rossetti, G.; Carloni, P.; Mandelli, D. Extreme scaling of the metadynamics of paths algorithm on the pre-exascale JUWELS Booster supercomputer. arXiv:2501.11962, 2025. arxiv.org/abs/2501.11962

[5] Tribello, G.A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604. linkinghub.elsevier.com/retrieve/pii/S0010465513003196

[6] Müllender, L.; Rizzi, A.; Parrinello, M.; Carloni, P.; Mandelli, D. Effective data-driven collective variables for free energy calculations from metadynamics of paths. PNAS Nexus 2024, 3, pgae159. doi.org/10.1093/pnasnexus/pgae159

go back