# Intramolecular and water mediated
tautomerism of solvated glycine

Pengchao Zhang1,2, Axel Tosello Gardini2,3, Xuefei Xu1,*, and Michele
Parrinello2,*. 1Center for Combustion Energy, Department of Energy and Power
Engineering, and Key Laboratory for Thermal Science and Power Engineering of
Ministry of Education, Tsinghua University, Beijing 100084, China 2Atomistic
Simulations, Italian Institute of Technology, Genova 16152, Italy 3Department
of Materials Science, Università di Milano-Bicocca, 20126 Milano, Italy *Co-
corresponding author e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The understanding of prototropic tautomerism in water and the characterization
of solvent effects on protomeric equilibrium pose significant challenges.
Using molecular dynamics simulations based on state-of-the-art deep learning
potential and enhanced sampling methods, we provide a comprehensive
description of all configurational transformations in glycine solvated in
water and determine accurate free energy profiles of these processes. We
observe that the tautomerism between the neutral and zwitterionic forms of
solvated glycine can occur by both intramolecular proton transfer in glycine
and intermolecular proton transfer in the contact ion pair (anionic glycine
and hydronium ion) or the separated ion pair (cationic glycine and hydroxide
ion).

## 1 Introduction

Prototropic tautomerism plays a crucial role in determining the thermodynamic
and kinetic properties of amino acids and their derivatives, and consequently
their reactivity and interaction with other biomolecules.1, 2, 3, 4, 5, 6, 7
The simplest amino acid, glycine, can exist in the zwitterionic [Z], neutral
[N], anionic [A] or cationic [C] form (Fig. 1) depending on the pH of the
solution8, 9, 10, providing a simple yet relevant model for the study of
tautomeric behavior. Here we will focus on the transformation of glycine in
water between the [Z] and [N] forms.

Figure 1: The schematic diagrams of four glycine states and a water molecule.
[Z], [N], [A], and [C] denote zwitterionic, neutral, anionic, and cationic
forms, respectively. Atomic notation is used to distinguish the atoms present
in the molecules of glycine and water. $\mathrm{X_{y}}$ represents the atom X
(H or O) that originates from molecule y (glycine or water).

Experimental studies have been carried out to determine the protonation states
of glycine under various conditions.11, 12, 9 It has been found that in the
gas phase the [N] tautomer is the more stable. Under microsolvation
conditions,13, 14 it has been demonstrated that proton transfer can occur and
the [Z] form is favoured. For fully solvated glycine,8, 11, 9, 10, 15 the [Z]
tautomer has also been shown to be more stable than the [N] tautomer. However,
little is experimentally known on the dynamic processes of transformation from
[N] to [Z], and these microscopic processes can be investigated by theoretical
simulations at the atomic level.

The gas phase experiments have been accompanied by a number of static
calculations.16, 17, 18 However, the glycine tautomerism is difficult to
achieve in the gas phase, but it is facilitated by the dynamics of the water
environment and the fluctuation of the hydrogen bond (HB) network. In an
effort to include water effects, various approaches have been taken, from
describing water as polarizable continuum models (PCMs),19, 20, 21, 22, 9, 10,
23 to quantum mechanical/molecular mechanical (QM/MM) calculations,24, 23 and
to the use of reactive force fields.25 Most studies have come to the same
conclusion as experiments, i.e., the [Z] tautomer is stable in water. However,
only a few simulations have explored the tautomeric dynamics of glycine, and
revealed the mechanism of short-range proton transfer, i.e., the transition
from [N] to [Z] occurs intramolecularly or sometimes via a nearby solvent
water molecule.25, 24, 18 Here we aim at describing the static and dynamic
behavior of glycine in water taking into account the full complexity of water
dynamics and chemistry. The most appropriate theoretical framework for such a
study is that of ab initio molecular dynamics (AIMD).26, 27, 28, 29, 18

Unfortunately, AIMD simulations are computationally expensive, and this has
limited the size of the system studied and the time scale of the simulations.
In pioneering work, Behler and Parrinello have shown that the cost of ab
initio-quality simulations can be greatly reduced without compromising
accuracy.30 Their solution was to represent the potential energy surface as a
suitably designed neural network whose inputs are a set of descriptors chosen
so as to best describe the local atomic environment, while enforcing the
symmetry of the problem and allowing scaling up to large systems. The neural
network is then trained on a set of DFT (density functional theory) energies
and forces performed on a carefully selected set of configurations. More
recently, many variants of this approach have been proposed, such as Gaussian
approximation potential,31 deep potential (DP),32, 33, 34 SchNet,35 and other
equivariant approaches.36, 37, 38 These developments benefit greatly from easy
access to efficient machine learning libraries and from efficient neural
network training strategies. Here, we use the DP model to perform molecular
dynamics (MD) simulations, which has proven its usefulness in many
applications, including the studies of water and reactions in water.39, 40,
41, 42, 43, 44

During the DP model training, it is important to present the network with an
appropriate set of configurations, especially when we want to study reactive
processes such as the protonation and deprotonation of glycine. Since these
processes take place on a time scale that is much larger than the one
accessible to standard MD simulations, the use of enhanced sampling methods is
necessary to sample such rare events. Here we employ the recently developed
on-the-fly probability enhanced sampling (OPES) method,45, 46 which is an
evolution of metadynamics.47, 48 Like many other enhanced sampling methods,
OPES uses collective variables (CVs) to accelerate the fluctuations that
eventually lead to the reactive process. In this work we adapt the Voronoi CVs
introduced by Grifoni et al.49, 50 so as to be able to simultaneously describe
the [Z], [N], [A], [C] forms of glycine and the intermediate structures
connecting these forms.

