Simulation-Based Defect Engineering in “α-Spodumene”

Naturally occurring lithium-rich α-spodumene (α-LiAlSi2O6) is a technologically important mineral that has attracted considerable attention in ceramics, polymer industries, and rechargeable lithium ion batteries (LIBs). The defect chemistry and dopant properties of this material are studied using a well-established atomistic simulation technique based on classical pair-potentials. The most favorable intrinsic defect process is the Al-Si anti-site defect cluster (1.08 eV/defect). The second most favorable defect process is the Li-Al anti-site defect cluster (1.17 eV/defect). The Li-Frenkel is higher in energy by 0.33 eV than the Al-Si anti-site defect cluster. This process would ensure the formation of Li vacancies required for the Li diffusion via the vacancy-assisted mechanism. The Li-ion diffusion in this material is slow, with an activation energy of 2.62 eV. The most promising isovalent dopants on the Li, Al, and Si sites are found to be Na, Ga, and Ge, respectively. The formation of both Li interstitials and oxygen vacancies can be facilitated by doping of Ga on the Si site. The incorporation of lithium is studied using density functional theory simulations and the electronic structures of resultant complexes are discussed.


Introduction
Escalating demand for mineral candidates, especially in industries of lithium ion batteries (LIBs), has fueled the transformational shift in mobile-electronics. In particular, lithium-rich minerals have drawn the attention of investigators to analyze their potential candidacy [1]. Lithium-rich minerals are significant in modern industries of LIBs, alloys, and pharmacy [2][3][4]. Spodumene is a lithium aluminosilicate system, which is categorized as a granitic-pegmatite in geo-mineralogy [5][6][7]. The theoretical Li content in the spodumene is reported to be 3.73% [8], where it offers a theoretical Li 2 O content of 8 wt % [9]. Mineralogists describe that the spodumene has some associations with quartz and albite [10] and as a three-phase system [11,12]. α, β, and γ are the three major phases of spodumene, where α-spodumene (LiAlSi 2 O 6 ) is the mined natural material [1]. In addition, α-spodumene has great potential to serve as a main lithium resource for lithium industries in comparison with the other two phases [13].
Previous computational investigations have shown the significance of atomistic scale computation to reveal the defect energetic features of the minerals for energy applications [14][15][16][17]. Spodumene has been experimentally studied in previous cases for its phase transformation mechanisms [18] and to clarify the density of defects in natural α-spodumene [19]. Systematic computational studies have not been carried out previously and this investigation presents the first atomistic modelling defect energetic study of αspodumene. Understanding the fundamental defect processes in α-LiAlSi 2 O 6 can be useful to optimize its performance in a wide industrial portrait. We have attempted to elucidate the energetics of the optimized theoretical model of α-spodumene and its intrinsic defect processes employing atomistic simulations based on classical pair-potentials. Density functional theory (DFT) simulations were performed to examine the stability of Li-incorporated α-spodumene and the electronic structures of the resultant complexes.

Computational Methods
Defect calculations were performed using a classical pair-wise potential simulation code GULP (general utility lattice program) in which lattice energy optimization procedures are concretized [20]. A summary of the aspects of these techniques will be presented as comprehensive reviews are given elsewhere [21]. Perfect and defective lattices are simulated on the basis of potential system attributed to the energy as attested by the atomic position coordinates within the crystal lattice. Inter-ionic interactions are described in terms of long-range Coulombic interactions and short-range Buckingham potential parameters (refer to Table 1) [22][23][24][25]. Ionic polarization is treated using the core-shell model where the instances with large spring values are approached with core alone treatment as the shell has no charge. The Mott-Littleton approach [26] was utilized to calculate defect formation energies in the lattice environment, in which the lattice is portioned as region I and region II. It is foreseen that the calculated defect process energies will be overvalued owing to the spherical treatment of ions at dilution limits. Although such an overestimation prevails in the present methodology, trends in relative energies will be congruous [27]. The contemporary atomistic scale simulations correspond to isobaric (i.e., constant pressure condition) parameters for the processes [28,29]. DFT simulation as implemented in the Vienna ab initio simulation program (VASP) code [30] was used to model the incorporation of Li. In this code, projected augmented wave (PAW) potentials [31] and plane wave basis sets are used. A plane wave basis set with a cut-off of 500 eV was used in all calculations. For bulk α-LiAlSi 2 O 6 , we used 4 × 4 × 8 Monkhorst-Pack [32] k-points. A 2 × 2 × 2 super cell containing 320 atoms was used to model Li incorporation. For this supercell, 2 × 2 × 4 Monkhorst-Pack k-points were used. The generalized gradient approximation (GGA) as parameterized by Perdew, Burke, and Ernzerhof (PBE) [33] was used to model exchange correlation. Full geometry optimizations were performed to relax both atomic positions and cell parameters using the conjugate gradient algorithm [34]. We used zero damping DFT+D3 as parameterized by Grimme et al. [35] to model short-range dispersive interactions.

