The surface reactivity and hydrophilicity of silicate materials are key properties for various industrial applications. However, the structural origin of their affinity for water remains unclear. Here, based on reactive molecular dynamics simulations of a series of artificial glassy silica surfaces annealed at various temperatures and subsequently exposed to water, we show that silica exhibits a hydrophilic-to-hydrophobic transition driven by its silanol surface density. By applying topological constraint theory, we show that the surface reactivity and hydrophilic/hydrophobic character of silica are controlled by the atomic topology of its surface. This suggests that novel silicate materials with tailored reactivity and hydrophilicity could be developed through the topological nanoengineering of their surface.

Understanding the chemical reactivity and hydrophilicity/hydrophobicity of silicate glasses in aqueous conditions is of great concern for several applications, including nuclear waste immobilization,1,2 laboratory glassware,3 bioactive glasses for bone regeneration and drug packaging,4 supplementary cementitious materials,5,6 or outdoor glasses exposed to atmospheric weathering (e.g., photovoltaic, automotive, and architectural applications).7,8 In particular, understanding the surface reactivity of glassy silica—an archetypal model for more complex multicomponent silicate glasses—is key for various technological applications.9,10

Various factors have been noted to potentially alter the chemical reactivity and hydrophilicity of silica. First, the presence of even small amounts of impurities can greatly affect silica’s surface reactivity.11,12 The hydrophilicity of silica can also be tailored via polymer deposition.13 Further, the roughness of silica surfaces can have a significant impact on their degree of hydrophilicity.14 In addition, although glassy silica is typically hydrophilic under ambient conditions, its surface becomes markedly hydrophobic at elevated temperature (e.g., around 700-800 °C).15 This has been interpreted from the fact that, at low temperature, glassy silica exhibits a high surface density of Si–OH silanol groups,16–18 which, in turn, enhances the physical adsorption of water molecules on the surface.19,20 In contrast, at elevated temperature, the surface silanol groups tend to react with each other to form Si–O–Si siloxane groups,16,17 which are intrinsically hydrophobic.15 This suggests that the structure of silica’s surface is closely linked to its hydrophilicity.

Although finely characterizing the atomic structure of silica is experimentally challenging, atomistic simulations have been thoroughly used to investigate the hydroxylation of silica.21–29 However, several questions remain unanswered.8,30,31 What is the relationship between the atomic structure of silica’s surface and its aqueous reactivity? How does the surface density of silanol groups control silica’s affinity for water? Can the hydrophilicity of silica be tailored by nanoengineering the atomic structure of its surface?

In glasses, understanding such linkages between structure and engineering properties is usually challenging due their complex disordered atomic structure. Specially, the complexity of silicate glasses renders it challenging to discriminate the structural features that control reactivity to the first order from those that only show a second order effect. This can be conveniently achieved by topological constraint theory (TCT or rigidity theory), which reduces complex atomic networks into simpler structural trusses,32–35 thereby capturing the important atomic topology while filtering out less relevant structural details that only weakly impact certain macroscopic properties. In this framework, atomic networks are described as nodes (the atoms) that are connected to each other through mechanical constraints (the chemical bonds). Such constraints comprise the radial two-body bond-stretching (BS) and angular three-body bond-bending (BB) bonds, which maintain the bond distances and bond angles fixed around their average values, respectively.36,37 As per Maxwell’s criterion, the mechanical stability of a network depends on the balance between the number of constraints and degrees of freedom.38 Namely, atomic networks can be flexible, stressed–rigid, or isostatic if the number of constraints per atom (nc) is lower, larger, or equal to 3 (i.e., the number of degrees of freedom per atom in three-dimensional networks). Interestingly, recent studies have suggested that the dissolution rate of silicate glasses and crystals is controlled to the first order by their atomic topology, that is, it exponentially decreases with nc.5–7,39–44

Here, based on reactive molecular dynamics simulations, we generated a series of artificial glassy silica surfaces annealed at different temperatures and, hence, exhibiting different structures. The surfaces were subsequently exposed to water to investigate the relationship between the structure and surface reactivity. Note that similar simulations (although using a different potential parametrization) were performed by Rimsza et al.28 However, here we use the present simulations to identify the linkages between the surface structure, thermodynamic stability, reactivity, and hydrophilicity. We show that higher annealing temperatures result in more stable and less reactive surfaces. This leads to a hydrophilic-to-hydrophobic transition that is driven by the surface density of silanol groups. By performing a topological analysis of the atoms present at the surface, we show that both reactivity and hydrophilicity are controlled by the atomic topology of the surface of the glass.

Investigating the chemical reaction between glassy silica and water requires the usage of a reactive potential. To this end, we used the ReaxFF potential45 implemented through the user-reaxc package in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).46 Unlike conventional classical force-fields, the interatomic energy terms used within ReaxFF depend on a dynamic bond order determined by the local environment of each atom.47 Hence, ReaxFF is able to capture the dynamics of bond-breaking and bond-forming, on which chemical reactions rely on. Additionally, unlike classical interatomic potentials, the charge of each atom is dynamically assigned via a charge equilibration (QEq) method.48 For instance, this permits to dynamically adjust the charge of O atoms based on their local environments, e.g., bridging- or non-bridging oxygen (BO or NBO). All of these features render ReaxFF ideal to simulate the chemical reaction between glassy silica and water, while avoiding the computational cost of ab initio simulations.49 More details about the ReaxFF framework can be found in previous references.45,50

