Simulation of Forces between Humid Amorphous Silica Surfaces: A Comparison of Empirical Atomistic Force Fields.

Atmospheric humidity strongly influences the interactions between dry granular particles in process containers. To reduce the energy loss in industrial production processes caused by particle agglomeration, a basic understanding of the dependence of particle interactions on humidity is necessary. Hence, in this study, molecular dynamic simulations were carried out to calculate the adhesion between silica surfaces in the presence of adsorbed water. For a realistic description, the choice of force field is crucial. Because of their frequent use and transferability to biochemical systems, the Clay and CWCA force fields were investigated with respect to their ability to describe the water-silica interface in comparison to the more advanced Reax force field, ab initio calculations, and experiments.

Introduction

Granular dry powders usually show poor flowability in process containers. To compensate for the strong shear stresses emerging between the particles, it is necessary to apply extra energy to the system. Therefore, in industrial production processes, one tries to find appropriate means to reduce the energy loss encountered in handling dry granular particles. A basic understanding of particle–particle interactions is essential for a realistic description of the friction between particles and a reliable prediction of possible particle agglomeration causing large clusters and, thereby, hindering flowability. Such interactions are influenced by particle roughness, by adsorbates or chemical impurities on the particle surfaces, and by plastic or elastic deformations in the particle surface layers.

Ultrafine powders consist of particles in the range from a few micrometers to several nanometers. In this study, powders consisting of amorphous silica were considered. Amorphous silica has the advantage that it is an experimentally well-examined substrate. Zhuravlev[1] provided a detailed description of the chemistry of amorphous silica, defining a maximum number of hydroxyl groups on the surface and categorizing them into subgroups according to their chemical bonding to the surface. Gösele and Tong[2] determined surface energies for bonding of silica wafers under moist and dry conditions by dehydroxylating the silica surface and, thereby, establishing siloxane bonds across the interface. Crack measurements were reported by Maszara et al.[3] Adhesion energies of contacting silica surfaces exposed to different relative humidities were determined by crack and atomic force microscopy (AFM) measurements by Wan et al.,[4] Michalske and Fuller,[5] and Chikazawa and co-workers,[6−8] who, additionally reported adsorption isotherms for water on amorphous silica surfaces. Under ambient conditions, water is the most frequently chemi- or physisorbed substance on amorphous silica surfaces. Moisture dramatically changes the adhesion of the particles, leading to enhanced agglomeration of powders, which, in most cases, is undesirable. Hence, it is essential to study particle adhesion for varying relative humidities.

For macroscopic bodies, continuum theories based on the Young–Laplace and Kelvin equations[9] are conventionally used to describe the capillary force that originates from water bridges formed between particles by the condensation of water vapor from the atmosphere. Macroscopic theories have been shown to be successful from the micrometer range[10] upward; however, when nanoparticles or particles exhibiting nanoscale roughness are treated, such theories, which usually rely on the chemical homogeneity and smoothness of the particle surfaces, suffer from several deficiencies.[11] The smallness of the investigated particles (in the range of a few nanometers) together with the requirement to determine highly accurate adhesion forces make such systems ideal for molecular dynamic simulations.

The choice of the atomic interaction potential is crucial to reproduce adhesion forces of realistic silica particles. Aside from computationally costly but very accurate quantum mechanical investigations of the water–silica interface, which are usually restricted to small crystalline silicon oxide clusters reacting with a few water molecules,[12−15] many empirical atomistic force fields have been developed in the past few decades. A review of atomistic bulk silica force fields (without water) was published, for instance, by Pedone.[16] A bulk silica force field that should be mentioned here, however, is the model of Goddard and co-workers[17] because of its later use in this work. It is based on a simple Morse-type short-range potential that includes long-range electrostatic interactions through a Coulomb term. Despite its simplicity, it reproduces the amorphous silica structure very accurately and, moreover, predicts the correct crystal-glass transition temperature, in contrast to the BKS (van Beest–Kramer–van Santen) model.[18]

In addition to the pure silica potentials, several more or less accurate force fields have been developed to study the interaction of silica surfaces with water. The wetting of crystalline silica plates was studied by Berendsen and co-workers[19] using molecular dynamic simulations with the GROMACS (Groningen Machine for Chemical Simulations) force field.[20] The applied force field provides only a crude description of the silica–water interactions; hence, they aimed only to obtain qualitative agreement of the simulated capillary forces with surface tension forces deduced from macroscopic theories. A more sophisticated force field was developed by Garofalini and co-workers.[21−23] They defined an empirical three-body potential for vitreous silica along with a dissociative water model to study hydrophilic wafer bonding and the chemisorption of water on the silica surface. Hassanali and Singer[24] extended the well-known BKS model to investigate the water–amorphous silica interface using classical molecular dynamic simulations. Cole and Payne[25] developed a classical force field based on ab initio calculations to describe the bonding of silica wafers subject to different amounts of adsorbed surface water. Puibasset and Pellenq[26] refined the PN-TrAZ (Pellenq–Nicholson-transferable to all zeolites) potential to study the absorption of water on mesoporous silica. Recently, van Duin and co-workers[27−29] extended their reactive (Reax) force field to describe the silica–water interface. In contrast to most conventional force fields, the Reax force field allows for chemical reactions as well as changes in the polarization of the water molecules and of the atoms on the silica surface, making simulations using it computationally very costly.

Many of the sophisticated force fields mentioned above involve complicated three-body potentials, bond-order terms to account for chemical bonding, or correction terms to screen out unphysical behavior, making them perfectly suitable for the study of pure silica–water systems but hardly transferable to systems containing chemical compounds other than water and silica. The force field of Cole and Payne,[25] the model used by Puibasset and Pellenq,[26] and conventional force fields such as the Clay[30] force field and the Charmm[31]-based CWCA (CHARMM water contact angle) force field of Schulten and co-workers,[32] however, use simple additive pair potentials to describe the atomistic interactions between the silica surface and the water molecules, along with conventional water models such as SPC (simple point charge)[33] and TIP3P (three-site transferable intermolecular potential).[34] Because the Clay, CWCA, and to a certain extent Reax force fields have been parametrized for biochemical simulations (although, because of their computational costs, simulations are restricted to small system sizes), these force fields provide a high degree of flexibility in terms of the incorporation of ions or the decoration of silica particles with organic molecules. This can be very important for the description of practical particle systems. The addition of ions to a colloidal suspension can lead to better dispersion of the particles in solution, for instance, and the coverage of particles with hydrophobic molecules can enhance the flowability of the particles by reducing humidity-induced agglomeration.

Thus, for the water–silica interface, we aimed to test conventional force fields such as the Clay,[30] CWCA,[32] and Reax force fields with respect to their ability to model the amorphous silica–water interface. In particular, the Clay force field has proven to be very successful[35,36] in describing the water structure on crystalline silica in comparison to ab initio simulations and scattering experiments. In the current work, the water structure and hydrogen-bond network on amorphous silica were investigated in comparison to ab initio calculations,[12−15,35−37] crack experiments,[2,4] infrared spectroscopy experiments, (38) and nuclear magnetic resonance (NMR) measurements.[39] In addition, we compared our results with the simulations performed by Cole and Payne[25] in a study of silica wafer bonding.

