Voronoi CVs for enhanced sampling autoionization and tautomerism

💭Catalogs

• Vorinoi-like polyhedra
• Water autoionization
• Glycine tautomerism
• Derivatives
• End

Vorinoi-like polyhedra

For acid-base equilibria in water, a crucial step is the identification of charge defects (H${}_3$O${}^{+}$ and OH${}^{–}$) within the hydrogen bond network of the system. The most intuitive and naive approach to detect such charge carriers is to partition the space with spheres around every water oxygen atom and to count how many hydrogen atoms fall within each sphere. Therefore, the oxygen atoms whose spheres count an excess or a defect of hydrogen atoms are identified as the charge carriers. This measurement can be easily defined mathematically by means of sum of smooth step-like functions, also known as a switching function $$\begin{array}{c}{c}_i=\sum _{j=1}^{\mathrm{N}}\frac{1-(\frac{|{\mathbf{R}}_i-{\mathbf{R}}_j{|}^n}{R_0})}{1-(\frac{|{\mathbf{R}}_i-{\mathbf{R}}_j{|}^m}{R_0})}\end{array}$$ reports an example of such measurement by means of a rational switching function counting how many atoms in position ${\mathbf{R}}_j$ are within the sphere of the atom $i$. Here $m$ and $n$ are two coefficients that control the smoothness of the function, $j$ runs over the hydrogen atoms and $R_0$ is the cutoff radius of the sphere.

Compared to standard chemical reactions, here the elusive nature of water ions makes their description through well-defined structures not possible and their identification only via hydronium and hydroxide ions leads to errors in the way these species are counted. In order to overcome these limitations and to build a generally applicable set of CVs we have made two conceptual steps: indistinguishable acid-base sites are considered as participating as a whole to the reaction and, once tessellated the entire space in Voronoi-like polyhedra centered on these sites, different protonation states are distinguished by the number of hydrogen atoms that contain.

The first point circumvents the problem of determining well defined geometries of the charged species. Reactants and products are in this way defined according to local anomalies in the density of hydrogen atoms without fixing beforehand any reference structures. The second one ensures a very precise and unambiguous way of evaluating this density. Figure below clearly shows how this way of partitioning the space does not present any of the previous limitations.

According to Voronoi principle, the total number of hydrogen atoms assigned to a given site $i$ is taken as the fraction of protons that are much closer to that site than from all the other ones. Then, the whole space is tessellated in polyhedra centered on any atom able to break or form bonds with an acid proton and every hydrogen atom is assigned unequivocally to one and only one of these polyhedra. This provides a very robust definition for assigning hydrogen atoms to the respective acid-base site without arbitrarily defining parameters and avoiding ambiguities when tiling the space.

