.gifpng.pngconvert gif:#1 png: .gif

# Learning a reactive potential for silica-water through uncertainty
attribution

Swagata Roy  Department of Materials Science and Engineering, Massachusetts
Institute of Technology  Johannes P. Dürholt  Evonik Operations GmbH  Thomas
S. Asche  Evonik Operations GmbH  Federico Zipoli  IBM Research  Rafael Gómez-
Bombarelli  Department of Materials Science and Engineering, Massachusetts
Institute of Technology  Corresponding Author. E-mail<EMAIL_ADDRESS>
###### Abstract

The reactivity of silicates in aqueous solution is relevant to various
chemistries ranging from silicate minerals in geology, to the C-S-H phase in
cement, nanoporous zeolite catalysts, or highly porous precipitated silica.
While simulations of chemical reactions can provide insight at the molecular
level, balancing accuracy and scale in reactive simulations in the condensed
phase is a challenge.

Here, we demonstrate how a machine-learning reactive interatomic potential can
accurately capture silicate-water reactivity. The model was trained on a new
dataset comprising 400,000 energies and forces of molecular clusters at the
$\omega$-B97XD
def2-TVZP level. To ensure the robustness of the model, we introduce a new and
general active learning strategy based on the attribution of the model
uncertainty, that automatically isolates uncertain regions of bulk simulations
to be calculated as small-sized clusters.

Our trained potential is found to reproduce static and dynamic properties of
liquid water and solid crystalline silicates, despite having been trained
exclusively on cluster data. Furthermore, we utilize enhanced sampling
simulations to recover the self-ionization reactivity of water accurately, and
the acidity of silicate oligomers, and lastly study the silicate dimerization
reaction in a water solution at neutral conditions and find that the reaction
occurs through a flanking mechanism.

Keywords: Neural Network potential; Chemical reactions; Active learning;
Uncertainty attribution; enhanced sampling; silica; water.

##  1 Introduction

One of the most abundant materials on the planet, silica[1] and its many
polymorphs are used in applications like catalysts[2], pharmaceuticals[3],
nanotechnology[4], and additives[5]. Amorphous silica, in particular, has
dozens of commercial applications and the global precipitated silica market is
expected to reach US$3.49 billion by 2023[6].

Precipitated silica or silica gel is produced by the polymerization of water-
soluble silicate salts in the presence of mineral acids, resulting in a
precipitate or gel[7]. Similarly, zeolites are crystallized from a precursor
gel in hydrothermal conditions[8]. While various experimental techniques can
characterize nanosized crystallites and precipitates, they cannot resolve the
individual silicate aggregates or the network of silicon connectivity [9, 10].

Molecular simulations of reactivity in the aqueous condensed phase can bridge
this knowledge gap on silica condensation. The key challenge is describing
covalent bond breaking and formation accurately while being able to simulate
large length- and time-scales. Numerous classical interatomic potentials for
anhydrous silica have been generated to simulate silica glasses or other
polymorphs[11, 12], including the well-known BKS[11] pair potential for
$\alpha$-quartz and silica glasses, and more recent reparametrizations for
nanoclusters and amorphous silica[12, 13]. These models, along with ReaxFF
models parameterized for water-silica[14, 15] cannot accurately replicate the
condensation reactions in aqueous solution. The recent advances in machine
learning (ML) techniques and the expansion of neural networks (NN) and
automatic differentiation routines[16, 17] have provided a branch between
physical sciences and statistical learning. NN-based interatomic potentials
(NNIP) built on a variety of molecular representations retain the accuracy of
the ab initio training data and can be executed at a computational cost much
lower than first principle calculations. Several NNIP have been constructed
for pure water[18, 19, 20], amorphous and liquid silica[21, 22] and siliceous
zeolites[23]. The authors of the latter study demonstrated superior
performance to ReaxFF. However, no NNIP exists that models reactive silica-
water interactions in acid or bases, which are required to understand the
fundamental steps of silica precipitation, crystallization and gel formation.

The initial step of silicate polymerization or oligomerization is dimer
formation, whose mechanism depends on pH and temperature[24, 25]. In neutral
conditions or in a slightly acidic medium, the dimerization can occur through
neutral monomers[24]. Aqueous reactions such as these are hard to understand
at an atomic level, given the need to simulate the role of explicit solvent in
the condensed phase, various bond breaking and formation, and proton-transfers
in the solvent and reactants. Quantum chemical simulations can provide
accurate energies for systems up to hundreds of atoms but are too costly to
address free energy questions, that arise in the presence of full explicit
solvation. A couple of mechanisms have been reported for the neutral
dimerization reaction[26]. One involves a $S_{N}2$-backside reaction[27]
whereas the other involves a lateral "flank-side" attack[28].

In this study, we trained a reactive water-silica NNIP through an active
learning loop. Our ground truth data consists of orbital-based hybrid DFT
calculations on molecular clusters of various sizes and the trained NNIP
generalizes to condensed phase simulation of thousands of atoms, both aqueous
and solid silica. Since production simulations contain many more atoms than we
can simulate with DFT - we introduce a new uncertainty attribution technique
that builds on adversarial attacks on uncertainty[29] to locate clusters of
uncertainty within large simulation boxes for annotating with DFT. The NNIP
was validated successfully against experimental properties of liquid water,
and silica, both physical (diffusivity, structure, vibrational) and chemical
(dissociation constants). We then used the potential to evaluate the
dimerization of orthosilicic acid with over 2 ns of enhanced-sampling
molecular dynamics (MD) simulations and resolve the mechanism, which had not
been resolved until now in an aqueous solution.

##  2 Results and discussions

###  2.1 Reference Data Set

A preliminary training data set was constructed containing molecular silica
clusters with explicit aqueous solvation. Molecular graphs for linear and
cyclic SixOyHz were generated with x < 10 programmatically, and embedded into
3D conformers with RDkit. Solvating water molecules were added at random with
geometrical constraints to avoid clashes and to maximize hydrogen bonding. A
subset of structures were deprotonated in one or more O-H bonds and calculated
as anion, or balanced with Na+ cations. We also added secondary and common
building units of zeolites from the international zeolite database[30]. All
structures were relaxed with xTB and then refined with $\omega$-B97xD3/TZVP
level of theory in Orca, which served as our ground truth forces and energies.
Reactive geometries for proton transfer and covalent Si-O-Si reactivity were
obtained using nudged elastic band at the xTB level, followed by single-point
DFT calculations. Figure 1a shows the distribution of the data set comprising
210K geometries based on sources as well as based on the stoichiometry of the
molecule clusters.

