Yi Mao		Paul D. Siders
Department of Chemistry		Department of Chemistry
Northwestern University		University of Minnesota - Duluth
Evanston IL 60208		Duluth MN 55812
ymao@chem.nwu.edu		psiders@d.umn.edu

Keywords for Indexing: calcite, carbonate, calcium carbonate, protonation

Abstract

Calculated properties of calcium carbonate monomers and dimers are presented. Properties of the dimer are compared to properties of calcite and of the cleavage surface of calcite. Protonation and hydration of calcium carbonate are discussed. The calcium carbonate dimer is a model for some features of the calcite surface.

Introduction

Calcite, one of the most ubiquitous materials at the earth's surface, is an environmentally important mineral. It plays a vital role in the chemical regulation of aqueous environments, including lakes, soils, sediments, hydrothermal systems and the ocean. Via precipitation, dissolution and sorption reactions in the solution, calcite affects the quality of water.

Calcite has been the object of a number of theoretical and computer simulation studies. Empirical methods were used to model calcite by two- and three- body interatomic potential energy functions with some parameters fitted to structural, elastic, and vibrational properties. [1-6] Energy, bonding and structure in crystalline calcite have been studied by several non-empirical methods. Calcite properties have been calculated for collections of independent calcium and carbonate ions. [7,8] Skinner, et al. [9], used density functional theory and the local density approximation, with a basis of linearized augmented plane waves, and included the crystal's spatial periodicity. Catti, et al. [10], used the CRYSTAL [11] code to do all-electron periodic Hartree-Fock calculations. However, none of these calculations treated calcite reactivity, especially the most interesting and challenging part of the geochemistry of calcite: the kinetics of calcite surface reactions.

Ab initio calculations on cluster models have been successfully used to study mineral surface reactions. The molecular modelling of solids and surfaces enables the extension of the usual treatment of molecules to solid problems. This approach has contributed to the knowledge of silicate and aluminosilicate dissolution, including detailed mechanistic proposals.[12-14] We believe that cluster models of the calcite surface will help explain chemical processes that occur at calcite interfaces, especially how interactions between the surface and molecules affect bonding and reactivity.

Structure and Bonding in the Calcite Crystal. The calcite(CaCO₃) crystal has the space group R3C with Z=6 (six molecules per unit cell), a=4.9894 Å and c=17.039 Å for the hexagonal cell, as determined by X-ray diffraction. [15] (Fig. 1) The crystal structure consists of alternating layers of Ca²⁺ and CO₃^2- groups oriented perpendicular to the C(z) axis, with interionic bonding between layers. [16] The carbonate groups are planar. Carbonate groups within a layer are all coplanar and have identical orientations, and alternate layers are rotated 180° with respect to the layers above and below. Each Ca atom is coordinated by six oxygen atoms from different carbonate groups in an octahedral configuration. Each carbon atom is bonded to three oxygen atoms. Each oxygen atom is connected to two calcium atoms and one carbon atom.

Calcite Surface Structure. In recent years, experimental techniques, such as X-ray photoelectron spectroscopy (XPS),[17] low-energy electron diffraction (LEED), [17] and atomic force microscopy (AFM),[18,19] have been used to investigate the calcite surface on a molecular level. The various experiments suggest a common conclusion: the calcite surface is ordered and crystalline in air, in solution, and in vacuum, and there is minimal reconstruction of the surface. Over time, however, especially when exposed to water the surface does rearrange. [20]

However, the surface structure cannot be exactly the same as that in the bulk crystal. Bond reconfiguration at the freshly fractured calcite surface and relative rotation of some surface CO₃ groups were observed by LEED experiments. [17] When a crystal is fractured and bonds are broken, the newly exposed surface is out of equilibrium and it reacts rapidly with the contacting medium. If there is scarcity of atoms and molecules available near the surface for reaction, (e.g. in ultra-high vacuum, the experimental condition of XPS and LEED), the dangling bonds are satisfied as much as is possible by rearrangement of bonding configuration and surface structure. After cleavage, surface Ca atoms share electrons with only five O atoms or fewer. This probably results in stronger Ca-O bonds at the surface. Some carbonate groups are rotated relative to their positions in the bulk crystal.

Computational Methods

