Ground-up generation of periodic slab models of dehydroxylated amorphous silica of varying roughness

Abstract

Simulation and systematic analysis of the surfaces of amorphous materials is a challenge for computational chemistry. For example, silica has found widespread industrial use as an adsorbent and catalyst support but available models for use with periodic DFT are limited in variety and representativeness of realistic materials. Herein we present a generic approach for the systematic construction of ensembles of amorphous materials surface models with varied roughness and termination characteristics. The power of the approach is shown with silica as the representative example. By combining MD simulations and Fourier-series-based randomization, bulk amorphous silica was modeled and cleaved to produce surfaces with systematically varied roughness and surface saturation. An automated saturation procedure resulted in surface models with silanol densities typical of high-temperature activation protocols in the range 0.35–2.00 OH nm⁻², in excellent agreement with the experimental data on surface chemistry of dehydroxylated silica materials.

---

1 Introduction

The knowledge of the molecular structure of materials is critical for understanding and predicting their properties through computational modeling. For crystalline materials, their regular, periodic structures simplify model construction.¹ However, many widely used materials lack this long-range periodic order,² such as amorphous silica (a-SiO₂), which finds diverse applications, including pharmaceuticals,^3,4 optics,⁵ and catalysis.^6–8

High-temperature (ca. 700 °C) dehydroxylation in vacuo is a common approach to reduce the heterogeneity of reactive sites on silica surfaces and produce model support materials well-suited for surface organometallic chemistry applications^9–11 and the generation of single-site catalysts.^6,12 For a-SiO₂, the resulting support can be referred to as SiO₂₋₇₀₀. This treatment results in surfaces with silanol group (Fig. 1) concentrations of ca. 1.15 OH nm⁻² (ref. 13 and 14) with a still debatable distribution of types.¹ The results of Zhuravlev¹³ suggest that after high-temperature dehydroxylation, silica surface features exclusively isolated (90%) and geminal (10%) silanols, while the study by Fleischman and Scott¹⁵ presented experimental evidence of the formation of also vicinal SiOH species (Fig. 1). Beyond the types of silanol groups present, the precise local structures and catalytic roles of the diverse sites on amorphous surfaces are also often unknown.^16–18

Fig. 1 Different types of silanol groups found on the surface of amorphous silica.

Computational modeling offers an opportunity to shed light on the structure and reactivity of such surface species. Modern approaches for constructing atomistic models of amorphous materials either focus on representing the local environment of the specific surface sites with finite cluster models^4,19–23 or approximate the amorphous surface with distorted periodic slab models.^14,24–26 While cluster models allow for detailed and highly accurate computations of local chemical behavior, they generally lack the high strain expected in highly dehydroxylated silica and a means of accounting for the inherent heterogeneity of the surface. Periodic slab models do not fall victim to these issues and provide a basis for systematically exploring the varied chemical environments of amorphous materials.

Fig. 2 schematically illustrates the common procedure for the construction of periodic slab surface models of amorphous materials (e.g. SiO₂₋₇₀₀). First, the crystal structure of SiO₂ (usually β-cristobalite) undergoes distortion via force-field MD simulation of melting and cooling processes. Next, a surface is created by cleaving the distorted structure with a plane at a fixed z-coordinate. The resulting dangling bonds are saturated by silanol groups, followed by successive condensation of neighboring SiOH pairs to ensure the overall charge neutrality and yield the dehydroxylated surface model.

Fig. 2 Overview of the general procedure commonly used for creating dehydroxylated periodic slab models of amorphous silica.

This protocol results in relatively flat surface models^14,24,25,27 and it does not allow for systematic variation and exploration of surface roughness, which plays a decisive role in surface properties^28,29 and reactivity^30,31 of a-SiO₂. The roughness of the resulting model can be increased through quicker rates of cooling, although to a relatively limited extent.³² Nguyen and Laird³³ proposed using a stochastic Fourier series in place of a flat plane at the cleaving step and showed that through their method, the roughness could effectively be controlled.

The surface resulting from the saturation step generally has silanol concentrations between 6 and 8 OH nm⁻².^14,24,25,27 For the reduction of silanol concentration to the dehydrated surface, the condensation step is ideally accompanied by relaxation and evaluation of the surface energetics. The consequences of a given condensation on surface stabilities cannot be known a priori,¹⁴ leading to noticeable computational and manual effort. Alternative approaches based on the addition of SiOH until the desired density is achieved have been proposed, efficient enough to produce a large number of models.^32,33 However, the implementations may result in charged surfaces³³ or remaining dangling bonds.^32,33

