Article pubs.acs.org/JPCC
How acidic is water on calcite? M. P. Andersson* and S. L. S. Stipp Nano-Science Center, Department of Chemistry, University of Copenhagen, Denmark ABSTRACT: We have used density functional theory in combination with the COSMO-RS implicit solvent model to calculate the pKa for water adsorbed on various sites on a calcite surface. The mineral surfaces were modeled using 80 atom clusters. We found that, on kink sites, water is acidic with pKa as low as 3.5, implying that, at equilibrium, the kink sites are associated with one hydroxyl ion in addition to water. We also calculated pKa for ethanol and phenol adsorbed on calcite and compared the free energies of adsorption for ethanol and water. Ethanol was found to bind more strongly on the {10.4} terrace and on the obtuse step, whereas on kink sites, water (in the form of hydroxyl ions) binds more strongly than ethanol (in the form of ethoxide). Our results are relevant for modeling the dissolution and precipitation of minerals and show how the OH functional group behaves in the interaction between organic molecules and mineral surfaces.
■
INTRODUCTION A detailed understanding of the water−mineral interface can help explain many fundamental and technological phenomena, such as wetting, retention, nutrient depletion, and pollutant migration through soil. It is particularly important to have as accurate a description as possible for mineral interface properties when studying mineral adsorption, desorption, dissolution, and precipitation. These phenomena form the base for geological processes, such as weathering and ore formation, for biological processes, such as bone growth and tooth decay, for chemical processes of relevance to industry, and in the environment, such as transport of contaminants in the soil. The presence of organic molecules also influences mineral growth and dissolution. One dramatic example is the interaction of calcite and polysaccharides produced by one celled marine algae called coccolithophorids. Their complex sugars control the growth of interlocking, submicrometer scale calcite crystals and protect them from dissolving in seawater.1 Calcite, CaCO3, and water is one of the most studied mineral−fluid interfaces, particularly the {10.4} plane, for example.2 This face is the most stable surface termination of calcite, which can be seen from rhombohedrally shaped particles formed in nature and in synthetic products (e.g.3 Figure 1a). Surface science experiments on calcite are often conducted on surfaces dominated by {10.4} terraces with only few steps, kinks and other defects. Atomic force microscopy (AFM) experiments show that the calcite {10.4} surface recrystallizes in less than a day, even in air.2a This happens because a thin surface layer of water is present under ambient conditions when the relative humidity is more than 4%.4 X-ray photoelectron spectroscopy (XPS)2b reveals evidence of bonding environments that reflect H associated with CO3 and OH associated with Ca from the surface. Simulations to date on the hydrated, defect-free {10.4} plane show no sign of dissociated water, © 2012 American Chemical Society
Figure 1. Cluster models (80 atoms) for calcite: (a) {10.4} terrace, (b) acute step, and (c) a kink on an acute step. The top layer is represented as balls and sticks and the lower layers as tubes. The number of Ca−O bonds for a Ca ion on each of the sites is shown with blue lines. Calcium is yellow, carbon is gray, and oxygen is red.
indicating that undercoordinated sites such as steps and kinks might play a large role. Metal surface step and kink sites are certainly more reactive than close-packed surfaces for heterogeneous catalysis, where the dissociation energy barrier for CO is lower by about 100 kJ/mol.5 This happens because the dissociated fragments are stabilized. Recently, undercoordinated sites on the calcite surface have attracted interest in the modeling community. Particularly, monatomic steps have been studied using molecular Received: May 14, 2012 Revised: July 4, 2012 Published: July 27, 2012 18779
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
mechanics2k,6 and ab initio methods.7 The effect of other surface defect sites, such as single ion vacancies,7 has also been examined and it was found that water dissociates spontaneously on single carbonate ion vacancies. The concentration of vacancies on a terrace is however small, because to remove an atom from a terrace, five Ca−CO3 bonds must be broken (Figure 1a). Creating a vacancy at a monatomic step requires breaking four bonds. Removing a Ca atom from a kink only requires three bonds to be broken. Dissolution produces etch pits a single molecular layer deep, and although etch pits originate from single vacancies, steady state dissolution and precipitation take place at steps, with kink sites being particularly important.6,8 This is consistent with experimental9 and theoretical10 studies. The properties of the most common undercoordinated surface sites thus play a very important role in the dissolution and precipitation rates of a mineral. Furthermore, it is not in vacuum, but in solution, that these processes take place so a realistic theoretical approach requires a system where the behavior of water can be considered. Including water introduces further difficulties for computational studies, such as the presence of ions and the influence of pH. Molecular dynamics offers possibilities, but a new simulation is required for each new pH. If bond breaking is important to the process, then it becomes more computationally demanding to model accurately with molecular mechanics because special reactive force fields are required.11 Ab initio molecular dynamics can also capture proton transfer events, which have been predicted for water on calcite monoion vacancies.7 These simulations are unfortunately much more computationally expensive, and predicting the properties of the system as a function of pH, even more so. The influence of hydroxyls on various calcite surfaces was studied by Kerisit et al.,2d who replaced carbonate groups by hydroxyl as part of the dissolution process. These vacancies resemble kink sites and the addition of hydroxyls maintained surface neutrality. In this paper, we approach the question in a novel way and investigate several sites on the calcite surface using a cluster approximation and the sophisticated implicit solvent model COSMO-RS.12 We have studied the {10.4} face, acute, and obtuse steps and two types of kinks on each step: one with long Ca−Ca distance from the kink Ca to the step and one with short Ca−Ca distance (Figure 2). The other types of kinks13
differ mostly in the second layer configurations and are not investigated here. Because the primary mode of binding for water to calcite requires forming a Ca−O bond, we have chosen to investigate only the kinks that expose Ca ions. CO32− terminated kinks are left for future work. COSMO-RS is able to predict pKa very accurately14 using a linear free energy relationship between the acid and its anion. The predictions for regular molecules have a root mean square error as low as 0.5 pKa units. We use this feature, not for a regular molecule, but for an adsorbate that attaches to a mineral cluster so we can explore how water behaves on the various surface sites. We show that on kink sites where the hydroxyl ion can bind by bridging two Ca atoms, water is quite acidic, with pKa as low as 3 to 4. This implies that in the pH range in Ca− CO3 solutions where calcite does not dissolve (pH > ∼7) these sites are covered with hydroxyl ions, not water, consistent with experimental evidence.2b Supplementary calculations for ethanol and phenol confirm that this behavior is not limited to water. pKa for water and the two alcohols adsorbed at the various surface sites are similar, compared with when they are free in solution, also consistent with experiment.2k We are also able to calculate free energies of adsorption from solution and to determine the relative adsorption strengths of ethanol and water on the various surface sites of calcite in a 1:1 mixture of ethanol and water. We find that ethanol adsorbs more strongly than water on terraces and obtuse steps, whereas hydroxyl adsorbs more strongly than ethoxide on kink sites.
■
COMPUTATIONAL DETAILS All quantum chemical calculations on mineral clusters were performed using density functional theory (DFT) with Gaussian type orbitals and used the program ORCA, v 2.8.0.15 We used the BP functional, which is composed of Becke exchange16 and Perdew correlation.17 All calculations used the Ahlrich TZVP basis set,18 which is a triple-ζ basis set with polarization functions. The resolution of identity (RI) approximation was used with the auxiliary basis set TZV/J. The COSMO implicit solvent model19 for a perfect conductor (infinite refractive index) was used in all calculations. The use of COSMO-RS requires using an infinite refractive index because this is how the model is parametrized. This enables the change of solvent without performing new DFT calculations with a different dielectric constant. More details are presented in ref 12. The method BP/TZVP/COSMO provides input files suitable for use with the most accurate parametrization of the COSMO-RS extension to COSMO12,20 and used the COSMOtherm program version C21_0111. All cluster representations of the calcite surface were made with Materials Studio and used experimental geometries, except where noted otherwise. In all cases, the adsorbate and some surface atoms were relaxed to default thresholds. A pKa calculation with COSMOtherm requires the relaxed geometries of both the protonated and the deprotonated form in implicit solvent. There is no need for calculations on the proton; its energy and solvation properties are handled internally in the program. The pKa values are based on a linear free energy relation between the free energy of protonated and deprotonated states of a large set of inorganic and organic acids, ranging from very weak to strong acids. More details are presented in ref 14. For our pKa predictions of water on calcite, we used the parametrization for acids in water in the COSMOtherm program.
Figure 2. Cluster models (80-atom) for kinks on the acute step. The Ca−Ca distance from the kink to the inner step is short (4.047 Å) in (a) and long (4.990 Å) in (b). The Ca atoms defining the distance are shown as balls and the rest of the cluster as a tube model. Calcium is yellow, carbon is gray, and oxygen is red. 18780
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
By using an implicit solvent model, we assume that the local behavior dominates our properties of interest. Any collective behavior such as water ordering is not taken into account. On the other hand, unless the change in the protonation state significantly changes any solvent ordering, most of these effects ought to cancel. Because we only calculate equilibrium properties of water and hydroxyl ions that are both present in bulk water, we do not need to consider any dissociation events on the surface. It is possible that the presence of bicarbonate ions close to our surface site of interest modifies the relative stability of water versus hydroxyl ions, but it is beyond the scope of the paper to include these effects. All gas phase energies corresponding to solvated structures that are required for COSMOtherm to calculate accurate solvation energies were obtained by deriving single point energies on the relaxed geometries from the implicit solvent calculations. To improve the accuracy of the calculated equilibrium constants, we performed periodic plane wave calculations for water and ethanol adsorbed on the calcite {10.4} surface using Quantum Espresso.21 We used the Perdew Burke Ernzerhof (PBE) functional22 together with ultrasoft pseudopotentials taken from the pseudopotential library at the quantum espresso web site: Ca.pbe-nsp-van.UPF, C.pberrkjus.UPF, O.pbe-rrkjus.UPF, and H.pbe-rrkjus.UPF. We used a kinetic energy cutoff of 35 Ry and a density cutoff of 350 Ry. We use the PBE functional for the periodic plane wave calculations to minimize pseudopotential errors because the pseudopotentials were created with this functional. Improving the gas phase reaction energy of a reaction in solvent using a different method than the one used for the solvation properties is standard procedure.23 Different DFT functionals give very similar gas phase adsorption energies, especially when coupled with empirical dispersion in the DFT-D scheme.24 Water and ethanol were adsorbed on a {10.4} slab of calcite four molecular layers thick, with the simulation cell equal to 1 × 2 unit cells (80 atoms). The lowest molecular layer in the slab was held fixed at bulk positions during the optimization to represent an essentially infinite solid. We adsorbed one water or one ethanol molecule per simulation cell, which gives a coverage of 0.25 counted with respect to the number of calcium atoms. Geometry optimizations were made using the gamma point. Single point energies with a 2 × 2 × 1 Monkhorst Pack k-point mesh confirmed that gamma point energies converged within 0.01 eV. We also included empirical dispersion25 because this increases the accuracy of gas phase adsorption energies using DFT.26,24
Table 1. Convergence of Calculated pKa Values for the Calcite {10.4} Surface as a (a) Function of Cluster Size (All Atoms Frozen) and (b) Function of the Number of Relaxed Atomsa (a) cluster size (Å)
No. of atoms
pKa
6×9×3 9×9×3 6×9×6 9 × 12 × 3 15 × 15 × 3 (b)
60 80 90 100 180
15.1 14.0 13.9 14.0 14.1
cluster size (Å) 6 6 9 9 9 6 6 9 a
× × × × × × × ×
9×3 9×3 9×3 9×3 9×3 9×6 9×6 12 × 3
No. of atoms
No. of atoms relaxed
pKa
60 60 80 80 80 90 90 100
1 5 1 5 10 5 10 15
13.9 12.5 12.7 12.8 12.7 12.6 12.3 12.4
Cluster sizes are rounded to the nearest Å for the sake of clarity.
unchanged upon solvation. For surface atoms, charge is transferred from Ca to CO3 upon solvation, ≤0.11 e. The effect is largest on corners and edges. In a pKa calculation, it is the energy difference between the adsorbed water and adsorbed hydroxyl ion that determines the pKa. When comparing adsorption energies for two molecules, pure DFT is reliable, whereas for accurate gas phase energies, a dispersion correction and hybrid DFT gives better results.24 Relaxing only the Ca atom(s) in contact with the adsorbed water and OH is enough to achieve convergence with respect to the number of relaxed surface atoms (Table 1b). We assume these findings can be generalized to the other cluster geometries used for representing steps and kinks. All our final results are thus obtained for 80 atom clusters with the only surface atoms relaxed being those bound to water and OH. This results in a single Ca ion being relaxed for most clusters and two Ca ions for clusters where OH and water binds to two Ca ions in a bridging geometry. If there is extensive restructuring of the kink sites that extends more than a single bond into the surface, our assumed bulk terminated structure with only a few atoms relaxed would be a somewhat inaccurate description of the site. However, molecular dynamics simulations of calcite steps in water show only small variations in atomic positions compared with the bulk positions.2k Additionally, we always compare two species adsorbing in a similar fashion, which makes the calculation less sensitive to such effects. How Many Water Molecules Have to Be Included Explicitly? Any pKa predictions using COSMO-RS involves a calculation of the free energy of solvation so we needed to check how well the model performs for the constituent ions of calcite. In Table 2, we present the calculated free energies of solvation for CO32− and Ca2+. Observe that the first solvation shell of Ca has to be included explicitly to get reasonable solvation energies. The comparison between theory and experiment is quite good and comparable to what is achieved for commonly used force fields for the calcite−water system.2m
■
RESULTS AND DISCUSSION pKa Convergence with Respect to Cluster Size and Number of Relaxed Surface Atoms. We have investigated the convergence of the calculated pKa value of adsorbed water with respect to cluster size and the number of relaxed surface atoms for calcite {10.4}. We tested a range of cluster sizes and found that an 80 atom cluster is sufficient for pKa to reach convergence (Table 1a). This is consistent with our previous sensitivity analysis on magnesite {10.4}, where 80 atom clusters gave converged adsorption energies for water.24 The shape of the clusters was rhombohedral at all times. Removing ions from the edge to create a charged cluster induced only small changes in the adsorption energy of ethanol on magnesite (∼0.05 eV for 80 atom clusters).24 A mulliken charge analysis of a 6 × 9 × 6 Å cluster in gas and implicit solvent showed that the charge on the bulk atoms not exposed to implicit solvent remains 18781
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
bond to a carbonate ion than a Ca ion in the {10.4} terrace. Molecular dynamics simulations of oxygen coordination to the barium ion in barite show that there is a linear relationship between a coordination number of 8 in aqueous solution to 12 in bulk barite with various surface sites in between the two extremes.28 In bulk calcite, calcium is coordinated with 6 O, whereas the free calcium ion in aqueous solution has a coordination number of ∼7.29 Thus, we approximate the behavior of Ca by assuming that its coordination number remains 6 for all surface sites. The difference in pKa for 6 or 7 coordinated waters is 1.7 pH units. Without any experimental input, we treat this difference as a measure of the accuracy of our pKa predictions. Our calculated pKa value for 7 coordinated water molecules, 13.8, agrees well with the experimental hydrolysis constant for Ca, which is 12.8.30 Finding a Useful Descriptor for Future Studies. A parameter to describe the deprotonation of water on various surface sites is of considerable interest if we are able to provide other researchers with ways to estimate pKa without performing full pKa calculations. Therefore, we have investigated several descriptors that could be used to estimate pKa without performing full calculations with explicit waters. We have plotted pKa derived for explicit water as a function of several proposed descriptors in Figure 3 and sorted them according to their performance. The best agreement (Figure 3a) is achieved with calculations using a single water molecule corrected for Ca ion coordination. The maximum deviation from the linear fit is 1.6 pKa units. This offers confidence that the COSMO-RS treatment for the surrounding water molecules is quite good once the known deficiencies in pKa calculations using a free Ca ion exposed to the implicit solvent are accounted for (Table 4). Only the behavior of a single water molecule needs to be simulated for accurate pKa predictions. Second best correlation (Figure 3b), and probably most useful for rapid estimates, is simply counting the number of unfilled coordination sites on Ca ions, assuming a coordination number of six. It is assumed that water (and OH) can easily bridge the short, but not the long, Ca−Ca distance (Figure 2). When a bridging geometry is the most stable one, the number of unfilled coordination sites is assumed to be the sum of unfilled coordination sites from both bridging Ca atoms. This results in the {10.4} terrace having one unfilled coordination site, the obtuse and acute steps have two, the obtuse kink with the long Ca−Ca distance has three, the acute kink with the long Ca−Ca distance has four, and the obtuse and acute kinks with the short Ca−Ca distance have five unfilled coordination sites. The maximum deviation from the linear fit is ±1.9 pKa units, and the good agreement allows for a very reasonable prediction of equilibrium species present on any kind of calcite surface site without any calculations. This relationship should be of particular interest for future modeling studies of kink sites on mineral surfaces. We would get a better fit with an S-shaped curve than a linear curve, but for the sake of consistency, we present the best linear fit. It is also not clear from a study of a single system if this behavior is general enough to justify a different type of fit. While lowering the coordination of surface Ca2+ ions decreases the pKa value, an isolated Ca2+ ion has a pKa value that is substantially higher. We believe that it is the asymmetric association of CO32− ions around the Ca that more strongly favors OH− adsorption over H2O. The third best descriptor (Figure 3c) is the enthalpy difference between the adsorption energies for adsorbed hydroxyl and adsorbed water in an implicit solvent. Maximum
We therefore include enough water molecules to keep the Ca− O coordination number to six for the Ca ion of interest. The pKa values from the calculations using explicit waters are listed in Table 3. These serve as reference values for attempts to find easy to use descriptors that provide a reasonable way to estimate pKa data. To improve the accuracy of our single water pKa predictions, we also investigated how the calculated pKa for water bound to a free Ca ion varies as a function of number of explicit water molecules in the complex (Table 4). This Table 2. Calculated Free Energy of Solvation for Ca2+ and CO32− Ions Using COSMO-RSa ion
dGsolv (calc)
dGsolv (exp)
CO32− Ca2+ Ca(H2O)62+ Ca(H2O)72+
−1290 −1087 −1560b −1537b
−131535 −150535,c −150535,c −150535,c
a Energies are in kJ/mol. bCalculated from Ca2+ + (H2O)n (aq) → Ca(H2O)n2+(aq), n = 6, 7; (aq) implies COSMO-RS implicit solvent. c The experimental solvation energy for Ca2+ varies substantially. We have chosen the median value from ref 2m.
Table 3. pKa Values for Water Adsorbed on Various Surface Sites of Calcitea adsorption site
No. of explicit water molecules in calculation
pKa
1 2 2 3
12.7 8.5 13.2 4.9
4
3.5
3
10.3
4
3.9
{10.4} terrace acute step obtuse step acute kink (long Ca−Ca distance) acute kink (short Ca−Ca distance) obtuse kink (long Ca−Ca distance) obtuse kink (short Ca−Ca distance) a
There are enough water molecules to keep six oxygen atoms coordinating the associated Ca ions.
Table 4. pKa Values for Water Bound to a Free Ca2+ Ion as a Function of the Number of Water Molecules Included Explicitly No. of water molecules
pKa
Δ(pKa − pKa (6 water)
1 2 3 4 5 6 7
8.4 10.1 11.6 13.6 14.2 15.5 13.8
−7.1 −5.4 −3.9 −1.9 −1.3 0 −1.7
provides a correction factor for our calculations, where fewer than six oxygens are coordinated to the Ca ion of interest. This correction is required because free energies of hydration calculated using implicit solvation models are not very accurate for free cations such as Ca2+ (Table 2). Similar to free Ca ions, Ca on steps and kinks are most likely coordinated to more than one water molecule. Molecular dynamics simulations of the acute and obtuse steps for calcite {10.4} faces show that the coordination number remains the same, i.e., six.27 A calcium ion on a step binds to two water molecules because it has one less 18782
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
Figure 3. pKa determined using explicit water molecules maintains octahedral Ca coordination (i.e., 6 O for each Ca) as a function of several possible descriptors: (a) pKa calculated for a single water molecule, corrected using results for calculations with a free Ca2+ ion; (b) the number of O that would be required to complete the six-fold coordination on a Ca ion that is bound to OH−. For bridge bound species, the sum of the two bridging Ca ions is used; (c) the deprotonation enthalpy change in the implicit solvent; (d) the Mulliken charge on the hydrogen atom with the highest atomic charge in the implicit solvent; (e) the Löwdin charge on the hydrogen atom with the highest atomic charge in the implicit solvent; and (f) the deprotonation enthalpy in the gas phase. The lines are best linear fits and the R2 values are shown for each.
deviation from the linear fit is ±2.1 pKa units. The fact that there is a reasonable linear fit shows that the enthalpy contribution dominates the pKa. On the other hand, three values have deviations close to two pKa units from the linear fit, which shows that details in the entropic contribution from solvation are important to the final result and considerably more accurate results can be obtained from the COSMO-RS treatment. The gas phase (Figure 3f) enthalpy difference between water and OH is the descriptor with the worst performance. Linear fit is poor and maximum deviation is ∼5 pKa units. The enthalpy in the gas phase is thus a poor descriptor, whereas the enthalpy in the implicit solvent is rather good. This shows that polarization of the electronic degrees of freedom is an important parameter to include for considering the relative stabilities of the two species. Thus, calculations in implicit solvent can be useful for pKa estimation, whereas calculations using gas phase species are not. We have also plotted pKa as a function of atomic charge for the most positively charged hydrogen of water adsorbed in implicit solvent. Figure 3d, e shows that the partitioning schemes according to Mulliken and Löwdin are not as reliable descriptors as simply counting the unfilled coordination sites on the Ca ions or using the enthalpy difference. Agreement is not nearly as good as was observed for substituted anilins and phenols.31 We ascribe this to the fact that changes induced by
switching between adsorption sites are not subtle enough to be described as a simple substitution effect on the electron density. The maximum error from the linear fit is about three pKa units for the Mulliken charge and about four pKa units for the Löwdin charge. Other parameter testing showed that there is no correlation between pKa and the atomic charge of the hydrogen on water for calculations in the gas phase (not shown). This emphasizes the necessity of using a solvent model to obtain a reasonable description of the deprotonation tendency in water. Equilibrium Coverage. In Figure 4, we show the optimized geometry of the predicted equilibrium species, that is, H2O or OH−, at pH = 8, along with the calculated pKa for single water molecule adsorption on each site. These pKa values are corrected for the coordination number of the adsorption site using calculations for a free Ca ion from Table 4. We have chosen pH 8, which is near that of sea water, pH 8.3, because many calcite biomineralization processes take place in the sea, such as coccolith growth. Even though the pKa for explicit waters is more accurate, we present the Ca ion corrected single water calculations here because they show more clearly the adsorption geometry of the key species involved in the equilibrium. At all sites where the water and OH can bind at a kink site by bridging a short Ca−Ca distance, the pKa is significantly lower than 8 and, consequently, a hydroxyl ion, 18783
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
Figure 4. pKa for water for each of the investigated sites is shown along with the predicted protonation state at pH = 8. The sites are presented in order of decreasing pKa. The most stable protonation state is water on (a) an obtuse step, (b) a {10.4} terrace, (c) a kink on the obtuse step with the long Ca−Ca distance, and (d) an acute step. The most stable protonation state is OH− on (e) a kink on the obtuse step with short Ca−Ca distance, (f) a kink on the acute step with the short Ca−Ca distance and (g) a kink on the acute step with the long Ca−Ca distance.
8. With pKa = 8.8, we still predict that, at equilibrium, most of the Ca sites are bonded with water and only a single OH group but that about 14% of the sites have two hydroxyl groups adsorbed. We can also compare our results to surface complexation models for the calcite/water interface. The pKa of water deprotonating to OH− on the calcite surface as a whole is 11.7.32 This fits very well with our site specific pKa. If we assume that the majority of the surface is the {10.4} plane, with minor contributions from other sites, our estimate of total surface pKa is slightly lower than our value for the {10.4} surface, 12.7. How much lower depends on how many steps and kinks are present but the two predictions agree very well. The surface complexation model predicts a low coverage of OH− at all pH because at high pH the carbonate species takes precedence over OH, an effect which is not taken into account in our model. Ethanol and Phenol. Although water is the main focus of the paper, it is of interest to compare the behavior of water with that of organic compounds such as alcohols. Water, as H−OH can be viewed as a special case, the simplest OH compound, to compare with CH3CH2−OH and (CH)5C−OH. We performed pKa calculations on some of the sites with adsorbed ethanol and phenol in a manner completely analogous to the water calculations, to explore if general conclusions could be drawn from our results. We chose ethanol because experimental results exist for CaCO3 precipitation in ethanol/water mixtures.2k Calcium carbonate precipitated by mixing solutions of sodium carbonate and calcium chloride follows the phase transition: amorphous calcium carbonate to vaterite to calcite, whereas in mixtures containing at least 50% ethanol, aragonite remains stable for several months.33 The recrystallization behavior of calcite is also quite different in a mixture of ethanol and water than in pure water. Ethanol is also useful for comparing with the behavior of water because in aqueous
rather than water, is adsorbed at equilibrium. The binding geometry of water on the {10.4} surface compares well to molecular dynamics studies, where the water molecule forms one Ca−O bond and a hydrogen bond to an adjacent CO32− ion. The pKa values for three kink sites (Figure 4e−g) are so low that we decided to investigate the deprotonation of a second water molecule in the presence of a hydroxyl ion. We performed a calculation of the second deprotonation for the calculation of explicit water molecules on all step and kink sites where two or more water molecules can bind to the same Ca ion. For the acute kinks, the 80 atom cluster introduces an artifact because the second deprotonating water bridges one of the edge ions of the cluster, where on a larger cluster, a terrace atom would have been. For these two cases, we used the relationship shown by Figure 3b, where a very good linear correlation exists between the coordination number and pKa. Thus, to counter the artifact, we increase the pKa for these calculations by the slope of the curve, which is 2.6 pH units/ coordination number. The results are listed in Table 5 and show that, only at the acute kink, with a short Ca−Ca distance, is the pKa of the second water deprotonation important at pH Table 5. Calculated pKa Values for the Second Deprotonation in the Presence of a Hydroxyl Ion adsorption site
pKa
acute step obtuse step acute kink (long Ca−Ca distance) acute kink (short Ca−Ca distance) obtuse kink (long Ca−Ca distance) obtuse kink (short Ca−Ca distance)
14.6 12.6 14.1a 8.8a 14.0 11.9
a
Values increased by 2.6 pH units to account for the artifact caused by the edge of the 80 cluster model (details explained in the text). 18784
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
and water. The results are presented in Table 7, expressed as the change in free energy for replacing an adsorbed water molecule by an ethanol molecule:
solution, its pKa is close to that of water. Phenol is a good example of an alcohol with a significantly lower pKa than water. Experimental and calculated pKa in water for the three molecules are presented in Table 6. Calculated pKa values for
EtOH(solv) + H 2O(ads) → EtOH(ads) + H 2O(solv)
Table 6. Calculated and Experimental pKa for Water, Ethanol, and Phenol in Aqueous Solution at 25 °C
a
molecule
pKa expa
pKa calc
water ethanol phenol
15.7 16.0 9.8
15.5 16.2 9.7
(1)
Table 7. Free Energies of Adsorption for Replacing an Adsorbed Water by Ethanol in a 1:1 Mixture of Ethanol/ Water at Room Temperaturea
Taken from ref 14.
adsorption site
ΔGads(cluster)
ΔGads(cluster + gas phase plane wave energies)
terrace {10.4} obtuse step acute step acute kink (long Ca−Ca distance) acute kink (short Ca−Ca distance)
8.6 15.1 23.8 42.7
−15.5 −9.0b −0.3b 18.6b
490 36 1.1 0.057c
54.1
30.0b
0.00060c
Kads (replacing water at surface by ethanol)
a
All energies are in kJ/mol. bThe difference in enthalpy between the PBE-D plane wave calculations and the BP/TZVP cluster calculations on the {10.4} surface is used for the other sites as well. cUses K2 from eq 3.
The associated equilibrium constant is denoted K1. We know that, on magnesite, the gas phase energies converge differently with respect to cluster size for ethanol and water, with the adsorption energy of ethanol being underestimated by about 0.14 eV for an 80 atom cluster and the adsorption energy of water already converged.24 We have therefore performed plane wave calculations for ethanol and water adsorbed on calcite {10.4} using the PBE functional and added empirical dispersion (-D). The PBE-D method has been found to give as good results as both the hybrid and double hybrid functionals in previous studies for gas phase adsorption energies.26,24 The calculated change in enthalpy for reaction 1 using plane waves is −13 kJ/mol. Using clusters and the BP/TZVP method, the change in enthalpy is +10 kJ/mol, which is significantly different from the plane wave calculations. This energy difference is mostly related to the fact that the adsorption energy of ethanol is not converged with respect to the cluster size. On 80 atom magnesite clusters, the adsorption energy of ethanol goes down by 0.13 eV upon increasing lateral cluster size and by 0.04 eV upon increasing vertical cluster size, whereas the adsorption energy for water remains unchanged.24 The total change is 16.4 kJ/mol, which accounts for the major part of the 23 kJ/mol energy difference between the periodic and cluster methods. The remaining change most likely comes from the use of empirical dispersion, where the CH3 group of ethanol also contributes more favorably to the binding than the H of water. Improving reaction free energies in this manner, by performing higher quality quantum chemistry calculations for the gas phase species, is standard procedure and has been successfully used to obtain accurate absolute rate constants for polymerization processes in various solvents.23 We have assumed that the difference in reaction enthalpy (−25 kJ/ mol) between the plane wave calculations and the 80 atom cluster can be applied to all sites. The correction is included for all entries in Table 7.
Figure 5. pKa as a function of the number of unfilled O sites on Ca ions, assuming six-fold coordination. Calcite terraces have one unfilled O site, step edges have two, and kink sites three unfilled O sites. Results for calculations for water are shown as black circles, for ethanol, as red squares and for phenol, as purple triangles. The pKa for ethanol and phenol follow very similar trends to that of water.
some surface sites are shown in Figure 5. The general trend of lower pKa for more undercoordinated sites holds for alcohols as it does for water but merely counting the number of nonbonded sites is a considerably worse descriptor for ethanol and phenol. When alcohols are adsorbed in a bridge bound geometry on the kink sites, the number of unfilled coordination sites is 4 or 5 depending on the site. Here a significant decrease of pKa and a deviation from a straight line is observed. The effect is only seen for ethanol and phenol and not for water. This could result from less steric hindrance for adsorption of the anions compared to the protonated acid. Free Energy of Adsorption for Ethanol Compared with Water in a 1:1 Mixture of Ethanol and Water. A very useful feature of COSMOtherm is that once the COSMO files are derived, a solvent of arbitrary composition can be modeled at no additional computational cost. Recent experiments and calculations2k,l show that ethanol binds more strongly than water at the {10.4} surface and it hinders recrystallization on this face. There is evidence from atomic force microscopy and molecular dynamics simulations that water is able to penetrate a structured ethanol layer only at steps and recrystallization takes place on undercoordinated sites.2k,l We have therefore calculated the free energies of adsorption for ethanol and water on terraces, steps, and kinks in a 1:1 mixture of ethanol 18785
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
low as 3.5 on kink sites where water and OH− can adsorb by bridging two Ca ions. At equilibrium, a hydroxyl ion is adsorbed on most kink sites, not only water. Water, H−OH can be considered the simplest alcohol. Other alcohols are even more strongly acidic on these kink sites but their free energy of adsorption compared with water (or hydroxyl) decreases the more undercoordinated the surface site is. Our results are likely to find welcome application in future modeling studies of calcite surface sites, particularly using molecular mechanics, because the presence of hydroxyl at equilibrium could potentially change the energy landscape.
For some of the kink sites, the pKa values for both water and ethanol are below 8 and hence, the deprotonated forms are adsorbed on these sites. Because the system is assumed to be in equilibrium, we can use our calculated pKa to obtain the equilibrium constant for the following reaction where the adsorbed species are deprotonated: EtOH(solv) + OH−(ads) → EtO−(ads) + H 2O(solv) (2)
In this case, the equilibrium constant, K2, is related to K1 and the pKa for the adsorbed ethanol and water as K 2 = K1
K a(EtOH) K a(H 2O)
■
AUTHOR INFORMATION
Corresponding Author
(3)
*E-mail:
[email protected].
The results are shown in Table 7, where the terrace and steps have used K1 and the kinks have used K2. We note that adsorption of ethanol, compared with water, is strongest on the {10.4} terraces and weakens the more undercoordinated the surface sites become, despite the fact that the pKa for ethanol is lower than that for water on these sites. The equilibrium constant for the replacement reaction predicts that ethanol is bonded on the obtuse steps, while hydroxyl ions are associated with the kink sites. This is quite consistent with the experimental results, showing that recrystallization is inhibited on terraces, while it does occur on step edges and kinks. It is worth noting that our implicit solvent model does not take into account entropy effects of solvent structuring other than the immediate bonding to the calcite surface of water and ethanol. However, as long as hydration water molecules for Ca2+ ions are included, free energies of solvation for both Ca2+ and CO32− are reasonable (Table 2). Implications. The surface behavior of something as ubiquitous as water depends critically on the local surface structure of the solid phase. We have shown that the deprotonation tendency, as a function of surface geometry, can be predicted quite effectively using a combination of density functional theory and implicit solvent models. Our results are highly relevant for studies of dissolution and precipitation, which should take into account the presence of hydroxyl ions on the very sites where the process takes place. This prediction is also of great interest to molecular dynamics simulations of mineral/liquid interfaces because these most often need to define the bond structure and maintain it throughout the simulation. If highly undercoordinated sites are simulated with only water present, it is not the equilibrium situation that is modeled. This introduces an error, not least because electrostatic interactions would be altered by transfer of protons to solution. One example where this should be relevant is the initial crystal nucleation of amorphous calcium carbonate.34 This process involves small clusters that necessarily present a high proportion of low coordination sites, similar to the kink sites we have investigated here.
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/I001514/1] which funds the Materials Interface with Biology (MIB) consortium and under the Nano-Chalk Venture, supported by Maersk Oil and Gas A/S and the Danish National Advanced Technology Foundation (HTF). Computational resources were provided by the Danish National Centre for Scientific Computing (DCSC).
■
REFERENCES
(1) Hassenkam, T.; Johnsson, A.; Bechgaard, K.; Stipp, S. L. S. Tracking single coccolith dissolution with picogram resolution and implications for CO2 sequestration and ocean acidification. Proc. Natl. Acad. Sci. U.S.A. 2011, 108 (21), 8571−8576. (2) (a) Stipp, S. L. S.; Gutmannsbauer, W.; Lehmann, T. The dynamic nature of calcite surfaces in air. Am. Mineral. 1996, 81 (1−2), 1−8. (b) Stipp, S. L.; Hochella, M. F., Jr. Structure and bonding environments at the calcite surface as observed with X-ray photoelectron spectroscopy (XPS) and low energy electron diffraction (LEED). Geochim. Cosmochim. Acta 1991, 55 (6), 1723−1736. (c) Wright, K.; Cygan, R. T.; Slater, B. Structure of the (10(1)overbar4) surfaces of calcite, dolomite, and magnesite under wet and dry conditions. Phys. Chem. Chem. Phys. 2001, 3 (5), 839−844. (d) Kerisit, S.; Parker, S. C.; Harding, J. H. Atomistic simulation of the dissociative adsorption of water on calcite surfaces. J. Phys. Chem. B 2003, 107 (31), 7676−7682. (e) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliott, J. A.; Lardge, J. S.; Duffy, D. M. New forcefields for modeling biomineralization processes. J. Phys. Chem. C 2007, 111 (32), 11943− 11951. (f) Harding, J. H.; Duffy, D. M.; Sushko, M. L.; Rodger, P. M.; Quigley, D.; Elliott, J. A. Computational techniques at the organic− inorganic interface in biomineralization. Chem. Rev. 2008, 108 (11), 4823−4854. (g) Villegas-Jimenez, A.; Mucci, A.; Pokrovsky, O. S.; Schott, J. Defining reactive sites on hydrated mineral surfaces: Rhombohedral carbonate minerals. Geochim. Cosmochim. Acta 2009, 73 (15), 4326−4345. (h) Villegas-Jimenez, A.; Mucci, A.; Whitehead, M. A. Theoretical insights into the hydrated (10.4) calcite surface: Structure, energetics, and bonding relationships. Langmuir 2009, 25 (12), 6813−6824. (i) Lardge, J. S.; Duffy, D. M.; Gillan, M. J. Investigation of the interaction of water with the calcite (10.4) surface using ab initio simulation. J. Phys. Chem. C 2009, 113 (17), 7207− 7212. (j) Yang, M. J.; Stipp, S. L. S. Molecular dynamics study of the interaction between water or ethanol and calcite. Geochim. Cosmochim. Acta 2009, 73 (13), A1478−A1478. (k) Sand, K. K.; Yang, M.; Makovicky, E.; Cooke, D. J.; Hassenkam, T.; Bechgaard, K.; Stipp, S. L. S. Binding of ethanol on Calcite: The role of the OH bond and its relevance to biomineralization. Langmuir 2010, 26 (19), 15239− 15247. (l) Cooke, D. J.; Gray, R. J.; Sand, K. K.; Stipp, S. L. S.; Elliott,
■
CONCLUSIONS We have investigated the calcite−water interface using density functional theory calculations for calcite cluster models in implicit solvent. The COSMO-RS extension of the COSMO model allowed us to calculate the pKa for water adsorbing on various surface sites. Our estimate for the dominating surface, the {10.4} face, is 12.7, which agrees very well with surface complexation models. We also showed that acidity depends critically on the local surface structure. The pKa of water is as 18786
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787
The Journal of Physical Chemistry C
Article
J. A. Interaction of ethanol and water with the {101̅4} surface of calcite. Langmuir 2010, 26 (18), 14520−14529. (m) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. Derivation of an accurate force-field for simulating the growth of calcium carbonate from aqueous solution: A new model for the calcite−water interface. J. Phys. Chem. C 2010, 114 (13), 5997−6010. (3) Dickinson, S. R.; McGrath, K. M. Switching between kinetic and thermodynamic control: calcium carbonate growth in the presence of a simple alcohol. J. Mater. Chem. 2003, 13 (4), 928−933. (4) Bohr, J.; Wogelius, R. A.; Morris, P. M.; Stipp, S. L. S. Thickness and structure of the water film deposited from vapour on calcite surfaces. Geochim. Cosmochim. Acta 2011, 74 (21), 5985−5999. (5) (a) Norskov, J. K.; Bligaard, T.; Logadottir, A.; Bahn, S.; Hansen, L. B.; Bollinger, M.; Bengaard, H.; Hammer, B.; Sljivancanin, Z.; Mavrikakis, M.; Xu, Y.; Dahl, S.; Jacobsen, C. J. H. Universality in heterogeneous catalysis. J. Catal. 2002, 209 (2), 275−278. (b) Andersson, M. P.; Abild-Pedersen, E.; Remediakis, I. N.; Bligaard, T.; Jones, G.; Engbwk, J.; Lytken, O.; Horch, S.; Nielsen, J. H.; Sehested, J.; Rostrup-Nielsen, J. R.; Norskov, J. K.; Chorkendorff, I. Structure sensitivity of the methanation reaction: H-2-induced CO dissociation on nickel surfaces. J. Catal. 2008, 255 (1), 6−19. (6) de Leeuw, N. H.; Parker, S. C.; Harding, J. H. Molecular dynamics simulation of crystal dissolution from calcite steps. Phys. Rev. B 1999, 60 (19), 13792. (7) Lardge, J. S.; Duffy, D. M.; Gillan, M. J.; Watkins, M. Ab initio simulations of the interaction between water and defects on the calcite (10(1)over-bar 4) surface. J. Phys. Chem. C 2010, 114 (6), 2664−2668. (8) (a) Liang, Y.; Baer, D. R.; McCoy, J. M.; Amonette, J. E.; LaFemina, J. P. Dissolution kinetics at the calcite−water interface. Geochim. Cosmochim. Acta 1996, 60 (23), 4883−4887. (b) De Yoreo, J. J.; Zepeda-Ruiz, L. A.; Friddle, R. W.; Qiu, S. R.; Wasylenki, L. E.; Chernov, A. A.; Gilmer, G. H.; Dove, P. M. Rethinking classical crystal growth models through molecular scale insights: Consequences of kink-limited kinetics. Cryst. Growth Des. 2009, 9 (12), 5135−5144. (9) (a) Macinnis, I. N.; Brantley, S. L. The role of dislocations and surface-morphology in calcite dissolution. Geochim. Cosmochim. Acta 1992, 56 (3), 1113−1126. (b) Davis, K. J.; Dove, P. M.; De Yoreo, J. J. The role of Mg2+ as an impurity in calcite growth. Science 2000, 290 (5494), 1134−1137. (10) (a) Lasaga, A. C.; Luttge, A. Variation of crystal dissolution rate based on a dissolution stepwave model. Science 2001, 291 (5512), 2400−2404. (b) Larsen, K.; Bechgaard, K.; Stipp, S. L. S. Modelling spiral growth at dislocations and determination of critical step lengths from pyramid geometries on calcite {10(1)over-bar4} surfaces. Geochim. Cosmochim. Acta 2010, 74 (2), 558−567. (11) Gale, J. D.; Raiteri, P.; van Duin, A. C. T. A reactive force field for aqueous-calcium carbonate systems. Phys. Chem. Chem. Phys. 2011, 13 (37), 16666−16679. (12) Klamt, A.; Jonas, V.; Burger, T.; Lohrenz, J. C. W. Refinement and parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102 (26), 5074−5085. (13) Kristensen, R.; Stipp, S. L. S.; Refson, K. Modeling steps and kinks on the surface of calcite. J. Chem. Phys. 2004, 121 (17), 8511− 8523. (14) Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. First principles calculations of aqueous pKa values for organic and inorganic acids using COSMO-RS reveal an inconsistency in the slope of the pKa scale. J. Phys. Chem. A 2003, 107 (44), 9380−9386. (15) Neese, F. ORCA, version 2.6; University of Bonn: Germany, 2008. (16) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38 (6), 3098. (17) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 33 (12), 8822. (18) (a) Schafer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted gaussian-basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97 (4), 2571−2577. (b) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence
quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7 (18), 3297−3305. (19) Klamt, A.; Schuurmann, G. Cosmo, a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, No. 5, 799−805. (20) Mu, T. C.; Gmehling, J. Conductor-like screening model for real solvents (COSMO-RS). Prog. Chem. 2008, 20 (10), 1487−1494. (21) Gianozzi, P.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condensed Matter 2009, 21 (39), 395502. (22) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77 (18), 3865− 3868. (23) Deglmann, P.; Müller, I.; Becker, F.; Schäfer, A.; Hungenberg, K.-D.; Weiß, H. Prediction of propagation rate coefficients in free radical solution polymerization based on accurate quantum chemical methods: vinylic and related monomers, including acrylates and acrylic acid. Macromol. React. Eng. 2009, 3 (9), 496−515. (24) Andersson, M. P.; Stipp, S. L. S. Sensitivity analysis of cluster models for calculating adsorption energies for organic molecules on mineral surfaces. J. Phys. Chem. C 2011, 115 (20), 10044−10055. (25) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27 (15), 1787−1799. (26) Tongying, P.; Tantirungrotechai, Y. A performance study of density functional theory with empirical dispersion corrections and spin-component scaled second-order Moller-Plesset perturbation theory on adsorbate-zeolite interactions. J. Mol. Struct.: THEOCHEM 2010, 945 (1−3), 85−88. (27) Spagnoli, D.; Kerisit, S.; Parker, S. C. Atomistic simulation of the free energies of dissolution of ions from flat and stepped calcite surfaces. J. Cryst. Growth 2006, 294 (1), 103−110. (28) Stack, A. G. Molecular dynamics simulations of solvation and kink site formation at the {001} barite−water interface. J. Phys. Chem. C 2008, 113 (6), 2104−2110. (29) Jiao, D.; King, C.; Grossfield, A.; Darden, T. A.; Ren, P. Y. Simulation of Ca2+ and Mg2+ solvation using polarizable atomic multipole potential. J. Phys. Chem. B 2006, 110 (37), 18553−18559. (30) Wulfsberg, G. Principles of Descriptive Chemistry; Brooks/Cole Publishing: Monterey, CA, 1987. (31) Gross, K. C.; Seybold, P. G.; Hadad, C. M. Comparison of different atomic charge schemes for predicting pKa variations in substituted anilines and phenols. Int. J. Quantum Chem. 2002, 90 (1), 445−458. (32) Pokrovsky, O. S.; Mielczarski, J. A.; Barres, O.; Schott, J. Surface speciation models of calcite and dolomite/aqueous solution interfaces and their spectroscopic evaluation. Langmuir 2000, 16 (6), 2677− 2688. (33) Sand, K. K.; Rodriguez-Blanco, J. D.; Makovicky, E.; Benning, L. G.; Stipp, S. L. S. Crystallization of CaCO3 in water−alcohol mixtures: spherulitic growth, polymorph stabilization, and morphology change. Cryst. Growth Des. 2012, 12 (2), 842−853. (34) (a) Tribello, G. A.; Bruneval, F.; Liew, C.; Parrinello, M. A molecular dynamics study of the early stages of calcium carbonate growth. J. Phys. Chem. B 2009, 113 (34), 11680−11687. (b) Raiteri, P.; Gale, J. D. Water is the key to nonclassical nucleation of amorphous calcium carbonate. J. Am. Chem. Soc. 2010, 132 (49), 17623−17634. (35) Marcus, Y. Thermodynamics of solvation of ions .5. gibbs freeenergy of hydration at 298.15 K. J. Chem. Soc., Faraday Trans. 1991, 87 (18), 2995−2999.
18787
dx.doi.org/10.1021/jp304671k | J. Phys. Chem. C 2012, 116, 18779−18787