Then we follow the deep Kohn-Sham (DeePKS)51, 52, 53, 54 strategy to
efficiently build a very accurate training dataset. In particular, we want to
perform simulations in which the energies and forces are as accurate as those
predicted by the hybrid meta-GGA functional M06-2X, which has been shown to
give an accurate representation of the thermochemistry and kinetics of main
group elements.55 However, the computational cost of performing single point
calculations at the M06-2X level for the current system is too high. In simple
terms, the way DeePKS gets around this difficulty is to train a neural network
that expresses the difference $E_{\delta}$ between the target energy and the
baseline energy. In our case the target energy is $E_{\mathrm{M06-2X}}$, and
the baseline energy $E_{\mathrm{PBE}}$ is from the cheaper GGA functional
PBE.56 The $E_{\delta}$ is calculated as the sum of atomic energy differences,
determined using a neural network that takes atomic coordinates, density
matrices, and orbitals as inputs. The remarkable finding is that to train
$E_{\delta}$ only a very small number of training data are needed. Once the
PBE-based DeePKS model is trained, energy is calculated from
$E_{\mathrm{PBE}}+E_{\delta}$, and forces can be driven from the energy. The
energies and forces generated using such a DeePKS model will have an accuracy
close to that of the full M06-2X calculation at a much lower cost, and will be
used to train the DP model.

We use the generated DP model to explore the tautomeric free energy landscape
of glycine in water, and find that transitions between the [N] and [Z] forms
of glycine can occur via multiple proton transfer (PT) processes.

## 2 Computational methods

### 2.1 Potential model generation

#### 2.1.1 Building training datasets

Key to the potential model training is the construction of appropriate
training datasets. As discussed in the Introduction section, we need to train
two potential models: One is the DP model to be used in molecular dynamics
simulations for investigating the glycine tautomerism in water; the other is
the PBE-based DeePKS model with the accuracy approaching the M06-2X functional
to be used in the labeling of the DP training datasets. In either case we use
pure water and solvated glycine configurations (see Table S1 in the SI). These
configurations are generated in multiple independent MD simulations with
enhanced sampling to generate as much uncorrelated data as possible and to
include reactive processes. The details of the MD simulations are given below.

The numbers of configurations in the training sets of the two models are much
different due to the difference between the descriptors. The similarities are
that the two model descriptors all satisfy the physical symmetry and the
locality. The differences are that the neural network descriptors in the DP
model contain only the angular and radial atomic environment, while the
descriptors of the DeePKS model additionally include density matrices
projected on the atomic orbitals and satisfy the gauge invariance symmetry33,
52, 54. This means that only a few hundred configurations are needed to train
the DeePKS model. In contrast, more than tens of thousands of configurations
are needed for training the DP model to have an accurate modeling of the
reactive process. Finally, a total of 300 and 55,498 configurations are
collected in the DeePKS and DP model training datasets, respectively.

#### 2.1.2 DFT calculations

To prepare DeePKS training datasets and all test datasets, the M06-2X55
energies and forces are calculated using the CP2K57 package. In the
calculations, the Goedecker-Teter-Hutter pseudopotentials58, 59 are used
together with a quadruple-zeta valence basis set with polarization60 (QZV3P).
The multi-grid level utilizes a plane-wave cutoff of 1000 Ry for the total
density and 70 Ry for the Kohn-Sham orbitals. To accelerate the convergence of
self-consistent field (SCF) iterations, the auxiliary density matrix method61
is utilized.

#### 2.1.3 DeePKS model training

The DeePKS model is then generated using an iterative approach, in which we
alternately train the network of the correction term from the PBE baseline to
the M06-2X target using the DeePKS-kit 51, 52, 53 package, and solve the
resulting DeePKS model in ABACUS62, 63, 54. In the ABACUS calculations, the
optimized norm-conserving Vanderbilt64 pseudopotentials are used together with
numerical atomic orbital65 basis. The kinetic energy cutoff is set at $100$
Ry, and the SCF convergence threshold for the density error is $1\times
10^{-7}$ Ry.

The above iterations are no longer needed once the DeePKS model has been
trained. One performs single point calculations based on the PBE functional
corrected with the $E_{\delta}$ term to obtain total energies and forces.
Relative to an independent set of testing data obtained by performing ordinary
M06-2X/QZV3P calculations, the PBE-based DeePKS model shows root mean square
errors (RMSEs) of 0.53 and 0.61 meV/atom for the energies and of 43 and 52
meV/Å for the forces, where we distinguish between errors relative to pure
water and solvated glycine systems. In terms of efficiency, the PBE-based
DeePKS model saves about an order of magnitude in time compared to using the
standard M06-2X functional (Table S3 in the SI).

#### 2.1.4 DP model training

The DP model is then trained in the DeePMD-kit34, 66 package using the
training dataset of M06-2X quality generated as described above. As discussed
earlier, at this stage we needed 55,498 configurations for an accurate result.
Then, the trained DP model is examined by testing the prediction on an
independently generated dataset containing 6,020 configurations for which
standard M06-2X energy and force calculations have been performed. The DP
model exhibits energy RMSEs of 0.79 and 1.72 meV/atom and atomic force RMSEs
of 58 and 63 meV/Å for water and solvated glycine systems, respectively,
suggesting that the DP model has achieved a level of precision comparable to
M06-2X. Compared to an M06-2X functional based simulation, the DP model is
five orders of magnitude faster (Table S3).

### 2.2 Molecular dynamics simulations

#### 2.2.1 MD simulation for the training dataset generation

The configurations required to train the DeePKS and DP models are sampled from
MD simulations using CP2K57 as the driver. To speed up the sampling, we use
the semi-empirical GFN1-xTB67 to derive potential energy surface. This
GFN1-xTB method has a lower computational cost, but is still accurate enough
to describe the structures of different glycine forms in water. The
simulations are performed in constant volume constant temperature (NVT)
ensemble with a time step of 1.0 fs. The temperature of 300 K is enforced
using the velocity rescaling thermostat68 with a damping time of 0.2 ps.

#### 2.2.2 DPMD simulation for studying glycine tautomerism

