Static lattice energy calculations of mixing and ordering enthalpy in binary carbonate solid solutions

Chemical Geology 225 (2006) 304 – 313 www.elsevier.com/locate/chemgeo

Static lattice energy calculations of mixing and ordering enthalpy in binary carbonate solid solutions Victor L. Vinograd a,*, Bjo¨rn Winkler a, Andrew Putnis b, Julian D. Gale c, Marcel H.F. Sluiter d a

d

Institute of Mineralogy, University of Frankfurt, Senckenberganlage 30, 60054 Frankfurt/Main, Germany b Institute of Mineralogy, University of Mu¨nster, Correnstrasse 24, 49149, Mu¨nster, Germany c Nanochemistry Recearch Institute, Curtin University of Technology, U1987, Perth, 6845, Western Australia Delft University of Technology, Department of Materials Science, Rotterdamseweg 137 2628 AL Delft, The Netherlands Accepted 16 August 2005

Abstract Enthalpies of mixing and ordering in the binaries of the quaternary carbonate (Ca,Mg,Fe,Mn)CO3 solid solution have been modeled using static lattice energy minimization calculations. A set of self-consistent empirical potentials has been specially developed for this task. The calculations illustrate the importance of size mismatch in determining magnitudes of the enthalpies of mixing and ordering in the binaries. The enthalpy effects show positive correlation with the predicted temperatures of the order/ disorder transitions in intermediate compounds with the dolomite structure. It is shown that the general behavior of the mixing enthalpy in high- and low-T limits can be constrained with a minimum set of sampled configurations. D 2005 Elsevier B.V. All rights reserved. Keywords: Carbonates; Mixing properties; Static lattice energy calculations

1. Introduction Activity–composition relations in the quaternary carbonates (Ca,Mg,Fe,Mn)CO3 play an important role in modeling phase equilibria in various geochemical systems. These relations at low temperatures are complicated due to the appearance of dolomite-like compounds, in which different cations alternate in layers parallel to [110]. The main direct experimental information on the thermodynamics of mixing and ordering in solid solutions is provided by solution calorimetry. * Corresponding author. Tel.: +49 69 79822764; fax: +49 69 798 22101. E-mail address: [email protected] (V.L. Vinograd). 0009-2541/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.chemgeo.2005.08.023

Experimental constraints are available for the Ca–Mg and Ca–Fe carbonates (Navrotsky and Capobianco, 1987; Chai et al., 1995; Navrotsky et al., 1999). The data on the other four binary systems are much less complete. Recently, Burton and Van de Walle (2003) have shown that the energies mixing and ordering in Ca–Mg carbonates can be constrained with state-of-theart ab initio calculations. Similar calculations could be applied to other carbonate binary systems. However, even with the present advancements in computer power, such calculations would require significant effort and, therefore, it is desirable to reduce the set of sampled configurations to minimum. Sluiter and Kawazoe (2002) and Sluiter et al. (2004) have shown that the mixing enthalpy in alloys and oxides in the disordered limit can be constrained with the impurity method. This

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

method is based on noting that the Margules parameters, Wij , in the subregular equation for the mixing enthalpy Hmix ¼ x1 x2 ðx1 W21 þ x2 W12 Þ;

ð1Þ

where x 1 and x 2 are mole fractions of the components of a binary system, are equal to the slopes of the enthalpy curve in the vicinity of the end-member compositions: W21 ¼

dHmix j ¼0 dx2 x2

¼ 

and W12 ¼

dHmix j ¼ 1: dx2 x2

dHmix j ¼0 dx1 x1 ð2Þ

The derivatives can be well estimated from finite changes in the static lattice energy and composition of an end-member after one exchangeable cation in its supercell is replaced with a cation corresponding to the opposite end-member. Thus, W 12 should be understood as a partial increase in the total enthalpy after addition of one mole of atoms of type 1 into the pure phase composed of atoms of type 2. Practically, one generates a supercell of a pure mineral phase, e.g. calcite, and substitutes one Ca2+ cation with Mg2+. The normalized difference in the enthalpy and composition before and after the incorporation of the impurity atom provides the enthalpy slope and the value of one Margules parameter. The second parameter is constrained with the analogous calculation applied to the other end-member, e.g., magnesite. In order to estimate the slope accurately, the supercell should contain a sufficiently large number of exchangeable sites. The simple subregular behavior of the enthalpy of mixing is perturbed, however, at low temperatures due to cationic ordering. Here we investigate the possibility of extending the method of Sluiter and Kawazoe (2002) to binary systems with ordered intermediate compound, e.g. calcite–dolomite–magnesite, and illustrate the new version of the method with the help of empirically constrained force-field calculations. The empirical approach is preferred here over a more rigorous ab initio calculation because our main aim is just to test the predictive power of the new sampling scheme. We show that the general shape of the low-temperature enthalpies in the calcite–magnesite, calcite–siderite and calcite–rhodochrosite systems can be well predicted with the help of the impurity calculations analogous to those applied to the end-members. Specifically, the low-temperature enthalpy can be described with the help of two subregular equations analogous to Eq. (1), where the composition variable x 3

