Table of Contents
Fetching ...

Constant pH Simulation with FMM Electrostatics in GROMACS. (B) GPU Accelerated Hamiltonian Interpolation

Bartosz Kohnke, Eliane Briand, Carsten Kutzner, Helmut Grubmüller

TL;DR

This article presents MAHI, a GPU-accelerated MAHI method that enables rigorous Hamiltonian Interpolation for constant-pH MD within GROMACS by leveraging Fast Multipole Method electrostatics. MAHI starts from a charge-scaled Hamiltonian and applies targeted P2P, lattice, and dipole corrections to recover HI forces, achieving near-linear scaling with the number of titratable sites and forms. Accuracy benchmarks show MAHI closely reproduces reference HI/PME results across random and biomolecular test systems, while performance evaluations reveal modest overhead (often under 10–25%) that remains nearly constant as system size grows; a hybrid PME+MAHI approach can yield further speedups. The work demonstrates practical, scalable constant-pH MD suitable for large biomolecular systems and opens avenues for exascale-ready constant-pH simulations with flexible electrostatics solvers.

Abstract

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or their complexes, are strongly influenced by protonation changes of their typically many titratable groups, which explains their pH sensitivity. In turn, conformational and environmental changes in the biomolecule affect the protonation state of these groups. With a few exceptions, conventional force field-based molecular dynamics (MD) simulations do not account for these effects, nor do they allow for coupling to a pH buffer. The $λ$-dynamics method implements this coupling and thus allows for MD simulations at constant pH. It uses separate Hamiltonians for the protonated and deprotonated states of each titratable group, with a $λ$ variable that continuously interpolates between them. However, rigorous implementations of Hamiltonian Interpolation (HI) $λ$-dynamics are prohibitively slow when used with Particle Mesh Ewald (PME). To circumvent this problem, it has been proposed to interpolate the charges instead of the Hamiltonians (QI). Here, we propose a rigorous yet efficient Multipole-Accelerated Hamiltonian Interpolation (MAHI) method to perform $λ$-dynamics in GROMACS. Starting from a charge-scaled Hamiltonian, precomputed with the Fast Multipole Method (FMM) or with PME, the correct HI forces are calculated with negligible computational overhead. We compare HI with QI and show that HI leads to more frequent transitions between protonation states, resulting in better sampling and accuracy. Our performance benchmarks show that introducing, e.g., 512 titratable sites to a one million atom MD system increases runtime by less than 20% compared to a regular FMM-based simulation. We have integrated the scheme into our GPU-FMM code for the simulation software GROMACS, allowing an easy and effortless transition from standard force field simulations to constant pH simulations.

Constant pH Simulation with FMM Electrostatics in GROMACS. (B) GPU Accelerated Hamiltonian Interpolation

TL;DR

This article presents MAHI, a GPU-accelerated MAHI method that enables rigorous Hamiltonian Interpolation for constant-pH MD within GROMACS by leveraging Fast Multipole Method electrostatics. MAHI starts from a charge-scaled Hamiltonian and applies targeted P2P, lattice, and dipole corrections to recover HI forces, achieving near-linear scaling with the number of titratable sites and forms. Accuracy benchmarks show MAHI closely reproduces reference HI/PME results across random and biomolecular test systems, while performance evaluations reveal modest overhead (often under 10–25%) that remains nearly constant as system size grows; a hybrid PME+MAHI approach can yield further speedups. The work demonstrates practical, scalable constant-pH MD suitable for large biomolecular systems and opens avenues for exascale-ready constant-pH simulations with flexible electrostatics solvers.