However, a function to be used in enhanced sampling method must be continuous and differentiable everywhere. Unfortunately, the standard Voronoi partition is described by a set of step functions, $$\begin{array}{c}{\mathrm{w}}_i({\mathbf{R}}_i)=\{\begin{array}{l}1,if|{\mathbf{R}}_i-{\mathbf{R}}_j|<|{\mathbf{R}}_m-{\mathbf{R}}_j|\forall m\ne i\\ 0,if|{\mathbf{R}}_i-{\mathbf{R}}_j|>|{\mathbf{R}}_m-{\mathbf{R}}_j|\forall m\ne i\end{array}\end{array}$$

so that ${\mathrm{w}}_i({\mathbf{R}}_j)=1$ if the site $i$ is the closest to the hydrogen $j$ and 0 otherwise. In proximity of the Voronoi boundaries the function derivatives, and therefore any force imposed by an external potential, are ill defined.

In order to reproduce the same behaviour of above equation without discontinuity or singularity, index functions ${\mathrm{w}}_i$ are replaced by soft-max functions $n_i$.

In the case of glycine tautomerism, several charge defects are possible: the -NH${}_3$${}^{+}$ and -COO${}^{-}$ groups in glycine, and the two water self-ions (hydronium H${}_3$O${}^{+}$ and hydroxide OH${}^{-}$). We then sum the charge contained in each polyhedron. If the charge in a polyhedron is different from zero, we attribute the charge defect to the atom N or O that is at the polyhedron center, and we distinguish the O defects according to whether they are centered on water or glycine oxygen (${\mathrm{O}}_{\mathrm{w}}$ and ${\mathrm{O}}_{\mathrm{g}}$).

To count the number of H atoms $n_i$ centered on O or N atom $i$, we use the formula $$\begin{array}{c}{n}_i=\sum _{j=1}^{\mathrm{Nu}{\mathrm{m}}_{\mathrm{H}}}\frac{e^{-\lambda |{\mathbf{R}}_i-{\mathbf{R}}_j|}}{{\sum }_{m=1}^{\mathrm{Nu}{\mathrm{m}}_{O\&N}}{e}^{-\lambda |{\mathbf{R}}_m-{\mathbf{R}}_j|}},\end{array}$$ where the first sum is over all H atoms with index $j$, and $m$ is the index of Voronoi centers. The parameter $\lambda$ regulates the smoothness of the function.

The charge defect number ${\delta }_i$ of the Voronoi center $i$ is calculated by subtracting the reference number $n_i^0$ (the original proton number connecting to the atom at the polyhedron center) from H number $n_i$. $$\begin{array}{c}{\delta }_i=n_i-n_i^0,\end{array}$$ where we take $n_i^0=2$ for the water oxygen (${\mathrm{O}}_{\mathrm{w}}$) and the glycine nitrogen (N), while for the carboxylic oxygens (${\mathrm{O}}_{\mathrm{g}}$) we set $n_i^0=\frac{1}{2}$ on account of the symmetry between two oxygen atoms.

Water autoionization

Since H${}_3$O${}^{+}$ and OH${}^{-}$ can be present in water, we use the CV ${\mathbf{s}}_a$ to promote the self-ionization processes, and its value varies from 0 to 2. $$\begin{array}{c}{\mathbf{s}}_a=\sum _{i=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}{\delta }_i^2.\end{array}$$

The CV ${\mathbf{s}}_a$ is supplemented by a CV ${\mathbf{s}}_t$ that represents the distance between H${}_3$O${}^{+}$ and OH${}^{-}$, $$\begin{array}{c}{\mathbf{s}}_t=-\sum _{i=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}\sum _{j>i}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}{r}_{i,j}{\delta }_i{\delta }_j.\end{array}$$

While ${\mathbf{s}}_t$ can distinguish between the pure water state (${\mathbf{s}}_t\approx 0$) and autoionization state (${\mathbf{s}}_t>0$), it is hard to identify the initial proton transfer which corresponds to a ${\mathbf{s}}_t$ value nearly concentrated around zero. Therefore, we use a piecewise logarithmic function $$\begin{array}{c}{\mathbf{s}}_t^{\prime }=\{\begin{array}{ll}\log ({\mathbf{s}}_t+ϵ),& 0\le {\mathbf{s}}_t<1,\\ {\mathbf{s}}_t-1+\log (1+ϵ),& {\mathbf{s}}_t\ge 1,\end{array}\end{array}$$ where $ϵ=0.03$ is a regularization parameter.

Glycine tautomerism

As the system is neutral overall, the [C] and [A] forms of glycine are compensated by OH${}^{-}$ and H${}_3$O${}^{+}$, respectively. Then, to distinguish between the [C]$-$OH${}^{-}$ pair, [Z]$\&$[N] states and [A]$-$H${}_3$O${}^{+}$ pair, we define a CV ${\mathbf{s}}_p$ that can identify these protonation states, by combining the charge defects of both glycine and water. $$\begin{array}{c}{\mathbf{s}}_p=\sum _{i=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}{\delta }_i+2(\sum _{j=1}^{\mathrm{Nu}{\mathrm{m}}_{\mathrm{N}}}{\delta }_j+\sum _{k=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{g}}}}{\delta }_k),\end{array}$$ where ${\mathbf{s}}_p\approx 1$ for [C]$-$OH${}^{-}$, ${\mathbf{s}}_p\approx -1$ for [A]$-$H${}_3$O${}^{+}$, and ${\mathbf{s}}_p\approx 0$ in the other two cases.

