The Zwitterion Isomer of Orthosilicic Acid and its Role in Neutral pH Dimerization from Density Functional Theory

The speciation and characteristics of silica in water at standard conditions are still poorly understood. Using density functional theory, the plausible existence of the zwitterion isomer of orthosilicic acid is proposed to account for some of the properties of silica in neutral pH water. Explicit hydration and explicit addition of salt are used in modeling the zwitterion and the dimerization reaction. Paths between orthosilicic acid, the zwitterion and the acid dissociation products are presented. Implicit solvation free energies and quasi-harmonic corrections are evaluated. The pK for the formation of the aqueous zwitterion species is calculated to be 6.5 in dilute silica solutions and the activation energy for the dimerization reaction ranges from 59.7 kJ/mol to 70.7 kJ/mol depending on salt concentration in agreement with experiment.


Introduction
Orthosilicic acid is the most fundamental building block of dissolved silica, which are themselves among the most ubiquitous classes of compounds found in waters of the Earth. A thorough understanding of the chemistry of orthosilicic acid is therefore vital to elucidate a wide variety of phenomena involving natural silica such as mineralization, nutrient bioavailability, and geothermal well scaling. Unfortunately, although great effort has been expended in experiments and simulation, the dimerization reaction, which is one of the simplest of its reactions, is still not very well understood; the same can be said about the reverse reaction, or the hydrolysis of the dimer.
Studies on the reaction order are conflicting and indicate that much is not known about the dynamics of silica in solution at the atomic scale. Studies at pH 4 by Alexander et al. [1] showed third order kinetics in agreement with Goto [2], who conducted work at pH 7 to 10. However, Rothbaum and Rohde [3] found the oligomerization at pH 7 to 8 to be second order and Icopini et al. [4] reported fourth order rates from pH 3 to 11. Gorrepati et al. [5] found second order rates below pH 0. The activation energies from these studies vary from 12.6 kJ/mol, to 58 kJ/mol, to 71 kJ/mol [3,6,7].
The dependence of the kinetics on ionic strength, or alternatively cation concentration, is just as puzzling [8]. Sodium ion has been shown to decrease the reaction rate of silica condensation at early stages [6] yet increase crystallization and growth at later stages [9]. Icopini et al. [4] found an exponential increase in oligomerization rates with ionic strength. Gorrepatti et al. [5] found rates that varied with the salt added according to the following relationship: AlCl 3 > CaCl 2 > MgCl 2 > NaCl > CsCl > no salt. How these phenomena mentioned above happen is unresolved. The complex relationships between rates and pH and ionic strength open the question of the existence of yet unidentified reactive species or class of species that favor neutral pH and high ionic strength.
Simulation has been providing immense insight to understanding the dimerization process, aiding in the discovery of species that would not otherwise have been easily identified by experiment, and allowing the examination of molecular scale local effects. Simulation methods have grown in sophistication in parallel with the development of computing memory and processing power. The primary emphasis in the present study is on density functional theory methods (DFT) specifically, although other schemes that are in active development and application such as molecular dynamics [10][11][12] and Monte Carlo [13][14][15] have also been applied to investigate these mechanism. These other schemes are able to model chemical systems up to the micrometer scale that would not be practicable with DFT. However, the elucidation of chemical reactivity is still principally a question that is best handled on the atomic-length scale by DFT and ab initio methods in general because of the changes in bonding, and therefore electronic structure, during the course of a reaction. At the same time, it is worth noting that there are methods that algorithmically combine the strengths of both ab initio and molecular dynamics aptly named ab initio molecular dynamics [16][17][18].
The impetus of these simulations is to model the real system as closely as computationally practical, and to incorporate all physical influences on the reaction center. In this regard, the nascent ab initio work by Xiao and Lasaga [19,20] on acid and base catalyzed hydrolysis of the dimer were attempts to draw silica solution chemistry from dry, gasphase simulations. Subsequently, the utilization of implicit solvation improved the energy calculations taking into account long range effects that the solvent imparts on the reacting ensemble [21][22][23] without significant computational overhead. However, the explicit addition of water molecules revealed that favorable reaction paths may be starkly different from the dry, gas-phase simulations, even with implicit solvation included [17,24,25] suggesting that long range interactions have a less influential role on the mechanisms than short range electrostatics and, more specifically, than hydrogen bonding. These effects are predominantly due to the ease of formation of hydrogen bonds by water molecules, which lowers energies significantly (12.6 kJ/mol - 33.4 kJ/ mol [26]), and the ability of water to form proton wires where hydrogens can concertedly move around and participate in the so-called Grotthus mechanism [27] (for example, the work of Liu et al. [28]). Furthermore, the hydroxides of H 4 SiO 4 themselves have been shown to participate in mobilizing hydrogen through similar mechanism in aqueous environments [29].
It is interesting to note that so far studies on the neutral pH dimerization have primarily focused on the reactivity of H 4 SiO 4 , and in particular on the nucleophilicity of hydroxide oxygens and the electro-positivity of the central silicon. For example, Trinh et al. [30] arrived at a one-step dimerization mechanism where the hydroxide oxygen of one orthosilicic acid attacks the silicon of another directly. However, this mechanism yielded a high activation energy of 127 kJ/mol that has been corroborated by succeeding simulation studies [28,31]. Furthermore, the mechanism does not explain why there is dependence of the kinetics on the ionic strength. It has increasingly become apparent that this mechanism cannot account for the empirical activation energies observed. This opens up the question of whether the canonical equations of orthosilicic acid speciation in water [32][33][34], or completely represent the relevant species in solution. In particular, some molecules that have appeared in simulations and that were considered reaction intermediates are more energetically stable than the acid dissociation products given by Eqs. (1) and (2). Also, these reaction intermediates are intriguing compounds in themselves from the standpoint of biochemical reactions, as well as their potential uses in industry and manufacturing.
The present work aims to accomplish three things. First, a zwitterion form of orthosilicic acid is examined in the neutral pH regime. This zwitterion was inspired by similar reaction intermediates reported by Mondal et al. [35] and Zhou et al. [36]. Second, its relevance in the dimerization reaction of silica through representative reactions is investigated. Third, the effects of explicitly adding sodium chloride to the system are evaluated. It is therefore envisaged that the dimerization proceeds as depicted in Fig. 1 and is positively influenced by the ionic strength.