We used the GAMESS [21] and Gaussian92 [22] computer programs to carry out ab initio Hartree-Fock calculations, and to find minimum-energy geometries.

To guide our choice of a basis set, we calculated energy, geometry and vibration frequencies of the CaCO₃ monomer using numerous basis sets. The smallest sets were STO-3G [23] and 3-21G [24], and the largest sets were Partridge's near-Hartree-Fock basis sets [25]. The minimum-energy geometry of the CaCO₃ monomer has C_2v symmetry. The geometry in Fig. 2 was calculated using Ahlrichs' double-zeta basis set plus polarization and diffuse functions, "Ahlrichs+*" [26]. Fig. 3 shows the total energy versus the number of primitive Gaussian-type orbitals in the CaCO₃ basis. The energy varies widely for small bases and may be converging with the larger basis sets. The basis-set labels in Figs. 3, 4, and 5 are explained in Appendix A. Fig. 4 shows the carbonate symmetric-stretching frequency, and Fig. 5 shows a Ca--CO₃ stretch. From the basis-set comparisons we conclude that polarization functions are important. We also want to include diffuse functions because of carbonate's anionic character.

To further test basis set quality we calculated basis set superposition error (BSSE) for the reaction of calcium ion Ca²⁺ with symmetric planar carbonate to form C_2v CaCO₃. We calculated BSSE using the simple Boys and Bernardi [27] superposition method. The reaction energies in Table 1 and Fig. 6 were calculated at geometries optimized for each basis set. The BSSE decreases greatly from 3-21G+* to Huzinaga's MIDI+* [28] to Ahlrichs+*. Most of the improvement is in the description of calcium.

The Ahlrichs+* basis and GAMESS triple-zeta valence basis [29] (TZV+*) give the best results. We chose the smaller Ahlrichs+* basis for most of the rest of our calculations because a smaller basis allowed us to do more extensive calculations.

Table 1. Basis Set Superposition Error for Ca²⁺ + CO₃^2- --> CaCO₃

basis set:			3-21+G*	MIDI+*	Ahlrichs+*	TZV+*
# of primitive GTO for CaCO3:			66	67	79	104
	ghost atom orbitals	calc'n at opt'd geometry for
CaCO3	none	CaCO3	-934.535021	-934.865449	-938.961503	-939.315685
CO32-	none	CO32-	-260.929822	-260.940853	-262.123523	-262.400009
Ca2+	none	atom	-672.815983	-673.140114	-676.058198	-676.130965
E (H)			-0.789216	-0.784482	-0.779782	-0.784711
Ca2+	CO32-	CaCO3	-672.822805	-673.147640	-676.058665	-676.131124
(Ca2+)= E(Ca2+)-E(Ca2+{CO32-})			0.006822	0.007526	0.000467	0.000159
CO32-	none	CaCO3	-260.914673	-260.927925	-262.110817	-262.387652
CO32-	Ca	CaCO3	-260.917810	-260.928300	-262.111103	-262.387897
(CO32-)= E(CO32-)- E({Ca2+}CO32-)			0.003137	0.000375	0.000286	0.000244
BSSE=SUM()			0.009959	0.007902	0.000753	0.000403
E corrected= E +BSSE			-0.779257	-0.776581	-0.779029	-0.784307

The (CaCO₃)₂ Dimer Model

Quantum chemical ab initio methods are a tool to relate molecular structure to reactivity, or the energetics and dynamics of reaction. The infinite dimension of solids complicates the application of quantum chemical methods to solid state problems. There exist ab initio methods which obtain the crystal wave function for a periodic lattice (e.g., CRYSTAL [11]). Such calculations are time-consuming and are ill suited to aperiodic situations such as surface defects and bonding sites. One way to overcome the difficulties is to adopt finite models such as clusters of atoms or ions, real or hypothetical molecules, that can be treated by the same ab initio methods applied to molecules. [13,30] This approach has gained popularity, not only because of the readily available computer programs for calculation, but also because it has the virtue of stressing the connections between structures and properties. It is also in line with chemists' intuition that local interactions dominate the energetics and dynamics of reactions. This is a key assumption of our approach. Cluster models have been widely applied to solid- state and surface problems, including kinetics of oxide mineral-fluid reactions. [12- 14] The topics of adsorption of CO₂ at surface and step sites of MgO,[31] the dissociation of H₂ on the Ti surface [32] and the mechanisms of silicate dissolution [33] have been addressed by this method. Carbonate minerals are in some respects more complicated than simple metallic and oxide crystals: the variety of atoms and bonding is greater. Our dimer model of calcite is the first one of carbonate minerals. Its high symmetry and closed structure allow most of the calcium-oxygen bonds to be included in a small model.