Here, we used the ReaxFF potential parametrized by Pitman et al.51,52 The potential parameters can be found in the supplementary material. Note that an alternative parametrization of ReaxFF from Yeon and van Duin53 has recently been shown to offer a more accurate prediction of the activation energy for the scission of Si–O bonds.49 However, we found this parametrization to yield a questionable description of glassy silica. First, we observed that Yeon’s parametrization results in the formation of a high fraction of 5-fold coordinated Si atoms (around 5%, in agreement with previous simulations28). Although the existence of 5-fold coordinated Si defects in silicate glasses has been suggested experimentally,54 such defects were noted to be extremely rare in most silicate glasses (with a maximum concentration of 0.16% in a fast-quenched potassium silicate glass54). Further, no 5-fold coordinated Si atoms were observed in glassy silica at ambient pressure.54 Hence, the defect concentration predicted by the parametrization from Yeon et al. appears unrealistically high. In addition, we observed that Yeon’s parametrization yields unrealistic stiffness values—we computed a bulk modulus of 602 GPa for α-quartz, which is about 20 times larger than the experimental value.55 In contrast, we recently showed that the parametrization from Pitman et al.51,52 yields an excellent description of glassy silica.50,56 The atomic structure was found to be in good agreement with the available experimental structure and did not comprise any 5-fold coordinated Si atoms (also see Fig. 1). In addition, the predicted stiffness was found to be realistic—with a computed Young’s modulus of 80 GPa (as compared to the experimental value of 73 GPa).50 Finally, the Pitman parametrization of the ReaxFF potential was found to yield realistic values for the diffusion constant of O atoms and the activation energy of diffusion, which suggests that the dynamics of silica melts are well reproduced.50 

FIG. 1.

Total correlation function of bulk glassy silica computed with ReaxFF. The results are compared with neutron diffraction data.64 

FIG. 1.

Total correlation function of bulk glassy silica computed with ReaxFF. The results are compared with neutron diffraction data.64 

Close modal

A bulk silica glass was generated based on the common melt-quench method,50,57–59 as detailed in the following. First, 3000 atoms were randomly placed in a large cubic simulation box (with a length of 50 Å) while ensuring the absence of any unrealistic overlap. The system was then relaxed at 4000 K for 100 ps to ensure the loss of the memory of the initial configuration. The melt was subsequently cooled from 4000 to 300 K with a cooling rate of 1 K/ps. Finally, the glass was further equilibrated at 300 K for 1 ns. The NPT ensemble, ensuring a zero pressure, was applied throughout the process. We obtained a final density of 2.18 g/cm3 at 300 K,50 which compares well with the experimental value (2.2 g/cm3).60 Note that, although the cooling rate used herein is much higher than those achieved experimentally, the structure of silica has been shown to weakly depend on the cooling rate.61–63 However, we recently noted that the use of cooling rates larger than 1 K/ps can yield unrealistic structural defects.61 As shown in Fig. 1, the structure of glassy silica predicted by ReaxFF is in excellent agreement with neutron diffraction data.50,64

Starting from the bulk glassy silica sample previously generated, two surfaces were created by instantly elongating the simulation box along the z-axis by 15 Å on both sides, thereby resulting in the creation of void spaces above and below the simulated system. Such a distance ensures that the two surfaces do not “see” each other. The system was subsequently equilibrated at 300 K for 1 ns in the NVT ensemble (i.e., constant number of atoms N, volume V, and temperature T). No significant variations in the density of silica were observed upon equilibration. Next, the system was independently annealed at 700, 1000, 1300, and 1600 K for 1 ns in the NVT ensemble—to preserve the void space in between the two surfaces. Note that these temperatures remain low as compared to the simulated glass transition temperature (around 3000 K58). In all cases, this duration was found to be long enough to achieve the convergence of the pressure and energy of the system.

For each annealed silica system, a silica–water interface was generated by inserting some water molecules into the void space with a density of 1 g/cm3 (see Fig. 2). The silica–water system was first subjected to an energy minimization and subsequently allowed to react at 300 K for 1 ns. During this process, the lateral dimensions (x and y) were kept fixed, while the NPT ensemble was used on the z-direction to keep a zero pressure.

FIG. 2.

Illustration of the atomic configuration of a silica–water system. Si, O, and H are represented in yellow, red, and white, respectively.

FIG. 2.

Illustration of the atomic configuration of a silica–water system. Si, O, and H are represented in yellow, red, and white, respectively.

Close modal

Starting from the basic structure of bulk silica—wherein tetrahedral Si units are connected to each other through their four corners, i.e., through BO atoms—the artificial dry silica surfaces simulated herein exhibit various structural defects, including 3-fold coordinated Si atoms (Siiii), terminating NBO, and edge-sharing (ES) Si units (i.e., two-membered Si ring),9,18,53 in agreement with previous simulations.28 It should be noted that, although they might exist in high-vacuum conditions, such completely dry surfaces (i.e., non-hydroxylated) are unlikely to be observed in realistic conditions of humidity since the silica surface will almost instantly get hydroxylated.