Methodology
Density functional theoretical calculations were performed on constructed model systems using the Gaussian 16 suite of programs [37] running in Dell PowerEdge M620 servers. The three-parameter hybrid functional B3LYP method [38] was utilized and two types of basis sets were used: one which purely uses 6-311++G(d,p) on all atoms, and another which is a customized combination of 6-311++G(d,p) and 6-31 + G(d) on specified atom centers, henceforth designated 6-311++G(d,p){36}/6-31 + G(d), where "{36}" represents the number of atom centers to which 6-311++G(d,p) are applied. The use of the B3LYP method and the 6-311++G(d,p) basis set together is mainly because they are known to produce good estimates of thermochemistry for silicate systems [39] and because their use would facilitate the ease of comparison with previous studies [40][41][42]. The use of 6-311++G(d,p){36}/6-31 + G(d) was used for the bigger chemical systems for computing tractability, and is similar to the method that was utilized by Kar et al. [43] and Yang et al. [44].
The initial model systems, or configurations of atoms, were prepared as follows: 1. The starting configuration of the "pure-water" model simulation was constructed first. Twenty-seven water molecules were added randomly and evenly around an H 4 SiO 4 molecule. The resulting system was optimized (1) H 4 SiO 4 ↔ H 3 SiO 4 − + H pK a1 = 9.472(in 0.6 M NaCl) (2) H 3 SiO 4 − ↔ H 2 SiO 4 2− + H + pK a2 = 12.65(in 0.6 M NaCl) at 6-31 + G(d) and subsequently to the desired level of 6-311++G(d,p). 2. Another H 4 SiO 4 was added to the system above and thirty-one water molecules were added randomly and evenly around. The resulting system was optimized at 6-31 + G(d) and subsequently to the desired level of 6-311++G(d,p){36}/6-31 + G(d). The choice of the thirty-six atoms that made use of the higher 6-311++G(d,p) basis sets was done from defining a reaction center -the set of molecules with the anticipated bond-breaking and bond formation and the immediately neighboring water molecules. This is the starting configuration of the pure-water model for dimerization. 3. Atoms of Na and Cl were added to the two systems above, leading to two additional systems. The initial placements of the two atoms were arbitrary but juxtaposed to each other at an interatomic distance of 3A.
These "salt present" systems were likewise fully optimized in their respective desired basis sets.
The four resulting systems from the procedure described above were clusters of molecules that had 90, 92, 195 and 197 atoms each. Constrained optimizations were performed on the structures to determine other minima and first order saddle-points in between them; these minima and first order saddle points are referred to as "stationary points". Transition state optimizations were done with the synchronous transit-guided quasi-Newton method at the respective basis sets mentioned above, followed by vibrational frequency calculations to establish that structures are true saddle points. All minima were fully optimized and analyzed for harmonic vibrational frequencies and zero-point energies. Optimization calculations yield stationary points that are internal energy minima along a potential energy surface (PES), which is a hypothetical surface where vibrations are absent. However, even in the ground state, a chemical system would have a ground state energy or zero-point energy. Hence, internal energies need to be corrected for zero-point energy producing another surface that is elevated in energy from the PES, or a surface of 0 K enthalpy. Finally, using the optimized stationary points obtained above, corresponding solvation free energies were calculated by implicit hydration using the solvent model density (SMD) method [45].
Atomic charge distribution calculations were accomplished by Hirshfeld population analysis [46], and charge model 5 (CM5) calculations [47,48]. The more common Mulliken population analyses [49] results were not reported as it is considered sensitive to basis set size and to the choice of basis [47] and deemed less reliable than either Hirshfeld or CM5.
Corrections to the free energy when modeling the solute-solvent system are necessary to minimize errors due to entropic contributions of low frequency vibrational modes using the harmonic oscillator model. The method used in this study [50] is described as quasi-harmonic as it modifies the vibrational modes of the harmonic oscillator before evaluating energies and entropies. A cutoff frequency is chosen and that frequency is assigned to modes with frequencies that are less than the cutoff. The recommended cutoff frequency is 100 cm −1 [50].
To evaluate the effect of dispersion on the reaction energies, the stationary points of the 90 and 92 atom clusters were reoptimized at the B97D3 [51] level and frequency calculations were performed on them. Solvation free energies using the SMD method were then evaluated for the optimized structures. The quasi-harmonic correction approach described above was then applied.

