Scalable Constant pH Molecular Dynamics in GROMACS

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 λ-dynamics
approach by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414–2423]
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 λ-coordinates, not only reduces the complexity
associated with parametrizing the intramolecular interactions between
the sites but also is 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.