Hybrid DFT is more accuracte for reactivity than the GGA functional, typically
used to train interatomic potentials on solids, but is restricted to a few
hundred atoms in molecule clusters. However, considering the local, atom-
centered architecture of most NNIPs, it is possible for them to learn from
cluster data and run predictions on periodic systems. Most of the interactions
in our concerning phase space are short-range coupled with electrostatic long-
range interactions which are relatively shielded in a high dielectric medium
like water[20, 31].

###  2.2 Active learning strategies

Our active learning methods rely on prediction uncertainty. We trained an
ensemble of three NNIP models with varying parameters initialized with
different random seeds. We used the variance of forces predictions as the
uncertainty metric. We implemented an uncertainty-based adversarial attack by
performing gradient-based optimization of the differentiable uncertainty
metric[29]. New molecular conformations from the data set were sampled by
back-propagating atomic displacements to find local optima that maximize the
uncertainty of the NNIP committee while balancing thermodynamic likelihood.
These new configurations are then evaluated using quantum mechanical
calculations and used to retrain the NNIPs in an active learning loop.

####  2.2.1 Differentiable uncertainty attribution

After early loops of adversarial attacks, the NNIP is stable enough to run MD
simulations of condensed phases for a decent amount of simulation time. An
adversarial attack on condensed phases will only generate condensed phase
configurations, while our data set comprises molecular clusters. Hence, we
adapted the idea of attribution. Attribution is a technique employed to test
how a change in a certain input neuron can impact an output neuron, thus
making it somewhat interpretable[32, 33]. Similarly, we calculate per-atom
attributes on uncertainty which qualitatively interpret which atoms will
contribute most to the uncertainty of the NNIP committee. It is calculated as
the derivative of the variance of forces or energy with respect to atomic
positions. Thus, from atoms with higher attributes, we can detect inter-atomic
interactions that the NNIP is highly uncertain about. A schematic of
attribution with an example is shown in Figure 1b.

We developed an active learning strategy around attribution that help us
choose uncertain atomic clusters from MD frames of condensed phases without
manual inspection. Atoms with attributes higher than two times the standard
deviation from the mean are chosen as atomic centers and all neighboring
molecules within a sphere of radius equal to the cutoff, the distance the NNIP
model uses to generate neighbors (6Å), are chosen as the molecule cluster to
add to training data. Thus a molecule in the structure remains intact in the
cluster. Only the first nearest neighbors are chosen as they have the major
contributions to the atomic energies or forces. Each MD frame of condensed
phases is used to obtain molecule clusters that can be put back in training
data after Quantum mechanical calculations. This active learning loop can work
automatically and thus provides an end-to-end framework for improved training
of NNIPs employed for condensed phase systems but trained on molecular cluster
data. Examples with active learning method are shown in Figure 1c. The NNIP
which is used to sample the first generation of molecular clusters using
attribution is tested on the same set of geometries after their DFT energies
and forces were calculated. It performs very poorly as shown in Figure 1d and
Figure 1e, thus proving the efficiency of attribution to choose atomic
environments with high errors. We further tested consecutive improvement of
attributes on a held-out set of condensed phases which was not put into
training data. Figure 1f shows a consecutive decrease in attributes over three
generations showing the generation-based improvement skills of our active
learning technique. Thus, the attribution-based technique helps sample
configurations with high error as well as improves the NNIP model on
consecutive generations, which proves the technique’s efficiency as an active
learning tool.

###  2.3 Properties of water

It has been seen based on a series of experiments that benchmarking force
error is not sufficient[34] and simulation-based statistics should be used.
Hence, it is imperative for us to begin by testing the properties of liquid
water.

####  2.3.1 Radial Distribution Functions

We compared the oxygen-oxygen (gOO) and oxygen-hydrogen (gOH) radial
distribution functions (RDFs) to experimental RDFs obtained from X-ray
diffraction measurements[35] and those obtained by reported Gaussian moment
neural network (GM-NN) potentials[20]. Thus, we compared our results with GM-
NNs trained on periodic systems of water as well as those trained on water
clusters similar to our NNIP. Figure 2a shows good agreement of gOO with the
experimental results with the first peak slightly overestimated similar to the
ones obtained by BLYP and B3LYP potentials trained on clusters. Both the peak
positions in gOH match experimental results as seen in Figure 2b but are
overestimated for all the NN potentials, which can be attributed to the lack
of Quantum nuclear effects (QNE)[36]. It is mentioned in the literature that
the O-O bond structure of water at ambient conditions is not as affected by
QNEs as the bonds involving hydrogen[18, 36] and we notice a similar behavior
for our NNIP.

####  2.3.2 Diffusion coefficient

The mean square displacements (MSDs) at different temperatures vs time are
shown in Figure 2c. At 300 K, we obtained our diffusion coefficient (D) as
0.21 $\pm$ 0.005 Å2/ps, close to an experimental value of 0.24 $\pm$ 0.015
Å2/ps[37]. Similar to RDFs, we compare our results to the GM-NN potentials
whose predicted values are 0.182 $\pm$ 0.006 Å2/ps, 0.215 $\pm$ 0.007 Å2/ps,
0.143 $\pm$ 0.004 Å2/ps, 0.139 $\pm$ 0.006 Å2/ps respectively. We plotted our
diffusion coefficients at different temperatures in Figue 2d and compare our
results with those obtained by another reported NN potential with MP2 accuracy
(DP-MP2)[19]. Our results are higher than those predicted by the DP-MP2 model
for most of the temperatures and are closer to experimental results. These
underestimated D values for the DP-MP2 model are attributed to the
overestimation of inter-molecular energy by MP2 with the aug-cc-pVDZ basis
set. Thus our potential can perform better at predicting the diffusion
properties of water at different temperatures and provide results comparable
to experimental results.

####  2.3.3 Vibrational Density of States

The vibrational density of states of water (VDOS) is shown in Figure 2e, where
we could observe that the peaks corresponding to O-H bond stretch, bending,
and libration of water at 300 K obtained using our NNIP are quite close to the
experimental observations[18]. We further employed Ring Polymer molecular
dynamics (RPMD) coupled with path integral Langevin equation thermostat
(PILE)[38] to introduce NQE. The RPMD simulations shifted the O-H stretch and
bending peaks towards the left. The O-H stretch peak is shifted from 3670 to
3600 $\mathrm{cm^{-1}}$ ( 3380 $\mathrm{cm^{-1}}$ from the experiment) and the
bending motion peak is shifted from 1664 to 1628 $\mathrm{cm^{-1}}$ ( 1637
$\mathrm{cm^{-1}}$ from the experiment). The libration has very minor changes
due to NQE, also observed for VDOS obtained by DP-MP2 model[19] and in
previous studies[18].