Computational Modelling of α-LiAlSi 2 O 6
The aluminosilicate α-spodumene crystallizes in the space group C 2/c with a monoclinic structure, where Al 3+ and Si 4+ exhibit octahedral and tetrahedral coordination, respectively [36]. The experimental crystal structure of α-LiAlSi 2 O 6 (lattice parameters: a = 9.461 Å, b = 8.395 Å, c = 5.217 Å, α = γ = 90.00 • , and β = 110.09 • ) [36] was reproduced by both classical and DFT simulations, where the calculated structural parameters are in good agreement with the experimental values (refer to Table 2). Figure 1 illustrates the structure and chemical environment of α-spodumene. This parameterized model of the α-LiAlSi 2 O 6 thus provides a valid fundamental base for subsequent defect modelling and calculation.

Intrinsic Defects
In this section, we report the results of intrinsic defect energies calculated using classical simulation. First, point defect (vacancy and interstitial) energies were calculated (refer to Table 3). They were then combined together with appropriate lattice energies to calculate Frenkel and Schottky defect energies. Anti-site defects in which cations exchange their atomic positions were also calculated in the forms of isolated and cluster. In the isolated form, the formation energies of impurities were calculated separately. In contrast, impurities were allowed to form next to each other in the cluster form and formation energy for the defect cluster was calculated simultaneously. Anti-site defect is an important defect that dominate the diffusion property of a material. The following reaction equations show the intrinsic defect processes (Schottky, Frenkel, and anti-site), as written using the Kröger-Vink notation [37].

Intrinsic Defects
In this section, we report the results of intrinsic defect energies calculated using classical simulation. First, point defect (vacancy and interstitial) energies were calculated (refer to Table 3). They were then combined together with appropriate lattice energies to calculate Frenkel and Schottky defect energies. Anti-site defects in which cations exchange their atomic positions were also calculated in the forms of isolated and cluster. In the isolated form, the formation energies of impurities were calculated separately. In contrast, impurities were allowed to form next to each other in the cluster form and formation energy for the defect cluster was calculated simultaneously. Anti-site defect is an important defect that dominate the diffusion property of a material. The following reaction equations show the intrinsic defect processes (Schottky, Frenkel, and anti-site), as written using the Kröger-Vink notation [37].
ChemEngineering 2021, 5, 57 4 of 12 Formation energies of all vacancies are endothermic. Formation energy increases with the increase in the charge of ion formed (Li > O > Al > Si). Conversely, all interstitial formation energies are exoergic and formation energy increases with the decrease in the charge of intertitial ion.
In Table 4, we report the calculated defect energies. The most favourable defect energy process is the Al-Si anti-site defect cluster Equation (14). This indicates that a small percentage of Al-Si cation mixing would be present in this material. The cluster form is energetically more favourable (by 0.31 eV) than its isolated form. This is owing to the isolated impurity defects preferring to form a cluster with exothermic binding energy of −0.31 eV. The second most favourable defect process is the Li-Al anti-site defect cluster. The defect energy is higher only by 0.09 eV than the Al-Si anti-site cluster. Anti-site defects have been found in many oxide materials experimentally and theoretically [38][39][40][41]. The Li-Frenkel is higher in energy by 0.33 eV than the Al-Si anti-site defect cluster. This Frenkel would ensure the formation of Li vacancies needed for the vacancy mediated Li-ion diffusion. Other Frenkel defect processes are highly endoergic, meaning that the Li 2 O partial Schottky is the lowest energy process compared with the other Schottky defect processes. This process will lead to the loss of Li 2 O in this material. However, this process would require energy in the form of heat. Table 4. Formation energies of intrinsic defect process in monoclinic α-LiAlSi 2 O 6 .