Results
The B3LYP level computations in this study resulted in a total of thirteen fully optimized energy minima and ten associated intervening fully optimized first-order saddle points that connect the paths of the reactions of interest. The fully optimized minima represent reactants, products, and intermediates, and the first-order saddle points represent transition states. These stationary points are in four distinct systems each with their own elemental compositions and assigned basis sets. Relative energy changes were calculated between stationary points along a reaction path within each system only. The complete set of configurations of stationary points are listed in the Online Resource. Note that these discovered pathways represent competing routes to dimerization. Among all the discovered pathways, the fastest mechanism would determine the current best upper bounds to activation energies or, equivalently, activation enthalpies at 0 K.

Pure Water: H 4 SiO 4 + 27H 2 O System
There were two species of silica that can be identified in this system through following reaction paths and exploring the PES. As expected from numerous previous DFT and ab initio studies, one of these species is orthosilicic acid. The other species is a zwitterion isomer of orthosilicic acid, which henceforth will be referred to as orthosilicic acid zwitterion (SiO − (OH) 2 (H 2 O) + , or OSAZ). Numerous attempts to optimize configurations representing acid dissociation of the orthosilicic acid to H 3 SiO 4 − and free H + or H 3 O + did not stabilize but optimized to either OSAZ or H 4 SiO 4 .
Because of the symmetry of OSAZ, there are twelve possible tetrahedral arrangements for two -OH groups, one -O − group and one -OH 2 + -group around a silicon atom. For computational practicality therefore, only one of these combinations was arbitrarily chosen for determining the isomer to zwitterion path. A transition state connecting the two minima was found and is shown in Fig. 2; the corresponding change in energies for the elementary and net reactions are listed in Table 1.
Units are kJ/mol. All are computed in the B3LYP/6-311++G(d,p) level. The dot (•) signifies hydrogen bonding. Elementary reaction energy changes are pseudo-thermodynamic and are denoted by a double-dagger superscript ‡ (e.g. ΔU ‡ ). The "s" subscript indicates that free energies of solvation from implicit hydration by SMD are included.
Presented are the internal energy changes (ΔU), 0 K enthalpy changes (ΔH 0 ), the Gibbs free energy changes at 298 K (ΔG 298 ), and the Gibbs free energy changes at 298 K with implicit solvation (ΔG 298,s ) included. Note that the equilibrium constant K, and pK are related to the ΔG of the net overall reactions by The computed net free energies indicate orthosilicic acid is thermodynamically more favored than OSAZ. Implicit The absence of local minima representing free protons (H + ), hydronium ion (H 3 O + ) or dihydroxonium ions (H 5 O 2 + ) in the pure water system may have several possible explanations. One is that the energy minima defining these positively charged aqueous species may be too shallow and therefore difficult to find. Another is that the limited size of the system constrains the separation of H + and H 3 SiO4 − producing edge effects and preventing the optimization to a stable configuration. Lastly, the symmetry of the system and lack of polarization may not be conducive to their formation.

