Protonation state sampling

Strategies for constant-pH simulation for atomistic potentials without explicit bonded terms

Theory

Consider a system at equilibrium in a semigrand ensemble relevant for constant-pH simulation.

The thermodynamic parameters $\theta \equiv \{ \beta, p, \mu_H \}$ denote the external inverse thermal energy $\beta = (k_B T)^{-1}$, external pressure $p$, and external proton chemical potential $\mu_H$.

Denote the instantaneous configuration of the system by $x$, where the relevant scalar properties conjugate to the thermodynamic parameters are the potential energy $U(x)$, the instantaneous volume $V(x)$, and the instantaneous number of protons $N_H(x)$.

We can write the reduced potential $u(x; \theta)$ for this ensemble as

$$ u(x; \theta) \equiv \beta U(x) + p V(x) + \mu_H N_H(x) $$

(NOTE: The inverse thermal energy $\beta$ has been absorbed into $p$ and $\mu_H$ in this equation, which differs from how we traditionally write the reduced potential.)

The probability density of interest is given by

$$ \pi(x; \theta) = Z(\theta)^{-1} q(x; \theta) $$

where the unnormalized probability density $q(x; \theta)$ is given by the Boltzmann law

$$ q(x; \theta) \equiv e^{-u(x; \theta)} $$

and the partition function $Z(\theta)$ is given by

$$ Z(\theta) \equiv \sum_{N_H=1}^{\infty} \int dx q(x_{N_H}; \theta) $$

where we have written the explicit subscript $x_{N_H}$ to specify that there are exactly $N_H$ protons in the system.

Algorithms

For simplicity, we first consider algorithms for dealign with tautomerization. These algorithms can be easily adapted to handle protonation states in a straightforward manner.

Consider the tautomerization of imidazole

(TODO: Replace this with a better figure)

Instantaneous Monte Carlo with atom-centered Gaussians

The simplest approach we will consider is an instantaneous Metropolis Monte Carlo proposal where we propose to "hop" the proton from one titratable heavy atom site to another.

• Step 1: Propose $x_{new} \sim p(x_{new} | x_{old})$ where the proton at position 1 is moved to a position that is a Gaussian random variable about the N at position 3.
• Step 2: Accept or reject the move with the Metropolis-Hastings criteria

$$ P_{accept} = \min \{ 1, \frac{p(x_{old} | x_{new})}{p(x_{new} | x_{old})} \frac{\pi(x_{new})}{\pi(x_{old})} \} $$

We expect the average acceptance rate $< P_{accept} >$ to be low in vacuum due to both the low proposal ratio $r \equiv \frac{p(x_{old} | x_{new})}{p(x_{new} | x_{old})}$ and the low probability ratio $s = \frac{\pi(x_{new})}{\pi(x_{old})}$ for this algorithm. In explicit solvent, the probability ratio $s$ may be especially low because protonation state

Useful experiments