305

(0 b x 3 b 1) corresponds to the mole fraction of the ordered phase in the two solid-solution segments: x3 ¼ 2x2 ; 0Vx2 V0:5 x3 ¼ 2ð1  x2 Þ; 0:5Vx2 V1

ð3Þ

The Margules parameters, which describe the lowtemperature enthalpy, are defined similarly to Eq. (2). Note that the differentiation is now performed with respect to the new composition variable. For example, in the segment 0 b x 2 b 0.5, the Margules parameters can be calculated as follows W31


dHmix

¼ dx3
x3 ¼0

and W13


dHmix

¼  : dx3
x3 ¼1

ð4Þ

The impurity calculations thus allow determination of the enthalpy of mixing both in the high and low temperature limits. The difference between the high and low-temperature curves gives the enthalpy effect of ordering. To show that the impurity calculations predict realistic mixing properties we perform here static lattice energy calculations not only for the impurity compositions, but for a larger set of supercell configurations. These calculations permit extracting of the cluster interaction energies, the Js, which in turn permit constraining the energy of any configuration. The mixing and ordering phenomena can be then studied using Monte Carlo simulations and the results of these calculations can be compared with the predictions of the impurity method. The idea of supercell calculations has been discussed in a number of recent publications (e.g. Dove et al., 1996; Myers, 1998; Becker et al., 2000; Becker and Pollok, 2002; Dove, 2001; Bosenick et al., 2001; Warren et al., 2001; Vinograd et al., 2004a,b). An effective Monte Carlo algorithm requires that the enthalpy of any configuration is calculated fast. This is best achieved with the help of the Js formalism (e.g. Dove, 2001). The idea of this approach is based on mapping the set of sampled configuration energies onto a set of the pairinteraction constants. Each constant characterizes the enthalpy of the exchange reaction AA þ BB ¼ 2AB

ð5Þ

for a pair of neighbors located at a distance r n . Accurate description of the enthalpy with the Js formalism requires consideration of a wide range of pair interactions. This might not be easily achieved due to limited size of the simulation cell. Recently, Vinograd et al. (2004b) have shown that the accuracy in the description

306

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

of ordering and mixing effects in Ca–Mg garnets significantly improves when a configuration-independent term is added to the Js expansion: X ðnÞ DE mix c1=2 zðnÞ PAB J ðnÞ n

þ x1 x2 ðA21 x1 þ A12 x2 Þ:

ð6Þ

Here z (n) is the coordination number with respect to a neighbor located at the n-th distance and PAB(n) is the probability to find an AB pair on the two neighboring sites. The first term in this equation corresponds to the standard Js formalism, while the second term describes the configuration independent term. The composition-dependent polynomial in Eq. (6) has the form of Eq. (1); however, the physical meaning of A ij parameters is different from that of Wij parameters in Eq. (1). It has been shown that A ij parameters in the garnet system reflect the strain energy arising in the response to the compression and expansion of the endmembers caused by their mutual accommodation to the common lattice geometry at an intermediate composition. In the case of cubic symmetry, the A ij parameters are simple functions of the bulk modulus of the endmembers and of the square of the difference between their volumes. In this study, we adopt the same form of the volume deformation energy as in the study of Vinograd et al. (2004a). We realize, however, that in rhombohedral carbonates, the A ij parameters depend on the elastic stiffness tensor in a more complex way than in cubic garnets. 2. The potential set Before we describe the interatomic potentials used here for constraining the force-field calculations, we note that this study is not aimed at a quantitative representation of phase relations in carbonate binaries. Our main objective is to illustrate the power of the impurity method and thus to provide methodology for future more precise ab initio studies. To illustrate the general applicability of the method, we will show that using the minimum set of calculations in the impurity range one can estimate correctly the enthalpies of mixing and ordering in any carbonate binary system. To show that this result is not necessarily restricted to carbonates, we will use a set of potentials, which is transferable within carbonates, oxides and silicates. Although many superior sets of empirical potentials for carbonate systems have been developed (e.g. Fisler et al., 2000; Archer et al., 2003), none of them have been shown applicable