To study the glycine tautomerism in water, deep potential molecular dynamics
(DPMD) simulations are performed using the LAMMPS 69 software. Before the
formal simulation, we have investigated the size effect by placing a glycine
molecule in varying amounts of water. Based on the convergence test result
(Fig. S9), the system C (Table S1) with one glycine molecule in 128 H2O
molecules is used in the final DPMD simulation to strike a balance between
efficiency and precision. In order to sufficiently sample all possible
tautomeric processes, simulations are run for 30 ns in the NVT ensemble with
an integration time step of 1.0 fs. The temperature is controlled at 300 K
using again the velocity rescaling thermostat with a damping time of 0.04 ps.

### 2.3 Enhanced sampling settings

#### 2.3.1 Bias potential

The OPES method using collective variables (CVs) is employed in both MD
simulations for constructing the datasets and the DPMD simulation of glycine
tautomerism. The CVs $\mathbf{{s}({R})}$, which are functions of atomic
coordinates ($\mathbf{{R}}$), skilfully capture the slow modes associated with
rare events. During MD simulations, the potential energy ($U(\mathbf{R})$) of
the system is modified by adding the external bias potential ($V(\mathbf{s})$)
through the utilization of the PLUMED70 plugin. Specifically, $V(\mathbf{s})$
at the $n$th step in OPES45, 46 is defined by:

$V_{n}(\mathbf{s})=(1-{\frac{1}{\gamma}})\frac{1}{\beta}\log(\frac{{P}_{n}{(\mathbf{s})}}{Z_{n}}+\epsilon)$
(1)

where $\beta=1/{k_{B}T}$ is the inverse Boltzmann factor. Probability
$P(\mathbf{s})$ is the unbiased marginal distribution, and parameter $Z$ is
the normalization factor. Bias factor $\gamma=\beta\Delta E_{\mathrm{bias}}$
ensures the reshaping of the original probability. The regularization term
$\epsilon=e^{-\gamma/(1-1/\gamma)}$ not only guarantees a positive argument
for the logarithm, but also imposes a bias constraint, thereby limiting
sampling within the specified region of interest. In particular, the bias
update is carried out every 500 steps with the adaptive kernel width, and the
value of $\Delta E_{\mathrm{bias}}$ is set at 35 kJ/mol.

Notably, the "explore" variant of OPES (OPES-explore) is used in the
configuration collection of datasets due to its ability to accelerate the
exploration of phase space, while the OPES is performed to generate well-
converged free energies during the final DPMD simulation. The reason is that
the ways of estimating the probability distribution are different in the OPES
and OPES-explore methods, although the idea of defining the bias potential is
similar. Specifically, unbiased probability $P(\mathbf{s})$ is estimated on-
the-fly using weighted kernel density estimation (KDE) in OPES, while the
well-tempered probability
$p^{\mathrm{WT}}(\mathbf{s})(\propto[P(\mathbf{s})]^{\frac{1}{\gamma}}$) is
estimated based on averaged KDE in OPES-explore.45, 46

#### 2.3.2 Collective variables

Proton can diffuse through water via the Grotthuss mechanism,71, 72, 44 in
which it is not a well-specified proton that moves, but rather a charge defect
that migrates through the water hydrogen network. Thus a CV that is used to
describe proton diffusion has to be able to identify charge defects without
making reference to a specific set of atomic coordinates. In the present case,
several charge defects are possible: the –NH3+ and –COO- groups in glycine,
and the two water self-ions (hydronium H3O+ and hydroxide OH-). To identify
automatically these charge defects we follow the strategy of References49, 50
and tessellate the space with Voronoi polyhedra centered on the O and N atoms.
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_{w}}$ and $\mathrm{O_{g}}$ in Fig. 1).

To count the number of H atoms $n_{i}$ centered on O or N atom $i$, we use the
formula

$n_{i}=\sum_{j=1}^{\mathrm{Num_{H}}}\frac{e^{-\lambda|\mathbf{R}_{i}-\mathbf{R}_{j}|}}{\sum_{m=1}^{\mathrm{Num_{O\&N}}}e^{-\lambda|\mathbf{R}_{m}-\mathbf{R}_{j}|}}$
(2)

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}$

$\delta_{i}=n_{i}-n_{i}^{0}$ (3)

where we take $n_{i}^{0}=2$ for the water oxygen ($\mathrm{O_{w}}$) and the
glycine nitrogen (N), while for the carboxylic oxygens ($\mathrm{O_{g}}$) we
set $n_{i}^{0}=\frac{1}{2}$ on account of the symmetry between two oxygen
atoms.

As the system is neutral overall, the [C] and [A] forms of glycine are
compensated by OH- and H3O+, respectively. Then, to distinguish between the
[C]$-$OH- pair, [Z]$\&$[N] states and [A]$-$H3O+ pair, we define a CV
$\mathbf{s}_{p}$ that can identify these protonation states, by combining the
charge defects of both glycine and water.

$\mathbf{s}_{p}=\sum_{i=1}^{\mathrm{Num_{O_{w}}}}\delta_{i}+2(\sum_{j=1}^{\mathrm{Num_{N}}}\delta_{j}+\sum_{k=1}^{\mathrm{Num_{O_{g}}}}\delta_{k})$
(4)

where $\mathbf{s}_{p}\approx 1$ for [C]$-$OH-, $\mathbf{s}_{p}\approx-1$ for
[A]$-$H3O+, 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_{w}}$ and N and between $\mathrm{O_{w}}$
and $\mathrm{O_{g}}$, as well as the average distance between N and
$\mathrm{O_{g}}$, as follows:

$\begin{split}\mathbf{s}_{d}=&-\sum_{i=1}^{\mathrm{Num_{O_{w}}}}\sum_{j=1}^{\mathrm{Num_{N}}}r_{i,j}\delta_{i}\delta_{j}-\sum_{i=1}^{\mathrm{Num_{O_{w}}}}\sum_{k=1}^{\mathrm{Num_{O_{g}}}}r_{i,k}\delta_{i}\delta_{k}\\\
&-\sum_{j=1}^{\mathrm{Num_{N}}}\sum_{k=1}^{\mathrm{Num_{O_{g}}}}r_{j,k}\delta_{j}\delta_{k}\end{split}$
(5)

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 ~3 Å in the other two cases.