We observe here an equal concentration of NBO and Siiii defects, as the breakage of a Si–O bond results in the simultaneous formation of both of these defects (i.e., one on the upper surface and the other on the lower surface). Note that, since no structural defects are observed in pristine bulk silica, the surface density of defects was calculated here from the total number of defects divided by the lateral cross section of the simulation box and, hence, does not depend on any arbitrary definition of the thickness of the surface. As shown in Fig. 3(a), the surface density of NBO decreases monotonically with increasing annealing temperature. We observe that the surface densities of NBO and 3-fold coordinated Si atoms remain equal to each other upon annealing. This suggests that thermal annealing partially allows the atoms to locally reorganize and relax some of the structural defects by mutual recombination of NBO and Siiii defects, in agreement with previous results.20,26,28,65 This reorganization partially induces the formation of ES units. We observe that such ES units may form when a 3-fold coordinated Si atom is initially directly connected to a Si atom comprising a NBO. The NBO then tends to create a bond with the undercoordinated Si, thereby resulting in the formation of a 2-membered ring structure. As shown in Fig. 3(b), the surface concentration of ES units is found to increase with increasing annealing temperature.

FIG. 3.

Surface density of (a) non-bridging oxygen (NBO) atoms (or 3-fold coordinated Si atoms) and (b) edge-sharing (ES) units as a function of the annealing temperature. The lines are guide to the eye.

FIG. 3.

Surface density of (a) non-bridging oxygen (NBO) atoms (or 3-fold coordinated Si atoms) and (b) edge-sharing (ES) units as a function of the annealing temperature. The lines are guide to the eye.

Close modal

The surface energy of the silica samples was estimated from the variation in the potential energy of the system upon the creation of the surface and subsequent thermal annealing (divided by the created surface area). In all cases, the potential energy was determined after energy minimization—to remove any thermal contribution to the computed energies. The obtained surface energy values are found to be within the range of available experimental66–68 and simulation results.28 As shown in Fig. 4, we observe that the surface energy decreases with increasing annealing temperature. This confirms that, when subjected to thermal annealing with increasing temperatures, the silica surface becomes more stable and achieves a more energetically favorable structure. Such relaxation-induced decrease in the energy is similar to that observed upon the thermal annealing of bulk samples.69,70 Such relaxation arises from the fact that glasses are by nature out-of-equilibrium materials.71,72 The degree of relaxation of a glass can be estimated from its fictive temperature, which is defined as the temperature of the metastable equilibrium supercooled liquid exhibiting an energy comparable to that of the glass.58,73 The fictive temperature of a glass typically decreases with decreasing cooling rate as the system is able to stay at the supercooled liquid metastable equilibrium down to lower temperatures.61 Upon annealing, bulk glasses are able to partially adjust their atomic structure in order to resemble glasses formed with lower cooling rates, that is, with lower fictive temperatures. As such, here, an increase in the annealing temperature results can be effectively described as a decrease in the fictive temperature of the glass surface.

FIG. 4.

Surface energy of silica as a function of the annealing temperature. The line is guide to the eye.

FIG. 4.

Surface energy of silica as a function of the annealing temperature. The line is guide to the eye.

Close modal

After contact with water, we observe that some H atoms diffuse within the silica glass and form a partially hydrated layer with a depth of around 4 Å [see Fig. 5(a)]. This induces a structural reorganization of the silica surface, in agreement with previous simulations.22,25,49,52 In detail, we observe that all free terminating NBO, 3-fold coordinated Si, and 2-membered ring ES units that are initially present in the dry silica surfaces disappear upon the first picoseconds of hydration and convert into Si–OH silanol groups. This suggests that these structural defects are highly energetically unfavorable. In the case of the as-cut silica, i.e., with no thermal annealing, the silanol surface density is found to plateau at around 5.8 nm−2. This agrees well with the range of values (from 4 to 6 nm−2) observed experimentally16,17,74 and in previous simulations.21,52

FIG. 5.

(a) Density profile of 4-fold Si (pristine Si, denoted Siiv), 3-fold Si (defected Si, denoted Siiii), and H atoms as a function of the distance from the silica–water interface (in the absence of annealing). (b) Silanol surface density as a function of the annealing temperature. The line is guide to the eye.

FIG. 5.

(a) Density profile of 4-fold Si (pristine Si, denoted Siiv), 3-fold Si (defected Si, denoted Siiii), and H atoms as a function of the distance from the silica–water interface (in the absence of annealing). (b) Silanol surface density as a function of the annealing temperature. The line is guide to the eye.

Close modal

We now study the effect of thermal annealing on the reactivity of silica surfaces. Note that the thermal annealing phase is initially performed on the dry silica surfaces to generate a series of surfaces with different structures and degrees of stability. However, the hydration of these surfaces is subsequently performed at 300 K. Our simulations do not intend to reproduce the effect of heating on the de-hydroxylation of silica surfaces. As shown in Fig. 5(b), we observe that the surface density of silanol groups decreases with increasing annealing temperature. This demonstrates that the reactivity of the artificial silica surfaces generated herein decreases upon annealing, which results from the relaxation of high-energy surface defects. This observation is in agreement with the outcomes of experiments performed on annealed silica samples.16,17,19,75

