Electronic and elastic properties of brownmillerite

The elastic and electronic properties of brownmillerite (Ca2AlFeO5), the fourth component in mass of Portland cement, have been determined using Density Functional Theory. The Young’s modulus obtained in this work (E = 169 ± 2 GPa) is much closer to the experimental value than all the previously reported classical calculations. The electronic structure has been analyzed by calculating the band structure, the density of states and the crystal orbital overlap population curves. Results show that there are anti-bonding bonds between the iron and oxygen atoms which may create instability in the structure. The presence of those anti-bonding bonds might explain the limitation of force-field simulations when modeling the elastic properties of brownmillerite. Such a limitation may be extended to other structures derived from perovskites that present anti-bonding states.


Introduction
The mineral brownmillerite, Ca 2 AlFeO 5 , was found naturally for the first time in 1964 by Hentschel [1] in the Ettringer-Belleberg volcano in Germany. The structure of brownmillerite derives from the perovskite prototype structure with the maximum concentration of oxygen vacancies [2]. An ideal ABO 3 perovskite crystalline structure is usually formed by octahedrally oxygen-coordinated B 4+ cations and more voluminous A 2+ cations placed between them with a variable coordination [3]. In brownmillerite, the B sites are occupied by Fe 3+ and Al 3+ atoms, so for every two cations an oxygen must be removed to maintain the unit cell charge-balance [2]. In the Kronger-Vink notation [4] the substitutions are written as 2B 4+ =Fe 3+ +Al 3+ +v O where v O represents an oxygen vacancy. The brownmillerite crystal is the limiting case in which all the B 4+ atoms have been substituted by Fe or Al atoms, and the high concentration of vacancies leads to the preferential oxygen/vacancy ordering shown in figure 1. As a consequence of these modifications, the originally octahedral sites where the oxygen vacancy has been created become tetrahedrally coordinated. Brownmillerite crystallizes in the orthorhombic I2mb (No. 46) space group. Experimental observations by x-ray diffraction [5] and computational simulations [6] show that Al atoms are preferentially placed in the tetrahedral sites while the Fe atoms maintain the octahedral coordination. Overall, the Ca 4 FeAlO 5 structure is composed of alternating layers of FeO 6 octahedra and AlO 4 tetrahedra in the b direction as shown in figure 1. The AlO 4 tetrahedra form chains along the a direction, and the Ca atoms are placed between the FeO 6 and AlO 4 layers.
Before its discovery in nature, brownmillerite was already known as a major constituent of the ferrite phase in Portland cement with the name of C 4 AF. In cement it is formed after a fast cooling of raw material (clays and limestone) melted at 1600 K, and reaches the 5%-15% of commercial Portland Cement [7,8]. C 4 AF is the cement clinker phase with lowest reactivity with water, and its hydration products do not play a key role on the strength development of cement [8]. However, unreacted brownmillerite grains remain as aggregates in the cement paste and contribute to the elasticity of concrete. The elastic properties of brownmillerite were measured using nanoindentation by Velez et al [9], obtaining a Young's modulus of 125±25 GPa. However, computational studies using empirical force-fields [10,11] reported values between 60% and 130% higher. Such a disagreement between experiments and simulations is not common when using force-field methods [12]. The disagreement on the values reported in [11] could be attributed to the choice of unsuitable force-fields, as CLAYFF [13], DREIDING [14] and the universal Force Field (UFF) [15] interatomic potentials are not parameterized for iron in the Fe 3+ oxidation state. However, the core-shell force-field used in [10] was parameterized specifically for the study of brownmillerite [6] and reproduces accurately the structure and preferential arrangement of Fe and Al within the octahedral and tetrahedral sites. In addition, the core-shell family of polarizable potentials have proved to reproduce elastic properties of all the main cement phases accurately except those of brownmillerite [12,16,17].
Despite the massive worldwide production of cement, and consequently of brownmillerite, its electronic structure and elastic properties have not been previously investigated by Density Functional Theory (DFT) methods. In this paper we first compute the full elastic tensor of brownmillerite, derive the elastic properties and compare them with previous inaccurate predictions using empirical potentials. Then we investigate its electronic structure in detail, discussing its implications in the reactivity of brownmillerite. Furthermore, given the wide range of industrial applications of perovskites, such as capacitors in electronic circuits [18], microwave dielectric applications [19], solar cells [20] and materials to store high level radioactive waste [21], we compare the electronic structure of brownmillerite with that of the common perovskite CaTiO 3 .