During the OPES simulations, both CVs $\mathbf{s}_{p}$ and $\mathbf{s}_{d}$
use the value of $\lambda=5$ so as to have a smoother definition of the
Voronoi polyhedra.

Since H3O+ and OH- can be present in water, we have added configurations
related to the autoionization process of water into our training datasets. In
the potential training process we use the CV

$\mathbf{s}_{a}=\sum_{i=1}^{\mathrm{Num_{O_{w}}}}\delta_{i}^{2}\ $ (6)

to promote the self ionization processes, and $\mathbf{s}_{a}$ varies from 0
to 2.

The CV $\mathbf{s}_{a}$ is supplemented by a CV $\mathbf{s}_{t}$ that
represents the distance between H3O+ and OH-:

$\mathbf{s}_{t}=-\sum_{i=1}^{\mathrm{Num_{O_{w}}}}\sum_{j>i}^{\mathrm{Num_{O_{w}}}}r_{i,j}\delta_{i}\delta_{j}\
$ (7)

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 hardly 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:

$\mathbf{s}^{\prime}_{t}=\begin{cases}\log(\mathbf{s}_{t}+\epsilon),&\
0\leq\mathbf{s}_{t}<1\\\ \mathbf{s}_{t}-1+\log(1+\epsilon),&\
\mathbf{s}_{t}\geq 1\end{cases}$ (8)

where $\epsilon=0.03$ is a regularization parameter.

#### 2.3.3 Free energy calculation

After the DPMD simulation, the free energy surface (FES) along a given CV can
be calculated as follows:

$F(\mathbf{s})=-\frac{1}{\beta}\log P(\mathbf{s})$ (9)

In the regime where the bias is quasi-static, $P(\mathbf{s})$ can be
reweighted45 as an average over the biased ensemble.

$P(\mathbf{s})=\frac{\langle\delta[\mathbf{s}-\mathbf{s}(\mathbf{R})]e^{\beta
V(\mathbf{s})}\rangle_{V}}{\langle e^{\beta V(\mathbf{s})}\rangle_{V}}$ (10)

See the supporting information for the full simulation workflow and more
computational details.

## 3 Results and discussion

### 3.1 Free energy surfaces

Figure 2: Free energy profiles. a. Two-dimensional FES as a function of the
glycine protonation coordinate $\mathbf{s}_{p}$ and the charge-charge distance
$\mathbf{s}_{d}$. b and c. One-dimentional projections of the FES along the
CVs $\mathbf{s}_{p}$ and $\mathbf{s}_{d}$. The standard deviations of one-
dimentional FESs are represented by transparent colors.

The converged FES is plotted in Fig. 2a as a function of the two CVs
$\mathbf{s}_{p}$ and $\mathbf{s}_{d}$, which reflect the glycine protonation
state and the charge-charge distance, respectively. The minima corresponding
to state [N] and state [Z] are easily identified. Less evident is the presence
of two other metastable states, [A] and [C]. Their existence can be made more
clear if we project the free energy along the CV $\mathbf{s}_{p}$ as shown in
Fig. 2b. In this one-dimensional representation, the [N] and [C] forms cannot
be resolved and are part of one single central minimum, but clearly two local
minima that correspond to [A] and [C] can be detected, and these states lie
higher in energy relative to the minimum by 44.1 kJ/mol ([A]) and 49.2 kJ/mol
([C]), respectively. To facilitate reading the result in a one dimensional
projection, we also plot the FES along the CV $\mathbf{s}_{d}$ (Fig. 2c). In
this projection, [A] and [C] cannot be identified, while the [Z] and [N]
states are now clearly visible. The energy difference between the lowest free
energy [Z] state and the [N] state is 32.6 kJ/mol, in good agreement with
experiments, which have reported values ranging from 30.4 kJ/mol to 32.1
kJ/mol.73, 74, 75

It is interesting to analyze in some detail the nature of the [A] and [C]
states. Because of the requirement that the system must be neutral, the
charged glycine protomers in these two states are accompanied by a counterion,
which is H3O+ in [A] and OH- in [C]. The behavior of these two ion pairs is
different; in [A] the hydronium remains close to the glycine, while in [C] the
ion pair can separate more easily. This different behavior is reflected in the
different $\mathbf{s}_{d}$ distributions. The $\mathbf{s}_{d}$ distribution in
[C]$-$OH- pair is much broader than that of [A]$-$H3O+ pair, and as a
consequence the FES exhibits a slight left-right asymmetry of Fig. 2a.

### 3.2 Prototropic tautomerism pathways

Figure 3: Three tautomeric pathways. a. Continuous sampling processes of
tautomeric pathways with the CVs $\mathbf{s}_{p}$ and $\mathbf{s}_{d}$, where
the color bar represents the relative simulation time. b. Trajectories of
[N]$-$[Z] tautomerism along the Intra-PT, H3O+-PT, and OH--PT pathways. The
molecular configurations involved in the reaction pathways are highlighted in
the sphere models, and other surrounding water molecules are shown in the
transparent stick models. c. The distribution of the CV $\mathbf{s}_{d}$,
i.e., intermolecular charge-charge distance (no intramolecular distance in the
[A] and [C] cases). d. The hydrogen-bond schematic diagrams of H3O+ and OH-
ions near [A] and [C].

We now analyze the pathways leading from [N] to [Z]. We have identified three
possible proton transfer pathways as shown in Fig. 3a,b. The most obvious is a
direct intramolecular proton transfer (Intra-PT), which is consistent with the
gas phase studies 76, 77, 78. The –NH2 group first rotates so that its lone
pair helps to accept the –COOH proton to form the ammonium group –NH3+, which
later rotates to assume a more stable [Z] conformation.

In addition to the direct route, there are two other pathways in which the
proton is transferred via a Grotthuss-like mechanism71, 72, 44 with the aid of
the water molecules in the solvent. These two mechanisms have either [A] or
[C] as an intermediate; in the first case, a contact ion pair
$[\mathrm{A}]-\mathrm{H}_{3}\mathrm{O}^{+}$ is formed, whereas in the second
one passes through the separated ion pair $[\mathrm{C}]-\mathrm{OH}^{-}$.