Li-Ion Diffusion
Ionic conductivity is an important parameter determining the property of a material. Here, we use classical simulation to examine the Li-ion diffusion pathways and their activation energies. In previous simulation studies [42][43][44][45], many oxide materials have been analysed to construct Li-ion migration pathways together with activation energies. For a promising Li-ion battery, electrode material with low activation energy of Li-ion migration is preferred. Computational simulations have been supportive to experimental studies in determining the diffusion pathways, as experimental investigation of diffusion is often challenging.
Four different Li local hops (A-D) were identified (see Figure 2). The Li hop distances and activation energies are reported in Table 5. The energy profile diagrams (see Figure 3) show the activation energies calculated for local Li hops. The local hop B exhibits the lowest activation energy of 2.61 eV. We constructed a long-range diffusion pathway (B→B→B→B) consisting of local hops B. Li-ion migrates in a zig-zag pattern in the ac plane (see Figure 2) in this diffusion pathway. Although this long-range pathway exhibits the lowest activation energy of 2.61 eV, diffusion of Li-ion is expected to be sluggish. For a promising electrode material, activation energy of Li-ion migration should be lower than 0.50 eV. Other local hops (A, C and D) have higher activation energies than that calculated for hop B. Although other long-range diffusion pathways [A→A→A→A, C→C→C→C and D→D→D→D] can be constructed, their overall activation energies are quite high. Migration calculations show that Li-ion conduction in this material is very slow. This is partly owing to the fact that Li-Li separation in this material is large. There is no experimental information available on the Li-ion migration pathways or activation energies in this material. Current simulation study can be useful to the future experimental verification. An efficient strategy to enhance the Li-ion diffusion would be introducing additional Li by appropriate doping. We discuss this strategy in the next section.
and D→D→D→D] can be constructed, their overall activation energies are quite high. Migration calculations show that Li-ion conduction in this material is very slow. This is partly owing to the fact that Li-Li separation in this material is large. There is no experimental information available on the Li-ion migration pathways or activation energies in this material. Current simulation study can be useful to the future experimental verification. An efficient strategy to enhance the Li-ion diffusion would be introducing additional Li by appropriate doping. We discuss this strategy in the next section.

Solution of Dopants
Cation doping is an important strategy to modify the properties of materials. Here, we consider isovalent and alovalent dopants at different cation sites. Charge compensating defects are not necessary in the case of isovalent doping. However, point defects (vacancies and interstitials) are used to compensate charges produced upon aliovalent doping. A range of cations were doped and their solution energies were calculated using classical simulation. We provide Buckingham potentials used for dopants in the Supplementary Information (see Table S1).

Solution of Dopants
Cation doping is an important strategy to modify the properties of materials. Here, we consider isovalent and alovalent dopants at different cation sites. Charge compensating defects are not necessary in the case of isovalent doping. However, point defects (vacancies and interstitials) are used to compensate charges produced upon aliovalent doping. A range of cations were doped and their solution energies were calculated using classical ChemEngineering 2021, 5, 57 7 of 12 simulation. We provide Buckingham potentials used for dopants in the Supplementary Information (see Table S1).

Isovalent Doping
First, monovalent dopants (M = Na, K, and Rb) were considered on the Li site and the following reaction equation was used to calculate the solution energy.
Solution energies are reported in Figure 4a. An exothermic solution energy of −0.24 eV is calculated for Na + , meaning that this dopant can be worth examining experimentally. The ionic radius of Li + is 0.76 Å [46]. Favourability of Na + on the Li site is partly owing to the fact that the ionic radius (1.02 Å) [46] of Na + is closer to that of Li + . Solution energy increases with increasing ionic radius. Doping of K and Rb is not favourable as solution energies are positive for both dopants. The largest solution of energy (5.23 eV) is calculated for Rb and this dopant is highly unlikely to be considered for doping at room temperature. The lowest solution energy is calculated for Ga with an exothermic solution energy of -0.14 eV (refer to Figure 4b). The preference of Ga is partly due to the fact that the ionic radius of Al 3+ (0.54 Å) matches reasonably with the ionic radius of Ga 3+ (0.62 Å) [46]. Other dopants have endoergic solution energies. Solution energy increases gradually with ionic radius. The La is the most unlikely dopant as it has the highest solution energy of 1.91 eV. Finally, the solution of tetravalent dopants (M = Ge, Ti, Sn, Zr, and Ce) was examined. Here, tetravalent dopants were substituted on the Si site. The following reaction equation was used to describe the doping process: All dopants exhibit endothermic solution energies (refer to Figure 4c). A promising dopant for this process is Ge, with a solution energy of 0.39 eV owing to the ionic radius of Si 4+ (0.26 Å) being closer to that of Ge 4+ (0.39 Å) [46]. There is an increase in solution energy for Ti (2.53 eV). Then, the solution energy decreases to 1.63 eV. Both Sn and Zr have lower solution energies than that of Ti. The highest positive solution energy is calculated for Ce.