To be able to judge the results, the nanoparticles were replaced by interacting silica plates, thus mimicking particles with infinitely large diameters. The Clay and CWCA force fields were used to describe the water–silica interface, whereas the bulk silica interactions were described by the model of Goddard and co-workers.[17] In addition, simulations using the Reax force field were carried out, with the force field applied to all interactions in the water–silica interface and in the silica bulk.

The simulations for calculating the interaction of silica plates in humid atmospheres were carried out using the grand canonical Monte Carlo (GCMC) molecular simulation approach.[40] The grand canonical ensemble is conventionally used when the number of molecules in the system is allowed to fluctuate, as is the case during the adsorption of water from the atmosphere. Manzano et al.[28] calculated adsorption isotherms for water on different faces of β-cristobalite, the filling of hydrophilic mesopores and the dynamics of water in confinement[41] using GCMC. Moreover, Ramachandran et al.[42] used GCMC to calculate adsorption isotherms for conventional water models in silicalite pores, Desbien et al.[43] did the same for zeolite and Shirono and Daiguji[44] investigated the phase behavior of confined water in silica nanopores.

This article is organized as follows: First, we describe the simulation methods and the preparation of the silica plates. In the next section, we compare the surface energies of plates subjected to a constant amount of adsorbed water using the Reax, CWCA, and Clay force fields and examine their ability to form hydrogen bonds at the water–silica interface in comparison to the results of ab initio calculations,[35,36] the simulations of Cole and Payne,[25] and experiments.[38,39] Finally, adsorption isotherms for water on silica and surface energies for silica wafers under humid conditions are calculated using GCMC simulations, and the results are discussed in comparison to experiments.[2,4,6]

Methodology and Simulation Details

The simulations were carried out using the open-source parallel code LAMMPS (Large-Scale Atomic/Molecular Massively Parallel Simulator).[45,46]

Preparation of the Amorphous Silica Surfaces

The structure of β-cristobalite served as the starting configuration for the creation of amorphous silica. The crystalline structure of β-cristobalite is conventionally chosen because it undergoes a direct phase transition to the liquid state as its temperature is raised to around 2000 K. In β-cristobalite, the silicon atoms occupy the lattice sites of the zincblende structure, and the oxygen atoms are located between the silicon atoms, forming tetrahedra around them. The cubic unit cell has a size of 7.16 Å and contains 24 atoms, 8 silicon and 16 oxygen, resulting in an oxygen-to-silicon ratio of 2:1.

To model the interactions of two atoms i and j in bulk silica, the Morse-style potential developed by Demiralp et al.[17]was applied using a cutoff of 10 Å with qSi = 1.3e and qO = −0.65e and the parameters given in Table 1. This potential reproduces well the melting temperature and the glass phase transition temperature of silica glass, as well as the density change from 2.33 to 2.2 g/cm3 caused by the transition from the crystalline to the amorphous phase.[47] Thus, it is more reliable than the Clay and CWCA force field, both of which poorly capture the bulk structure of silica at room temperature[35] and are definitely unsuitable at temperatures far above it. The β-cristobalite structure consisted of 5184 atoms in a 42.96 Å cubic simulation box with periodic boundary conditions. The amorphous or glassy samples were prepared by melting the crystal at 7000 K, equilibrating the liquid, and finally quenching it to room temperature at a rate of 4 K/ps using a simulation step of 1.6 fs. These simulations were executed first in the NPT ensemble[48] at a pressure of 1 bar to allow for slight adaptions of the simulation cell volume during the crystal/liquid and liquid/glass phase transitions. Finally, the amorphous sample was equilibrated for 200 ps in the NVT ensemble before plates with a thickness of 20 Å were cut out of the bulk. To minimize the number of dangling silicon and oxygen atoms at the freshly cleaved surfaces, the plates were annealed for 200 ps. This gave amorphous structures with dangling bond densities of less than 1/nm2.

Table 1

Morse Potential Parameters[47]

interaction	R0 (Å)	D0 (kcal/mol)	γ
O–O	3.791	0.5363	10.4112
Si–Si	3.7598	0.17733	15.3744
Si–O	1.628	45.997	8.6342

The partial radial distribution functions and bond-angle distributions (not shown) were found to agree very well with those calculated by Hoang[47] using the same silica potential. Also, the coordination numbers of silicon and oxygens and the fractions of silicon and oxygen atoms found in over-, under-, and chemically saturated coordination match well with the values given in ref (47). Finally, the dangling silicon and oxygen atoms were saturated with hydroxyl groups and hydrogens, respectively. To obtain a silanol density of 3/nm2 for the wafers, siloxane bonds (so-called bridging oxygens) were randomly broken. A silanol density of 3/nm2 was chosen to allow the comparison of our results with those of Cole and Payne[25] who used the same density. The silanols on the silica surface were equilibrated at room temperature for 200 ps, where the hydroxylated structures can be seen in Figure 1. According to Zhuravlev’s model,[1] the silanol density ranges from 2.6 to 4.6/nm2 depending on the preparation process. The silanol groups can be divided into isolated (single silanols), geminal (two hydroxyl groups per silicon), and vicinal or bridged OH groups. These three types of silanols show different hydrogen-bonding affinities with water molecules. Isolated silanols, for instance, have a smaller adsorption energy for water than vicinal or geminal hydroxyls.[49] The type of silanols, not only their density, therefore determines the hydrophilicity of the silica surface. For a hydroxyl density close to 3/nm2, Zhuravlev found that almost all of the silanol groups on the surface are isolated, with only a few being present in geminal or vicinal configuration, in good agreement with our created samples.

Figure 1

Side and top views of a hydroxylated silica plate of thickness 20 Å and silanol density 3/nm2. Silicon atoms in yellow, oxygen atoms in red, and hydrogen atoms in cyan.

Preparation of Wetted Silica Plates

The force fields tested in this work for the water–silica interface are the CWCA force field of Cruz-Chu et al.,[32] the Clay force field,[30] and the reactive (Reax) force field of van Duin and co-workers.[27] The CWCA force field was specially designed to reproduce experimental water contact angles on silica, from which the parameters for the bulk silica atoms were taken, while the ones for the hydroxyl groups were taken from ref (35). Along with the CWCA force field,[32] the TIP3P potential was used to model the water interactions, whereas with the Clay force field,[30] the SPC water model was chosen (see Table 2). The parameters for bond stretching and angle bending interactions are listed in Table 3. For the Clay force field, only hydroxyl (Oh–Hh) bonds and silanol (Si–Oh–Hh) angles were defined, so for the other bonds and angles, the parameters from the CWCA force field were used. For the Lennard-Jones interactions, the standard Lorentz–Berthelot combining rules were used for unlike pairs, where the cutoff was set to 10 Å with an integration step of 2 fs. All hydrogen bonds and angles were frozen through application of the SHAKE algorithm.[50] Moreover, in the course of the Clay and CWCA simulations, the bulk of the plates was kept rigid, and only the silanol groups and their neighboring bulk oxygen atoms were allowed to move.

Table 2

Nonbonded Force Field Parameters for the Clay and CWCA Force Fields from Ref (35)a

	Clay			CWCA
