Article pubs.acs.org/JPCC
Is the Calcite−Water Interface Understood? Direct Comparisons of Molecular Dynamics Simulations with Specular X‑ray Reflectivity Data Paul Fenter,*,† Sebastien Kerisit,‡ Paolo Raiteri,§ and Julian D. Gale§ †
Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne Illinois 60439, United States Chemical and Materials Sciences Division, Pacific Northwest National Laboratory, Richland Washington 99352, United States § Nanochemistry Research Institute, Department of Chemistry, GPO Box U1987, Curtin University, Perth, WA 6845, Australia ‡
S Supporting Information *
ABSTRACT: New insights into the understanding of calcite− water interface structure are obtained through direct comparisons of multiple classical molecular dynamics (MD) simulations with high-resolution specular X-ray reflectivity (XR) data. This set of comparisons, with four different stateof-the-art force fields (including two nonpolarizable, one polarizable, and one reactive force field), reveal new insights into the absolute accuracy of the simulated structures and the uniqueness of the XR-derived structural results. These four simulations, though qualitatively similar, have visibly distinct interfacial structures and are distinguished through a quantitative comparison of the XR signals calculated from these simulations with experimental XR data. The results demonstrate that the simulated calcite−water interface structures, taken as a whole, are not consistent with the XR data (i.e., within the precision and accuracy of the XR data). This disagreement is largely due to the simulated calcite interfacial structure. The simulated interfacial water profiles show a higher level of consistency with the XR data, but with substantially different levels of agreement, with the rigid-ion model (RIM) simulations showing semiquantitative agreement. Further comparisons of the structural parameters that describe the interfacial structure (derived from both the MD simulations and the XR data) provide further insight into the sources of differences between these two approaches. Using the new insights from the RIM simulations, new structures of the calcite−water interface consistent with both the experimental data and the simulation are identified and compared to recent results.
1. INTRODUCTION Robust knowledge of the structure and interactions of matter derives from the interplay of precise experimental observations (e.g., by X-ray crystallography) and computational studies (e.g., high-level quantum mechanics). The accuracy of our understanding is tested, ultimately, by the consistency of observations using different and complementary approaches. The absolute accuracy of the derived equilibrium structure (corresponding to their free energy minima) is important as it forms the basis for understanding excitations and reaction kinetics in these systems (which depend, respectively, on the curvature of the free energy about the minimum location, and the barriers between distinct minima). Although this approach has been well-developed for understanding bulk crystalline structures, there are a number of challenges with respect to establishing the accuracy of our understanding of liquid−solid interfaces: (1) The presence of an interface leads to strong gradients in the molecular-scale structure and composition. This can lead to significant changes in the local atomic structure and coordination relative to those that are normally present in homogeneous media, leading perhaps to the creation of unique local environments. © 2013 American Chemical Society
(2) Solid−liquid interfaces inherently are extended structures. Substrate surface displacements can extend many unit cells below the solid surface, and the interfacial solution structure may be perturbed many nanometers into the solution.1 This is inherently challenging to both experimental and computational analyses of interfacial structure. (3) There are few experimental probes that can measure the interfacial structure with the accuracy and precision that can be achieved for bulk structures (e.g., X-ray crystallography). X-ray-based probes (including various X-ray scattering and spectroscopic probes), offer numerous approaches to probing the structure and bonding of interfaces through direct observations in the environment of interest. X-ray reflectivity (XR),2,3 which is the direct analogue of X-ray crystallography for interfaces, offers perhaps the best opportunity to achieve both accuracy and precision in Received: November 5, 2012 Revised: January 16, 2013 Published: February 13, 2013 5028
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
defining the structural chemistry of liquid−solid interfaces. Classical molecular dynamics simulations based on interatomic potentials (also referred to as force fields) are attractive as a means to assess our understanding of such interfaces because of their ability to simulate systems having a large number of atoms over long simulation periods, with statistically significant results that can be compared directly to experimental observations.4−6 Yet a general challenge for performing such simulations is to evaluate the uncertainties and absolute accuracy of the results. In particular, the relationship between both the choice of functional form and parametrization of the force field and the accuracy can be difficult to determine. An additional challenge lies in comparing computational and experimental results directly. In part, this is because experimentally derived structures are often obtained with some assumptions (either explicit or implicit) that may bias the interpretation of the results (Figure 1). For instance, XR data
Figure 2. Schematic of the calcite−water interface. The calcite lattice is shown as red, black, and blue spheres (for oxygen, carbon, and calcium, respectively), and the water profile near the interface is shown as large blue circles whose size indicates the rms width of each water layer. The laterally averaged density profile is shown schematically on the right. Also shown are the structural parameters that are used to describe the interfacial structure, including the vertical displacements of the Ca (Δz, with a similar parameter for the CO3 ion), and the tilt and twist of the carbonate groups (θCO3 and ϕCO3, respectively). These parameters, shown in the top layer, are used to describe the structure in the top six layers.
the calcite−water interface performed here for several different descriptions of the interatomic forces, taken in total, do not provide a quantitative description of the XR data. Comparisons that use portions of the simulated structure (e.g., only the calcite or water portion of the simulation cell) show that this disagreement derives largely from the inability of the different force fields to reproduce the measured calcite surface structure to within the precision and sensitivity of the XR data. Direct comparisons of the measured and simulated interfacial water profiles reveal variable agreement for the different interatomic potentials, with semiquantitative agreement in some cases. The variations in the quality of agreement for the different simulations correlate with visible differences in the simulated interfacial water structures. More generally, these results demonstrate an approach to quantify each interatomic potential’s ability to reproduce the experimentally observed XR data (and by implication, the actual interfacial structure). The differences in the simulated and measured water structures allows for the accuracy of the XR results to be assessed, leading to the identification of revised structures that reproduce the observed XR data.
Figure 1. Conceptual comparison of XR and MD results. Black arrow indicates traditional comparison of XR structure (derived from parameter-dependent least-squares fitting of XR data) and timeaveraged simulated structure. Red dashed arrow indicates the present approach where the MD simulation results are compared directly to the XR data without any intervening interpretation.
are normally analyzed by comparing the data to a structural model described by a set of parameters, and the “best-fit model” is one in which the parameters give the best agreement.2,3 However, such approaches can be blind to whether another model might lead to better agreement with the data. There is always the possibility that a particular structure might disagree with the structure obtained from model-dependent fits, even if that structure were fully consistent with the experimental data (e.g., if there were multiple structures consistent with the data). One way to address this issue is through the use “modelindependent” methods to interpret the experimental data.7−9 Such approaches are powerful but are not necessarily precise, because they often are susceptible to the noise and systematic errors in real data. Instead, they are best-suited to providing a conceptual picture of the structure that is being optimized, so that an appropriately parametrized model can be built. Here, we demonstrate the utility of a new approach to evaluate the accuracy of simulated structures to reproduce interfacial properties, quantitatively and directly, without any assumptions or parameter-dependent analysis of the interfacial structure. This comparison is performed for the calcite−water interface (Figure 2). The structure of the calcite−water interface has been well-studied with both computational10−18 and experimental probes.19−25 The method that we use here provides a new route to test our understanding of this wellstudied system. The results show that the MD simulations of
2. EXPERIMENTAL AND COMPUTATIONAL PROCEDURES 2.1. XR-MD Comparisons. The X-ray reflectivity signal, R(Q), is the fraction of the incident beam that is specularly reflected from an interface as a function of vertical momentum transfer, Q = 4π/λ sin(θ), where λ is the X-ray wavelength and θ is the angle of incidence of the X-ray beam with respect to the surface plane. The XR signal has a typical shape known as a crystal truncation rod (CTR).26 Direct comparisons of simulated structures with experimental XR data can be obtained using a recently developed formalism.2,3 This approach makes use of the simple and direct relationship between the interfacial 5029
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
structure, as expressed as a laterally averaged electron density, ρ(z), and the specular reflectivity signal, R(Q): R(Q ) = (4πre/AUCQ )2 |
∫ ρ(z) exp(iQz) dz|2
= (4πre/AUCQ )2 |F(Q )|2
Here, AUC is the surface unit cell area, re = 2.818 × 10−5 Å is the classical electron radius, and F(Q) is the interfacial structure factor. The structure factor can be separated into distinct components: F(Q ) = Fcalcite bulk + Fcalcite int + Fwater int + Fwater bulk
These correspond to the structure factors of the bulk calcite lattice, the distorted interfacial calcite lattice, the layered interfacial water, and the bulk water. The bulk structure factor is in the form of a crystal truncation rod.26 The MD simulated profile can be used as the basis for calculating the contribution of the calcite interface and interfacial water to the measured XR signals,27 whereas the contributions to the structure factor from the bulk media are known. This process allows the MD simulation to be compared directly to the XR data without any adjustable parameters, at least in principle. The MD-simulated profile is in the form of a nuclear density distribution of each atom in the system, ρMD, whereas XR probes the electron density profile, ρe. These are related by the convolution of the nuclear density profile with the spatial distribution of the electron cloud surrounding each nucleus, ρatom; that is, ρe = ρatom ⊗ ρMD. The structure factor is the Fourier transform of the electron density profile which, by the convolution theorem, is the product of the Fourier transforms (FT) of ρMD and ρatom. Therefore, the interfacial structure factor can be calculated directly from the MD simulation by
Figure 3. (a) Simulated RIM-1 profiles (red line) embedded between a bulk calcite lattice and bulk water media (black line). (b) Calcite portion of the same RIM-1 simulation embedded between a bulk calcite and a parameter-dependent- fit to the interfacial water profile. (c) Water portion of the RIM-1 simulation embedded between a bulk calcite lattice, six layers of parameter-dependent fit of the near surface calcite structure and bulk water.
so that the interface between the MD-simulated and bulk water densities becomes invisible to the X-rays. The height of the MD simulated profile is adjusted initially so that the lowest calcite layer (here at z ∼ −15 Å) is one lattice spacing above the first bulk calcite layer (black line on left-hand side of Figure 3a). A fitting parameter is used to optimize the relative location of the MD profile with respect to the calculated structure because this is not known a priori. Portions of the simulation can be embedded, as shown for the MD simulation of the calcite lattice (Figure 3b) and the interfacial water portion of the lattice (Figure 3c). In these cases, the nonsimulated portion of the interfacial structure can be calculated on the basis of a parametrized model of the interfacial structure. The accuracy of the numerical Fourier transformation, due to the discrete “binning” of the simulated profiles, is determined by a combination of the bin size used in the simulation, the spatial resolution of the data, π/Qmax (where Qmax is the maximum momentum transfer that is measured), and the intrinsic width of the spatial distribution of any element in the structure.27 On the basis of these considerations, the comparison should have an absolute error in the calculated structure factors of 10 μm). Using the MD simulation, alone, to calculate the XR signal therefore leads to substantial extrinsic errors, including oscillations in the reflected intensity (due to the finite size of the MD simulation cell), and broadening of the substrate diffraction peaks due to Debye−Scherrer broadening associated with the finite crystal thickness.27 To eliminate these artifacts, we select a portion of the MD simulation to account for the interfacial structure (i.e., the portion of the calcite and water profiles that are not “bulk-like”) and embed this portion into the ideal structures of the respective bulk media. This is illustrated in Figure 3 using the data from the RIM-1 simulation (see detailed description of this simulation, below). Briefly, a portion of the MD simulation is selected by multiplying the data with a window function: W = [1 + erf(zmin,σmax)] × [1 − erf(zmax,σmax)]/4, where erf(zmin,σmax) is the error function located at zmin and having width σmax. We choose zmin to be in the zero density region between the calcium carbonate layers with a narrow width, whereas zmax is located in the water media with σmax = 0.4 Å. A complementary error function with the same width, σmax = 0.4 Å, corresponding to the bulk water media (black line on right-hand side of Figure 3a) is included 5030
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
between the data and the calculation. These qualities of fit can be compared to the values achieved with the traditional modeldependent fitting to directly and quantitatively evaluate the relative ability of the simulated structures to explain the XR data. This comparative approach opens up a new window through which we can evaluate the accuracy and precision of the experimental and computational approaches. 2.2. Experimental Details. The details of the XR measurements have been described elsewhere.25 Briefly, specular XR measurements were performed at sector 6-ID of the Advanced Photon Source using an incident photon energy of 16.0 keV, corresponding to an X-ray wavelength of λ = 0.7748 Å. The specular reflectivity signal and background are measured simultaneously by an X-ray CCD detector as a function of the momentum transfer, Q. The calcite unit cell is defined by a vertical layer spacing, c = 3.035 Å, and a surface unit cell area of AUC = 20.198 Å2. The calcite (104) surface was cleaved with a razor blade, quickly transferred (within seconds) to a calcite-saturated aqueous solution, and mounted in the thin-film cell for X-ray characterization. The XR measurements were performed in a “thin-film” cell3 with a beam cross-section of 0.04 mm × 0.5 mm. The data can be understood by comparing the measured reflectivity with that calculated for an atomistic model of the calcite surface. The model is parametrized with rotations of the individual carbonate groups (including a tilt angle, θCO3, and a twist angle, ϕCO3) along with independent vertical displacements of Ca ions (e.g., ΔCa and ΔCO3) and CO3 groups. The structure is allowed to distort within the top six CaCO3 layers. We include the presence of specifically adsorbed water layers in the present analysis. These are described by their position, zH2O_i, occupation, OH2O_i, and rms width, uH2O_i, for each layer, i, whose values are unconstrained. A multilayer water model is used to describe the weakly modulated fluid water density with a layer spacing, dwater, and a parameter, ubar, that specifies the increase in the rms width of water layers at increasing distance from the surface, un = [u02 + (n − 2)ubar2]1/2. The occupation factor in the multilayer water model, OH2O, is a calculated quantity chosen to fix the bulk water density to its known value, based on the water molecular volume, Vwater = 30.3 Å3, the surface unit cell area, AUC, and dwater through the relation: OH2O = AUCdwater/Vwater. 2.3. Computational Details. Atomistic simulations have been used extensively to study the bulk and surface structures of calcite, both in vacuo and in the presence of solvents or adsorbates, as well as the other crystalline polymorphs of calcium carbonate (for example, refs 10, 13, 16, 18, and 29−42). Correspondingly, there have been many different interatomic force fields created as the number of properties included in the fitting process has been increased or the objective of the simulation model changed.13,29,31,34,35,37,42,43 Initially, the focus was primarily on reproducing the structure, elastic, dielectric, and vibrational properties of calcite.31,34,37,43 However, the best current models also aim to describe the relative thermodynamics of polymorphism.13,42 Similarly, numerous models have been proposed to describe the interaction of calcium carbonate with water.12,14,17,32,44 It has recently become apparent that careful consideration of the free energies of solvation of the ions and the solubility of calcium carbonate is important in determining the strength of such interactions, especially for studies of crystal growth.13,42
Aside from the details of how the force field for calcium carbonate and its interaction with water is parametrized, there is also the question of the functional form. In essence, there are three broad types of interatomic potential models being used to describe the aqueous calcium carbonate interface. The simplest of these is the so-called rigid ion model (RIM) in which all atoms are point particles with both Coulombic and short-range interactions between them. Here, it is common to describe the carbonate anion as a molecular species such that the intramolecular interactions are typically harmonic restraints while the Coulomb interactions are excluded, though this is not universally true. Formal charges of +2 and −2 are applied to the calcium cation and carbonate anion, respectively, with a set of partial charges assigned to individual atoms within the latter. The second type of force field includes atomic polarizability through the use of the shell model (SM) to more accurately describe the electrostatic response properties. Although extensively used for lattice dynamics studies, the use of a shell model increases the cost of molecular dynamics substantially through the need either to adiabatically relax the shells to the ground state for each configuration or to reduce the time step to integrate the shell positions having assigned them a small, fictitious mass. The third type of model is one that is based on a bond order potential, typically with charge equilibration on the fly to determine the electrostatics, and thereby allows chemical reactivity to be included. Here the ReaxFF force field of van Duin and co-workers is the most widely employed model for aqueous species.45,46 All of the three types of force field described above represent different compromises between computational expediency and capturing the underlying physical interactions. Hence, they all have the potential to yield distinct results for the aqueous calcite interface. In the present work, we examine the results of representative models of each type. The details of the molecular-dynamics simulations have been described previously,10,13,14,42 whereas the key attributes of each of the force fields are summarized below to highlight the specific models and parametrizations used. 2.3.1. RIM-1. For calcite, model RIM-1 uses a rigid-ion model version of the calcite model of Pavese et al.31 whereby the shells are removed and the oxygen charge is the sum of the core and shell charges of the oxygen ion in the original Pavese model. All other parameters remain the same. Although this model gives an incorrect structure for aragonite, calcite is correctly reproduced, albeit with hexagonal lattice parameters that are ∼3% too short in the ab plane, whereas the c axis is correspondingly too long. For water, model RIM-1 uses the extended simple point charge model (SPC/E).47 The carbonate ion is fully flexible in the Pavese model whereas the water geometry is kept fixed in the SPC/E model. For the calcite− water interactions, RIM-1 uses three potentials: the calcium− water oxygen potential of de Leeuw and Parker;32 the carbon− water oxygen potential of Kerisit et al.;12 and the SPC/E water oxygen−water oxygen potential for the carbonate oxygen− water oxygen interactions. The calcite slab had a surface area of 23.98 × 24.12 Å2 and a thickness of 12 atomic layers for a total of 360 CaCO3 units. The calcite slab was placed in contact with 630 water molecules. Data were collected from a 1 ns NVT simulation at 300 K with a time step of 1 fs. The MD simulation was performed using the computer code DL_POLY48 with the same simulation parameters as described in reference 10. 5031
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
Figure 4. Distribution of oxygen atoms belonging to water during the MD simulations of the calcite−water interface for (a) RIM-1, (b) RIM-2, (c) SM, and (d) ReaxFF. Here the positions of the oxygen atoms of water (hydrogen atom positions are not shown, reflecting the sensitivity of the XR measurements) are cumulative over 500 frames sampled at regular intervals throughout the simulation. The spheres representing the oxygen positions of water are colored according to the different components to the water structure as a guide to the eye (first layer, yellow; second layer water, orange; third and higher layers, green), as indicated in (a). The averaged positions of the atoms in the calcite lattice are shown with Ca (green), C (cyan), and oxygen (red).
2.3.2. RIM-2. During the last two years a new set of rigid ion models for aqueous calcium carbonate systems has evolved that are derived to capture the thermodynamics of solvation, as well as calcium carbonate polymorphism, at least to the extent that is possible within the limits of the functional form. In the original model,13 water was treated as a rigid entity using the TIP4P-Ew potential.49 This was later revised42 to include fully flexible water, based on the SPC/Fw model,50 and then finally higher order coupling terms were introduced within the carbonate group to capture bond-angle coupling effects, thereby leading to accurate reproduction of the experimental vibrational spectrum.51 In the present study, we employ this final fully flexible rigid ion model to simulate the calcite−water interface. For comparison against RIM-1, this potential shows the same qualitative trend for the calcite structure (i.e., the a and b cell parameters are underestimated whereas c is too large) though the errors are now less than ∼1%. All simulations were run using the LAMMPS code.52 The calcite slab surface area was 62.44 × 76.59 Å2 and 12 molecular layers thick, for a total of 2,880 CaCO3 formula units. The system was placed in contact with 5,862 water molecules and the c axis was equilibrated for 1 ns with an NPT simulation at ambient pressure. The average value of the c axis was then used in a 6 ns long NVT simulation for the data collection. The temperature was set at 300 K and controlled with a Nosé−Hoover chain thermostat of length 5 and a relaxation time of 0.1 ps. The time step was always kept
to 1 fs. The trajectory of the atoms was saved every 100 fs for analysis. To verify that there were no finite size effects due to the slab thickness, we repeated the calculations using a calcite slab with 22 molecular layers. The surface water was not affected by the increased calcite layer thickness and only minor changes were observed in the slab itself. 2.3.3. SM. In the SM model, the model of Pavese et al.31 is used to simulate calcite and the SWM4-NPD model of Lamoureux et al.53 is employed for water. Both models use a shell model to simulate the ionic polarizability of the oxygen ions. As for SPC/E, the water geometry is kept fixed in the SWM4-NPD model, except for the position of the oxygen shell. For the calcite−water interactions, the same potentials used in the RIM-1 model for the calcium−water oxygen and carbon− water oxygen interactions are also employed in this model. In addition, the potentials of de Leeuw and Parker32 and Kerisit and Parker10 are used to account for the carbonate oxygen− water oxygen and carbonate oxygen−water hydrogen interactions. The size of the simulation cell was the same as that used for RIM-1. The MD simulation was carried out as for the RIM-1 model, except for the time step, which was reduced to 0.2 fs. The oxygen shells were given a finite mass of 0.2 au, with the corresponding amount subtracted from the mass of the core. 2.3.4. ReaxFF. Because of the interest in being able to simulate solution speciation as part of the crystal growth 5032
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
Figure 5. (a) Time averaged density profiles obtained by MD simulations (RIM-1, RIM-2, SM, and ReaxFF; see text for description of the details for each simulation) compared with density profiles obtained from the parameter-dependent analysis of XR data (XR). The origin in height is chosen such that the top layer Ca position is at z = 0. (b) Detail of the top two calcite layers for each simulation. (c) Detail of the interfacial water structure within 12 Å of the calcite surface.
process, there has been a need to develop reactive force fields for calcium carbonate in contact with water. This avoids the requirement to make assumptions as to the protonation state of species because proton transfer should be able to occur spontaneously during the dynamical simulation. To date, the most widely applied reactive model is that of ReaxFF and this has been adapted to the simulation of calcium carbonate within the approximation of calcium being a formally charged cation.14 All reactive force field simulations have been performed in LAMMPS with in-house modifications to allow for the constraint on the calcium charge. The simulations were performed with the same system used for the RIM-2 model using a 0.5 fs time step and analogously the temperature was controlled at 300 K with a Nosé−Hoover chain thermostat of length 5 and 0.1 ps relaxation time. Due to the higher complexity of the interactions the computational cost was such that the overall simulation time was limited to 1 ns. The trajectory was saved every 100 fs for analysis.
layer of Ca ions located at the height, z = 0. All four sets of results have the general characteristics of the calcite−water interface that have been observed previously, through both computational and experimental results. These include (1) the bulk calcite structure exhibiting a triple layer structure within each layer, with coplanar C, Ca, and O located on the lattice plane and two additional oxygen atoms located above and below the Ca−C−O plane, (2) the distortion of the calcite structure in the near-surface layers as seen in the profiles through the reduction in the peak density of the central Ca− C−O layer due to slight displacements of individual atomic species, (3) an interfacial water structure that exhibits one to two distinct water layers followed by more weakly modulated water layers, and (4) the electron density that tends to the bulk water value of 0.33 e−/Å3 at heights sufficiently high above the interface so that the influence of the interface is negligible (interfacial water profiles are known to transition to a bulk-like water without any residual layering within ∼20 Å of the mineral surface1). For this system, we choose the upper window function to be at a height of ∼13 Å from the interface, with a relatively broad width so that this interface is effectively invisible to the XR measurement. From the general trends that appear in these results, it is apparent that there is good consensus on the qualitative aspects of the calcite−water interface. A closer view of the structures shows, however, that there are many fine-scale details that distinguish these results. For example, the interfacial density within the top two calcium carbonate layers (Figure 5b) shows that there are visible differences in the second layer Ca−C−O layer position, suggesting differences in the separation between the first and second layers. Also apparent, from the vertical oxygen atom positions above and below the Ca−C−O layer, is that there are significant differences in the carbonate ion orientation within the top two layers.
3. RESULTS 3.1. Comparison of MD Simulated Structures. The time averaged equilibrium structures of the calcite−water interface (in the form of the laterally averaged electron density as a function of height, z) were obtained from the molecular dynamics simulations for four different choices of interaction parameters (referred to as “RIM-1”, “RIM-2”, “SM”, and “ReaxFF”). These simulations are shown in Figure 4, revealing qualitatively similar, but yet distinct, water profiles at the calcite−water interface showing two distinct layers of water organized at the interface. Differences in the water structures of these simulations include the degree of separation between the first two water layers and their lateral distribution widths. These same simulations are shown as laterally averaged vertical electron density profiles, along with a comparable interfacial density profile obtained from a parameter-dependent analysis of XR data, in Figure 5a. These results are presented with the top 5033
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
Figure 6. (a) XR reflectivity data and (b) normalized X-ray reflectivity data of the calcite−water interface (blue circles) are shown with the calculated reflectivity signals derived from the RIM-1 simulation. (c) Residual deviations between the measured and calculated reflectivity data. Calculations include the best fit structure derived from the XR data (red lines), the embedded RIM-1 simulation (illustrated in Figure 3) including the full RIM-1 simulation (green line), the calcite portion of the RIM-1 simulation (blue line), and the water portion of the XR simulation (black line).
(red and black lines, respectively, Figure 3a) so that the calculated XR signal is directly comparable to that measured for the calcite−water interface. The reflectivity signal from the full RIM-1 structure, embedded in the semi-infinite calcite and water media, is shown in Figure 6a (green line, labeled “Full”). The same calculations and data are shown after dividing out the rapid variation in the XR signal due to the CTR shape (Figure 6b) so that the differences between the calculated and measured signals are more evident. These calculations include the numerically calculated contributions from the MDsimulated profile, along with the closed form calculation from the bulk calcite and water components. For comparison, the best-fit XR result is also shown (red line). Here, we use two fitting parameters to optimize the relative height of the MD simulated structure with respect to the bulk calcite lattice (to eliminate any systematic errors associated with the placement of the MD simulation) and an overall intensity scale factor. We did not include any additional contribution from a surface roughness factor26 because the previous parameter-dependent fitting results showed that there was negligible contribution from this in the measured reflectivity signals. The agreement between the calculated reflectivity for this structure and the XR data is visibly poor (green lines, Figure 6). Note the substantial oscillations (e.g., near Q = 5 Å−1) in the calculated reflectivity, with a period that is significantly smaller than that seen in the data. This is highlighted by the substantial oscillations in the Qdependent variation in the fractional residuals between the calculated and measured reflectivity signals (i.e., residual(Qi) = (Rdata(Qi) − Rcalc(Qi))/Rcalc(Qi)) (Figure 6c), with residual values exceeding 2 (corresponding to 200% fractional error, which can be compared to the typical experimental uncertainty of ∼2%). The quality of fits for this structure are χ2 = 176.5 and R = 0.293. These values can be compared to those obtained for the “best-fit” structural model (χ2 = 1.67, R = 0.033).25 That is,
A more dramatic difference is seen in comparison of the derived water structures above the calcite surface (Figure 5c). The first water layer (i.e., corresponding to the waters coordinating with surface calcium ions) is found at heights ranging from 2 Å to 2.5 Å, and the second water layer (i.e., that coordinating with surface carbonate ions) is found at heights ranging from ∼3 to 3.4 Å. A third distinct water layer is found in the simulations at heights of 4.3 to 5.3 Å, whereas the corresponding layer for the XR-derived structure is at a height of 6.7 Å. Although the MD simulations show as many as three distinct water layers, the XR results found one distinct water layer with a weakly modulated water density beyond the first layer. It is challenging to make truly meaningful comparisons of these various results to the derived XR results, in part due to the complexity of these interfacial structures that contain both crystalline and liquid components. In particular, there is no reason that a specific derived structural parameter (e.g., the first water layer height) obtained from fitting the XR data will be accurate if any other aspect of the parametrized model (i.e., the calcite interfacial displacements) does not reproduce the actual structure. In this sense, it is useful to compare the results against a common benchmark using a single well-defined metric. For this purpose, we make direct comparisons between calculated specular XR signals derived directly from the MD simulations to the experimental data using the χ2 and R-factor quality of fits. We can also compare the optimized structural parameters to comparable quantities derived from the MD simulations to gain insights into these qualities of fit. These two approaches will be discussed separately. 3.2. Direct Comparisons of MD Simulations to XR Data. 3.2.1. Comparison with RIM-1. The RIM-1 results are used to illustrate the comparison process in detail. We begin by embedding the full MD-simulated profile into the bulk media 5034
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
the RIM-1 simulation has a χ2 quality of fit that is >100 times worse than the best fit model, demonstrating that there are serious quantitative discrepancies between the RIM-1 model (as a whole) and the XR data based on the quality of fit. See Table 1 for a listing of quality of fits for all of the simulated structures.
3.2.2. Comparison of XR Data with RIM-2, SM, and ReaxFF. The same approach as above has been used to evaluate the other three simulations of the calcite−water interface, to test the generality of the observations and the sensitivity of the XR data to the differences in the simulated structures. The calculated reflectivity profiles and residual plots for each of the three additional simulations (RIM-2, SM, and ReaxFF) are shown in the Supporting Information (Figures S1−S3), and the resulting qualities of fit are summarized in Table 1. We begin with a discussion of the relative qualities of fits for the “full” simulated structures and focus specifically on the χ2 quality of fit because that measures the agreement with respect to the experimental uncertainties (and is the guide for the leastsquares structural optimizations). The quality of fit for fullRIM-2 (χ2 = 163) simulation is slightly better than that of RIM1 (χ2 = 176.5). The SM simulation has the best agreement for the full simulation (χ2 = 118). A surprising aspect of these comparisons is that the most sophisticated force field shows the poorest agreement with the data: ReaxFF shows a quality of fit (χ2 = 504) that is more than twice as large as for RIM-1. The poor agreement between ReaxFF and the data is not entirely surprising given that the observed density profile for the ReaxFF simulation shows a top calcite surface layer structure that is the most highly distorted of the simulated structures and that is the least similar to that derived from the traditional XR analysis. Furthermore, the ReaxFF simulation yields a qualitatively distinct structure for the first hydration layer above the surface. Here the water dipole is strongly oriented so as to be nearly parallel with the surface normal, with oxygen directed toward the calcium ions of the calcite surface. All other simulations give a first water layer in which the dipole is much closer to being parallel to the surface plane instead. Although the SM simulation has the “best” overall agreement, which is substantially better than either RIM-1 or RIM-2, this level of agreement is nevertheless poor in any quantitative sense (a value of χ2 ∼ 1 would indicate quantitative agreement). It cannot be concluded, therefore, that SM “better” represents the calcite−water interface structure, because all of these simulations disagree with the experimental data when taken as a whole. The “calcite-only” comparisons show that the RIM-2 simulation has the best quality of fit (χ2 = 33.63), which is nearly twice as good as that of the RIM-1 simulation (χ2 = 57.7). The SM simulation has a slightly better quality of fit (χ2 = 49.7) than RIM-1 but is inferior to RIM-2 in this regard. Again, we find that the agreement for ReaxFF is substantially worse than the others (χ2 = 372). These results suggest that the SM and RIM-2 simulations of the calcite portion of the interface agree best with the XR data. The level of agreement is, again, quantitatively poor with respect to the uncertainties of the experimental data. The disagreement with ReaxFF suggests that the structural details of the calcite portion of the interface are not well reproduced by this simulation. Comparison of only the interfacial water portion of the simulations (“water-only”) shows that the best agreement is obtained for RIM-1 (χ2 = 3.27). A somewhat worse agreement is found for RIM-2 (χ2 = 6.92). This is nominally surprising because the actual water layer heights for RIM-2 match that derived from the XR data reasonably well, as will be discussed in more detail below. The water structures for SM and ReaxFF have substantially worse quality of fits (χ2 = 19.36 and 20.24, respectively). That is, without any reference to parametrized models of the interfacial water structure, we find that RIM-1
Table 1. Direct Comparisons of the Quality of Agreement between the XR Data and Each MD Simulation (Indicated by χ2/R-Factor Qualities of Fit for Each Comparison)a quality of fits χ2/R-factor
experiment XR simulation model
full
RIM-1 RIM-2 SM ReaxFF
176.5/0.29 163.4/0.33 118.0/0.26 504.0/0.54
1.67/0.033 calcite only
water only
57.7/0.159 33.6/0.139 49.7/0.154 372.8/0.405
3.3/0.043 6.9/0.073 17.9/0.097 20.2/0.089
a
Each comparison is made for three scenarios illustrated in Figure 3: the full simulation embedded in the bulk calcite and water phases (“full”), the calcite portion of the simulation along with a parameterized water profile (“calcite-only”), and the water portion of the simulation along with a parameterized calcite profile (“water only”). The χ2 and R-factor for each fit are listed in each box. The fits for RIM-1 to the XR data are shown in Figure 6, whereras those for RIM-2, SM, and ReaxFF are shown in Figures S1−S3, respectively, in the Supporting Information. Also included are the results from the parameter-dependent fits to the XR data.
This lack of agreement does not provide insight into the source of these discrepancies. These insights can be identified by using a hybrid approach that evaluates which aspects of the model, if any, are in agreement with the data. We start by using only the calcite portion of the MD-simulated structure, which is combined with a parameter-based fit of the interfacial water structure. In this approach, we use the same parametrization that was used to optimize the water profile for the previously obtained best-fit structure (i.e., a layered water model for the interfacial water layer; the associated structural parameters are shown schematically in Figure 2). The calculated reflectivity for the MD-simulated “calcite-only” structure (from RIM-1) with a parameter-optimized water profile (as shown in Figure 3b) is shown in Figure 6 (blue lines). The agreement between this calculation and the data is visibly better (and with smaller residual errors) than for the full RIM-1 calculation, but the residuals are still large with respect to the experimental uncertainties. The quality of fits for this structure are χ2 = 57.7 and R = 0.159. Finally, we use only the MD-simulated water profile from the RIM-1 simulation and fit the interfacial calcite structure (as shown in Figure 3c). In this case, the calculated reflectivity data show very good visual agreement with the data (black lines, Figure 6, “water only”), with consistently small residuals (i.e., errors between the measured and calculated reflectivity signals). The quantitative measures of agreement for this result are χ2 = 3.27 and R = 0.0427. These qualities of fit are, nevertheless, quantitatively slightly worse than those obtained from the traditional parameter-based fits to the XR data (χ2 = 1.67, R = 0.033). From this set of comparisons we conclude that there are significant discrepancies between the calcite portion of the MD-simulated RIM-1 structure and the experimental data, whereas the water profile from the same simulation provides semiquantitative agreement with the experimental data. 5035
dx.doi.org/10.1021/jp310943s | J. Phys. Chem. C 2013, 117, 5028−5042
The Journal of Physical Chemistry C
Article
Table 2. Comparison of Parameters Describing the MD-Simulated Interfacial Water Structures and That Derived from Analysis of XR Data, Including the Water Layer Height, z, Its rms Width, σ, and Its Occupation, Occa result XR
RIM-1
RIM-2
SM
ReaxFF
XR-revised (RIM-1)
XR-revised (RIM-2)
parameter
W1
W2
W3
W4
W5
z (Å) σ (Å) Occ z (Å) σ (Å) Occ z (Å) σ (Å) Occ z (Å) σ (Å) Occ z (Å) σ (Å) Occ z (Å) σ (Å) Occ z (Å) σ (Å) Occ
2.07(0.02) 0.27 (0.04) 1.12 (0.1) 2.37 0.17 1.1 2.01 0.22 1.5 2.21 0.22 1.22 2.39 0.18 1.2 2.36 (0.35) 0.24 (0.2) 1.2 (0.65) 2.08 (0.02) 0.22 (--) 1.32 (0.14)
3.37 (0.11)† 1.19 (0.06) 2.26* 3.4 0.17 0.9 3.3 0.25 0.6 3.15 0.25 0.85 2.95 0.18 0.2 3.27 (0.26) 0.34 (--) 1.00 (0.35) 3.54 (0.20) 0.33 (--) 0.66 (0.7)
6.7 (0.18)* 1.38 (0.10)* 2.26* 5.1† 0.7 1.41* 4.7 0.6 1.6 5.10† 0.69 1.68* 4.3 0.5 1.9 5.24 (0.34)† 0.73 (0.67) 0.56* 4.31 (0.74) 0.5 (--) 0.64 (0.31)
10.1 (0.24)* 1.54 (0.13)* 2.26* 7.2* 0.82* 1.41* 7.0† 1.3 1.55* 7.6* 1.06* 1.68* 6.45† 0.78 1.21* 6.07 (1.2)* 1.60 (1.1)* 0.56* 6.01 (0.31)† 1.15 (0.6) 1.55*
13.4 (0.30)* 1.70 (0.16)* 2.26* 9.3* 0.93* 1.41* 9.3* 1.84* 1.55* 10.1* 1.33* 1.68* 8.25* 1.05* 1.21* 6.90 (2.1)* 2.15 (1.4)* 0.56* 8.30 (2.2)* 1.74 (0.4)* 1.55*
Layers that are described through the “layered water model” have an occupation factor that is controlled by the layer spacing, an rms width σn = [σ02 + (n − 1)σbar2]1/2, and a constant layer spacing, z = z0 + (n − 1)dwater. The values derived from this model are indicated by “*”, and the height of the first water layer of this structural model is indicated by “†”. Uncertainties indicated by “--” denote values that were fixed during the optimization. a
parameter-based analysis of the XR data. This is valuable for understanding the source of the discrepancies and consistencies between the XR-derived and MD-simulated structures. Consequently, we describe the MD-simulated structures using the same parameters that are used to fit the XR data of the calcite−water interface. For the calcite lattice, this includes the vertical displacements of the Ca and CO3 ions, as well as the tilt, θCO3, and twist, ϕCO3, of the CO3 ion, each defined in the top six layers (Figure 7). For the interfacial water layer, the structure is defined by the height, occupation, and rms width of each water layer (Table 2). 3.3.1. Calcite Interfacial Structure. The comparisons of the near-surface calcite structural parameters obtained from the XR data and the simulations are shown in Figure 7. Generally speaking, the simulated Ca and CO3 displacements are substantially smaller than those observed by XR. For instance, the XR sees typical displacements of ∼0.05 Å near the surface that decay into the crystal. The RIM-1, RIM-2, and SM simulation show much smaller displacements (