To assess the affinity of the hydroxylated silica surfaces with water, the undissociated water molecules were manually separated from the surface by inserting a void space (20 Å) between the silica surface and the water. This spacing was chosen to be larger than the ReaxFF potential cutoff (i.e., 10 Å) so that the atoms of silica and water effectively do not “see” each other any longer. Note that the silanol groups formed during hydration were still attached to the surface. The difference of potential energy, before and after separation of the water from the silica surface, was used to calculate the silica–water interface energy (or binding energy), which offers a direct estimation of the level of hydrophilicity of the silica surfaces. As shown in Fig. 6, we note that the interface energy decreases with increasing annealing temperature. This shows that, although it is largely hydrophilic without any thermal treatment, silica becomes significantly more hydrophobic as its surface becomes more stable and exhibits fewer structural defects.

FIG. 6.

Silica–water interface (binding) energy as a function of the annealing temperature.

FIG. 6.

Silica–water interface (binding) energy as a function of the annealing temperature.

Close modal

We now investigate the relations between structural defects, silanol density, and hydrophilicity. First, as shown in Fig. 7, we observe a positive linear correlation between the silica–water interface energy and the silanol surface density. This strongly supports the idea that the hydrophilicity of silica is primarily controlled by the number of silanol groups per unit of surface. This can be understood from the fact that silanol groups can interact with polar water molecules through the formation of hydrogen bonds (H-bonds), whereas BO (siloxane groups) are intrinsically hydrophobic.19,20,26 Namely, if the silanol groups are, on average, close enough to each other, each water molecule will create two H-bonds with the silica surface (hydrophilic surface), whereas, if the silanol groups are far from each other, only one H-bond per water molecule will be created (hydrophobic surface).76 By extrapolation, we find that silica would become fully hydrophobic (i.e., negative silica–water interface energy) for a silanol surface density of around 0.6 nm−2. This is in agreement with previous simulation data,20 suggesting a hydrophilic-to-hydrophobic transition occurring between 0 and 3.7 OH nm−2.

FIG. 7.

Silica–water interface energy as a function of the surface density of silanol groups. The grey area indicates the range of surface density of silanol groups for which a hydrophobic silica surface is expected. The line is a linear fit.

FIG. 7.

Silica–water interface energy as a function of the surface density of silanol groups. The grey area indicates the range of surface density of silanol groups for which a hydrophobic silica surface is expected. The line is a linear fit.

Close modal

Finally, we investigate how the silanol surface density and, consequently, silica’s hydrophilicity are controlled by the atomic topology of the surface. In the absence of any defects, Si and O atoms create 2 and 1 BS constraints, respectively [note that each BS constraint is shared by two atoms so that the number of BS constraints created by a given atom is equal to half of its coordination number (CN)]. In addition, Si atoms create 5 BB constraints, which correspond to the 5 independent angles that need to be fixed to define the Si tetrahedra. In contrast, O atoms do not create any BB constraints due to the large angular flexibility of the Si–O–Si bond angles.37,70,77,78 This leads to an isostatic structure, i.e., nc = 3 (see enumeration in Table I). In turn, the presence of structural defects (e.g., Siiii or NBO) tends to decrease the rigidity of the surface as compared to that of the pristine bulk.

TABLE I.

Summary of the constraint enumeration for Q4 (i.e., 4-fold coordinated Si atoms comprising no terminal oxygen atom), Siiii (i.e., 3-fold coordinated Si atoms comprising 3 BO atoms), and Q3 units (i.e., 4-fold coordinated Si atoms comprising one terminal oxygen atom). For each type of unit, the table summarizes the coordination number (CN) of Si, the number of bond-stretching (BS) and bond-bending (BB) constraints per unit, the total number of constraints (BS+BB), and the number of silanol groups per Si unit. The last row indicates the variation in these quantities when a Si–O bond is broken, which results in the transformation of two Q4 units into a Siiii and a Q3 unit (see text).

UnitCNBSBBBS+BBNo. of silanol
Q4 
Siiii 
Q3 
2Q4 → Siiii + Q3  −1 −2 −3 +2 
UnitCNBSBBBS+BBNo. of silanol
Q4 
Siiii 
Q3 
2Q4 → Siiii + Q3  −1 −2 −3 +2 

To understand the linkages between surface topology and reactivity, we computed a “surface number of topological constraint per atomnc*, that is, the average number of topological constraints per atom of the surface before exposure to water. Note that the surface is defined here as a slab ranging from z = 0 to 4 Å, where z is the distance from the silica–water interface [see Fig. 5(a)]. Rather than only considering the rigidity of the “skeleton network,” that is, wherein terminating atoms are ignored, we explicitly accounted for all the atoms of the network herein (i.e., by including terminating O atoms in the constraint enumeration). Indeed, it was shown that, although the two enumeration schemes predict the same isostatic threshold (nc = 3), the explicit incorporation of 1-fold coordinated atoms into the constraints enumeration offers a more accurate description of the glass rigidity in the flexible (nc < 3) and stressed–rigid domains (nc > 3).79,80