In the anionic pathway, the carboxylic proton is transferred to a nearby
solvent water, which transfers a positive charge to another water molecule via
a Zundel intermediate to end up on the glycine, forming the ammonium group of
the [Z] state. This $\mathrm{H_{3}O^{+}}$ mediated proton transfer
($\mathrm{H_{3}O^{+}}$-PT) process is also found at the microhydration limit
of glycine-water clusters.76, 18, 77

Similar but different is the new pathway that leads from [N] to [Z] via an
intermediate [C]$-$OH- pair mediated proton transfer (OH--PT). The first step
is the transfer of a proton from a water molecule to the amine group of
glycine to form an ammonium cation and an OH- ion. The OH- ion is then
solvated in water and diffuses a relatively long distance via a Grotthuss
mechanism, eventually ending up back at the glycine to abstract the carboxylic
proton to form the [Z] protomer and a water molecule.

The reason for the distinct charge-charge distance behavior comes from the
different amphipathy of $\mathrm{H_{3}O^{+}}$ and $\mathrm{OH^{-}}$ ions as
depicted in Fig. 3c,d, where the distance between the contact ion pair
$[\mathrm{A}]-\mathrm{H}_{3}\mathrm{O}^{+}$ is in the short range (< ~4.5 Å),
while the distance between the separated ion pair
$[\mathrm{C}]-\mathrm{OH}^{-}$ is in the relatively long range (> ~4.5 Å).
Since glycine disrupts the HB network of water and provides a hydrophobic
environment, $\mathrm{H_{3}O^{+}}$ tends to stay in proximity to glycine due
to its hydrophobic O atom. Compared to $\mathrm{H_{3}O^{+}}$, the O atom in
$\mathrm{OH^{-}}$ is hydrophilic and can easily form HB with surrounding water
molecules as HB acceptors,79, 44 resulting in a more extensive HB network
around $\mathrm{OH^{-}}$. This indicates that $\mathrm{OH^{-}}$ can diffuse
further away from glycine towards the outer water solvation shell compared to
$\mathrm{H_{3}O^{+}}$, and thus the outer solvation shell of glycine is
involved in the prototropic tautomerism and cannot be ignored.

Note that experimental investigations75 on the glycine protonation reactions
were mostly conducted in the water with glycine concentrations falling within
the interval of 0.03-0.25 M at room temperature. The present simulation
results are obtained based on a neutral aqueous simulation system with a
little bit higher glycine concentration of 0.43 M. Nevertheless, as documented
in the literature, the solubility of glycine in water spans a range of
0.36-0.56 M,80 and therefore the present simulation ensures a solvated glycine
molecule. In addition, it is noteworthy that the discovery of reaction
pathways involving ion pairs in the simulation, accompanied by the presence of
$\mathrm{H_{3}O^{+}}$ or $\mathrm{OH^{-}}$ in the aqueous environment,
indicates an instantaneous change of the pH, leading to values of 0.4 or 13.6,
respectively. Under the present simulated condition, we have clearly
delineated three distinct [N]$-$[Z] tautomerization pathways for the first
time, elucidated the roles and characters of [C]$-$OH- and [A]$-$H3O+ ion
pairs in these processes, and indicated a way by which the pH of the solution
can influence the protomeric equilibrium. This potentially have broader
applicability in chemistry, biology, engineering and other proton transfer
processes. We look forward to the upcoming advanced experiments that will
serve to validate our findings. It is also important to emphasize that the
quantitative results of these experiments will depend on the glycine
concentration and, more interestingly, the pH of the system.

## 4 Conclusion

In summary, the present study provides novel insights into the prototropic
tautomerism of glycine in water and covers all possible configurational
transformations using an accurate description of the interaction potential and
a thorough sampling of the potential energy surface. We discover three
pathways for tautomerization between the neutral and zwitterionic forms of
solvated glycine; one is via intramolecular proton transfer in glycine, the
second one involves short-range intermolecular proton transfer in the contact
ion pair between anionic glycine and hydronium ion, and the third one has the
aid of long-range intermolecular proton transfer in the separated ion pair
between cationic glycine and hydroxide ion. In the two intermolecular proton
transfer pathways, the observed remarkably distinct charge-charge distance of
the intermediate ion pairs is attributed to the different amphipathy of water
self-ions.

The combination of our computational technologies, including DeePKS, DeePMD,
OPES, and Voronoi CVs, not only deepens our understanding of glycine
tautomerism in water, but also provides a comprehensive framework and
methodology for facilitating further research into the intricate dynamics of
proton transfer.

## Acknowledgments

This work received partial support from the National Natural Science
Foundation of China under Grant No. [21973053]. The authors would like to
express their gratitude to Umberto Raucci, Enrico Trizio, Sudip Das, Linfeng
Zhang, Qi Ou, Andrea Rizzi and Francesco Mambretti for their valuable
discussions. Computational resources were provided by the High Performance
Computing (HPC) platform at Tsinghua University and the HPC Franklin at
Fondazione Istituto Italiano di Tecnologia.

## Author contributions

All authors contributed to the conception and design of the study. P.Z. and
A.T.G. conducted the testing of the enhanced sampling method. P.Z. performed
the simulation workflow and data analysis, with the assistance of all authors.
X.X. and M.P. provided supervision throughout stages of the project. The
initial draft of the manuscript was prepared by P.Z., and all authors
participated in editing and reviewing the final version.

