Merlin Meheut
email: merlin.meheut@get.omp.eu

Vasileios Mavromatis

Magali Benoit

Jacques Schott

Merlin Méheut

Bασίλǫιoς Mαυρoµάτης

Magali Benoît

Calcium isotope fractionation in solids (silicates, carbonates, oxides) from Density Functional Theory: superiority of the PBEsol functional, and vibrational dragging by strong anionic groups

Keywords:

come

Calcium isotope fractionation in solids (silicates, carbonates, oxides) from Density Functional Theory: superiority of the PBEsol functional, and vibrational dragging by strong anionic groups

INTRODUCTION

Calcium is the 5

th most abundant element in Earth's crust. It is also the 6 th most important macronutrient in plants, and one of the most important elements in human physiology. Variations of its stable isotope composition, measurable by recent progresses in mass spectrometry and analytical chemistry, bring major insights into the geochemical and physiological processes involving Ca, with application to e.g. climatic reconstruction, bone cancer diagnosis or design of efficient agricultural nutrients [START_REF] Gussone | Calcium stable isotope geochemistry[END_REF]. Ca isotopes have shown their potential to unravel processes occurring at the water-soil-plant interface [START_REF] Schmitt | Calcium stable isotope geochemistry[END_REF], as well as during its transfer to oceans and within oceans [START_REF] Gussone | Calcium stable isotope geochemistry[END_REF].

To extract the information contained in those signatures, constraints on isotopic fractionation induced by individual physico-chemical processes are mandatory. Equilibrium isotope fractionation is a key parameter and can be quantified from theoretical considerations based on statistical physics and electronic structure calculations. Efficient electronic structure calculations are made possible by the Density Functional Theory (DFT). Calcium isotope equilibrium fractionation has been explored by [START_REF] Antonelli | Kinetic and equilibrium Ca isotope effects in high-T rocks and minerals[END_REF] and [START_REF] Huang | First-principles calculations of equilibrium Ca isotope fractionation: Implications for oldhamite formation and evolution of lunar magma ocean[END_REF] for high temperature silicate systems, by [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF] and Wang et al. (2017b) at low temperature, between carbonates (calcite, dolomite, magnesite) and aqueous Ca 2+ , and by [START_REF] Colla | Calcium-isotope fractionation between solution and solids with six, seven or eight oxygens bound to Ca(II)[END_REF] and [START_REF] Moynier | Calcium isotope fractionation between aqueous compounds relevant to low-temperature geochemistry[END_REF] between various dissolved Ca species. However, it is difficult to infer the reliability of those studies. For example, [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF] and Wang et al. (2017b) provide very different estimate of the calcite-dolomite equilibrium (2.1 vs 1.2 at 300K). A key aspect of validation is the comparison to experimental data (structure, vibrations, and possibly isotopes fractionation). In particular, calculated frequencies should be properly compared to experimental measurements. This requires exhaustive compilations of vibrational properties (from Raman and Infrared spectroscopies).

The materials in this study were chosen as the most exhaustively characterized among the carbonate, silicate and oxide families. They all correspond to ionic, anhydrous, compact solids. The important cases (see below) of hydrous minerals (such as ikaite or Ca-oxalates)

or dissolved species cannot be considered as they are insufficiently well characterized for their vibrational properties.