The first step is to replace the infinite, periodic structure of the crystal surface by a finite set of atoms in a cluster. Intuitively, cutouts from the bulk are made first and regarded as molecules. Commonly, the dangling bonds are then treated at the boundary of such a cluster. There are various conceptual approaches to the termination of a cluster. The most common way is to use hydrogen atoms as terminators. Hydrogen termination has been widely applied to oxide minerals, [12- 14] in which covalent bonds dominate. The nuclear charge on the hydrogen atoms can be adjusted to mimic the Coulomb potential resulting from the rest of the crystal. The arbitrariness of termination of a cluster results in the necessity of using more than one model. First because the results of calculations on a cluster depend on its size, a hierarchy of cluster sizes should be used, e.g. from a one- layer model to a multi-layer model, and results are expected to converge as cluster size increases. Second, for a given size (if not very large), complementary structures arise from choosing either a metal or an oxygen atom as the center. [30] For a simple mineral like silica (SiO₂), this causes little trouble. However, in calcite the high coordination number of calcium and the existence of carbonate groups would lead to a large cluster even for a one layer model if the usual termination were employed. Moreover, calcite is an ionic crystal, [34] so the nonlocal aspects of the bonding (e.g. electrostatic interactions or long-range nonbonded interactions such as Van der Waals interactions) should be taken into account, which requires a multi-layer model. Directly extending the surface model from (hydro)oxide to carbonate mineral is difficult.

Another common way to terminate a cluster is to embed it in a model of the bulk ionic crystal. [30] This is the usual treatment of metallic and ionic crystals. The charges, fractional or integer, are determined to reproduce the long- range potential. Embedding is a good method to characterize the electronic states of a crystal and investigate impurity ions in a host crystal. [35] There are problems with charge embedding for studying chemical reactions on the crystal surface. Fixed point charges do not respond to geometry and charge changes on the surface. They also may not reflect the geometry and charge relaxation expected at a cut crystal surface even in the absence of surface reactions.

To avoid large models and multiple models and the lengthy calculations associated with them, we have chosen a small cluster of high symmetry and closed structure. We use neither terminating atoms nor embedding. We arrived at the model as follows: we cut two CaCO₃ molecules from a slice (which is the cleavage plane) through the calcite crystal. This surface pair is shown in Fig. 7. Full optimization of this dimer yielded our model shown in Fig. 8. In the bulk crystal, all oxygen atoms are bonded to two different calcium atoms, but in the cut-out molecules shown in Fig. 7 some calcium-oxygen bonds are not present. The two oxygen atoms for which no bonds to calcium appear in Fig. 7 form the bridging oxygen atoms by bonding to two calcium atoms in the optimized dimer.

Geometry. Fig. 8 shows the structure of our model, which is supposed to represent a cleaned, unreacted calcite surface. The geometry was fully optimized at different levels. The dimer has C_2h symmetry: an inversion center, a C₂ axis passing through two calcium atoms and a symmetry plane containing two carbon atoms and two oxygen atoms that connect with two Ca each. Each Ca atom is bonded to four oxygen atoms, two from the same carbonate group, which constitute a closed ring structure. There are two kinds of oxygen atoms: one (nonbridging O) is connected to one C and one Ca only, the situation on the surface; the other (bridging O) is bonded to one C and two Ca atoms, the same environment as in the crystal. In all the optimizations, the carbonate groups remain nearly planar, showing a rigid structure as indicated in other studies. Table 2 shows the model geometry obtained with different basis sets, compared with the structure in the crystal.

Table 2. Optimized geometry and energy of (CaCO₃)₂

	MIDI+*	Ahlrichs+*	TZV+*	calcite data