####  2.3.4 Equilibrium Density

Our NNIP predicted an equilibrium density of 1.08 g/cc. The GM-NN potentials
produced liquid water with densities 0.86 g/cc, 1.02 g/cc, 1.1 g/cc, and
1.12g/cc respectively[20]. Our NNIP performed equally well as the GM-NN
potentials trained on clusters. Considering our results we can claim that our
NNIP though trained on molecular clusters can predict properties of bulk
periodic water systems.

###  2.4 Crystalline Silica

We next tested our ability to predict the relative formation energies of
crystalline silica; pure siliceous zeolites to take it a notch further and
establish the success of our NNIP on periodic systems. We compared the
relative formation energies of 236 zeolites with respect to $\alpha$-Quartz
predicted by our NNIP to the ones calculated using a periodic PBE-D3 level of
DFT theory. Figure 3a depicts our NNIP performing quite closely to periodic
DFT results. We also collected experimental relative transition enthalpies for
fifteen siliceous zeolites[23] and compared them to the relative formation
energies calculated using periodic PBE-D3 DFT theory and our NNIP in Figure
3b. The NNIP performed better than the periodic DFT. This better prediction
can be due to the fact that our NNIP is trained on data obtained using a
higher level of hybrid functional DFT theory than PBE-D3. Thus, training on
cluster data has a lesser impact than the theory of DFT employed to obtain the
training data and our NNIP serves better in predicting the relative formation
energies of crystalline zeolites.

####  2.4.1 Elastic constants

We also calculated elastic constants of $\alpha$-Quartz and siliceous FAU
zeolite framework to further solidify our claim on NNIP’s capability to
predict properties of periodic systems involving silicates. We then compared
the elastic constants with those from published literature[39]. The calculated
elastic constants compared to results found in the literature are shown in
Table 1. It can be seen the NNIP ensemble is quite close to obtaining the
elastic constants for both $\alpha$-Quartz and FAU zeolite framework. Thus,
our NNIP can capture the symmetry in periodic silicates adequately and can be
used for predicting their properties.

Crystalline silica |  |  Bulk Modulus (GPa) |  Young Modulus (GPa) |  Shear Modulus (GPa) |  Poisson’s ratio
---|---|---|---|---|---
$\alpha$-Quartz |  NNIP ensemble |  51.92$\pm$10 |  112.16$\pm$8 |  49.20$\pm$3 |  0.14$\pm$0.06
literature |  38.23 |  101.41 |  47.93 |  0.06
FAU |  NNIP ensemble |  44.34$\pm$5 |  42.11$\pm$4 |  15.69$\pm$5 |  0.34$\pm$0.02
literature |  61.37 |  50.24 |  18.42 |  0.36
Table 1: Elastic constants of $\alpha$-Quartz and FAU siliceous zeolite
framework calculated by our NNIP ensemble and compared to results published in
the literature.

###  2.5 Chemical reactions

####  2.5.1 p$\mathrm{K_{a}}$ of water

To obtain the p$\mathrm{K_{a}}$ of water or p$\mathrm{K_{w}}$, we chose a
single $\mathrm{H_{2}}$O molecule in a cubic box of 100 water molecules and
varied one of the O-H bond lengths as our reaction coordinate from 0.8 to 1.9
Å with an increment of 0.05 Å. The deprotonation reaction is shown in equation
1

$\mathrm{{H_{2}O(aq)+99H_{2}O(l)\longleftrightarrow
OH^{-}+H_{3}O^{+}+98H_{2}O(l)}}$ (1)

Where $\mathrm{H_{2}O}$ (aq) is the water molecule undergoing dissociation.
The potential of mean force (PMF) vs O-H bond length ($\mathrm{r_{OH}}$) is
shown in Figure 4a. As seen from the figure, the PMF falls into a minimum
between $\mathrm{r_{OH}}$ = 0.95 and 1 Å. This is the equilibrium O-H bond
length in the neutral water molecule, thus signifying the reactant state. We
find another equilibrium point at $\mathrm{r_{OH}}$ = 1.85 Å. Here, the
$\mathrm{H^{+}}$ has left the dissociated water molecule and has been
successfully transferred to $\mathrm{H_{3}O^{+}}$. The only interaction
between $\mathrm{OH^{-}}$ and the $\mathrm{H_{3}O^{+}}$ is long-range
Coulombic now and hence this state signifies the product state where the
dissociation is complete. From the difference in PMF, we obtain our
p$\mathrm{K_{w}}$ as 15.4. The reference p$\mathrm{K_{w}}$ value is 14[40] but
considering the NQE the expected p$\mathrm{K_{w}}$ becomes 17[41]. Our result
is also comparable to a reported value of 16 obtained by first principle
calculations using SCAN functional[42]. They used the same procedure as us to
find the p$K_{a}$ of water.

####  2.5.2 p$\mathrm{K_{a}}$ of silicate oligomers

We also calculated the p$\mathrm{K_{a}}$ of orthosilicic acid, silicate dimer,
and trimer. Orthosilicic acid and dimer have all identical O-H bonds. Hence
dissociating any one of them is fine. However, for trimer two O-H bonds
attached to the middle Si are different than the end ones. We chose one of the
end O-H bond dissociations for trimer to get its p$\mathrm{K_{a}}$. Similar to
water, we took a single orthosilicic acid, dimer, and trimer in a cubic box of
100 water molecules and varied their chosen O-H bond lengths. The individual
PMFs vs $\mathrm{r_{OH}}$ are shown in Figure 4b-d. The deprotonation reaction
of orthosilicic acid looks as shown in equation 2.

$\mathrm{Si(OH)_{4}(aq)+100H_{2}O(l)\longleftrightarrow
SiO_{4}H_{3}^{-}+H_{3}O^{+}+99H_{2}O(l)}$ (2)

We found the equilibrium reactant state at $\mathrm{r_{OH}}$ between 0.95 and
1 Å and product state at $\mathrm{r_{OH}}$ = 1.55 Å for orthosilicic acid and
between 1.5 and 1.55 Å for the dimer and the trimer. From the PMF, we obtain
the p$\mathrm{K_{a}}$ of orthosilicic acid as 10.45. Orthosilicic acid is a
weak acid with a known p$\mathrm{K_{a}}$ of 9.8[43]. Further, for dimer and
trimer, we obtained a p$K_{a}$ of 10.43 and 10.42 respectively, whose
p$\mathrm{K_{a}}$ values have been reported to vary between 9.5 and 10.7[44].
Our reactive NNIP is thus shown to be capable of replicating complex reactions
involving pure water as well as amorphous silicates in water.

####  2.5.3 Silicate dimerization