Salt-Present: H 4 SiO 4 + 27H 2 O + NaCl System
There were two species of silica that can be identified in this system through following reaction paths and exploring the PES. Again and as expected, orthosilicic acid is one of them. − , a representative transformation path from one dipolar complex I configuration to another was computed. The results are shown in Fig. 3 and the corresponding energies are listed in Table 1. Attempts to optimize OSAZ to an energy minimum were not successful in this system. Instead, OSAZ-like structures manifest as transition states between different dipolar complex I configurations.
To estimate the concentration of NaCl in the system being modeled, gross volume calculations were performed by defining the surface as the contour of 0.001 electrons/Bohr3 density and using a Monte Carlo integration method [52]. This method yielded volumes of 490 cc/mol and 460 cc/mol for stationary points C and D of Fig. 3, respectively, corresponding to concentrations of 2.0 M and 2.2 M of NaCl; the ionic strength is likewise 2.0 M and 2.2 M. In comparison, saturated brine is 6.14 M NaCl and seawater is 0.469 mol/kg NaCl [53]; the ionic strength of saturated brine is therefore >6.14 M while seawater has a known ionic strength of about 0.7 M. There were small but noticeable migrations of sodium and chlorine from the original 3A during optimization; optimized Na-Cl distances were 5.125A for the cluster with orthosilicic acid (Fig. 3 C) and 5.251A for the cluster with dissociation ( Fig. 3  D). This relative ease of sodium and chlorine migration during optimization is consistent with the interpretation that they do not hydrogen bond with water. Note that the concentration of H 4 SiO 4 are also 2.0 M and 2.2 M.
The absence of local minima corresponding to OSAZ suggests that the dipole field introduced by NaCl discourages its formation and encourages acid dissociation instead. Relatedly, the uneven coulombic field and the general asymmetry of the system could be factors contributing to this absence.
The negative activation enthalpy in Eq. (B2) of Table 1 suggests that there are precomplexes that have not been found and that therefore the reaction can likely be broken down into more elementary steps. The free energy governs reactions in solution and the free energy surface may have a very different topology from enthalpy surface at 0 K, especially when solvent is playing a role.

Dimerization in Pure Water: 2H 4 SiO 4 + 60H 2 O System
Four species of silica were identified in this system through following the reaction paths and exploring the PES. Namely,  7 . Hence, there are potentially at least six reaction paths that can be explored in this system. The focus however is to find the relevance of OSAZ to the dimerization reaction and to the other three species, and therefore not all reaction paths were investigated. The results are shown in Table 2.
Units are kJ/mol. All are computed at the B3LYP/6-311++G(d,p){36}/6-31 + G(d) level. The dot (•) signifies hydrogen bonding. The asterisk '*' indicates a reassignment of the 36 atoms in the basis set. The "s" subscript indicates free energies of solvation from SMD have been included. Table 2 and Fig. 4 show the possible mechanism to dimerization under neutral pH in pure water. The process begins with two H 4 SiO 4 molecules sitting side-by-side. Next, one of the H 4 SiO 4 molecules transforms to OSAZ. Subsequently, H-bonding draws the two different silica molecules closer together. Finally, the H 4 SiO 4 molecule simultaneously donates a proton, performs a nucleophilic attack on the silicon of OSAZ, and causes the H 2 O group of OSAZ to eject. Figure 5 shows the formation of H 3 SiO 4 − from OSAZ. Note that the H 3 SiO 4 − species is charge balanced by the immediately neighboring H 3 O + formed from the freed proton (stationary point E of Fig. 5) and will be referred to as "dipolar complex II" from hereon. A second configuration where H 3 O + is more distant (stationary point F of Fig. 5) but much higher in energy was also found. The energy profiles corresponding to the stationary points in Figs. 4 and 5 are shown in Fig. 6.
Because the bond breaking and bond formation during H 3 SiO 4 − production (stationary points C, E and F of Fig. 6) involve atoms that were not among the thirty-six that were initially assigned the high basis set, a test on the reliability of the results was performed by reassignment of basis sets on atom centers. This was also an opportunity to see if the placement of diffuse functions on hydrogens that were available in 6-311++G(d,p), but not in 6-31 + G(d), would change outcomes significantly. Optimization and frequency calculations were performed producing a second, generally similar energy curve, shown in blue (Eqs. A7 to A10 in Table 2) and to nearly identical configuration geometries, based on RMSD of 0.020A-0.028A; here, the RMSD was determined by Kabsch algorithm [54,55]. Eq. (A9) in Table 2 was the most sensitive to this reassignment and change in basis, but exhibiting a free energy difference of only 3.1 kJ/mol.
Net reactions and thermodynamic quantities are shown in Table 3. It should be emphasized that because of the   The formation of either the acid dissociation products and OSAZ were less favored than orthosilicic acid. In summary, the free energies reveal the following neutral pH favorability relationships within the vicinity of the reaction zone: H 6 Si 2 O 7 > H 4 SiO 4 > OSAZ > H 3 SiO 4 − . The negative activation enthalpies again suggests that there are precomplexes that have not been found and that therefore the reaction can be broken down into more elementary steps. Note that these negative values occur where barriers are often low and the energy surface is locally relatively flat.