Figure 8 shows the surface number of topological constraints per atom, nc*, as a function of the annealing temperature. Overall, the obtained nc* values are found to be in the 2.25–2.75 range. This range of values is located between two notable extrema, namely, nc = 2, which is the minimum number of constraints per atom needed to obtain a rigid 1-dimensional structure, and nc = 3, which is the minimum number of constraints needed to obtain a rigid 3-dimensional structure.34 We observe that the rigidity of the surface (i.e., as captured by nc*) increases upon annealing. This arises from the transformation of undercoordinated 3-fold Si and NBO atoms, which results in an increase in the overall connectivity of the surface.

FIG. 8.

Surface number of constraints per atom (nc*) as a function of the annealing temperature.

FIG. 8.

Surface number of constraints per atom (nc*) as a function of the annealing temperature.

Close modal

Figure 9(a) shows the number of silanol groups per Si atom (limited to the Si atoms belonging to the surface) as a function of the surface number of constraints per atom nc*. We observe that the number of silanol groups per Si atom decreases fairly linearly with nc* and eventually becomes zero when nc*=3. This can be understood as follows. The nc*=3 threshold corresponds to the value obtained for pristine silica, that is, when the atomic network is fully polymerized and only comprise hydrophobic siloxane groups. In this situation, silica does not exhibit any “open” chemisorption site for water molecules to dissociate into hydroxyl groups. Starting from this ideal network, the breakage of a Si–O bond can be expressed as the following reaction 2Q4 → Siiii + Q3, where a Qn unit is defined as a 4-fold coordinated Si atom connected to n other Si units (i.e., comprising 4 − n terminating O atoms) and Siiii being a 3-fold coordinated Si (E′ center) connected to 3 other Si units (i.e., with 3 BO atoms). As summarized in Table I, this reaction induces the loss of 3 constraints (1 BS and 2 BB). In turn, the resulting Siiii and Q3 defects offer two potential additional chemisorption sites for silanol groups. After normalization by the number of atoms in the surface, one gets the following expression for the number of potential silanol groups per Si atom in the surface OH/Si:

OH/Si=2(3nc*).
(1)

As shown in Fig. 9(a), the computed number of silanol groups (after hydration) agrees very well with the prediction based on the surface number of constraints per atom offered by Eq. (1) (i.e., relying on the atomic topology of the surface before hydration). This suggests that the short-term reactivity of silica (as measured by the surface density of silanol groups) can be fully predicted from the atomic topology of the unreacted surface (with no fitting parameter).

FIG. 9.

(a) Number of silanol groups per surface Si atom as a function of the surface number of constraints per atom nc*. The line is a theoretical prediction given by Eq. (1) (no fitting parameters, see text). (b) Silica–water interface energy as a function of the surface number of constraints per atom nc*. The line is a linear fit. The grey area indicates the range of surface number of constraints per atom for which a hydrophobic silica surface is expected.

FIG. 9.

(a) Number of silanol groups per surface Si atom as a function of the surface number of constraints per atom nc*. The line is a theoretical prediction given by Eq. (1) (no fitting parameters, see text). (b) Silica–water interface energy as a function of the surface number of constraints per atom nc*. The line is a linear fit. The grey area indicates the range of surface number of constraints per atom for which a hydrophobic silica surface is expected.

Close modal

Since the silanol surface density offers a good estimation of silica’s affinity for water (see Fig. 7), the correlation between silanol surface density and hydrophilicity suggests that the atomic topology of the surface controls the hydrophilicity of silica. The close relationship between surface hydrophilicity and surface topology is demonstrated in Fig. 9(b), which shows the water–silica interface (binding) energy as a function of the number of surface topological constraints per atom, nc*. We observe that the affinity of the silica surface for water (i.e., as captured by the water–silica interface energy) decreases linearly with nc*. This suggests that the more rigid the surface is (i.e., more connected, with higher nc* values), the more hydrophobic it becomes. Namely, flexible (low nc* values) and more rigid (nc* values close to 3) surfaces are associated with hydrophilic and hydrophobic behaviors, respectively. Based on our results, we identify that the hydrophilic-to-hydrophobic transition occurs around nc*=2.82.

Reactive molecular dynamics simulations were used to artificially generate a series of dry (de-hydroxylated) silica surfaces exhibiting varying atomic structures. It was observed that, upon thermal annealing, the surface density of structural defects decreases, thereby increasing the overall thermodynamic stability of the surface. The surfaces were subsequently exposed to water to investigate the relationship between the structure and surface reactivity. It was observed that the more stable the surface is, the lower the aqueous reactivity is. In particular, the surface density of silanol groups was found to be lower on pre-annealed surfaces—this effect increases upon increasing annealing temperatures. In turn, it was shown that the hydrophilicity of the surface is directly controlled by the surface density of silanol groups. Based on topological constraint theory, it was demonstrated that both surface reactivity and hydrophilicity are determined by the atomic topology of the surface.