atomb	q	σ (Å)	ε (kcal/mol)	q	σ (Å)	ε (kcal/mol)
Oh	–0.95	3.1655	0.1554	–0.66	3.1553	0.1521
Hh	0.425	0.0	0.0	0.43	0.4	0.046
Si	2.1	3.302	1.84 × 10–6	0.9	3.8264	0.3
Ob	–1.05	3.1655	0.1554	–0.45	3.118	0.1521
Ow	–0.82	3.1655	0.1554	–0.834	3.15061	0.1521
Hw	0.41	0.0	0.0	0.417	0.4	0.046

ULJ = 4ε[(σ/r)12 – (σ/r)6], UCoul = (qq)/(4πεr).

Subscripts b, h, and w stand for bridging oxygens, for atoms part of hydroxyls and water, respectively.

Table 3

Bond Stretching and Angle Bending Parameters for the Clay and CWCA Force Fields from Refs (30 and 51)

bond parametersa					angle parametersb
	CWCA		Clay			CWCA		Clay
bond	Kb	R0	Kb	R0	angle	Ka	r0	Ka	r0
Si–Ob	885.10	1.61	885.10	1.61	Ob–Si–Oh	153.26	111.09	153.26	111.09
Si–Oh	428.0	1.61	428.0	1.61	Oh–Si–Oh	89.62	116.26	89.62	116.26
Oh–Hh	545.0	0.96	554.13	1.0	Si–Oh–Hh	57.50	106.0	30.5	109.47
Ow–Hw	450.0	0.95	554.13	1.0	Hw–Ow–Hw	55.4	104.52	45.76	109.47

Bond function is harmonic with force constant Kb in kcal/(mol Å2) and equilibrium distance R0 in Å.

Angle function is harmonic with force constant Ka in kcal/(mol rad2) and equilibrium angle r0 in degrees.

The Reax force field provides a more advanced description of the silica–water interface than the Clay and CWCA force fields. It uses distance-dependent bond-order functions to represent the contributions of chemical bonding to the potential energy and is thus able to model bond formation and breakage. Moreover, at each simulation step, the partial charges of the atoms of water and silica are adjusted by minimizing the system's electrostatic energy accounting for the polarization of the water molecules and the silica surface. The Reax force field for the silica–water interface reproduces the correct bulk water and bulk silica structures, the binding of a single water molecule to a Si(OH)4 molecule, and the reaction energies for Si(OH)4 polymerization.[27] Thus, the Reax force field is expected to give reliable results for the liquid and vapor bulk phases of water, as well as for the water–silica interface. Recent simulations using the Reax force field[28,29] gave oxygen–hydrogen radial pair distribution functions and X-ray diffraction patterns for calcium silicate hydrate gel in excellent agreement with ab initio calculations. Thus, the hydrogen-bond structure of water on the silica surface seems to be reasonably well captured by this force field. The simulation time step using Reax force field was set to 0.5 fs, leaving the silica plate completely flexible.

In all simulations, the long-range electrostatic interactions were taken into account by applying the particle–particle particle-mesh solver,[52] which is faster than traditional Ewald summation. The equations of motion were integrated using the velocity Verlet algorithm. At the beginning of the simulations, two silica plates were aligned parallel to the x–y plane and placed at a distance of 50 Å from each other along the z axis, thus forming a slit pore. The long-range electrostatic interactions were treated in the slab configuration, that is, periodic in the x and y directions and nonperiodic with respect to z. At the beginning, the pore was completely filled with water, to establish an equilibrated water profile between the plates from which the number of water molecules adsorbed in the first monolayer can be deduced. All filled-pore simulations were carried out in the canonical (NVT) ensemble[48] at a temperature of 298 K using the Clay and CWCA force fields. During equilibration for 500 ps, the width of the slit pore was adjusted in the course of the simulations to guarantee that the pressure did not exceed 1 bar on average and to establish bulk water density in the center of the pore. Finally, a production run for 1 ns was performed to obtain the water profile in the slit pore. The resulting configurations were then used to prepare wetted silica plates by removing the water molecules in the center of the pore, leaving sets with a predefined number of adsorbed water molecules.

Finally, the wetted plates were brought into contact to measure the forces between the silica surfaces as a function of their distance. For this purpose, the silica surfaces were kept at fixed center-of-mass distances along the z axis (perpendicular to the plates), and the mean force between them was measured. The simulations started at a surface-to-surface distance of around 35 Å that was successively reduced by 0.5 Å until the repulsive forces on the silica surfaces became dominant. At each constant surface-to-surface distance, the system was equilibrated for 1 ns, and a production run of 3 ns was conducted to collect the center-of-mass forces on each silica plate. The simulation runs with the Reax force field used 20 ps of equilibration, followed by a production run of 100 ps for each fixed surface-to-surface distance. The shorter simulation times for the Reax force field in comparison to the CWCA and Clay force fields were necessary because of the computational costs involved in charge equilibration. The mean forceacting between plates A and B at a certain separation along the z axis was then deduced as the arithmetic average over the total center-of-mass forces F⃗A and F⃗B collected at a single simulation step.[53] Assuming ergodicity, the arithmetic average corresponds to the canonical ensemble average denoted by the brackets ⟨···⟩ in eq 2. The change in free energy, also called potential of mean force (PMF), encountered by moving the silica plates from an infinitely large separation to the distance of closest approach z0, was calculated by integration of the mean forces along the z axiswhere U is the potential energy of the water–silica system under the constraint that the plates' centers of mass are fixed and F is the force defined in eq 2. In the limit of two interacting plates (studied here), the PMF corresponds to the effective pair potential between the plates.

In addition to systems with fixed amounts of adsorbed water, interplate forces were calculated under ambient conditions using the Clay force field. In the presence of water vapor in the surrounding atmosphere, certain amounts of water molecules are chemisorbed (by the formation of silanol groups) and physisorbed (adsorbed through hydrogen bonds) on the silica surface. The number of adsorbed water molecules depends not only on the partial water vapor pressure or relative humidity, but also on the pore width (plate separation). Thus, the canonical ensemble is inappropriate for modeling the forces between silica plates under humid atmosphere. For these systems, the grand canonical ensemble was applied because it allows for fluctuations in the number of water molecules. Along with the grand canonical ensemble, the inhomogenous phase given by the adsorbed water molecules on the silica surfaces was in equilibrium with an infinite reservoir of water vapor molecules imposing its chemical potential and temperature. The chemical potential is directly related to the pressure of the gas phase, such that, with the choice of a fixed chemical potential, a constant partial water vapor pressure and, with that, a constant relative humidity is set.

Grand canonical Monte Carlo (GCMC) molecular simulations were used to study the influence of humidity on silica–wafer bonding. First, the water vapor pressure as a function of a given chemical potential was simulated. Because of the very low water vapor density, very large simulation boxes with side lengths of 750 Å had to be used for the vapor pressure simulations to obtain system sizes between 50 and several hundred water molecules. Moreover, the adsorption isotherm of SPC water on silica was determined. These simulations consisted of 300000 cycles for equilibration followed by a production run of an additional 100000 cycles to deduce the ensemble average for the number of water molecules. At each cycle, N displacement, insertion, and deletion attempts were executed, where N corresponds to the number of water molecules in the simulation box. The adsorption isotherm provided the starting configurations for wetted noninteracting silica plates at a certain relative humidity. When the silica surfaces were brought into closer contact for force calculations, capillary condensation set in, and the number of water molecules was allowed to adjust for at least 300000 cycles for each pore width. After a stable water meniscus and a constant number of water molecules had been established, the simulations were continued in the canonical ensemble to save computational time during the collection of the forces.