Dimerization with Salt Present: 2H 4 SiO 4 + 60H 2 O + NaCl System
As with the dimerization in pure water system, four species of silica were identified in this system through following the reaction paths and exploring the PES: orthosilicic acid, OSAZ, H 3 SiO 4 − , and H 6 Si 2 O 7 . Similar to the previous section, only the dimerization reaction of OSAZ is pursued. The calculated paths are shown in Table 2. Table 2 and Fig. 7 show one possible mechanism to dimerization in neutral pH with NaCl present. The process begins with two H 4 SiO 4 molecules sitting side-by-side. Next, one of the H 4 SiO 4 molecules transforms to OSAZ. Finally, the H 4 SiO 4 molecule simultaneously performs a nucleophilic attack on the silicon of OSAZ, transfers a proton via concerted reaction to the OSAZ oxygen, and causes the H 2 O group of OSAZ to eject. Figure 8 shows the formation of H 3 SiO 4 − from OSAZ. Note that the H 3 SiO 4 − species is charge balanced by a distant H 3 O + formed by concerted transfers from the freed proton (stationary point J of Fig. 8) and is analogous to stationary point F of Fig. 5. The energy profiles corresponding to the stationary points in Figs. 7 and 8 are shown in Fig. 9. Neither dipolar complex I nor dipolar complex II was successfully optimized in this system.

Fig. 5 Acid dissociation in pure water
Although the production of H 3 SiO 4 − involved bond breaking and formation between atom centers that were not initially assigned the higher basis set, and no reassignment of atom centers to the higher basis set was deemed necessary. This was decided because the reassignment procedure in the previous section (pure water case) did not significantly modify the reaction direction or the optimized geometries and, therefore, it is surmised that the additional computational effort of altering the atom centers of the basis sets was unnecessary as it would not change the general findings regarding OSAZ reactivity.    Computing the molar volumes using the method of Parsons and Ninham [52] yielded volumes of 1033 cc/mol and 1069 cc/mol for stationary points G and I of Fig. 7 [56]. Note that although NaCl concentration is approximately half that of the system in Section 3.2, the shorter optimized interatomic distances seem to suggest that there is an influence being exerted by the higher concentration of H 4 SiO 4 on the sodium and chlorine atoms. This needs a more thorough investigation, however.
The free energies indicate that the dimerization reaction (Eq. B1 in Table 3) is thermodynamically favored but inclusion of implicit hydration decreases the favorability. In summary, the free energies reveal the following neutral pH favorability relationships within the vicinity of the reaction

The Orthosilicic Acid Zwitterion
Configurationally, OSAZ is what would result if one of the hydrogen atoms from a hydroxide group was moved to another hydroxide in orthosilicic acid resulting in a molecule with an -O − group and a -H 2 O + group. This is essentially the reaction in Fig. 2, where the movement of proton occurs through concerted transfers facilitated by surrounding water molecules.
A quick examination of the Lewis formal charges of OSAZ (Fig. 2) yields a positive charge on the -OH 2 group and a negative charge on the -O group. This doubly charged molecule is the zwitterion state. To confirm the zwitterion characteristics of OSAZ, Hirshfeld and CM5 calculations were conducted on stationary points A and B of Fig. 2. The results are shown in Table 4. Upon transformation to OSAZ, the -OH 2 group assumes a positive charge and the -O group takes on a relatively high negative charge.
Note that the magnitude of the charges differs between the Hirshfeld and CM5 results for the -OH 2 group. The CM5 method is considered a more reliable estimate than the Hirshfeld method [47] and what is crucial is the direction of change. The sums of the charges are non-zero and slightly negative, whereas the molecule is expected to be neutral. This is likely due to a numerical artifact because the electron density distribution disperses the charges throughout the 90-atom system.
While the recognition of the zwitterion characteristics of OSAZ is new, similar structures have been seen in previous work [35,36]. Table 5 2 ) with key differences. In OSAZ, four oxygens are covalently bonded around a silicon atom in a distorted tetrahedron, whereas in metasilicic acid, three covalently bonded oxygen atoms are trigonal planar with silicon [57]. Interestingly, metasilicic acid is observed in gas phase ab initio studies, whereas OSAZ-like structures are observed when water molecules are explicitly added to simulations [35,36]. Ab initio studies of gas phase metasilicic acid show a H 2 O-Si bond that is significantly longer (~1.9A) compared to OSAZ (~1.7A) [35]. These seem to suggest that relaxation of the structure to a more tetrahedral geometry from a trigonal planar one is facilitated by the hydrogen bonds provided by water molecules surrounding the molecule. Furthermore, metasilicic acid is not known to experimentally exist in the aqueous phase, although it has been identified in steam [58].