Aliovalent Doping
Trivalent dopants (M = Ga, Sc, In, Y, Gd, and La) were substitutionally doped on the site. This process required charge compensating defects. Two possible such defects (Li interstitial or oxygen vacancy) were considered (refer to Equations (18) and (19)). Additional Li + ions in this materials would increase its capacity. The formation of oxygen vacancies would enhance the vacancy mediated oxygen ion diffusion. Furthermore, this doping process can lead to the loss of Li2O as Li interstitials and oxygen interstitials as facilitated by the formation oxygen vacancies can aggregate. Next, trivalent dopants (M = Ga, Sc, In, Y, Gd, and La) were doped on the Al site. We describe the reaction process using the following equation: The lowest solution energy is calculated for Ga with an exothermic solution energy of −0.14 eV (refer to Figure 4b). The preference of Ga is partly due to the fact that the ionic radius of Al 3+ (0.54 Å) matches reasonably with the ionic radius of Ga 3+ (0.62 Å) [46]. Other dopants have endoergic solution energies. Solution energy increases gradually with ionic radius. The La is the most unlikely dopant as it has the highest solution energy of 1.91 eV.
Finally, the solution of tetravalent dopants (M = Ge, Ti, Sn, Zr, and Ce) was examined. Here, tetravalent dopants were substituted on the Si site. The following reaction equation was used to describe the doping process: All dopants exhibit endothermic solution energies (refer to Figure 4c). A promising dopant for this process is Ge, with a solution energy of 0.39 eV owing to the ionic radius of Si 4+ (0.26 Å) being closer to that of Ge 4+ (0.39 Å) [46]. There is an increase in solution energy for Ti (2.53 eV). Then, the solution energy decreases to 1.63 eV. Both Sn and Zr have lower solution energies than that of Ti. The highest positive solution energy is calculated for Ce.

Aliovalent Doping
Trivalent dopants (M = Ga, Sc, In, Y, Gd, and La) were substitutionally doped on the site. This process required charge compensating defects. Two possible such defects (Li interstitial or oxygen vacancy) were considered (refer to Equations (18) and (19)). Additional Li + ions in this materials would increase its capacity. The formation of oxygen vacancies would enhance the vacancy mediated oxygen ion diffusion. Furthermore, this doping process can lead to the loss of Li 2 O as Li interstitials and oxygen interstitials as facilitated by the formation oxygen vacancies can aggregate. Figure 5a reports the solution energies for the trivalent doping in which Li interstitials are the charge compensating defects (refer to Equation (18)). Ga is the most prominent dopant for this process, although its solution energy is 2.31 eV. The solution energy gradually increases with the increasing ionic radius. The dopant La has the highest solution energy of 5.29 eV. A similar trend is observed for the formation of oxygen vacancies, but with slightly high energies (refer to Figure 5b). Overall, Ga is the promising dopant on the Si site for both processes. dopant for this process, although its solution energy is 2.31 eV. The solution energy gradually increases with the increasing ionic radius. The dopant La has the highest solution energy of 5.29 eV. A similar trend is observed for the formation of oxygen vacancies, but with slightly high energies (refer to Figure 5b). Overall, Ga is the promising dopant on the Si site for both processes.

Incorporation of Li into LiAlSi2O6
Li incorporation was examined using DFT simulations. Here, we report the geometries, electronic structures, and volume change upon incorporation. The calculated total DOS plot of bulk LiAlSi2O6 is shown in Figure 6. The calculation shows that LiAlSi2O6 is an insulator with a band gap of 5.50 eV, which is in agreement with the band gap value of 5.575 eV calculated using DFT simulation by He et al. [47]. The total magnetic moment of the bulk LiAlSi2O6 is zero, as evidenced by the total DOS plot, in which spin up states are the mirror image of spin down states about the x-axis.