number of primitive GTO	134	158	208	n.a.
Hartree-Fock E (H)	-1870.0037	-1878.1902	-1878.8947
monomer E(H)	-934.8654	-938.9615	-939.3157
E for dimerization (H)	-0.2728	-0.2672	-0.2634
RCa-C (Å)	2.6493	2.6766	2.6745	3.21
RCa-O bridging	2.3392	2.3670	2.3709	2.35
RCa-O nonbridging	2.3500	2.3754	2.3750	n.a.
RC-O bridging	1.3213	1.3088	1.3099	1.28
RC-O nonbridging	1.2623	1.2551	1.2554	n.a.
MP2 corrections	__________	__________	__________	______
E(2) dimer (H)	-1.2677	-1.2706		n.a.
E(2) monomer (H)	-0.6889	-0.6853	-0.7216	n.a.
MP2 Corrected E for dimerization	-0.1628	-0.1676		n.a.

As indicated in Table 2, the calculated geometries are similar to the crystal structure. The carbonate groups have almost the same bond lengths and angles as in the bulk. The orientation of carbonate relative to calcium gives Ca-O distances as in the bulk. The favorable Ca-O interactions lead to short Ca-C distances. This is consistent with LEED evidence of Ca-C bonds in calcite near-surface layers. [17] The bond angles formed by the carbonate group and Ca differ significantly from those in the bulk. This agrees qualitatively with the proposed rotation of carbonate group relative to Ca at the surface. We propose to use the optimized-geometry calcium carbonate dimer as a model for local bonding on the calcite surface.

Atomic charges. Atomic charge is an important parameter to characterize the chemical bonding in ionic crystals. The charge distribution on individual atoms can be assigned to reproduce the electrostatic potential calculated from the Hartree-Fock wavefunction at a grid of points surrounding the atoms. [37-39] Common methods of choosing grid points were recently reviewed by Spackman. [39] The potential determined charges (PDC) in Table 3 were determined on Spackman's {3,5+}_3,0 grid, as implemented in GAMESS. Charges thus calculated are within 0.03 of those from the CHELPG and MKS-Connolly methods (also as implemented in GAMESS) for (CaCO₃)₂ and the Ahlrichs+* basis. For comparison, atomic charges calculated by means of a Mulliken population analysis of the wavefunction are included in Table 3, for the Ahlrichs+* basis set. Mulliken and potential-derived charges are nearly the same for calcium but not within the carbonate ion, presumably reflecting the ionic nature of the calcium-carbonate interaction and the delocalized bonding within the carbonate ion. In Table 3 we compare charges to those obtained previously for bulk calcite by different methods (Table 3). Models a-d in Table 3 are for atoms in the bulk crystal. We know of no previous theoretical studies that yielded atomic charges for atoms on the crystal surface.

Table 3. Comparison of atomic charges

Ca	C	O	References and Methods
1.83	0.88	-0.80(Onbr) -1.10(Obr)	(CaCO3)2 Mulliken with Ahlrichs+* basis
1.80	1.53	-1.05(Onbr) -1.24(Obr)	(CaCO3)2 PDC with Ahlrichs+* basis
1.63	1.38	-0.96(Onbr) -1.09(Obr)	(CaCO3)2 PDC with basis sets from ref's 36 and 40
1.865	1.075	-0.980	Model a, Periodic ab initio method, ref.10
2.00	0.95	-0.98	Model b, empirical method, ref 41
2.00	0.985	-0.995	Model c, empirical method, ref 1
1.642	1.041	-0.894	Model d, empirical method, ref 2

For both basis sets, the dimer's nonbridging O atom has a lower charge than the bridging O atom. This result supports the belief that the atomic charge (for an anion, the absolute value of the charge) on an atom increases with coordination number. The same trend is observed in X-ray studies. In Table 3, a similarity appears among charge distributions on the CO₃ group from different sources. Model a used ab initio periodic Hartree-Fock calculations (CRYSTAL code), and Ca,C and O atoms were represented by 22, 18 and 14 atomic orbitals, respectively, in the form of contracted Gaussian-type functions. Models b, c, and d are all empirical-potential approaches. All three empirical potentials use two-body exponential terms to express short-range repulsions. A three-body harmonic angular O-C-O potential was included in c and d only. In b and c the calcium charge was fixed at the formal +2e value. The charges of the rest of the atoms in b and c and all the atomic charges in d were obtained by a fitting to a variety of physical properties of calcite. The similarity of the charges from our dimer model and those from periodic Hartree-Fock and empirical bulk calcite models suggests that the dimer's atomic charges are reasonable for calcite.