Altogether, these results suggest that tuning the atomic topology of the surface of silica is a promising route toward the design of novel phases with tailored reactivity and hydrophilicity. Besides thermal treatments,36 the atomic topology can be altered by compositional design37,77 (e.g., by inserting various fractions of network-modifying and network-forming species), pressure effects,81,82 exposure to irradiation,43,44,59,83–85 or ion exchange treatments.86,87 Hence, various routes are available—individually or in combination—to finely control the atomic architecture of the surface of disordered networks in order to optimize their reactivity and hydrophilicity.

More generally, this study extends for the first time the topological constraint theory to investigate the rigidity of surfaces—as opposed to bulk atomic networks. It is notable that, although silicates interact with their aqueous environment through their surface, their dissolution rate was shown to be controlled by the topology of their bulk atomic structure,7 whereas both their surface reactivity and hydrophilicity appear here to be controlled by the topology of their atomic structure. Clearly, further studies are needed to better understand the linkages, if any, between the topology of the bulk and that of the surface, and to discern which properties depend on the topology of the surface from those that are controlled by the topology of the bulk.

See supplementary material for the list of ReaxFF parameters used for the present molecular dynamics simulations.

This work was supported by the National Science Foundation under Grant No. 1562066.

1.
C. M.
Jantzen
,
K. G.
Brown
, and
J. B.
Pickett
,
Int. J. Appl. Glass Sci.
1
,
38
(
2010
).
3.
M. M.
Smedskjaer
,
R. E.
Youngman
, and
J. C.
Mauro
,
Appl. Phys. A
116
,
491
(
2014
).
5.
T.
Oey
,
J.
Timmons
,
P.
Stutzman
,
J. W.
Bullard
,
M.
Balonis
,
M.
Bauchy
, and
G.
Sant
,
J. Am. Ceram. Soc.
100
,
4785
(
2017
).
6.
T.
Oey
,
A.
Kumar
,
I.
Pignatelli
,
Y.
Yu
,
N.
Neithalath
,
J. W.
Bullard
,
M.
Bauchy
, and
G.
Sant
,
J. Am. Ceram. Soc.
100
,
5521
(
2017
).
7.
I.
Pignatelli
,
A.
Kumar
,
M.
Bauchy
, and
G.
Sant
,
Langmuir
32
,
4434
(
2016
).
8.
J. C.
Mauro
and
E. D.
Zanotto
,
Int. J. Appl. Glass Sci.
5
,
313
(
2014
).
9.
A.
Burneau
,
J.
Gallas
, and
A.
Legrand
,
Vibrational Spectroscopy
(
John Wiley
,
Chichester
,
1998
), p.
147
.
10.
A.
Rimola
,
D.
Costa
,
M.
Sodupe
,
J.-F.
Lambert
, and
P.
Ugliengo
,
Chem. Rev.
113
,
4216
(
2013
).
11.
M. L.
Hair
,
J. Non-Cryst. Solids
19
,
299
(
1975
).
12.
S.
Haastrup
,
M. S.
Bødker
,
S. R.
Hansen
,
D.
Yu
, and
Y.
Yue
,
J. Non-Cryst. Solids
481
,
556
(
2018
).
13.
B.
Akhavan
,
K.
Jarvis
, and
P.
Majewski
,
Powder Technol.
249
,
403
(
2013
).
14.
N.
Giovambattista
,
P. G.
Debenedetti
, and
P. J.
Rossky
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
15181
(
2009
).
15.
J.
Laskowski
and
J. A.
Kitchener
,
J. Colloid Interface Sci.
29
,
670
(
1969
).
16.
L. T.
Zhuravlev
,
Colloids Surf., A
173
,
1
(
2000
).
17.
L. T.
Zhuravlev
,
Langmuir
3
,
316
(
1987
).
18.
A. S.
D’Souza
and
C. G.
Pantano
,
J. Am. Ceram. Soc.
82
,
1289
(
1999
).
19.
G.
Ori
,
C.
Massobrio
,
A.
Bouzid
, and
B.
Coasne
,
Molecular Dynamics Simulations of Disordered Materials
(
Springer
,
Cham
,
2015
), pp.
345
365
.
20.
B.
Siboulet
,
B.
Coasne
,
J.-F.
Dufrêche
, and
P.
Turq
,
J. Phys. Chem. B
115
,
7881
(
2011
).
21.
J.
Du
and
A. N.
Cormack
,
J. Am. Ceram. Soc.
88
,
2532
(
2005
).
22.
S. H.
Garofalini
,
J. Non-Cryst. Solids
120
,
1
(
1990
).
23.
E. A.
Leed
and
C. G.
Pantano
,
J. Non-Cryst. Solids
325
,
48
(
2003
).
24.
B. P.
Feuston
and
S. H.
Garofalini
,
J. Chem. Phys.
91
,
564
(
1989
).
25.
B. P.
Feuston
and
S. H.
Garofalini
,
J. Appl. Phys.
68
,
4830
(
1990
).
26.
V. A.
Bakaev
and
W. A.
Steele
,
J. Chem. Phys.
111
,
9803
(
1999
).
27.
J. C.
Fogarty
,
H. M.
Aktulga
,
A. Y.
Grama
,
A. C. T.
van Duin
, and
S. A.
Pandit
, “
A reactive molecular dynamics simulation of the silica-water interface
,”
J. Chem. Phys.
132
,
174704
(
2010
).
28.
J. M.
Rimsza
,
R. E.
Jones
, and
L. J.
Criscenti
,
Langmuir
33
,
3882
(
2017
).
29.
F. S.
Emami
,
V.
Puddu
,
R. J.
Berry
,
V.
Varshney
,
S. V.
Patwardhan
,
C. C.
Perry
, and
H.
Heinz
,
Chem. Mater.
26
,
2647
(
2014
).
30.
31.
J. C.
Mauro
,
C. S.
Philip
,
D. J.
Vaughn
, and
M. S.
Pambianchi
,
Int. J. Appl. Glass Sci.
5
,
2
(
2014
).
32.
J. C.
Phillips
,
J. Non-Cryst. Solids
34
,
153
(
1979
).
33.
J. C.
Phillips
,
J. Non-Cryst. Solids
43
,
37
(
1981
).
34.
J. C.
Mauro
,
Am. Ceram. Soc. Bull.
90
,
31
(
2011
).
35.
M.
Bauchy
,
Am. Ceram. Soc. Bull.
91
,
34
(
2012
).
36.
M.
Bauchy
and
M.
Micoulaut
,
J. Non-Cryst. Solids
357
,
2530
(
2011
).
37.
M.
Bauchy
,
M.
Micoulaut
,
M.
Celino
,
S.
Le Roux
,
M.
Boero
, and
C.
Massobrio
,
Phys. Rev. B
84
,
054201
(
2011
).
38.
J. C.
Maxwell
,
Philos. Mag. Ser. 4
27
,
294
(
1864
).
39.
N.
Mascaraque
,
M.
Bauchy
, and
M. M.
Smedskjaer
,
J. Phys. Chem. B
121
,
1139
(
2017
).
40.
N.
Mascaraque
,
M.
Bauchy
,
J. L. G.
Fierro
,
S. J.
Rzoska
,
M.
Bockowski
, and
M. M.
Smedskjaer
,
J. Phys. Chem. B
121
,
9063
(
2017
).
41.
I.
Pignatelli
,
A.
Kumar
,
R.
Alizadeh
,
Y. L.
Pape
,
M.
Bauchy
, and
G.
Sant
,
J. Chem. Phys.
145
,
054701
(
2016
).
43.
I.
Pignatelli
,
A.
Kumar
,
K. G.
Field
,
B.
Wang
,
Y.
Yu
,
Y.
Le Pape
,
M.
Bauchy
, and
G.
Sant
,
Sci. Rep.
6
,
20155
(
2016
).
44.
Y.-H.
Hsiao
,
E. C.
La Plante
,
N. M. A.
Krishnan
,
Y.
Le Pape
,
N.
Neithalath
,
M.
Bauchy
, and
G.
Sant
,
J. Phys. Chem. A
121
,
7835
(
2017
).
45.
A. C. T.
van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
,
J. Phys. Chem. A
105
,
9396
(
2001
).
46.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
47.
48.
A. K.
Rappe
and
W. A.
Goddard
 III
