Acid-base dissociation mechanisms and energetics at the silica–water interface: An activationless process Journal of Colloid and Interface Science

Hypothesis: Silanol groups at the silica–water interface determine not only the surface charge, but also have an important role in the binding of ions and biomolecules. As the pH is increased above pH 2, the silica surface develops a net negative charge primarily due to deprotonation of the silanol group. An improved understanding of the energetics and mechanisms of this fundamentally important process would further understanding of the relevant dynamics. Simulations: Density Functional Theory ab initio molecular dynamics and geometry optimisations were used to investigate the mechanisms of surface neutralisation and charging in the presence of OH (cid:2) and H 3 O þ respectively. This charging mechanism has received little attention in the literature. Findings: The protonation or deprotonation of isolated silanols in the presence of H 3 O þ or OH (cid:2) , respectively, was shown to be a highly rapid, exothermic reaction with no signiﬁcant activation energy. This process occurred via a concerted motion of the protons through ‘water wires’. Geometry optimisations of large water clusters at the silica surface demonstrated proton transfer to the surface occurring via the rarely discussed ‘proton holes’ mechanism. This indicates that surface protonation is possible even when the hydronium ion is distant (at least 4 water molecules separation) from the surface. (cid:2) 2015 The Authors. Published by Elsevier Inc. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Silica and water represent two of the most abundant chemical systems, and therefore it is unsurprising that understanding the interface between them is relevant to a wide variety of systems. Surface charging is fundamental to a range of phenomena such as dissolution rates [1][2][3] and the surface adsorption of ions and molecules [4,5]. Chemical reactions of reactive silanol groups ðSi-OHÞ with H þ =OH À are thought to be the primary surface charging mechanism for silica [6], with electrolyte effects having a measurable but less significant effect [7].