We next simulated a dimerization reaction of two silicate monomers in a
neutral solution with 100 water molecules at 300K. We ran umbrella sampling on
a reaction coordinate, $\mathrm{d_{1}-d_{2}}$ varying it from -5 Å to 5 Å with
an increment of 0.05 Å. $\mathrm{d_{2}}$ signifies the bond distance between
attacking oxygen and the attacked silicon, whereas $\mathrm{d_{1}}$ signifies
the distance between the leaving oxygen and the attacked silicon. Our
attacking oxygen, attacked silicon, and the leaving oxygen atoms for the
$S_{N}2$ and the lateral attack mechanisms are shown in Figure S1. The
reported $S_{N}2$ mechanism and our observed one are shown in Figure S1a-b. In
the case of the lateral attack mechanism, we saw the reaction happening as
explained in literature[25, 24, 26], where the leaving oxygen leaves with the
hydrogen from the attacking oxygen (Figure S1c). From the free energy profile
in figure 4e, we see that at the coordinate close to 5 Å, the product dimer is
formed. The reactant state with two monomers is at the coordinates between -3
and -5 Å where the free energy is almost flat. The reaction energy was
obtained by subtracting these two energies which shows that the product phase
is only 0.2 kJ/mole higher in energy than the reactant phase which can be
noise. The activation energy is obtained as 103 kJ/mole which is closer to 129
kJ/mole found in a DFT study, with solvent effect incorporated and product
phase found less stable with free energy 9 kJ/mole higher[24]. From the plot
of the coordination between the attacking oxygen with the transferring
hydrogen and that of the leaving oxygen with the transferring hydrogen in
Figure 4e, we observed that at the transition state, the hydrogen leaves the
attacking oxygen with a drop in their coordination value and moves to the
leaving oxygen with an increase in their coordination value almost immediately
as expected in the lateral flank-side mechanism. Hence, we can conclude that
silicate dimerization in an aqueous solution at neutral conditions occurs
through a lateral flank-side attack mechanism with pentavalent silicon as a
transition state.

##  3 Conclusion

Neural Network potentials have proven to work as a bridge between first
principle calculations and classical potentials several times. In this work,
we have proved the same by developing a reactive equivariant NNIP trained on a
complex domain of silicate-water interactions. Our NNIP though trained on
molecule clusters can predict properties of bulk periodic systems including
liquid water as well as crystalline solids. We further proposed a new active
learning strategy based on attribution of differentiable uncertainty which
further displays highly uncertain inter-atomic interactions in amorphous
systems. This method not only serves as a qualitative measure for unsure
atomic environments but also serves to minimize computational costs by
extracting only the uncertain atomic environments. The NNIP trained on this
strategy then proved to be successful in predicting the properties of bulk
water as well as crystalline silicates. The reactive potential is further
adept at predicting the reaction path of deprotonation of water, small
silicate oligomers, and silica dimerization reaction in a water solution. This
NNIP can further be used to probe into other silicate polymerization reactions
leading to precipitated silica. We can also study the impact of pH and
temperature on these reactions in the future as we have sodium in our data set
to maintain charge neutrality or we can add aluminum in our data set to train
a potential on preliminary stages of the synthesis of zeolites.

##  4 Methodology

###  4.1 Model Architecture

Our NN potential is based on PaiNN architecture[45], which is an equivariant
neural network with message passing as its backbone. Equivariant models[45,
46, 47] can act on non-invariant inputs like displacement vectors in a
symmetry-respecting way and thus offer good accuracy with a low amount of
training data for properties that are equivariant in Euclidean space like
forces, thus being more data sufficient than invariant potentials like
Schnet[48]. In this method, an invariant feature vector for each atom is
generated with their atomic numbers which are then updated through
convolutions with "messages" from neighboring atoms which consist of distance,
orientation, and features of nearest neighbor atoms within a cutoff distance.
Each atomic feature is updated with information from neighbors and through
convolutions obtains information from atoms far away as well, thereby
generating a representation vector as a function of atomic positions and
orientations, which is then mapped to atomic energies through feed-forward NN
and summed to obtain the energy of the system. Through automatic
differentiation, the forces on each atom can also be obtained. We further
modified the model to predict stresses of a system based on the virial
theorem[49] as given in equation 3 for pair potentials

${\sigma_{ij}^{V}}=\frac{1}{V}\sum_{\alpha}\left[\sum_{\begin{subarray}{c}\beta=1\\\
\beta>\alpha\end{subarray}}^{N}\left(r_{i}^{\beta}-r_{i}^{\alpha}\right)f_{j}^{\alpha\beta}-m^{\alpha}v_{i}^{\alpha}v_{j}^{\alpha}\right],\qquad
f_{j}^{\alpha\beta}=\frac{\partial{E}}{\partial{\left(r_{j}^{\beta}-r_{j}^{\alpha}\right)}}$
(3)

where (i,j) signifies x,y, and z directions, $\beta$ varies from 1 to N
neighbors of atom $\alpha$, $r_{i}^{\beta}$, and $r_{i}^{\alpha}$ are
positions of atom $\beta$ and $\alpha$ along direction i respectively,
$f_{j}^{\alpha\beta}$ is the force on $\alpha$ by $\beta$ along direction j
and E is total energy, V is total volume, $m^{\alpha}$ is the mass of atom
$\alpha$ and $v^{\alpha}$ is the thermal velocity of atom $\alpha$.

###  4.2 Training details

Molecular geometries of the same stoichiometry are split at a ratio of 3:1:1
into training, validation, and test sets. We employed a mini-batch gradient
descent optimization with Adam optimizer. The learning drops from 10-4 to 10-6
at a rate of 0.5 when the validation loss hits a plateau for 20 epochs. The
loss function is the same as described elsewhere[45], with a loss coefficient
of 0.95 for forces and 0.05 for energies. We fixed the number of convolutional
layers at 4 and the cutoff for nearest neighbors at 6 Å. We optimized the
other hyperparameters with Sigopt on validation loss. The optimized dimension
of the feature vector for each atom is 310; the number of radial basis
functions used are 20 and the batch size is chosen as 9.

###  4.3 Active learning

We ran adversarial attacks on molecular geometries randomly chosen from the
preliminary data set for three generations. Then we used our model to run MD
simulations on 10-100 orthosilicates with different ratios of water molecules
at different temperatures of 300, 500, 750, and 1000 K and ran attribution-
based active learning on them. We also ran simulations replicating the
deprotonation of water, orthosilicic acid, and other silicates. We further ran
attributions on zeolite crystals and dimerization reactions to extract
uncertain molecule clusters from them. Our active learning consists of nine
cycles, involving three from adversarial attacks and then six generations of
attribution-based active learning.