Zwitterion and Dipolar Complex Formation: H 4 SiO 4 + 27H 2 O ± NaCl Systems
The models in these systems, where a lone H 4 SiO 4 interacts with its aqueous environment, are intended to represent silica in dilute concentrations. Hence, the results of the calculations (Tables 1 and 2) suggest that, in dilute silica concentrations, the zwitterion is more likely to be observed as a stable species where salts are in low concentration while dipolar complex I is more likely to be observed in elevated salt concentrations. Furthermore, the inclusion of implicit solvation increased the free energy to form zwitterion in pure water by 5.2 kJ/mol, while it increased the free energy  to form dipolar complex I by only 0.3 kJ/mol. These suggest that increasing the polarization of the environment, by explicit ions present and by implicit solvation, decreases the favorability of OSAZ formation. (The preference is reversed in environments of higher concentrations of silica, and will be discussed in the next section.) Observe that the formation of OSAZ should not depend on pH ( Table 1).
The results predict that the acid dissociation products at 2 M NaCl and very dilute silica concentrations is dipolar complex I, where the hydronium ion H 3 O + sits one water molecule away from the H 3 SiO 4 − ion. This configuration is consistently reproducible, for example, in stationary points D and D′ of Fig. 3. The existence of H 3 O + is corroborated by theoretical [59] and empirical studies [60,61].
The pK a of 6.7 for the Eq. (B "Net reaction") in Table 1 is much lower than values for the acid dissociation of H 4 SiO 4 reported by Sjöberg et al. [34] (pK a1 = 9.472) and similarly by Fleming and Crerar [62] (pK a1 = 9.687). However, the pK a value is in close agreement to that derived from amorphous silica titration experiments by von Schindler and Kamber [63] (pK a1 = 6.8) and near the calculated values of -OH groups on the surface of quartz [64] (pK a = 5.5). These relationships need further investigation especially since different sites in silica appear to have different acidities [65,66]. The dipolar complex may be another silica species in dilute solutions whose concentration is three orders more than the fully ionized H 3 SiO 4 − . Observe that the pK of the formation of the zwitterion (pK = 7.8-8.7) is within these ranges and that therefore the concentration of zwitterion is only about one to two orders of magnitude less than dipolar complex I at these concentrations of silica.
The kinetics of the reactions from the Gibbs free energies of activation can be estimated by using transition state theory (TST). The TST rate constant is given by where k B is Boltzmann's constant, R is the gas constant, T is the temperature and h is Planck's constant. The rate constant corresponding to the forward reaction, Eq. (A1) in Table 1, with a free energy of 56.8 kJ/mol at 298 K, is 680 s −1 indicating an equilibration time of 1.5 milliseconds. Likewise, the rate constant of the reverse reaction, Eq. (A2) in Table 1, with a free energy of 12.5 kJ/mol, is 4.0 × 10 10 s −1 , indicating an almost instantaneous equilibration time. Rate constant calculations for the acid dissociation reactions, Eqs. (B1) and (B2) in Table 1, are of similar magnitudes and therefore the reaction completion times are expected to be very rapid. It would be worthwhile to mention here, however, that although TST is useful for estimating rates and for drawing comparisons, TST was developed as a gas phase theory and extending its application to systems in solution is not as straightforward. Truhlar and Pliego [67], and references found therein, discuss this topic and adaptations of TST that are more suitable for reactions in the liquid phase.