Oxide surface charging
One of the most popular approaches to modelling surface proton reactions has been Surface Complexation Models (SCM). One such example is the 2-pK model which assumes that the surface state can be modelled as two consecutive protonation reactions, the equilibrium constants of which are often obtained empirically from acid-base titration data [8,9]. This methodology has the advantage of being highly generalisable, but neglects any direct information obtained about the oxide-water interfacial structure and dynamics obtained experimentally or via simulation. While the thermodynamic description of the 2-pK model is suitable for describing certain processes such as ion complexation and surface dissolution [8,10], it provides little insight into the atomic-scale interactions which are present at water-surface interfaces and are required to understand many dynamic interfacial processes, such as double layer formation, solvent structure, surface-charging kinetics and non-equilibrium interfacial processes.
As the pH is increased above 2 ± 1, the silica surface becomes increasingly negatively charged [11]. Silicon surface atoms (S surf ) are present as a neutral oxide that can dissociate or protonate according to the following chemical equilibria: X-ray photoelectron spectroscopy measurements have shown that significant quantities of Si-OH þ 2 are only present at extremely low pH ($2 or lower) [12]. High level ab initio calculations have demonstrated the chemical instability of this species in neutral water, supporting the notion that at conditions relevant to most applications (pH [2][3][4][5][6][7][8][9][10][11][12][13][14] this species is an insignificant component of the surface composition [13]. Given this observation, Borkovec explained how the 2-pK model should be interpreted as reactions between a pair of neighbouring ionizable groups, as opposed to a single site as has often been assumed [14]. Although the 2-pK model can be used to describe the average surface charge from a system-scale perspective by considering only the reaction with hydronium ions, the surface at an atomistic scale presents a far more complex environment in which Eq. (2) is competing with the analogous hydroxyl reaction show in Eq. (3).

Charge transport in pure water
Given the chemical similarity between proton transfer of water-water proton transfers and of silanol-water proton transfers, the mechanism of proton transport in pure water is relevant. In pure water, solvated hydronium ions are transported via the Grotthuss mechanism [15] which involves interconversion between the symmetric Zundel cation H 5 O þ 2 and the triply-hydrogen bonded Eigen cation H 9 O þ 4 [16,17]. The rate limiting step for proton transfer is believed to be the reorientation of water molecules, which necessarily involves breaking hydrogen bonds [16]. The concerted motion of protons along a chain of water molecules is sometimes referred to as a 'water wire'. Solvated hydroxide ions are thought to have a different mechanism for transport than hydronium ions and correspondingly demonstrate a lower ionic mobility [18]. This mechanism has been suggested to involve interconversion between the square planar H 9 O À 5 anion and the tetrahedral H 7 O À 4 anion, with the rate limiting step being the formation of the latter [16,19,20].

Modelling the silica-water interface
With regard to the silica-water interface, both ab initio and classical molecular modelling have been used to describe atomistic surface charge. Many classical forcefields capable of representing negative charges have been developed, the majority of which require a priori knowledge of the surface charge, treating surface charges as predefined and fixed throughout the simulation [21][22][23][24][25][26][27][28][29][30]. In order to study time-varying surface charge using molecular dynamics, a force field must be used which can allow bond breaking and formation. Such a forcefield is usually referred to as reactive or dissociative. The reactive force field of Rustad et al. incorporated water dissociation and led to a fully hydroxylated surface with no surface charge [31], but was not designed to accurately represent surface charging and only considered a short timeframe of 10 ps. The dissociative force field of Mahadevan and Garofalini [32] was used to study formation of silanols and transfer of protons at the surface, similarly, the Hybrid-QM/MM study of Du et al. [33] was used to model formation of silanols at the surface, however none of these models were designed to investigate protonation/deprotonation dynamics of surface silanols and do not discuss this aspect of their model. The reactive forcefield 'reaxFF', developed by Goddard et al. [34] and applied by Fogarty et al. [35], was utilised to study the silica-water interface. Their 600 ps dynamics implied that a fully hydroxylated surface was produced from a freshly cleaved slab after approximately 250 ps, the concentration of silanolate groups at the surface was not explicitly stated.
Ab initio molecular dynamics (AIMD) studies based on periodic Density Functional Theory (DFT) have investigated protonation/deprotonation involving hydronium ions in order to calculate pK a values [13,36,37] and in order to investigate dissolution mechanisms involving hydronium ions [32,38,39]. DFT has shown that hydronium ions can facilitate transfer of negative charge across the surface via the Grotthuss mechanism [13]. Mahadevan and Garofalini have noted that hydronium ions are important in short lived proton transfer processes at the surface [32], which has been supported by experimental observations [40][41][42].
Reaction of hydronium ions with silanolate groups at the silicawater interface has been shown by Leung et al. to have no barrier along the reaction coordinate based on Potential of Mean Force calculations [13,36]. In agreement with this result, Liu et al. have shown that there is no energetic barrier to acid dissociation of orthosilicic acid (SiðOHÞ 4 ) [36]. In contrast to silanol acid dissociation reactions, the counterpart basic reaction shown in Eq. (3) has received much less attention within the literature using ab initio methods for silica surfaces.
It has long been known that the dynamics and energetics of hydroxyl-based proton transfer in pure water differs from hydronium-based proton transfers [16], and therefore Eq. (3) must be considered separately from hydronium reactions (Eq. (2)) for an accurate representation of atomic interfacial proton transfer reactions. The hydroxyl transfer reaction shown in Eq. (3) has received little attention with respect to ab initio studies, one of the key exceptions being the work of Xiao and Lasaga who have investigated this reaction as a possible precursor to silica dissolution [43]. Xiao and Lasaga used SiH 3 OH and ðHOÞ 3 Si-O-SiðOHÞ 3 cluster models for the silica surface at the HF/6-31G ⁄ level (for geometries), and their results indicated that this reaction was activationless and exothermic.

Motivation
Understanding both Eqs. (1) and (3) is important from two perspectives. Firstly, adsorption of H þ and OH À is thought to be important in determining the rate of dissolution of silicates, as indicated by Xiao and Lasaga [39]. One aim of this work is to improve on the limitations of the work of Xiao and Lasaga by exploring the hydroxyl-transfer reaction (Eq. (3)) using a more representative model of the silica (periodic slab), water (solvated hydroxide molecule) and a higher level of theory via Second Order Møller-Plesset (MP2) and DFT calculations. Secondly, an atomistic model that has been shown to accurately describe both water dynamics and surface protonation/deprotonation kinetics as a function of pH does not exist. Such a model is required for fundamental understanding of double layer dynamics, for an improved understanding of the geochemical properties of oxides [44] and in order to interpret the response of charge-sensitive silica-water nanodevices such as silica nano-pore Field Effect Transistors (FETs) [45] and Ion-Sensitive and Biologically-Sensitive FETs [46,29]. A basic level of understanding of the available energetic barriers and mechanisms involved in proton transfer reactions at the surface is required before such a model can be considered. Thirdly, calculation of the transition states of these important reactions may allow a Transition State Theory description of the kinetics of the system which can be used to interpret the results of, for example, titration experiments. Understanding these reactions is important in the context of empirical macroscopic models such as SCMs in interpreting the physical significance of the reactions being modelled with empirical equilibrium constants.
In this work, DFT simulations in the form of AIMD and geometry optimisations have been used to investigate both the acid association mechanism (Eq. (1)) and the base-dissociation reaction (Eq. (3)) between silica and water. To our knowledge, this work represents the first ab initio simulation of the base-dissociation reaction which goes beyond a simple cluster-model for the silica surface, and the first ab initio simulation which explicitly explores the effect of differing solvent structure on both Eqs. (1) and (3).

Computational methods
Calculations were performed using DFT with the PBE-GGA exchange-correlation functional [47]. The PBE-GGA functional has been shown to produce accurate structures for crystalline silica [48]. The PBE functional is known to over-structure liquid water in dynamic simulations [49,50] however it can provide reasonable geometries for optimised water clusters as compared to MP2 calculations [51,52] and it has been used in various studies of the silica-water interface [13,[53][54][55].
Periodic boundary condition calculations were performed using the linear-scaling pseudopotential DFT software ONETEP version 3.5.9.8 [56,57]. PBE OPIUM [58] norm conserving pseudopotential (NC-PPs) bundled with Accelerys Material Studio 6.0.0 were utilised in all ONETEP calculations. An effective kinetic energy cutoff of approximately 800 eV was used for the psinc basis set [59], which is equivalent to the energy cutoff used in conventional plane-wave DFT codes. DFT in ONETEP was performed using self-consistent field convergence criteria whereby the RMS gradient of the NGWFs must be equal to or less than 1.8375 Â 10 À6 E h a À3=2 0 . Geometry optimisations proceeded using the BFGS algorithm until the difference in energy between iterations was equal to or less than 1 Â 10 À5 eV, 0.03 eV Å À1 and 0.001 Å for the energies, forces and maximum atomic displacement respectively. Unless otherwise stated, all calculations were performed using these settings. NWChem software version 6.3 [60] was used to perform all-electron calculations. Unless otherwise specified, these calculations were performed using the driver module, DFT and the PBE-GGA functional. All calculations used a total energy SCF tolerance of 10 Â 10 À8 E h and the aug-cc-pvtz basis set. For the data presented in Fig. 5, NWChem geometry optimisations were performed using the stepper module and a 0.05 Å maximum displacement per iteration. Example input files used and unit conversions can be found within the Supplementary Information Section 1.
In order to validate the pseudopotential used, the geometry, deprotonation energy and adsorption energy of monomeric silanol-water and silanolate-water systems were compared with allelectron calculations, the results of which are presented in Tables  1 and 2 respectively. It can be seen that there is good agreement (1-4% difference) in the calculated ONETEP NC-PPs and all-electron energies. Optimised geometries (not shown) demonstrated excellent agreement, with bond lengths within 0.01 Å and angles within 0.1°.
Born-Oppenheimer AIMD simulations were performed using ONETEP to investigate the proton transfer dynamics of three water clusters (H 3 O þ , OH À or H 3 O À 2 ). Simulations were performed using the same electronic and simulation cell settings as the slab geometry optimisations, but without any geometry constraints. A 0.5 fs AIMD timestep and the Nosé-Hoover thermostat (one chain, 8.8 fs relaxation time). The water cluster was placed above the silica surface, and 200 fs of molecular dynamics was performed. A temperature of 300 K was utilised, consistent with the AIMD of Musso et al. [53] however it should be noted that the properties of bulk water are known to be poorly reproduced without the use of elevated temperature using this functional [49].
Implicit solvation calculations were performed using ONETEP, using a self-consistent cavity and a fine grid scale of 3.0 in a 47.5 Å cubic simulation cell with open boundary conditions [62].
Similarly to the work of Leung et al. [13], the calculations reported herein treat the nuclei classically and it is assumed that the effects of zero point motion and tunnelling do not affect the qualitative nature of proton transfer mechanisms. This has been shown to be the case for electron transfer and pure water proton transport [63,64]. Table 1 Deprotonation energies DE d;gas (kJ/mol) calculated as the total energy of the deprotonated system minus the total energy of the protonated system. AE = Allelectron. NC-PP = norm-conserving pseudopotential.  Visualisation was performed using the Visual Molecular Dynamics software [65] with O-H bonds and Si-O bonds drawn of internuclear separations of less than 1.1 Å and 1.7 Å respectively. Bond distances are given in Å. In some figures, hydrogen bonds have been drawn as unlabelled dotted lines using a cutoff of 3 Å and 20°angles between hydrogen bond acceptors and donors.

Models
Musso et al. have performed a series of DFT studies on various silica polymorphs [53,66,67], including an investigation of the dynamics of water upon a (1 0 1) cleaved plane of the silica polymorph a-cristobalite. This surface is attractive from the point of view of modelling amorphous silica surface due to its surface hydroxyl density of approximately 5 OH per nm À2 , similar to fully hydroxylated silica [68]. Furthermore a-cristobalite itself has a bulk density (2.23 g cm À3 [69]) close to amorphous silica (2.20 g cm) [70]. Therefore, a-cristobalite was chosen as a model crystal structure for the DFT calculations.
The initial silica structure of a-cristobalite was obtained from the structures bundled with Accelrys Material Studio 6.0.0 which was itself generated based on a paper by Dollase [71] (primitive tetragonal P4 1 2 1 2 space group, a = b = 4.978 Å, c = 6.948 Å). A variable-cell geometry optimisation using the CASTEP software [72] was performed on the primitive cell in order to obtain relaxed unit cell-parameters for use in future calculations (a = b = 5.075 Å, c = 7.085 Å). For this calculation, a 900 eV kinetic energy cutoff was utilised with a 4 Â 4 Â 4 k-point grid and the aforementioned NC-PPs.
Using the CASTEP relaxed crystal geometry, a supercell was created from these coordinates with doubled lattice parameters, and this was optimised using the ONETEP software, which is a fixedcell dimension calculation. This produced no significant change in molecular geometry of the crystal. A (1 0 1)-plane slab of 14 Å thickness was cleaved from this crystal and passivated with a layer of hydrogen on both top and bottom, resulting in a system of 168 atoms of isolated silanol groups. The resulting lattice parameters were a = 17.431 Å, b = 10.150 Å, c = 105.929 with approximately 90 Å of this being vacuum padding. The slab was relaxed using ONETEP, with no significant rearrangement of the bulk. The optimisation resulted in a contraction of approximately 0.1 Å slab thickness. The final coordinates are shown in Fig. 1.
The optimised silica slab (Fig. 1) demonstrated isolated silanols with an OÁ Á ÁO distance of 4-5 Å and the closest OÁ Á ÁH approach distance of 4.8 Å. This result deviates from that reported by Musso et al. [53], who reported a zig-zag pattern of hydrogen bonds. However, Musso et al. comment that these hydrogen bonds are weak and disrupted at room temperature and entirely broken in the presence of water [53]. This was investigated by repeating the geometry optimisation using CASTEP, (1000 eV kinetic energy cutoff, C-point sampling of the Brillouin zone and ultrasoft pseudopotentials of Civalleri and Harrison [73]). This resulted in the same geometry as the previous ONETEP optimisation. This indicates that the deviation in structure between this work and that of Musso et al. is a result of the latter being in a different local minimum. The local minimum obtained herein provides an idealised model of a silica surface composed of isolated silanols.
Taking the neutral slab, a proton was removed from a surface silanol (indicated with a circle in Fig. 1)) to create a negatively charged silanolate group and the system was geometry optimised in ONETEP. In the protonated system the in-plane Si-O bond length was 1.638 Å and the out-of-plane (Si)-(OH) bond length was 1.643 Å, in the deprotonated system the in-plane Si-O bond was slighlty stretched (1.692 Å) and the out-of-plane Si-O À bond was shortened (1.547 Å). The geometry of the bulk and the other surface silanols were not significantly affected by the deprotonation, demonstrating that the silanols are truly isolated even in the deprotonated system. The surface charge density used in this work ($0.05 Si-O À per nm 2 ) is similar to that calculated by Behrens and Grier for a silica plate in deionised water [74], however the surface charge density of silica is highly variable depending on surface preparation, ionic strength and pH.
Unless otherwise specified, calculations were performed using 3D periodic boundary conditions using ONETEP and a vacuum gap with a neutralising background charge to minimise periodic interactions. For systems with a net charge which are also orthorhombic, it has been shown that there will be some uncompensated neutralising background charge [75] that leads to a divergent system energy, though the forces remain convergent. The ONETEP implementation of DFT has the advantage that there is little computational cost to using a large vacuum gap, therefore allowing the forces within this work to be well converged with respect to the simulation cell size to within $0.005 eV/Å.
For explicit solvent calculations, the isolated water clusters were initially optimised in vacuum. The PBE-GGA functional was found to provide a reasonable description of the ground-state geometry of simple water clusters (see Supplementary Information Section 2). The initial coordinates for the 11-water (H 3 O þ ðH 2 OÞ 11 ) and 20-water (H 3 O þ ðH 2 OÞ 20 ) hydronium clusters In order to generate a charged silica surface slab model, the highlighted silanol (black circle) was deprotonated and the system geometry optimised as described in the main text (Section 2.1). The surface has been hydrogen passivated. Silicon atoms are shown as yellow, oxygen atoms as red and hydrogen atoms as grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) mentioned in this work were obtained from the work of Hodges and Wales [76] using the Kozack-Jordan potential [77] and are candidate global minimum for a solvated hydronium ion in 11and 20-water molecules respectively. Experiments have shown the 20-water cluster is an unusually stable water cluster [78]. The optimised water geometries were placed approximately 2.3 Å distant from the optimised silica surface, and relaxation energies listed in the main text are simply calculated as the energy of the entire system (water cluster + silica) after optimisation minus the configuration they were initially placed. The bottom half of the slab was constrained during geometry optimisations, as shown in Fig. 1 using stick representation.
A summary of all model systems geometry optimised in this work can be found in Table 3.