The core approximation that permits electronic structure calculations based on the density functional theory resides in the approximate functional used to describe the way that electronic energy depends on electronic density. In that respect, an issue insufficiently addressed by the above mentioned studies is the effect of the theoretical scheme, and, central to it, the choice of the density functional. The simple local density approximation (LDA) is considered as suitable for silicates (e.g. [START_REF] Tsuchiya | Phase transition in MgSiO3 perovskite in the earth's lower mantle[END_REF], but inadequate for materials with H-bonds, for which gradient-corrected approximation (GCA) are generally preferred [START_REF] Lee | Ab initio studies on highpressure phases of ice[END_REF][START_REF] Hamann | H 2 O hydrogen bonding in density-functional theory[END_REF]. The hybrid B3LYP functional has also been used by several studies to compute carbonate solids and dissolved species (e.g Rustad et al., 2010, Moynier and[START_REF] Moynier | Calcium isotope fractionation between aqueous compounds relevant to low-temperature geochemistry[END_REF]. Although not the main focus of this paper, Ca in solution is an important system to understand fractionations observed in nature, in particular for low-temperature systems, in which isotopic signatures generally result from precipitation in solution. Ca in solution seems better described (at least for its structure) by theoretical schemes including

Van der Waals corrections [START_REF] Baer | Local Aqueous Solvation Structure Around Ca2+ During Ca2+•••Cl-Pair Formation[END_REF].

Fractionation properties are occurring between phases, that can belong to different families, better described by different theoretical schemes. However, because they ultimately originate from small differences in vibrational properties, it is very important to compute those properties as consistently as possible. One (quite satisfying) possibility is to keep the same theoretical scheme for all phases. In this goal, we evaluate the efficiency of six different theoretical schemes, among those described above. Another possibility is to correct the results of a particular scheme based on experiment. The most common approach (e.g. [START_REF] Schauble | Preferential formation of 13 C -18 O bonds in carbonate minerals, estimated using first-principles lattice dynamics[END_REF] is to consider that a given theoretical scheme biases frequencies by a constant relative amount at first order, and to take a single rescaling factor determined from the comparison with experiment to correct theoretical frequencies. The efficiency of such corrections has to be evaluated.

Various solids were considered in this work (oxides: lime CaO; carbonates: calcite CaCO 3 , dolomite CaMg(CO 3 ) 2 , aragonite CaCO 3 ; silicates: grossular Ca 3 Al 2 Si 4 O 12 , diopside CaMgSi 2 O 6 ), chosen for their structural variety, and most importantly for the existence of exhaustive quantification of their vibrational properties at the experimental level. Structures with H-bonds are absent from this list, as the more suitable ones (such as ikaite or Ca-oxalates) remain partly characterized for their vibrational properties. The structural, vibrational and fractionation properties of these materials are calculated with various theoretical schemes (PZ, PBE, vdW-DF2, BLYP, BLYP+D2, PBEsol). The code that is used here (Quantum Espresso) cannot implement the B3LYP functional efficiently, due to the use of plane wave basis sets, and the efficiency of this functional was only estimated based on literature. The PBEsol functional [START_REF] Perdew | Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces[END_REF] previously showed efficient in the challenging case of the aragonite-calcite equilibrium [START_REF] Demichelis | Ab Initio Modelling of the Structure and Properties of Crystalline Calcium Carbonate[END_REF]. For the sake of a precise comparison with the work of [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF], we also computed calcite properties within the BP86 functional [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Perdew | Density-functional approximation for the correlation energy of the inhomogeneous electron gas[END_REF], that is similar to PBE or BLYP .

METHODS

2.1. The Isotopic Fractionation Factor α.

β 44/40 Ca A is the isotopic fractionation factor of the element Ca between the phase A and a perfect gas of Ca atoms, having the natural mean isotopic concentration. The isotopic fractionation factor α 44/40 Ca A-B relative to an atom Ca, between two phases A and B can be written as the ratio of the β-factors relative to this atom and to each phase separately [START_REF] Richet | A review of hydrogen, carbon, nitrogen, oxygen, and chlorine stable isotope fractionation among gaseous molecules[END_REF].

If A is a crystalline solid, β 44/40 Ca A can be computed from the harmonic vibrational properties of A using

β 44/40 Ca A = Q( 44 Ca N ) Q( 40 Ca N ) 1/N m 40 m 44 3/2 (1)
where Q( 44 Ca N ) and Q( 40 Ca N ) are the partition functions of the system having all the N Ca atoms substituted by 44 Ca and 40 Ca, respectively. The harmonic partition function of a crystalline solid is:

Q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 3N at i=1 * {q} e -hν q,i /(2kT )
1 -e -hν q,i /(kT ) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ 1/(N q )

(2)

where ν q,i are the frequencies of the phonon with wavevector q and branch index i=1,3N at .

N at is the number of atoms in the unit cell, N is the number of Ca atoms in the unit cell,

T is the absolute temperature, k is the Boltzmann constant and h is Plank's constant. The product is performed on a sufficiently large grid of N q q-vectors in the Brillouin zone. The ⋆ symbol above the product of Eq. ( 2) indicates that the three translational modes with ν 0,i = 0 are not considered.

Vibrational Analysis.

In order to identify the relative contributions of each vibrational modes to the overall βfactor, [START_REF] Méheut | Structural control over equilibrium silicon and oxygen isotopic fractionation: A first-principles density-functional theory study 258, 28 -37[END_REF] developed an approach called "vibrational analysis". In this approach, we consider the β-factor computable from a linear development in µ = 1 m :

ln β 44/40 Ca A ≈ 1 N q N q,i ∂ ln(f (ν µ q,i )) ∂µ µ=µ ∆µ, (3)
where µ = (µ 44 + µ 40 )/2 , ∆µ = µ 44 -µ 40 and

f (ν) = ν e -hν 2kT
1 -e -hν kT .

Eq. ( 3) expresses the β-factor as a product of the contributions from each vibrational mode.

Instead of running over all phonon modes, Eq. ( 2) can be limited to the Γ-point (q=0) [START_REF] Dupuis | Fractionation of silicon isotopes in liquids: The importance of configurational disorder[END_REF], giving:

β 44/40 Γ Ca A = 3N at i=3 f (ν 44,i ) f (ν 40,i ) (4)
where ν 44/40,i refers to the frequency associated with the system with all sites occupied by 44 Ca or 40 Ca isotopes. The product runs over all Transversal Optic modes.

By opposition to β Γ , the β-factor defined by Eqs. ( 1) and (2) will be further referred as "exact".

Then, we derive the equivalent of Eq. (3) for Γ-point sampling:

ln β 44/40 Γ Ca A ≈ 1 N i ∂ ln(f (ν µ i )) ∂µ µ=µ ∆µ = i h i , (5)
where we have defined the individual contribution h i (in ) of mode i to the β-factor at Γ.

Error estimation based on experimental frequencies

To quantitatively assess the error of our calculated fractionation properties based on existing experimental frequencies, we have set up an original approach. First, the experimental frequencies are usually only available for limited q-points in the First Brillouin Zone (FBZ), and for the 40 Ca material. Except for CaO, they are here originating from IR and Raman spectroscopy, and correspond to the Γ-point (q=0) of the FBZ. For these materials, we have computed the vibrational analysis only at q=0 (Eq. ( 5)), giving for each frequency its contribution (in ) to the β-factor at q=0. Limiting the sampling to the Γ-point (Eq. ( 4))

shows valid for large systems [START_REF] Dupuis | Fractionation of silicon isotopes in liquids: The importance of configurational disorder[END_REF], for Si isotopes).

In the case of the present study, Table 1 shows the β-factors calculated with the PBEsol functional, exactly (Eqs ( 1) and ( 2)) and with Γ-point sampling (Eq (4)). For grossular, lnβ Γ shows negligibly different from the exact calculation. For dolomite, calcite, diopside and aragonite, probably due to the small size of some cell parameters, the difference between lnβ exact and lnβ Γ is not negligible in absolute value (around 1 ). Given the limited relative error, we however expect the Γ sampling as adequate to evaluate the relative contributions of the different parts of the vibrational spectrum to fractionation properties. This is necessary for our correction procedure presented below.

For CaO, since we have dispersion curves giving frequencies on the whole Brillouin zone, and since considering the frequencies at q=0 only gave a result significantly different from the well converged one, we considered the frequencies on a 2x2x2 grid of the Brillouin zone (centered on zero, see Table 1).

For each frequency we can now assess what is its contribution to Ca isotope fractionation properties (h i ) and its offset-vs-experiment (χ i ). To correct for this offset and obtain the "true" β-factor, we will follow the lines exposed in Appendix B of [START_REF] Méheut | Structural control over equilibrium silicon and oxygen isotopic fractionation: A first-principles density-functional theory study 258, 28 -37[END_REF] We assume that a calculated frequency 40 ν i and its corresponding isotopic frequency shift 44/40 ∆ i = 40 ν i -44 ν i are inaccurate by the same relative amount χ i :

40 ν i,calc = χ i × 40 ν i,exp and 44/40 ∆ i,calc = χ i × 44/40 ∆ i,exp
In this case, this will induce an error χi on h i equal to:

χi ≡ h i,calc h i,exp = χ i 1 + f 2 (u i ) -1 1 -f (u i )cosh(u i /2) (6) with u i = h 40 ν i kT and f (u) = u e u/2 -e -u/2 .
As will be shown below, the frequencies contributing to Ca isotope fractionation at equilibrium extend typically from 150 to 350cm -1 . For this range, and for T>300K, equation ( 6)

gives χ very close to 2χ (more precisely, χ will be between 1.98χ -for f=150cm -1 -and 1.91χ

-for f=350cm -1 -) . We will therefore consider the following equation to correct the β-factor for frequency errors:

ln β 44/40 Γ,corr Ca = i (h i × (1 -2 * χ i )) (7)
This equation is strictly equivalent to the following one:

ln β 44/40 Γ,corr Ca = ( i h i ) × (1 -2 * χ h ) = ln β Γ × (1 -2 * χ h ), (8)
where χ h is the relative frequency error weighted by the contribution h:

χ h = ∑ i χ i * h i ∑ h i (9)
If this reasoning is correct, the error on the logarithmic β-factors is therefore twice the error χ h . χ h measures the frequency error of the frequencies important for Ca isotope fractionation properties. Practically, χ h was estimated based on existing experimental data, i.e. the sum in Eq. ( 9) was realized on the modes for which experimental frequencies are available. The accuracy of this estimate is therefore depending on the completeness of the experimental data.

In addition, if the exact β-factor is different from β Γ , and if χ h is known, the same correction can be applied to the exact β-factor:

ln β 44/40 corr Ca = ln β 44/40 exact Ca × (1 -2 * χ h ) (10)
In the literature (e.g. [START_REF] Schauble | Preferential formation of 13 C -18 O bonds in carbonate minerals, estimated using first-principles lattice dynamics[END_REF], the correction for functional error on vibrational frequencies is realized through a "scaling factor" (SF), equivalent to multiplying all frequencies by the same χ factor, that we will determine here by directly averaging the χ i :

ln β 44/40 corr-SF Ca = ln β 44/40 exact × (1 -2 * χ) (11)
with:

χ = ∑ i χ i nf req (12)
with nf req the number of available experimental frequencies.

DFT calculations

The phonon frequencies are computed from first principles using density functional theory (DFT) [START_REF] Hohenberg | [END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects 140[END_REF].

The core approximation that permits electronic structure calculations based on DFT resides in the building of an approximate functional that describes the relation between the electronic density and electronic energy.

The choice of the functional can be rooted on various considerations. The simple local density approximation (LDA), such as parametrized by Perdew and Zünger (PZ, Perdew and Zunger, 1981), is numerically very efficient. It has been shown reliable for silicate systems to predict their structural, vibrational and thermodynamical properties such as state diagrams (e.g. [START_REF] Tsuchiya | Phase transition in MgSiO3 perovskite in the earth's lower mantle[END_REF]. It was chosen by Fang Huang et al. to compute Ca isotope fractionation in anhydrous silicates and carbonates (Wang et al., 2017b;[START_REF] Huang | First-principles calculations of equilibrium Ca isotope fractionation: Implications for oldhamite formation and evolution of lunar magma ocean[END_REF]. The PZ functional is however inadequate for materials containing hydrogen.

Gradient-corrected functionals such as PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]; BLYP [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] are prefered in this case. These functionals were used in our preceeding works on Si, O and H fractionation properties [START_REF] Méheut | Structural control over equilibrium silicon and oxygen isotopic fractionation: A first-principles density-functional theory study 258, 28 -37[END_REF][START_REF] Méheut | Silicon isotope fractionation in silicate minerals: Insights from first-principles models of phyllosilicates, albite and pyrope[END_REF][START_REF] Füger | Oxygen isotope fractionation during smithsonite formation from aqueous solutions[END_REF]. On the other hand, [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF][START_REF] Colla | Calcium-isotope fractionation between solution and solids with six, seven or eight oxygens bound to Ca(II)[END_REF][START_REF] Moynier | Calcium isotope fractionation between aqueous compounds relevant to low-temperature geochemistry[END_REF] focus on the more evolved B3LYP functional [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF]. Several properties of the Ca atom appear to be very sensitive to the choice of the functional, but so far little is known about their reliability. An illustrative example of this is the number of water molecules in the coordination sphere of the Ca 2+ aquo ion as inferred from molecular dynamics simulations, that vary from 5 to 10 depending on the way the atomic interactions are taken into account (considering both electronic structure and calculations based on empirical potentials, [START_REF] Chizhik | Structure of hydration shell of calcium cation by NMR relaxation, Car-Parrinello molecular dynamics and quantum-chemical calculations[END_REF]. Experimental results (X-ray, Neutron Diffraction or EXAFS studies) are, in this regard, no more conclusive (see discussion in [START_REF] Chizhik | Structure of hydration shell of calcium cation by NMR relaxation, Car-Parrinello molecular dynamics and quantum-chemical calculations[END_REF]. Considering only methods based on first-principles calculations, the range of coordinations found for Ca

2+ (aq.)
is smaller (6-7, Baer and Mundy, 2016 and references therein), but their reliability remains questioned, and a limited number of functionals have been tested so far (BLYP, PBE, HCTH, PBE96, [START_REF] Chizhik | Structure of hydration shell of calcium cation by NMR relaxation, Car-Parrinello molecular dynamics and quantum-chemical calculations[END_REF]. [START_REF] Baer | Local Aqueous Solvation Structure Around Ca2+ During Ca2+•••Cl-Pair Formation[END_REF] argued on the validity of a calculation based on the BLYP functional in addition to the dispersion correction (D2) put forth by Grimme [START_REF] Grimme | Accurate description of van der Waals complexes by density functional theory including empirical corrections[END_REF]. This correction, intended to remedy to the incomplete account of weak dispersion forces by gradient-corrected functionals, is generally considered for aqueous systems. Theoretical schemes accounting for dispersion forces (vdW-DF2 functional, [START_REF] Lee | Higher-accuracy van der Waals density functional[END_REF] were tested to compute Si fractionation properties [START_REF] Stamm | Extreme silicon isotope fractionation due to Si organic complexation: Implications for silica biomineralization[END_REF], but did not show any net advantage in that particular case. As illustrated above, the choice of the functional is generally system dependent. However, as fractionation properties result from the thermodynamic equilibrium between various phases, it is utterly important to find a theoretical scheme that is able to properly model various kinds of systems with the same level of accuracy. As such, the efficiency of our theoretical scheme is primarily estimated as its consistency, i.e. its capacity to induce errors that are as much system-independent as possible.

In this work we have used the local approximation of Perdew and Zunger (PZ) [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF], the generalized-gradient approximations to the exchange-correlation functional of Perdew, Burke and Ernzerhof (PBE) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] and of Becke, Lee and Par (BLYP) [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF]. We have also considered calculations with the vdW-DF2 non-local functional [START_REF] Lee | Higher-accuracy van der Waals density functional[END_REF], and calculations with the BLYP functional with the Grimme-D2 correction (BLYP-GD2) [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF][START_REF] Barone | Role and effective treatment of dispersive forces in materials: Polyethylene and graphite crystals as test cases[END_REF]. Those two theoretical frameworks are considered here for their interest to treat aqueous systems [START_REF] Baer | Local Aqueous Solvation Structure Around Ca2+ During Ca2+•••Cl-Pair Formation[END_REF]Mundy, 2016, Ducher et al., 2017), as the fractionation properties of aqueous Ca are of high interest for Ca isotopes geochemistry. Here, the materials in consideration should not require such theoretical frameworks. However, as isotopic fractionation is measured between two different phases, it is of primary importance to treat both systems consistently, to avoid systematic errors. As a consequence, if the treatment of Ca in solution requires the use of BLYP-GD2 or vdW-DF2 frameworks, it is essential to check for their efficiency for crystalline solids as well. A last functional, PBEsol [START_REF] Perdew | Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces[END_REF] was also tested here.

This functional previously showed efficient as an alternative to B3LYP in the challenging case of the aragonite-calcite equilibrium [START_REF] Demichelis | Ab Initio Modelling of the Structure and Properties of Crystalline Calcium Carbonate[END_REF].

Ionic cores of each atom are described by norm-conserving pseudopotentials [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF] in the Kleinman-Bylander form [START_REF] Kleinman | Efficacious form for model pseudopotentials[END_REF]. The pseu-dopotentials used for Al, Si and O are described in the electronic annexes of [START_REF] Méheut | Equilibrium isotopic fractionation in the kaolinite, quartz, water system: predictions from first-principles density-functional theory[END_REF]. The pseudopotential used for C is described in [START_REF] Füger | Oxygen isotope fractionation during smithsonite formation from aqueous solutions[END_REF] . The pseudopotentials used for Ca and Mg were taken from the PSlibrary (Dal Corso, 2014).

Computational details for calcite are the same as in [START_REF] Füger | Oxygen isotope fractionation during smithsonite formation from aqueous solutions[END_REF] . For all other minerals (lime CaO, dolomite, aragonite, diopside, grossular), the electronic wave-functions are expanded in plane-waves up to an energy cut-off ǫ cut =80Ry, and the charge density cut-off is 4ǫ cut .

The electronic integration is performed by sampling the Brillouin zone with a 4×4×4 k-points grid for CaO, a 3×3×3 k-point grid for dolomite, a 3×2×3 k-points grid for aragonite, and a 2×3×2 k-points grid for diopside according to the Monkhorst-Pack scheme [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF]. Due to the large unit-cell of grossular, the electronic sampling can be restricted to a single k-point (1/4, 1/4, 1/4), chosen according to [START_REF] Baldereschi | Mean-Value Point in the Brillouin Zone[END_REF].

Phonon frequencies are computed using linear response theory [START_REF] Baroni | Phonons and related crystal properties from density-functional theory[END_REF], with the Quantum-Espresso package [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF]. Interatomic force-constants are obtained from the dynamical matrices computed exactly (within DFT) on a n × m × p grid of q-vectors (4×4×4, 2×3×2, 3×2×3 and 1×1×1 for CaO, diopside, aragonite and grossular, respectively)

Long-range effects are taken into account by computing Born effective-charges and static dielectric constants [START_REF] Baroni | Phonons and related crystal properties from density-functional theory[END_REF]. Dynamical matrices and thus phonon frequencies, can then be obtained in any point of the reciprocal space by Fourier-interpolation of the force constants. For all materials, the vibrational partition function (Eq. ( 2)) is converged with a 5×5×5 interpolation grid.

RESULTS

Structural properties

The relaxed structures are reported in In terms of the overall variation of the error of a given theoretical scheme from one material to the other (hereafter called spread), the PBEsol calculation shows also superior to the others (except vdW-DF2, cf Table EA-1), with a spread of 0.25%, versus 0.4% to 0.7% for the other calculations. It means that, as fas as we can judge, PBEsol (together with vdW-DF2) is more consistent, or more transferable than the other theoretical schemes.

Lastly, it should be noted that the structural parameters of CaO, calcite, dolomite and aragonite are usually more underestimated (or less overestimated), whereas the structural parameters of grossular and diopside are generally less underestimated (or more overestimated). There is however no systematic trend.

Vibrational properties

Tables EA-2 to EA-7 present, for each material considered here and the 6 considered theoretical schemes (functionals+correction), the complete vibrational properties, and compare them with experimental frequencies. For the frequencies that could be compared to their measured counterpart, we computed the offset χ i = ν th ν exp .

Figure 2, corresponding to our PBEsol calculation, shows typical patterns for the χ(ν exp )

diagram: for CaO, calcite and aragonite, χ is highly variable below 400 cm -1 , and much less for higher frequencies. In comparison, the two silicates grossular and diopside show much less scatter below 400 cm -1 . Note that, as discussed in 3.3, the 0-400cm -1 interval is associated with the largest contribution to Ca isotope fractionation properties at equilibrium.

This can be seen also based on statistical analysis (see Fig. 3). For the PBEsol functional, the variability of χ (σ(χ)) is 3.7%, 3.6% and 4.0% for CaO, calcite and aragonite, respectively, and only 1.3% and 1.7% for grossular and diopside, respectively.

Note also that for aragonite with the BLYP functional (with or without Grimme-D2 correction), two calculated frequencies are negative (B 2u (1) and A u (1), Table EA-5). Although this may appear as an argument to disqualify this functional, those frequencies contribute negligibly to the β-factor (see discussion below), so we decided to simply exclude them from our estimations.

Based on statistical analysis, Fig. 3 shows that, for a given material, the variability of the frequency error χ (σ(χ)) can be quite important, up to 12.3% for the BLYP calculation of calcite. Not surprisingly, this variability is generally smaller for the two silicates considered here, as underlined above for PBEsol. Noticeably, PBEsol shows the lowermost variability σ(χ) compared to other theoretical schemes for all materials, with a few cases where PZ or PBE show equivalent results. Interestingly, even B3LYP calculations found in the literature do not show less variability than PBEsol. Looking at the spread (total domain of variation of χ between the different materials), PBEsol shows also more consistent than other schemes except PZ. However, to conclude on the consistency of a given theoretical scheme, it is better to take into account the contribution of the various frequencies to the fractionation, as Ca fractionation is related to some particular modes only. This will be discussed below.

β-factors and correction

The results will be discussed in terms of the logarithmic β-factors (lnβ) and logarithmic fractionation factors (lnα) expressed in parts per thousand ( ). Table 2 summarizes all the fractionation properties calculated at 300K within all the theoretical frameworks considered.

Figure 4 represents all the mineral pair fractionations calculated here.

Tables EA-2 to EA-7 present, for each material considered here and the 6 considered theoretical schemes (functionals+correction), the complete vibrational properties at the Γ-point, the calculated exact logarithmic β-factors at 300K, and the contribution h i of each Γ-point vibrational frequency to the β-factor (except for CaO, for which a 2x2x2 Brillouin Zone Sampling is considered). do not depend on the theoretical scheme). The main contributing frequencies are in the range 100-400cm -1 . Note, however, that the silicates (grossular and diopside) show highfrequency contributions, if not dominating, at least much higher than for the other materials, in particular the carbonates. Those high vibrational frequencies correspond to modes (i.e., collective atomic displacements) implying the displacements of the more strongly bonded atoms (SiO, AlO). The fact that they are contributing to the fractionation of Ca isotopes means that the Ca atom is in fact displaced in those modes, even if slightly. In other words,

the Ca atom appears dragged by the high-frequency displacement of other structural units (in particular, SiO 4 or AlO 6 ), resulting in a contribution of those high-frequency modes to fractionation.

Lastly, Tables EA-2 to EA-7 also compare calculated and experimental frequencies, giving the offset χ i and the corrected β-factor at 300K following Eqs. ( 11) or ( 9). Table 2 summarize these fractionation properties (raw and corrected) at 300K.

DISCUSSION

Uncertainty on Ca fractionation

As seen on Figure 4, the calculated fractionation properties are varying significantly, typically by 1 , between the different theoretical frameworks. More precisely, the variation due to the theoretical framework, hereafter called "functional spread" goes from 0.4 for the grossular-aragonite pair (fractionation between 7.3 for BLYP and 6.9 for vdW-DF2) to 1.7 for the grossular-diopside pair (fractionation between 4.5 for PZ and 2.8 for BLYP+D2). Trying to acknowledge if some minerals are more problematic than others, one can observe that diopside appears the most (three times over four) in pairs showing the largest functional spreads. Mineral pairs involving only carbonates show on average less functional spread (from 0.5 to 0.9), but this is clearly not the case among silicates, as exemplified by the grossular-diopside pair. In any case, this variation due to theoretical framework appears significant with respect to the corresponding fractionations, and well above the level of analytical accuracy of isotopic measurements.

Trying to acknowledge which functionals behave worse, we can observe that PZ, and functionals accounting for van der Walls interactions (BLYP+D2, vdW-DF2) behave particularly badly. Indeed, from our 15 mineral pair calculations, we have 30 "extreme" calculations (i.e.

theoretical schemes representing either the larger or the smaller value for a given fractionation). These extreme values are coming from BLYP-D2 (10 times), PZ (7 times), vdW-DF2 (6 times), PBEsol (4 times), BLYP (2 times) and PBE (once).

Correction procedure

The interest of the above discussion is that it is independent of any correction. Its weakness is that we don't know what the correct fractionation could be, in order to determine which functional may approximate it best.

To improve the accuracy of those calculations, several works (e.g. [START_REF] Schauble | Preferential formation of 13 C -18 O bonds in carbonate minerals, estimated using first-principles lattice dynamics[END_REF][START_REF] Méheut | Structural control over equilibrium silicon and oxygen isotopic fractionation: A first-principles density-functional theory study 258, 28 -37[END_REF] have attempted to correct for the error on frequencies, based on experiment.

Scaling factor

The common approach assumes that a given theoretical scheme systematically underestimates or overestimates all vibrational frequencies by the same relative amount. Likewise, the average relative error on frequencies (χ, Eq. ( 12)) is computed. χ is then used to correct all frequencies, leading to Eq. ( 11) for the β-factor.

If the starting hypothesis was correct, this procedure should give the same result for any functional, as it corrects for its systematic error. Figure 6 illustrates this procedure for the β-factors of calcite. The black filled symbols represent the β-factors computed within the different theoretical frameworks. Then, from the comparison of the calculated vibrational properties with experiment, an average scaling factor is obtained (χ), for each theoretical scheme (for calcite, χ=+2.1%, -8.2%, , -6.3%, -9.1%, -7.3% and -3.4% respectively for PZ, PBE, vdW-DF2, BLYP+D2, and PBEsol, cf Table EA-3). This unique factor is then used to correct fractionation properties (Eq. ( 12), blue empty symbols on Figure 6). The outcome of this approach is disappointing: indeed, the corrected β-factors show the same dispersion as the raw ones, before correction.

Table 1 further illustrates the functional spread on β-factors for each mineral, before (first line for each mineral) and after rescaling (second line). In general, the rescaling procedure indeed reduces the spread, but with very variable efficiencies (almost nul for calcite -see Figure 6 -for which the spread goes from 3.6 to 3 , to quite high for CaO, for which the spread goes from 3.11 to 0.59). However, when mineral pair fractionations are now considered instead (lower part of Table 1), the spread always increases upon rescaling. This means that rescaling increases the uncertainty of calculated fractionations, instead of decreasing it.

It stems from the high variability of the frequency error χ, and from the very different contribution of each modes to Ca fractionation properties. Figure 2 shows typical plots of the error as a function of the frequency for different materials, for PBEsol calculations. The frequencies more important for Ca fractionation properties are essentially below 400cm -1 (see Fig 5 andsection 3.3). Their relative error with respect to experiment χ i is highly variable from one material to the other. χ i is also varying a lot between different modes, particularly for calcite, CaO and aragonite. Statistically speaking, the average rescaling factor presents a high uncertainty, as shown on Figure 3. For example, for the PBE calculation we have χ = -2.0 ± 8.2% for calcite, and χ = -8.2 ± 8.8% for aragonite. Correcting uniformly for these errors without further caution (e.g. for PBE, taking scaling factors of 1.02 for calcite and 1.082 for aragonite) leads to fractionations even more variable from one theoretical framework to the other, than without correction.

Improved procedure

To correct for theoretical errors on fractionation properties based on experiment, it is necessary to account for the relative contributions of each frequencies. The most important modes are those for which the Ca atom will move the most, as it is the condition for the Ca mass to affect the frequency. These important modes can be identified based on vibrational analysis. Fig. 5 shows vibrational analyses for the materials considered here. The important frequencies are concentrated in a narrow range of low frequencies, typically between 150 and 350cm -1 . Note that this frequency range corresponds also to largest uncertainties on χ (Figure 2).

The procedure consists in calculating a "h-weighed" scaling factor χ h (Eq. ( 9), see section 2.3), h denoting the contribution of a particular mode to ln β (Eq. ( 5)) . Table 2 and Figure 7 present the results of this procedure, (see also Figure 6 illustrating the particular case of calcite). On the β-factors, this improved correction significantly reduces the functional spread, typically around 0.5 , with larger values remaining for grossular and dolomite (≈ 1 ).

For the calculated fractionation factors relative to calcite, shown on Table 2 (and Fig. 8), it also permits an improvement of the spread by a factor of approx. 2 except for grossularcalcite and dolomite-calcite pairs.

A worse spread after correction could be due to an incomplete (and not representative) sampling of the vibrational properties of a mineral by experiment. This is particularly expected for grossular. Indeed, among the 97 different vibrational frequencies of grossular corresponding to Γ-point modes, 55 are silent (meaning that they are neither IR nor Raman-active), and therefore impossible to sample experimentally (except possibly by neutron inelastic scattering). Four other modes could not be associated with an experimental frequency, possibly due to insufficient spectroscopic intensity. In total, these 59 "orphan" modes contribute for more than 50% of the β-factor (∑ orphan h in Table EA-6). The estimate of χ h for grossular, that is based upon existing experiment, might therefore be less accurate than for minerals with a more complete set of active vibrational frequencies. Other important orphan contributions are for calcite (25% of the total), whereas it is ≈ 10% for aragonite and diopside. CaO lime and dolomite show no orphan contribution, and are therefore the better determined minerals. They should be preferentially used as reference for fractionation in mineral pairs.

On Table 3, we represent the functional spread before and after correction for each of the 15 independent mineral pairs. We can observe that (i) correction generally reduces the spread, except for calcite-dolomite, grossular-calcite, grossular-aragonite and CaO-diopside. (ii) the pairs not implying grossular or calcite show generally lower spreads after correction (the worse being 0.8 for dolomite-diopside), whereas the presence of grossular or calcite among the pair tends to increase the spread after correction.

Overall, if the improvement with the standard rescaling procedure is real, the improved procedure does not significantly reduce the functional variability. It does not permit to attain uncertainty ranges comparable to analytical uncertainty. Keeping only mineral pairs involving CaO, aragonite, diopside and dolomite (green and black pairs on Table 3), we obtain spreads of 0.3 to 0.8 after correction, which could be taken as representative of the uncertainty attached to our correction procedure.

Other origins for this variability after correction could be errors on experimental frequencies, anharmonicity (unaccounted for by theory, implicitly present in experimental frequencies),

or the difference between the displacements of the calculated modes and the real ones. In fact, looking at the contributions (h values) of similar modes in different calculations of the same material (important modes are emphasized in bold character in Tables EA-2 to EA-7 ), we can see that similar modes contribute very differently. For example, for dolomite (Table EA-4), the A u (2) mode contributes for 0.62 when calculated with BLYP, but for 0 when calculated with PZ. It means that the mode that is assimilated to the same displacement, associated with the experimental frequency at 361 cm -1 , in fact implies the displacement of the Ca atom with BLYP, but not with PZ. This likely reflects some complex effects related to diagonalization of slightly different matrices holding significantly different eigenvectors with similar eigenvalues, and puts some limits to a proper comparison with experiment.

In the following, this correction procedure will be further used to determine the efficiency of the different functionals. Efficient functionals should give similar results whether corrected or not.

Evaluation of the different functionals

Ideally, the efficiency of a functional should be estimated towards the exact isotope fractionation values (α-factors between two minerals). We do not know of accurate and reliable estimates of Ca isotope mineral-pair equilibrium fractionation. However, our corrected isotope fractionation properties may serve as estimates for those exact values. On Figure 9, we chose to estimate the error of a scheme by the difference between the raw and corrected isotope fractionation properties for this scheme. An efficient functional should change little with correction.

Figure 9 demonstrates the general superiority of the PBEsol functional in that matter (PBE being almost as effective). Its estimated error remains below 0.67 for all mineral pairs not including grossular (for which the correction is very uncertain, see 4.2.2), whereas all other schemes except PBE present errors larger than 1 in some cases. We have also considered the functional spread after correction (Table 3) as an estimate of the uncertainty of our procedure (dashed line on Figure 9). As seen on Figure 9, the error estimated for PBEsol remains below the inferred procedure uncertainty, except for grossular-dolomite and CaOdolomite isotope fractionations (for which PBEsol presents the smaller estimated error).

In a different attempt to characterize the error of a theoretical scheme, we considered the difference of a raw calculation with the average isotope fractionation value after correction (=average of the corrected isotope fractionation for a given mineral pair over the different theoretical schemes), as it could be considered a better estimate of the "exact" value. This alternative estimate of the error do not show very significantly different than the former one.

Our first approach therefore appears satisfactory to estimate the efficiency of a functional. In this approach, following Eq. ( 8), 2χ h appears naturally as an estimator of the relative error of our procedure on the ln β values (in %). For this error to cancel out on isotope fractionation factors, it should be independent on the mineral. In this regard, the better accuracy of PBEsol is a consequence of its higher consistency on errors of vibrational frequencies, as visualized on Figure 3.

Note also that, as discussed in the results section, the PBEsol functional showed superior to the other frameworks regarding (1) the absolute accuracy on structural parameters and average Ca-O distances in particular (Figure 1); (2) the amplitude of the spread on structural parameters, depending on the material (Figure 1); (3) the total χ variability for the frequencies of a given material (σ(χ), Figure 3). Considering that our correction procedure is only possible when experimental frequencies are available and exhaustive, which is far from systematic, and that our correction appears to present some uncertainty (around 0.8 , see 4.2.2), we argue that the PBEsol functional gives the most reliable results without correction, at least for materials comparable to those present in this study.

In the same order of idea, the estimated error for the PBEsol functional (red squares on Fig9) is essentially comparable to the assessed uncertainty of this error estimation (dashed line on Fig9).

Following these remarks, we recommend the raw PBEsol calculations as most reliable for estimating isotope fractionation properties. Our recommended isotope fractionation laws are listed in Table 4 .

Comparison with previous calculations

Previous works [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF]Wang et al., 2017a;Wang et al., 2017b;[START_REF] Antonelli | Kinetic and equilibrium Ca isotope effects in high-T rocks and minerals[END_REF][START_REF] Huang | First-principles calculations of equilibrium Ca isotope fractionation: Implications for oldhamite formation and evolution of lunar magma ocean[END_REF] have estimated Ca isotope fractionation properties between the same mineral pairs considered here. Table 5 compares To evaluate the role of the functional in these discrepancies, we compare in Table 5 calculations realized with the same functionals.

Regarding [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF], unfortunately, we cannot use B3LYP in our theoretical framework. Instead, we computed the calcite β-factor with the BP86 functional, that can be compared with their work ( difference, still far from the 25% observed between our BP86 calculations of calcite and the one of [START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF]. Further investigations will be necessary to clarify this matter.

Regarding the works of Fang Huang & al. (Wang et al., 2017a, Wang et al., 2017b[START_REF] Huang | First-principles calculations of equilibrium Ca isotope fractionation: Implications for oldhamite formation and evolution of lunar magma ocean[END_REF], excluding the role of the functional is quite straightforward, since we realized PZ calculations. The comparison (Table 5) shows the β-factors that we calculate with the same functional to be systematically lower than theirs, by various amounts, from 0.5-0.8 at 300K for calcite, dolomite, diopside, aragonite (by order of increasing difference), and up to 1.75 for grossular. As for the works of Rustad & co-workers, the functional cannot be claimed as the sole parameter responsible for the differences between the two calculations. To further document the discrepancies between both studies, we compared the relaxed structural parameters on Table 6. On the whole, our calculated cell parameters and Ca-O distances are lower (by ≈1%) than Fang Huang & al's. We looked at other studies using periodic boundary conditions, plane-wave basis sets and pseudopotentials (or,

for CaO, all-electron calculations) and the PZ functional. For CaO lime, calcite, dolomite and aragonite, we found several works fairly consistent with our findings. We hypothesize that the discrepancies of Fang Huang & al originate from inequivalences posed by using different pseudopotentials and other convergence parameters.

Comparison to isotope fractionations measured in experiment and Nature

The Ca isotope fractionation factor between calcite and fluid is the most well studied example. For this mineral-fluid pair however, an unambiguous isotope fractionation factor under isotopic equilibrium conditions does not exist. At present the general consensus is that no isotope fractionation occurs between calcite and the forming fluid. This argument is based on the observations by [START_REF] Fantle | Ca isotopes in carbonate sediment and pore fluid from ODP Site 807A: The Ca2+(aq)-calcite equilibrium fractionation factor and calcite recrystallization rates in Pleistocene sediments[END_REF] that no difference occurs between carbonitic nannofossil ooze and chalk and pore fluid samples from ODP Site 807A

in Ontong Java Plateau. This negligible isotope fractionation between solution and calcite at low precipitation rates was further confirmed by aquifer studies (Jacobson and Holm-den, 2008). Note however that this value has not been confirmed to date from experimental studies. To model the dependence in growth rate of experimentally defined Ca isotope fractionation between nucleated calcite and fluid, DePaolo (2011) suggest that slightly larger values (∆ 44/40 Ca eq calcite-f luid = -0.5 to -0.2 ) should be considered at equilibrium. More recently [START_REF] Oelkers | The rapid resetting of the Ca isotopic signatures of calcite at ambient temperature during its congruent dissolution, precipitation, and at equilibrium[END_REF] studied isotopic fractionation between calcite and Ca (aq) at near equilibrium conditions using a pH-jump technique, where calcite initially dissolved at pH 6.2 and precipitated at pH 7.5. They reported calcite fluid isotope fractionation factors at equilibrium of -0.8 for the precipitation phase and -1.6 for the dissolution phase of the experimental runs, arguing that the difference likely stems from the speciation of Ca in the fluid phase, which exhibits differences due to the prevailing pH. As we did not compute a dissolved Ca β-factor, we cannot compare our calculations directly to mineral-fluid experiments. In the following, to compare our calculated ∆ 44/40 Ca eq mineral-calcite with existing ∆ 44/40 Ca eq mineral-f luid , we will combine them to obtain the value of ∆ 44/40 Ca eq calcite-f luid that would be consistent if all values were representing equilibrium, and discuss this value.

For the dolomite-calcite isotope fractionation, Holmden (2009) report a difference in δ 44/40 Ca of -0.61 between dolomite and calcite in natural systems (sediments). This value compares well with our estimated -0.9 estimate of the dolomite-calcite isotope fractionation at 300K, and even more if we consider that the equilibration temperature might be larger (∆ 44/40 Ca eq dolomite-calcite = -0.5 at 125

• C, see below). Note however that [START_REF] Holmden | Ca isotope study of Ordovician dolomite, limestone, and anhydrite in the Williston Basin: Implications for subsurface dolomitization and local Ca cycling[END_REF] considered as a starting hypothesis that the dolomite-calcite Ca isotope fractionation should be null at equilibrium ( ∆ 44/40 Ca eq dolomite-calcite = 0). They chose to interpret their results as a consequence of the isotopic signature of the involved reacting fluids.

For the dolomite-fluid system, experimental works under hydrothermal conditions (125 A particularly striking feature of our results is the very high β-factor of grossular, although Ca in this mineral presents a VIII coordination (versus VI for CaO, calcite, dolomite), and its Ca-O distance is longer than in calcite, and similar to that in CaO. The vibrational analysis (reported in Fig. 5) shows that for the silicates grossular and diopside, the contributions of high frequencies, although not the main ones, are not negligible, and "make the difference" with e.g. carbonates, for which the highest frequencies contributing are around 300cm -1 . Those high-frequency modes are mainly modes implying Si-O or Al-O bonds, and the contribution of those modes to Ca isotope fractionation means that those modes "drag" the Ca atom. This coupling, or this interference of the Ca atom with those high frequency atomic displacements, seems able to explain that the β-factor of grossular is finally larger than carbonates, when all frequencies have been considered. On the contrary, carbonates do not show this high-frequency contribution, as the CO stretching modes appear uncoupled with the Ca displacements. Such a control of strong anionic groups on the isotopic fractionation of cations has already been proposed to account for Zr isotope fractionation in minerals [START_REF] Méheut | Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects[END_REF]. The question that remains to answer is why some modes are coupled whereas some others are not.

Apart from this control of strong anionic groups coupled to the atom of interest, it cannot be excluded that Ca-O distance is correlated with Ca isotope fractionation properties, for similar materials at least. On Figure 10, we show in particular that the isotope fractionation properties of silicates (diopside, grossular) and carbonates (at least dolomite and calcite) evolve on straight lines of similar slopes, as a function of Ca-O distances. Aragonite, however, does not follow the trend of dolomite and calcite, suggesting that coordination might also influence this relationship. On another hand, CaO lime, which presents a perfectly regular CaO octahedron, is much heavier than calcite and dolomite, of the same coordination but much more distorted. In this case, it is tempting to invoke a potential control of polyhedral distortion.

Although practical to predict isotope fractionation properties, the correlation with CaO distances asks for tangible explanations. In the case of dolomite versus calcite, one has to recognize that other parameters, such as those measuring polyhedral distortion, do not make sense: dolomite is less distorted than calcite but shows shorter CaO distances whereas we would have expected CaO distances to decrease with decreasing distortion, not increase.

In those conditions, and in line with our study on the structural controls of Si isotope fractionation [START_REF] Méheut | Silicon isotope fractionation in silicate minerals: Insights from first-principles models of phyllosilicates, albite and pyrope[END_REF], it is tempting here to invoke the different electronegativity of Ca versus Mg to explain the increasing CaO bond-length -and decreasing isotope fractionation properties -of dolomite versus calcite.

CONCLUSION

• In this work, we have estimated Ca isotope fractionation properties between carbonate, silicate and oxide minerals presenting accurate experimental estimates of their vibrational frequencies. Six theoretical frameworks were considered, showing results varying from single to double.

• To evaluate the efficiency of those theoretical frameworks, we tried to correct our calculated isotope fractionations, accounting for the error shown between calculated and experimental frequencies. The classical approach, assuming an identical error for all frequencies, appeared quite disappointing, in that corrected isotope fractionations did not show less variable than raw ones. We have set up an original approach for a more accurate correction. It relies on the estimate of the contribution of individual frequen-cies to isotope fractionation properties. To be efficient, it requires to have exhaustive experimental estimates for the modes contributing the most to isotope fractionation properties. The case of grossular, for which more than half of the contributive modes are silent, is most illustrative.

• At the term of this correction, PBEsol and PBE functionals appear the most efficient, gives better results (closer to experiment) for these properties, but it also gives results with a more systematic error with respect to experiment. For these reasons, we judge this theoretical scheme as the most reliable to predict the fractionation of Ca isotopes at equilibrium.

• The comparison with other estimates in the literature reveals the crucial importance of all numerical aspects of those calculations, such as electronic basis sets, or modelling approach (periodic boundary versus cluster modelling).

• This sensitivity to numerical details may seem specific to Ca isotopes. In fact, for Si or O for example, calculated properties appeared influenced only slightly by the considered numerical schemes [START_REF] Méheut | Equilibrium isotopic fractionation in the kaolinite, quartz, water system: predictions from first-principles density-functional theory[END_REF] and within the various theoretical schemes tested in this study, showing their variability. For each material, the first line gives the raw calculated β-factor, the second gives the β-factor corrected for the frequency error using the average relative error χ (Eq.( 11)), the third gives the β-factor corrected for the frequency error using Eq.( 10) . The functional spread represents the width of the interval obtained with the results of the different theoretical schemes for a given mineral, either raw or corrected. • C for calcite and within the various approximate functionals tested in this study. Black filled symbols : raw calculated properties; Blue empty symbols: properties corrected for the frequency error using the average relative error χ (Eq.( 11)). Red: properties corrected for the frequency error using Eq.( 10 Errors on calculated isotope fractionation properties (ln α 44/40 Ca), as estimated by the absolute difference between the "raw" calculation for a given mineral pair and a given theoretical scheme, and the corrected calculation for the same pair and scheme. Symbols indicate PBEsol (■), BLYP-D2 ( ), BLYP ( ), vdW-DF2 (△), PBE (▲), and PZ (◆). The estimated errors appear by increasing order of uncertainty, as measured by the functional spread after correction (dashed line, see • C for calcite and within the various approximate functionals tested in this study. Black filled symbols : raw calculated properties; Blue empty symbols: properties corrected for the frequency error using the average relative error χ (Eq.( 11)).

Red: properties corrected for the frequency error using Eq.( 10) properties corrected for the frequency error using Eq.( 10) Fig. 8: Calculated fractionation factors relative to calcite for the different materials and within the various approximate functionals tested in this study. Black: "raw" calculated properties; Red:

properties corrected for the frequency error using Eq.8 (see text)

Fig. 9: Errors on calculated isotope fractionation properties (ln α 44/40 Ca), as estimated by the absolute difference between the "raw" calculation for a given mineral pair and a given theoretical scheme, and the corrected calculation for the same pair and scheme. Symbols indicate PBEsol (■), BLYP-D2 ( ), BLYP ( ), vdW-DF2 (△), PBE (▲), and PZ (◆). The estimated errors appear by increasing order of uncertainty, as measured by the functional spread after correction (dashed line, see Table 3 and discussion in 4.2.2). See Table 2, for raw data. ) of lime CaO at Γ, X and L, and comparison to experiment. Measurements (Exp) are compiled from [START_REF] Saunderson | The lattice dynamics of calcium oxide[END_REF], [START_REF] Vijayaraghavan | Proceedings of Symposium on Neutron Inelastic Scattering[END_REF][START_REF] Rieder | Measurement and Comparative Analysis of the Second-Order Raman Spectra of the Alkaline-Earth Oxides with a NaCl Structure[END_REF] Values of lnβ and h (in ) are computed at 300K for the 44/40 Ca pair. χ is the deviation (in %) between the theoretical frequency and its experimental counterpart.

PZ PBE vdW-DF2 BLYP BLYP+GD2 PBESOL exp ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) Γ T O
291.1 -3.6 0.51 254.5 -15.7 0.39 261.6 -13.4 0.42 245.2 -18.8 0.37 260.2 -13.8 0.41 272.6 -9.7 0.45 302 Γ LO 569.4 -1.0 0.90 536.2 -6.7 0.81 537.9 -6.4 0.81 531.5 -7.6 0.80 547.5 -4.8 0.84 551.3 -4.1 0.85 575 X T A 205.3 -4.9 2.33 197.5 -8.6 2.33 197.4 -8.6 2.33 195.7 -9.4 2.36 199.0 -7.9 2.31 201.5 -6.7 [START_REF] Gillet | Raman spectroscopic studies of carbonates part I: Highpressure and high-temperature behaviour of calcite, magnesite, dolomite and aragonite[END_REF], [START_REF] Deines | Carbon isotope effects in carbonate systems1 1Associate editor[END_REF] for Raman, and [START_REF] Gillet | Thermodynamic properties and isotopic fractionation of calcite from vibrational spectroscopy of 18Osubstituted calcite[END_REF] for Infrared. Values of lnβ and h (in ) are computed at 300K for the 44/40 Ca pair. χ is the deviation (in %) between the theoretical frequency and its experimental counterpart. .

PZ: PBE vdW-DF2 BLYP BLYP+GD2 PBESOL BP Exp B3LYP (a) sym. ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) Raman E g
166.9 7.0 0.00 151.5 -2.9 0.00 160.7 3.0 0.00 151.7 -2.8 0.00 155.9 -0.1 0.00 156.9 0.6 0.00 152.8 -2.1 0.00 152.9 -2.0 E g 297.8 6.0 0.00 263.8 -6.1 0.00 277.6 -1.2 0.00 259.8 -7.5 0.00 270.8 -3.6 0.00 277.6 -1.2 0.00 265.2 -5.6 0.00 276.8 -1.5 E g 695.9 -2.0 0.00 679.7 -4.3 0.00 680.8 -4.1 0.00 676.3 -4.7 0.00 679.7 -4.3 0.00 685.0 -3.5 0.00 678.4 -4.4 0.00 705.0 -0.7 E g 1427.9 -0.4 0.00 1372.0 -4.3 0.00 1314.4 -8.3 0.00 1329.3 -7.2 0.00 1344.5 -6.2 0.00 1402.7 -2.1 0.00 1363.1 -4.9 0.00 1438.5 0.4 A 1g 1082.6 = 0.00 1045.0 -3.5 0.00 1017.4 -6.1 0.00 1023.3 -5.5 0.00 1032.0 -4.7 0.00 1063.2 -1.8 0.00 1040.2 -3.9 0.00 1081.7 -0. sym.

ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%)
Table EA-4: Modelled vibrational modes (frequencies ν in cm -1 ) for dolomite and comparison to experiment. Measurements (Exp) are compiled from [START_REF] Matas | Thermodynamic properties of carbonates at high pressures from vibrational modelling[END_REF], [START_REF] Hellwege | Zwei-phononen-absorptionsspektren und dispersion der schwingungszweige in kristallen der kalkspatstruktur[END_REF] and [START_REF] Böttcher | Characterization of inorganic and biogenic magnesian calcites by Fourier Transform infrared spectroscopy[END_REF] (see Table 6 of Valenzano et al., 2007). Values of lnβ and h (in ) are computed at 300K for the 44/40 Ca pair. χ is the deviation (in %) between the theoretical frequency and its experimental counterpart. .

PZ: PBE vdW-DF2 BLYP BLYP+GD2 PBESOL Exp B3LYP (a) sym. ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) Raman E g
187.1 6.3 0.00 169.1 -3.9 0.00 175.9 -0.0 0.00 167.5 -4.8 0.00 179.1 1.8 0.00 175.9 -0.1 0.00 177.0 E g 316.3 5.1 0.00 283.3 -5.9 0.00 292.4 -2.9 0.00 277.7 -7.7 0.00 292.6 -2.8 0.00 296.8 -1.4 0.00 295.5 E g 710.1 -1.9 0.00 692.4 -4.4 0.00 693.0 -4.3 0.00 688.3 -4.9 0.00 694.7 -4.0 0.00 698.2 -3.6 0.00 722.5 E g 1437.0 -0.5 0.00 832.3 -5.3 0.00 826.7 -6.0 0.00 818.2 -6.9 0.00 823.9 -6.3 0.00 824.9 -6.2 0.00 826.7 -5.9 0.00 877.6 A u LO 858. [START_REF] Carteret | The vibrational spectrum of CaCO3 aragonite: a combined experimental and quantummechanical investigation[END_REF]. Note that, contrary to the structure (Pmcn setting), the symmetry assignments are defined based on the Pnma spacegroup, consistently with the more general choice made in the literature. Values of lnβ and h (in ) are computed at 300K for the 44/40 Ca pair. χ is the deviation (in %) between the theoretical frequency and its experimental counterpart. sym.

PZ

PBE vdW-DF2 BLYP BLYP+GD2 PBESOL Exp B3LYP (a) sym. ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (
ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) A g
1084.8 -0.1 0.00 1047.5 -3.5 0.00 1018.6 -6.2 0.00 1027.1 -5.4 0.00 1030.7 -5.0 0.00 1065.3 -1.9 0.00 1085.5 1095.3 0.9 A g 1455.7 -0.5 0.03 1406.4 -3.9 0.04 1350. -7.7 0.04 1367.5 -6.5 0.04 1377.4 -5.9 0.04 1433.9 -2.0 0.03 1463.0 1473.9 0.8 sym. sym. 327.9 -0.6 0.14 309.1 -6.3 0.10 313.5 -5.0 0.13 307.5 -6.8 0.11 313.5 -5.0 0.08 317.1 -3.9 0.11 330 332.9 0.9 T 2g 349.1 0.0 0.00 331.9 -4.9 0.00 336.4 -3.6 0.00 331.7 -5.0 0.00 340.2 -2.5 0.00 338.9 -2.9 0.26 349 355.9 2.0 T 2g 387.3 1.1 0.19 362.0 -5.5 0.11 366.9 -4.2 0.16 358.9 -6.3 0.16 368.1 -3.9 0.15 374.0 -2.3 0.17 383 390.0 1.8 T 2g

B 1g 120
ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (
ν χ h ν χ h ν χ h ν χ h ν χ h ν χ h ν ν χ (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (%) ( ) (
471.3 -1.4 0.02 448.5 -6.2 0.01 460.5 -3.7 0.01 450.0 -5.9 0.01 460. 529.4 539.2 1.9 T 1u TO 537.9 -0.3 0.15 514.0 -4.8 0.12 515.8 -4.4 0.13 510.2 -5.5 0.12 523.7 -3.0 0.14 524.9 -2.8 0.14 539.8 546.4 1.2 T 1u 558.9 -3.0 576.2 583.1 1.2 T 1u TO 609.4 -1.4 0.09 587.3 -5.0 0.10 594.8 -3.7 0.08 587.4 -5.0 0.09 598.1 -3.2 0.09 596.1 -3.5 0.10 618.0 625.9 1.3 T 1u 623.8 634.4 1.7 T 1u TO 835.4 -0.3 0.01 793.2 -5.3 0.01 778.8 -7.0 0.01 780.9 -6.8 0.01 797.8 -4.8 0.01 814.9 -2.7 0.01 837.7 829.6 -1.0 T 1u LO 849.0 -0.3 806.7 -5.2 793.5 -6.8 795.0 -6.6 812.2 -4.6 828.3 -2.7 851.4 845.6 -0.7 T 1u TO 854.1 -0.3 0.02 811.2 -5.3 0.02 797.4 -6.9 0.03 799. 1 -6.7 0.03 816.6 -4.6 [START_REF] Etchepare | Amorphous materials international conference on noncrystalline solids[END_REF], [START_REF] Swamy | Hightemperature Raman spectroscopy and quasi-harmonic lattice dynamic simulation of diopside[END_REF], [START_REF] Chopelas | Spectroscopic evidence for pressure-induced phase transitions in diopside[END_REF] and [START_REF] Prencipe | Simulation of vibrational spectra of crystals by ab initio calculations: an invaluable aid in the assignment and interpretation of the Raman signals. The case of jadeite NaAlSi 2 O 6[END_REF] for Raman, and are taken from [START_REF] Zulumyan | Study of calculation of the vibrational spectrum of a crystal with complex polyatomic anions[END_REF] 652.5 -2.2 0.29 626.9 -6.0 0.25 635.9 -4.7 0.27 627.0 -6.0 0.25 641.6 -3.8 0.29 636.7 -4.5 0.27 671 A g 849.0 -0.8 0.00 813.9 -4.9 0.00 794.1 -7.2 0.00 810.4 -5.3 0.00 817.1 -4.5 0.00 830.6 -3.0 0.00 850 A g 998.0 -1.6 0.10 958.1 -5.5 0.10 956.6 -5.7 0.10 950.3 -6.3 0.10 966.4 -4.7 0.10 976.5 -3.7 0.10 1007 A g 1031.7 -1.6 0.02 988.0 -5.7 0.02 985.3 -6.0 0.03 976.5 -6.8 0.02 992.4 -5.3 0.02 1009.0 -3.7 0.02 1038

Fig 5

 5 Fig 5 visualizes these contributions h as a function of theoretical frequencies, and their integrated value (for the PBEsol calculations only, since thhe general features of these curves

  the fractionations obtained by these authors for the dolomite-calcite, aragonite-calcite, grossular-diopside and grossularlime mineral pairs. Our results show a fair agreement with Antonelli et al. (2019), and with the aragonite-calcite and dolomite-calcite calculations by Wang et al. (2017b), but our results significantly differ from Rustad et al. (2010) for dolomite-calcite (-0.9 vs -2.3 ), and with Huang et al. (2019) for grossular-diopside (4.2 vs 5.8 ). Our calculations differ from these authors by the choice of functional: Rustad et al. (2010) used B3LYP, whereas Huang et al. (2019) uses PZ.

  BLYP and BLYP+GD2 being approximately equivalent, and PZ and vdW-DF2 being significantly worse. From the point of view of several other criteria, PBEsol appears superior for the considered materials: (1) the consistency of the frequency error χ i for different modes of a given material, (2) the consistency of the average frequency error χ, or of the average error weighted by the contribution to the β-factor χ h for different minerals, (3) the amplitude of this average frequency error (4) the amplitude and consistency of the error on average CaO distances. Said differently, not only PBEsol

  (a) CaO: Speziale et al. (2006) (b) calcite: Effenberger et al. (1981) X-ray diffraction at ambient temperature (c) dolomite: Althoff (1977) X-ray diffraction (d) aragonite: Caspi et al. (2005) Synchrotron Hi-Res X-Ray (e) grossular: Geiger and Armbruster (1997) X-ray diffraction at 293K (f ) diopside: Cameron et al. (1973) X-ray diffraction at ambient temperature

- 1

 1 (in %) between calculation and experimental frequencies. Error bars represent the χ±σ(χ) uncertainty interval for χ, whereas symbols indicate χ h , for CaO lime ( ), calcite (+), dolomite(▲),aragonite (■), grossular (◆), and diopside (⋆). The four B3LYP uncertainty intervals correspond to calcite, dolomite, aragonite and grossular, from top to bottom. See Tables EA-2 to EA-7 for raw data. Note the smaller uncertainty on χ for diopside and grossular. Note also the good efficiency of the PBEsol calculation, giving smaller uncertainty on χ (smaller σ(χ) for a given material), and smaller variability on χ h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4 Fractionation calculated at 300K for all the 15 mineral pairs considered in this work. The pair order is chosen (between A-B and B-A) to give a positive fractionation in PBEsol. The pairs are ordered by increasing value of the fractionation calculated in PBEsol. . . . . . . . . . . . . . . . . . . . . . . . 48 5 Plot of the contribution h (bars, in ) of individual frequencies to lnβ Γ at 300K, for the PBEsol calculations. Dashed lines are integrated contributions, converging towards lnβ Γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Logarithmic β-factors calculated at 25

Figures 894 895 Fig. 1 :

 8951 Figures 894

Fig. 2 :

 2 Fig. 2: Theoretical frequency offset χ (in %) as a function of experimental frequency, for calculations with the PBEsol functional. As seen on Fig.5, the range of frequencies most important for Ca isotope fractionation properties are between 100 and 400 cm -1 (double arrow) .

Fig. 3 :Fig. 4 :

 34 Fig. 3: Frequency offset χ = ν calc ν exp -1 (in %) between calculation and experimental frequencies. Error bars represent the χ ± σ(χ) uncertainty interval for χ, whereas symbols indicate χ h , for CaO lime ( ), calcite (+), dolomite(▲),aragonite (■), grossular (◆), and diopside (⋆). The four B3LYP uncertainty intervals correspond to calcite, dolomite, aragonite and grossular, from top to bottom. See Tables EA-2 to EA-7 for raw data. Note the smaller uncertainty on χ for diopside and grossular. Note also the good efficiency of the PBEsol calculation, giving smaller uncertainty on χ (smaller σ(χ) for a given material), and smaller variability on χ h .

Fig. 7 :

 7 Fig. 7: Logarithmic β-factors calculated for the different materials and within the various approximate functionals tested in this study. Black filled symbols : raw calculated properties; Red:

Fig. 10 :

 10 Fig. 10: Calculated isotope fractionation properties (logarithmic β-factors at 300K) as a function of mean CaO distance for the various minerals considered here. Coordination is given in Roman number.

  (a) CaO: Speziale et al. (2006) (b) calcite: Effenberger et al. (1981) X-ray diffraction at ambient temperature (c) dolomite: Althoff (1977) X-ray diffraction (d) aragonite: Caspi et al. (2005) Synchrotron Hi-Res X-Ray (e) grossular: Geiger and Armbruster (1997) X-ray diffraction at 293K (f ) diopside: Cameron et al. (1973) X-ray diffraction at ambient temperature (g) CaO: Calculated from cell volumes at 0K extracted from Figure 4 of Erba et al. (2015) (h) calcite: Valenzano et al. (2007), BSD basis set (i) dolomite: Valenzano et al. (2007), BSD basis set (h) aragonite: Carteret et al. (2013) (j) grossular: Zicovich-Wilson et al. (2008) (k) diopside: Prencipe (2012)

Table EA

 EA

	on average Ca-O distances < Ca -O >.
	It shows that the PZ functional tends to underestimate CaO distances (and cell parameters,
	cf Table EA-1), whereas PBE and BLYP functionals tend to overestimate them, consistently
	with previous studies (for PZ functional: see e.g. Huang et al., 2019; for PBE and BLYP
	functionals: see e.g. Méheut and Schauble, 2014, Dupuis et al., 2015). The vdW-DF2 and
	BLYP+GD2 calculations, that have been implemented to effectively correct for the short-
	comings of GGA functionals in accounting for the weak Van der Waals forces (Kristyán
	and Pulay, 1994), show contrasting results. Adding the Grimme-D2 correction to the BLYP
	functional (BLYP+GD2 scheme) leads to smaller (and closer to experiment) structural pa-
	rameters, which seems reasonable considering that the Grimme-D2 correction corresponds
	to adding an attractive force. On the contrary, calculations with the vdW-DF2 functional
	lead to elongate Ca-O distances (+1.7±0.3%) and cell parameters, in particular relative
	to PBE (+1.0±0.4% for < Ca -O >). This is contradictory with our results in Stamm
	et al., 2020, which were dealing with organic materials containing Si-O bonds. In that case,
	the vdW-DF2 lead to cell parameters smaller than PBE, and closer to experiment, whereas
	Si-O bond lengths were the same with PBE and vdW-DF2, as expected (Supplemental Table
	A.1 of Stamm et al., 2020). For the materials considered in this study, the improvement of
	vdW-DF2 over PBE is not clear. Lastly, in terms of the general agreement with experiment,
	the PBEsol calculation shows the best result, with a difference of -0.45 ± 0.25% relative to
	measured bond lengths.
	-1, for each material and within each of
	the 6 theoretical schemes considered here. Average Ca-O distances are also reported. Since
	those distances are often correlated with equilibrium fractionation properties, an efficient
	way to comment on the efficiency of the various theoretical schemes is to look a their error
	on average Ca-O distances. On Figure 1, we report the relative error of each calculations

  Table5). Whereas their BP86 β calcite calculation is very similar to their reference B3LYP one (16.01 vs 16.18 ), our BP86 calculation yields a much lower value (12.26 ), similar to our PBE result (12.05 ). This suggests that other details of[START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF] calculations may have significant effects on the final result.These calculations differ from ours by two aspects: first,[START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF] use gaussian basis sets for electronic wave-function description, and second, they use cluster-like models when we use periodic boundary conditions.[START_REF] Rustad | Isotopic fractionation of Mg 2+ (aq) , Ca 2+ (aq) , and Fe 2+ (aq) with carbonate minerals[END_REF] did not see any significant effect of the choice of gaussian basis set (6-31G* vs 6-311++G(2d,2p) on the mineral β-factors. On another hand, based on a periodic boundary calculation with Gaussian basis sets,Valenzano et al., 2007 and De La Pierre et al., 2016 report different values for the Eu(2)TO frequency determined experimentally at 223 cm -1 : respectively 219.6 [-1.5%] and 225.2 [+1.0%]. In terms of its contribution to isotope fractionation, this makes a 5%

  For the aragonite-fluid system mineral growth experiments at low degrees of saturation of the fluid do not exist in the literature. Both the experimental works by[START_REF] Gussone | Model for kinetic effects on calcium isotope fractionation (δ 44 Ca) in inorganic aragonite and cultured planktonic foraminifera[END_REF] and AlKhatib and Eisenhauer (2017) exhibit Γ aragonite values >25. In both these studies however aragonite is enriched in the lighter 40 Ca isotope, by approximately -1.7 at the lower precipitation rates, closest to equilibrium. At the same temperature, we calculate:

	∆	44/40 Ca	eq aragonite-calcite = -3 , therefore ∆	44/40 Ca	eq calcite-f luid = -1.3 would be necessary
	to reconcile our calculations with those estimates. This is surprisingly similar to what we
	find based on dolomite-fluid estimates of Perez-Fernandez et al. (2017), and suggests again
	that those experiments might not be fully equilibrated.
							• C)
	report a ∆	44/40 Ca dolomite-f luid = -1.2 during dissolution, under conditions interpreted
	as close to equilibrium (Perez-Fernandez et al., 2017). At the same temperature, we find
	∆	44/40 Ca	eq dolomite-calcite = -0.5 . Our calculation and Perez-Fernandez et al.'s study would
	be consistent for ∆	44/40 Ca	eq calcite-f luid = -0.7 at 125	• C. Considering that ∆	44/40 Ca	eq values
	are essentially proportional to 1 T 2 , they should be twice as large at 300K compared to 125	• C
	( 400 2 300 2 =1.8). This would translate into ∆	44/40 Ca	eq calcite-f luid = -1.3 at 300K. This value
	seems incompatible with the general consensus (∆	44/40 Ca	eq calcite-f luid = 0 ), suggesting that
	Perez-Fernandez et al's experiments might not be fully equilibrated. Only the larger estimate
	of Oelkers et al. (2019) (∆	44/40 Ca

eq calcite-f luid = -1.6 ) would be consistent.

4.6. Controlling parameters of Ca isotope fractionation: anionic group dragging

In Table

7

, we compare the calculated β-factors with various structural properties. The two parameters that have been the most discussed in the literature are coordination (e.g.

[START_REF] Colla | Calcium-isotope fractionation between solution and solids with six, seven or eight oxygens bound to Ca(II)[END_REF]

and CaO distances (e.g.

[START_REF] Huang | First-principles calculations of equilibrium Ca isotope fractionation: Implications for oldhamite formation and evolution of lunar magma ocean[END_REF]

. Figure

10

reports our calculated isotope fractionation properties as a function of Ca-O distance.

Table 1 :

 1 • In the systems considered, and independently of CaO bond lengths, the fractionation of Ca isotopes appears strongly affected by anionic group dragging, i.e. the presence of strong anionic groups (such as SiO 4 or AlO 6 groups in silicates) vibrating at high frequency and which atomic displacements are "dragging" the Ca atom. This was previously observed for predicted Zr isotope fractionation properties[START_REF] Méheut | Drivers of zirconium isotope fractionation in Zr-bearing phases and melts: The roles of vibrational, nuclear field shift and diffusive effects[END_REF]. β-factors (in ) calculated at 300K with the PBEsol functional and various samplings of the Brillouin zone. mineral lnβ exact lnβ Γ lnβ 2×2×2

	: PZ vs PBE for oxygen isotopes;

Table 2 :

 2 Fractionation properties (raw and corrected) calculated at 300K for the different materials

Table 3 :

 3 Functional spread obtained for each mineral pair before/after correction. Colors emphasize minerals for which correction is assumed less precise due to inhomogeneous experimental sampling of their vibrational properties. green: mineral fully characterized (CaO, dolomite). Orange: calcite,

	PZ	PBE	vdW-DF2 BLYP BLYP+D2 PBEsol spread
		ln β mineral ( )
	CaO 16.97 14.38 14.56 13.86 14.96 15.60 3.11

Table 4 :

 4 Recommended fits of 1000ln β 44/40 Ca mineral and of 1000ln α 44/40 Ca mineral-calcite for 0-1200• C, based on ax 2 +bx 3 with x= 10 2 T (K) . We used the raw PBEsol results (see text for discussion).

	System Fit parameters
			a	b
	1000ln β	44/40 Ca mineral
	calcite 127.9	-23.9
	dolomite 119.5	-23.1
	aragonite 99.1	-16.4
	CaO 148.2	-23.6
	grossular 171.2	-38.1
	diopside 129.8	-26.3
	1000ln α	44/40 Ca mineral-calcite
	dolomite -8.4	0.8
	aragonite -28.8	7.5
	CaO 20.3	0.3
	grossular 43.3	-14.2
	diopside 1.9	-2.4

Table 5 :

 5 Comparison of the isotope fractionation properties calculated here with literature. R10:

	Rustad et al. (2010); FH: Fang Huang et al. (Wang et al., 2017a; Wang et al., 2017b; Huang
	et al., 2019); A19: Antonelli et al. (2019)
	property	this work R10 FH A19
			preferred isotope fractionation values at 300K
	ln α	44/40 Ca dolomite-calcite	-0.9 -2.3 -1.1
	ln α	44/40 Ca aragonite-calcite	-2.9	-2.8
	ln α	44/40 Ca grossular-diopside	4.2	5.8 4.4
	ln α	44/40 Ca grossular-lime	2.0	2.2
	isotope fractionations at 300K calculated with BP86
	ln β calcite	12.26 16.01
	isotope fractionations at 300K calculated with PZ
	ln α	44/40 Ca dolomite-calcite	-1.2	-1.1
	ln α	44/40 Ca aragonite-calcite	-3.1	-2.8
	ln α	44/40 Ca grossular-diopside	4.5	5.8
	ln β	44/40 Ca calcite	15.34	15.81
	ln β	44/40 Ca dolomite	14.14	14.74
	ln β	44/40 Ca aragonite	12.21	13.01
	ln β	44/40 Ca diopside	14.86	15.51
	ln β	44/40 Ca grossular	19.36	21.31

Table 6 :

 6 Lattice parameters (in Å and

	• ) and aver-

Table 7 :

 7 Comparison of the ln β 44/40 Ca calculated at 300K (PBEsol functional) with structural parameters of the different structures (estimated on experimental structures). Ca-O: average Ca-O distances. Vol: volume of the coordination polyhedra. D: Baur's distortion index on distances. Ca -O > th ) and measurements (< Ca -O > exp ), for

	Qel: quadratic elongation of the coordination polyhedra. BAV: bond angle variance. Coord:
	coordination.
	mineral	Ca-O	Vol	D	Qel BAV coord ln β
		Å	Å3			•2
	calcite	2.3574 17.414	0 1.0020 7.4660 6 13.33
	dolomite 2.3815 17.960	0 1.0017 6.1194 6 12.43
	CaO lime 2.4075 18.605		1.0000	0 6 15.60
	aragonite 2.5283 29.808 0.0258		9 10.41
	grossular 2.4019 23.765 0.0335		8 17.61
	diopside 2.4976 25.757 0.0567		8 13.45

-1 (double arrow) . . . . . . . . . . . . . . . . . . . . . . . . 46 3 Frequency offset χ = ν calc

  ) . . . . . . . . . . . . . . . . . . . . . . 50 7 Logarithmic β-factors calculated for the different materials and within the various approximate functionals tested in this study. Black filled symbols : raw calculated properties; Red: properties corrected for the frequency error using Eq.(10) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 8 Calculated fractionation factors relative to calcite for the different materials and within the various approximate functionals tested in this study. Black:

"raw" calculated properties; Red: properties corrected for the frequency error using Eq.8 (see text) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 9

  , for the materials considered here. The numbers in parenthesis refer to uncertainties on the last significant digit. For each distance, ∆(%) is the relative difference between the experimental and calculated value.

				Table EA-1: Experimental (exp, from literature) and calculated lat-
				tice parameters and Ca-O distances (this work, and from literature for B3LYP calculation) in Å and
						lime CaO
		Exp	(a-f )	PZ:	PBE	vdW-DF2	BLYP	BLYP+GD2	PBESOL	B3LYP	(g-k)
				value ∆(%) value ∆(%) value ∆(%) value ∆(%) value ∆(%) value ∆(%)
	a	4.8115[5]	4.7157	4.8382	4.8822	4.8781	4.8180	4.7746	4.858
	Ca-O 2.4058[3]	2.358 -2.0	2.419 +0.6 2.441 +1.5 2.439 +1.4 2.409 +0.2 2.387 -0.7	2.429 +1.0
						calcite
	a	4.9896[2]	4.9614 -0.6	5.0602 +1.4 5.0930 +2.1 5.0929 +2.1 5.0643 +1.5 5.0099 +0.4 5.0373 +1.0
	c	17.0610[11] 16.4176 -3.8 17.2059 +0.8 17.1122 +0.3 17.3614 +1.8 17.0348 -0.2 16.8516 -1.2 17.3304 +1.6
	Ca-O 2.3598[6]	2.3121 -2.0	2.3871 +1.2 2.3941 +1.5 2.4072 +2.0 2.3792 0.8	2.3504 -0.4	2.3907 +1.3
						dolomite
	a	4.8033[9]	4.7714 -0.7	4.8662 +1.3 4.8894 +1.8 4.8943 +1.9 4.8626 +1.2 4.8216 +0.4 4.8376 +0.7
	c	15.984[4]	15.5651 -2.6 16.2415 +1.6 16.2287 +1.5 16.4176 +2.7 16.0957 +0.7 15.9324 -0.3 16.2756 +1.8
	Ca-O 2.378[1]	2.3315 -2.0	2.4049 +1.1 2.4167 +1.6 2.4272 +2.1 2.4001 +0.9 2.3682 -0.4	2.4099 +1.3
						aragonite
	a	4.96183[1] 4.8863 -1.5	5.0173 +1.1 5.0543 +1.9 5.0581 +1.9 5.0677 +2.1 4.9517 -0.2	5.008 +0.9
	b	7.96914[2] 7.8346 -1.7	8.0433 +0.9 8.0513 +1.0 8.1074 +1.7 8.0516 +1.0 7.9406 -0.4	8.029 +0.8
	c	5.74285[2] 5.5381 -3.6	5.8239 +1.4 5.8540 +1.9 5.9162 +3.0 5.8192 +1.3 5.6768 -1.2	5.861 +2.1
	Ca-O 2.5272 [5] 2.4692 -2.3	2.5562 +1.1 2.5689 +1.7 2.5838 +2.2 2.5679 +1.6 2.5120 -0.6	2.560 +1.3
						grossular
	a	11.847[1]	11.7000 -1.2 11.9923 +1.2 12.0535 +1.7 12.0651 +1.8 11.9986 +1.3 11.8471 =	11.9368 +0.8
	Ca-O 2.405[1]	2.3711 -1.4	2.4332 +1.2 2.4530 +2.0 2.4538 +2.0 2.4383 +1.4 2.4007 -0.2	2.4180 +0.5
						diopside
	a	9.746[4]	9.6478 -1.0	9.9031 +1.6 9.9232 +1.8 10.0017 +2.6 9.8583 +1.2 9.7752 +0.3 9.8931 +1.5
	b	8.899[5]	8.7739 -1.4	9.0448 +1.6 9.0594 +1.8 9.1242 +2.5 9.0183 +1.3 8.9151 +0.2 9.0199 +1.4
	c	5.251[6]	5.1996 -1.0	5.3375 +1.6 5.3430 +1.8 5.3926 +2.7 5.3254 +1.4 5.2661 +0.3 5.3265 +1.4
	β	105.63[6]	105.77	106.46	105.86	107.01	105.68	106.05	106.37
	Ca-O 2.4976	2.4547 -1.7	2.5335 +1.4 2.5418 +1.8 2.5662 +2.7 2.5325 +1.4 2.4930 -0.2	2.5350 +1.5
	∆ CaO			-1.85±0.45	+1±0.4	+1.75±0.25	+2.05±0.65	+0.9±0.7	-0.45±0.25	+1±0.5
	interval
	(%)

•

Table EA -

 EA 2: Modelled vibrational frequencies ν (in cm-1

  Table EA-3: Modelled vibrational modes (frequencies ν in cm -1 and intensities) for calcite, and comparison to experiment. PBE refers to the present calculations. Measurements (Exp) are compiled from

									2.33 216
	X LA	311.1 5.8 2.90 276.6 -5.9 2.17 284.3 -3.3 2.36 269.7 -8.3 2.04 288.1 -2.0 2.43 292.9 -0.4 2.49 294
	X T O	315.6 0.5 0.98 283.0 -9.9 0.39 293.0 -6.7 0.44 277.5 -11.6 0.23 292.5 -6.8 0.56 298.0 -5.1 0.68 314
	X LO	393.2 -1.7 0.27 381.9 -4.5 0.52 387.6 -3.1 0.42 383.8 -4.1 0.58 386.2 -3.5 0.38 385.7 -3.6 0.42 400
	L T A	189.0 -3.1 3.07 167.2 -14.2 2.41 168.9 -13.4 2.46 161.4 -17.2 2.25 172.1 -11.8 2.55 178.0 -8.7 2.73 195
	L LA	382.5 3.9 6.06 361.5 -1.8 5.43 359.9 -2.2 5.39 356.9 -3.0 5.31 365.3 -0.7 5.54 371.6 1.0 5.73 368
	L T O	323.9	0.00 289.6	0.00 295.2	0.00 283.5	0.00 298.0	0.00 305.7	0.00
	L LO	485.7 2.7 0.00 457.6 -3.3 0.00 461.7 -2.3 0.00 456.0 -3.6 0.00 468.8 -0.9 0.00 469.2 -0.8 0.00 473
	∑ h = ln β 2×2×2	17.02	14.45	14.63	13.94	15.02	15.68
	ln β exact	16.97	14.38	14.56	13.86	14.96	15.60
	χ				-0.2	-7.8	-6.6	-9.3	-5.8	-4.2
	σ(χ)	3.7	4.7	4.4	5.7	4.7	3.7
	χ	h †		1.0	-6.5	-6.0	-8.0	-4.8	-3.0
	†χ	h	=	∑ χ×h ∑ h

  [START_REF] Zicovich-Wilson | Ab initio simulation of the IR spectra of pyrope, grossular, and andradite[END_REF] for infrared,[START_REF] Dovesi | Ab initio quantum-mechanical simulation of the Raman spectrum of grossular[END_REF] for Raman Table EA-7: Modelled vibrational modes (frequencies and intensities) for diopside, and comparison to experiment. Measurements (Exp) correspond to typical values compiled from

												Table EA-6: (continued)
						PZ:			PBE		vdW-DF2		BLYP		BLYP+GD2		PBESOL		Exp. B3LYP	(a)
	sym.	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	ν	χ
						(%) ( )		(%) ( )		(%) ( )		(%) ( )		(%) (	)		(%) (	)	(%)
	χ					-0.7			-6.4			-5.4			-7.1			-5.0			-3.9		0.5
	σ(χ)		1.3			1.8			2.1			2.3			2.4			1.3		1.5
	χ	h †			-0.1			-7.0			-5.5			-7.8			-5.7			-3.9
	†χ	h	=	∑ χ×h ∑ h
	(a) Zicovich-
																		0.03 833.2 -2.7 0.02 856.3 849.6 -0.8
	T 1u LO 882.4 -0.3		839.0 -5.2		824.6 -6.8		826.7 -6.6		844.2 -4.6		861.0 -2.7		884.7 879.8 -0.6
	T 1u TO 902.8 -0.5 0.08 860.9 -5.2 0.08 846.3 -6.8 0.09 848.9 -6.5 0.09 864.9 -4.7 0.08 882.3 -2.8 0.08 907.8 901.8 -0.7
	T 1u LO 991.2 -1.6		948.5 -5.8		940.8 -6.6		939.7 -6.7		954.3 -5.3		969.0 -3.8		1007.4 1002.5 -0.5
														Silent
	T 2u		162.1		0.39 147.1		0.37 147.7	0.38 140.4	0.36 144.6	0.37 153.8	0.38
	T 2u		191.7		0.39 182.1		0.32 184.3	0.33 182.1	0.31 180.2	0.34 185.4	0.34
	T 2u		249.6		0.53 232.0		0.46 236.8	0.45 230.7	0.45 236.0	0.43 238.9	0.48
	T 2u		295.9		0.14 283.6		0.10 291.2		0.09 283.9	0.07 293.0	0.06 287.8	0.13
	T 2u		314.9		0.29 290.6		0.21 298.1	0.22 289.1	0.20 297.1	0.20 300.6	0.24
	T 2u		329.0		0.06 313.0		0.02 315.0		0.01 311.5	0.00 321.9	0.02 319.6	0.04
	T 2u		354.0		0.30 327.3		0.22 330.3	0.17 321.0	0.17 333.7	0.17 339.5	0.27
	T 2u		401.3		0.20 377.0		0.16 377.8	0.19 372.3	0.17 383.0	0.15 387.9	0.17
	T 2u		439.6		0.03 412.0		0.04 409.2		0.08 405.3	0.06 422.0	0.06 425.2	0.03
	T 2u		465.7		0.07 437.2		0.09 427.5		0.11 426.7	0.11 446.9	0.12 451.7	0.07
	T 2u		514.5		0.06 490.9		0.05 500.0		0.05 492.3	0.05 502.4	0.06 499.8	0.06
	T 2u		587.0		0.07 560.4		0.05 565.2		0.04 558.5	0.04 571.5	0.05 571.8	0.06
	T 2u		600.9		0.09 577.5		0.09 589.3		0.10 580.9	0.10 592.2	0.10 586.0	0.09
	T 2u		826.1		0.02 788.8		0.02 777.5		0.02 778.2	0.02 794.6	0.02 807.7	0.02
	T 2u		869.6		0.02 829.2		0.02 822.5		0.02 820.2	0.02 837.5	0.02 848.9	0.02
	T 2u		889.6		0.06 847.6		0.06 839.3		0.07 838.0	0.07 856.3	0.07 868.1	0.06
	T 1g		185.0		0.63 171.0		0.45 171.8	0.53 166.6	0.57 173.5	0.62 177.6	0.63
	T 1g		194.3		0.69 174.4		0.55 178.9	0.56 170.7	0.45 175.7	0.51 183.9	0.63
	T 1g		235.6		0.24 214.5		0.19 213.0	0.26 208.6	0.23 215.0	0.20 224.6	0.20
	T 1g		288.0		0.13 271.9		0.03 276.0		0.04 269.6	0.05 277.2	0.10 278.6	0.08
	T 1g		300.6		0.53 279.3		0.58 284.7	0.46 276.8	0.44 281.1	0.44 289.0	0.56
	T 1g		330.0		1.25 307.2		1.12 315.6	1.36 306.7	1.28 310.0	1.25 316.6	1.15
	T 1g		380.0		0.07 360.4		0.05 369.1		0.04 362.3	0.04 368.0	0.05 367.3	0.06
	T 1g		408.0		0.32 383.1		0.21 394.5	0.22 384.7	0.18 390.8	0.18 392.4	0.26
	T 1g		460.5		0.20 442.4		0.14 454.8	0.15 445.7	0.13 452.8	0.13 448.5	0.17
	T 1g		488.3		0.04 470.8		0.01 482.7		0.02 473.4	0.01 480.8	0.01 476.3	0.02
	T 1g		560.9		0.40 535.5		0.30 545.8	0.30 536.8	0.28 546.4	0.29 544.9	0.35
	T 1g		810.0		0.00 765.7		0.00 748.4		0.00 750.9	0.00 768.7	0.00 789.1	0.00
	T 1g		833.9		0.11 799.6		0.10 795.2		0.11 792.6	0.10 808.1	0.10 815.8	0.10
	T 1g		921.2		0.10 879.9		0.10 870.4		0.10 870.8	0.10 888.0	0.10 900.0	0.09
	E u		234.9		0.07 220.3		0.07 222.2		0.13 217.9	0.08 223.3	0.09 226.8	0.07
	E u		300.0		0.82 276.5		0.70 283.9	0.70 275.8	0.67 281.4	0.73 286.6	0.76
	E u		330.0		0.03 310.2		0.03 312.9		0.06 307.3	0.04 319.3	0.01 318.5	0.02
	E u		390.1		0.31 366.7		0.22 374.0	0.22 364.1	0.22 371.4	0.21 376.5	0.26
	E u		419.8		0.03 395.7		0.06 399.4		0.14 395.7	0.10 407.4	0.08 405.4	0.04
	E u		454.1		0.01 428.6		0.02 423.9		0.01 420.4	0.02 434.6	0.02 441.3	0.01
	E u		550.5		0.26 525.2		0.21 536.4	0.18 526.8	0.18 538.4	0.19 534.5	0.24
	E u		630.1		0.14 605.2		0.11 617.2		0.12 608.1	0.11 617.5	0.11 614.2	0.12
	E u		851.9		0.10 816.7		0.09 811.0		0.10 809.1	0.09 824.9	0.09 833.5	0.09
	E u		884.1		0.00 843.0		0.00 826.8		0.00 830.4	0.00 845.5	0.00 864.4	0.00
	A 1u		300.6		0.00 285.0		0.00 288.1		0.00 283.0	0.00 295.1	0.00 291.8	0.00
	A 1u		378.6		0.00 361.0		0.00 370.1		0.00 362.9	0.00 371.9	0.00 367.4	0.00
	A 1u		460.8		0.00 435.9		0.00 439.7		0.00 432.0	0.00 444.1	0.00 446.6	0.00
	A 1u		652.7		0.00 624.6		0.00 627.4		0.00 621.6	0.00 633.5	0.00 637.1	0.00
	A 1u		976.7		0.00 933.1		0.00 926.3		0.00 924.9	0.00 940.4	0.00 953.9	0.00
	A 2u		189.9		0.04 179.2		0.04 187.3		0.04 179.8	0.04 180.0	0.04 182.8	0.03
	A 2u		294.3		0.39 274.7		0.36 282.2	0.37 274.3	0.35 276.4	0.36 282.8	0.37
	A 2u		440.1		0.31 413.2		0.23 420.1	0.26 412.1	0.23 421.6	0.22 423.9	0.26
	A 2u		482.4		0.00 460.5		0.01 460.6		0.00 457.4	0.00 470.8	0.01 470.5	0.01
	A 2u		845.3		0.00 808.9		0.00 804.6		0.00 801.8	0.00 819.8	0.00 826.2	0.00
	A 2g		201.0		0.14 190.9		0.13 196.1		0.14 191.5	0.14 195.2	0.14 194.7	0.13
	A 2g		376.8		0.38 346.3		0.24 357.5	0.36 344.3	0.32 348.8	0.32 359.3	0.35
	A 2g		440.0		0.20 417.2		0.14 428.9		0.14 418.3	0.14 424.5	0.14 425.9	0.16
	A 2g		562.3		0.00 544.5		0.00 555.2		0.00 548.3	0.00 556.0	0.00 550.2	0.00
	A 2g		968.7		0.03 927.7		0.03 916.0		0.03 918.0	0.03 934.2	0.03 947.8	0.03
	∑ h = ln β Γ		19.25			16.08		16.89		15.79		16.19		17.44
	∑ orphan h		11.56			9.58		10.30		9.60		9.84		10.54
	lnβ exact		19.36			16.28		17.06		15.93		16.36		17.61

  for Infrared. . Values of lnβ and h (in ) are computed at 300K for the44/40 Ca pair. χ is the deviation (in %) betwenn the theoretical frequency and its experimental counterpart. -1.0 0.01 121.9 -12.9 0.02 127.7 -8.8 0.00 115.5 -17.5 0.03 123.1 -12.1 0.00 129.6 -7.5 0.02 137 A g 182.5 0.3 0.85 170.7 -6.2 0.77 173.6 -4.6 0.85 168.0 -7.7 0.71 172.6 -5.2 0.84 175.7 -3.4 0.80 180 A g 234.0 0.0 0.05 221.2 -5.5 0.06 234.8 0.4 0.11 222.4 -4.9 0.09 234.1 0.0 0.00 224.8 -3.9 0.05 235 A g 250.6 -2.1 0.38 237.5 -7.2 0.44 242.6 -5.2 0.39 235.5 -8.0 0.47 245.7 -4.0 0.47 242.6 -5.2 0.42 253 A g 314.6 1.5 0.87 293.6 -5.3 0.41 299.6 -3.4 0.32 283.0 -8.7 0.08 302.5 -2.4 0.48 301.6 -2.7 0.66 307 A g 325.9 -0.0 0.06 304.0 -6.8 0.11 309.8 -5.0 0.19 299.8 -8.0 0.38 315.9 -3.1 0.10 314.9 -3.4 0.04 321 A g 369.3 2.6 0.01 335.3 -6.9 0.01 337.9 -6.1 0.02 328.5 -8.7 0.00 342.6 -4.8 0.04 351.4 -2.4 0.01 351 A g 396.9 1.8 0.04 368.6 -5.5 0.02 375.7 -3.7 0.01 365.1 -6.4 0.02 379.6 -2.7 0.03 379.9 -2.6 0.03 387 A g 498.0 -2.0 0.00 481.7 -5.2 0.01 497.5 -2.1 0.01 485.9 -4.4 0.00 494.6 -2.6 0.01 486.1 -4.3 0.01 518 A g 526.3 -0.7 0.04 507.0 -4.3 0.02 517.2 -2.4 0.01 510.3 -3.7 0.02 518.3 -2.2 0.03 513.0 -3.2 0.03 535 A g

											Raman
			PZ:			PBE		vdW-DF2		BLYP		BLYP+GD2		PBESOL		Exp B3LYP	(a)
	sym.	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	χ	h	ν	ν	χ
			(%) ( )		(%) ( )		(%) ( )	(%) ( )		(%) (	)		(%) ( )		(%)
	A g	138.6

Acknowledgments

This research was founded by l'Agence Nationale de la Recherche (ANR) through grant FILi-CaBeSo ANR-19-CE01-0016 (Fractionation of Isotopes in Liquids -Calcium Behaviour in Soils, PI Merlin Méheut). This work was performed using HPC resources from CALMIP (Calcul en Midi-Pyrénées; grant 2020-P1037).

698.9 -1.4 0.10 673.7 -5.0 0.09 676.9 -4.5 0.10 673.4 -5.0 0.09 685.2 -3.4 0.10 683.8 -3.5 0.10 715 B g 897.7 -2.1 0.00 856.4 -6.6 0.00 853.9 -6.9 0.00 845.8 -7.8 0.00 862.6 -5.9 0.00 876.4 -4.4 0.00 901 B g 958.3 -0.9 0.06 924.2 -4.4 0.07 913.0 -5.6 0.07 923. sym. 493.7 -2.2 0.00 472.6 -6.4 0.00 488.4 -3.3 0.01 473.0 -6.3 0.00 484.9 -4.0 0.02 479.6 -5.0 0.00 A u 507.2 0.00 487.8 0.00 502.0 0.00 490.1 0.00 503.9 0.00 493.7 0.00 A u 663.1 -1.0 0.02 644.4 -3.8 0.03 639.6 -4.5 0.02 646.4 -3.5 0.03 648.7 -3.2 0.02 651.7 -2.7 0.02 A u 899.0 -1.7 0.00 855.7 -6.5 0.00 852.4 -6.8 0.00 843.6 -7.8 0.00 862.5 -5.7 0.00 877.1 -4.1 0.00 A u 945.9 -1.5 0.00 912.4 -5.0 0.00 902.6 -6.0 0.00 911.9 -5.0 0.00 920.1 -4.2 0.00 927. 1 -3.4