Results and Discussion

Validation of the CWCA, Clay, and Reax Force Fields

In Figure 2, the density profiles of water in an amorphous silica slit pore obtained using the Clay and CWCA force fields are depicted. To clarify the onset of the silica surfaces indicated by the vertical lines in the left plot, the number densities of the atomic contributions are given on the right in the same figure. The origin of the coordinate system was shifted along the z axis in all figures, such that it coincides with the surface of the left silica plate. As can be seen in the figure, the first density peak obtained using the CWCA force field is higher and broader than that obtained using the Clay force field, resulting in a slightly larger number of adsorbed water molecules in the first monolayer (10.5 instead of 9.6/nm2), counting all the waters starting inside the silica bulk up to the first minimum in the density profile. A small amount of water molecules diffuses up to 5 Å into the amorphous silica structure, which explains the higher values in comparison to early assumptions of 9.52 water molecules per square nanometer made by Whalen[54] for the quartz surface. Cole and Payne,[25] for instance, obtained a much higher water number density of 10.7/nm2. As a consequence of the atomistic surface roughness of about 2 Å, on the amorphous plate, the first density peak was hardly elevated with respect to the bulk value in both cases.

Figure 2

Left: Density profiles of water between the amorphous silica plates calculated with the CWCA and Clay force fields. Right: Number densities of the different atomic contributions.

In Figure 3, 3 ML of water was adsorbed on each silica plate. The mean forces calculated for 3 ML using the CWCA and Clay force fields show very similar results as do the potentials of mean force. The PMF was divided by the surface area of the plates, such that its minimum denoting the surface energy can be compared with experiments. In the current case, the CWCA force field gives surface energies in better agreement with the experimental water/vapor surface energy of 72 mJ/m2 that should be encountered when a high amount of water is adsorbed between the plates.

Figure 3

Left: Mean force over plate distance for silica plates covered with 3 ML of water obtained using the CWCA and Clay force fields. Right: Corresponding potential of mean force.

Comparison of Mean Forces and PMF at Low Water Coverage

In empirical force fields such as the CWCA and Clay force fields, a static polarization of the molecules attributable to their neighboring atoms is assumed. This is done, for instance, by assigning higher fixed atomic partial charges (or van der Waals well depths) to the atoms than are observed in a highly diluted gas phase. The CWCA and Clay force fields were parametrized for the bulk water environment. This works fine as long as the chemical surrounding of a certain water molecule or surface atom do not change much in the course of the simulation, that is, far away from the water/vapor interface. Thus, when low numbers of water molecules condense from the atmosphere onto the silica surface, it is questionable whether the static polarization used to mimic the bulk water environment is able to reproduce correctly water–water and water–silica interactions.

To examine this issue further in the limit of low amounts of adsorbed water, the reactive (Reax) force field[27] was chosen in addition to conventional (i.e., nonpolarizable) force fields. In contrast to ab initio simulations, the use of the Reax force field provides the advantage of examining exactly the same surface as was used for the Clay and CWCA force fields, in terms of size, silanol density, and number of adsorbed water molecules. The same argument is also valid in comparison to experiments, where surface roughness, silanol density and type, impurities on the surface, and so on, can hardly be tuned, but considerably change the water adsorption capability. Nevertheless, comparisons with ab initio simulations and experiments were made whenever data were available.

The definitely shorter simulation times used for the Reax force field can be justified because, in the current low-water systems, no large water bridges requiring long production runs were expected to form between the plates. Also, conventional simulations almost converge to the displayed mean forces after 100 ps for each plate separation, so that the nonpolarizable and Reax force field runs are comparable. Moreover, for short Reax force field runs, it is not expected that chemical reactions between the surface silanols and the water molecules take place, which could lead to hydroxylation or dehydroxylation of the surfaces, making the samples less comparable. Around neutral pH, silica is known to exhibit a negative surface charge in water environments. The surface charge can be devoted to a certain number of deprotonated silanol groups. Ab initio calculations[49] have shown that, even at neutral pH, the number of deprotonated silanol hydroxyls and, thus, the surface charge are very low. Bolt[55] derived the surface charge in deionized water at pH 7 to be on the order of −1 μC/cm2. This corresponds to just one deprotonated hydroxyl group per simulation cell in our systems. It is therefore not expected that the neglect of the surface charge strongly influences the calculated forces. In fact, within the short simulation times no deprotonation was observed using the Reax force field.

In Figure 4, the mean forces and surface energies of silica wafers covered with 0.5 ML of water are displayed. For plate distances smaller than 6 Å, the shapes of the mean forces for all three force fields are in quantitatively good agreement. All mean forces in Figure 4 show two more or less pronounced minima, one at around 4 Å and a second at around 6 Å from the surface of the left plate. The global minimum in the mean forces from the CWCA force field is lowest, whereas those deduced from the Clay and Reax force fields nearly coincide. The minimum in the PMF was at a separation close to 3.5 Å in all cases studied. However, the mean forces between the silica plates in the conventional simulations were of considerably longer range than those arising from the Reax force field run. Thus, the surface energies obtained using the Clay and CWCA force fields (105 mJ/m2) are much higher than those for the Reax force field (80 mJ/m2). The longer range of the interactions could be an artifact of the higher water polarization (partial charges) in the Clay and CWCA force fields (see Table 2) in comparison to the Reax force field. On average, the Reax force field charges were found to be 0.38e for the water hydrogen atoms and −0.76e for the water oxygen atoms close to the silica surface.

Figure 4

Left: Mean forces over plate distance for hydroxylated silica plates with 0.5 ML of adsorbed water obtained using the CWCA, Clay, and Reax force fields. Right: Corresponding PMFs.

The surface energies measured[4] on silica glass scatter around 105 mJ/m2, irrespective of the relative humidity between 30% and 95%, showing the same trends as observed in the simulations. However, no information about surface roughness and chemistry was provided in ref (4). A higher surface roughness than in the simulated samples and the presence of hydrophobic adsorbates would diminish the adhesion energy, whereas a higher silanol density or surface defects would increase it, such that the values deduced from the Reax force field and the conventional force fields all lie in the range of possible surface energies.

In Figure 5, molecular water was completely neglected to determine the interactions between dry hydroxylated plates to judge the polarization of the silica surfaces in a vacuum. The mean forces down to a surface separation of 4 Å agree nearly perfectly with each other. The deviations between conventional and Reax forces at closer approach come from the fact that the silica plates are completely flexible using the Reax force field. Thus, the corresponding mean forces are much less repulsive for small distances than for the conventional force fields. Deformations of the silica plates result in a decrease in the equilibrium plate distance for the Reax force field in comparison to the Clay and CWCA force fields. Moreover, this is responsible for the smaller surface energies given by the minimum in the PMF using the conventional force fields. Experiments[5,9] gave surface energies of 25 mJ/m2 for dry surfaces in good agreement with our simulations. Thus, amorphous silica has a polarization that is nearly independent of the surrounding medium.

Figure 5