###  4.4 Simulation details

We used the atomic simulation environment (ASE) package coupled with our NNIP
to run simulations. All periodic structures with silicates and water molecules
are generated using Packmol and rendered using OVITO. The simulation boxes
were generated with the equilibrium density of water predicted by our NNIP of
1.08 g/cc. All MD simulations were run with canonical (NVT) ensemble with
Nosé-Hoover thermostat and timestep of 0.5 fs.

####  4.4.1 Radial distribution functions of water

We obtain the oxygen-oxygen (gOO) and oxygen-hydrogen (gOH) RDFs from 100 ps
long MD simulations of a periodic box containing 400 H2O molecules in
canonical (NVT) ensemble at 300 K after 20 ps equilibration period. We divide
the whole trajectory into ten parts and take the average RDF curve as the
estimated one.

####  4.4.2 Diffusion coefficients of water

The diffusion coefficients were calculated from the MSDs by Einstein’s
diffusion equation as given in equation 4

${D=\frac{MSD}{6\Delta t}=\frac{\langle{|x(t+\Delta
t)-x(t)|}^{2}\rangle}{6\Delta t}}$ (4)

Where D is the diffusion coefficient, $\Delta t$ is the time interval used for
measuring MSDs. We used 10 ps as the time interval. We further obtained the
slope of ln(MSD) vs ln(t) and choose the region for linear fitting where the
slope is within 1 $\pm$ $10^{-4}$. The MSDs were calculated based on the
relative squared displacement of oxygen atoms and averaged over trajectories
for all water molecules. We further needed to extrapolate the diffusion
coefficient to infinite-system size owing to the finite-size effect[50]. The
extrapolated diffusion coefficient (D($\infty$)) is obtained from finite-
system size diffusion coefficient (D(L)) with an extra correction term as
shown in equation 5

${D(\infty)=D(L)+\frac{\zeta}{6\pi\beta\eta L}}$ (5)

Where L is the simulation box size, $\zeta$ depends on the geometry of the
simulation box (2.837297 for a cubic box), $\beta$ is $\frac{1}{{k_{B}}T}$,
and $\eta$ is experimental shear viscosity = 0.8925 x $10^{-3}$ Pa s. We
calculated the diffusion coefficients by first running an NVT simulation for
50 ps to equilibrate the system and then ran five independent NVE simulations
for 150 ps at different temperatures of 260, 280, 300, 320, 340, and 360 K and
chose the average MSDs.

####  4.4.3 Vibrational density of states of water

The VDOS of water can be procured from the Fourier transform of the velocity-
velocity autocorrelation function as given in equation 6

${C_{vv}(\omega)=\int\langle{v(\tau)v(\tau+t)}_{\tau}\rangle e^{-i\omega
t}dt}$ (6)

Where v($\tau$) is the centroid velocity of the atom and $\omega$ is the
vibrational frequency. The RPMD simulations were implemented using
SchnetPack[51] with 5 beads and 0.25fs timestep.

####  4.4.4 Equilibrium density of water

We ran MD simulations in NPT (constant pressure) ensemble using the isobaric-
isothermal form of the Nosé–Hoover dynamics at a temperature of 300K and a
pressure of 1 atm. We chose a cubic box of 400 water molecules and the box was
allowed to relax in a hydrostatic manner. A 1 ns long trajectory after 100 ps
of equilibration period was divided into ten parts and the time-averaged
volume of the simulation box was taken as equilibrium volume. The equilibrium
density was then calculated from it.

####  4.4.5 Elastic constants of silicates

We provided normal strains of -1% to 1% in steps of 0.5% and shear strains of
-4% to 4% in steps of 2% to the equilibrium structures of $\alpha$-Quartz and
FAU zeolite and obtained their deformed structures using Pymatgen. We then
optimized each deformed structure by the BFGS algorithm implemented in the ASE
package. We then obtained the stresses at different strains and their
corresponding slopes to measure the elastic tensors for each system. We then
used ELATE software[39] to obtain elastic constants from the tensors. We used
our NNIP ensemble to procure the average elastic constants and their standard
deviation.

####  4.4.6 Chemical reactions

We performed constrained MD simulations with umbrella sampling varying our
reaction coordinates chosen for different reactions. We then obtained a smooth
reaction path from reactant to product state by calculating the PMF as a
function of our chosen reaction coordinates. The PMF was calculated using the
Multistate Bannett’s Acceptance Ratio (MBAR)[52] on the sampled trajectories.
The MBAR analysis was carried out using codes implemented in adaptive sampling
package[53, 54]. We can use the PMF of the deprotonation reactions of water
and silicates to obtain their p$K_{a}$ values. The difference in the PMF of
reactant and product phase ($\Delta$F) is related to the acid dissociation
constant ($K_{c}$) as shown in equation 7 and the relation between p$K_{a}$
and $K_{c}$ is shown in equation 8, Combining both we can get p$K_{a}$ from
$\Delta$F as shown in equation 9.

${\Delta F=-k_{B}T\ln k_{c}}$ (7) ${pK_{a}=-\log k_{c}}$ (8)
${pK_{a}=\frac{\Delta F}{k_{B}T\times 2.3026}}$ (9)

Where $k_{B}$ is the Boltzmann constant.

##  5 Acknowledgement

R.G-B. and S.R. acknowledge funding support from Evonik AG and the MIT-IBM
Watson AI Lab. We acknowledge the MIT Engaging cluster at the Massachusetts
Green High-Performance Computing Center (MGHPCC) for providing high-
performance computing resources.

##  6 Data Availability

The silica-water data set can be obtained from the authors on request and the
NNIP model generated during this study is available at Materials Cloud Archive
under accession code https://doi.org/10.24435/materialscloud:61-x8. The codes
used for this study can be downloaded from https://github.com/learningmatter-
mit/NeuralForceField.

##  7 References