Results and discussion
3.1. Acid-base dissociation reactions of single hydronium or hydroxide molecules with silica surface models In order to investigate the proton transfer described in Eq. (1), 200 fs of AIMD were performed on the S surf O À þ H 3 O þ system. The same initial configuration was also used in a geometry optimisation. Similarly, in order to investigate the proton transfer reaction described in Eq. (3), the model system S surf OH þ OH À was considered. As model systems for S surf , both an isolated silanol cluster (SiH 3 OH/SiH 3 O À ) and a periodic silica slab (S surf OH/ S surf O À ) was considered.
The initial and final coordinates of the geometry optimisation of the H 3 O þ system are shown in Fig. 2(a) for the slab system, and Supplementary Information Section 3 for the cluster system. Snapshots of the AIMD are shown in Fig. 3 and a video 1  The initial and final coordinates of the geometry optimisation of the hydroxide system is shown in Fig. 2(b) for the slab, and in Supplementary Information Section 3 for the cluster model. The optimisation showed a proton transfer from the Si-OH to the OH À , resulting in a Si-O À Á Á Á H 2 O hydrogen bonded system as described in Eq. (3). Snapshots of the AIMD are shown in Fig. 4 and a video 2 of the AIMD trajectory can be found within the Supplementary Information. Reorientation of the OH À occurred for the first 75 fs, followed by rapid proton transfer over the next $50 fs, after which the H 2 O diffused 4 Å away from the now negatively charged silanolate group.
Proton transfer during a geometry optimisation indicates that the initial encounter-pair is energetically unstable and that there is no activation energy to the proton transfer process. Fig. 5 shows the energy profile for a geometry optimisation performed upon a cluster system at both the PBE-GGA and MP2 level of theory, and on the periodic slab model of the silica surface (PBE). It can be seen that the total energy of the system decreases smoothly and monotonically. Consistent with this observation, using both the ONETEP and NWChem transition state search functionality, no transition state could be identified for these proton transfer coordinates.
It is possible that hydrogen bonding in geminal or vicinal silanols would introduce energetic barriers, for example Sulpizi et al. recently published a study indicating that silanols with in-plane hydrogen bonds of a hydroxylated quartz surface are 3 pK a units more acidic than those forming hydrogen bonds out-of-plane with water [37]. Leung et al. have shown that highly strained sites can be significantly more acidic [13]. Geometry optimisations using geminal and vicinal silanol cluster models (shown in Supplementary Information Section 3) indicated that these proton transfers remain activationless. This work will be restricted to the study of isolated silanols.
The results thus far presented have not taken into account the effect of solvation on hydronium and hydroxide ions, and therefore it is possible that the instability of the reactants in the above DFT studies may be a result of neglecting these interactions. It was found that the incorporation of implicit solvent also demonstrated activationless proton transfer during geometry optimisations of either the hydronium or hydroxide (Supplementary Information Section 3). The implicit solvent model cannot explicitly incorporate the effects of water cooperativity [79] and Grotthuss proton transport [15], and therefore solvation of the periodic silica slab model was investigated via explicit solvation, in terms of water clusters of increasing size placed at the silica surface.