Vibration frequencies. Vibration frequencies from Hartree-Fock calculations using the Ahlrichs+* basis set are reported in Table 4, and compared to experimental Raman data.[ 1, 2, 42, 43] The calcite crystal has point symmetry D_3d and five Raman active vibrations. Two of the E_g vibrations, as well as the A_1g vibration, have been assigned to internal vibrations of the CO₃^2- ion.[43] The carbonate groups in the surface model have almost the same structure as in the crystal, thus the comparison of the vibrational frequencies of internal modes is feasible. Our dimer has C_2h symmetry, so the symmetry labels for the corresponding frequencies are different. In Table 4, C_2h frequencies that correlate with the crystalline frequencies are listed, and averaged. The calculated data are direct Hartree-Fock results with no scaling corrections. The errors are as expected for Hartree-Fock calculations.

Table 4. CO₃^2- vibration frequencies (cm^-1)

Experimental	D3d symmetry	Calculated Hartree-Fock w. Ahlrichs+*	Dimer C2h symmetry	Average of calc'd frequencies	Percent error
714	Eg	753 756 808 817	Bu Ag Bg Au	784	10
1088	A1g	1170 1175	Ag Bu	1172	8
1432	Eg	1542 1554 1752 1773	Ag Bu Bg Au	1655	16

Protonation. As mentioned before, in the dimer model nonbridging O atoms are representative of O atoms on the surface, each bonded to only one calcium. Bridging oxygens represent oxygens in the crystal, each bonded to two calcium atoms. It is difficult to distinguish surface atoms from bulk atoms experimentally. However, theoretical studies reveal different reactivity for these two types of oxygen atoms.

Protonation at the calcite-water interface has been examined experimentally and theoretically. [44] It is a typical interfacial reaction. Protons are transported from the solution to the crystal surface, and adsorbed onto the surface. To model protonation, we attached one and two protons to the dimer. Table 5 shows attachment energies. The single-proton results in Table 5 are for optimized H⁺ coordinates on the frozen C_2h dimer geometry. Fully optimized single-proton calculations are queued for future calculation. The two-proton results in Table 5 are for fully optimized geometries, at least at local energy minima. Figs. 9, 10, 11, and 12, show the protonated calcium carbonate dimers.

Table 5. Protonated Calcium Carbonate Dimer

Hartree-Fock w. Ahlrichs+* Basis	1H+ on bridging O	1H+ on nonbridging O	2H+, 1 on each bridging O	2H+, 1 on each of two opposite nonbridging O
dimer geometry	frozen	frozen	locally optimized*	locally optimized
E (H)	-1878.515	-1878.549	-1878.845	-1878.857
E - E(dimer)	-0.325	-0.359	-0.655	-0.667
RO-H (Å)	0.959	0.958	0.947	0.956
RCa-O protonated O to nearest Ca	2.367 as in dimer	2.375 as in dimer	2.311	2.693
*dissociated to Ca(OH)2 and 2CO2 upon full optimization

Protons attach preferentially to the nonbridging oxygens. The bonding of a single proton at bridging and nonbridging oxygens is qualitatively similar. Fig. 13 shows the attachment energy versus O-H distance for both types of oxygens. The curves are alike except that the bridging O-H has higher energy, by 90 kJ/mol. Attaching protons to both bridging oxygens has a dramatic effect on geometry: the protonated dimer dissociates. Protonating two opposite nonbridging oxygens also weakens the dimer (the Ca-to-protonated O distances increase 0.3 Å) but not to the extent of dissociation. The distinction between protonation of bridging and nonbridging oxygen atoms should carry over to the calcite surface, where nonbridging oxygens are those most accessible to incoming water and cations.

Hydration. Fig. 14 shows the fully optimized geometry of water adsorption. Weak hydrogen bonding appears between the water and a nonbridging carbonate oxygen: the H-O_carbonate distance is 1.97 Å. The bond lengths at the non-hydrated Ca are almost unchanged: 2.37 Å to the carbonate oxygens. The Ca-O bond lengths at the hydrated calcium are slightly lengthened: 2.42 Å to the water O, 2.48 Å to the oxygen associated with water's H, and 2.40 Å to the other nonbridging carbonate oxygen. That Ca-O bond lengths are little changed by water adsorption at calcium suggests that hydration does not directly lead to the dissolution of calcite. Protonation has a greater effect on the Ca-O bond lengths.