Herein we present a method for the systematic generation of surfaces of amorphous materials enabling systematic variation of surface roughness and algorithmic surface saturation. The method is demonstrated for a-SiO₂, for which a comprehensive set of models was generated. The resulting silanol concentrations and compositions were evaluated and compared to experiment and previously created models. To further validate models, atomistic thermodynamics was used to compute stability diagrams over temperature and P_H2O for three select surfaces functionalized with additional silanol groups.

2. Methods

A model of bulk a-SiO₂ was generated using the classical molecular dynamics carried out using the Large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)³⁴ software package. To describe the interaction between atoms the forcefield of van Beest, Kramer, and van Santen (BKS)³⁵ was used overlayed with the Lennard-Jones potential described by Wimalasiri et al.³² (parameters found in ESI† Table S1). An Ewald summation³⁶ to a tolerance of 10⁻⁴ was used to account for long-range coulomb interactions. Using a time step of 0.5 fs, a 3 × 3 × 2 supercell of β-cristobalite³⁷ (144 units of SiO₂) was heated to 8000 K at constant pressure (1 atm) and subsequently cooled to 298 K at a rate of 1 K ps⁻¹. This resulting bulk silica sample was heated to 4000 K and allowed to equilibrate for 200 ps then cooled to 298 K at a rate of 1 K ps⁻¹. During these simulations, the x and y dimensions were held constant at 21.5 Å to simplify analysis. To control temperature and pressure, a Berendsen thermostat (damping constant of 0.1 ps) and barostat³⁸ (damping constant of 0.1 ps, modulus of 360 000 atm) were employed. The resulting amorphized bulk was cleaved using the method of Nguyen and Laird³³ using 10 different stochastic Fourier series with roughness being controlled through the empirical parameter α, which is essentially an attenuation factor controlling the decay of Fourier coefficients based on frequency magnitude. In simple terms, lower values of α result in surfaces of higher roughness, on average. For more details, the reader is referred to the work of Nguyen and Laird.³³ For this study, α 's equal to 0.01, 0.02, 0.03, 0.04, 0.05, and 0.1 were used, resulting in a total of 60 different structures. To cleave the surface, atoms were not removed but rather shifted downwards by the z length of the simulation box, schematically shown in Fig. 3 (blue). Each structure was subsequently re-equilibrated for 2 ns at 298 K under an NVT ensemble with a 15 Å vacuum in the z-dimension.

Fig. 3 General overview of the method implemented to create surfaces and saturate dangling bonds of the models.

The saturation of models was achieved through a genetic algorithm, summarized in Fig. 3. The initial, dry, structure was relaxed and dangling bonds were saturated in a step-wise manner. At each step, a single H₂O molecule was added to the structure which was then relaxed. 10 possible structures were considered per step and partially optimized, after which the one of lowest electronic energy was carried forward. This continued until there were no remaining dangling bonds as they are uncommon at the high-temperature activation conditions (low humidity, 700 °C).³⁹ Si and O atoms were considered bound to each other if they were within 2 Å of each other. This assumption is based on a geometric constraint that atoms are in contact if they are closer than 62% of the sum of their van der Waals radii.⁴⁰ If necessary, dangling atoms were induced by shifting a Si atom within 2.0 Å of a 3-coordinate oxygen defect (³O). To achieve higher silanol densities, the models were functionalized through probing for surface atoms, as described in ref. 32, and water molecules were added across Si–O bonds, which were selected at random. Probing was carried out through splitting the x–y plane into a square grid with a spacing of 0.15 Å between points and a probe of 1.4 Å descending from the vacuum slab until contact is made with an atom. Contact was determined when the radius of the probe intersects the contact radius of an atom, defined as 1.04 Å for Si atoms and 0.76 Å for O atoms.