,
J. Phys. Chem.
95
,
3358
(
1991
).
49.
J. M.
Rimsza
,
J.
Yeon
,
A. C. T.
van Duin
, and
J.
Du
,
J. Phys. Chem. C
120
,
24803
(
2016
).
50.
Y.
Yu
,
B.
Wang
,
M.
Wang
,
G.
Sant
, and
M.
Bauchy
,
J. Non-Cryst. Solids
443
,
148
(
2016
).
51.
M. C.
Pitman
and
A. C. T.
van Duin
,
J. Am. Chem. Soc.
134
,
3042
(
2012
).
52.
J. C.
Fogarty
,
H. M.
Aktulga
,
A. Y.
Grama
,
A. C. T.
van Duin
, and
S. A.
Pandit
,
J. Chem. Phys.
132
,
174704
(
2010
).
53.
J.
Yeon
and
A. C. T.
van Duin
,
J. Phys. Chem. C
120
,
305
(
2016
).
54.
J. F.
Stebbins
and
P.
McMillan
,
J. Non-Cryst. Solids
160
,
116
(
1993
).
55.
P.
Heyliger
,
H.
Ledbetter
, and
S.
Kim
,
J. Acoust. Soc. Am.
114
,
644
(
2003
).
56.
Y.
Yu
,
B.
Wang
,
M.
Wang
,
G.
Sant
, and
M.
Bauchy
,
Int. J. Appl. Glass Sci.
8
,
276
(
2017
).
57.
M.
Bauchy
,
J. Chem. Phys.
141
,
024507
(
2014
).
58.
N. M. A.
Krishnan
,
B.
Wang
,
Y.
Yu
,
Y.
Le Pape
,
G.
Sant
, and
M.
Bauchy
,
Phys. Rev. X
7
,
031019
(
2017
).
59.
N. M. A.
Krishnan
,
B.
Wang
,
Y.
Le Pape
,
G.
Sant
, and
M.
Bauchy
,
J. Chem. Phys.
146
,
204502
(
2017
).
60.
N. P.
Bansal
and
R. H.
Doremus
,
Handbook of Glass Properties
(
Elsevier
,
2013
).
61.
X.
Li
,
W.
Song
,
K.
Yang
,
N. M. A.
Krishnan
,
B.
Wang
,
M. M.
Smedskjaer
,
J. C.
Mauro
,
G.
Sant
,
M.
Balonis
, and
M.
Bauchy
,
J. Chem. Phys.
147
,
074501
(
2017
).
62.
K.
Vollmayr
,
W.
Kob
, and
K.
Binder
,
Phys. Rev. B
54
,
15808
(
1996
).
63.
J. M. D.
Lane
,
Phys. Rev. E
92
,
012320
(
2015
).
64.
D. I.
Grimley
,
A. C.
Wright
, and
R. N.
Sinclair
,
J. Non-Cryst. Solids
119
,
49
(
1990
).
65.
S. M.
Levine
and
S. H.
Garofalini
,
J. Chem. Phys.
86
,
2997
(
1987
).
66.
S. M.
Wiederhorn
,
J. Am. Ceram. Soc.
52
,
99
(
1969
).
67.
Y. K.
Shchipalov
,
Glass Ceram.
57
,
374
(
2000
).
68.
J.
Mizele
,
J. L.
Dandurand
, and
J.
Schott
,
Surf. Sci.
162
,
830
(
1985
).
69.
M. M.
Smedskjaer
,
M.
Bauchy
,
J. C.
Mauro
,
S. J.
Rzoska
, and
M.
Bockowski
,
J. Chem. Phys.
143
,
164505
(
2015
).
70.
M.
Wang
,
B.
Wang
,
T. K.
Bechgaard
,
J. C.
Mauro
,
S. J.
Rzoska
,
M.
Bockowski
,
M. M.
Smedskjaer
, and
M.
Bauchy
,
J. Non-Cryst. Solids
454
,
46
(
2016
).
71.
Y.
Yu
,
M.
Wang
,
D.
Zhang
,
B.
Wang
,
G.
Sant
, and
M.
Bauchy
,
Phys. Rev. Lett.
115
,
165901
(
2015
).
72.
Y.
Yu
,
M.
Wang
,
M. M.
Smedskjaer
,
J. C.
Mauro
,
G.
Sant
, and
M.
Bauchy
,
Phys. Rev. Lett.
119
,
095501
(
2017
).
73.
A. K.
Varshneya
,
Fundamentals of Inorganic Glasses
(
Academic Press, Inc.
,
1993
).
74.
R.
Mueller
,
H. K.
Kammler
,
K.
Wegner
, and
S. E.
Pratsinis
,
Langmuir
19
,
160
(
2003
).
75.
S. J.
Gregg
and
K. S.
Sing
, “
Adsorption, Surface Area and Porosity
” (
Auflage, Academic Press
,
London
,
1982
), http://onlinelibrary.wiley.com/doi/10.1002/bbpc.19820861019/abstract.
76.
V.
Bolis
,
B.
Fubini
,
L.
Marchese
,
G.
Martra
, and
D.
Costa
,
J. Chem. Soc., Faraday Trans.
87
,
497
(
1991
).
77.
M.
Bauchy
,
M. J. A.
Qomi
,
C.
Bichara
,
F.-J.
Ulm
, and
R. J.-M.
Pellenq
,
Phys. Rev. Lett.
114
,
125502
(
2015
).
78.
M.
Zhang
and
P.
Boolchand
,
Science
266
,
1355
(
1994
).
79.
P.
Boolchand
and
M. F.
Thorpe
,
Phys. Rev. B
50
,
10366
(
1994
).
80.
P.
Boolchand
,
M.
Zhang
, and
B.
Goodman
,
Phys. Rev. B
53
,
11488
(
1996
).
81.
M.
Bauchy
and
M.
Micoulaut
,
Phys. Rev. Lett.
110
,
095501
(
2013
).
82.
M.
Bauchy
and
M.
Micoulaut
,
Nat. Commun.
6
,
6398
(
2015
).
83.
B.
Wang
,
N. M. A.
Krishnan
,
Y.
Yu
,
M.
Wang
,
Y.
Le Pape
,
G.
Sant
, and
M.
Bauchy
,
J. Non-Cryst. Solids
463
,
25
(
2017
).
84.
N. M. A.
Krishnan
,
B.
Wang
,
Y.
Le Pape
,
G.
Sant
, and
M.
Bauchy
, “
Irradiation-driven amorphous-to-glassy transition in quartz: The crucial role of the medium-range order in crystallization
,”
Phys. Rev. Mater.
1
,
053405
(
2017
).
85.
N. M. A.
Krishnan
,
B.
Wang
,
G.
Sant
,
J. C.
Phillips
, and
M.
Bauchy
,
ACS Appl. Mater. Interfaces
9
,
32377
(
2017
).
86.
M.
Wang
,
B.
Wang
,
N. M. A.
Krishnan
,
Y.
Yu
,
M. M.
Smedskjaer
,
J. C.
Mauro
,
G.
Sant
, and
M.
Bauchy
,
J. Non-Cryst. Solids
455
,
70
(
2017
).
87.
M.
Wang
,
M. M.
Smedskjaer
,
J. C.
Mauro
,
G.
Sant
, and
M.
Bauchy
,
Phys. Rev. Appl.
8
,
054040
(
2017
).

Supplementary Material