Abstract

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or their complexes, are strongly influenced by protonation changes of their typically many titratable groups, which explains their pH sensitivity. In turn, conformational and environmental changes in the biomolecule affect the protonation state of these groups. With a few exceptions, conventional force field-based molecular dynamics (MD) simulations do not account for these effects, nor do they allow for coupling to a pH buffer. The -dynamics method implements this coupling and thus allows for MD simulations at constant pH. It uses separate Hamiltonians for the protonated and deprotonated states of each titratable group, with a variable that continuously interpolates between them. However, rigorous implementations of Hamiltonian Interpolation (HI) -dynamics are prohibitively slow when used with Particle Mesh Ewald (PME). To circumvent this problem, it has been proposed to interpolate the charges instead of the Hamiltonians (QI). Here, we propose a rigorous yet efficient Multipole-Accelerated Hamiltonian Interpolation (MAHI) method to perform -dynamics in GROMACS. Starting from a charge-scaled Hamiltonian, precomputed with the Fast Multipole Method (FMM) or with PME, the correct HI forces are calculated with negligible computational overhead. We compare HI with QI and show that HI leads to more frequent transitions between protonation states, resulting in better sampling and accuracy. Our performance benchmarks show that introducing, e.g., 512 titratable sites to a one million atom MD system increases runtime by less than 20% compared to a regular FMM-based simulation. We have integrated the scheme into our GPU-FMM code for the simulation software GROMACS, allowing an easy and effortless transition from standard force field simulations to constant pH simulations.
Paper Structure (32 sections, 52 equations, 12 figures, 1 table)

This paper contains 32 sections, 52 equations, 12 figures, 1 table.

Figures (12)

  • Figure 1: FMM far field calculation. The five individual steps and operators involved in the far field calculation, shown for the lowest two levels of the octree: /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $1$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $1$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $1$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $1$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) 1 P2M: At the lowest level, the individual charges (yellow dots) are combined into a multipole representation. /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $2$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $2$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $2$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $2$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) 2 M2M: The multipoles of the higher levels are derived from those of the lower levels (blue). /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $3$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $3$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $3$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $3$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) 3 M2L: The multipoles (blue) are transformed into local moments (pink) at each level of the tree. /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $4$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $4$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $4$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $4$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) 4 L2L: The local moments are propagated down the tree to the deepest level. /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $5$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $5$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $5$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) $5$ /csteps/inner xsep /csteps/inner ysep /csteps/inner ysep (0,) (0,0)*(0,0)(0,) (-.50,0) 5 : The local moments are used to calculate the far field contribution to the forces on the particles.
  • Figure 2: Sketch of the four different types of interactions that occur in an MD system with titratable sites. Each gray box illustrates one term of eqs \ref{['eqn:iter_ham_lambda']} and \ref{['eqn:iter_ch_ham_lambda']}, with particles as circles and interactions as lines. The first three terms (top three boxes) are calculated from scaled charges ($\tilde{q}_i$, orange circles) and are identical for Hamiltonian (HI) and charge interpolation (QI). HI differs from QI for the intra-site interactions (gray boxes in the middle), which are calculated from scaled charges for QI (left), but from pure charges (yellow and red) for HI (right). Scaled charges are obtained by weighing form 0 (yellow) and form 1 (red), as seen at the bottom.
  • Figure 3: Starting from a charge-scaled Hamiltonian, the MAHI scheme calculates the correct (periodic) $\mathcal{H}_\text{form-form}$ interactions for Hamiltonian interpolation (HI). The central magenta box shows the actual simulation volume containing a two-atomic site (black/green dots), while the surrounding boxes are periodic images. A. First, FMM calculates the interactions for the scaled charges $\tilde{q}$ using multipole expansion in the yellow areas. B. Corrections are then computed so that HI is retrieved for a site with charges $\text{Q}_{\text{P2P}}$ and $\text{Q}_\mathcal{L}$ (green). Here, in contrast to a regular FMM, all corrections to interactions coming from the first layer around the central box are computed directly (blue), while corrections from distant boxes are handled by a lattice operator (yellow).
  • Figure 4: Benzene ring solvated in water. The ball-and-stick drawing shows hydrogen atoms in white, carbon atoms in cyan, and oxygen atoms in red.
  • Figure 5: Quantification of the accuracy of the MAHI scheme. The plots show the relative deviation of the MAHI forces from the reference forces for the 1,010 particle test system with typical particle distribution (a) and (b), and with hypothetical worst-case particle distribution (c) and (d).
  • ...and 7 more figures