⚗️ biochemistry

Gaussian Accelerated Molecular Dynamics in GROMACS

This paper presents the first GPU-enabled implementation of Gaussian accelerated molecular dynamics (GaMD) in GROMACS 2025.4, demonstrating its effectiveness for enhanced sampling of protein folding and ligand binding through successful free-energy recovery across four benchmark systems without requiring predefined collective variables.

Original authors: Yang, Y.

Published 2026-08-10

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Yang, Y.

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ⚕️ This is an AI-generated explanation of a preprint that has not been peer-reviewed. It is not medical advice. Do not make health decisions based on this content. Read full disclaimer

Imagine the microscopic world inside your body as a bustling, chaotic city made entirely of tiny, vibrating Lego bricks. These bricks are atoms, and when they snap together in specific patterns, they form proteins—the molecular machines that keep you alive. To understand how these machines work, scientists use a powerful tool called Molecular Dynamics (MD). Think of MD as a super-fast movie camera that records every single move these atoms make. However, there's a catch: the city is so huge and the atoms move so slowly that the camera often gets stuck in one neighborhood, unable to film the exciting events happening across town, like a protein folding into its final shape or a drug molecule finding its target. It's like trying to find a specific key in a dark, giant attic by only looking at one corner for a few seconds.

To solve this, scientists invented "enhanced sampling" methods. One popular technique is called Gaussian Accelerated Molecular Dynamics, or GaMD. If the standard movie camera is too slow, GaMD is like giving the atoms a gentle, smooth push whenever they get stuck in a low-energy valley. It doesn't force them in a specific direction (which would ruin the movie); instead, it just makes it easier for them to hop over the hills and explore the rest of the attic. The magic is that this push is designed to be mathematically predictable, so scientists can later "rewind" the movie and calculate exactly what the un-pushed, natural world looks like. For a long time, this clever trick was missing from GROMACS, one of the most popular and powerful cameras used by scientists to film these atomic movies.

This paper introduces a brand-new version of GROMACS (specifically version 2025.4) that finally has GaMD built right into its engine. The author, Yuefeng Yang, didn't just add a plugin; they rebuilt the engine to handle this "gentle push" directly on the computer's graphics card (GPU), making it incredibly fast and efficient. To prove this new feature works, the author ran four different "test drives" using simulated molecules. First, they tested a tiny molecule called alanine dipeptide, showing that the new method could map its energy landscape in just 100 nanoseconds of simulation time, a result that matched what took 1,000 nanoseconds using the old, un-pushed method. Next, they watched two small proteins, chignolin and TC5b, fold themselves from a messy, stringy state into their perfect, functional shapes. In these simulations, the proteins folded in 300 nanoseconds and 1 microsecond respectively, capturing the entire process from start to finish. Finally, they simulated a benzene molecule trying to find its way into a protein pocket (T4 lysozyme). In two out of five attempts, the benzene successfully swam into the pocket and found the exact spot it occupies in real life, matching the crystal structure with a tiny error of just 0.06 nm.

The results suggest that this new GROMACS-GaMD tool is a reliable way to speed up the filming of molecular movies without losing the accuracy of the story. The simulations showed that the "push" given to the atoms followed a predictable, bell-curve pattern, which allowed the scientists to mathematically correct the footage and recover the true energy maps of the molecules. While the method worked beautifully for these specific test cases, the paper notes that it is currently limited to three specific modes of boosting and doesn't yet include some specialized variants designed for very specific biological tasks. However, by successfully integrating this method into a widely used engine, the work provides a practical, high-speed toolkit for researchers to study how proteins fold and how drugs bind, potentially helping them solve complex biological puzzles much faster than before.

Technical Summary: Gaussian Accelerated Molecular Dynamics in GROMACS

Problem Statement
Atomistic molecular dynamics (MD) simulations are limited by restricted accessible time scales, often failing to adequately sample rare biologically relevant events such as conformational transitions, protein folding, and ligand binding due to high free-energy barriers. While Gaussian accelerated molecular dynamics (GaMD) offers a solution by enhancing sampling without predefined collective variables, its software support has been uneven across major MD engines. Although AMBER, NAMD, and OpenMM offer GaMD implementations, GROMACS—one of the most widely used high-performance MD engines—lacked a native, engine-integrated GaMD module. Existing enhanced-sampling tools in GROMACS (e.g., AWH, replica exchange, PLUMED interface) do not provide a dedicated, collective-variable-free GaMD action.

Methodology and Implementation
This work presents the first engine-level implementation of GaMD in GROMACS 2025.4. Unlike external wrappers, the method is integrated directly into the simulation engine to ensure strict consistency between bias forces, reported energies, and statistical updates.