Computational methods
All the simulations were done using the density functional theory (DFT) method as implemented in the SIESTA code [22]. For all the DFT calculations we used the general gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) parametrization [23] for the exchange-correlation functional. A double polarized basis set (DZP) was used to describe the valence electrons, and pseudopotentials in the Troullier-Martins method [24] flavour to mimic the effect of the core electrons. The Al, Ca, and O basis sets and pseudopotentials were previously validated for CaO and Al 2 O 3 [25], and the Fe for Fe 3 O 4 [26]. All the calculations were performed using a Monkhorst-Pack [27] k point spacing of 0.024Å −1 , which gives a mesh of 8×3×8 for the experimental structure. The energy minimization was done using the conjugate gradient method until the maximum difference in the density matrix was below 10 −4 or the maximum force was below 10 −4 eV/Å (for cell variable minimization a maximum value of 0.006 eV/Å 3 for the stress components was also considered).
Empirical simulations were done using GULP [28] and the core-shell potentials used by Zacate et al [6] to model the structure of brownmillerite. The energy minimizations were done using the conjugate gradient method until the gradient norm was smaller than 10 −3 . Then, the elastic constants were calculated using the second derivatives of the energy with respect to external stress.

Spin ordering
In order to analyze the spin ordering, the Fe atoms were spin-polarized. The ferromagnetic ordering and all the antiferromagnetic possible orderings compatible with the lattice periodicity (see figure 2) were relaxed and their energies were calculated. It is already known that Ca 2 Fe 2 O 5 brownmillerite is a G-type antiferromagnet [29], so these calculations serve as a test of our methodology. According to the results, the energies of the two structures with alternated intralayer spin polarization (G-type and C-type antiferromagnetics) are almost 1 eV lower than the other two magnetic orderings. Between the two low energy structures, the G-type ordering is slightly more stable due to the presence of both intralayer and interlayer alternated spin polarization. Therefore, having two atoms with the same spin polarization next to each other is penalized energetically, specially when they are located in the same Fe layer. The result matches perfectly with the experimental observation [9] and hereafter the G-type antiferromagnetic configuration was used for the rest of the calculations.

Crystalline structure
The crystal structure of brownmillerite is made of alternated layers of FeO 6 octahedra and AlO 4 tetrahedra. It's unit cell contains 4 formula units and it is orthorhombic with space grup I2mb (No. 46). The comparison between the experimental and the computed atomic positions is show in table 1. In table 2 the optimized orthorhombic lattice parameters of the G-type antiferromagnetic structure can be compared with the experimental values obtained by Readhammer et al [5] by x-ray powder diffraction for the Ca 2 (Fe 2−x Al x )O 5 structure with x=0.986. There is a good agreement between the experimental results and the computed values as the difference in every cell parameter is less than 3%, which is within the standard DFT-GGA error [30]. All the cell parameters are overestimated, a common trend of the DFT-GGA approximation [31,32]. The cell parameters calculated with the core-shell empirical potentials are also included in table 2. It can be seen that this force-field describes the lattice constants at the same level of accuracy as DFT calculations, which indicates its suitability for structural studies of brownmillerite despite its poor performance in elasticity. Proton transfer in perovskites [33,34] are highly influenced by the magnitude of the tilts of the oxygen octahedra at the B sites. In brownmillerite, the octahedra are tilted around the b direction in direct correlation with the Fe-O-Fe angle between two adjacent octahedra of the same layer. The lower the Fe-O-Fe angle, the bigger the tilting is. The computed values from both DFT and the core-shell force-field simulations present deviations within 2°from the experimental value [5]. Comparing this results with the common CaTiO 3 perovskite, which has an experimental Ti-O-Ti angle of 156°, we can observe that brownmillerite has a bigger tilt angle. Consequently, it is closer to the perfect cubic perovskite structure which presents no tilting.

