ARTICLE pubs.acs.org/JPCC
Simulations of the Quartz(1011)/Water Interface: A Comparison of Classical Force Fields, Ab Initio Molecular Dynamics, and X-ray Reflectivity Experiments A. A. Skelton,*,† P. Fenter,‡ J. D. Kubicki,§ D. J. Wesolowski,|| and P. T. Cummings||,† †
Department of Chemical Engineering, Vanderbilt University, Nashville, Tennessee 37240, United States Chemical Sciences and Engineering, Argonne National Laboratory, Argonne, Illinois 60439, United States § Department of Geosciences, The Pennsylvania State University, University Park, Pennsylvania 16802, United States Chemical Sciences Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, United States
)
‡
bS Supporting Information ABSTRACT: Classical molecular dynamics (CMD) simulations of the (1011) surface of quartz interacting with bulk liquid water are performed using three different classical force fields, Lopes et al., ClayFF, and CHARMM water contact angle (CWCA), and compared to ab initio molecular dynamics (AIMD) and X-ray reflectivity (XR) results. The axial densities of the water and surface atoms normal to the surface are calculated and compared to previous XR experiments. Favorable agreement is shown for all the force fields with respect to the position of the water atoms. Analyses such as the radial distribution functions between water and hydroxyl atoms and the average cosine of the angle between the water dipole vector and the normal of the surface are also calculated for each force field. Significant differences are found between the different force fields from such analyses, indicating differing descriptions of the structured water in the near vicinity of the surface. AIMD simulations are also performed to obtain the water and hydroxyl structure for comparison among the predictions of the three classical force fields to better understand which force field is most accurate. It is shown that ClayFF exhibits the best agreement with the AIMD simulations for water hydroxyl radial distribution functions, suggesting that ClayFF treats the hydrogen bonding more accurately.
’ INTRODUCTION The interaction of quartz (R-SiO2) with water is of interest in geological,1-4 biological,5,6 and technological7,8 contexts. Quartz is a significant component of soils, clays, and rocks constituting about 20% of the Earth’s exposed crust.9,10 Moreover, the properties of amorphous silicas involving high surface area gels and colloidal particles11,12 depend on the silica surface chemistry.13 The interaction of quartz and other silica phases with a range of biomolecules is relevant to diseases such as silicosis14 and in the controlled assembly of nanocomposites.15 Water is important as it is ubiquitous in natural environments, and it mediates the interaction of biomolecules with surfaces.16 Quartz is challenging to study both experimentally and theoretically because, unlike many other metal oxides, it exhibits significant solubility in water. To add to the complexity, it has been found that the dissolution rate of quartz depends on the solution pH and dissolved salt concentration in aqueous solutions.17-21 Quartz also exhibits surface acid/base properties; it has been suggested that the pH at the point of zero charge r 2011 American Chemical Society
(PZC) for quartz is approximately 2.22 At higher pH, deprotonaton of surface silanol groups (dSiOH) leads to the development of negative surface charge.23-27 In the present study, classical molecular dynamics (CMD) simulations will be compared to previous X-ray reflectivity (XR) experiments28 as well as ab initio molecular dynamics (AIMD). In the aforementioned XR study, two surfaces of R-quartz, (1010) and (1011), were investigated. The study provided information about the atomic-scale structure of the quartzwater interface, wherein the electron density normal to the surface in the interfacial region is converted to the axial density of atoms at the surface and the distribution of water in direct contact with the surface, which indicates interfacial water layering (we use the term “axial” to refer to profiles normal to the surface). The atomic distribution and electron density can be extracted Received: October 1, 2010 Revised: December 10, 2010 Published: January 18, 2011 2076
dx.doi.org/10.1021/jp109446d | J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C from CMD, provided valid force fields and structural models are available. A direct comparison between the XR results and molecular dynamics will therefore be made for the (1011) surface. The present work follows on from previous studies involving the oxide surfaces, titania (TiO2) and cassiterite (SnO2), which used a similar combined CMD/ab initio/XR approach29-33 to study the interfacial water structure. Previous studies involving CMD simulation of the silicawater and quartz-water interfaces using various classical force fields have been reported.34-43 ClayFF, a force field for studying the interaction of clays with water, was parametrized by Cygan et al.44 The water structure, orientation, and dielectric properties were investigated for the (0001) surface of quartz using ClayFF, and it was shown that simulation using ClayFF compared favorably with dielectric spectroscopy data.45 A force field for quartz-water interactions was developed and tested by Lopes et al.46 The water structure and dynamics were analyzed for the (1010) and (0111) surfaces. The vibrational density of states (VDOS) was calculated from simulations of water confined between two quartz surface slabs in close proximity using the Lopes et al. force field. A comparison was made with VDOS obtained by neutron scattering experiments involving water confined within various types of porous silica, Vycor glass,47 and sol-gel silica glass.48 Agreement between experiment and simulation indicated the validity of the Lopes et al. force field for describing the vibrational structure of water at silica surfaces. Further water structuring at the (1010), (0001), and (0111) surfaces was investigated using the Lopes et al. force field.49 A separate force field was parametrized for modeling amorphous silica,50 called the CHARMM water contact angle (CWCA), and was used to study the permeation of water through silica nanopores. CWCA was used to calculate electrokinetic properties of electrolyte solutions in a silica nanochannel.51 Ab initio molecular dynamics (AIMD) simulations can be especially important for modeling the quartz-water and silicawater interfaces as a test of CMD simulations, for dynamically modeling chemical events, and for providing key data to calibrate classical force fields.52 A few studies have been made for water on quartz surfaces. AIMD (BLYP, 952 eV cutoff) was used53 to simulate the extent of hydroxylation of the nonhydroxylated (0001) surface when exposed to bulk water, and free energy calculations of the dissolution of a Q2 silica unit (one connected to the surface via 2 Si-O-Si bonds) were performed by constrained molecular dynamics.54 A detailed mechanism of hydroxylation of a nonhydroxylated (0001) surface was illuminated by AIMD (PW91, 680 eV cutoff), and the hydrogenbonding pattern of hydroxyls on (0001) for the clean surface as compared to water-covered surface was studied (PW91, 300 eV cutoff).55 As described above, considerable information has already been gathered about the quartz-water interface via simulation. As far as CMD is concerned, however, the robustness of the force fields is not well-defined, that is, the degree to which the results depend on the choice of the force field. This can be assessed by comparing different available force fields for a given system. This is the focus of the present study, where we compare structures obtained with three force fields (Lopes et al.,46 ClayFF,44 and CWCA50,51,56) to AIMD simulations. The quartz (1011) surface was chosen for this study because the structure of this interface has been obtained from XR measurements of this system for comparison and testing of the computational results. Because the force fields available are currently limited to modeling neutral
ARTICLE
quartz surfaces, throughout the present study, the quartz/water systems will be modeled at the PZC of quartz (acidic pH).
’ METHODOLOGY The Quartz (r-SiO2) (1011) Surface. Figure 1 shows a side and top view of the quartz (1011) surface. The surface was obtained by cleaving the bulk crystal at the corresponding crystallographic plane. It was hydroxylated by adding hydroxyl groups to silicon dangling bonds to account for dissociation of water.53,57 It should be noted that this is an idealized structure and the structure of the real surface structure is not fully known. The quartz (1011) surface contains only Q3 silica functional groups, which are connected to the surface via three siloxane (Si-O-Si) bonds. These Q3 groups are expressed as two distinct silanol (Si-O-H) groups, distinguished by their height above the surface (indicated in Figure 1). The surface siloxane linkages (also indicated in Figure 1) are significant because they are the first bonds that must be broken for dissolution to occur (i.e., release of silicic acid into solution). Force Fields. Molecular dynamics of the quartz (1011) surface was performed, testing three different force fields. These are the Lopes et al. designated here as LFF,46 ClayFF,44 and CWCA50,51,56 force fields. LFF was designed to model quartz and quartz-water interactions specifically. It was parametrized from ab initio calculations (model silica compounds with water) and involves explicit bonds, angles, and dihedrals. Specifically, the Lennard-Jones (LJ) parameters and partial atomic charges were optimized by reproducing the interaction energies between water and the model compounds. LFF uses the modified Transferable Intermolecular Potential 3P (TIP3P) water model.58 ClayFF was designed to model clays and contains no special bonding terms except for those between the hydrogen and oxygen in hydroxyls. LJ parameters were obtained by fitting to structures of model compounds quartz and kaolinite, although all oxygens were assigned the same LJ parameters as the oxygen of TIP3P water.44 Partial charges were initially assigned by reference to Mulliken and electrostatic charge determinations obtained from single point electronic structure calculations from periodic functional theory of the model compounds. Quantum calculations of small silica clusters were then performed, and electronic distributions were analyzed to further refine the partial charges as well as examine the variation of charge with structure. ClayFF uses the flexible single point charge (flexible SPC) water model.59 The CWCA force-field derived its parameters from many sources, and the parameters are given in Lorenz et al.51 The partial charges of the bulk oxygen and silicon and LJ epsilon of the bulk silicon were refitted to reproduce contact angle experiments on amorphous silica surfaces.50 The partial charges and LJ parameters of the hydroxyl hydrogen and hydroxyl oxygen were obtained from the hydroxyl of serine in the CHARMM protein force field.60 CWCA also contains bond and angle parameters taken from Hill and Sauer.56 The silicon/hydroxyl oxygen bond parameters of Lorenz et al.51 were given with an equilibrium value, r0, of 1.42 Å where the bonding potential Ubond is:
Ubond ¼ Σkbond ðri - r0 Þ2 i
ð1Þ
kbond is the bond spring constant, ri is the bond distance, and r0 is i the equilibrium bond distance. However, considering the bonding parameters were based on the potential of Hill and Sauer,56 2077
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
ARTICLE
Figure 1. Side (left) and top (right) views of the quartz (1011) surface.
Table 1. Intermolecular Parameters of the Three Force Fields, LFF, ClayFF, and CWCAa LFF
CWCA
q
σ
ε
q
σ
Oh
-0.54
3.1538
0.1521
-0.95
3.1655
0.1554
Hh
0.32
0.4
0.046
0.425
0.0
0.0
0.43
0.4
0.046
Si
1.08
3.91996
0.6
2.1
3.302
1.84 10-6
0.9
3.8264
0.3
-0.53 -0.834
3.1538 3.15061
0.1521 0.1521
-1.05 -0.82
3.1655 3.1655
0.1554 0.1554
-0.45 -0.834
3.118 3.15061
0.1521 0.1521
0.4
0.32
0.0
0.0
0.4
0.32
Obulk Ow a
ClayFF
Hw
0.417
0.41
ε
q
σ
ε
-0.66
3.1553
0.1521
0.417
q is the Coulombic partial charge in electrons; σ and ε are the LJ collision diameter (in Å) and well depth (kcal mol-1), respectively. Oh are hydroxide oxygens, Hh are hydroxide hydrogens, Ow are water oxygens, Hw are water hydrogens, and Obulk are oxygens in the bulk of the surface slab.
which gives the r0 as 1.6128 Å (a higher order harmonic was used but is effectively equivalent), the r0 was changed to 1.61 Å, equal to that of the silicon-bulk oxygen value in the Lorenz et al. work.51 Table 1 shows the intermolecular parameters (partial charges and LJ) for the different force fields. It is interesting to compare the parameters because they will dictate the water structuring at the surface. First, it can be seen that the LJ parameters are quite similar for all force fields, the main differences being in the hydroxyl hydrogens (Hh), which have nonzero values of σ and ε (describing the LJ collision diameter and well-depth, respectively) for the LFF and CWCA, but are zero for ClayFF. This mimics the water hydrogen (Hw) of TIP3P used with CWCA and LFF and Hw of SPC used with ClayFF. Also, the Si LJ parameters are different between force fields, but it is not expected that this would have a significant effect considering the water does not directly come in contact with the silicons. The partial charges are different between models, and this could be expected to have a significant impact because the water partial charges are the same or similar for the water models (SPC and TIP3P). This is especially relevant for the hydroxyls, which will come in direct contact with water. For the hydroxyl oxygen (Oh), ClayFF has the largest partial charge, followed by CWCA, with LFF having the smallest partial charge. ClayFF and CWCA have similar Hh partial charges but quite different Oh charges. It is
possible that the different relative magnitudes of the Hh and Oh charges may have consequences for the water structure. The bulk silicon (Si) and bulk oxygens (Obulk) are also different between force fields with ClayFF having the largest partial charge (about double that of LFF), followed by LFF, followed by CWCA. The charges of Si and O atoms in the bulk crystal, despite not coming into direct contact with water, may still significantly influence the water structuring due to long-range interactions. Ab Initio Molecular Dynamics (AIMD). All AIMD calculations reported here were performed using the Vienna ab initio Simulation Package (VASP)61 using an energy cutoff of 400 eV. The Perdew-Burke-Ernzerhof (PBE) exchange and correlation functionals62 and the projector-augmented wave (PAW) method63 were used in all density functional theory (DFT) energy calculations. The initial configurations for the VASP quartz-water model systems were obtained as follows: the bulk quartz structure64 was cleaved at the (1011) surface. The surface was manually converted to a silanol (Si-O-H) terminated surface by adding and dissociating H2O molecules to underbonded Si and O atoms on the cleaved surface. An energy minimization was performed using VASP. A periodic box of H2O molecules was created and inserted into the vacuum space of the quartz (1011) surface model. Another energy minimization was then performed for the surface and water model before the AIMD simulations were begun. Cell dimensions were 2078
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
ARTICLE
Figure 3. A snapshot showing the two types of interactions (a) Hw-Oh and (b) Ow-Hh. The ball and stick representation indicates the surface and two specific waters. The stick representation indicates all other waters. The dotted lines indicate hydrogen bonds.
Figure 2. The configuration of the cell used for AIMD and corresponding CMD calculations under the same conditions. The balls show the atoms that are allowed to move in the CMD calculations, and the sticks indicate the fixed atoms.
a = 13.755 Å, b = 9.829 Å, c = 41.5 Å, and contained 60 Si atoms, 236 O atoms, and 232 H atoms, 108 water molecules, and 204 atoms in the surface. Reciprocal space was sampled with a Monkhorst-Pack 1 1 1 grid65 for all supercells. Lattice parameters for these calculations were fixed at the experimental values for R-quartz. Simulations were performed with a time step of 1 fs, run for 5 ps. During all minimizations and MD simulation, the central atoms were held fixed to maintain the perfect crystal structure.64 Figure 2 shows which atoms were allowed to be flexible (balls) and fixed (stick representation). Classical MD Details. The CMD calculations were performed using the same conditions as those used for the AIMD calculations to make a direct comparison of the results, because the AIMD cell was much smaller and less real time was simulated than is typically simulated using CMD. This involved using the same initial configuration, time step, and length of simulation time. The degree of noise due to sample size will, therefore, be the same in both systems. In the case of the small scale CMD simulations, a cutoff of 4.8 Å (one-half the box) was used for electrostatic and van der Waals interactions. Radial distribution functions (RDFs) of Hh and water oxygen (Ow) and the RDF of Oh and Hw were used to compare the simulated water structures. The radial distribution functions measure the number of particles of a particular type as a function of distance from a reference atom of a particular type per
configuration (e.g., Oho-Hw indicates the distribution of Hw around the outer Oh site). This number is normalized by the volume of the radial interval with which it is measured and is divided by 2 to account for the surface that limits the accessibility of water. Finally, the number is normalized by the density of water in g cm-3. CMD simulations were performed using the classical MD code, LAMMPS.66 A 1 fs time step was used with the velocity Verlet integration, Nose/Hoover thermostat, in the canonical (constant number of molecules N, constant volume V, and constant temperature T) ensemble. For the larger scale simulations, there were 2754 atoms in the slab representing the solid and 2320 water molecules. An orthorhombic cell was used with dimensions a = 41.2587 Å, b = 43.9618 Å, and c = 54.8425 Å. The water density was fixed at 1 g cm-1. The simulations were run for 4 ns, and particle-particle-particle mesh (PPPM) was used to treat the long-range electrostatics.67 The initial configurations of the classical simulations were produced by cleaving the surface slab from a perfect crystal, adding hydroxyl groups to the surface. Water was added at 1 g cm-3 density (and was equilibrated by running a separate simulation of only SPC water). For the larger CMD calculations, the atoms in the bulk part of the surface slab were held fixed (as for the small CMD calculations), but the surface atoms were allowed to relax (see Figure 2). Comparison to X-ray Reflectivity. One of the main objectives of this Article is the testing of the classical force fields through comparison of MD simulations with results of XR measurements. From the simulations, we can compute axial densities, ni(z), where z is the distance perpendicular to the surface. XR probes the vertical interfacial electron density profile. This can be obtained by summing the individual number densities for each element obtained from simulation (laterally averaged over the surface unit cell) after scaling for each element’s atomic number, Zi, giving the weighted electron density, F(z) = ΣiZini(z). The Schlegel et al.28 results present the density profiles after taking into account the finite experimental resolution of the XR measurement. This leads to differences in the widths of the individual features in the electron density profiles as defined by XR and for those obtained by computational approaches. Because of these differences in the 2079
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
ARTICLE
Figure 4. Radial distribution function of Oh and Hw for outer (Oho-Hw) and inner (OhI-Hw) hydroxyls plotted for CMD simulations using force fields ClayFF (green), CWCA (black), LFF (red), and AIMD simulation (blue dashed).
Figure 5. Radial distribution function of Hh and Ow for outer (Hho-Ow) and inner (HhI-Ow) hydroxyls plotted for CMD simulations using force fields ClayFF (green), CWCA (black), LFF (red), and AIMD simulation (blue dashed).
two sets of profiles, the comparison of MD- and XR-based profiles in this Article will focus on the derived positions of the interfacial species obtained by the two techniques, which are quantities that can be compared directly.
followed by CWCA, with LFF being the smallest. The larger charge on the Oh will create a stronger Coulombic interaction leading to greater structuring. The lowest Hh partial charge is for LFF, with ClayFF and CWCA exhibiting similar values, yet the first peak height of the Hh/Ow RDFs follows the same trend as the Oh/Hw RDFs. However, because the Hw atoms are connected to Ow atoms, the RDFs are highly correlated, and it appears that the RDFs are dominated by the Oh charges. It should be noted that the simulations involve many waters and many competing interactions; the differing bulk charges as well as the Coulombic repulsions between oxygens in the surface and water and hydrogens in the surface and in the water may have an effect. The analysis of the different interactions for explaining the differences in force fields is therefore nontrivial, and further speculations will be avoided. It could be argued that this comparison could be biased by the starting conditions because of the short time scale and small system size. To address this, the RDFs were compared to larger simulations for nanosecond time scales described in the methods section. The first peaks of RDFs for the small and larger simulations are shown in the Supporting Information (Figures S1 and S2 and accompanying text). Overall, qualitative agreement is shown between the first peaks of the small cell and large cell simulations, adding support to the comparison between small scale simulations for evaluation of RDFs. Other preliminary simulations have been performed (results not shown), where the AIMD was performed without any constraints on center positions and comparisons were made to CMD. The center of the slab was held fixed at the AIMD configuration for CMD simulations, and a completely different starting water and hydroxyl structure were used. The same order of peak heights and distance of closest approach were observed, indicating that sampling errors will not significantly affect the result. It should be noted that, although ClayFF gave the closest agreement to AIMD for water structure, calculations performed when the whole slab was allowed to be flexible visibly showed that the bulk quartz structure was not reproduced by ClayFF or CWCA but LFF did a better job (data not shown). MD simulations of bulk quartz were performed for the three different force fields, and some results are presented in the Supporting Information. Table S0 shows LFF was able to reasonably predict the bulk quartz structure. ClayFF, on the other hand, gives a large discrepancy in both the average Si-O bond distance (underpredicts by 0.7 Å) and the average O-Si-O angle (overpredicts by 14) when compared to experimental data.64 When the RDFs were repeated for fully flexible surfaces, the underlying structure was modified causing the water structure to change. This did not occur for LFF because it was able to
’ RESULTS Comparison of CMD to AIMD. AIMD simulations were performed and used as one test of the accuracy of the three force fields; comparison with XR is considered in the next section. When considering the interaction between water and quartz, the dominating interactions are the hydrogen bonds between the water and surface hydroxyls. Figure 3 shows a snapshot of the AIMD run, indicating the main two forms of the interaction, a Oh and Hw hydrogen bond (a) and a Hh and Ow hydrogen bond (b). Radial distribution functions were calculated for the AIMD simulation to quantify how each force field performs in relation to the AIMD. As stated earlier, there are two different types of hydroxyl at different heights above the surface. It might be expected that the hydrogen bonding at these surfaces should be different. The RDFs for the individual hydroxyl types are therefore shown in Figure 4 for Oh-Hw and Figure 5 for Hh-Ow. Inner hydroxyl oxygens and hydrogens, respectively, will be labeled OhI and HhI, and outer hydroxyl oxygens and hydrogens respectively will be labeled Oho and Hho. First, peak heights in both O and H RDFs are greater for outer hydroxyl profiles than for inner hydroxyls. This is not surprising because the accessibility for water to the outer hydroxyl should be greater since it protrudes further into the solution. For three of the interactions (OhI-Hw, Hho-Ow, and HhI-Ow), the first peak is more pronounced for the AIMD than for classical MD using all three force fields. The Oho-Hw RDF for ClayFF has a larger peak than AIMD, however. For all four RDF plots, the heights of the first peaks for ClayFF are closest to the corresponding AIMD peaks, and the distance of closest approach is similar to AIMD; LFF and CWCA differ from AIMD quite significantly. In fact, LFF gives the poorest agreement with AIMD, yielding a minimal first peak for all four interactions. The inner hydroxyl for ClayFF seems to match AIMD very closely for both Oh and Hh RDFs. However, the outer hydroxyls differ slightly with the distance of closest approach and size of the peak being underrepresented for Hho-Ow and the peak shape being different for Oho-Hw. Overall, the heights of the peaks decrease in the order AIMD ≈ ClayFF > CWCA > LFF, and the distance of closest approach increases in the order AIMD ≈ ClayFF > CWCA > LFF. These trends are not surprising when the O and H partial charges are considered (Table 1). The Oh/Hw RDFs follow the same trend as the Oh charges. That is, ClayFF is the greatest,
2080
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
Figure 6. The distribution of inner silanol angles (Si-OhI-Hh) and outer silanol angles (Si-Oho-Hh); comparisons between classical force fields ClayFF (green solid - outer, green dashed - inner), CWCA (black solid - outer, black dashed - inner), LFF (red solid - outer, red dashed - inner), and AIMD simulation (blue solid - outer, blue dashed - inner).
reproduce the bulk quartz structure adequately. Therefore, the results presented are given with the proviso that when modeling quartz surfaces with ClayFF, the bulk of the crystal must be held fixed. ClayFF, however, has been parametrized for silicate/water interactions, and there is no a priori (to this study) reason for considering that ClayFF cannot reliably treat the surface structure when the underlying bulk structure is fixed at the correct configuration. One aspect of the surface structure was tested by comparison of the Si-Oh-Hh angle for the three force fields with AIMD. Figure 6 shows the distribution of these angles for the two types of hydroxyl throughout the simulation. The first noticeable feature is that the outer and inner hydroxyl angles for AIMD have different preferred angles (outer hydroxyl ∼121 and inner hydroxyl ∼114). Hydrogen bonding between Hh and Ow apparently leads to a straighter hydroxyl angle due to the attractive interactions between the Hh and the water. The inner hydroxyl angle is more acute than the outer hydroxyl, probably because the smaller degree of hydrogen bonding with water (see Figure 5) means that it is straightened less. All the classical simulation results, however, show similar peak positions with differences between outer and inner hydroxyl that are less that 2. This suggests an inherent weakness in classical simulation for describing more subtle effects. The difference between the average Si-O-H angle for different types of silanol surface hydroxyls is likely due to the environment of the particular hydroxyl. This would involve the position of hydroxyl atoms in relation to the surface and the extent of hydrogen bonding to water. For the classical simulations, however, the environment seems to be less important than for AIMD, and the angle is probably dominated by the angle bending terms. The angles for both hydroxyls for CWCA are very different from those for AIMD, ClayFF, and LFF; both peaks for CWCA are offset by 15-20. It will be shown later that this more acute angle explains structural differences in large scale simulations. The outer hydroxyl peaks for ClayFF and AIMD are at similar positions (although AIMD shows a broader peak). The inner hydroxyl, however, has an angle greater than AIMD because there is not the reduction of angle as compared to the outer hydroxyl (in fact, there is a slight increase in angle). The inner hydroxyl for LFF has a peak position similar to that of AIMD, but the outer hydroxyl is greater than that observed for AIMD. This is because LFF has similar values for outer and inner hydroxyls but AIMD does not. A simulation was performed where the CWCA Si-Oh-Hh angle potential was modified to reproduce the angle given by
ARTICLE
Figure 7. The distribution of Hh-Oh bond lengths (left) and Oh-Si (right) bond lengths obtained from CMD simulations using different classical force fields: ClayFF (green), CWCA (black), LFF (red), and AIMD simulation (blue dashed).
AIMD to test if this might affect the water structure. The hydroxyl-water RDFs, however, were not affected by such changes, indicating the central importance of intermolecular parameters for dictating water structure. The Hh-Oh bond lengths and Oh-Si bond lengths for the different force fields and AIMD were also investigated to give further insight into differences in hydroxyl structure and their effects on the water structuring. Figure 7 shows the distribution of Hh-Oh bond lengths and Oh-Si bond lengths. LFF best reproduced the AIMD results for both distances, whereas ClayFF gave the poorest agreement with a larger distance for Hh-Oh (about þ0.05 Å) and a smaller distance for Oh-Si (about -0.1 Å). CWCA gave good agreement with AIMD for Hh-Oh, but was off by about -0.04 Å for Oh-Si. It is not surprising that CWCA and ClayFF underestimated the Si-Oh bond distance, but LFF gave a comparable value to AIMD because CWCA and ClayFF also underestimate the Si-O bond distance in bulk quartz (see the Supporting Information). Despite the disagreements between ClayFF and the AIMD results, it is not expected that differences of less than 0.1 Å will make a significant difference to the water structure. There were no significant differences between bond lengths of inner and outer hydroxyls. Comparison of X-ray Reflectivity Results to Classical MD. One goal of this Article is to compare XR experimental results to molecular dynamics using different force fields to gain an understanding of how the force-field parameters (especially the intermolecular parameters) affect the structure of water at the interface. It will become evident later that the axial density profiles are generally similar for the different force fields. Details of the local water structure can still be quite different, and this may have profound effects on some of the other properties (e.g., the angular distribution of water). Figure 8 shows the axial electron density of the (1011) surface in contact with bulk water for the three force-fields as compared to the interfacial density profiles obtained previously from XR results.28 The quartz structure is held fixed to the bulk crystal structure for the second layer and all deeper layers and is the same for all force fields. The axial densities overlap perfectly in the bulk. The differences for the surface atom distributions (Si, siloxane oxygens, and hydroxyls) are due to the differences in force field because the surface atoms are allowed to relax freely (see Figure 2). Because XR is sensitive to the electron density, elements with a greater atomic number will yield a greater density. The electron density of water is seen for heights above 21 Å with a bulk value of 0.33 e-/Å. 2081
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
ARTICLE
Figure 8. Axial density plot for all atoms in the (1011) surface of Rquartz. These involve plots obtained from CMD simulations using three force fields: ClayFF (green), CWCA (black), LFF (red). Each atom is scaled by its atomic number. The XR data are taken from Schlegel et al.28
The first noteworthy feature is that the hydroxyl peaks (located between ∼14 and ∼16 Å) are distinct when comparing the XR results with all of the force fields. The experimental result has one clearly resolved hydroxyl peak (at ∼15.5 Å) corresponding to the outer hydroxyl. A second surface inner hydroxyl peak is found at ∼14.2 Å but is only partially resolved from the SiO layer below it due to the finite resolution of the data. Two peaks between 14.5 and 15 Å and between 15.2 and 16 Å for MD plots correspond to separate hydroxyl types (see Figure 1). The force fields also give peaks at different positions. There are discrepancies in the interfacial water profiles between the simulations and XR results. The first layer of water from XR is broader and higher than for the simulations showing a higher density of water for this layer. However, the XR results are closest to the LFF result, with a similar layer position and density, but with the LFF results showing a reduced density. CWCA and ClayFF, however, predict a sharper interfacial water profile, with CWCA rising higher than for XR, indicating a higher degree of order. The ClayFF result is intermediate between the LFF and CWCA. The first layer in the CWCA profile is closer to the surface, probably because the Oh are also closer. The ClayFF water peak is closer to the surface than the LFF water peak, and this is probably because of the stronger hydroxyl-water interactions. All simulations have a second water peak (located at ∼19 Å for ClayFF and CWCA and a smaller peak at ∼20 Å for LFF), although the experimental results did not show any evidence for additional ordering in the hydration structure at these heights. To gain an understanding of the total axial density shown in Figure 8, axial density profiles of the separate atom types are shown in Figure 9. Also shown are the densities of surface atoms for an AIMD simulation (here the water axial density was not included because it was not considered that it was sampled to convergence; see the Discussion) . These results show that there are distinct differences in the surface structure between the different force fields and AIMD. First, the LFF and CWCA have two roughly equal density siloxane oxygen (Osil) peaks (although for CWCA they are of different widths and heights), but the ClayFF profile for these atoms has one large and one small peak. This is due, for ClayFF, to a relaxation of the Osil, which tends to flatten this layer. This does not occur for CWCA and LFF because the Osil are distributed between two resolved heights. The AIMD result shows a broad double peak rather than two defined peaks. There are differences in the surface silicons (Si1, Si2, Si3), and this is significant because Si2 and Si3 are the ones that are connected to the Oh. Some of the separation distances
Figure 9. The axial densities of separate atom types obtained from CMD simulations using the three force fields (bottom): Hh - hydroxyl hydrogen (black), Hw - water hydrogen (red), Oh - hydroxyl oxygen (green), Obulk - bulk oxygen (blue), Ow - water oxygen (purple), Si bulk and surface silicon (orange). Also indicated are various distances (Δz1-4). The top shows the quartz surface with atom labels indicated.
Table 2. Various Distances in z Corresponding to the Axial Density Plot of Separate Atoms in Figure 9a
a
Δz1
Δz2
Δz3
Δz4
ClayFF
1.221
1.5634
0.807
0.4646
LFF
1.3303
1.5927
0.68
0.4536
CWCA
1.2893
1.5525
0.6865
0.4233
AIMD
1.302
1.626
0.781
0.458
Measured in angstroms.
between different types of silicon (Δz4) and Oh (Δz3) and between Si and Oh (Δz1 and Δz2) are indicated in Figure 9 and shown in Table 2. The distances were measured using the peak average z-position (Avz) calculated by integrating the peak from the start of the peak (pst) to the end (pe) and dividing by the integral of densities F(z) within the peak: R pe pst zFðzÞ dz ð2Þ Avz ¼ R pe pst FðzÞ dz For CWCA, the outer silicons are close together (Δz4 = 0.42 Å). They are, however, further apart for LFF (Δz4 = 0.45), ClayFF (Δz4 = 0.46), and AIMD (Δz4 = 0.46). Distances between different types of Oh (Δz3) show larger values for AIMD and ClayFF than for CWCA and LFF, and this is related to Δz4 because the Oh are connected to Si2 and Si3. The similarity of Δz3 values for AIMD and ClayFF may be important for the water structure, that is, that water “sees” the same surface topography. The distances between Si and Oh (both Δz1 and Δz2) are greater for LFF and AIMD than for ClayFF and CWCA, and this is because of differences in Si-O distances (as shown in Figure 7). It should be noted that the distance between Si2 and 2082
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
Figure 10. The relaxation in the z-direction of various surface atoms from their position in the ideal bulk crystal obtained from CMD simulation using different classical force fields: ClayFF (green 0), CWCA (]), LFF (red þ), and AIMD simulation (blue *). XR data (purple O) taken from Schlegel et al.28 are also shown. A positive value indicates a relaxation toward the solution. A negative value indicates a relaxation toward the bulk.
OhI (Δz1 = 1.2-1.3 Å) is lower than Si3 and Oho (Δz2 = 1.51.6 Å) because the first Oh does not protrude straight into the solution but at an angle. One of the most striking differences between force fields (and AIMD) is in the Hh peak. For LFF, there are two main peaks, one at the same distance as OhI and one between Oho and the solution. For ClayFF, there is a very small peak close to OhI, one centered over the Oho peak, and one between Oho and the solution. For CWCA, there is a very small peak close to OhI and a broad peak (with a shoulder) centered over the Oho peak. There is not therefore a separate peak extending toward the water. This can be related to the Si-Oh-Hh angle for CWCA, which has an angle of about 20 lower than the other force fields, causing the Hh to be closer to the surface than for ClayFF, LFF, and AIMD. The Hh peaks for AIMD appear to be most similar to ClayFF, that is, two peaks, one at the level of Oho and one between Oho and Ow. This similar profile could be another reason for the similarity in water structure between AIMD and ClayFF; that is, a similar position of Hh will again yield a similar topography. Despite the similarities, there is no small peak close to OhI for AIMD. The Hw comes closer to the surface than the Ow for all force-fields; the Hw density first becomes nonzero at the same distance from the surface as the Oho, indicating that water can locate between Oh atoms. In the work by Schlegel et al.,28 the analysis of the XR data was achieved by a flexible model of the interfacial structure that allowed for relaxed heights for surface Si and Oh (with the location of the bridging oxygen sites being defined by the Si atom positions). This was performed by an optimization starting with the perfect quartz bulk structure terminated by Oh. The positions relative to the perfect crystal of the various atoms were therefore obtained. These data are shown in Figure 10 and are compared to results from CMD (large simulations) and AIMD. In Figure 10, a negative value shows a displacement toward bulk quartz, and a positive value is a displacement toward the solution. The atom labels relate to those in Figure 9; layer 1 is the layer of interest, the interface atoms. Layers 2-4 correspond to the atoms that were fixed at their bulk positions in the simulations. The displacements for these atoms obtained by XR are small, typically less than 0.05 Å. For atoms of layer 1, the simulations give results that were distinct from that obtained by XR. The most notable differences are that OhI, Oho, and Si1 heights obtained by XR are displaced
ARTICLE
Figure 11. The axial density profile of Hw for the immediate vicinity of the interface obtained from CMD simulations using three force fields: ClayFF (green), CWCA (black), LFF (red).
Figure 12. The axial density profile of Ow for the immediate vicinity of the interface obtained from CMD simulations using three force fields: ClayFF (green), CWCA (black), LFF (red).
toward quartz but are displaced toward the solution for all simulations. Only the Si3 has a positive displacement as seen by XR, and the Si2 site shows no measurable displacement. The main trend for the simulation methods is an incremental increase as we move closer to the solution, that is, an expansion as compared to the perfect structure. This is certainly the case for AIMD and LFF, but there are deviations from this trend for CWCA and ClayFF. The largest interfacial displacements are in the AIMD, corresponding to the simulation with the highest level of theory. Possible explanations for the discrepancies between the AIMD, CMD, and XR results will be given below in the “Issues Affecting the Comparison between CMD, AIMD, and XR”. The Hw and Ow peaks differ between force fields, and closeups of these in the immediate vicinity of the interface are shown in Figure 11 (Hw) and Figure 12 (Ow). ClayFF and CWCA show a closer approach to the surface than LFF for both the Hw and the Ow. In the case of ClayFF, this probably because of the stronger interactions of water for hydroxyl atoms as seen for the RDFs in Figures 4 and 5; for CWCA, it is because the Oh are closer to the bulk crystal meaning closer approach of water. The structure of the Hw predicted by ClayFF and CWCA differs significantly from that of LFF. Specifically, LFF predicts a single first peak, while ClayFF and CWCA predict a split first peak. CWCA and ClayFF show a somewhat sharper first Ow peak than LFF. Again, this could be explained by a stronger hydroxyl-water interaction. It would be expected that the first peaks in the Hw and Ow axial density profile are related to the first peaks of the RDFs (shown above). The waters involved with hydrogen bonding with the hydroxyls should be the same waters occupying the first water layer of the axial density profiles. One question is how the first axial density peaks are partitioned into waters hydrogen bonded 2083
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
Figure 13. The contribution of Hw that are hydrogen bonded to Oho (black) and OhI (red) to the total Hw axial density profile (green). “Oho-Hw other” (dashed black) and “OhI-Hw other” (red dashed) involve the Hw that are on the same water as Hw that are hydrogen bonded to Oho and OhI, respectively.
with each type of hydroxyl. Figure 13 shows such contributions for all three force fields. This analysis is a combination of RDFs with axial density profiles because the RDFs are calculated and the axial density of any atoms within the first RDF layer is plotted. Figure 13 shows axial densities for the first Hw RDF peak. In this case, the axial density is counted if the distance is less than 2.5 Å to completely incorporate all atoms in the first RDF layer. Oho-Hw and OhI-Hw involve any Hw that is hydrogen bonding with outer and inner hydroxyls, respectively. “Oho-Hw other” and “OhI-Hw other” are the other Hw on the same molecule that were hydrogen bonded to outer and inner hydroxyls, respectively. Figure 13 provides insights into the reasons for the differences in axial density between force fields. The first noticeable feature is that the peaks for ClayFF are the largest because the first peak in the RDFs contains more waters. For ClayFF, it is obvious how the Hw axial density is made up. The Hw that hydrogen bond with inner hydroxyls (OhI-Hw) make up the density closest to the surface (at ∼15.7 Å, this peak exactly lines up). The Hw that hydrogen bond with the outer hydroxyls (Oho-Hw) make up the first shoulder at ∼16.5 Å. The other Hw within the same molecule as the hydrogen that hydrogen bonds with the outer hydroxyl (Oho-Hw other) makes up much of the second shoulder at ∼17.5 Å. For ClayFF, most of the axial density up to 17 Å can accounted for in water that is hydrogen bonded to hydroxyls. For LFF, this is not the case; there is density that cannot be accounted for. There are Hw that are close to the surface not because of surface-induced structuring but probably due to indirect affects (e.g., because the Hw are attached to structured Ow). This indicates a massive difference in behavior among force fields, which is due to the differences in partial charges of surface hydroxyls. It explains the double Hw peak for ClayFF but single one for LFF; each part of the double peak relates to specific hydrogen bonding between the water and the surface. CWCA shows a less defined double peak than does ClayFF because of the weaker hydrogen bonding. We notice that the first Hw axial density for the LFF model consists of waters that are not hydrogen bonded to hydroxyls. What waters does it contain? To address this question, a similar analysis was performed showing the axial density of waters in the second hydration shell of Oh (data not shown), revealing that these waters filled up the first Hw axial peaks. It would be expected that the Ow that are hydrogen bonded to the different types of Hh will make up different parts of the first peak of the Ow axial density; that is, a separation of the first Ow
ARTICLE
Figure 14. The contribution of Ow that are hydrogen bonded to Hho (black) and HhI (red) to the total Ow axial density profile (green).
Figure 15. A snapshot showing a water hydrogen bonded to both inner and outer hydroxyls. The angles with respect to the plane of the surface (θ1 and θ2) are indicated to explain how the water is able to hydrogen bond to both hydroxyls despite HhI being closer to the surface than Hho.
axial density peak where Ow hydrogen bonded to HhI will form the part of the peak closest to the surface and Ow hydrogen bonded to Hho will form the part of the peak furthest from the surface. Figure 14 shows the contribution of water oxygens that are hydrogen bonded to the surface to the Ow axial density for all force fields: HhO-Ow and HhIOw are for outer and inner hydroxyls, respectively. It can be seen that CWCA has the narrowest peaks and LFF the broadest peaks following the broadness of the Oh axial density peak. More importantly, for all force fields, the peaks for inner and outer hydroxyls are at the same axial density, and this is surprising considering the different contributions to the axial density of the Hw. This can be explained by considering the angle at which hydrogen bonding occurs in relation to the plane of the surface. Figure 15 shows a snapshot indicating an example where a HhO-Ow and HhI-Ow hydrogen bond occurs for the same water molecule. HhO is above the HhI as expected, yet both are able to hydrogen bond to the same Ow. Moreover, the Oh-Hh-Ow angle is approximately the same (there is an equal distribution, data not shown), and the Si-Oh-Hh angles are similar (see Figure 6). HhI is able to hydrogen bond to the same water as HhI because the inner hydroxyl is bent with respect to the surface, whereas the outer hydroxyl extends straight into the solution (θ3 ≈ 90, θ4 > 90). This means that the angle of the hydrogen bond with respect to the plane of the surface is less for the outer hydroxyl than for the inner hydroxyl (θ1 < θ2). This difference in angle means that the distance in z between Hh and Ow is less and explains Figure 14. As for the Hw axial density, the rest of the first Ow axial density peak involves waters in the second peak of the Hh-Ow RDF. 2084
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C
Figure 16. The distribution of the cosine of the angle between the water dipole vector and the normal of the surface as a function of distance from the surface. A positive value indicates the waters on average sit with their oxygens facing the surface, and a negative value indicates waters on average sit with their hydrogens facing the surface.
Water Dipole Distribution. To gain more insight into the differences between the different force fields, the distribution of the cos(θ), where θ is the angle between the water dipole vector and the surface normal vector, is shown in Figure 16. This tells us whether the water, on average, has its hydrogens facing toward the surface (negative cos(θ)) or away from the surface (positive cos(θ)). It can be seen that this distribution is different for different force fields, showing that this analysis is more sensitive to changes in force field than the axial density. It is likely that the distribution of the cos(θ) is a result of a balance between the water-water and hydroxyl-water interactions. LFF and ClayFF have a negative value for the waters closest to quartz, but CWCA is positive at the closest separation (z ≈ 16 Å). It should be noted, however, that because this analysis gives us an average over the available waters at a particular position, at close distances, the value is a result of densities lower than that of bulk water (from a density of 0.1 to ∼1 g cm-3), meaning that such differences are not so significant as differences in high density regions relating to the peaks in axial density. The start of the first axial density peaks, above bulk water densities, therefore lines up with the peaks in Æcos(θ)æ at ∼16.5 Å for CWCA and ∼17 Å for ClayFF and LFF. The Æcos(θ)æ for all force fields then undulates from positive to negative, approximately lining up with the density of Ow. This value becomes less ordered further from the surface as the influence of the surface is diminished. LFF, however, exhibits a smaller amplitude of undulation than do CWCA and ClayFF, indicating a lower influence of the surface on the water. This is probably because of the lower hydroxyl partial charges for LFF as compared to the CWCA and ClayFF.
’ DISCUSSION This work has shown the importance of testing available force fields for describing the quartz/water interface. All of the force fields evaluated have been parametrized using available experimental information, and they all predict a total axial density that is generally similar to, but not the same as, that obtained by X-ray reflectivity. For predictions of more sensitive properties such as hydrogen bonding (in comparison to AIMD result), it was shown that ClayFF was significantly closer than the other two force fields. The results reported here are a precursor for further work for understanding the interaction of water and ions at the quartzwater interface at different solution chemistries. Now that we
ARTICLE
have a basis for believing that ClayFF is better than the other two force fields at treating water structure at the quartz-water interface, we will be more confident in its potential use in the future. That is, ClayFF will be used as a starting point for the development of a force field used for studying negatively charged (deprotonated) surfaces. This leads to the question of why the treatment of water structure is so important in the properties of the quartz/water interface? And why have such great lengths been made in the present study to ensure that the best force field be chosen for describing the water structure? One of the reasons for studying the interaction of quartz with ions and water is the charging behavior of silica. It is known that silica surfaces develop a pH-dependent surface charge that depends on the choice of dissolved cation.68 Specifically, the surface charge increases with ionic radius for monovalent ions and decreases with ionic radius for divalent ions. It was postulated that the reason for this behavior is the structure modifying effects of cations on the interfacial water. Ions will have their own solvation shell whose size and structure vary with ion size and charge. Interfacial water will inevitably be affected by the ions solvation shell. Furthermore, the approach of ions themselves will be affected by the interaction between the ion solvation shell and interfacial water. Therefore, the correct treatment of water is crucial for understanding charging phenomena. If we use a force field that does not describe the interfacial water structure correctly, the interaction between the ion solvation shell and interfacial water will not be described properly, and it will be impossible to gain an understanding about the charging behavior. It was shown that the orientation of water dipoles near the surface was highly dependent on the force field. Nonlinear optical spectroscopies (e.g., sum frequency generation studies)69-74 are sensitive to the orientation of water near the quartz surface. It will be important to compare these results to simulation results from this Article and the proposed ClayFFbased force field for negatively charged surfaces. Because the Si-Oh bond length is underestimated for ClayFF, the force-field may need to be modified before moving on further. There have been many ab initio calculations involving the calculation of the reaction pathway for breaking the surface Si-O bond and hence obtaining estimates for the activation energy of dissolution.75-81 The activated complexes were found in situations where both the water Ow-Si and the Hw-Osil distances were much smaller than the corresponding weakly absorbed states. The absorbed states, involving water interacting with the surface in a specific configuration, are attractive and therefore have a negative energy. The activation energy was consequently reduced by the attractive interaction energy of the absorbed state because the absorbed state is the starting point for dissolution. This suggests that water has to be in a specific configuration, a particular orientation above the surface, for dissolution to occur. It is proposed that CMD can be used as a tool for deciding how reactive to water dissolution a quartz surface will be because MD will show the frequency with which a water molecule will approach the surface in the specific configuration identified by the ab initio calculations. This approach can then be used to evaluate the possible difference in reactivity for surfaces in different crystallographic directions, roughness (e.g., steps, etch pits),28,82 and changes in solution chemistry because all these factors will inevitably affect the water structuring. However, before MD can be used for this purpose, it must be verified that the force field is able to treat the water structure correctly. The simulations in the present study have provided valuable insight, 2085
dx.doi.org/10.1021/jp109446d |J. Phys. Chem. C 2011, 115, 2076–2088
The Journal of Physical Chemistry C showing the viability of ClayFF for treating hydrogen bonding correctly and therefore for evaluating accurate frequencies of water approach. Issues Affecting the Comparison between CMD, AIMD, and XR. We should exercise some caution in reaching strong conclusions based on comparing the classical MD results to AIMD because sampling and reaching equilibrium is a major issue for AIMD. Nevertheless, the RDFs show profiles similar to CMD over system sizes and times much larger than AIMD. The consistency at least shows that the lack of sampling has not really affected the comparative sizes of the first RDF peaks. RDFs are a good analysis for comparison because RDF space is better sampled than axial density profile space since RDFs take into account many water-hydroxyl interactions (i.e., 16 in the case of these). The axial density profiles do not reveal as much structural detail as the RDFs because the hydroxyl-water distributions are a direct effect of the respective interactions, but the axial density profile is a result of a combination of many interactions. The angular distribution of water dipoles is even more sensitive to sampling issues. As far as the axial density of surface atoms is concerned, the time needed for sampling is lower than for water because the temporal variation of the structure is substantially reduced. Therefore, the axial density of surface atoms can be evaluated for AIMD simulations. A separate issue is the accuracy of the AIMD simulations. For example, DFT lacks dispersion interactions. The water structure can be expected to be dominated by electrostatic and short-range repulsive interactions, which DFT handles quite well. The RDFs computed from AIMD for bulk water are in good agreement with neutron scattering, except for the peaks being slightly narrower.83 Thus, we expect the AIMD results for quartz-water reported here to be similarly accurate, and if anything slightly overstructured. To be more confident about AIMD results, different starting configurations and longer simulations will be used in future work to increase the sampling. It may also be helpful to increase the size of the system to account for boundary effects. Because the charges are key to the interfacial water and surface structure, a discussion of the various methods of charge determination should be made (see Methodology for further details). ClayFF appears to have produced the most accurate charges, probably because the most consistent method of charge determination was employed. CWCA was put together from other force fields, which may not be the best way of producing force fields. LFF developers fitted the charges using ab initio clusters with isolated gas-phase molecules (single waters), which could result in the condensed phase properties of water being neglected. LFF, however, gave the best bulk quartz structure (see the Supporting Information), and it is not surprising that the structure produced from ClayFF was less reliable because there are no special bonding or angle terms within the bulk crystal for this model. Much effort was put into parametrizing LFF to reproduce the bulk structure, that is, fitting to ab initio geometries and vibrations. It should be noted that (currently) for the sake of understanding the quartz/water interface, there is little requirement to explicitly treat the vibrations in the bulk part of the slab, holding the bulk atoms fixed is adequate, and ClayFF is probably the best force field to use. However, for future studies where there is a focus on the bulk quartz structure, ClayFF would not be suited, and LFF may be a better choice. In this case, LFF would have to be modified to better describe the water/hydroxyl hydrogen bonding.
ARTICLE
In the results of the Schlegel et al. study,28 the measured vertical relaxation was small but nonzero (