Since ${\mathbf{s}}_p$ cannot distinguish the [Z] and [N] forms, we introduce another CV ${\mathbf{s}}_d$ to estimate the charge-charge distance that is capable of separating the two cases. The ${\mathbf{s}}_d$ is measured in terms of the distances between ${\mathrm{O}}_{\mathrm{w}}$ and N and between ${\mathrm{O}}_{\mathrm{w}}$ and ${\mathrm{O}}_{\mathrm{g}}$, as well as the average distance between N and ${\mathrm{O}}_{\mathrm{g}}$, as follows $$\begin{array}{c}{\mathbf{s}}_d=-\sum _{i=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}\sum _{j=1}^{\mathrm{Nu}{\mathrm{m}}_{\mathrm{N}}}{r}_{i,j}{\delta }_i{\delta }_j-\sum _{i=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{w}}}}\sum _{k=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{g}}}}{r}_{i,k}{\delta }_i{\delta }_k-\sum _{j=1}^{\mathrm{Nu}{\mathrm{m}}_{\mathrm{N}}}\sum _{k=1}^{\mathrm{Nu}{\mathrm{m}}_{{\mathrm{O}}_{\mathrm{g}}}}{r}_{j,k}{\delta }_j{\delta }_k,\end{array}$$ where $r$ is the modulus distance between the two Voronoi centers. The CV ${\mathbf{s}}_d$ will approximately be 0 and 3 Å in the [N] and [Z] forms, respectively, and will become larger than \textasciitilde 3 Å in the other two cases.

Based on OPES with well-defined Voronoi CVs, we can obtain accurate free energy surfaces that distinguish all possible glycine configurations and their transition processes.

Derivatives

When adapting Voronoi CVs ($\mathbf{s}(\mathbf{R})$) to distinguish glycine states, it is important to note that these CVs can be derived with respect to the atomic coordinates ($\mathbf{R}$), as indicated in the following equation. $$\begin{array}{c}\frac{\partial V(\mathbf{s})}{\partial \mathbf{R}}=\frac{\partial V(\mathbf{s})}{\partial \mathbf{s}(\mathbf{R})}\frac{\partial \mathbf{s}(\mathbf{R})}{\partial \mathbf{R}},\end{array}$$ where the implementation of the first term $\frac{\partial V(\mathbf{s})}{\partial \mathbf{s}(\mathbf{R})}$ has been completed within the PLUMED plugin, we will proceed to describe the second term $\frac{\partial \mathbf{s}(\mathbf{R})}{\partial \mathbf{R}}$ in detail and make corresponding changes to the code. For clarity and ease of understanding, we will focus on the relatively complicated part of the Voronoi CVs, as represented by the following equation. If there is a need to obtain derivatives for other CVs (${\mathbf{s}}_p$, ${\mathbf{s}}_d$, ${\mathbf{s}}_a$, and ${\mathbf{s}}_t$), such calculations can be easily achieved by applying the derivative chain rule. $$\begin{array}{c}\mathbf{s}(\mathbf{R})=\sum _{i\in \text{group}}{\delta }_i.\end{array}$$ Next, we will illustrate the process using the example of distinguishing autoionization in pure water without any prior information. In this scenario, all O atoms act as Voronoi centers and form a defined group. To identify the protonation states of all O atoms, all nearby H atoms around the target O${}_i$ atom (labeled as O${}_i\to$H${}_{\mathrm{all}}$) and all nearby O atoms around the target H${}_j$ atom (labeled as O${}_i\to$H${}_j\to$O${}_{\mathrm{all}}$) are searched.