Elastic properties
In order to calculate the elastic tensor coefficients, the strain-stress approach was used [12,25]: where ε and σ are the strain and stress tensors respectively and C ij are the elastic tensor coefficients in the Voigt notation. Equation (1) constitutes a set of six linear equations that can be simplified by deforming the unit cell in such a way that there is only one non-zero strain component. Then, we can relate each of the stress tensor components σ i with each of the strain tensor components ε j by just one elastic tensor component C ij : Strain was applied to each component of the strain tensor in steps of 1% from −3% (contraction) to +3% (expansion) by modifying the lattice vectors matrix following the relationship e = + L L 0 · ( ) where L 0 and L are the equilibrium and deformed lattice vectors matrices respectively, and e is the strain tensor.
Once all the elastic tensor components C ij are calculated from linear fit of the computationally obtained data, the bulk modulus (K ), shear modulus (G), Young's modulus (E) and the Poisson's ration (μ) can be calculated on the Voigt-Reuss-Hill (VRH) [35] approximation: where Voight's and Reuss's expressions for the bulk and shear modulus are given by: where C ij are the elastic tensor components calculated with the strain-stress method, and S ij are the components of the elastic compliance tensor C −1 . The calculated elastic constants are given in table 3. The obtained Young's modulus in the Voigt-Reuss-Hill approximation is E=169±2 GPa, 35.2% higher than the experimental result of 125±25 GPa obtained by Velez et al [9]. The disagreement might be due in part to the DFT methodology but also to the experimental conditions and approximations [9]. Nevertheless, the error is much smaller than that one obtained by the core-shell force-field method and our results almost overlap with the experimental ones. This discrepancy can be explained by the presence of antibonding bonds as discussed later in section 3.4. It can be noticed that most of the elastic tensor coefficients calculated with the core-shell force-field method are bigger than those calculated with DFT, resulting in large errors of the elastic coefficients. Moreover, those differences are notably more significant in some directions. With the aim of studying this anisotropy, the Young's modulus was represented as a three dimension parametric surface using the following definition [36,37]:   Table 3. Elastic constants of brownmillerite calculated with DFT and with the core-shell force-field method together with the previous experimental and computational results. All the values are given in GPa except ν (which is dimensionless). The percentage error of the core-shell values are given with respect to the DFT results.

DFT
Core-shell(error%) [6] CLAYFF [11] DREIDING [11] UFF [11] Exp [9] directions correspond respectively to the a, b and c crystallographic axes. The core-shell force-field method presents a more anisotropic result than DFT. In the XZ projection plane both approaches result in a similar shape. In the YZ plane DFT calculation has a rounded shape while the empirical potential displays a starlike shape which implies a more anisotropic behavior. Finally, the XY plane clearly shows that the X direction is overestimated in the core-shell force-field method. Overall, the Young's modulus calculated with the core-shell method is bigger than that calculated with DFT specially along the X and Z directions. The fact that the Y direction has the smallest Young's modulus in our DFT calculations matches with the structural discontinuity due to the oxygen vacancies. On the other hand, the X and Z directions have bigger values as there is always a linked octahedra and tetrahedra. The main drawback of the core-shell force-field method is that the elastic properties calculated along those directions that have continuous chains are greatly overestimated. Comparing the Young's modulus with the CaTiO 3 perovskite calculated with DFT [38] (E=263.84 GPa), brownmillerite has a lower Young's modulus. This is probably due to the oxygen vacancies which weaken the structure as well as the difference in chemical composition.

