J. Phys. Chem. C 2009, 113, 14507–14511
14507
Andradite-Uvarovite Solid Solutions. An ab Initio All-Electron Quantum Mechanical Simulation with the CRYSTAL06 Code A. Meyer,† Ph. D’Arco,‡ R. Orlando,§ and R. Dovesi*,† Dipartimento di Chimica IFM, UniVersita` di Torino and NIS -Nanostructured Interfaces and Surfaces - Centre of Excellence, Via P. Giuria 7, 10125 Torino, Italy, UniVersite´ Pierre et Marie Curie, Lab. PMMP, Tour 44-00, case 110, 4 Place Jussieu, 75005 Paris, France, and Dipartimento di Scienze e Tecnologie AVanzate, UniVersita` del Piemonte Orientale, Viale T. Michel 11, I-15121 Alessandria, Italy ReceiVed: April 21, 2009
Andradite-uvarovite (Ca3Fe2Si3O12-Ca3Cr2Si3O12) solid solutions have been investigated at an ab initio quantum-mechanical level by using an all-electron Gaussian-type basis set and the hybrid B3LYP functional in its unrestricted formulation. Only ferromagnetic phases have been considered. All possible nonequivalent geometrical configurations resulting from the substitution of Cr atoms for x ) 1-8 Fe atoms in the 16a site in the garnet primitive cell have been fully optimized (cell parameters and fractional coordinates of 80 atoms). As the lattice parameters of the two end-members are very similar (11.99 and 12.05 Å for uvarovite and andradite, respectively), geometry rearrangements at the various x are extremely small, the largest excess volume being 0.15 Å3 and the largest excess energy 3.68 kJ/mol. Thermodynamic functions are calculated from the configurational contribution to entropy and disregarding the vibrational contribution, which is expected to be negligible. Almost ideal miscibility is predicted. I. Introduction X32+Y23+Si3O12
Garnets are important rock-forming silicates, as major constituents of the Earth’s upper mantle and relevant phases of high-pressure metamorphic rocks in the Earth’s crust.1 Building blocks are SiO4 tetrahedra sharing corners with YO6 octahedra, whereas X2+ cations are in dodecahedral coordination. The structural,2 elastic,3,4 vibrational,5-8 thermodynamic,9,10 electronic,11 and magnetic12 properties of garnets have been extensively investigated with various experimental techniques. In regard to simulation, ab initio quantum mechanical methods have been used to investigate the vibrational properties and the dielectric and Born charge tensors of pyrope, grossular, andradite (Ca3Fe2Si3O12, indicated in the following as And), spessartine and uvarovite (Ca3Cr2Si3O12, Uva).13-16 Calculated and experimental IR and Raman5,17,18 spectra are found to be in excellent agreement. More recently the magnetic properties of Uva and And, two of the most abundant members of the family, have been investigated.19 Natural garnets never appear as pure end-members. Instead, they form solid solutions extending over a broad chemical range involving up to 12 end-members.20 The most common cases refer to substitutions of divalent cations in the X sites, and of trivalent cations in the Y sites, respectively. A miscibility gap is observed in mineral Ca-rich (pyralspite) and Ca-depleted (ugrandite) garnet solutions, which is not supported by any experimental evidence and its origin is not fully understood. Simultaneous substitutions at X and Y sites should then also be considered. Intermediate compositions, as well as the end-members, have been studied extensively. Several hundreds of structure refinements have been carried out,21 usually with no evidence of * To whom correspondence should be addressed. E-mail: roberto.dovesi@ unito.it. † Universita` di Torino and NIS. ‡ Universite´ Pierre et Marie Curie. § Universita` del Piemonte Orientale.
reduction in the cubic symmetry. Almandine-grossular22 and almandine-pyrope22,23 have been investigated with IR spectroscopy. Other properties, such as thermal expansion and elasticity,24-26 electron density,27 and optical and electronic structure28 have been explored, too. Substitutions at Y sites seem responsible for anomalous birefringence in some calcium rich garnets. Most data are available for the grossular-andradite series,29,30 but optically anisotropic grossular-uvarovite garnets have also been described.31,32 Optical anisotropy has been related to various features such as loss of the end-members cubic symmetry, secondary plastic deformation,33 early lattice mismatch at compositional twin boundaries,29,30 formation of defects such as the substitution of (OH)4 for SiO4 groups, which is accompanied by symmetry reduction,34 and others.33,35-37 In regards to simulation, semiclassical models38-41 have been for a long time the highest level of theory applicable to these complicated problems. More recently, two papers about pyropegrossular solid solutions have been published based on ab initio methods (LDA42 or GGA,43 plane-waves, pseudopotentials). Because of high computational cost, only a small subset of the many possible configurations has been explored in ref.42 by using the interaction parameters obtained within the cluster expansion scheme.44,45 To the authors’ knowledge, no ab initio study tackles solid solution involving transition metals (TM) in garnets in any of the two sites. Ab initio modeling of TM garnets, such as And and Uva, faces additional difficulties such as functionals of the LDA and GGA form fail in properly localizing d electrons of the TMs (d3 for Cr and d5 for Fe), as documented by the many studies on magnetism46-49 and spin polarized defects such as trapped holes.50 The open shell treatment increases the computational cost; the high atomic number of TM often requires the use of pseudopotentials, particularly when associated with plane-wave basis sets, that may be inappropriate for the description of the relatively deep d levels of TMs.
10.1021/jp903654w CCC: $40.75 2009 American Chemical Society Published on Web 07/17/2009
14508
J. Phys. Chem. C, Vol. 113, No. 32, 2009
Meyer et al.
TABLE 1: And-Uva Solid Solutionsa x
L
M
Ndf
D
a
0 1 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
1 8 4 12 12 8 24 24 2 6 6 8 12 12 24 8 24 24 4 12 12 8 1
4 38 19 56 57 38 114 114 9 28 27 38 57 56 114 38 114 114 19 56 57 38 4
1 1-5 3-4 1-6 5-6-7 1-2-5 1-2-7 5-6-7-8 1-2-5-6 1-2-7-8 1-6-7-8 1-2-5-8 1-2-5-7 1-2-3-5 5-6-7 1-2-5 1-2-7 1-5 3-4 1-6 1
12.193 12.185 12.177 12.176 12.177 12.168 12.170 12.169 12.160 12.159 12.160 12.160 12.161 12.159 12.160 12.152 12.153 12.150 12.144 12.144 12.143 12.135 12.127
3 4
5 6 7 8
b
12.176 12.168 12.168
c
12.178
12.169 12.162
12.159 12.161
12.161 12.161 12.161
12.151 12.151
12.152 12.154
12.143
12.144 12.144
R 90.000 89.970 89.942 90.000 90.000 90.033 90.098 90.028 90.000 90.000 90.000 90.000 89.998 90.001 89.964 89.975 89.899 89.949 90.049 90.000 90.001 90.024 90.000
β
γ
90.005 89.946 89.967 90.019 89.971
89.967 89.970 89.897
90.003 89.982 90.055 90.025 89.979 90.017
89.998 89.982 89.978 90.028 90.004
89.999
89.995 90.049
V
E
P
0.000 0.100 0.094 -0.049 0.109 0.019 0.084 0.059 -0.152 0.072 0.030 0.007 0.043 -0.066 -0.009 -0.022 0.129 -0.105 0.103 0.038 0.007 -0.012 0.000
0.000 1.045 2.176 2.557 1.835 3.647 2.274 2.035 3.248 3.274 -0.097 1.586 2.838 2.859 2.707 3.421 2.704 1.441 2.019 2.418 1.916 1.554 0.000
1.00 1.00 0.14 0.37 0.49 0.08 0.44 0.48 0.02 0.06 0.22 0.15 0.14 0.13 0.29 0.34 0.09 0.57 0.15 0.38 0.47 1.00 1.00
a For each configuration (labeled by L), its composition x (number of Cr atoms at the Y site), multiplicity M and the number of structural variables to be optimized Ndf are reported. The octahedral sites (see Figure 1) occupied by the minority atoms (Cr of Fe) are listed in column D. a, b, c (Å) and R, β, γ (in degrees) denote the pseudo-cubic cell parameters (if omitted, refer to the value reported on the left). V (Å3) and E (kJ/mol) are the excess volume and energy defined as the difference between the calculated values for the actual configuration and the weighted average of end-member values. P is the occupation probability of the Lth configuration among the different configurations of identical composition, estimated at 300 K.
In the present paper, we investigate the Uva-And series using an all electron Gaussian type basis set and the DFT-hybrid B3LYP functional. These ingredients, as implemented in the code CRYSTAL06,51 mostly overcome the difficulties mentioned above. The paper is organized as follows. The computational details for the ab initio calculations (basis set, density functional theory functional, optimization technique, computational parameters) are given in Section II A. The model for describing solid solutions is illustrated in Section II B. Sections III A and III B provide the results concerning equilibrium geometries, energies, and miscibility. II. Computational Details A. Method and Basis Set. The CRYSTAL06 program51 has been used for the present calculations. As in previous studies on end-member garnets,13-15 the hybrid B3LYP Hamiltonian52 has been employed, here in its unrestricted formulation. B3LYP is used widely and successfully in molecular quantum chemistry,53 as well as in solid-state calculations, where it has been shown to provide equilibrium geometries, vibrational frequencies,13-16,54-58 and magnetic properties such as the superexchange coupling constants47-49,59,60 in good agreement with experimental data. The same computational conditions (tolerances for the truncation of the infinite Coulomb and exchange sums, SCF convergence criteria, grid size for the integration of the DFT exchange and correlation contribution, number of points in the reciprocal space) as in our previous studies on garnets13-15,57 have been used. An all electron basis is used for all atoms. Oxygen, silicon, and calcium are described by (8s)(411sp)(1d), (8s)(6311sp)(1d), (8s)(6511sp)(21d), contractions; iron and chromium are described by (8s)(64111sp)(41d) contractions as used in previous papers.13-15 Full details of the basis sets and computational parameters can be found at the CRYSTAL Web site,61 along with the input and output files used for these calculations.
Structure optimizations were performed by use of analytical energy gradients with respect to atomic coordinates and unitcell parameters,62-64 and of the BFGS scheme for Hessian updating.65-68 The Hessian matrix of one end member, computed by numerical differentiation of the gradient as for the calculation of frequencies,13-15,57 is used as an initial guess for determining the step in the minimal energy search. Convergence was checked on both gradient components and extension of nuclear displacements (TOLDEG and TOLDEX51 were set to 3.0 × 10-5 Ha/ Bohr and 1.2 × 10-4 Bohr, respectively). B. The Model for the Solid Solutions. Solid solutions have been mimicked by using the 80 atoms primitive cell that contains eight equivalent trivalent cations in the And and Uva endmembers. Eight different compositions have been studied, characterized by x, the number of chromium atoms in the cell (correspondingly there are 8 - x Fe atoms). For every x, the 8!/x!(8 - x)! possible geometrical configurations can be grouped into symmetry equivalent classes, all members in a class being energetically degenerate, so that structure optimization is required only for one member per class. Such a classification is performed automatically by the CRYSTAL code. In the present case 23 classes exist, each labeled by L. As an example, consider the x ) 2 case. There are 8!/(2! 6!) ) 28 possible configurations grouped into 3 classes with multiplicity M equal to 4, 12, and 12 (12 + 12 + 4 ) 28). The unit cell representative of each class has been fully optimized, by relaxing the cubic symmetry constraint when necessary. The 48 point operators characterizing the original Oh symmetry reduce to 48/M, and the total number of atomic coordinates plus cell parameters to be optimized is indicated as Ndf in Table 1. Ndf increases from 4 for the end-members (one lattice parameter and the three fractional coordinates of oxygen) to 114 for the lowest symmetric cases (L ) 7, 8, 15, 17, 18). Column D in Table 1 describes the 23 nonequivalent And-Uva configurations; atomic labels refer to Figure 1. In the L ) 8 configuration, for example,
Andradite-Uvarovite Solid Solutions
Figure 1. Location of octahedrally coordinated atoms (Fe and Cr in this study) in the conventional cell of garnets. Atoms 1-8 are found within the primitive cell. Primed and unprimed atoms are related by primitive translations.
Cr atoms 1, 2, and 7 are substituted for Fe atoms. Note that 2 and 7 are nearest and next nearest neighbors of 1, respectively. Atoms 2′ and 7′, that are translationally equivalent to 2 and 7, are obviously substituted, too. In this preliminary work, only the highest spin ferromagnetic solutions have been considered. This choice is dictated by both practical reasons (the number of spin configurations for each substitutional configuration can be relatively high) and the fact that magnetic interactions entail energy differences of the order of two kJ/mol or smaller, that might affect our results only at very low temperature. III. Results A. Excess Energy and Volume. The calculated cubic cell parameters of the end-members are 12.197 (And) and 12.123 Å (Uva) to be compared with 12.0518-12.05969 Å (And), and 11.9886-11.99969 Å (Uva). Thus our calculations overestimate the experimental lattice parameters by about 1% in line with previous B3LYP results. The difference between end-member cells is, however, about the same for experiment (0.060 Å)69 and calculation (0.066 Å), indicating a sort of systematic shift. Upon substitution, symmetry reduces from cubic (L ) 1, 9, 23) to tetragonal (L ) 11), orthorhombic (L ) 10), monoclinic (L ) 4, 5, 13, 14), triclinic (L ) 7, 8, 15), and trigonal (L ) 2, 3, 6, 12). Cell parameters for each configuration are reported with reference to the (pseudo)-cubic lattice.70 Deviation from the cubic geometry is always very small. The largest angular deviation from 90° is 0.1° (L ) 7), and the largest difference among a, b, and c is 0.004 Å (L ) 18). In many cases (L ) 9, 11, 12), these differences are smaller than 10-4 Å and are not reported in Table 1. The deviation from cubic metric reported for uvarovite-grossular32 or grossular-andradite36 birefringent garnets are of the same order of magnitude as the ones found here. The excess volume (V ) and energy (E ), defined as the difference between the calculated values and the ones obtained by linear interpolation between end-members, are also reported in Table 1 for every configuration and represented as crosses in Figure 2. Molar volumes of the two end-members differ only by 29.28 Å3 (