In O${}_i\to$H${}_{\mathrm{all}}$ stage, if H${}_j$ (selected in H${}_{\mathrm{all}}$) is a neighbor of a given target O${}_i$, the derivative for H${}_j$ coordinate is governed by the following equation. For simplicity, some small terms are ignored. $$\begin{array}{c}\frac{\partial \mathbf{s}(\mathbf{R})}{\partial {\mathbf{R}}_{{\mathrm{H}}_j}}\approx \sum _{i\in \mathrm{O}}\sum _{j\in \mathrm{H}}-\lambda \frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\left[1-\frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\right]\frac{{\mathbf{R}}_{{\mathrm{H}}_j}-{\mathbf{R}}_{{\mathrm{O}}_i}}{\left|{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}\right|}.\end{array}$$ If H${}_j$ (selected in H${}_{\mathrm{all}}$) is a neighbor of a non-target Voronoi center O${}_n$ ($n\ne i$), the derivative for H${}_j$ coordinate is given by the following equation. For simplicity, some small terms are ignored. $$\begin{array}{c}\frac{\partial \mathbf{s}(\mathbf{R})}{\partial {\mathbf{R}}_{{\mathrm{H}}_j}}\approx \sum _{i\in \mathrm{O}}\sum _{j\in \mathrm{H}}-\lambda \frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_n}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\left[0-\frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\right]\frac{{\mathbf{R}}_{{\mathrm{H}}_j}-{\mathbf{R}}_{{\mathrm{O}}_n}}{\left|{\mathbf{R}}_{{\mathrm{O}}_n}-{\mathbf{R}}_{{\mathrm{H}}_j}\right|}.\end{array}$$

In O${}_i\to$H${}_j\to$O${}_{\mathrm{all}}$ stage, if O${}_i$ (selected in O${}_{\mathrm{all}}$) is a neighbor of H${}_j$, the derivative for O${}_i$ coordinate is governed by the following equation $$\begin{array}{c}\frac{\partial \mathbf{s}(\mathbf{R})}{\partial {\mathbf{R}}_{{\mathrm{O}}_i}}=\sum _{i\in \mathrm{O}}\sum _{j\in \mathrm{H}}-\lambda \frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\left[1-\frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\right]\frac{{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}}{\left|{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}\right|}.\end{array}$$ If O${}_k$ (selected in O${}_{\mathrm{all}}$, $k\ne i$) is a neighbor of H${}_j$, the derivative for O${}_k$ coordinate is given by the following equation $$\begin{array}{c}\frac{\partial \mathbf{s}(\mathbf{R})}{\partial {\mathbf{R}}_{{\mathrm{O}}_k}}=\sum _{i\in \mathrm{O}}\sum _{j\in \mathrm{H}}-\lambda \frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_k}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\left[0-\frac{e^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_i}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}{{\sum }_{m\in \mathrm{O}}{e}^{-\lambda |{\mathbf{R}}_{{\mathrm{O}}_m}-{\mathbf{R}}_{{\mathrm{H}}_j}|}}\right]\frac{{\mathbf{R}}_{{\mathrm{O}}_k}-{\mathbf{R}}_{{\mathrm{H}}_j}}{\left|{\mathbf{R}}_{{\mathrm{O}}_k}-{\mathbf{R}}_{{\mathrm{H}}_j}\right|}\approx 0.\end{array}$$ At this point, all derivatives with respect to the coordinates have been determined. In addition, we incorporate the neighbor list function, which efficiently searches for nearby atoms, reducing the computational complexity from $\mathcal{O}(N[{\mathrm{O}}_{\mathrm{all}}\times {\mathrm{H}}_{\mathrm{all}}\times {\mathrm{O}}_{\mathrm{all}}])$ to $\mathcal{O}(N[{\mathrm{O}}_{\mathrm{all}}\times {\mathrm{H}}_{\mathrm{neighbor}}\times {\mathrm{O}}_{\mathrm{neighbor}}])$. In particular, a neighbor list with a cutoff $2.4$ Å is set, updated at each step, to search for atoms in groups.

End