Surface protonation in the presence of explicitly solvated hydronium
Proton transport was investigated for hydrogen-bonded water clusters at the surface via geometry optimisations of water clusters in contact with the silica surface in vacuum. The following systems were investigated: H 3 O þ ; H 5 O þ 2 (''Zundel cation''), H 9 O þ 4 (''Eigen cation''), a hydronium ion solvated in 11 water molecules (H 3 O þ ðH 2 OÞ 11 ), and a hydronium ion solvated in 20 water molecules (H 3 O þ ðH 2 OÞ 20 ). See Supplementary Information Section 2 for images of these structures in isolation.
The initial structures and geometry optimised structures of Eigen cation and Zundel cation systems are shown in Figs. 6 and 7 respectively. Both of these simulations demonstrated proton transfer via the Grotthuss mechanism. A hydrogen bonded network between the water-cluster and the surface was formed prior to the proton transfer.
The optimisation for the larger H 3 O þ ðH 2 OÞ 11 system is presented in Fig. 8. Similarly to the previous hydronium systems, the surface was protonated within the first few iterations which was followed by a rearrangement of the protons in the system to stabilise the hydroxide ion produced. Unlike the previous simulations, the resulting water cluster demonstrated some structural character of both a Zundel cation and a H 3 O À 2 dimer (Fig. 8c), indicating that there is some charge separation within the water cluster and that the surface has been protonated via the H 3 O þ anion stabilising the deprotonation of a water molecule.
A Natural Population Analysis (NPA) ( Supplementary  Information Section 4) confirmed that the Zundel-like substructure had a Natural Charge significantly more positive that the other oxygen atoms in the simulation, however the H 3 O 2 substructure did not show a particularly negative charge relative to other oxygens ( Supplementary Information Section 4). Prior to optimisation, the silanolate group surface-terminal oxygen had a Natural Charge of À1.17 and all other surface-terminal oxygens showed a Natural Charge of À1:04 AE 0:01. After optimisation all surface-terminal oxygens showed a Natural Charge of À1:04 AE 0:01, which indicates that the silanolate had been neutralised. NPA analysis showed that water cluster itself is almost neutral, with a net natural charge of +0.04 relative to the net Natural Charge of À0.04 for the silica slab. It can be concluded that the silica surface has been protonated, and the system neutralised. However the water cluster has distributed itself so as to retain a structural defect analogous to a Zundel cation. The next section uses a larger water cluster to investigate how the distance of the hydronium ion from the surface may affect this mechanism of surface protonation.
The H 3 O þ ðH 2 OÞ 20 water cluster system was studied, in which the H 3 O þ could be placed initially distant from the surface or close to the surface depending on the orientation of the cluster. Figs. 9-11 present the results of a geometry optimisation of this cluster in three different initial orientations respectively. Conformation A, shown in Fig. 9, initialised the H 3 O þ ion at a distance from the surface with the shortest path between the silanolate and H 3 O þ being four water molecules. Using chemical notation, the initial structure can be described as S surf -O À Á Á Á ðH 2 OÞ 20 Á Á Á H 3 O þ . This conformation had four water molecules in the shortest path between the silanolate and H 3 O þ . Rotation of this water cluster relative to the surface resulted in Conformation B (Fig. 10) and Conformation C (Fig. 11)) which have three and two intervening water molecules, respectively.
Geometry optimisation of this cluster in all three Conformations A, B and C showed activationless protonation of the silanolate surface. In Conformation A (Fig. 9), protonation of the surface silanolate occurred via water dissociation resulting in a substructure of   3. AIMD simulation of a neutral silica slab with a H3O þ above a silanolate group at the silica surface (Eq. (1)). The silanolate group is rapidly protonated by the H3O þ within the first 50 fs. From left to right, snapshots are shown at every 25 fs. Fig. 4. AIMD simulation of a neutral silica slab with an OH À above a silanol at the silica surface (Eq. (3)). The silanol is rapidly deprotonated by the OH À within $50 fs after an initial period of $75 fs. From left to right, snapshots are shown at every 25 fs. Proton transfer occurs rapidly within approximately 25 fs.
water near the surface similar to a solvated hydroxide ion H 3 O À 2 À Á near the silica-water interface. NPA showed that the H 3 O þ remained positive compared to the rest of the water cluster, but the H 3 O À 2 -like substructure was not negative. Using chemical notation, the optimised structure could be schematically drawn as Interestingly, the H 3 O þ ion in the cluster was unperturbed by the silanolate environment, indicating that, at least at 0 K, there is a distance beyond which H 3 O þ will not recombine with the silanolate group directly (herein referred to as the 'basin of attraction') which in this case occurs at 4 water molecules separation from the silanolate. In contrast, for the optimisation of Conformation C (Fig. 11), the initial proximity of the H 3 O þ facilitated complete proton transfer to the silanolate, resulting in a neutral system S surf -OH Á Á Á ðH 2 OÞ 21 . Conformation B showed a mechanism in between these two extremes: the optimised structure (Fig. 10)) can be seen schematically as S surf -OH Á Á Á OH Á Á Á H 5 O þ 2 Á Á Á ðH 2 OÞ 18 , and again showed a Natural Charge (Supplementary Information Section 4) which was positive for the H 5 O þ 2 substructure but not negative for the OH À -like substructure.
By comparing these three conformations, the ground-state basin of attraction for activationless proton transport for isolated silanols on silica surfaces is seen to be 2-3 water molecules. The proton transfer mechanism observed in Figs. 8-10 involves initial deprotonation of the mediating water molecules stabilised by the hydronium ion, and resembles the 'proton holes' transport discussed by [80]. That is to say, these results suggest that even if the H 3 O þ is distant from the surface, its presence is enough to stabilise the surrounding waters such that the surface can be protonated without the H 3 O þ ion losing its localised proton, as shown most clearly in Fig. 9.
Thermal fluctuations might be expected to reduce the basin of attraction by breaking the hydrogen-bonded water-wires required for Grotthuss mechanism-like proton transport, however, the thermal energy would also allow activated proton transfer, thereby increasing the size of the basin of attraction. It is interesting to compare these results with the results of the AIMD simulations  1) and (3) and the respective system energy versus geometry optimisation step is shown in figure (a) and figure (b) respectively. A minimal cluster silanol model was used, and the chemical system is drawn as insets within each figure. The initial energy (y-axis) is normalised to zero. As each optimisation took a different number of steps, for comparison the optimisation progress (x-axis) is presented, in which the optimisation has been scaled to range from the initial structure (left of x-axis) to the fully optimised structure (far right of the x-axis). It can be seen that in both optimisations there was a smooth, monotonic decrease in energy upon optimisation, indicating an activationless proton transfer. Images of each geometry optimisation can be found in Fig. 2 within the main text for the ONETEP optimisations, and in Section 3 of the Supplementary Information for the NWChem optimisations. Fig. 6. Geometry optimisation of a Zundel cation (H5O þ 2 ) above a silanolate group at the silica surface. As the geometry optimisation proceeds, the silanolate group is protonated by the cation (b), resulting in two water molecules hydrogen bonded to a silanol group (c). of the water-silica interface by Leung et al., which indicated that once the H 3 O þ was further than 2-3 water molecules from the surface, then the H 3 O þ would diffuse away without protonating the surface. This suggests that the thermal contribution and/or screening from the bulk prevents proton transfer to the surface, however we note that these AIMD simulations were run at elevated temperature in order to preserve the liquid dynamics of water under the PBE functional [13].

