Theoretical and Computational Chemistry

# Scalable Constant pH Molecular Dynamics in GROMACS

## Abstract

Molecular dynamics (MD) computer simulations are used routinely to compute atomistic trajectories of complex systems. Systems are simulated in various ensembles, depending on the experimental conditions one aims to mimic. While constant energy, temperature, volume, and pressure are rather straightforward to model, pH, which is an equally important parameter in experiments, is more difficult to account for in simulations. Although a constant pH algorithm based on the $\lambda$-dynamics approach by Brooks and co-workers was implemented in a fork of the GROMACS molecular dynamics program, uptake has been rather limited, presumably due to the poor scaling of that code with respect to the number of titratable sites. To overcome this limitation, we implemented an alternative scheme for interpolating the Hamiltonians of the protonation states that makes the constant pH molecular dynamics simulations almost as fast as a normal MD simulation with GROMACS. In addition, we implemented a simpler scheme, called multisite representation, for modeling side chains with multiple titratable sites, such as imidazole rings. This scheme, which is based on constraining the sum of the $\lambda$-coordinates, not only reduces the complexity associated with parameterizing the intra-molecular interactions between the sites, but is also easily extendable to other molecules with multiple titratable sites. With the combination of a more efficient interpolation scheme and multisite representation of titratable groups, we anticipate a rapid uptake of constant pH molecular dynamics simulations within the GROMACS user community.

## Version notes

Funding information updated

## Supplementary material

Fitting code
Mathematica notebook with the implemented algorithm for fitting polynomials to gradients obtained from TI runs of reference state
Zip file with all the parameters used in the simulation
In this file, all the information needed to reproduce our simulations is provided: force field parameters, input files, etc.
Scalable constant pH in GROMACS SI
The Supporting information file, where the description of lambda-potentials, effect of neglecting the interpolation of Lennard-Jones interactions, titration results for Asp and Glu within the single-site representation, a comparison of pKa values for HEWL obtained with various lambda-dynamics-based constant pH methods, demonstration that charge interpolation requires a single evaluation of the electrostatic potential for both single- and multisite representations, and introduction into routines to fit V^MM.