The silanol density of the given sample can be assessed experimentally through physisorption of liquid Kr or N₂.^13,14 Taking inspiration from this, the solvent-accessible surface area for Kr was calculated as to estimate the surface area of generated models. To get this estimate, a polygon mesh of the structure was generated, like that shown in Fig. 4, using the marching cubes algorithm⁴¹ as implemented in scikit-image.⁴² There are two independent surfaces in the slab model. Both surfaces could be a valid model, so both sides were analyzed as individual data points. The surface area was taken as the sum of the area of all tetrahedrons that are part of the mesh. Starting from the top-most and bottom-most vertices of the mesh, the polygons belonging to either surface of the model were identified through a breadth-first search for all connecting vertices. The roughness of the model was calculated as the mean-square deviation of these vertices. To identify which silanol groups belong to which side of the model, k-means clustering⁴³ based on the z-coordinate of the H atoms was used. It was assumed that there would be two distinct clusters corresponding to each side of the model and that silanol groups would find themselves close to the surfaces of the model.

Fig. 4 Representation of the meshes generated used to estimate surface areas of a model separated into the (a) “top” of the model and (b) “bottom” of the model.

For the entirety of the algorithm outlined in Fig. 3, structure relaxations were performed using the semi-empirical SCC-DFTB⁴⁴ method as implemented in CP2K.⁴⁵ The parameter set of Guimarães et al.⁴⁶ was benchmarked against PBE-D3 results of Comas-Vives²⁴ (ESI† Fig. S1). A plane-wave energy cutoff of 450 Ry and a relative cutoff of 60 Ry were selected. Dispersion effects were accounted for using the D3BJ dispersion correction scheme of Grimme et al.⁴⁷ To increase throughput, only during algorithmic saturation, the structures were partially optimized over 25 steps, resulting in forces generally below 1 × 10⁻³ Hartree Bohr⁻¹. Otherwise, structures were optimized to tolerances of 4.5 × 10⁻⁴ and 3 × 10⁻³ Bohr for force and displacement, respectively.

Stability diagrams were constructed using atomistic thermodynamics.⁴⁸ The Gibbs free energy (ΔG) was estimated using the equation

ΔG = E_SiO2,ρhigh − E_SiO2,surf − n_H2Oμ_H2O(T,P_H2O) (1)

where E_SiO2,ρhigh is the electronic energy of the surface of the highest silanol density, E_SiO2,surf is the electronic energy of the silica surface for which ΔG is to be calculated, and n_H2O is the difference in the number of water molecules between the surfaces. μ_H2O(T,P_H2O) is the chemical potential of water which can be approximated by the ideal gas law through which the energy of the system at varying temperature (T) and partial pressure of water (P_H2O) can be determined:

μ_H2O(T,P_H2O) = E_H2O + μ_H2O(T,P°) + k_BT ln(P_H2O/P°) (2)

where E_H2O is the electronic energy of water, and μ_H2O(T,P°) is the chemical potential of water at a given temperature T and standard pressure P°. This term contains all contributions from rotations, vibrations, and ideal gas entropy at 1 atm. While this term can be calculated from first principles, the values stated in the JANAF thermodynamic tables⁴⁹ were used instead. This method assumes that the phonon modes do not change during the functionalization process; even so, the error in this assumption is known to be smaller than the error inherent to DFT⁵⁰ for a surface made with a perfect β-cristobalite. For evaluation of this assumption in the context of amorphous surfaces, see Section S2 in ESI.†

3. Results and discussion

In preparation for algorithmic saturation, 60 slabs of different roughness had the coordination of each of their atoms calculated after the step “Relax initial structure” in Fig. 3. The saturation occurs at dangling atoms, namely singly bound O atoms (¹O) and 3-coordinated Si atoms (³Si). Fig. 5 summarizes the distributions of such atoms among all slabs. One can see that the distributions of ³Si and ¹O do not overlap substantially. Naively saturating all dangling bonds would therefore lead to models of net formal charge, non-physical considering the system we are aiming to model is a solid surface in contact with a gaseous phase. Thus, 3-coordinate O atoms (³O) were used to create additional ³Si to give formally-neutral slabs. The distributions of ¹O and ³O overlap significantly, resulting in an excess of ³O that can be used. As the energetic effect of removing silanol groups cannot be known a priori, the impact of adding groups using one ³O over another was screened before choosing which structure to carry forward, in hopes of avoiding unnecessary additional instability to the model.

Fig. 5 Count of a given number of defects occurring within a model for (blue) under-coordinated Si atoms (³Si), (red) under-coordinated O atoms (¹O), and (green) over-coordinated O atoms (³O).