## References

  * [1] Peter J Heaney, Charles T Prewitt and Gerald V Gibbs “Silica: Physical behavior, geochemistry, and materials applications” Walter de Gruyter GmbH & Co KG, 2018
  * [2] Jens Weitkamp “Zeolites and catalysis” In _Solid state ionics_ 131.1-2 Elsevier, 2000, pp. 175–188
  * [3] Yating Gao et al. “Multifunctional Role of Silica in Pharmaceutical Formulations” In _AAPS PharmSciTech 2022 23:4_ 23.4 Springer, 2022, pp. 1–18 DOI: 10.1208/S12249-022-02237-5
  * [4] Naomi J. Halas “Nanoscience under glass: The versatile chemistry of silica nanostructures” In _ACS Nano_ 2.2 American Chemical Society, 2008, pp. 179–183 DOI: 10.1021/NN800052E/ASSET/IMAGES/LARGE/NN-2008-00052E{\\_}0005.JPEG
  * [5] HE Bergna and WO Roberts “Colloidal silica: fundamentals and applications”, 2005 URL: https://books.google.com/books?hl=en&lr=&id=d0huBwAAQBAJ&oi=fnd&pg=PP1&ots=uYjxTawhd_&sig=Pf92yxgRxSVaJ1on0L2VdXsvb_c
  * [6] Nadeem Raza et al. “Synthesis and characterization of amorphous precipitated silica from alkaline dissolution of olivine” In _RSC Advances_ 8.57 The Royal Society of Chemistry, 2018, pp. 32651–32658 DOI: 10.1039/C8RA06257A
  * [7] R. Dewati et al. “Precipitated Silica from Pumice and Carbon Dioxide Gas (Co2) in Bubble Column Reactor” In _Journal of Physics: Conference Series_ 953.1 IOP Publishing, 2018, pp. 012226 DOI: 10.1088/1742-6596/953/1/012226
  * [8] Colin S. Cundy and Paul A. Cox “The hydrothermal synthesis of zeolites: History and development from the earliest days to the present time” In _Chemical Reviews_ 103.3 American Chemical Society, 2003, pp. 663–701 DOI: 10.1021/cr020060i
  * [9] Durgesh K. Rai et al. “In situ study of aggregate topology during growth of pyrolytic silica” In _Journal of Aerosol Science_ 118 Elsevier Ltd, 2018, pp. 34–44 DOI: 10.1016/J.JAEROSCI.2018.01.006
  * [10] Manuel Meier et al. “Multi-scale characterization of precipitated silica” In _Powder Technology_ 354 Elsevier B.V., 2019, pp. 45–51 DOI: 10.1016/J.POWTEC.2019.05.072
  * [11] B..H. Van Beest, G.. Kramer and R.. Van Santen “Force fields for silicas and aluminophosphates based on ab initio calculations” In _Physical Review Letters_ 64.16 American Physical Society, 1990, pp. 1955 DOI: 10.1103/PhysRevLett.64.1955
  * [12] A. Carré, J. Horbach, S. Ispas and W. Kob “New fitting scheme to obtain effective potential from Car-Parrinello molecular-dynamics simulations: Application to silica” In _Europhysics Letters_ 82.1 IOP Publishing, 2008, pp. 17001 DOI: 10.1209/0295-5075/82/17001
  * [13] E. Flikkema and S.. Bromley “A new interatomic potential for nanoscale silica” In _Chemical Physics Letters_ 378.5-6 North-Holland, 2003, pp. 622–629 DOI: 10.1016/J.CPLETT.2003.07.017
  * [14] Joseph C. Fogarty et al. “A reactive molecular dynamics simulation of the silica-water interface” In _The Journal of Chemical Physics_ 132.17 American Institute of PhysicsAIP, 2010, pp. 174704 DOI: 10.1063/1.3407433
  * [15] J.. Rimsza, Jejoon Yeon, A..T. Van Duin and Jincheng Du “Water Interactions with Nanoporous Silica: Comparison of ReaxFF and ab Initio based Molecular Dynamics Simulations” In _Journal of Physical Chemistry C_ 120.43 American Chemical Society, 2016, pp. 24803–24816 DOI: 10.1021/ACS.JPCC.6B07939/SUPPL{\\_}FILE/JP6B07939{\\_}SI{\\_}004.TXT
  * [16] Felix Musil et al. “Physics-inspired structural representations for molecules and materials” In _Chemical Reviews_ 121.16 American Chemical Society, 2021, pp. 1–60 DOI: 10.1021/acs.chemrev.1c00021
  * [17] Simon Axelrod et al. “Learning Matter: Materials Design with Machine Learning and Atomistic Simulations” In _Accounts of Materials Research_ 3.3, 2022, pp. 343–357 DOI: 10.1021/accountsmr.1c00238
  * [18] Yi Yao and Yosuke Kanai “Nuclear Quantum Effect and Its Temperature Dependence in Liquid Water from Random Phase Approximation via Artificial Neural Network” In _Journal of Physical Chemistry Letters_ 12.27 American Chemical Society, 2021, pp. 6354–6362 DOI: 10.1021/ACS.JPCLETT.1C01566/SUPPL{\\_}FILE/JZ1C01566{\\_}SI{\\_}003.ZIP
  * [19] Jinfeng Liu, Jinggang Lan and Xiao He “Toward High-level Machine Learning Potential for Water Based on Quantum Fragmentation and Neural Networks” In _Journal of Physical Chemistry A_ 126.24 American Chemical Society, 2022, pp. 3926–3936 DOI: 10.1021/ACS.JPCA.2C00601/SUPPL{\\_}FILE/JP2C00601{\\_}SI{\\_}002.ZIP
  * [20] Viktor Zaverkin, David Holzmüller, Robin Schuldt and Johannes Kästner “Predicting properties of periodic systems from cluster data: A case study of liquid water” In _Journal of Chemical Physics_ 156.11 AIP Publishing LLCAIP Publishing, 2022, pp. 114103 DOI: 10.1063/5.0078983
  * [21] Linus C. Erhard, Jochen Rohrer, Karsten Albe and Volker L. Deringer “A machine-learned interatomic potential for silica and its relation to empirical models” In _npj Computational Materials 2022 8:1_ 8.1 Nature Publishing Group, 2022, pp. 1–12 DOI: 10.1038/s41524-022-00768-w
  * [22] I.. Balyakin, S.. Rempel, R.. Ryltsev and A.. Rempel “Deep machine learning interatomic potential for liquid silica” In _Physical Review E_ 102.5 American Physical Society, 2020, pp. 52125 DOI: 10.1103/PhysRevE.102.052125
  * [23] Andreas Erlebach, Petr Nachtigall and Lukáš Grajciar “Accurate large-scale simulations of siliceous zeolites by neural network potentials” In _arXiv preprint arXiv:2102.12404_ , 2021, pp. 1–33 URL: http://arxiv.org/abs/2102.12404
  * [24] Thuat T. Trinh, Antonius P.J. Jansen and Rutger A. Van Santen “Mechanism of oligomerization reactions of silica” In _Journal of Physical Chemistry B_ 110.46 American Chemical Society, 2006, pp. 23099–23106 DOI: 10.1021/JP063670L/ASSET/IMAGES/MEDIUM/JP063670LU00008A.GIF
  * [25] Xue Qing Zhang, Thuat T. Trinh, Rutger A. Van Santen and Antonius P.J. Jansen “Mechanism of the initial stage of silicate oligomerization” In _Journal of the American Chemical Society_ 133.17 American Chemical Society, 2011, pp. 6613–6625 DOI: 10.1021/JA110357K/SUPPL{\\_}FILE/JA110357K{\\_}SI{\\_}001.PDF
  * [26] Charles L. Schaffer and Kendall T. Thomson “Density functional theory investigation into structure and reactivity of prenucleation silica species” In _Journal of Physical Chemistry C_ 112.33 American Chemical Society, 2008, pp. 12653–12662 DOI: 10.1021/JP066534P/SUPPL{\\_}FILE/JP066534P-FILE002.PDF
  * [27] J..G. Pereira, C..A. Catlow and G.. Price “Silica condensation reaction: an ab initio study” In _Chemical Communications_ The Royal Society of Chemistry, 1998, pp. 1387–1388 DOI: 10.1039/A801816B
  * [28] Mohamed Elanany et al. “A quantum molecular dynamics simulation study of the initial hydrolysis step in sol-gel process” In _Journal of Physical Chemistry B_ 107.7 American Chemical Society, 2003, pp. 1518–1524 DOI: 10.1021/JP026816Z/ASSET/IMAGES/LARGE/JP026816ZF00008.JPEG
  * [29] Daniel Schwalbe-Koda, Aik Rui Tan and Rafael Gómez-Bombarelli “Differentiable sampling of molecular geometries with uncertainty-based adversarial attacks” In _Nature Communications_ 12.1 Springer US, 2021, pp. 1–12 DOI: 10.1038/s41467-021-25342-8
  * [30] Ch. Baerlocher and L.B. McCusker “Database of Zeolite Structures” In _http://www.iza-structure.org/databases/_ , 2021 URL: http://www.iza-structure.org/databases/
  * [31] Tobias Morawietz, Andreas Singraber, Christoph Dellago and Jörg Behler “How van der waals interactions determine the unique properties of water” In _Proceedings of the National Academy of Sciences of the United States of America_ 113.30 National Academy of Sciences, 2016, pp. 8368–8373 DOI: 10.1073/pnas.1602375113
  * [32] Avanti Shrikumar, Peyton Greenside and Anshul Kundaje “Learning important features through propagating activation differences” In _34th International Conference on Machine Learning, ICML_ 7, 2017, pp. 4844–4866
  * [33] Mukund Sundararajan, Ankur Taly and Qiqi Yan “Axiomatic Attribution for Deep Networks” In _34th International Conference on Machine Learning, ICML_ 7 International Machine Learning Society (IMLS), 2017, pp. 5109–5118 URL: https://arxiv.org/abs/1703.01365v2
  * [34] Xiang Fu et al. “Forces are not Enough: Benchmark and Critical Evaluation for Machine Learning Force Fields with Molecular Simulations”, 2022 DOI: 10.48550/arxiv.2210.07237
  * [35] Lawrie B. Skinner et al. “Benchmark oxygen-oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range” American Institute of PhysicsAIP, 2013, pp. 074506 DOI: 10.1063/1.4790861
  * [36] Ondrej Marsalek and Thomas E. Markland “Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects” In _Journal of Physical Chemistry Letters_ 8.7 American Chemical Society, 2017, pp. 1545–1551 URL: https://pubs.acs.org/doi/full/10.1021/acs.jpclett.7b00391%20https://pubs.acs.org/doi/abs/10.1021/acs.jpclett.7b00391
  * [37] M. Holz, S.. Heil and A. Sacco “Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1H NMR PFG measurements” In _Physical Chemistry Chemical Physics_ 2.20 The Royal Society of Chemistry, 2000, pp. 4740–4742 DOI: 10.1039/B005319H
  * [38] Michele Ceriotti, Michele Parrinello, Thomas E. Markland and David E. Manolopoulos “Efficient stochastic thermostatting of path integral molecular dynamics” In _The Journal of Chemical Physics_ 133.12 American Institute of PhysicsAIP, 2010, pp. 124104 DOI: 10.1063/1.3489925
  * [39] Romain Gaillac, Pluton Pullumbi and François Xavier Coudert “ELATE: an open-source online application for analysis and visualization of elastic tensors” In _Journal of Physics: Condensed Matter_ 28.27 IOP Publishing, 2016, pp. 275201 DOI: 10.1088/0953-8984/28/27/275201
  * [40] Todd P. Silverstein and Stephen T. Heller “PKa Values in the Undergraduate Curriculum: What Is the Real pKa of Water?” In _Journal of Chemical Education_ 94.6 American Chemical Society, 2017, pp. 690–695 DOI: 10.1021/acs.jchemed.6b00623
  * [41] Michele Ceriotti et al. “Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges” In _Chemical Reviews_ 116.13 American Chemical Society, 2016, pp. 7529–7550 DOI: 10.1021/ACS.CHEMREV.5B00674/ASSET/IMAGES/CR-2015-00674D{\\_}M006.GIF
  * [42] Ruiyu Wang, Vincenzo Carnevale, Michael L. Klein and Eric Borguet “First-Principles Calculation of Water p Ka Using the Newly Developed SCAN Functional” In _Journal of Physical Chemistry Letters_ 11.1 American Chemical Society, 2020, pp. 54–59 DOI: 10.1021/ACS.JPCLETT.9B02913/ASSET/IMAGES/LARGE/JZ9B02913{\\_}0006.JPEG
  * [43] Carole C. Perry “Biogenic Silica: A Model of Amorphous Structure Control” In _Growth, Dissolution and Pattern Formation in Geosystems_ Springer, Dordrecht, 1999, pp. 237–251 DOI: 10.1007/978-94-015-9179-9{\\_}11
  * [44] David J. Belton, Olivier Deschaume and Carole C. Perry “An overview of the fundamentals of the chemistry of silica with relevance to biosilicification and technological advances” In _The FEBS journal_ 279.10 FEBS J, 2012, pp. 1710–1720 DOI: 10.1111/J.1742-4658.2012.08531.X
  * [45] Kristof T. Schütt, Oliver T. Unke and Michael Gastegger “Equivariant message passing for the prediction of tensorial properties and molecular spectra” In _arXiv preprint arXiv:2102.03150_ , 2021 URL: http://arxiv.org/abs/2102.03150
  * [46] Simon Batzner et al. “E(3)-Equivariant Graph Neural Networks for Data-Efficient and Accurate Interatomic Potentials” In _Nature Communications_ 13.1 Nature Publishing Group, 2021, pp. 1–11 DOI: 10.1038/s41467-022-29939-5
  * [47] Albert Musaelian et al. “Learning Local Equivariant Representations for Large-Scale Atomistic Dynamics”, 2022 URL: http://arxiv.org/abs/2204.05249
  * [48] K.. Schütt et al. “SchNet: A continuous-filter convolutional neural network for modeling quantum interactions” In _Advances in Neural Information Processing Systems_ 2017-Decem.1, 2017, pp. 992–1002
  * [49] Arun K. Subramaniyan and C.. Sun “Continuum interpretation of virial stress in molecular simulations” In _International Journal of Solids and Structures_ 45.14-15 Pergamon, 2008, pp. 4340–4346 DOI: 10.1016/j.ijsolstr.2008.03.016
  * [50] In Chul Yeh and Gerhard Hummer “System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions” In _Journal of Physical Chemistry B_ 108.40 American Chemical Society, 2004, pp. 15873–15879 DOI: 10.1021/JP0477147/ASSET/IMAGES/LARGE/JP0477147F00008.JPEG
  * [51] K.. Schütt et al. “SchNetPack: A Deep Learning Toolbox for Atomistic Systems” In _Journal of Chemical Theory and Computation_ 15.1 American Chemical Society, 2019, pp. 448–455 DOI: 10.1021/acs.jctc.8b00908
  * [52] Michael R Shirts and John D Chodera “Statistically optimal analysis of samples from multiple equilibrium states”, 2008 DOI: 10.1063/1.2978177
  * [53] Johannes C.. Dietschreit et al. “From Free-Energy Profiles to Activation Free Energies”, 2022 URL: http://arxiv.org/abs/2206.02893
  * [54] Andreas Hulm, Johannes C.B. Dietschreit and Christian Ochsenfeld “Statistically optimal analysis of the extended-system adaptive biasing force (eABF) method” In _Journal of Chemical Physics_ 157.2 American Institute of Physics Inc., 2022, pp. 24110 DOI: 10.1063/5.0095554/2841454

