There is growing interest to discover suitable calcium containing oxides that can be used as electrode materials in calcium ion batteries. A comprehensive computational investigation of ionic defects and Ca-ion diffusion in Ca-bearing oxide materials at the atomic level is important so as to predict their suitability for use in Ca-ion batteries. In this study, we apply atomistic simulation techniques to examine the energetics of defects, dopants, and Ca-ion diffusion in Ca3Fe2Si3O12. The calculations suggest that the Ca/Fe anti-site defect is the most favorable intrinsic defect causing such significant disorder, which would be sensitive to synthesis conditions. Diffusion of Ca2+ ions within Ca3Fe2Si3O12 is three-dimensional, with the activation energy of migration of 2.63 eV inferring slow ionic conductivity. The most favorable isovalent defects are Mn2+, Sc3+, and Ge4+ on Ca, Fe, and Si, respectively, for this process. The formation of extra calcium was considered to increase the capacity and diffusion of Ca in this material. It is found that Al3+ and Mn2+ are the candidate dopants on the Si and Fe sites, respectively, for this process and there is a reduction observed in the activation energies. The electronic structures of favorable dopant configurations are discussed using density functional theory simulations.
Although lithium-ion batteries have dominated for years,1–4 calcium-ion batteries (CIBs) are gaining attraction for large-scale storage applications due to cost, abundance, and enhanced safety.5,6 Compared to other multivalent battery systems such as Mg-ion batteries,7,8 the research activity of CIBs is in the early stage. CIBs are also expected to provide high energy and power density due to their two-electron reaction during intercalation. Furthermore, the electrode potential of Ca/Ca2+ is lower by 0.5 eV than that of Mg/Mg2+ and higher only by 0.1 eV than that of Li/Li+, implying higher cell voltage.9
Identification of suitable electrode materials is important in the development of CIBs. There are only a few published studies on electrode materials including layered compounds (e.g., V2O5),10 chevrel phases [e.g., CaMo6X8 (X = S, Se, and Te)],11 and Prussian blue analogs [e.g., NaMnFe(CN)6]9 for CIBs. Arroyo-de Dompablo et al.12 recently carried out a joint experimental and theoretical study on CaMn2O4 polymorphs and concluded that full Ca extraction from these polymorphs introduces a substantial volume change with an average voltage of 3.1 V and a high Ca-ion migration barrier. In a recent theoretical study, Torres et al.13 considered a few Ca based minerals including Ca3Mn2(SiO4)3, CaMn(SiO3)2, and CaMn(CO3)2 to calculate their diffusion barriers, average voltages, and volume changes over Ca de-insertion.
Materials consisting of iron with structures formed by linking through polyhedral units such as PO4, SiO4, and SO4, to provide structural stability, are of considerable interest in the development of electrode materials.14–16 Andradite (Ca3Fe2Si3O12)17 is a rock-forming silicate garnet and consists of three Ca2+ ions per formula unit, leading to high theoretical capacity, and SiO44− units offering structural stability via strong Si–O bonds. Furthermore, the low cost and availability of iron and silicon coupled with their non-toxicity also make this material ideal.
Computer modeling techniques based on pair-wise potentials are powerful tools to predict and understand the defects, ion transport, and dopant properties in solid-state oxide materials. A variety of oxide materials including battery materials have been examined using this technique.18–22 In the present study, using inter-atomic pair potentials, we characterize the defect, diffusion, and dopant properties of Ca3Fe2Si3O12 on an atomistic level. Furthermore, density functional theory (DFT) allowed us to study the electronic structures of doped configurations.
II. COMPUTATIONAL METHODS
Defect, diffusion, and dopant calculations were performed using the classical pair potential method as implemented in the GULP code,23 which describes ionic interactions in the form of long-range Coulomb attractions and short-range attractions (electron–electron repulsion and attractive dispersion). Buckingham potentials (refer to Table S1) were used to describe short-range interactions. Cell parameters and ionic positions were relaxed using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.24 Defects were modeled using the Mott–Littleton method.25 Ca-ion migration was calculated by considering two adjacent Ca vacancy sites as initial and final configurations. Seven interstitial Ca ions were considered in a direct linear route, and they were allowed to relax. The difference between the vacancy formation energy and the maximum energy along the diffusion path is defined as the activation energy of migration. Relative energies and trends will be consistent although the present model assumes that ions are fully charged and defect concentrations are at a dilute limit. Thermodynamically, the defect parameters (i.e., migration energies) can be defined through the comparison of the defective crystal to an isobaric (or isochoric) non-defective crystal. These defect formation parameters are interlinked by thermodynamic relations.26,27 The present calculations correspond to the isobaric parameters for the formation and the migration processes.28,29
Electronic structure calculations are based on DFT using plane wave basis sets and projected augmented wave (PAW) potentials. A plane wave DFT code VASP (Vienna Ab initio Simulation Package)30 which diagonalizes the Kohn–Sham equations was used. The calculations were performed using the generalized gradient approximation (GGA) as proposed by Perdew, Burke, and Ernzerhof (PBE).31 In all cases, a plane wave basis set with a cut-off of 500 eV and a 2 × 2 × 2 Monkhorst and Pack32 k-point mesh were used. Energy minimization calculations were performed using the conjugate gradient (CG) algorithm33 and a force tolerance value of 0.001 eV/Å. Short-range attractive interactions (dispersion) were modeled using a DFT + D3 method as implemented by Grimme et al.34
III. RESULTS AND DISCUSSION
A. Crystal structure of Ca3Fe2Si3O12
The crystal structure of Ca3Fe2Si3O12 at room temperature is cubic (space group , no. 230).17 Experimental lattice parameters reported by Novak and Gibbs.17 are a = b = c = 12.058 Å, α = β = γ = 90°, and V = 1753.18 Å3. An eight-fold coordination is observed for Ca2+ cations. The Fe2+ ions reside in the FeO6 octahedra. Tetrahedral (SiO4) units share their corners with FeO6 units in a three-dimensional network, as shown in Fig. 1. Energy minimization calculation was performed to relax both ionic positions and cell dimensions in order to obtain equilibrium lattice parameters. The results allowed us to check the quality of the Buckingham potentials used in the classical simulation and pseudopotentials together with basis sets used in the DFT simulation. There is good agreement between calculated and experimental values (refer to Table I).
|.||Calculated .||.|||∆| (%) .|
|Parameter .||Force field .||DFT .||Experiment17 .||Force field .||DFT .|
|.||Calculated .||.|||∆| (%) .|
|Parameter .||Force field .||DFT .||Experiment17 .||Force field .||DFT .|
B. Intrinsic atomic defects
Diffusion of ions in an ionic material is typically influenced by Schottky and Frenkel disorders, which can be calculated by combining point defects (vacancies and interstitials). These point defect energies can be difficult to determine via experimentation. The current classical simulation technique allowed us to calculate a series of point defects. The lowest energy point defects were then combined to calculate Frenkel and Schottky defect processes. Anti-site defects in which cations exchange their positions were also calculated as these defects have been observed in many oxide materials both experimentally and theoretically.35–40 Here, we write reaction equations for all possible Schottky, Frenkel, and anti-site defects in Ca3Fe2Si3O12 using the Kröger–Vink notation,41
Figure 2 shows the resulting defect energies. The most favorable intrinsic disorder type is the Ca–Fe anti-site defect cluster (0.83 eV/defect). In this defect process, it is anticipated that there would be some Ca2+ ions on the Fe site and Fe3+ ions on the Ca site. In the isolated form of the anti-site defect, defects were modeled individually, and their energies were combined. Both defects were considered simultaneously in the cluster form. The energy difference between the isolated and cluster forms is considered as the binding energy (−0.19 eV). Notably, other anti-site defect cluster energies (Ca–Si and Fe–Si) are lower than the Schottky and Frenkel defect energies. Among other defect processes, the CaO Schottky is the lowest energy process, with a defect energy of 3.73 eV/defect. This process will facilitate the formation of and at high temperatures. The oxygen Frenkel energy is calculated to be 3.98 eV/defect, with a deviation of only 0.25 eV from the CaO Schottky. However, the Ca Frenkel energy is calculated to be ∼5.00 eV, inferring that the formation of Ca vacancies is not significant. Other Frenkel and Schottky energies are highly endoergic, suggesting that they will not be present at significant concentrations.
C. Self-diffusion of Ca2+ ions
The diffusion of Ca2+ ions is discussed in this section. A material with high ionic conductivity generally exhibits a high rate performance and is suitable as an electrode for use in devices such as batteries. Classical simulation allowed us to calculate the Ca-ion diffusion pathways and the corresponding activation energies in Ca3Fe2Si3O12. In general, diffusion pathways are difficult to examine in the experiment, particularly for complex oxides such as Ca3Fe2Si3O12. The current classical simulation technique has been successfully applied in previous studies to calculate ion migration pathways together with activation energies.42,43 For example, the Li-ion migration pathway calculated in LiFePO4 by Fisher et al.44 was in excellent agreement with the neutron diffraction experiment.45
A potential Ca–Ca hop with a jump distance of 3.70 Å was identified. Many of these identical hops were then connected to form a three-dimensional long-range migration path, as shown in Fig. 3. The activation energy for this hop is calculated to be 2.63 eV, and Ca-ions migrate in a curved pathway. High activation energy of Ca2+ ion migration in Ca3Fe2Si3O12 shows that the Ca-ion diffusion is slow in this material. A possible solution to increase the ion diffusion is by reducing the Ca–Ca distance in this material via synthesizing nanoparticles.46 Table II lists the activation energies calculated for Ca-ion migration in different Ca-based oxide materials.12,13,47 Torres et al.13 recently performed DFT calculations on a variety of Ca based minerals and reported the activation energies of Ca-ion migration. Interestingly, Ca3Cr2Si3O12 and Ca3Mn2Si3O12 materials, isostructural with Ca3Fe2Si3O12, were also considered, and their activation energies of migration were reported to be 2.07 eV and 2.09 eV, respectively.13 Activation energies calculated for Ca-ion migration are lower by ∼0.50 eV than those calculated in a similar Ca3Fe2Si3O12 structure. This is due to different methodologies considered and different transition elements present in the crystal structures. Activation energies are greater than 2 eV in all cases, inferring poor ionic conductivity. A possible reason can be due to the strong electrostatic interaction between the migrating Ca2+ cations and the other ions in the lattice.
|Ca-based oxide material .||Activation energy (eV) .|
|Ca-based oxide material .||Activation energy (eV) .|
D. Solution of dopants
This section examines the introduction of cations via divalent, trivalent, and tetravalent doping into Ca3Fe2Si3O12. Doping can result in beneficial or detrimental effects on the material. In particular, suitable dopants can stabilize a material in order to prevent unfavorable phase transformations. The aim here is to provide guidance to the dopants that require the least energy to integrate into the lattice, the effects of which can be studied in subsequent experimental testing.
1. Divalent dopants
Divalent cations (M = Ni, Mg, Co, Mn, Sr, and Ba) were first considered for doping on the Ca site. The doping mechanism can be seen in the following equation:
The calculated solution enthalpies are reported in Fig. 4. The lowest solution enthalpy (0.21 eV) is predicted for Mn2+, inferring the highest solubility. The lowest solution enthalpy suggests that Mn2+ would be an ideal candidate to stabilize Ca3Fe2Si3O12. The ionic radius of Ca2+ in an eight coordination environment is 1.00 Å. Favorability of the Mn2+ dopant can be due to its ionic radius (0.96 Å), which is close to the ionic radius of the Ca2+ ion. Both Mg2+ and Co2+ exhibit an identical solution enthalpy of 0.59 eV. This is due to a very small difference in their ionic radii. Solution enthalpy increases with the difference in the ionic radius between dopants and the Ca2+ ion. The solution enthalpies of Sr2+ and Ni2+ were calculated to be 0.64 eV and 0.77 eV, respectively. The highest solution enthalpy of 2.20 eV is calculated for Ba2+, suggesting that this dopant is very unfavorable for doping on the Ca site.
To increase the capacity of Ca3Fe2Si3O12, divalent dopants were considered on the Fe site. The following reaction equation [Eq. (16)] shows that this doping strategy can introduce additional Ca2+ ions in the form of Ca interstitials in the lattice. A previous experimental study considered Co3+ doping on the Ru site in Li2RuO3 to increase the Li+ concentration, and doping resulted in an enhancement in the electrochemical reversibility of Li+ ions and the concentration of Li+ ions in the lattice,48
The favorable dopant for this process is Mn2+, and its solution energy is calculated to be 3.36 eV. Solution energies calculated for Co2+, Ni2+, Mg2+, and Zn2+ are very close to the value calculated for Mn2+; therefore, they are also candidate dopants for experimental testing (refer to Fig. 5). The ionic radius of Fe3+ is 0.65 Å. Favorability of the five dopants (Mn, Co, Ni, Mg, and Zn) is due to their ionic radii (0.65 Å–0.74 Å) matching closely with the ionic radius of Fe3+. The ionic radius of Sr2+ (1.18 Å) deviates significantly from that of Fe3+, reflecting in solution enthalpy. The highest positive solution enthalpy of 6.17 eV is calculated for Ba2+, suggesting that this dopant is highly unfavorable on the Fe site.
2. Trivalent dopants
Here, we consider a range of trivalent cations (M = Al, Ga, In, Sc, Y, Gd, and La) on the Fe and Si sites. Doping on the Fe site produced no charge compensating defects as the charge on Fe is +3 [refer to Eq. (17)]. Solution enthalpies are plotted as a function of ionic radii of dopants in Fig. 4(b),
The most favorable dopant is Sc3+, with a solution enthalpy of 0.09 eV. The solution enthalpy of In3+ is higher only by 0.02 eV than that calculated for Sc3+, meaning that In3+ is also a promising dopant. Favorability of these two dopants can be partly due to their ionic radii being closer to the ionic radius of Fe3+ (0.65 Å). The solution enthalpy of Ga3+ is 0.29 eV although its ionic radius (0.62 Å) is very close to the ionic radius of Fe3+. Other dopants exhibit high solution enthalpies due to their ionic radii deviating from the ionic radius of Fe3+ [refer to Fig. 4(b)].
The formation of additional Ca2+ ions in the lattice was considered by doping trivalent cations on the Si site. As we discussed earlier, this process can also introduce Ca interstitials, as shown in the following reaction,
The solution enthalpies are reported in Fig. 5(b). In all cases, solution enthalpies are highly positive (>5 eV), suggesting that this process required high energy possibly in the form of heat. The most favorable dopant is Al3+, with a solution enthalpy of 5.49 eV. The lowest solution enthalpy for Al3+ can be due to its ionic radius (0.39 Å) being closer to the ionic radius of Si4+ (0.26 Å). High endoergic solution enthalpy for this mechanism can be due to the high defect formation energy of quadruply charged Si.
3. Tetravalent dopants
Finally, tetravalent dopants (M = Ge, Ti, Sn, Zr, and Ce) were considered on the Si site. Solution enthalpies are reported in Fig. 4(c). The following defect equation explains the mechanism involved in this process,
The most favorable solution energy (0.43 eV) is calculated for Ge4+. This is due to the ionic radius of Ge4+ (0.39 Å) being closer to that of Si4+ (0.26 Å). The solution enthalpy for Ti4+ is 3.86 eV, meaning that it is an unfavorable dopant. The second favorable dopant is Sn4+, with a solution enthalpy of 2.10 eV. Solution enthalpy gradually increases with an increase in the ionic radius. The highest enthalpy of solution for Ce4+ (4.12 eV) suggests that doping would require external heat energy.
E. Ca-diffusion in the presence of dopants
Here, we calculate the activation energies for Ca-ion diffusion in the presence of Mn2+ on the Fe site and Al3+ on the Si site. Figure 6 shows the energy profile diagrams for the local Ca hops. There is a reduction (by 0.45 eV) in the activation energy for the doping of Mn2+ on the Fe site [refer to Fig. 6(a)]. In the case of doping of Al3+ on the Si site, activation energy is reduced only by 0.09 eV [refer to Fig. 6(b)]. In both cases, Ca–Ca distances have been slightly elongated compared to the distances present in the un-doped crystal structure. This perturbation in distances is due to the charge and ionic radius mismatch between Mn2+ and Fe3+ and Al3+ and Si4+. The current simulation shows that doping of Mn2+ and Al3+ would simultaneously increase the concentration of Ca2+ ions in the lattice and reduce the activation energy of the Ca-ion migration.
F. Electronic structure of doped Ca3Fe2Si3O12
DFT simulations were used to examine the electronic structures of doped configurations. Here, we discuss the results of the most favorable dopants (Mn2+, Ga3+, Al3+, and Ge4+), as discussed in the previous sections.
The DOS plot calculated for bulk Ca3Fe2Si3O12 is shown in Fig. 6(a). Metallic behavior is observed for the ferromagnetic configuration of bulk Ca3Fe2Si3O12. Doping of Mn alters the electronic structure slightly, as shown in the DOS plot [refer to Fig. 7(b)]. Peaks arising from the 3d states of Mn appear near the Fermi level [refer to Fig. 7(c)]. Calculated Mn–O bond distances (2.20 Å–2.50 Å) [refer to Fig. 7(d)] are slightly shorter than the Ca–O (2.37 Å–2.54 Å) bond distances due to the ionic radius of Mn2+ being smaller than that of Ca2+. The constant charge density plot of the Mn-doped configuration is shown in Fig. 7(e).
In the case of Sc, the Fermi level does not change significantly, and the doped configuration maintains its metallic behavior [refer to Fig. 7(f)]. The Fermi level does not have any contribution from Sc. States associated with Sc appear at ∼8 eV in the conduction band [refer to Fig. 7(g)]. Sc–O bond distances are slightly longer (by ∼0.04 Å) than the Fe–O bond distances. This is due to the ionic radius of Sc3+ being larger than that of Fe3+, as discussed in Sec. III D 2.
Doping of Ge on the Si site has a little effect on the DOS structure [refer to Fig. 7(j)]. The doped configuration exhibits metallic behavior. Small peaks that belong to the 3s and 3p states of Ge appear just below the Fermi level [refer to Fig. 7(k)]. The Si–O bond distance in pristine Ca3Fe2Si3O12 is 1.66 Å. The Ge–O bond distances are longer than (by 0.15 Å) the Si–O bond distances. This is because of the ionic radius of Ge4+ (0.39 Å) being larger than that of Si4+ (0.26 Å).
Doping of Al on the Si site shifts the Fermi energy level only by 0.02 eV [refer to Fig. 8(b)]. The doped configuration is metallic. Additional peaks arising from 3s and 3p states of Al appear near the Fermi level [refer to Fig. 8(c)]. The calculated Al–O bond distances (1.77 Å) are longer than the Si–O bond distances (1.66 Å). This is because of the ionic radius and charge density of Al3+ being larger and smaller, respectively, than that of Si4+.
There is a slight shift in the Fermi energy level by 0.14 eV upon doping of Mn on the Fe site. Peaks associated with Mn appear near the Fermi level and in the valence band. Metallic behavior is still maintained. A very small perturbation is observed in the Mn–O bond distances compared to Fe–O bond distances. This is due to the small difference in the ionic radius and charge mismatch. The constant charge density plot of Mn doped on the Fe site in Ca3Fe2Si3O12 is shown in Fig. 8(h).
Computational simulation was applied to analyze the behavior of the defects, diffusion pathways, and dopant properties in Ca3Fe2Si3O12. The Ca/Fe anti-site defect was predicted to be the most energetically favorable intrinsic defect process. The activation energy for three dimensional Ca vacancy migration in Ca3Fe2Si3O12 was calculated to be 2.63 eV, inferring low calcium mobility. Low solution energies were calculated for Mn2+, Sc3+, and Ge4+ on the Ca, Fe, and Si sites, respectively. Here, we show that Al3+ and Mn2+ are promising dopants on the Si and Fe sites, respectively, to introduce additional Ca in the form of Ca interstitials in Ca3Fe2Si3O12 for improving its capacity and Ca-ion diffusion. The difference in the electronic properties between un-doped and doped configurations was explained using electronic structures calculated using DFT.
See the supplementary material for the Buckingham potentials used in the classical simulation of the current study.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Computational facilities and support were provided by the High Performance Computing Centre at Imperial College London.
This research was financially supported by the European Union’s H2020 Programme under Grant Agreement No. 824072–HARVESTORE.
The authors declare that there is no competing financial interest.