Theoretical Framework: The implementation supports three modes: total-boost (applying to the total potential energy), dihedral-boost (applying to dihedral energy), and dual-boost (simultaneously applying to both). The boost potential, $\Delta V$ , is a smooth harmonic function applied when the instantaneous potential energy $V(r)$ falls below a threshold $E$ . The parameters (threshold $E$ and force constant $k$ ) are dynamically determined from potential energy statistics (minimum, maximum, average, and standard deviation) collected during a preparatory stage to ensure the boost potential remains approximately Gaussian, enabling free-energy recovery via cumulant expansion.
Technical Architecture: The implementation utilizes a staged workflow: (1) conventional MD (cMD) for statistics collection, (2) boost equilibration for parameter refinement, and (3) production GaMD.
- GPU Acceleration: Dedicated GPU kernels evaluate boost factors and apply force corrections. Nonbonded, PME, bonded, and update calculations remain offloaded to the GPU, minimizing host-device transfers and enabling CUDA Graph acceleration.
- Data Management: The system handles restarts reproducibly across parallel and heterogeneous environments. It outputs necessary state files for reweighting, including the evolution of energy statistics, boost potentials, and force-weight terms.
Validation Protocols: The implementation was benchmarked on four systems using the CHARMM36m force field and TIP3P water:
1. Alanine Dipeptide: To test free-energy reweighting accuracy (100 ns GaMD vs. 1000 ns cMD).
2. Chignolin: To test protein folding sampling (3 independent 300 ns trajectories).
3. TC5b: To test folding in a larger, more complex system (3 independent 1 $\mu$ s trajectories).
4. Benzene–T4 Lysozyme: To test ligand binding/unbinding pathways (5 independent 500 ns trajectories).

Key Results

Alanine Dipeptide: A 100 ns reweighted GaMD trajectory successfully recovered the major free-energy basins and rotational barriers on the $(\phi, \psi)$ surface, showing overall agreement with a 1000 ns conventional MD reference. The boost potential distribution was unimodal and Gaussian (anharmonicity $8.91 \times 10^{-3}$ ), validating the second-order cumulant reweighting.
Chignolin: All three independent 300 ns trajectories sampled native-like folded states from extended conformations. The best folded structure achieved a backbone RMSD of 0.07 nm from the experimental NMR structure. The average folding time was reduced to ~100 ns compared to ~600 ns in previous cMD studies. The reweighted PMF clearly resolved folded, intermediate, and unfolded basins.
TC5b: All three independent 1 $\mu$ s trajectories reached the low-RMSD folded region. The best structure had a backbone RMSD of 0.03 nm from the experimental crystal structure. The method captured folding events within 1 $\mu$ s, significantly faster than the experimental folding time (~~3.6 $\mu$ s) and previous cMD estimates (~~5.5 $\mu$ s).
Benzene–T4 Lysozyme: Two of five independent 500 ns trajectories captured both ligand binding and dissociation events. The bound pose achieved a minimum ligand RMSD of 0.06 nm from the crystal structure. The reweighted PMF distinguished unbound, intermediate, and bound states, revealing a binding barrier of ~3.5 kcal/mol.
General Observations: Across all systems, the boost-potential distributions were approximately Gaussian. Second-order cumulant reweighting successfully resolved expected conformational and binding free-energy basins. Local anharmonicity remained low (<0.1) in major energy wells, supporting the reliability of the reweighting in thermodynamically dominant regions.

Significance and Claims
The author claims that this work fills a critical gap in the GROMACS ecosystem by providing a native, GPU-enabled, collective-variable-free enhanced-sampling framework. The implementation is presented as a practical tool for:

Biomolecular Free-Energy Calculations: Accurately recovering free-energy landscapes (e.g., alanine dipeptide) with sampling times an order of magnitude shorter than conventional MD.
Protein Folding: Efficiently sampling folding pathways and native states for peptides (chignolin) and larger proteins (TC5b).
Ligand Binding: Characterizing binding pathways and recovering bound poses for protein-ligand systems.

The paper emphasizes that the implementation satisfies rigorous engineering requirements, including reproducible restarts, consistent bias forces, and complete output for post-processing. While the current version implements the standard total, dihedral, and dual-boost modes, the author notes that future work could extend this to selective variants (e.g., LiGaMD, Pep-GaMD). The study concludes that GROMACS-GaMD is a verifiable and reusable method that expands the capabilities of GROMACS for studying complex biomolecular processes without the need for predefined reaction coordinates.

More like this