The relative stabilities of the saturated slabs were analyzed. Energies were normalized per SiO₂ unit relative to that of bulk β-cristobalite and H₂O. For comparison, the readily available models of Comas-Vives²⁴ and Rozanska et al.⁵⁰ were relaxed using the same semi-empirical SCC-DFTB parameterization used during algorithmic saturation and the same relative stability was calculated, the results of which are summarized in Fig. 6. As expected, surface structures are less stable than bulk β-cristobalite, however, this is not true for all the previous models. All models made by Rozanska et al.⁵⁰ and the model representing 2.0 OH nm⁻² of Comas-Vives²⁴ are more stable compared to the chosen reference states. We ascribe this over-stabilization to the asymmetric hydroxylation of the slabs. One side of those models remains at its initial silanol density and is not dehydroxylated in tandem with the other; thus, a network of hydrogen bonds remains. Our models are 1.5–27 kJ mol⁻¹ (16 kJ mol⁻¹ per SiO₂ unit on average) less stable than the cristobalite reference and have two highly hydroxylated sides instead of one. These stability trends are in line with the expected lower intrinsic stability of the amorphous structures than that of the crystalline ones.

Fig. 6 Count of models of given relative energy per unit SiO₂ (reference states of β-cristobalite and H₂O) of this work, Rozanska et al.,⁵⁰ and Comas-Vives.²⁴

The algorithmic saturation gives rise to versatile silanol environments on the surface, which together with the variation in the roughness establish a wide range of chemical and topological compositions. Fig. 7a presents the probability of a given number of geminal groups forming on a surface using the presented approach. About 60% of surfaces had no geminal pairs, about 30% had one pair of geminal –OH groups and the remaining had two pairs. A surface with an odd number of geminal groups was found and is an artifact due to insufficient separation between the two surfaces of the slab. For simplicity, these models were not included in further analysis. Experimentally, it is postulated that 10% of silanol groups present on SiO₂₋₇₀₀ are geminal.¹³ For the generated ensemble of surfaces, geminal silanol groups make up 14% of all SiOH.

Fig. 7 Count of surfaces within the data set with a given number of (a) –OH groups part of a geminal pair (b) vicinal silanol groups.

Fig. 7b reports the number of surfaces with given numbers of vicinal silanol pairs. Surfaces of 0, 2, or 4 vicinal groups were predominantly generated with 2 groups accounting for 35% of the dataset and those of 0 or 4 accounting for 20%. Other periodic models for SiO₂₋₇₀₀ have either one or no vicinal groups present. This is in line with the experimental findings that vicinal groups may be present even after SiO₂ dehydroxylation at 800 °C.¹³ Furthermore, models with vicinal groups have been able to accurately reproduce experiments such as IR spectra^24,25 or experimental silanol density as a function of temperature.¹⁴

To account for the generated models having increased roughness, the area of each surface, shown in Fig. 8a, was estimated to determine the silanol density more accurately. This resulted in values of surface area between 5.8 nm² and 11.6 nm², an average of 7.8 nm². Using the previous methodology, the flat cleaving plane allows for the surface area to be estimated through the xy cross-sectional area of the periodic unit cell. This results in a surface area of 4.6 nm² for all generated surfaces. In the worst case, this is a difference of a factor of 3. Compared to the average, this is underestimated by a factor of 1.7. Applying this method to the readily available models in literature, similar results were achieved (ESI† Table S2). In all cases, the calculated surface area of the dehydroxylated side of the model is approximately 50% larger than the cross-sectional area of the simulation box. For the models of Comas-Vives,²⁴ this could be a fair approximation, for those of Rozanska et al.⁵⁰ it could be as well. Visually inspecting the latter models, there are slight cavities that could lead to this greater-than-expected surface area. In any case, for models of amorphous nature, these estimates should be a sufficiently close approximation to the surface area of a model.

Fig. 8 (a) Count of surfaces with a given solvent-accessible surface area. (b) Count of surfaces with given silanol densities. (c) Relationship between surface roughness and surface area with the densities of the points given through the color of the scatter.

Based on the calculated surface areas, silanol densities ranging from 0.35 OH nm⁻² to 2.00 OH nm⁻² were generated, with an average density of 0.96 OH nm⁻² (Fig. 8b). The range of silanol densities represents surfaces treated between 1000 °C and approximately 400 °C, respectively. The average almost perfectly matches the expected silanol density of SiO₂₋₇₀₀.^13,14 In line with prior studies, the roughness of our surfaces correlates strongly with their surface area (Fig. 8c).³² That said, there is no simple relationship between the resulting silanol density or roughness ultimately implying that the proposed method samples a wide range of topological and unique silanol group environments.