Dimerization by Zwitterion: 2H 4 SiO 4 + 60H 2 O (± ) NaCl Systems
The models in these systems, where two H 4 SiO 4 molecules are mutually interacting in an aqueous environment, are intended to explore the dimerization reaction, and therefore to investigate energetics of configurations where the local concentration of silica is sufficient to enable encounters between H 4 SiO 4 molecules. This scenario can occur at super-saturation, at the water-mineral interface, or during chance encounters due to diffusion.
The results of the calculations show that, OSAZ, dipolar complexes and H 3 SiO 4 − simultaneously occur with or without the presence of NaCl in the presence of a neighboring H 4 SiO 4 . This is in contrast to the results of the "dilute" system in the previous section where the occurrence of OSAZ is exclusive of the occurrence of dipolar complex I, and vice versa. Another, and more significant, difference is that in the dilute system, acid dissociation is favored over zwitterion formation, whereas the reverse is true in the presence of another H 4 SiO 4 . As seen in Table 3, the addition of implicit solvation in the pure water case makes the formation of H 6 Si 2 O 7 , OSAZ and H 3 SiO 4 − more favorable than without implicit solvation. On the other hand, the addition of implicit solvation in the salt present case makes their formation less favorable. Hence, at least for the formation of OSAZ and H 3 SiO 4 − , the presence of another H 4 SiO 4 seems to reverse the free energy trend seen in Table 1 on adding implicit solvation.
The dimerization mechanism in this study proceeds differently from previous mechanisms proposed. Instead of forming a chemically reactive nucleophile such as H 3 SiO 4 − , or forcing two monomers to bond together, the conversion to the zwitterion isomer turns one monomer into a target for nucleophilic attack by a neighboring H 4 SiO 4 . This mechanism can be more easily extended to oligomerization since it only requires a dangling -OH group from the other reactant (Fig. 1). The relative rates between pure water and salt present cases that are expected from their activation energies are in agreement with the observed enhancement of rates with ionic strength.
The rate determining step for the neutral pH dimerization process has a 0 K activation enthalpy of 70.7 kJ/ mol for the pure water system and 59.7 kJ/mol for the salt-present system (Eqs. A1 and B2 in Table 2), which strongly suggests that the zwitterion pathway is a good candidate for the dimerization mechanism in neutral pH. According to experimental studies [3,7], the activation energy for this reaction should be between 12.6 kJ/mol and 71 kJ/mol and therefore, the results of this study are in excellent agreement with experiment. In comparison, the mechanism previously reported by Trinh et al. [30] has an activation energy of 127 kJ/mol for the rate-determining step, while that of Liu et al. [28] has an activation energy of 133 kJ/mol.
It might be argued that OSAZ is an insignificant species by virtue of its low concentration upon equilibrium. However, its pH neutrality, its dependence on ionic strength, its dimerization activation energy, and its thermodynamic favorability in the reaction zone all strongly suggest that OSAZ is the key reactive species in the dimerization mechanism and therefore conditions affecting its formation influences the rate of dimerization. Hence, regulating its concentration influences the dimerization mechanism.
It is yet to be demonstrated that this can be extended to oligomerization in general.
Note that in the pure water case, the occurrence of OSAZ is thermodynamically slightly more preferred over H 3 SiO 4 − (Eqs. A2 and A4 in Table 3) and this preference is enhanced when NaCl is present (Eqs. B2 and B3 in Table 3). These local thermodynamic preferences are crucial in determining which mechanism predominates because while the slow step in the OSAZ mechanism is the formation of OSAZ itself, the slow step in the next competing mechanism is the attack by H 3 SiO 4 − on another monomer [19,30] and is therefore directly dependent on the concentration of H 3 SiO 4 − . This is why the latter mechanism predominates in higher pH ranges.
Equations (A1) and (B1) in Table 3 indicate that the dimer is slightly thermodynamically more favored than two monomers. The direction of reaction favorability is similar to that found by Noguera et al. [68], who computed pK 303K = −0. 69. Some workers [34,69] have however implied that the monomer should be slightly more preferred. This deserves further exploration in the future.
Comparison with the experimental values of the reverse mechanism, or hydrolysis, is also interesting. The activation energies of the rate-determining step in this study are 126 kJ/mol for the pure water system and 111 kJ/mol for the salt-present system (Eqs. A6 and B4 in Table 2), and are slightly lower than the results of Trinh et al. [30] and Xiao and Lasaga [20], who computed activation energies of 137 kJ/mol and 119 kJ/mol, respectively. However, according to experimental studies [70][71][72][73], the activation energy for the hydrolysis reaction should be between 67 and 92 kJ/mol, and therefore, the related activation energies in this study (Eqs. A6 and B4 in Table 2) are still too high by at least 19 kJ/mol. Thus, it is for future work to show that the D → [CD] ‡ (pure water) step and the I → [HI] ‡ (salt-present) step can be broken up into smaller lower activation energy steps.

Improvements to the Free Energy
The previous sections discuss energy calculations at the B3LYP level from frequency calculation within the harmonic-oscillator model. The current section discusses refinements to the free energy by quasi-harmonic treatments and using a density funtional theory method that includes dispersion. This will be an opportunity to evaluate the qualitative findings of this study, and examine possible ways to enhance the quantitative results.