Figure 1:  Training dataset and active learning method. a) Distribution of
data based on data source - Cluster: silicate and water clusters based on
random position guesses. Opt: molecule clusters optimized at a lower level of
DFT theory or X-TB. AL: data produced through active learning and NEB:
transition states using NEB; and based on stoichiometry - Si$\mathrm{O_{n}}$:
silicate clusters in a vacuum. Si$\mathrm{O_{n}}$-$\mathrm{H_{2}}$O: silicate
clusters surrounded by water molecules, $\mathrm{H_{2}}$O: water molecules
only and with Na: molecular clusters of silicate or water with Na in it. b)
Schematic diagram of attribution. The potential energy and forces and their
uncertainty are calculated for a molecular system using a NN ensemble. The
uncertainty is then differentiated with respect to atom positions to obtain
individual per-atom uncertainty attributes. The neighborhood of the hydrogen
atom in orthosilicic acid which is attracted by the adjacent water molecule
has high attributes. This shows that the schematic potential is bad at
hydrogen transfer from orthosilicic acid to adjacent water. c) Two examples of
attribution. The color bar is based on attributes. The red circles show the
atomic environments added for active learning with atoms with high attributes
at the center. The $\mathrm{H_{3}O^{+}-H_{2}O}$,
$\mathrm{H_{3}O^{+}-H_{3}O^{+}}$, $\mathrm{H_{3}O^{+}-OH^{-}}$,
$\mathrm{Si(OH)_{4}-H_{2}O}$ interactions and $\mathrm{Si(OH)_{3}OH_{2}^{+}}$
complex lead to high attributes as they were not present in the training data.
Parity plot of d) forces and e) energies of a set of molecular clusters
obtained through attribution-based active learning. The first-generation NNIP
performs poorly on them. f) Consecutive generation-wise improvement of
attributes on a held-out condensed phase system. Attribution-based active
learning thus improves our model.  Figure 2: Properties of water. a) O-O and
b) O-H radial distribution functions at 300 K calculated by NVT simulations
with our NNIP and compared to experimental and GM-NN potentials’ results.NNIP:
our Neural network inter-atomic potential; EXP: experimental data; revPBE-
bulk: GM-NN potential trained on periodic water box calculated at the
revPBE-D3 level of theory; revPBE-cluster: GM-NN potential trained on water
clusters calculated at the revPBE-D3 level of theory; B3LYP: GM-NN potential
trained on water clusters calculated at the B3YLP-D3 level of theory; BLYP:
GM-NN potential trained on water clusters calculated at the BYLP-D3 level of
theory. c) Mean squared displacements vs time at different temperatures from
260K to 360K. d) diffusivity vs temperature obtained using our NNIP compared
with results from a DP-MP2 potential. e) Vibrational density of states of
liquid water obtained from classical (MD) and Ring polymer (RPMD) molecular
simulations with our NNIP. The dashed vertical lines represent the
experimental infrared spectrum. Figure 3: Properties of crystalline silica. a)
Relative formation energies of pure siliceous zeolites with respect to
$\alpha$-Quartz calculated using our NNIP compared to those obtained by PBE-D3
level of DFT theory. b) Relative formation energies of fifteen zeolites
obtained by NNIP and PBE-D3 compared to their experimental transition
enthalpies. Figure 4: Simulations of chemical reactions. a) PMF of water
deprotonation. The insets show reactant state (0.95 Å $<$ $\mathrm{r_{OH}}$
$<$ 1 Å) and product state ($\mathrm{r_{OH}}$ = 1.85 Å). The free energy of
the reactant state is set to 0. b) PMF of orthosilicic acid deprotonation. The
insets show reactant state (0.95 Å $<$ $\mathrm{r_{OH}}$ $<$ 1 Å) and product
state ($\mathrm{r_{OH}}$ = 1.55 Å). c) PMF of dimer deprotonation. d) PMF of
trimer deprotonation. e) Free energy profile of silica dimerization,
coordination between attacking O and transferring hydrogen and that between
leaving oxygen and the transferring hydrogen along the reaction coordinate
$\mathrm{d_{1}-d_{2}}$ (Å). The insets show the reactant, product, and
transition states.

##  8 Supplementary Materials

Figure S1: Schemes of Silicate dimerization reactions. a) Expected scheme of
$S_{N}2$ back side attack mechanism. The hydrogen from the attacking oxygen
goes to the lateral oxygen, which then loses hydrogen to the leaving oxygen.
b) Observed scheme of back-side attack mechanism. When we chose the back side
oxygen as leaving oxygen to replicate $S_{N}2$ mechanism, a surrounding water
molecule interfered and took up hydrogen from the attacking oxygen forming a
$\mathrm{H_{3}O^{+}}$ and $\mathrm{OH^{-}}$. c) Observed lateral flank-side
attack mechanism. This mechanism was not hampered by surrounding water
molecules and occurred as reported in the literature. The transition state has
a pentavalent Silicon atom.