To establish whether the generated models could be representative of reaction conditions, atomistic thermodynamics was applied to calculate the relative stabilities. Surfaces of low, medium, and high roughness (surfaces noted in Table S3, ESI†) were further functionalized through the simple addition of H₂O across siloxane bonds to get a range of silanol densities. Fig. 9 reports the silanol density of the most stable surface resulting from this additional functionalization as a function of T and P_H2O. For the high-temperature dehydroxylation conditions employed for the generation of SiO₂₋₇₀₀, the most stable medium- and high-roughness surfaces show SiOH densities of 1.1 OH nm⁻² or 1.0 OH nm⁻² in a perfect agreement with the experiment.¹³ Following the range of surface silanol densities generated, the stable surfaces at 400 °C, 500 °C, and 600 °C under vacuum are also denoted. Experimentally, they are determined to be 2.35 OH nm⁻², 1.8 OH nm⁻², and 1.5 OH nm⁻², respectively.¹³ At 400, all three model surfaces are in line with these results, showing stable densities of 2.2 OH nm⁻², 2.0 OH nm⁻², and 1.9 OH nm⁻², for low, medium, and high roughness. Similar results are also seen for 500 °C. At 600 °C, the stable surface is the same as that stable at 700 °C.

Fig. 9 Stability diagrams as a function of log(P_H2O/P°) and temperature (°C) constructed through further functionalizing a surface of (a) low, (b) medium, and (c) high surface roughness. Scatter points are colored according to experimental silanol density at thermodynamic conditions for 1.15 OH nm⁻² (700 °C), 1.5 OH nm⁻² (600 °C), 1.8 OH nm⁻² (500 °C), 2.35 OH nm⁻² (400 °C). Colors are chosen according to the color bar. Scatter points are labeled with the calculated stable silanol density for the surface at given thermodynamic conditions.

As the number of Si–O bonds that can be randomly broken for any one addition H₂O is so large within the ensemble, resulting in errors in stability diagrams due to non-exhaustive search of configurational space and, thus, neglect of configurational entropy may be significant. This should be the most noticeable at higher silanol densities and still difficult to quantify even at low silanol densities. The process of further functionalizing the surfaces was repeated to assess whether results can be replicated. This was done for the high and low roughness surfaces (ESI† Fig. S3 and S4). As silanol density increases, deviations in the stable densities at given T and P_H2O do indeed arise. For the low roughness surface, the stable surface at 400 °C is different between each iteration is different. For high roughness, slight deviation even at 600 °C. Still, the surfaces of highest stability at 700 °C remain the same across iterations, specifically those around 1.1 OH nm⁻².

4. Conclusion

Within this study, a systematic and general approach for the creation of period slab models of amorphous materials is proposed and was used to generate 120 surfaces. These models will help simulate and understand the reactivity of amorphous materials. Dehydroxylated amorphous silica was chosen as a representative example, where we were able to generate 120 surfaces with a variety of roughness, silanol compositions, and silanol densities, aiming to represent the material activated at 700 °C.

The silanol densities of the surfaces were on average 0.96 OH nm⁻² after estimating the surface area through a polygon mesh of the solvent-accessible surface area of Kr. Without calculating this estimate, surface areas would have errors around a factor of 1.7, showing that the previous method of using the xy-cross-sectional area quickly underestimates this property even after the smallest amount of surface roughness is introduced. Not only is the initial silanol density almost exactly that of SiO₂₋₇₀₀, but the numbers of geminal and vicinal silanol groups are in good agreement with experimental results. While vicinal groups are occasionally in excess, the removal of 1 or 2 H₂O molecules is minimal compared to the number that has to be removed following older methodologies.

The plausibility of the models was checked by comparing them against experimentally measured silanol density at given T and P_H2O. For three surfaces of varying roughness, silanol groups were naively added to the surface and the preferred stability was evaluated using atomistic thermodynamics. At 700 °C in vacuum, the preferred silanol density of the model surfaces was between 1.0 OH nm⁻² and 1.1 OH nm⁻². Decreasing the temperature to 600 °C, the preferred stability did not change. Moving to 500 °C and subsequently 400 °C, the preferred stability followed experimental results of 1.8 OH nm⁻² and 2.35 OH nm⁻². Repeating the process for the same models twice more, results remained consistent. Thus, with not only surface silanol characteristics but also thermodynamic properties being accurate for that of SiO₂₋₇₀₀, this ground-up method can be a more computationally efficient alternative for the generation of highly dehydroxylated periodic slabs of amorphous silica. The efficiency of this method, and algorithm, allows one to generate numerous surface models for an amorphous material at a minimal computational cost.