Left: Mean forces over plate distance for hydroxylated silica plates without the presence of molecular water obtained using the CWCA, Clay, and Reax force fields. Right: Corresponding PMFs.

For comparison with the simulated surface energies of Cole and Payne,[25] the surface energies for a predefined amount of adsorbed water are summarized in Table 4, where it is assumed that the missing CWCA energies perform similar to those deduced from the Clay force field. The simulated surface energies are largest close to a water coverage of 1 ML (compare Table 4), in agreement with the findings of Cole and Payne.[25] The force field of Cole and Payne,[25] however, gives considerably smaller surface energies below a coverage of 0.75 ML than the other empirical force fields. In the dry case, the attractions even vanish, in contradiction with the experimental results.[4,5,9]

Table 4

Calculated Surface Energies (mJ/m2) in Comparison to Those Obtained with the Force Field of Cole and Payne[25]

	coverage/plate
force field	dry	0.1 ML	0.25 ML	0.5 ML	0.75 ML	1.25 ML	3 ML
Clay	27	69	80	105	112	112	95
CWCA	27	–	–	105	–	–	83
Reax	32	–	–	78	–	–	–
Cole and Payne	0	–	19	45	80	90	80

Hydrogen Bonds

For the adhesion energies, the CWCA and Clay force fields gave very similar results with a slight advantage for the CWCA force field at saturated water vapor. This is not surprising because the CWCA force field was parametrized to reproduce experimental contact angles of water droplets on amorphous silica surfaces and, thus, experimental water–silica surface energies at high relative humidity. The Reax force field gives lower surface energies in the presence of water, which might be an artifact of the static polarization assumed for the Clay and CWCA force fields. However, the simulated energies of all three force fields lie within the range of experimental values.

To go into more detail, we aimed to investigate the water structure on the silica surface in a manner similar to that used by Cummings and co-workers[35,36] for the quartz–water interface. For the following investigations, the filled-pore configuration from Figure 2 was taken for the CWCA and Clay force fields, whereas for the Reax force field, a silica plate was covered by 3 ML of water corresponding to a water-layer thickness of around 10 Å to save computation time.

On the left in Figure 6, the radial pair distribution functions for the hydroxyl hydrogens with the water oxygens can be seen, and those for the hydroxyl oxygens with the water hydrogens are depicted on the right. The hydroxyl donor–water acceptor correlations for the CWCA and Clay force fields show nearly the same peak, with the peak maximum from the Clay force field shifted to slightly lower values in better agreement with the Reax force field. The pair distribution function for the Reax force field is more localized than those for the other force fields and coincides nearly perfectly with that found in ab initio calculations by Skelton et al.[35] for the quartz (1011) surface consisting of Q3 silica functional groups only (i.e., distinct silanol groups connected to the surface through three siloxane groups), which makes the system comparable to our amorphous one. The area under the first peak of the radial pair distribution functions (between 1.5 and 2.5 Å) is a measure for the number of established hydrogen bonds. In Table 5, the average numbers of hydrogen bonds per water molecule in the first monolayer from the silica surface are listed. A hydrogen bond is considered to be present between two oxygen atoms when the O–H···O angle is larger than 140° and the O–O separation is smaller than 3.5 Å. It can be seen in Table 5 that the numbers of water oxygens involved in hydrogen bonds with silanol groups are practically the same for all force fields treated here, whereas the simulations of Cole and Payne[25] give a clearly higher number of hydrogen bonds per water molecule.

Figure 6

Left: Radial pair distribution functions of silanol hydrogens with water oxygens comparing the Clay, CWCA, and Reax force fields. Right: Radial pair distribution functions of silanol oxygens with water hydrogens.

Table 5

Average Number of Hydrogen Bonds Per Water Molecule in the First Monolayer for the Three Force Fields and the Force Field of Cole and Payne[25]

force field	Ow–Hw···Ow	Ow–Hw···Ob	Ow–Hw···O	O–H···Ow
Clay	2.4	0.16	0.23	0.21
CWCA	2.5	0.07	0.18	0.20
Cole and Payne	2.54	0.14	0.30	0.32
Reax	2.4	0.37	0.29	0.22

The plot on the right in Figure 6 shows pronounced differences in the trends of the radial pair distribution functions. The Reax force field gives by far the highest number of hydrogen bonds between silanol oxygens and water hydrogens (see Table 5), in good agreement with Cole and Payne. The hydrogen bond length is shifted to higher values as compared to the value originating from the Clay force field, indicating stretched and weaker hydrogen bonds as compared to the hydroxyl donor–water acceptor ones (because the silanol oxygens are not as exposed to water as the silanol hydrogens), which is consistent with ab initio calculations.[35] However, the peak height in the simulations using the Reax force field is somewhat overestimated. The ab initio radial distribution function has a broader, less localized maximum, in better agreement with that obtained using the Clay force field. Finally, the CWCA force field gives the lowest number of hydroxyl acceptor–water donor hydrogen bonds.

On the left of Figure 7, the radial pair distribution functions of water with the siloxane (bridging) oxygens are presented. Close to the silica surface the radial pair distribution functions are declined in comparison to the bulk value which reflects the hydrophobic character of the bridging oxygens. Whereas the radial pair distribution functions from the Reax and Clay force fields are well-structured, the CWCA force field gives hardly any evidence for ordering of the water molecules close to the bridging oxygens. This is also confirmed by the calculation of the number of corresponding hydrogen bonds. Again, the first peak in the radial pair distribution function for the Clay force field is shifted to lower distances compared to that for the Reax force field, indicating a closer approach of the water molecules to the silica surface and, thus, stronger but considerably fewer hydrogen bonds.

Figure 7

Left: Radial pair distribution functions of bridging oxygens with water hydrogens comparing the Clay, CWCA, and Reax force fields. Right: Radial pair distribution functions of vicinal hydroxyls.

As can be seen in Table 5, the numbers of water donor to water acceptor hydrogen bonds with silanol are nearly the same for the Clay and CWCA force fields. Moreover, hydrogen bonding with siloxane bridges is less preferred, both in agreement with results of Cole and Payne.[25] For the Reax force field, generally, water donor hydrogen bonds are formed more frequently than acceptor ones. If one compares the snapshots showing the hydrogen-bonding network on the silica surfaces for the Clay, CWCA, and Reax force fields in Figure 8, one can see that, for the Reax force field, a considerable number of silanols have three hydrogen-bonded waters, two in donor configurations and one in an acceptor configuration, whereas especially for the CWCA force field, at most two waters hydrogen bond to a single silanol group, balancing donor and acceptor contributions. NMR experiments of Chuang and Maciel[39] on silica gel surfaces represented by a β-cristobalite model stated that a single silanol group can form up to three hydrogen bonds, one in donor configuration and two in acceptor configurations, which could explain the observed imbalance in the two kinds of water–silanol hydrogen bonds obtained using the Reax force field. The water for the Reax force field, moreover, has a very high affinity to hydrogen bond to bridging oxygens. The reason for this behavior might lie in the fact that the silica surface is not saturated with hydroxyl groups (where the silanol density usually should amount to approximately 4.6/nm2) such that many bridging oxygens are exposed to the water molecules. It is known[39] that especially pentasiloxane (five Si–O units) and trisiloxane (three Si–O units) rings existing on the silica surface are energetically unfavorable at room temperature. To release strain, these rings are broken by chemisorption of water molecules from the atmosphere to form hydroxyl groups. Because the Reax force field is able to account for bond formation, the high number of water–siloxane hydrogen bonds could indicate the onset of hydroxylation taking place later in the simulation to relax the strained silica surface.