The binding energy of one water molecule onto the surface is estimated as 110 kJ/mol, simply as the difference of Hartree-Fock energies. Although this calculated energy is subject to systematic error, it does indicate that the binding energy is not small. The experimental hydration energy of a calcium ion is 141 kJ/mol [32], for about six water molecules, so our calculated binding energy is too large.

Stipp, et al., [17] used X-ray photoelectron spectroscopy, a technique that probes local bonding environments at the mineral surface, to examine the calcite surface freshly cleaved and after hydration. Their study provides direct evidence for the formation of CaOH and HCO₃ on the calcite surface upon hydration. It also suggests that the hydrated surface is not significantly restructured.

Fig. 14 indicates that water adsorption does not involve the cleavage of H₂O into H bonded to carbonate and OH to calcium. To confirm this suggestion, we calculated the optimized geometry and energy of the dimer with OH on the Ca site and H on the carbonate site, as shown in Fig. 15. The geometry in Fig. 15 is optimized to a local energy minimum for H and OH, but it is not a global minimum for even the H and OH positions, because the energy of H₂O on a frozen dimer is lower by .048 Hartree. However, protonation of carbonate is known to elongate Ca-O bonds, as discussed above, so the dissociative adsorption of water yields a less stable but more reactive structure. This is contrary to the X-ray results, but may be consistent with Stipp's suggestion that "H₂O might sit with its O attached to underbonded surface Ca atoms and with one H directed toward underbonded O from CO₃ groups."[19]

Conclusions

The Ahlrichs+* and TZV+* basis sets are suitable for calculating CaCO₃ geometries and energies. We used mainly the smaller Ahlrichs+* basis and intend to continue using that basis for calculating properties of carbonate clusters. The energy-optimized calcium carbonate dimer incorporates features of bulk and surface calcite. Bond lengths, charges, and frequencies, and the carbonate geometry are reasonable. Nonbridging oxygens can represent dangling surface oxygens, while bridging oxygens represent oxygens bound as in bulk calcite.

Protonation of the dimer is energetically preferred on nonbridging oxygens. Protonation weakens the cluster, reflecting the role of protonation in dissolution. Water attached to the dimer model has lower energy when not dissociated. The water attaches with oxygen bonded to calcium, and a weak water-carbonate hydrogen bond. The calculations suggest that water should attach nondissociatively to the calcite cleavage surface in vacuum. Further calculations will investigate the effects of solvation and of calcite's Madelung potential on hydration and protonation.

Acknowledgments

We are grateful for financial support from the Petroleum Research Fund of the American Chemical Society and for a copy of GAMESS provided by the Quantum Chemistry Group of Iowa State University.

References are currently in a separate file

Appendix A Basis sets for CaCO₃.

Tables are in the text.

1. Table 1
2. Table 2
3. Table 3
4. Table 4
5. Table 5

Illustrations

1. Figure 1: CaCO₃ unit cell
2. Figure 2: Planar C_2v CaCO₃
3. Figure 3: CaCO₃ Energy versus number of primitive GTOs
4. Figure 4: CaCO₃ symmetric stretch versus number of primitive GTOs
5. Figure 5: Ca -- CO₃ stretch versus number of primitive GTOs
6. Figure 6: BSSE for Ca²⁺+CO₃^2- --> CaCO₃
7. Figure 7: CaCO₃ pair cut from the surface
8. Figure 8: (CaCO₃)₂ C_2h dimer
9. Figure 9: H⁺ on nonbridging oxygen of frozen (CaCO₃)₂ dimer
10. Figure 10: H⁺ on bridging oxygen of frozen (CaCO₃)₂ dimer
11. Figure 11: Optimized 2H⁺ on nonbridging oxygens of calcium carbonate dimer
12. Figure 12: Optimized 2H⁺ on bridging oxygens of calcium carbonate dimer
13. Figure 13: Energy of (CaHCO₃)⁺(CaCO₃)
14. Figure 14: Optimized H₂O plus calcium carbonate dimer
15. Figure 15: Optimized H⁺ and OH^- on frozen calcium carbonate dimer