Author contributions

M. P. Klein: investigation; acquire data; coding; formal analysis of data; writing original draft; writing – review & editing A. A. Kolganov: conceptualization, supervision, writing – review & editing E. A. Pidko: writing – review & editing, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

All supporting data, including trajectories, output files, surface structures, and code, will be available at 4TU. Research Data at DOI 10.4121/b8b88018-26a2-41cf-b5cc-10627cb94c0 at the time of publication.

Acknowledgements

The authors thank Nuong P. Nguyen and Prof. Dr Brian B. Laird for providing the force field parameters for silica materials, which greatly facilitated our calculations. The use of supercomputer facilities was sponsored by NWO Domain Science.

Notes and references

1. B. R. Goldsmith, B. Peters, J. K. Johnson, B. C. Gates and S. L. Scott, ACS Catal., 2017, 7, 7543–7557.
2. H. Fritzsche, S. R. Ovshinsky and B. B. Schwartz, New Science, World Scientific Publishing Co. Pte Ltd, 2008, ch. 2, pp. 5–46.
3. P. G. Jeelani, P. Mulay, R. Venkat and C. Ramalingam, Silicon-neth, 2019, 12, 1337–1354.
4. Z. Xu, C. Zhu, H. Shen, Y. Liu, J. Tang, D. A. Weitz and L. Xu, Proc. Natl. Acad. Sci. U. S. A., 2025, 122, e2423426122.
5. W. Qarony, M. I. Hossain, M. K. Hossain, M. J. Uddin, A. Haque, A. Saad and Y. H. Tsang, Results Phys., 2017, 7, 4287–4293.
6. S. Soled, Science, 2015, 350, 1171–1172.
7. Q. Sun, N. Wang, Q. Xu and J. Yu, Adv. Mater., 2020, 32, 2001818.
8. N. Martín and F. G. Cirujano, Org. Biomol. Chem., 2020, 18, 8058–8073.
9. N. LiBretto, Y. Xu, A. Quigley, E. Edwards, R. Nargund, J. Vega-Vila, R. Caulkins, A. Saxena, R. Gounder, J. Greeley, G. Zhang and J. Miller, Nat. Commun., 2021, 12, 2322.
10. A. G. B. Getsoian, B. Hu, J. T. Miller and A. S. Hock, Organometallics, 2017, 36, 3677–3685.
11. A. Panda, A. Alotaibi, I. Kozhevnikov and S. Raveendran, Waste Biomass Valori., 2020, 11, 6337–6345.
12. M. C. Capel-Sanchez, G. Blanco-Brieva, J. M. Campos-Martin, M. P. de Frutos, W. Wen, J. A. Rodriguez and J. L. G. Fierro, Langmuir, 2009, 25, 7148–7155.
13. L. Zhuravlev, Colloids Surf., A, 2000, 173, 1–38.
14. C. S. Ewing, S. Bhavsar, G. Veser, J. J. McCarthy and J. K. Johnson, Langmuir, 2014, 30, 5133–5141.
15. S. D. Fleischman and S. L. Scott, J. Am. Chem. Soc., 2011, 133, 4847–4855.
16. B. Peters and S. L. Scott, J. Chem. Phys., 2015, 142, 104708.
17. M. B. J. Roeffaers, G. De Cremer, J. Libeert, R. Ameloot, P. Dedecker, A.-J. Bons, M. Bückins, J. A. Martens, B. F. Sels, D. E. De Vos and J. Hofkens, Angew. Chem., Int. Ed., 2009, 48, 9285–9289.
18. I. F. Reis and M. H. Gehlen, J. Phys. Chem. Lett., 2024, 15, 2352–2358.
19. M. Caricato, J. Phys. Chem. C, 2021, 125, 27509–27519.
20. B. Rhers, A. Salameh, A. Baudouin, E. A. Quadrelli, M. Taoufik, C. Copéret, F. Lefebvre, J.-M. Basset, X. Solans-Monfort, O. Eisenstein, W. W. Lukens, L. P. H. Lopez, A. Sinha and R. R. Schrock, Organometallics, 2006, 25, 3554–3557.
21. I. Del Rosal, I. C. Gerber, R. Poteau and L. Maron, J. Phys. Chem. A, 2010, 114, 6322–6330.
22. Y. Bouhoute, D. Grekov, K. C. Szeto, N. Merle, A. De Mallmann, F. Lefebvre, G. Raffa, I. Del Rosal, L. Maron, R. M. Gauvin, L. Delevoye and M. Taoufik, ACS Catal., 2016, 6, 1–18.
23. P. Patel, R. H. Wells, D. M. Kaphan, M. Delferro, R. T. Skodje and C. Liu, ACS Catal., 2021, 11, 7257–7269.
24. A. Comas-Vives, Phys. Chem. Chem. Phys., 2016, 18, 7475–7482.
25. P. Ugliengo, M. Sodupe, F. Musso, I. J. Bush, R. Orlando and R. Dovesi, Adv. Mater., 2008, 20, 4579–4583.
26. S. A. Khan, S. M. Godahewa, P. N. Wimalasiri, W. H. Thompson, S. L. Scott and B. Peters, Chem. Mater., 2022, 34, 3920–3930.
27. F. Tielens, C. Gervais, J. F. Lambert, F. Mauri and D. Costa, Chem. Mater., 2008, 20, 3336–3344.
28. A. Galarneau, N. Calin, J. Iapichella, M. Barrande, R. Denoyel, B. Coasne and F. Fajula, Chem. Mater., 2009, 21, 1884–1892.
29. Y. Niu, M. Yu, J. Zhang, Y. Yang, C. Xu, M. Yeh, E. Taran, J. J. C. Hou, P. P. Gray and C. Yu, J. Mater. Chem. B, 2015, 3, 8477–8485.
30. M. A. Smith, A. Zoelle, Y. Yang, R. M. Rioux, N. G. Hamilton, K. Amakawa, P. K. Nielsen and A. Trunschke, J. Catal., 2014, 312, 170–178.
31. D. B. van den Heuvel, E. Gunnlaugsson and L. G. Benning, Geothermics, 2020, 87, 101835.
32. P. N. Wimalasiri, N. P. Nguyen, H. S. Senanayake, B. B. Laird and W. H. Thompson, J. Phys. Chem. C, 2021, 125, 23418–23434.
33. N. P. Nguyen and B. B. Laird, J. Phys. Chem. A, 2023, 127, 9831–9841.
34. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
35. B. W. H. van Beest, G. J. Kramer and R. A. van Santen, Phys. Rev. Lett., 1990, 64, 1955–1958.
36. P. P. Ewald, Ann. Phys., 1921, 369, 253–287.
37. K. Persson, Materials Data on SiO₂ (SG:227) by Materials Project, 2014, An optional note.
38. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
39. B. Siboulet, B. Coasne, J.-F. Dufrêche and P. Turq, J. Phys. Chem. B, 2011, 115, 7881–7886.
40. S. Wells, S. Menor, B. Hespenheide and M. Thorpe, Phys. Biol., 2005, 2, S127–36.
41. W. E. Lorensen and H. E. Cline, SIGGRAPH Comput. Graph., 1987, 21, 163–169.
42. S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, T. Yu and The scikit-image contributors, PeerJ, 2014, 2, e453.
43. S. Lloyd, IEEE T. Inform. Theory, 1982, 28, 129–137.
44. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai and G. Seifert, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 7260–7268.
45. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
46. L. Guimarães, A. N. Enyashin, J. Frenzel, T. Heine, H. A. Duarte and G. Seifert, ACS Nano, 2007, 1, 362–368.
47. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
48. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 65, 035406.
49. D. R. Stull and H. Prophet, JANAF thermochemical tables, 2nd edn, 1971.
50. X. Rozanska, F. Delbecq and P. Sautet, Phys. Chem. Chem. Phys., 2010, 12, 14930–14940.

---

Footnote

† Electronic supplementary information (ESI) available: F-MD parameters, compared stability diagram using SCC-DFTB, Table summarizing surface properties of all generated models, stability diagrams of repeated saturation, structure before and after application of algorithm. See DOI: https://doi.org/10.1039/d5cp01570g

---