Figure 8

Simulation snapshot of the water hydrogen-bond network after 200 ps on part of the silica surface for the Clay, CWCA, and Reax force fields (from top to bottom). For clarity, only those water molecules involved in hydrogen bonds with the silica surface are shown.

In Figure 8, one can see that the hydrogen-bond networks obtained using the Reax and Clay force fields are much denser than that obtained using the CWCA force field. Whereas there exist only a few inner silanol groups (silanols sitting in the valleys caused by the surface roughness) without hydrogen-bonded water molecules for the Reax and Clay force fields, for the CWCA force field, there are several non-hydrogen-bonded silanols on the silica surface. Moreover, there is a much higher amount of single-coordinated silanols for the CWCA force field than observed for the other two force fields.

In good agreement with the model of Zhuravlev[1] and the simulations of Nangia et al.,[56] the hydrophobic region in the center of the plate in Figure 8 is nearly free of hydrogen-bonded waters for all force fields studied, showing that a closed hydrogen-bonded water layer is not able to form on corrugated not fully hydroxylated amorphous silica surfaces. Close to the silanol groups, a more or less dense hydrogen-bond network builds up, where hydrogen bonds to siloxane oxygens are preferentially formed in the vicinity of destabilizing silanol groups.[39] Only in the case where the hydrophobic regions are small enough are they bridged by hydrogen-bonded water chains that link water molecules involved in hydrogen bonds with surface silanols in different hydrophilic domains.

On the right in Figure 7, the radial pair distribution functions for the correlations between the hydroxyl groups can be seen. They provide a measure for the number of vicinal hydroxyls on the silica surface. All force fields in Table 6 show a small number of hydrogen bonds between neighboring silanol groups, indicating that the surface mainly consists of isolated hydroxyls that are hydrogen-bonded to each other through connecting water molecules. The average number of hydrogen bonds per silanol group using the Reax force field is highest.

Table 6

Average Number of Hydrogen Bonds Per Silanol Group

	hydrogen-bond type
force field	H2O	vicinal	total
Clay	1.45	0.07	1.52
CWCA	1.31	0.10	1.41
Reax	1.85	0.035	1.89

In Table 7, the average numbers of hydrogen bonds per water molecule in the bulk and the first monolayer are listed. The Clay and CWCA force fields show very similar images. The water hydrogen-bond coordination in the bulk is around 3.1, whereas at the surface, it decays to around 3, in agreement with what has been found for the empirical OPLS force field.[57] However, the force field of Cole and Payne[25] and experiments[38] provide evidence that the coordination number of water increases at the silica–water interface in accord with the findings obtained using the Reax force field (see Table 7). Infrared spectroscopy experiments[38] have shown that monolayers of water are present in a quasi-ice-like state on the fully hydroxylated amorphous silica surface with a water–water hydrogen-bond coordination of 3 and a water–silica coordination of 1. If one takes into account that the surfaces in the simulations are not fully hydroxylated in contrast to the experimental ones, the Reax, Cole and Payne, and to a certain extent Clay force fields quantitatively reproduce the observed silica–water hydrogen-bonding affinity, whereas the CWCA force field considerably underestimates it.

Table 7

Total Average Number of Hydrogen Bonds Per Water Molecule in the Bulk and the First Monolayer for the Three Force Fields and the Force Field of Cole and Payne[25]

force field	surface	1 ML	bulk
Clay	0.60	3.0	3.12
CWCA	0.45	3.0	3.15
Cole and Payne	0.76	3.3	3.13
Reax	0.88	3.3	3.14

In Figure 9, on the left, we have plotted the water hydrogen-bond coordination together with its water–water and water–silica contributions for the three force fields with respect to the distance from the silica surface. To see how much the coordination numbers at the respective distance from the silica surface (indicate by z = 0) contribute to the averages for 1 ML given in Table 7, the probability density of water is drawn on the right of the same figure. As already shown in Tables 5 and 7, the water coordination number for the Reax force field increases, whereas those for the two other force fields decay with closer approach to the silica surface. Whereas the water–water hydrogen-bond coordinations in all three systems are nearly the same, the water–silica hydrogen-bond coordination is definitely higher for the Reax force field than for the CWCA and Clay force fields. For the conventional force fields, the loss of water–water coordination at the silica surface is not compensated to the same extent by the formation of water–silica hydrogen bonds.

Figure 9

Left: Total hydrogen-bond coordination (solid line), water–water coordination (dashed line), and water–silica coordination (dotted line) for water at the silica surface indicated by z = 0. Right: Corresponding probability densities.

To summarize, our investigations have shown that, regarding surface energies (see Table 4), the CWCA and Clay force fields perform very similarly, whereas the Reax force field gives lower surface energies in the presence of water. However, the structuring of water at the silica surface obtained using the CWCA force field is less well captured in comparison to that of the Clay and Reax force fields. Cummings and co-workers[35,36] came to very similar results describing the water–quartz interface, judging the Clay and CWCA force fields with respect to quantum calculations. Well-structured water with the ability to form a dense hydrogen-bond network on the silica surface seems to be the more physical description of the system. Thus, the Clay force field was used for a final verification. Because we were interested in the behavior of wetted silica surfaces, we aimed to judge simulated surface energies under ambient conditions in comparison to experiments.

Interaction of Humid Plates

On the left in Figure 10, the water vapor pressure and density are depicted as functions of the applied chemical potential. For the ideal-gas approximation, the expressionswere taken to calculate pressure P and density ρ from the chemical potential μ, where Λ = (h/2πmkBT)1/2 is the de Broglie wavelength and h is Planck's constant. In addition to the analytical values, the gas density and pressure were adjusted by grand canonical Monte Carlo molecular simulations given a certain chemical potential. As can be seen in Figure 10, the simulated values did not deviate strongly from those deduced with the ideal-gas assumption. From simulated phase diagrams,[58] the saturated vapor pressure of SPC water at T = 300 K was found to be P0 = 0.044 bar. Given the saturated vapor pressure, the partial vapor pressure P/P0 or the relative humidity in dependence of the chemical potential can easily be obtained.

Figure 10

Left: Density and pressure of SPC water vapor for varying chemical potentials. Right: Adsorption isotherm of SPC water on the silica surface as a function of relative humidity (where 16.0 μmol/m2 corresponds to 1 ML of adsorbed water) in comparison to that simulated in ref (26) (magenta spheres).

On the right in Figure 10, the adsorption isotherm for SPC water on a single amorphous silica plate is depicted with respect to the partial vapor pressure (relative humidity). The adsorption isotherm is in considerably good agreement with that determined by the GCMC simulations of Puibasset and Pellenq,[26] who approximated the adsorption isotherm of an amorphous silica plate by the superposition of adsorption isotherms calculated for different faces of the β-cristobalite crystal. Moreover, they[26] showed that their adsorption isotherm coincides nicely with experiments of Takei et al.[8] The number of water molecules per square nanometer in the first monolayer, found to be 9.6 in the last section, corresponds to the amount of 16.0 μmol/m2.

