- © 2010 E. Schweizerbart’sche Verlagsbuchhandlung, D-70176 Stuttgart
Structure and vibrational frequencies (at the Γ point) of CaCO3 aragonite have been calculated from first principles, by using the hybrid Hartree-Fock/DFT B3LYP Hamiltonian, at different unit-cell volumes in the 185–242 Å3 range. By using the frequencies evaluated at such different volumes, the mode-γ Grüneisen’s parameters were estimated for each vibrational mode, and the zero point and thermal pressure contributions to the total pressure, at each volume and temperature, have then been determined by means of standard thermodynamic formulas, within the limits of the quasi-harmonic approximation. This allowed for the determination of (i) the equation of state at different temperatures; (ii) the thermal expansion as a function of pressure and temperature, and (iii) the evaluation of some thermodynamic properties (entropy and specific heat) together with their temperature dependences. Results were directly compared with relevant experimental data.
The agreement of the calculated frequencies with the experimental data, at variable pressure, shows that the ab initio simulation can reproduce, at a relatively low computational cost, the full vibrational spectra of crystalline compounds of mineralogical interest. Moreover, the elastic properties (bulk modulus in particular), thermal expansion and thermodynamic properties, which play an important role in the characterization and in the understanding of the stability relations between carbonate phases, at high-pressure and high-temperature conditions, can be satisfactorily estimated. Precisely, at room temperature and pressure conditions, the calculated bulk modulus was 64.7 GPa, to be compared with an experimental value of 65(1) GPa (mean of three different experimental determinations); the estimated thermal expansion was about 6.1 · 10−5 K−1, which is only slightly underestimated with respect to the experimental datum [6.3(2) · 10−5 K−1]; the calculated entropy (S) and the constant-pressure specific heat (CP) were 87.5 and 83.1 J mol−1 K−1 respectively, which are in close agreement with the experimental data [84(6) and 82.6 J mol−1 K−1, for S and CP respectively].
Like many other simple chemical compounds, CaCO3 is able to crystallize in a variety of polymorphic forms. Aragonite is one of the five known polymorphs of calcium carbonate; the other ones include vaterite and three forms of calcite. Aragonite, vaterite and one form of calcite occur naturally and are stable phases; the other phases are known to be metastable (e.g., Carlson, 1983).
Aragonite is by far the most common orthorhombic carbonate in its occurrence as it is the inorganic constituent of many invertebrate skeletons and sediments derived from them. It also occurs as a primary phase in high-pressure metamorphic rocks (e.g., Carlson, 1983). Its structure was first determined by Bragg (1924) and Wyckoff (1925); refinements were made by Dickens & Bowen (1971), by De Villiers (1971) and by Dal Negro & Ungaretti (1971). The space group of aragonite is Pmcn (a = 4.960, b = 7.964, c = 5.738 Å; Deer et al., 1996) and there are 4 CaCO3 formula units in the primitive cell.
In previous studies, calcite and aragonite phases have been studied and thermodynamic data were retrieved from spectroscopic measurements (Salje & Viswanathan, 1976), the vibrational frequencies were calculated by ab initio methods (Prencipe et al., 2004; Valenzano et al., 2006, 2007) and thermodynamic properties have been computed by using empirical interatomic potentials (Catti et al., 1993). As a further contribution to the understanding and to the interpretation of the physical and thermodynamic properties of the complex calcium carbonate system, in this work the vibrational spectrum at the Γ point of aragonite has been calculated by means of the periodic ab initio CRYSTAL06 code (Dovesi et al., 2006), by employing a high-quality basis set and the hybrid HF/DFT B3LYP Hamiltonian (Becke, 1988). The vibrational frequencies were computed at different volumes of the unit cell, which correspond to values of the pressure in the 0–30 GPa range. On the basis of the calculated vibrational frequencies, total pressures, bulk modulus, thermal expansion, heat capacity and entropy were obtained in the limit of the quasi-harmonic approximation (Anderson, 1995), through the evaluation of the dependence of the frequencies of the vibrational normal modes, at the Γ point, by the volume of the unit cell (mode-γ Grüneisen’s parameters). The dispersion effects in the phonon spectra at various pressures have not been considered.
In a first section, a brief description of the computational method is given, concerning geometry, basis set, Hamiltonian and algorithms used in the calculation of the equilibrium geometries and phonon frequencies. A second section is devoted to the description of the calculation of some thermodynamic properties of aragonite. The calculated results are finally compared with experimental data.
2. Computational details
Static energy, geometry and spectra calculations were performed by the ab initio CRYSTAL06 code (Dovesi et al., 2006), which implements the Hartree–Fock (Pisani et al., 1988) and Kohn–Sham (Kohn & Sham, 1965) self consistent field (SCF; Roothaan, 1960) method for the study of periodic systems, by using localized basis sets of Gaussian functions (GTO).
2.1. Basis set
In CRYSTAL, the multi-electronic wave function is constructed as an antisymmetrized product (Slater determinant) of mono-electronic crystalline orbitals (CO) that, in turn, are linear combinations of local functions (atomic orbitals; AO) centered on each atom of the crystal (Pisani et al., 1988). The basis set of aragonite was the same already used in the case of calcite (Valenzano et al., 2006); it consists of a 8-6511(21) contraction for Ca, where the (21) symbol indicates the presence of 2 d-type shells containing two Gaussian functions and a single one, respectively, and 6–311(11) and 8–411(11) contractions for C and O, respectively (Valenzano et al., 2006).
2.2. SCF Hamiltonian calculation
The hybrid HF/DFT B3LYP Hamiltonian (Becke, 1988, 1993), which contains a fraction (20 %) of the exact non local Hartree–Fock exchange, has been used; such Hamiltonian is widely and successfully employed in molecular quantum chemistry (Koch & Holtausen, 2000), as well as in solid-state calculations. B3LYP is known to provide very accurate values of vibrational frequencies (Catti et al., 1993; Prencipe et al., 2004; Valenzano et al., 2006, 2007). The DFT exchange and correlation terms of the Hamiltonian are numerically evaluated by integrating, over the cell volume, the relevant DFT functionals of the electron density and its gradient. In the present case, a (99,974)p pruned grid (corresponding to the keyword XXLGRID in the CRYSTAL input; Dovesi et al., 2006) was employed. The thresholds for the evaluation of the bielectronic integrals (ITOL1 to ITOL5; Dovesi et al., 2006) were set to 7 (ITOL1 to ITOL4) and 16 (ITOL5). The Hamiltonian matrix was diagonalized at 8 k points in the reciprocal space, which were chosen by setting to 3 the shrinking factor (IS) in the CRYSTAL input (Dovesi et al., 2006).
2.3. Geometry optimization
Geometry (cell parameters and fractional coordinates) was optimized at the equilibrium volume and at a set of fixed unit-cell volumes (185.46, 191.49, 194.54, 197.64, 200.75, 203.91, 208.74, 213.59, 220.20, 226.94, 236.44, 240.70, 242.15 Å3), by using the conjugate gradient method as implemented in CRYSTAL06 (Civalleri et al., 2001; Dovesi et al., 2006). Geometry optimization was considered converged when each component of the gradient (TOLDEG parameter in CRYSTAL06) was smaller than 0.00001 har-tree/bohr and displacements (TOLDEX) with respect to the previous step were smaller than 0.00004 bohr.
2.4. Frequency calculation
Vibrational frequencies and normal modes at the Γ point were calculated at the different cell volumes reported in Table 1⇓ within the limit of the harmonic approximation, by diagonalizing a mass-weighted Hessian matrix, whose elements are the second derivatives of the full potential of the crystal with respect to mass-weighted atomic displacements (Pascale et al., 2004). The threshold for the convergence of the total energy, in the SCF cycles, was set to 10−10 hartree (TOLDEE parameter in CRYSTAL06; Dovesi et al., 2006). Results are provided as supplementary material (Table S1) linked to this article and freely available on the GSW website of the journal at http://eurjmin.geoscienceworld.org/.
2.5. Thermodynamic properties
Bulk modulus, thermal expansion, specific heat and entropy were obtained in the limit of the quasi-harmonic approximation (Anderson, 1995), through the evaluation of the unit-cell volume dependence of the frequencies of the vibrational normal modes (mode-γ Grüneisen’s parameters) at the Γ point. Dispersion effects in the phonon spectra at various pressures were not taken into account as the large volume of the unit cell makes it impossible the required calculation of the vibrational frequencies in the case of super-cells, due to limitation of the available computational resources. On the other hand, the Grüneisen’s parameters corresponding to the zone-centre vibrational modes can reasonably be considered representative of the whole set of parameters, due to their large number which depends upon the number of atoms in the unit cell.
On the other hand, it has been demonstrated that highly reliable estimations of thermoelastic and thermodynamic quantities can be obtained even if dispersion effects are neglected, and even in the case of systems having relatively small unit cells (Ottonello et al., 2009a and b). The total pressure P at each cell volume (V) and temperature (T) is given by (Anderson, 1995):
where k and h represent Boltzmann’s and Planck’s constants respectively; EST is the static energy; νj is the frequency of the jth vibrational normal mode; nj(νj, T) is the number of phonons of frequency νj at the temperature T (according to the Bose-Einstein statistics; Equation 2) and γj is the mode-γ Grüneisen’s parameter as defined in Equation (3).
The total pressure P(V, T) (Equation 1) is a sum of three terms, that is, static pressure PST(V), zero point pressure PZP(V) and thermal pressure PTH(V, T).
For each normal mode, the relative γj Grüneisen’s parameter was determined through the analytic first derivative at each volume (pressure) of the second order polynomial resulting from the fitting of the numerical νj(V) curves.
Static pressures PST were determined by numerical derivatives, with respect to V, of a fourth order polynomial interpolating the EST(V) curve; zero point (PZP) and thermal (PTH) pressures were obtained by direct application of Equations (1)–(3).
2.6. Equation of state
The bulk modulus (KT) at a given cell volume (pressure) and temperature has been calculated through: (i) Birch–Murnaghan third order EOS (BM3) (Birch, 1966) fittings of the P(V, T) curves, and (ii) numerical derivatives of the P(V, T) curves, at the corresponding volume (NumDeriv); when using the latter method, K′ was obtained from the numerical derivative of KT with respect to P. More precisely, by using isothermal P(V, T) curves calculated at different T, from BM3 EOS we were able to obtain, for each temperature T, the bulk modulus (KT), the first derivative of it with respect to P(K′) and reference volume (V0), at any fixed pressure and temperature.
2.7. Thermal expansion
For any given volume V and temperature T, the thermal expansion α has been calculated either by: (i) numerical derivative, with respect to T, of the total pressure P(V, T) calculated by Equation (1):
where KT is the bulk modulus estimated at a given pressure and temperature; or (ii) using the VT values obtained by fitting isothermal curves, at different temperatures, by means of BM3 EOS:
2.8. Heat capacity and entropy
The heat capacity at constant volume (CV), as a function of T, has been calculated at zero pressure through the formula:
where Xj = hνj/kT and N0 is Avogadro’s number.
The specific heat at constant pressure (CP) was calculated at zero pressure by the standard formula:
where VT represents the cell volume at the temperature T. Depending upon the method used for calculating VT, KT and αT (see above) we obtained two different values for CP which, in the following, are indicated by the symbols CPBM3 and CPNumDeriv.
Entropy is given by:
3. Results and discussion
Starting from the optimized equilibrium geometry of the orthorhombic aragonite primitive cell (a = 5.0129 Å, b = 8.0587 Å, c = 5.8529 Å, V = 236.44 Å3, E = −3766.021280 hartree, at the B3LYP level), obtained in a previous work, new geometries have been obtained for reduced volumes of the unit cell which correspond to pressures higher than zero. Precisely, these new optimized geometries (Table 1⇑) correspond to static pressures (first term of formula 1) ranging from PST = 0 GPa at V = 236.44 Å3 (equilibrium geometry) to PST = 28.62 GPa at V = 185.46 Å3. Linear fits of parameter ratios b/a and c/a, as functions of volume, show a higher slope of c/a ratio with respect to the b/a one (Fig. 1⇓).
3.1. Grüneisen’s parameters, zero point and thermal pressures
A very good agreement has been obtained between the calculated Γ point vibrational frequencies at various pressures and the experimental data by Kraft et al. (1991). For example, at the equilibrium volume, we obtain a Δ̄ = 5.5 cm−1 mean difference between the calculated (B3LYP level) and experimental (Raman) frequencies (Carteret et al., 2009), which is in line with the average discrepancy observed in the analogous calculation on calcite (Valenzano, 2006).
In the low-frequency region of the vibrational spectrum (ranging from 50 to about 300 cm−1) the majority of the vibrational modes are associated to large positive Grüneisen’s parameters (≥1), whereas for frequencies above 700 cm−1 all of such parameters are in 0.1 – 0.2 positive range (Fig. 2⇓). The average Grüneisen’s parameter (γA) is positive and equal to 1.637. According to Equation (1), modes having a large γj produce major contributions to the vibrational pressure: in our case, the zero-point pressure is positive, since γA is positive, and it amounts to 0.713 GPa at the static equilibrium cell volume (VST 0 ), where PST = (∂EST/∂V )VST0 = 0. The zero total pressure at the athermal limit (where the thermal pressure is zero: T = 0 K) corresponds to an equilibrium volume (V0 = 238.41 Å3), which is larger than the static one (VST 0 = 236.44 Å3; at the B3LYP level).
At higher temperatures, a thermal pressure contribution depending on (i) the average number of phonons, (ii) the frequencies and (iii) the Grüneisen’s parameters associated to each normal mode (Equation 1), must be taken into account. Due to the Bose-Einstein distribution of phonons, at 400 K the thermal pressure is mainly determined (99.2 %) by phonon frequencies below 300 cm−1 (external modes); at 700 K the relative contributions to the thermal pressure are 53, 20 and 27 %, for external modes, bending and stretching modes respectively (Fig. 3⇓).
In Table 2⇓ some Grüneisen’s parameters, calculated at P = 0 GPa and T = 298 K, are compared with the available experimental ones (Gillet et al., 1993) corresponding to the six Raman active modes reported by Gillet et al. (1993). As it can be seen, a good agreement between calculated and experimental data was found. All of the Grüneisen’s parameters, calculated at the same P-T conditions, are reported in Table S2 (supplementary material).
In the investigated pressure range, all the frequency values increase with pressure (i.e., γj > 0) by following nearly linear trends: the only exceptions are the B1u and Ag modes at respectively νj = 850.7 and 851.6 cm−1, having γj = −0.07, and the B2g and B3u modes at respectively νj = 913.4 and 914.7 cm−1, having γj ? 0; see Table S2.
The calculated total pressure curve at room temperature (Equation 1), as a function of the scaled volume (according to the well known volume over-estimation provided by the B3LYP functional; see e.g., Prencipe et al., 2004), is shown in Fig. 4⇓ and compared with one experimental isothermal curve at 298 K (Martinez et al., 1996). A very good agreement between our calculated isothermal curve and the experimental one (Martinez et al., 1996) in the 298–573 K temperature range is observable.
By using Equations (2)–(5), it was possible to calculate thermodynamic properties as a function of T and P in the investigated volume range.
3.2. Bulk modulus
The values of the bulk modulus, calculated by means of the two different methods which have been previously described [third order Birch-Murnaghan BM3-EoS, and numerical derivative of the P(V, T) curve (NumDeriv)] are reported in Table 3⇓ as a function of T (298 – 973 K) at fixed P = 0 GPa.
By expanding the cell at volumes larger than the one corresponding to the static limit, we observed a mode whose frequency was rapidly decreasing as the volume was increased (soft mode); this is the Au mode at 35.2 cm−1, for the 242.15 Å3 cell volume (Table S1). Such mode reaches a zero value of the frequency, possibly corresponding to a zone-centre soft-mode phase transition, at a volume of about 245 Å3 (at B3LYP level). For this reason it was not possible to use our calculations to predict properties of aragonite at P = 0 GPa and T >600 K (since, at those conditions, thermal expansion would produce calculated cell volumes larger than 245 Å3). Indeed, by observing the experimental values of the bulk modulus measured at different temperatures (Martinez et al., 1996), which are reported in Fig. 5⇓ (and in the fourth column of Table 3⇑), two different slopes of the KT(T) curves are seen (precisely, before and after 673 K). This feature could be the sign of the presence of a phase transition which is theoretically predicted by our calculation.
For temperatures lower than 673 K, a very good agreement is observed between the calculated KNumDeriv and the experimental values (DiffRMS – 2 GPa). The agreement is somewhat less satisfactory if the calculated BM3 values are considered.
In Table 4⇓, the values of the bulk modulus (KT = K0T + K′P, where K0T is the bulk modulus at zero pressure, and K′ is the first derivative of K with respect to P), calculated by means of the two different methods (BM3 and NumDeriv) as a function of P (0 – 5 GPa), at the fixed 298 K temperature, are reported. At P = 0 GPa, the KT values obtained from the application of the NumDeriv (third column of Table 4⇓) and of the BM3 methods (second column of Table 4⇓) are very close; however, they become significantly different at higher pressures. By comparing the calculated KT values with the experimental ones (fourth column of Table 4⇓ and Fig. 6⇓; Salje & Viswanathan, 1976), it appears that a good agreement is reached between BM3 and the experimental data, in the 0 – 3.8 GPa pressure range.
The calculated BM3 and NumDeriv bulk moduli at zero pressure (K0T), at 298 K, were found to be in close agreement with respect to the experimental data (see Table 7a⇓). The bulk modulus in the 0 – 10 GPa pressure range and at higher temperatures (in the 298 – 773 K range) has been evaluated. Values of the coefficients (I and s) of the KT(T) = I + s · T linear fits of the bulk modulus as a function of temperature, at different pressures, are reported in Table 7b⇓; the negative slope (s) of the KT(T) curve decreases, as the pressure is increased, from −0.022 GPa/K, at P = 0 GPa, to −0.017 GPa/K, at P = 10 GPa (Table 7b⇓).
3.3. Thermal expansion
In Table 5⇓, the thermal expansion coefficient is reported, as a function of T, at a fixed P = 0 GPa, calculated by using either the VT values from the BM3 fittings of isothermal EoS (Equation 4b; second column of Table 5⇓), or numerical derivation (NumDeriv method; Equation 4a; third column of Table 5⇓), (see also Fig. 7⇓). In this case the NumDeriv values diverge from a linear trend above 473 K (Fig. 7⇓) so that our calculated values stops as 573 K. BM3 values seems more stable above 753 K, likely because the isothermal P(V, T) curves above 473 K (used by the BM3 method to get the values of the equilibrium volume at P = 0 GPa, at each fixed temperature) contain (P, V) points corresponding to P higher than 1 – 2 GPa, which are far from the soft-mode phase transition, which is predicted to occur at a temperature close to 600 K at 0 GPa (this likely causes the failure of the quasi-harmonic approximation which is at the basis of our calculations). Thus, at variance with the NumDeriv method, the BM3 values of the equilibrium volume at each temperature are obtained by using high-pressure data which are not affected by the error related to the presence of the soft mode. The comparison with the experimental values (fourth column of Table 5⇓; Salje & Viswanathan, 1976) gives similar trends for both methods in the 298 – 573 K range (Fig. 7⇓). The relative DiffRMS between the calculated and the experimental values in all the considered temperature range is smaller for BM3 method (4.4 %) than for NumDeriv method (7.4 %).
Moreover, the values of thermal expansion, obtained at temperatures above 600 K (at P = 0 GPa), are somewhat uncertain for the reason reported above related to the presence of the soft mode.
We also calculated the thermal expansion coefficient, α, at higher pressures (in the 0 – 10 GPa pressure range, and in the 298 – 773 K temperature range; BM3 method): the positive slope of the α(T) curve decreases, with increasing pressure, from 5.40 · 10−8 K−2, at P = 0 GPa, to 1.54 · 10−8 K−2, at P = 10 GPa (Table 7b⇓).
3.4. Heat capacity and entropy
As entropy S is concerned (Fig. 8⇓), a small DiffRMS ≈ 1 J mol−1 K−1 is generally observed along the whole T = 0 – 300 K range considered; the calculated S values reach a maximum of over-estimation, with respect to the experimental data, of + 1.7 J mol−1 K−1 at 100 K. In Table 6⇓, heat capacity values at constant volume (CV) calculated by Equation (3) (second column), CP values at constant pressure calculated by using the BM3 method (third column), together with the experimental data (fourth column; Staveley & Linford, 1969) are reported (see also Fig. 9⇓). All these properties have been calculated as a function of T (in the 20 – 300 K range) at a fixed P = 0 GPa. Table 6⇓ and Fig. 9⇓ show a very good agreement between the calculated CP data and the experimental ones from Staveley & Linford (1969), in all the investigated T range, with a mean DiffRMS <1 J mol−1 K−1 (maximum overestimation of +2.2 J mol−1 K−1, at 40 K). The agreement is particularly good in the 120–300 K temperature range, where DiffRMS < 0.2 J mol−1 K−1.
The calculated thermodynamic properties are summarized in Table 7a and b⇓. Table 7a⇓ (referring to data evaluated at T = 298 K) contains calculated values obtained by means of the BM3 (column 2) and NumDeriv (column 3) methods for the following properties: bulk modulus (K0T), thermal expansion (α), entropy (S), heat capacity (CV and CP); for comparison, some experimental values (from Staveley & Linford, 1969; Salje & Viswanathan, 1976; Martinez et al., 1996) are reported in columns 4–6.
Table 7b⇑ contains linear fittings (intercept and slope) of the data calculated by means of the BM3 (columns 2 – 3) and NumDeriv methods (columns 4 – 5) of the following quantities: bulk modulus as a function of temperature (298 – 773 K) at different pressures (0–10 GPa), thermal expansion as a function of temperature (298 – 773 K) and at different pressures (0 – 10 GPa). For comparison, in columns 6 – 7 some linear fits of the relevant experimental values (from Salje & Viswanathan, 1976; Martinez et al., 1996) are reported.
By using the α, KT and CV values of Table 7a⇑, it was possible to calculate the following Grüneisen ratio γ = αKTV/CV = 1.898 (Anderson, 1995), at P = 0 GPa and T = 298 K, which appears to be in good agreement (over-estimation of 3 %) with respect to the experimental room temperature value of 1.843 (Salje & Viswanathan, 1976).
The calculated BM3 values of the bulk modulus are in good agreement with the experimental ones, the difference between the calculated and the experimental data being less than 2 GPa in the whole pressure and temperature range considered (respectively, 0–10 GPa; 298–973 K).
The agreement of the estimated thermal expansion coefficient [α(T)] and its temperature dependence with the experimental data is reasonably good. In particular, the calculated value at room temperature and pressure is under-estimated by 2.3 · 10−6 K−1 with respect to the experimental datum (0.925 %), whereas the slope of the α(T) curve appears to be overestimated by 1.8 %; the magnitude of such a slope decreases as the pressure increases. Concerning heat capacity (CV and CPBM3), a very good agreement with the experimental values in the 120 – 300 K temperature range is obtained, and the average absolute difference |CPBM3 − CPEXP| is about 0.7 J mol−1 K−1.
The Au mode at 35.2 cm−1 (at a unit-cell volume of 242.15 Å3) reaches a zero value of the frequency for volumes above 245 Å3; this could possibly correspond to a zone-centre soft-mode phase transition (P = 0 GPa and T > 600 K conditions), which requires further investigations from both the theoretical and experimental sides.
- Received 4 December 2009.
- Modified version received 20 April 2010.
- Accepted 18 May 2010.