Surface deprotonation in presence of explicitly solvated hydroxide
In this section, AIMD simulation of the H 3 O À 2 cluster at the silica interface is presented, followed by geometry optimisations of the H 3 O À 2 ; H 9 O À 5 and H 7 O À 4 clusters at the silica interface. As discussed in the introduction, several of these clusters are of particular interest as they have been presented in the literature as important in the transfer of protons within pure water.
Snapshots of the H 3 O À 2 system AIMD simulation are shown in Fig. 12 and a video 3 is included in the Supplementary Information. From the AIMD simulation it can be seen that during the first $50 fs the symmetric structure of the H 3 O À 2 anion is broken as a H 2 O . . . OH À structure is formed, after which the silanol is deprotonated within the next $25 fs.
The geometry optimisation of the same H 3 O À 2 system is shown in Fig. 13. Similarly to the previously shown isolated hydroxyl system in Fig. 2(b), geometry optimisation of the H 3 O À 2 system Fig. 7. Geometry optimisation of an Eigen cation (H9O þ 4 ) above a silanolate group at the silica surface is shown in (a). As the optimisation proceeds a proton transfer occurs such that the silanolate is protonated by the Eigen cation (b), forming a water cluster which is hydrogen bonded to the silanol surface (c). Fig. 8. Geometry optimisation of a hydronium ion (shown in blue) solvated in 11 water molecules H3O þ ðH2OÞ 11 À Á above a silanolate group on the silica surface is shown in (a). The surface is protonated within the first few iterations, as shown in (b), forming a water network shown in (c). The water network contains a substructure which is structurally similar to a Zundel cation (blue) and H3O À 2 dimer (green) near the surface, a shown in (c). Natural Population Analysis of the optimised system shown in (c) can be found in the Supplementary Information Section 4, and demonstrates that the blue-highlighted atoms are more positively charged than the rest of the water cluster. The green-highlighted oxygens showed a similar Natural Charge to the rest of the water cluster. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) resulted in deprotonation of the silanol surface by the hydroxide cluster and the formation of a pair of water molecules hydrogen bonded to a silanolate group on the surface.
The optimisation for the H 7 O À 4 system is shown in Fig. 14   À Á above a silanolate group at the silica surface. The H3O þ is initially separated from the silanol by 3 water molecules (a). This is the same water cluster as Fig. 9, but rotated such that the hydronium ion (blue) is closer to the surface silanolate. As the optimisation proceeds the silanolate is protonated by a nearby water molecule (b) and the water cluster rearranges forming a substructure which contains what could be described as a solvated hydroxide ion (green) and a solvated hydronium ion (blue). Natural Population Analysis showed that the blue-highlighted oxygens in the optimised structure are more positive that the other oxygen atoms, as shown in Supplementary Information Section 4. The green-highlighted oxygen showed a similar Natural Charge to the rest of the water cluster. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) redistribution of electrons within the water molecules occurs which can reuslt in a reduction of hydrogen bond length with increasing cluster size [81]. It should be noted that this correlation is not guaranteed as many other factors can affect hydrogen bond lengths, for example, if the hydrogen bond conformation is nonoptimal then its strength will decrease (and thereby length increase). Hydrogen bonds at the periphery of a cluster may be expected to show a smaller change in hydrogen bond length [79]. This increased stabilisation of the SiO À Á Á Á HO-O hydrogen bond with increasing water cluster size is significant as this effect is often neglected in the parameterisation and validation of classical force-fields, which are often constructed using a single water molecule interacting with the surface [25,82]. For the square planar H 9 O À 5 system (Supplementary Information Section 3), geometry optimisation did not deprotonate the silanol group at the surface and the hydroxyl cluster remained stable and relatively unperturbed by the silica surface environment, with only a slight distortion (0.1-0.2 Å) of the hydrogen bonds normal to the surface. The lack of proton transfer indicates that a significant activation barrier to deprotonation was present for this system. This result is consistent with the AIMD of pure water by Tuckerman et al., who observed no proton transfer for the more stable square planar H 9 O À 5 complex [16,20], and the experimental and theoretical study of Cwiklik et al. on pure water, who observed that this structure is more stable than the H 7 O À 4 tetrahedral cluster [83].
These simulations extend the work presented by Xiao and Lasaga [43] and, to our knowledge, represent some of the first dynamic and mechanistic ab initio descriptions of surface charging due to solvated hydroxide at the silica/water interface. An interesting find of this study is that the SiO À . . . :HO-O hydrogen bond length is strongly dependent upon the degree of solvation. Furthermore, this study indicates that, similarly to proton transfer in pure water, deprotonation of SiOH in the presence of OH À demonstrates no significant activation energy, except in the case of the highly stable H 7 O À 4 solvated cluster.