In Figures 11–13, the mean forces and corresponding PMFs of wetted wafers are depicted for partial pressures P/P0 of 0.8, 0.55, and 0.3, respectively, corresponding to amounts of adsorbed water molecules on the noninteracting plates of around 1, 0.75, and 0.5 ML, respectively. We denote these amounts in the following discussion as equilibrium surface coverage. If the plate separation becomes close enough that the water molecules at the different plates can interact with each other, capillary condensation sets in, which is characterized by the formation of a water meniscus between the contacting surfaces. The water in the water bridge condenses out from the surrounding atmosphere, where its surface tension and the adsorbed water layers on the silica surfaces are in equilibrium with the partial vapor pressure. The higher the relative humidity, the more water condenses into the slit pore. A stable water bridge is formed for distances between zmax, marking the onset of a small connecting water path, and zmin, where the pore starts to be continuously filled with water. The capillary excess giving the number of condensed water molecules in the bridge (excluding those of the equilibrium surface coverage) rises slowly as soon as zmax is reached, runs through a maximum close to zmin, decays for smaller pore width, and becomes zero around the equilibrium plate distance, which roughly corresponds to twice the thickness of the equilibrium surface coverage. This can be seen in Figure 14, where the number of adsorbed water molecules is presented as a function of plate separation and relative humidity. The dotted lines indicate the equilibrium coverage, whereas the dashed lines mark the positions of zmax and zmin. As one would expect, for high relative humidity, the pore starts to fill with water already at larger separations (location of zmin) than for low relative humidity. So, for instance, at 30% relative humidity, it fills at a pore width of around 7 Å, whereas for intermediate relative humidity, this happens at 12 Å, and for high relative humidity, it occurs already at around 16 Å. Also, the onset of the water meniscus (location of zmax) is shifted to larger distances with rising relative humidity. As soon as zmin is reached, water condensation is restricted by the available volume between the plates, such that the numbers of water molecules in the slit for small distances coincide, irrespective of the relative humidity. Finally, when the pore becomes so narrow that its width falls below the thickness of the equilibrium coverage, water is squeezed out. However, because of the amorphous character and surface roughness of the silica plates, water is not completely removed from the gap between the plates. Even at closest approach, traces of water remain in the slit.

Figure 11

Left: Mean force split into its contributions from the water–plate and plate–plate interactions at 80% relative humidity. Right: Corresponding PMFs.

Figure 13

Left: Mean force split into its contributions from the water–plate and plate–plate interactions at 30% relative humidity. Right: Corresponding PMFs.

Figure 14

Number of water molecules as a function of pore width and relative humidity. Dotted lines mark the equilibrium coverage, whereas dashed lines indicate the positions of zmax and zmin.

Left: Mean force split into its contributions from the water–plate and plate–plate interactions at 55% relative humidity. Right: Corresponding PMFs.

The onset of the water bridge at large plate separations and high relative humidity can be seen by a step in the mean forces in Figure 11, indicating a sudden increase in the attractive forces, resulting from the negative capillary pressure inside the meniscus, pulling the surfaces together. For low relative humidity, one cannot strictly speak of a continuous capillary bridge, but rather of single hydrogen-bonded water chains linking the silanol groups at opposite plates to each other. The step in the onset of the mean forces at large distances disappears, and the curve becomes smooth. The mean forces and PMFs in Figures 11–13 were split into their contributions arising from capillary and silica–silica interactions. The water-mediated or capillary forces dominate the interactions between the plates, whereas the range of direct silica–silica interactions, which originate mainly from van der Waals forces, is limited to a plate separation of around 6 Å. Thus, for high and intermediate relative humidities, the water layers formed between the plates prevent the surfaces from interacting directly with each other. When the surface distances fall below zmin, where the water bridge becomes a continuous water film, the mean forces oscillate, showing minima close to 4, 6, 8, and 12 Å. The oscillations are caused almost entirely by the water-mediated forces alone. The oscillations in the mean forces indicate the formation of water layers between the plates. At a distance of around 3.5 Å, the plates are separated by a single water layer; at 6 Å, by two water layers; at 8 Å, by three water layers; and so on. What is meant by water layering is illustrated more clearly in Figure 15. For the system at 30% relative humidity, at a larger distance of 8 Å (right figure), water chains consisting of, on average, three water molecules connect the surfaces by forming hydrogen bonds with silanol groups on opposite plates and among each other. Then, these hydrogen bonds are broken by moving the plates closer together, and others are established at closer approach. At the equilibrium distance of around 3.5 Å, the surfaces are preferentially bridged by single water molecules hydrogen-bonded to silanols, seen in the snapshot on the left in the same figure. Moreover, a considerable number of silanols directly form hydrogen bonds across the plates.

Figure 15

Snapshots of the hydrogen-bond network for surface–surface distances of 3.5 and 8 Å at 30% relative humidity after 1 ns of simulation. (For clarity, only slices of the silica slit are shown.)

For a better comparison, the mean forces subject to different amounts of adsorbed water are drawn in a single plot on the left in Figure 16, where, with decreasing water-layer thickness and relative humidity, the minima in the forces indicating the pull-off forces are shifted to lower plate separations. The striking difference in the force curves is that the mean force becomes less attractive with rising relative humidity (excluding the dry case). This effect is especially visible in the maximum of the mean force at around 5 Å, where, for high and intermediate relative humidities, the mean force becomes repulsive, whereas for low relative humidity, it stays attractive. For high relative humidity, water vapor has a strong affinity to condense out between the hydrophilic plates. For small plate separations, the silica–silica interaction attempts to pull the plates together. If the equilibrium coverage on both plates becomes broader than the gap between the plates, the silica–silica force acts against water condensation. This causes a repulsive capillary force (seen for relative humidities of 55% and 80%) that keeps the plates at larger equilibrium separations than observed for low relative humidity.

Figure 16

Left: Mean forces over distance for varying relative humidity. Right: PMFs for varying relative humidity.

If one examines the local minima of the mean forces for dry and humid plates, one can see that the pull-off forces run through a maximum with respect to the relative humidity, in agreement with experiments. The maximum of the pull-off forces for our four selected systems lies at around 30% relative humidity or 0.5 ML of equilibrium surface coverage. For this relative humidity, the agglomeration between infinitely large particles would be strongest, stronger than for dry powders and much stronger than for completely wetted particles. However, to draw a picture of pull-off forces against relative humidity and to provide a statement about the absolute maximum of the pull-off force, more systems would have to be calculated.

On the right in Figure 16, the PMFs for wafers are shown at varying relative humidities. The location of the minimum in the PMF (surface energy) moves from a plate–plate separation of near 6 Å at 80% relative humidity down to 3.5 Å for dry surfaces. So, the plates at low relative humidity have a distance of 3.5 Å separated by a single water layer; those at intermediate relative humidity are 5 Å from each other, separated by around 1.5 ML of water; and the plates for high relative humidity have their equilibrium position at 6 Å from each other, separated by around 2 ML of water that have been merged into a single broad layer. The surface energies given by the local minima in the PMF are 100, 90, and 82 mJ/m2 for 30%, 55%, and 80% relative humidities, respectively. Crack experiments by Wan et al.[4] gave nearly constant surface energies scattering around 105 mJ/m2 between 30% and 95% relative humidity, which is in the range of the simulated values.