## References

  * Elguero et al. 2000 Elguero, J.; Katritzky, A. R.; Denisko, O. V. Prototropic tautomerism of heterocycles: Heteroaromatic tautomerism-General overview and methodology. _Advances in Heterocyclic Chemistry_ 2000, _76_ , P1–P84
  * Singh et al. 2015 Singh, V.; Fedeles, B. I.; Essigmann, J. M. Role of tautomerism in RNA biochemistry. _Rna_ 2015, _21_ , 1–13
  * Pospisil et al. 2003 Pospisil, P.; Ballmer, P.; Scapozza, L.; Folkers, G. Tautomerism in computer-aided drug design. _Journal of Receptors and Signal Transduction_ 2003, _23_ , 361–371
  * Raczyńska et al. 2005 Raczyńska, E. D.; Kosińska, W.; Ośmiałowski, B.; Gawinecki, R. Tautomeric equilibria in relation to pi-electron delocalization. _Chemical reviews_ 2005, _105_ , 3561–3612
  * Douhal et al. 1995 Douhal, A.; Kim, S.; Zewail, A. Femtosecond molecular dynamics of tautomerization in model base pairs. _Nature_ 1995, _378_ , 260–263
  * Rodriquez et al. 2001 Rodriquez, C. F.; Cunje, A.; Shoeib, T.; Chu, I. K.; Hopkinson, A. C.; Siu, K. M. Proton migration and tautomerism in protonated triglycine. _Journal of the American Chemical Society_ 2001, _123_ , 3006–3012
  * Nemeria et al. 2009 Nemeria, N. S.; Chakraborty, S.; Balakrishnan, A.; Jordan, F. Reaction mechanisms of thiamin diphosphate enzymes: defining states of ionization and tautomerization of the cofactor at individual steps. _The FEBS journal_ 2009, _276_ , 2432–3446
  * Sheinblatt and Gutowsky 1964 Sheinblatt, M.; Gutowsky, H. A nuclear magnetic resonance study of the protolysis kinetics of glycine. _Journal of the American Chemical Society_ 1964, _86_ , 4814–4820
  * Ottosson et al. 2011 Ottosson, N.; Børve, K. J.; Spangberg, D.; Bergersen, H.; Sæthre, L. J.; Faubel, M.; Pokapanich, W.; Öhrwall, G.; Björneholm, O.; Winter, B. On the Origins of Core- Electron Chemical Shifts of Small Biomolecules in Aqueous Solution: Insights from Photoemission and ab initio Calculations of Glycineaq. _Journal of the American Chemical Society_ 2011, _133_ , 3120–3130
  * Hernández et al. 2013 Hernández, B.; Pflüger, F.; Kruglik, S. G.; Ghomi, M. Protonation–deprotonation of the glycine backbone as followed by Raman scattering and multiconformational analysis. _Chemical Physics_ 2013, _425_ , 104–113
  * Locke and McIver Jr 1983 Locke, M. J.; McIver Jr, R. T. Effect of solvation on the acid/base properties of glycine. _Journal of the American Chemical Society_ 1983, _105_ , 4226–4232
  * Kumar et al. 2005 Kumar, S.; Rai, A. K.; Singh, V.; Rai, S. Vibrational spectrum of glycine molecule. _Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy_ 2005, _61_ , 2741–2746
  * Alonso et al. 2006 Alonso, J. L.; Cocinero, E. J.; Lesarri, A.; Sanz, M. E.; López, J. C. The glycine–water complex. _Angewandte Chemie_ 2006, _118_ , 3551–3554
  * Schwaab et al. 2022 Schwaab, G.; de Tudela, R. P.; Mani, D.; Pal, N.; Roy, T. K.; Gabas, F.; Conte, R.; Caballero, L. D.; Ceotto, M.; Marx, D.; others Zwitter ionization of glycine at outer space conditions due to microhydration by six water molecules. _Physical Review Letters_ 2022, _128_ , 033001
  * Sun et al. 2014 Sun, J.; Niehues, G.; Forbert, H.; Decka, D.; Schwaab, G.; Marx, D.; Havenith, M. Understanding THz spectra of aqueous solutions: Glycine in light and heavy water. _Journal of the American Chemical Society_ 2014, _136_ , 5031–5038
  * Pérez de Tudela and Marx 2016 Pérez de Tudela, R.; Marx, D. Water-induced zwitterionization of Glycine: stabilization mechanism and spectral signatures. _The journal of physical chemistry letters_ 2016, _7_ , 5137–5142
  * Bachrach 2008 Bachrach, S. M. Microsolvation of glycine: a DFT study. _The Journal of Physical Chemistry A_ 2008, _112_ , 3722–3730
  * Tripathi et al. 2021 Tripathi, R.; Durán Caballero, L.; Pérez de Tudela, R.; Hölzl, C.; Marx, D. Unveiling Zwitterionization of Glycine in the Microhydration Limit. _ACS omega_ 2021, _6_ , 12676–12683
  * Wood et al. 2008 Wood, G. P.; Gordon, M. S.; Radom, L.; Smith, D. M. Nature of Glycine and Its $\alpha$-Carbon Radical in Aqueous Solution: A Theoretical Investigation. _Journal of Chemical Theory and Computation_ 2008, _4_ , 1788–1794
  * Senn et al. 2003 Senn, H. M.; Margl, P. M.; Schmid, R.; Ziegler, T.; Blöchl, P. E. Ab initio molecular dynamics with a continuum solvation model. _The Journal of chemical physics_ 2003, _118_ , 1089–1100
  * Aikens and Gordon 2006 Aikens, C. M.; Gordon, M. S. Incremental solvation of nonionized and zwitterionic glycine. _Journal of the American Chemical Society_ 2006, _128_ , 12835–12850
  * Kayi et al. 2012 Kayi, H.; Kaiser, R. I.; Head, J. D. A theoretical investigation of the relative stability of hydrated glycine and methylcarbamic acid—from water clusters to interstellar ices. _Physical Chemistry Chemical Physics_ 2012, _14_ , 4942–4958
  * Valverde et al. 2018 Valverde, D.; da Costa Ludwig, Z. M.; Da Costa, C. R.; Ludwig, V.; Georg, H. C. Zwitterionization of glycine in water environment: Stabilization mechanism and NMR spectral signatures. _The Journal of Chemical Physics_ 2018, _148_ , 024305
  * Choi et al. 2012 Choi, C. H.; Re, S.; Feig, M.; Sugita, Y. Quantum mechanical/effective fragment potential molecular dynamics (QM/EFP-MD) study on intra-molecular proton transfer of glycine in water. _Chemical Physics Letters_ 2012, _539_ , 218–221
  * Rahaman et al. 2011 Rahaman, O.; Van Duin, A. C.; Goddard III, W. A.; Doren, D. J. Development of a ReaxFF reactive force field for glycine and application to solvent effect and tautomerization. _The Journal of Physical Chemistry B_ 2011, _115_ , 249–261
  * Car and Parrinello 1985 Car, R.; Parrinello, M. Unified approach for molecular dynamics and density-functional theory. _Physical review letters_ 1985, _55_ , 2471
  * Born and Heisenberg 1985 Born, M.; Heisenberg, W. Zur quantentheorie der molekeln. _Original Scientific Papers Wissenschaftliche Originalarbeiten_ 1985, 216–246
  * Sun et al. 2010 Sun, J.; Bousquet, D.; Forbert, H.; Marx, D. Glycine in aqueous solution: solvation shells, interfacial water, and vibrational spectroscopy from ab initio molecular dynamics. _The Journal of chemical physics_ 2010, _133_ , 09B609
  * Leung and Rempe 2005 Leung, K.; Rempe, S. B. Ab initio molecular dynamics study of glycine intramolecular proton transfer in water. _The Journal of chemical physics_ 2005, _122_ , 184506
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. _Physical review letters_ 2007, _98_ , 146401
  * Bartók et al. 2010 Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. _Physical review letters_ 2010, _104_ , 136403
  * Zhang et al. 2018 Zhang, L.; Han, J.; Wang, H.; Saidi, W.; Car, R.; others End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems. _Advances in neural information processing systems_ 2018, _31_
  * Zhang et al. 2018 Zhang, L.; Han, J.; Wang, H.; Car, R.; Weinan, E. Deep potential molecular dynamics: a scalable model with the accuracy of quantum mechanics. _Physical review letters_ 2018, _120_ , 143001
  * Wang et al. 2018 Wang, H.; Zhang, L.; Han, J.; E, W. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. _Computer Physics Communications_ 2018, _228_ , 178–184
  * Schütt et al. 2018 Schütt, K. T.; Sauceda, H. E.; Kindermans, P.-J.; Tkatchenko, A.; Müller, K.-R. Schnet–a deep learning architecture for molecules and materials. _The Journal of Chemical Physics_ 2018, _148_ , 241722
  * Schütt et al. 2021 Schütt, K.; Unke, O.; Gastegger, M. Equivariant message passing for the prediction of tensorial properties and molecular spectra. International Conference on Machine Learning. 2021; pp 9377–9388
  * Batzner et al. 2022 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. _Nature communications_ 2022, _13_ , 2453
  * Hoogeboom et al. 2022 Hoogeboom, E.; Satorras, V. G.; Vignac, C.; Welling, M. Equivariant diffusion for molecule generation in 3d. International conference on machine learning. 2022; pp 8867–8887
  * Zhang et al. 2021 Zhang, L.; Wang, H.; Car, R.; Weinan, E. Phase diagram of a deep potential water model. _Physical review letters_ 2021, _126_ , 236001
  * Piaggi et al. 2022 Piaggi, P. M.; Weis, J.; Panagiotopoulos, A. Z.; Debenedetti, P. G.; Car, R. Homogeneous ice nucleation in an ab initio machine-learning model of water. _Proceedings of the National Academy of Sciences_ 2022, _119_ , e2207294119
  * de la Puente et al. 2022 de la Puente, M.; David, R.; Gomez, A.; Laage, D. Acids at the Edge: Why Nitric and Formic Acid Dissociations at Air–Water Interfaces Depend on Depth and on Interface Specific Area. _Journal of the American Chemical Society_ 2022, _144_ , 10524–10529
  * Galib and Limmer 2021 Galib, M.; Limmer, D. T. Reactive uptake of N2O5 by atmospheric aerosol is dominated by interfacial processes. _Science_ 2021, _371_ , 921–925
  * Yang et al. 2022 Yang, M.; Bonati, L.; Polino, D.; Parrinello, M. Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water. _Catalysis Today_ 2022, _387_ , 143–149
  * Zhang et al. 2023 Zhang, P.; Feng, M.; Xu, X. Double-layer distribution of hydronium and hydroxide ions in the air-water interface. _ChemRxiv_ 2023,
  * Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking metadynamics: From bias potentials to probability distributions. _The journal of physical chemistry letters_ 2020, _11_ , 2731–2736
  * Invernizzi and Parrinello 2022 Invernizzi, M.; Parrinello, M. Exploration vs convergence speed in adaptive-bias enhanced sampling. _Journal of Chemical Theory and Computation_ 2022, _18_ , 3988–3996
  * Laio and Parrinello 2002 Laio, A.; Parrinello, M. Escaping free-energy minima. _Proceedings of the national academy of sciences_ 2002, _99_ , 12562–12566
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. _Physical review letters_ 2008, _100_ , 020603
  * Grifoni et al. 2019 Grifoni, E.; Piccini, G.; Parrinello, M. Microscopic description of acid–base equilibrium. _Proceedings of the National Academy of Sciences_ 2019, _116_ , 4054–4057
  * Grifoni et al. 2020 Grifoni, E.; Piccini, G.; Parrinello, M. Tautomeric equilibrium in condensed phases. _Journal of Chemical Theory and Computation_ 2020, _16_ , 6027–6031
  * Chen et al. 2020 Chen, Y.; Zhang, L.; Wang, H.; E, W. Ground state energy functional with Hartree–Fock efficiency and chemical accuracy. _The Journal of Physical Chemistry A_ 2020, _124_ , 7155–7165
  * Chen et al. 2020 Chen, Y.; Zhang, L.; Wang, H.; E, W. DeePKS: A comprehensive data-driven approach toward chemically accurate density functional theory. _Journal of Chemical Theory and Computation_ 2020, _17_ , 170–181
  * Chen et al. 2023 Chen, Y.; Zhang, L.; Wang, H.; Weinan, E. DeePKS-kit: A package for developing machine learning-based chemically accurate energy and density functional models. _Computer Physics Communications_ 2023, _282_ , 108520
  * Li et al. 2022 Li, W.; Ou, Q.; Chen, Y.; Cao, Y.; Liu, R.; Zhang, C.; Zheng, D.; Cai, C.; Wu, X.; Wang, H.; others DeePKS+ABACUS as a Bridge between Expensive Quantum Mechanical Models and Machine Learning Potentials. _The Journal of Physical Chemistry A_ 2022, _126_ , 9154–9164
  * Zhao and Truhlar 2008 Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. _Theoretical chemistry accounts_ 2008, _120_ , 215–241
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Physical review letters_ 1996, _77_ , 3865
  * Kühne et al. 2020 Kühne, T. D.; Iannuzzi, M.; Del Ben, M.; Rybkin, V. V.; Seewald, P.; Stein, F.; Laino, T.; Khaliullin, R. Z.; Schütt, O.; Schiffmann, F.; others CP2K: An electronic structure and molecular dynamics software package-Quickstep: Efficient and accurate electronic structure calculations. _The Journal of Chemical Physics_ 2020, _152_ , 194103
  * Goedecker et al. 1996 Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. _Physical Review B_ 1996, _54_ , 1703
  * Hartwigsen et al. 1998 Hartwigsen, C.; Gœdecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. _Physical Review B_ 1998, _58_ , 3641
  * VandeVondele et al. 2005 VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. _Computer Physics Communications_ 2005, _167_ , 103–128
  * Guidon et al. 2010 Guidon, M.; Hutter, J.; VandeVondele, J. Auxiliary density matrix methods for Hartree- Fock exchange calculations. _Journal of chemical theory and computation_ 2010, _6_ , 2348–2364
  * Chen et al. 2010 Chen, M.; Guo, G.; He, L. Systematically improvable optimized atomic basis sets for ab initio calculations. _Journal of Physics: Condensed Matter_ 2010, _22_ , 445501
  * Li et al. 2016 Li, P.; Liu, X.; Chen, M.; Lin, P.; Ren, X.; Lin, L.; Yang, C.; He, L. Large-scale ab initio simulations based on systematically improvable atomic basis. _Computational Materials Science_ 2016, _112_ , 503–517
  * Schlipf and Gygi 2015 Schlipf, M.; Gygi, F. Optimization algorithm for the generation of ONCV pseudopotentials. _Computer Physics Communications_ 2015, _196_ , 36–44
  * Lin et al. 2021 Lin, P.; Ren, X.; He, L. Strategy for constructing compact numerical atomic orbital basis sets by incorporating the gradients of reference wavefunctions. _Physical Review B_ 2021, _103_ , 235131
  * Zeng et al. 2023 Zeng, J.; Zhang, D.; Lu, D.; Mo, P.; Li, Z.; Chen, Y.; Rynik, M.; Huang, L.; Li, Z.; Shi, S.; others DeePMD-kit v2: A software package for deep potential models. _The Journal of Chemical Physics_ 2023, _159_
  * Grimme et al. 2017 Grimme, S.; Bannwarth, C.; Shushkov, P. A robust and accurate tight-binding quantum chemical method for structures, vibrational frequencies, and noncovalent interactions of large molecular systems parametrized for all spd-block elements (Z= 1–86). _Journal of chemical theory and computation_ 2017, _13_ , 1989–2009
  * Bussi et al. 2007 Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. _The Journal of chemical physics_ 2007, _126_ , 014101
  * Plimpton 1995 Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. _Journal of computational physics_ 1995, _117_ , 1–19
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Computer physics communications_ 2014, _185_ , 604–613
  * von Grotthuß 1805 von Grotthuß, T. _Mémoire sur la décomposition de l’eau et des corps qu’elle tient en dissolution à l’aide de l’électricité galvanique_ ; 1805
  * Marx 2006 Marx, D. Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. _ChemPhysChem_ 2006, _7_ , 1848–1870
  * Wada et al. 1982 Wada, G.; Tamura, E.; Okina, M.; Nakamura, M. On the ratio of zwitterion form to uncharged form of glycine at equilibrium in various aqueous media. _Bulletin of the Chemical Society of Japan_ 1982, _55_ , 3064–3067
  * Haberfield 1980 Haberfield, P. What is the energy difference between $\mathrm{H_{2}NCH_{2}CO_{2}H}$ and $\mathrm{{}^{+}H_{3}NCH_{2}CO_{2}^{-}}$? _Journal of Chemical Education_ 1980, _57_ , 346
  * Slifkin and Ali 1984 Slifkin, M.; Ali, S. Thermodynamic parameters of the activation of glycine zwitterion protonation reactions. _Journal of molecular liquids_ 1984, _28_ , 215–221
  * Jensen and Gordon 1995 Jensen, J. H.; Gordon, M. S. On the number of water molecules necessary to stabilize the glycine zwitterion. _Journal of the American Chemical Society_ 1995, _117_ , 8159–8170
  * Tang et al. 2018 Tang, W.; Cai, C.; Zhao, S.; Liu, H. Development of reaction density functional theory and its application to glycine tautomerization reaction in aqueous solution. _The Journal of Physical Chemistry C_ 2018, _122_ , 20745–20754
  * Tuñón et al. 2000 Tuñón, I.; Silla, E.; Ruiz-López, M. F. On the tautomerization process of glycine in aqueous solution. _Chemical Physics Letters_ 2000, _321_ , 433–437
  * Chen et al. 2018 Chen, M.; Zheng, L.; Santra, B.; Ko, H.-Y.; DiStasio Jr, R. A.; Klein, M. L.; Car, R.; Wu, X. Hydroxide diffuses slower than hydronium in water because its solvated structure inhibits correlated proton transfer. _Nature chemistry_ 2018, _10_ , 413–419
  * Bowden et al. 2018 Bowden, N. A.; Sanders, J. P.; Bruins, M. E. Solubility of the proteinogenic $\alpha$-amino acids in water, ethanol, and ethanol–water mixtures. _Journal of Chemical & Engineering Data_ 2018, _63_ , 488–497

## Table of contents (TOC)