Quasi-Harmonic Correction
The computed [50,74] quasi-harmonically corrected free energies of all the previously discussed reactions are shown in Table 6. Included are their values with their corresponding implicit solvation free energy added.
Reactions in A and B were computed at the B3LYP/6-311++G(d,p) level and C and D at the B3LYP/6-311++G(d,p){36}/6-31 + G(d) level. The dot (•) signifies hydrogen bonding. The asterisk '*' indicates a reassignment of the 36 atoms in the basis set. The "s" subscript indicates free energies of solvation from SMD have been included. The "t" and "o" subscript stand for inclusion of t-correction and o-correction, respectively.
Observe that the correction changes the former values (Tables 1, 2 and 3) only slightly. The maximum difference in free energies is a decrease of 8.5 kJ/mol (Eq. D3 of Table 6), which corresponds to the dimerization step. Almost all corrections decreased the former values except Eqs. (B3, C1, C2, C8, C12, and D6).
If it is to be assumed that negative free energies of activation, apart from already discussed cases in Section 4.3, are the result of overcorrection, then it can be seen that overcorrection is evident in the following applications: implicit solvation (Eq. A2 of Tables 2); quasi-harmonic correction (Eqs. B2 of Table 6); and implicit solvation combined with the quasi-harmonic correction (Eqs. C2 of Table 6). It is conceivable but yet to be demonstrated that optimizing the explicitly hydrated geometries within implicit hydration constraints would minimize or eliminate these presumed overcorrections.

Dispersion
The B3LYP optimized stationary points in the "dilute" cases were re-optimized to the B97D3 level. Subsequent frequency calculations and implicit solvation was applied to the re-optimized geometries. The configurations of stationary points are listed in the Online Resource. The  58.5 49.9 quasi-harmonic correction was applied and the resulting free energies are shown in Table 7.
All were computed at the B97D3/6-311++G(d,p) level. The dot (•) signifies hydrogen bonding. The "s" subscript indicates free energies of solvation from SMD have been included. The "t" subscript stands for inclusion of the quasi-harmonic correction.
All activation and equilibrium free energy changes decreased relative to those in B3LYP. The maximum difference was a decrease by almost 40 kJ/mol (Eq. B1 of Table 7). Implicit solvation decreased activation free energies with the exception of the same reaction.
If it is again to be assumed that negative free energies of activation are a result of overcorrection, then it can be seen that overcorrection happens in the application of implicit solvation on the explicitly hydrated clusters (Eq. B2 of Table 7). Again, it is conceivable that optimizing within implicit hydration constraints would minimize or eliminate the presumed overcorrections. These issues would be left for future studies to resolve.
Taking 37 kJ/mol as an estimate for the free energy of forming the zwitterion from H 4 SiO 4 in dilute conditions yields a pK of 6.5. Quasi-harmonic corrections yield values lower than these. In a similar manner, taking 33.5 kJ/ mol and 18.7 kJ/mol as estimates of the free energy of acid dissociation of H 4 SiO 4 in dilute conditions yield pK of 5.9 and 3.2, respectively; the former value is in good agreement with some experiments [63,65] and the latter value has been observed in simulations [75]. It is yet to be investigated at this level of density functional theory if the zwitterion can be found in the dilute condition with NaCl present and, likewise, if the acid dissociation products could be found in the pure water condition.

Conclusion
The speciation of silica dissolved in water is almost universally acknowledged to be predominantly defined by Eqs. (1) and (2), and hence, the recognized main species of monomeric silica are H 4 SiO 4 , H 3 SiO4 − and H 2 SiO4 2− . The results of this study suggest that there is at least a third reaction that needs to be considered, and that therefore there is a fourth monomeric silica species, SiO − (OH) 2 (H 2 O) + or OSAZ, which is present in solution in amounts significant enough to influence silica chemistry.
The kinetic and thermodynamic simulations suggest that the dimerization process in neutral pH regime is facilitated by this fourth species. In the reaction zone, the formation of OSAZ was found to be energetically more favored than H 3 SiO 4 − . Thus, it is proposed that the neutral pH dimerization rate depends on the formation of OSAZ, whose concentration is influenced mainly by ionic strength. Only upon the elevation of pH where the concentration of H 3 SiO 4 − attains a sufficient level that the basic mechanism [17,19] takes over.
Acknowledgements All computations in this study were performed on Yale high-performance computing servers and the author wishes to thank Mark Gerstein and Yale University for the use of these. The author thanks Prashant Emani and Declan Clarke for proofreading the text and giving invaluable comments. The author also thanks Catherine Dubois for help in editing the references. The author wants to express his appreciation to the two anonymous reviewers for their helpful comments and insightful suggestions.   18.7