Incorporation of Li into LiAlSi 2 O 6
Li incorporation was examined using DFT simulations. Here, we report the geometries, electronic structures, and volume change upon incorporation. The calculated total DOS plot of bulk LiAlSi 2 O 6 is shown in Figure 6. The calculation shows that LiAlSi 2 O 6 is an insulator with a band gap of 5.50 eV, which is in agreement with the band gap value of 5.575 eV calculated using DFT simulation by He et al. [47]. The total magnetic moment of the bulk LiAlSi 2 O 6 is zero, as evidenced by the total DOS plot, in which spin up states are the mirror image of spin down states about the x-axis.
Li incorporation was examined using DFT simulations. Here, we report the geometries, electronic structures, and volume change upon incorporation. The calculated total DOS plot of bulk LiAlSi2O6 is shown in Figure 6. The calculation shows that LiAlSi2O6 is an insulator with a band gap of 5.50 eV, which is in agreement with the band gap value of 5.575 eV calculated using DFT simulation by He et al. [47]. The total magnetic moment of the bulk LiAlSi2O6 is zero, as evidenced by the total DOS plot, in which spin up states are the mirror image of spin down states about the x-axis. Incorporation of up to four Li atoms was considered. Figure 7 shows the relaxed structures. Incorporation of first two Li atoms is endoergic with respect to both the gas phase and bulk Li (see Table 6). Incorporation energy calculated using gas phase bulk Li metal is more positive than that calculated using the gas phase Li atom. This is because of the extra energy required to dissociate bulk Li to form Li atoms. Incorporation becomes exothermic for three and four Li atoms with respect to gas phase Li atom. This can be partly owing to the volume expansion introduced by the previous incorporation. Bader Incorporation of up to four Li atoms was considered. Figure 7 shows the relaxed structures. Incorporation of first two Li atoms is endoergic with respect to both the gas phase and bulk Li (see Table 6). Incorporation energy calculated using gas phase bulk Li metal is more positive than that calculated using the gas phase Li atom. This is because of the extra energy required to dissociate bulk Li to form Li atoms. Incorporation becomes exothermic for three and four Li atoms with respect to gas phase Li atom. This can be partly owing to the volume expansion introduced by the previous incorporation. Bader charge analysis [48] confirms that all Li atoms exhibit a +1 charge. Positively charged Li atoms strongly interact with negatively charged oxygen atoms in the lattice.
ChemEngineering 2021, 5, x FOR PEER REVIEW 9 of 12 charge analysis [48] confirms that all Li atoms exhibit a +1 charge. Positively charged Li atoms strongly interact with negatively charged oxygen atoms in the lattice.     Figure 8 shows the total DOS plots calculated for the Li-incorporated structures. Fermi energy level shifts towards conduction band upon incorporation. This is because of the presence of electrons that were produced during the Li incorporation. The resultant complexes exhibit a metallic character and Fermi levels are occupied by the states associated with Li according to the atomic DOSs plotted for incorporated Li (see Figure 8e-h).

Conclusions
Intrinsic disorder mechanisms are investigated in the simulated crystal model of α-LiAlSi2O6 utilizing static atomistic simulation techniques. The Al/Si cation anti-site disorder is found to be the most dominant defect in this material. The formation of Li vacancies is ensured by the Li-Frenkel defect process, which is higher in energy by 0.33 eV than the Al/Si anti-site defect. A slow diffusion of Li-ion is noted with high activation energy of 2.62 eV. Na, Ga, and Ge are the favorite dopants on the Li, Al, and Si sites, respectively. Here, we suggest that the concentration of both Li interstitials and O vacancies can be facilitated by doping of Ga on the Si site. Incorporation of lithium was studied using density functional theory simulations and it is shown that, at a high concentration, incorporation becomes easier and the insulating nature of pure LiAlSi2O6 becomes metallic upon incorporation.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1: Interatomic potential parameters used for dopant in the atomistic simulations of LiAlSi2O6.

Conclusions
Intrinsic disorder mechanisms are investigated in the simulated crystal model of α-LiAlSi 2 O 6 utilizing static atomistic simulation techniques. The Al/Si cation anti-site disorder is found to be the most dominant defect in this material. The formation of Li vacancies is ensured by the Li-Frenkel defect process, which is higher in energy by 0.33 eV than the Al/Si anti-site defect. A slow diffusion of Li-ion is noted with high activation energy of 2.62 eV. Na, Ga, and Ge are the favorite dopants on the Li, Al, and Si sites, respectively. Here, we suggest that the concentration of both Li interstitials and O vacancies can be facilitated by doping of Ga on the Si site. Incorporation of lithium was studied using density functional theory simulations and it is shown that, at a high concentration, incorporation becomes easier and the insulating nature of pure LiAlSi 2 O 6 becomes metallic upon incorporation.