Electronic structure
The band structure shows that brownmillerite is an insulator with a band gap of 2.04 eV in the R point of the Brillouin zone [see figure 4(a)]. Comparing the band structure with the common CaTiO 3 perovskite the main difference is that in CaTiO 3 the gap is direct in the Γ symmetry point, yet the band gap is nearly the same (2.003-2.55 eV) [38][39][40]. The projected density of states of figure 4(b) shows that the valence bands correspond mostly to the 2p orbitals of the oxygen atoms, as usual in oxides and perovskites [39,41,42]. There is a significant overlap between these oxygen orbitals and the Fe 3d orbitals while it is negligible with other atom types (Ca and Al). From the analyzed valence and conduction bands and by analogy with the frontier orbital theory in molecules, we can draw conclusions on the reactivity of brownmillerite. Under attack by water molecules, as in cement hydration, the valence band regions would be the preferential points for the attack of water hydrogen atoms or hydronium cations. As expected for an oxide, those regions are located in the oxygen atoms, yet the reaction would be a priori more favourable with the oxygen atoms in the Fe octahedral plane. On the other hand, water oxygen atoms and hydroxyl anions would react preferentially with the conduction band regions, located on the Fe atoms.
In addition to the analysis of the spatial location of the valence and conduction bands, we have analyzed the bonding/antibonding character of the cation-oxygen bonds computing the crystal orbital overlap population (COOP) [see figure 4(b)]. The COOP can be understood as an extension of the Mulliken population analysis formalism to the periodic crystalline case [44,45], and are calculated multiplying the DOS by the overlap population of the desired orbitals. When the overall overlap is positive the bonds between those orbitals are bonding and when the overall overlap is negative the bonds between those orbitals are antibonding. Figure 4(b) shows that in the conduction region all the bonds between Al, Fe and Ca atoms with the O atoms are antibonding. However, below the Fermi level, the Ca-O and Al-O bonds are bonding, while the Fe-O are antibonding in the valence region between −3.6 eV and the Fermi energy. That indicates an instability of these particular bonds which would be prone to major reactivity. The spatial location of the antibonding Fe-O bonds can be seen looking at the valence local density of states [see figure 5(a)]. It must be noticed that the unstable Fe-O bonds are mostly located in the XZ plane. As discussed previously, these directions correspond to the maximum discrepancy between the core-shell force-field and DFT calculations. Therefore, the antibonding character of the Fe-O bonds might imply an instability that cannot be captured by empirical force-fields and this could be the reason of the huge disagreement in the Young's modulus.

Conclusion
The value for the Young's modulus of brownmillerite calculated with DFT is E DFT =169±2 GPa. The 3-dimensional representation of the Young's modulus shows that the b crystallographic direction presents a considerable smaller than a and c. This is attributed to the linked octahedra-tetrahedra. The obtained result is much closer to the experimental one (E exp =125±25 GPa) [9] than all the previous calculations [11]. The remaining disagreement may be due to the difference between the experimental sample (compressed powder) and the present calculation (single crystal).
The calculated band gap of 2.04 eV is very similar to that of the CaTiO 3 perovskite [40]. Analyzing the COOP curves, it is found that the Fe atoms have an antibonding bond with the oxygen atoms creating instability within the structure. Those anti-bonding bonds are mostly located in the ac plane, precisely the plane in which the Young's modulus is overestimated with the core-shell force-field calculation. Therefore, it could be stated that force-fields methods cannot reproduce correctly these antibonding bonds.
This paper shows that certain structural properties, the elasticity in the present study, might depend intrinsically on the electronic structure, and therefore cannot be evaluated with force-field methods. The present case might be a particular example of a more general trend regarding the weakness of force-field approaches when anti-bonding bonds are present, which should be considered when testing the accuracy of new forcefields.