Summary

PMFs and MFs calculated for dry and humid silica wafers subjected to different amounts of adsorbed water were found to show close agreement between the Clay and CWCA force fields. Thus, a closer inspection of the water structure in the water–silica interface should shed more light on the reliability of the investigated force fields. It is well-known[1] that hydroxyl groups on the silica surface show a high affinity to form hydrogen bonds with neighboring water molecules. The ability to reproduce the hydrogen-bonded water structure close to the silica surface as realistically as possible is crucial for a physical description of the capillary forces.

To investigate this issue, partial radial distribution functions that describe the correlations of the atoms on the silica surface with water were calculated, together with the average numbers of hydrogen bonds per water molecule using the Clay, CWCA, and Reax force fields. The water acceptor–silanol donor radial pair distributions are nearly the same for the Clay and CWCA force fields, and that for the Reax force field agrees almost perfectly with quantum calculations done by Skelton et al.[35] Despite small differences in the radial distribution functions, the number of formed silanol–water hydrogen bonds is rather the same for all three force fields. The water donor–silica acceptor radial pair distribution functions, however, showed dedicated differences. Whereas the radial distribution function for the Clay force field agrees nicely with ab initio ones,[35] the peak height of the radial distribution function was overestimated for the Reax force field and underestimated for the CWCA force field. NMR experiments[39] on different β-cristobalite faces have pointed out that surface silanols can encounter up to three hydrogen bonds, with one in a donor configuration and two in acceptor configurations. This could explain the higher portion of water donor–silanol acceptor hydrogen bonds (in comparison to the water acceptor–silanol donor ones) observed for the Reax force field. The radial distribution function between bridging oxygens and water hydrogens shows good structuring of water on the silica surface for the Clay and Reax force fields, whereas the CWCA force field gives hardly any evidence of hydrogen bonding between the mentioned atoms. In the case of the Reax force field, those hydrogen bonds are weak and elongated. The larger number of water–siloxane hydrogen bonds in comparison to the other force fields could be attributed to unreleased strain in the silica surface, where the hydrogen bonds mark the onset of hydroxylation taking place later in the simulation.

In agreement with the simulations of Cole and Payne,[25] the average hydrogen-bond coordination per water molecule for the Reax force field rises, approaching the silica surface coming from the bulk. However, the CWCA and Clay force fields indicate a decay of the water coordination in the vicinity of the silica surface, similar to what was found for the empirical OPLS force field.[57]

Whereas the average number of water–water hydrogen bonds is nearly the same for all three force fields, differences in the total water hydrogen-bond coordination arise from the number of water–silica hydrogen bonds. The average hydrogen-bond coordination of water with the atoms on the silica surface in the first monolayer is higher for the Clay force field (0.60) than for the CWCA force field (0.45) and in better agreement with the findings using the Reax force field (0.88) and the Cole and Payne force field (0.76), which are supported by infrared spectroscopy experiments[38] that give a water–silanol coordination of 1 on fully hydroxylated silica surfaces. Also, in comparison to quantum calculations[35,36] and NMR measurements[39] the water structuring on the silica surface is better captured by the Reax and Clay force fields than by the CWCA force field.

The water–silica interface is better described by the Clay force field than by the CWCA force field; hence, for further analysis, the Clay force field was chosen to calculate surface energies between humid wafers in comparison to experiments[2,4] and recent simulations.[25] In first instances, we restricted the mean force calculation to interacting silica plates subjected to a fixed amount of adsorbed water, to compare the calculated surface energies with those found by Cole and Payne.[25] For higher water coverages, the simulations of Cole and Payne are in good qualitative agreement with our results. Cole and Payne[25] found that the adhesion energy between the wafers reaches a maximum of 90 mJ/m2 for one monolayer of adsorbed water molecules, whereas our simulations gave a maximum of around 112 mJ/m2 using the Clay force field at the same water-layer thickness. The Cole and Payne force field, however, gave considerably lower surface energies in the case of lower relative humidity than obtained for all three empirical force fields used in this study. For dry wafers, for instance, the adhesion energy vanished, in contradiction with experiments giving around 25 mJ/m2,[5] which is nicely reproduced by our simulations. This result provides evidence that the force field of Cole and Payne underestimates the wafer interactions at low water coverage.

To study the interaction of humid silica plates, a GCMC molecular simulation was performed, to account for the fluctuations in the number of water molecules when capillary condensation set in. The adsorption isotherm calculated for SPC water on amorphous silica is in good agreement with those deduced by Puibasset and Pellenq[26] using the PN-TrAZ potential.[41] Moreover, the agreement with experimental adsorption isotherms[8] is reasonable. Also, the surface energies show considerably good agreement with experiments.[2,4] The experiments of Gösele and Tong[2] on fully hydroxylated amorphous silica gave surface energies between 100 and 200 mJ/m2 depending on storage time of the interacting humid silica plates. The experiments of Wan et al.[4] showed nearly constant adhesion energies of around 105 mJ/m2 over a wide range of relative humidities between 30% and 95%. The simulated surface energies have values of 100, 90, and 82 mJ/m2 for 30%, 55% and 80% relative humidity, respectively, and, thus, are close to the measured values, especially if one takes into account that the simulated silica surface has a lower silanol density and, consequently, a lower surface energy than the surfaces used in most experiments. Whether the surface energy reaches a value of 72 mJ/m2 for saturated water vapor corresponding to the surface tension of water at the liquid–vapor interface cannot be answered; however, the surface energy of 82 mJ/m2 for 80% relative humidity comes close to this value.

Conclusions

Because of their widespread use and transferability to a variety of chemical compounds (apart from water and silica), the Clay[30] and CWCA[32] force fields have been tested with respect to their applicability to describing the silica–water interface in comparison to the more advanced Reax force field,[27] ab initio calculations,[35] and infrared spectroscopy experiments.[38] Although a static polarization was used with the Clay and CWCA force fields, which leads to the known underestimation of the water-vapor surface tension for SPC and TIP3P water, the Clay and CWCA force fields gave surface energies for humid silica plates in reasonable quantitative agreement with crack measurements[4] and simulations using the Reax force field. Also, adsorption isotherms deduced by GCMC molecular simulations using the Clay force field showed reasonable consistency with experiments.[8]

Moreover, the ability to form a highly coordinated water hydrogen-bond network in the silica–water interface was tested for the three force fields. The Reax force field gives a silica–water hydrogen bond number in good agreement with experiments by forming, however, a large number of siloxane–water hydrogen bonds, which are energetically less preferred. For the Clay and CWCA force fields, the number of water–silica hydrogen bonds was underestimated in comparison to the force field of Cole and Payne[25] and infrared spectroscopy experiments. However, the water–silica radial distribution functions obtained using the Clay and Reax force fields agree considerably better with ab initio calculations than those obtained using the CWCA force field. Taking into consideration that the surfaces in the simulations were not fully hydroxylated, the Clay force field gives a physically meaningful description for the ordering of the water molecules on the silica surface in qualitative agreement with the Reax force field. Thus, despite its simplicity, the Clay force field provides feasible accuracy for the investigation of humid amorphous silica surfaces, leaving space for further refinement.