Timescales and energetics
All proton transfer events observed in geometry optimisations occurred within the first few optimisation steps, indicating a strong energy gradient driving the reaction. AIMD results already presented (Figs. 3, 4 and 12) have shown these proton transfer events occur on a femtosecond timescale ($25-100 fs). This result is consistent with Car-Parinello molecular dynamics studies on pure water, which have shown that once a water wire is formed, Fig. 11. Conformation C: Geometry optimisation of a hydronium ion solvated in 20 water molecules H3O þ ðH2OÞ 20 À Á above a silanolate group at the silica surface. The H3O þ is initially separated from the silanol by 2 water molecules (a). This is the same water cluster as Fig. 9 and 10, but rotated such that the hydronium ion (blue) is closer to the surface silanolate. The closer proximity of the hydronium ion relative to Conformation A results in proton transfer from the hydronium ion to the silanolate group of the surface as shown in (b), ultimately producing a neutral water cluster hydrogen bonded to a neutral silanol group, as shown in (c). The intermediate hydroxide-like substructure is shown in green. Natural Population Analysis showed that oxygens in the optimised structure have similar Natural Charges, as shown in Supplementary Information Section 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 12. AIMD simulation of a neutral silica slab with the H3O À 2 anion above a silanol at the silica surface (Eq. (3)). Reorientation of the system occurs for the first $50 fs from which the silanol is deprotonated by the H3O À 2 within another $25 fs resulting in a water molecule hydrogen bonded to the silanolate group. From left to right, snapshots are shown at every 25 fs. OH À =H 3 O þ recombination occurs extremely rapidly (on a femtosecond timescale) [84,85].
In order to quantify the energetics of this proton transfer, allelectron calculations using NWChem were performed upon a silica cluster model of an isolated silanol molecule (SiH 3 OH) and an orthosilicic acid molecule (SiðOHÞ 4 ). Three different reaction schemes were considered for each reaction in which a single additional water molecule stabilised the reactants, the full details of these calculations can be found within the Supplementary Information Section 5. Each different reaction scheme considers a different combination of hydrogen bonding between the products, which can lead to large differences in the reaction energies. Reaction energies of between À637 and À682 kJ/mol for orthosilicic acid, and of between À655 to À693 kJ/mol for silanol were calculated for Reaction (1). Reaction energies of between À43.5 and À105 kJ/mol were calculated for orthosilicic acid, and of between À25.9 and À142 kJ/mol for silanol were calculated for Reaction (3). For Reaction (3), using a ðHOÞ 3 Si-O-SiðOHÞ 2 ðOHÞ 3 model of the surface, Xiao and Lasaga calculated a reaction energy of À232.6 kJ/mol at the MP2/6-31G ⁄ level, The resulting reaction energy is likely more exothermic than the silicic acid and silanol due to the formation of multiple hydrogen bonds in the resulting complex [43].
These values indicate that hydronium based protonation events are significantly more exothermic than their counterpart hydroxide deprotonation events, suggesting that the chemistry of silica deprotonation cannot be treated as simply the reverse of protonation at an atomistic scale. These calculations also demonstrate that both protonation and deprotonation reactions are highly exothermic; by comparison k b T is approximately $2.48 kJ/mol at 298 K. The fast time scale and high exothermicity of this reaction has significance for building dynamic models of surface charging, indicating that these reaction coordinates might be modelled as diffusionlimited once a hydrogen-bonded encounter pair has been formed.  13. Geometry optimisation of H3O À 2 silanol group at the silica surface. Initial structure is shown in (a). As the optimisation proceeds proton transfer from the silanol to the H3O À 2 anion occurs, as shown in (b). Image (c) shows the optimised structure, in which the silanol has been deprotonated.
Fig. 14. Geometry optimisation of the H7O À 4 anion above a silanol group at the silica surface. The initial structure is shown in (a). As the optimisation proceeded, Grotthuss transfer of the proton from the surface occurs (b) resulting in a deprotonated silanolate group in the optimised structure (c).