outside the limits of the carbonate group. Therefore, we use here a transferable set. This set has fewer internal parameters in comparison to the sets of Fisler et al. (2000); Archer et al. (2003), but permits sufficiently accurate description of structural and elastic properties of carbonates. The main requirement to a transferable set is that the potentials fitted to physical properties of pure oxides such as periclase and stishovite should automatically predict the corresponding properties of a complex oxide such as MgSiO3 perovskite. Previous attempts to develop such sets were based on the assumption of formal charges on cations and anions (e.g. Price et al., 1987; Winkler et al., 1991; Bush et al., 1994) and were only partially successful. Vinograd et al. (2004b) suggested recently that the lack of transferability could be related to the difficulty in fitting short cation–cation distances in dense oxides: at short distances, the repulsion of formally charged cations becomes too strong and the fitted cation–cation distances markedly exceed experimental values. It was shown that the repulsion can be effectively suppressed by scaling the charges of all cations and anions down. It was shown that the scaling

Table 1 The interatomic potentials used in the present study Buckingham Atom

Atom

A (eV)

˚) B (A

Si (core) Al (core) Ca (core) Mg (core) Fe (core) Mn (core) C (core) O (shell)

O O O O O O O O

1096.42 1226.70 2881.04 925.135 1242.278 1427.82 3517.73 614.709

0.299937 0.289261 0.281475 0.297097 0.289621 0.287978 0.187705 0.301585

(shell) (shell) (shell) (shell) (shell) (shell) (shell) (shell)

˚  6) C (eV A 0 0 0 0 0 0 0 27.0667

Three-body potential Atom

Atom

Atom

K (eV deg2)

G (deg)

O O O O O O O O

Si [4] Si [6] Al [4] Al [6] Fe [6] Mn [6] Mg [6] C [3]

O O O O O O O O

3.79 3.77 0.6686 1.8864 0.4497 0.8330 1.4132 2.1964

109.47 90.00 109.47 90.00 90.00 90.00 90.00 120.00

(shell) (shell) (shell) (shell) (shell) (shell) (shell) (shell)

(shell) (shell) (shell) (shell) (shell) (shell) (shell) (shell)

Spring Atom

Atom

˚ 2) k (eV A

O (core)

O (shell)

54.7039

For the three-body potentials, the coordination number of the pivot atom is given in parenthesis.

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

307

Table 2 ˚ ) of carbonates The observed and predicted lattice constants (A Mineral name

a (experimental)

c (experimental)

a (calculated)

c (calculated)

Calcite Dolomite Magnesite Rhodochrosite Siderite

17.061 16.002 15.018 15.635 15.380

4.988 4.807 4.634 4.768 4.692

17.081 16.154 15.031 15.647 15.341

5.042 4.814 4.618 4.747 4.694

factor of 0.85 is sufficient for bringing the predicted cell parameters and elastic stiffness constants of mono oxides such as stishovite and corundum in good agreement with the experiment. When this task is achieved, the potentials become automatically transferable to complex oxides and only slight refitting is needed to achieve an overall consistency. Here we present the new scaled-charge potentials and demonstrate their performance in predicting relative magnitudes of the enthalpies of mixing and ordering in different carbonate binary systems. The new potentials (Table 1) adopt the mathematical form of the THB set (Price et al., 1987; Winkler et al., 1991), but assume that the charges on cations and anions are lower than the formal ones by the factor of 0.85. Oxygen is modeled as a two-particle (core-shell) system with the charges 0.746527 (core) and  2.446527 (shell). The cutoff radius for the Bucking˚ . The cation–oxygen (M–O) ham potential is set at 12 A and oxygen–oxygen (O–O) cutoffs of the three-body ˚ , respectively. The potentials are set at 2.3 and 3.4 A core and shell charges, the coefficients of the Buckingham M–O and O–O potentials, the O (core)–O (shell) and the three-body O–M–O potentials have been fitted to structural and elastic stiffness constants of a large number of oxides and silicates in the system Si4+, Al3+, Mg2+, Ca2+, Fe2+, Mn2+, O2. The C4+–O potential was fitted to the structural and elastic constants of calcite, magnesite, siderite and rhodochrosite (in this fit, the other potentials were kept fixed). The structural and the elastic stiffness constants have been taken from the American Mineralogist Crystal Structures Data Base

(R.T. Downs, http://www.geo.arizona.edu/AMS) and the review of Bass (1995), respectively. The fitting was performed with the relax-fitting algorithm as implemented in the GULP program (Gale, 1996). Table 2 compares the predicted and observed lattice constants of the carbonates related to this study. 3. The impurity calculations of the high- and low-temperature enthalpy Table 3 lists the predicted Margules parameters for the six binary carbonate systems. The first two columns give the parameters for the high-temperature enthalpy and the last four columns give the values for the two branches of the low temperature enthalpy. For Mn–Fe, Fe–Mg and Mn–Mg binaries, we list only the hightemperature enthalpy. Significantly ordered states in these systems are probably inaccessible due to slow kinetics at the very low temperatures, at which they are stable. The high-temperature enthalpies are compared in Fig. 1. The mixing enthalpies of Ca–Mg, Ca– Fe and Ca–Mn systems are separately shown in Fig. 2a and c. The results of the impurity method calculations are shown as solid lines. In all calculations, we have

Table 3 The predicted Margules parameters (kJ/mol cations) for the six binary carbonate systems System

W 12

W 21

W 13

W 31

W 32

W 23

Ca–Mg Ca–Fe Ca–Mn Mn–Mg Fe–Mg Mn–Fe

35.03 22.07 17.87 5.13 1.67 0.82

51.60 32.58 23.16 5.40 1.86 0.93

34.68 23.58 16.51

41.19 26.79 19.11

42.86 27.01 18.22

50.51 31.11 23.40 Fig. 1. The high-temperature enthalpy of mixing in the six binaries of the quaternary carbonate solid solution predicted with the impurity method of Sluiter and Kawazoe (2002).

308

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

used 4  4  1 supercell containing 96 exchangeable sites. 4. The supercell calculations of the mixing and ordering enthalpies 4  4  1 supercell has been also used to calculate static lattice energies of a large number of randomly generated Ca–Mg, Ca–Fe and Ca–Mn configurations. These energies (after subtraction of the linear combination of the end-member energies) are also plotted in Fig. 2a and c. Each binary system was sampled with 25, 100 and 25 randomly generated configurations at the compositions of 0.25, 0.5, and 0.75, respectively. The calcite–magnesite system has been additionally studied at 0.33 and 0.66 compositions. In this system, we have also calculated static energies of several ordered configurations with 0.25, 0.33, 0.66 and 0.75 compositions. These configurations were constrained to have each second layer completely filled either with Ca or with Mg, while random mixing of Ca and Mg was allowed in the other layers. In carbonates, one can also assume the existence of ordered schemes in which Ca and Mg layers alternate in sequences and ratios different from that of dolomite. Our simulation cell containing 6 layers permitted to sample only 1:5, 2:4, 3:3, 4:2 and 5:1 ratios between Ca- and Mg-filled layers, but, evidently, larger supercells would have permitted simulation of virtually any ratio. Our calculations show that such mixed-layered configurations have very low enthalpies. We have also performed impurity calculations for such mixed-layer configurations to see that these configurations represent local minima in the energy. (The introduction of the impurities caused a significant increase in the enthalpy with respect to the stoichiometric compositions). 5. Monte Carlo simulations of the temperature-dependent disorder

Fig. 2. The results of the static lattice energy calculations for (a) Ca– Mg, (b) Ca–Fe and (c) Ca–Mn binaries. The dashed and solid lines show the results of the impurity calculations. (These lines are predicted based on the enthalpies of only seven configurations). The enthalpies of all sampled configurations within 4  4  1 supercell are shown as crosses.

The supercell calculations have permitted us to predict sets of the pair ordering constants, J n , and A kl parameters of the elastic energy term, which gave best least-squares fits to the sampled SLEC energies. These parameters are given in Table 4. Fig. 3 illustrates the quality of the fit to the SLEC results for the Ca–Mg system. The estimated J n constants have been used further to predict the enthalpies of mixing at different temperatures and the temperature dependent disorder in dolomite-like compounds. For these calculations, we have created a 16  16  2 unit cell containing 3072

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313 Table 4 Ordering interaction parameters, J n (kJ/mol cations), and the elastic term coefficients, A ij (kJ/mol cations), determined by least-squares fit to the SLEC energies

309

where PAB is the probability to find AB (e.g. Ca–Mg) ˚ , which is equal to half pair at the distance of 38.486 A of the supercell size in the a direction. This probability

The ordering energies N

˚) Distance (A

J n (Ca–Mg)

J n (Ca–Fe)

J n (Ca–Mn)

1 2 3 4 5 6 7 8 9 10 11

3.8659 4.8108 6.0528 6.1716 7.3718 7.8251 8.0669 8.3325 9.1062 9.3925 9.6216

3.3946 7.8486 0.0026 7.2188 1.3940 2.1136 0.4354 0.2478 0.5092 0.0980 1.8162

2.6748 5.7996 0.3640 4.5730 0.9398 1.6920 0.2262 0.3648 0.3374 0.3324 1.1994

1.9542 4.0428 0.2878 3.2238 0.5946 1.2356 0.1368 0.2566 0.2344 0.2338 0.8210

The configuration independent term ij

A ij (Ca–Mg)

A ij (Ca–Fe)

A ij (Ca–Mn)

12 21

47.3338 31.9834

23.0868 14.8484

15.6635 10.4524

exchangeable atoms with periodic boundary conditions. Each state different in temperature and composition was sampled with 2 million Monte Carlo steps using Metropolis algorithm. It was observed that this number of steps was enough to achieve homogeneously ordered states at low temperatures. The LRO parameter was calculated using the equation of Myers (1998) Qod ¼

PAB ðsampleÞ  PAB ðrandomÞ PAB ðorderedÞ  PAB ðrandomÞ

ð7Þ

Fig. 3. The correlation between calculated (SLEC) and fitted excess enthalpies in the calcite–magnesite system. The symbols correspond to the same set of configurations as in Fig. 2a.

Fig. 4. The results of the Monte Carlo simulation of the enthalpy of mixing in (a) Ca–Mg, (b) Ca–Fe and (c) Ca–Mn binaries. The black lines show the enthalpy isotherms between 0 and 2200 8C with the step of 100 8C. The dotted lines show the enthalpy at an infinitely high temperature.

310

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

Fig. 5. The dependence of the long-range order parameter in Ca–Mg, Ca–Fe and Ca–Mn dolomites on temperature according to the results of Monte Carlo simulation.

is zero in the ordered dolomite and is equal to 0.25 in the disordered phase. This parameter is less sensitive to fluctuations in the degree of LRO within the unit cell in comparison to the usual Q od defined as the relative difference in the atomic fractions within the sublattices. Figs. 4-6 show the results of Monte Carlo simulations for Ca–Mg, Ca–Fe and Ca–Mn systems. 6. Discussion The main result of this study is that the magnitudes of ordering and mixing effects in binary carbonate systems can be predicted on the basis of a very limited set of sampled configurations. Fig. 2a,b and c show that the majority of the SLEC energies corresponding to randomly generated configurations cluster along the curves predicted with the impurity method. The hightemperature enthalpies predicted with the Monte Carlo and the impurity methods are also in good agreement: the dotted curves in Figs. 2 and 4 are nearly coincident. This shows that the Margules equation indeed works well for the high-temperature mixing enthalpy both in the diluted limit and at intermediate compositions. This is probably because the elastic energy of volume deformation comprises the largest part of the total enthalpy of mixing, while this elastic energy is well described by the subregular Margules equation. The Js-expansion in the high temperature limit also converges to the Margules polynomial. The high-temperature impurity curve thus shows the enthalpy of a completely disordered solid solution. The low-temperature impurity curves (at the intermediate range of compositions) path through the average energy of configurations

with random distribution of cations within one sublattice of the dolomite structure and with the complete order within the other sublattice. Thus, the low-temperature curves approximately limit the range of existence of the dolomite structure. The difference between the curves characterizes the total enthalpy of dolomite ordering. It is clear that the impurity calculations done on only seven anchor configurations can predict the enthalpy difference between the two very important states avoiding their random sampling. Accurate characterization of these states with the traditional random sampling approach would require a 20–30 times greater computational effort with nearly the same thermodynamic result. The present results show, however, that the low-temperature impurity curves does not necessary coincide with the minimum enthalpy in the system. Other ordering schemes, e.g. mixed-layered structures formed in the calcite–magnesite system, can be present at compositions deviating from the 50:50 ratio. The impurity curves correspond to the minimum enthalpy only if the 50:50 phase is the only ordered compound in the system. Our results permit several observations, which could be relevant not only for the carbonates, but also for other solid solutions with size mismatch. First there is positive correlation between the magnitude of the hightemperature mixing enthalpy and the temperature of the order/disorder transition in the dolomite-type compounds. This correlation becomes clearer when one observes that the magnitude of the mixing enthalpy is nearly equal to the ordering enthalpy at dolomite composition. This in turn suggests that the mixing and ordering enthalpies have the same origin: the size

Fig. 6. The enthalpy effect of disorder in Ca–Mg, Ca–Fe and Ca–Mn dolomites predicted with Monte Carlo simulation. Values are per 2 cations in the formula unit.

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

mismatch between the mixing compounds. A correlation between the mixing enthalpy and volume difference of the end-members is well known in silicate solid solutions (e.g. Geiger, 2001). Here we observe the same correlation between the size-mismatch and the ordering enthalpy. Obviously, the ordering enthalpy will coincide with the mixing enthalpy when the enthalpy of the ordered compound is equal to the mean enthalpy of the end-members. The results of assessment of mineral equilibrium data (Berman, 1988; Holland and Powell, 1998) suggest that this relationship is nearly fulfilled for dolomite. The calorimetric data of Navrotsky and Capobianco (1987), in fact, offer significantly more negative enthalpy of the dolomite phase. Modeling studies, e.g. Davidson (1994), suggest that the enthalpy of dolomite should not be that low. The present calculations suggest a slightly larger enthalpy of dolomite with the reference to the equal mixture of the end-members. Although this result is apparently wrong qualitatively (dolomite is a stable compound according to all available experimental information), the quantitative error in the relative enthalpy of ordering in dolomite does not appear to be large. It appears that if the enthalpy of ordering in dolomite was significantly larger than that predicted here, the temperature of the order/disorder transition would be also consistently larger. The temperature of 1460 8C predicted here already exceeds the value of 1200 8C of Goldsmith and Herd (1961), which is based on the observations of thermal disordering in dolomite. The transition in Ca–Mn dolomite (kutnahorite) predicted here at 500 8C is consistent with phase relations discussed by Goldsmith and Graf (1957) and Goldsmith (1983). The experimental data show that the solvus in Mn-rich system closes at 550 8C. The analogy with the Ca–Mg system suggests that the closure of the solvus occurs in the same temperature range as the order/ disorder transition. In the Ca–Fe system, we predict the ordering to ankerite at 850 8C. This result cannot be easily compared with experiment because ankerite is unstable with respect to the mixture of calcite-rich and siderite-rich phases below 900 8C (Goldsmith, 1983). We should note that our order/disorder calculations have been done under the conditions where the phase separation was artificially suppressed. The composition independent term, which presumably governs the phase separation, has not been directly included in our Monte Carlo simulations, but added to the final result. Direct comparison with the phase relations is not possible within the limits of the present method. Our calculations ignore magnetic interactions, which might be important in this system (Burton, personal communication). We

311

also ignore the effects of the excess vibrational entropy, which could stabilize intermediate compositions. The observed correlation of the enthalpies of ordering and mixing with the size-mismatch suggests that both these effects are related through the common factor, the elastic strain. The strain accumulates in the solid solution at the volume deformation stage, and releases at the intermediate composition via ordering. It is conceivable that the elastic strain concentrates within AA- and BB-type pairs representing the mixture of the end-members, while the relaxation occurs via the preferential formation of AB-type (A–O–B, A–C–B) contacts. Note that Ca–Ca and Mg–Mg pairs are either too long or too short to fit intermediate bond distances in the solid solution, while Ca–Mg pairs are consistent with the intermediate spacing. Table 4 shows that the maximum strain, which is presumably reflected in the values of A ij parameters, is the strongest in the calcite– magnesite system. The magnitudes of the ordering interactions, which on total balance the positive elastic effect, are also the largest in the Ca–Mg system. All this again suggests (see the discussion section in Vinograd et al., 2004a) that the addition of the configuration independent term to the Js expansion is an important requirement for the magnitudes and signs of the ordering interactions were predicted correctly. Reasonable values of the calculated transition temperatures (Figs. 5 and 6) confirm this statement. The importance of introducing the elastic term is also reflected in Fig. 3, which shows that the J n –A ij expansion permits an excellent fit to the SLEC energies of virtually all compositions and all ordering states. The data set is well described with the assumption that the asymmetry is fully reflected in the elastic term. Contrary to expectation, the strength of the Js does not correlate with the distance between the cations. Table 4 shows that the strongest negative interaction is that of J 4, but not J 1. J 4 acts exactly across the CO3 group. Apparently, the presence of a rigid CO3 group limits the space available for the bond relaxation: the strongly negative J 4 value thus shows that the compression or extension across the CO3 group is extremely unfavorable. The A–C–B arrangement is the best possibility to release the strain accumulated in A–C–A and B–C–B pairs. Similar effect has been observed in Ca–Mg garnets (Bosenick et al., 2000; Vinograd et al., 2004a) where the strongest J is that acting across the SiO4 group. The separation of the excess enthalpy into the configuration dependent and configuration independent terms leads to certain methodical difficulties, however. One can see in Fig. 3 that the fit does not work well in the range of 2–3 kJ/mol. This range of energies corresponds

312

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313

to the mixed-layered ordered configurations. These configurations are built of regions resembling pure calcite or magnesite and ordered dolomite and thus have significant chemical heterogeneity. The chemically different layers within the supercell should be certainly modeled with different sets of the composition dependent A ij parameters, but the present formalism does not allow this. Due to this problem configurations with marked tendency to phase separation should not be included in the fit. One can argue that the including of the configuration independent terms in Monte Carlo simulations does not permit correct determination of miscibility gaps. In fact, when the simulations are performed within a single supercell, there is no possibility to make the elastic term consistent with the compositional heterogeneity. However, the way out of this problem is to determine phase boundaries not from Monte Carlo results, but to calculate them indirectly from the free energies. The thermodynamic integration should be then performed with the configuration-dependent Hamiltonian only, while the configuration-independent terms should be added to the final free energy result. Such calculations are out of the scope of the present study. We think that the future progress in solid solution modeling will be achieved with a combination of the empirical SLEC and ab initio calculations. The impurity method could serve as a tool permitting to tune SLEC energies and make them consistent with more accurate ab initio energies. By comparing the SLEC and ab initio calculated energies of impurity configurations one can work out correction factors and adjust the values of the Js and A ij parameters. Alternatively, ab initio energies can be used as the constraints in the potential derivation procedure. When both the empirical and ab initio impurity results are shown consistent, one can proceed with large-scale SLEC and Monte Carlo simulations. Acknowledgements We thank Deutsche Forschungsgemeinschaft for support. JDG would like to thank the Government of Western Australia for funding under the Premier’s Research Fellowship program. Comments and suggestions from two anonymous reviewers were very helpful. [LMW] References Archer, T.D., Birse, S.E.A., Dove, M.T., Redfern, S.A.T., Gale, J.D., Cygan, R.T., 2003. An interatomic potential model for carbonates allowing for polarization effects. Phys. Chem. Miner. 30, 416 – 424.

Bass, J.D., 1995. Elasticity of minerals, glasses, and melts. In: Ahrens, T.J. (Ed.), Mineral Physics and Crystallography. AGU Reference Shelf, vol. 2. American Geophysical Union, pp. 45 – 63. Burton, B.P., Van de Walle, A., 2003. First-principles-based calculations of CaCO3–MgCO3 and CaCO3–MgCO3 subsolidus phase diagrams. Phys. Chem. Miner. 30, 88 – 97. Becker, U., Fernandez-Gonzalez, A., Prieto, M., Harrison, R., Putnis, A., 2000. Direct calculation of thermodynamic properties of the barite/celestite solid solution from molecular principles. Phys. Chem. Miner. 27, 291 – 300. Becker, U., Pollok, K., 2002. Molecular simulations of interfacial and thermodynamic mixing properties of the grossular–andradite garnets. Phys. Chem. Miner. 29, 291 – 300. Berman, R.G., 1988. Internally-consistent thermodynamic data for minerals in the system Na2O–K2O–CaO–MgO–FeO–Fe2O3– Al2O3–SiO2–TiO2–H2O–CO2. J. Petrol. 29, 445 – 522. Bosenick, A., Dove, M.T., Geiger, C.A., 2000. Simulation studies of pyrope–grossular solid solutions. Phys. Chem. Miner. 27, 398 – 418. Bosenick, A., Dove, M.T., Myers, E.R., Palin, E.J., Sainz-Diaz, C.I., Guiton, B.S., Warren, M.C., Craig, M.S., 2001. Computational methods for the study of energies of cation distributions: applications to cation-ordering phase transitions and solid solutions. Mineral. Mag. 65, 193 – 219. Bush, T.S., Gale, J.D., Catlow, C.R.A., Battle, P.D., 1994. Selfconsistent interatomic potentials for simulation of binary and ternary oxides. J. Mater. Chem. 4, 831 – 837. Chai, L., Navrotsky, A., Reeder, R.J., 1995. Energetics of calciumrich dolomite. Geochim. Cosmochim. Acta 59, 939 – 944. Davidson, P.M., 1994. Ternary iron, magnesium, calcium carbonates: a thermodynamic model for dolomite as an ordered derivative of calcite-structure solid solutions. Am. Mineral. 79, 332 – 339. Dove, M.T., 2001. Computer simulations of solid solutions. In: Geiger, C.A. (Ed.), Solid Solutions in Silicate and Oxide Systems of Geological Importance. EMU Notes in Mineralogy, vol. 3. Eo¨tvo¨s University Press, Budapest, pp. 225 – 249. Dove, M.T., Thayaparam, S., Heine, V., Hammonds, K.D., 1996. The phenomenon of low Al/Si ordering temperatures in aluminosilicate framework structures. Am. Mineral. 81, 349 – 362. Fisler, D.K., Gale, J.D., Cygan, R.T., 2000. A shell model for the simulation of rhombohedral carbonate minerals and their point defects. Am. Mineral. 85, 217 – 224. Gale, J.D., 1996. Empirical derivation of interatomic potentials for ionic materials. Philos. Mag., B 73, 3 – 19. Geiger, C.A., 2001. Thermodynamic mixing properties of binary oxide and silicate solid solutions determined by direct measurements: the role of strain. In: Geiger, C.A. (Ed.), Solid Solutions in Silicate and Oxide Systems of Geological Importance. EMU Notes in Mineralogy, vol. 3. Eo¨tvo¨s University Press, Budapest, pp. 71 – 100. Goldsmith, J.R., 1983. Phase relations of rhombohedral carbonates. In: Reeder, R.J. (Ed.), Carbonates: Mineralogy and Chemistry, Reviews in Mineralogy, vol. 11. Mineralogical Society of America, Washington, D.C., pp. 145 – 190. Goldsmith, J.R., Graf, D.L., 1957. The system CaO–MnO–CO2 solid solutions and decomposition relations. Geochim. Cosmochim. Acta 11, 310 – 334. Holland, T.J.B., Powell, R., 1998. An internally consistent thermodynamic data set for phases of petrological interest. J. Metamorph. Geol. 16, 309 – 343.

V.L. Vinograd et al. / Chemical Geology 225 (2006) 304–313 Myers, E.R. 1998. Al/Si ordering in silicate minerals. Ph.D. Dissertation, Newham College, Cambridge. Navrotsky, A., Capobianco, C., 1987. Enthalpies of formation of dolomite and magnesian calcite. Am. Mineral. 79, 782 – 787. Navrotsky, A., Dooley, D., Reeder, R., Brady, P., 1999. Calorimetric studies of the energetics of order–disorder in the system Mg1x Fex Ca(CO3)2. Am. Mineral. 84, 1622 – 1626. Price, G.D., Parker, S.C., Leslie, M., 1987. The lattice dynamics and thermodynamics of the Mg2SiO4 polymorphs. Phys. Chem. Miner. 15, 181 – 190. Sluiter, M.H.F., Kawazoe, Y., 2002. Prediction of mixing enthalpy of alloys. Europhys. Lett. 57, 526 – 532. Sluiter, M.H.F., Vinograd, V.L., Kawazoe, Y., 2004. Intermixing tendencies in garnets: pyrope and grossular. Phys. Rev., B 70, 184120 – 184124. Vinograd, V.L., Sluiter, M.H.F., Winkler, B., Putnis, A., Ha˚lenius, U., Gale, J.D., Becker, U., 2004a. Thermodynamics of mixing and

313

ordering in pyrope–grossular solid solution. Mineral. Mag. 68, 101 – 121. Vinograd, V.L., Sluiter, M.H.F., Winkler, B., Putnis, A., Gale, J.D., 2004b. Thermodynamics of mixing and ordering in silicates and oxides from static lattice energy and ab initio calculations. In: Warren, M., Oganov, A., Winkler, B. (Eds.), First-Principles Simulations: Perspectives and Challenges in Mineral Sciences (CECAM workshop. Lion, September 2004), Berichte aus Arbeitskreisen der DGK vol. 14. Deutsche Gesellschaft fu¨r Kristallographie, pp. 143 – 151. Warren, M.C., Dove, M.T., Myers, E.R., Bosenick, A., Palin, E.J., Sainz-Diaz, C.I., Guiton, B.S., 2001. Monte Carlo methods for the study of cation ordering in minerals. Mineral. Mag. 65, 221 – 248. Winkler, B., Dove, M.T., Leslie, M., 1991. Static lattice energy minimization and lattice dynamics calculations on aluminosilicate minerals. Am. Mineral. 76, 313 – 331.