Conclusions
Using AIMD and geometry optimisations, we have investigated the acid-base dissociation mechanism for protonation and deprotonation events of isolated silanol and silanolate groups by hydronium ions and hydroxide ions in solution. The reaction mechanism was observed to be rapid, highly exothermic and predominantly activationless. We believe that this work is the first to go beyond simple cluster models of the surface to study surface charging due to adsorption of solvated hydroxide at the silica-water interface.
The acid dissociation of isolated silanols (Eq. (1)) did not demonstrate an energetic barrier to the proton transfer, whether in the gas phase, implictly solvated or explicitly solvated using water clusters. Simulations of the H 3 O þ ðH 2 OÞ 11 and H 3 O þ ðH 2 OÞ 20 clusters demonstrated proton transfer via the 'proton holes' mechanism in which the hydronium ion stabilises water-dissociation which, in turn, protonates the silanolate group. This mechanism has rarely been considered in the literature, but could indicate that surface protonation is possible even when the hydronium ion is distant (4 water molecules at least) from the surface.
The dissociation of isolated silanols in the presence of hydroxide (Eq. (3)) was also found to behave as an activationless process for the cases of both the gas phase hydroxide ions and the implicitly solvated hydroxide ions. For the case of explicitly solvated hydroxide ions, the local environment of hydrogen bonded silanols and waters was shown to be capable of creating an energetic barrier to deprotonation in the case of the H 9 O À 5 anion, but showed complete or partial deprotonation for the H 3 O À 2 and H 7 O À 4 hydroxide clusters. This energetic barrier to proton transfer for the H 9 O À 5 is consistent with the pure water simulations of Tuckerman et al. [20]. The Si-O À -Á Á Á H 2 O hydrogen bond length was found to be strongly dependent upon the degree of solvation, which could have significant implications for the accurate parameterisation of this bond in molecular dynamics force fields.
This work suggests that proton transfer events at the isolated silanol-water interface often do not exhibit a well-defined transition state and therefore Transition State Theory is likely inapplicable. Furthermore, the fast time scale and high exothermicity of the reactions discussed herein shows these hydronium/hydroxide systems will readily transfer protons, therefore surface scientists should be cautious when simulating such systems in the context of non-reactive forcefields.
Proton transfer was shown to involve first a reorientation of the water into hydrogen bonds with the surface group followed by proton transfer along a hydrogen bonded network of water to/from the surface group. This indicates that proton transfer for both hydroxide and hydronium ions is likely to be limited by the rate of reorientation of solvated hydroxide/hydronium clusters at the surface, as is thought to be the case for proton transfer in pure water [16,86].
In this work, the silica geometry was assumed to be regular and composed of isolated silanols, however it is possible that the silica structure (e.g. geminal or vicinal silanols) could affect energetic barriers. Furthermore, this study does not incorporate the electrostatic effect of counterions in the double layer which will introduce stabilisation to negatively charged sites and can stabilise hydroxide molecules at the surface [87].
As this study has been primarily focused upon static geometry optimisation calculations and short-timescale AIMD, it is the hope of the authors that this work will stimulate further work towards a model of the system which is capable of accurately describing the complex and dynamic nature of surface charging at an atomistic scale, without the need for surface-specific empirical parametrization or computationally expensive ab initio calculations.