Computational Investigation of a NASICON-Type Solid Electrolyte Material LiGe 2 (PO 4 ) 3

: Phosphate-based electrolyte materials are of great interest in the ﬁeld of Li-ion batteries due to their rigid structural integrity. LiGe 2 (PO 4 ) 3 is a NASICON-type phosphate material with high thermal and electrochemical stability. Computational simulation techniques were employed to study the defects, diffusion, and dopant properties of LiGe 2 (PO 4 ) 3 . Furthermore, the reaction energies for the formation of LiGe 2 (PO 4 ) 3 and the incorporation energies for the insertion of additional Li into this material were calculated. The calculations revealed that the Li Frenkel is the lowest-energy defect. The second most favorable defect is the Ge-P anti-site defect cluster. A low Li migration energy of 0.44 eV implies high Li ionic conductivity. The most favorable isovalent dopants on the Li and Ge sites are Na and Si, respectively. The formation of Li interstitials and oxygen vacancies can be facilitated through the doping of Ga on the Ge site. The doping of Ga slightly enhances the Li ionic conductivity. Li incorporation (up to four Li) is thermodynamically feasible.

Lithium germanium phosphate [LiGe 2 (PO 4 ) 3 ] is one of the NASICON-type solid electrolyte materials that has many advantages, such as high electrochemical and thermal stability [13,14]. The high stability is attributed to the stable oxidation state of germanium (IV) in LiGe 2 (PO 4 ) 3 . Practical application of this material is limited by its low ionic conductivity. Many experimental studies considered single or dual doping of trivalent cations on the Ge site to improve its Li-ion conductivity [15][16][17][18][19]. Such improvement is assisted via the formation of additional Li + ions in the lattice. Kotobuki et al. [20] used a co-precipitation method to prepare Li 1 . 5 Al 0 . 5 Ge 1 . 5 (PO 4 ) 3 at high temperatures and observed an enhancement in the Li-ion conductivity. A similar conclusion was reached in other experimental studies, although their methodologies of preparation were different. Nikodimos et al. [15] studied the Li-ion conductivity in Al-Sc dual-doped LiGe 2 (PO 4 ) 3 both experimentally and computationally and concluded that a significant enhancement in the ionic conductivity was achieved with the reduction of the activation energy.
Little is known about the defect chemistry of LiGe 2 (PO 4 ) 3 . Computer modeling techniques are useful for finding the possible intrinsic defects present in this material, Liion diffusion pathways and their activation energies, and solutions of a variety of dopants, as discussed in previous simulation studies [21][22][23][24][25][26]. To this end, we used a classical simulation based on the pair-wise potentials. Furthermore, density functional theory (DFT) simulations were employed to elucidate the electronic structures of LiGe 2 (PO 4 ) 3 , doped-LiGe 2 (PO 4 ) 3 , and Li-incorporated LiGe 2 (PO 4 ) 3 .

Computational Methods
The calculations were carried out by using classical pair potentials as implemented in the GULP (General Utility Lattice Program) code [27] and projected augmented wave (PAW) pseudo-potentials [28] as implemented in the VASP (Vienna ab initio Simulation Package) code [29]. The GULP code enabled us to examine the most probable defect energy processes, Li-ion migration pathways, and solutions of dopants. The electronic structures of pristine LiGe 2 (PO 4 ) 3 and doped-LiGe 2 (PO 4 ) 3 were modeled by using the VASP code.
The GULP code uses ionic interactions, such as long-range attraction (Coulomb) and short-range (Pauli repulsion and van der Waals attraction) forces. The short-range interactions were described using Buckingham potentials (see Table 1) [30,31]. Both lattice constants and ionic positions were relaxed with the aid of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [32]. Point defect relaxation was modeled using the Mott-Littleton method [33]. Li-ion diffusion pathways were calculated by considering seven interstitial positions between two local Li vacancy sites. The activation energy was calculated as the energy difference between the local maximum energy along this diffusion pathway and the initial vacancy defect energy. The VASP code uses PAW potentials and plane wave basis sets. The exchangecorrelation term was described using the generalized gradient approximation (GGA) as formulated by Perdew, Burke, and Ernzerhof (PBE) [34]. In all cases, a plane wave basis set with a cut-off value of 500 eV was used. A 4 × 4 × 4 Monkhorst-Pack [35] k-point mesh was used to model bulk LiGe 2 (PO 4 ) 3 . A 2 × 2 × 1 supercell consisting of 432 atoms (Li: 24, Ge: 48, P: 72, and O: 288) was used to model defects. For defect calculations, we used a 2 × 2 × 2 Monkhorst-Pack k-point mesh. The conjugate gradient algorithm [36] was applied to optimize both bulk and defect configurations. The force tolerance was set to be 0.001 eV/Å. The Bader charge analysis [37] was used to estimate the charges on the dopants and incorporated Li atoms.
The two-body expression [Φ ij (r ij ) = A ij exp (−r ij/ ρ ij ) − C ij /r ij 6 ] was used, where A, ρ, and C are parameters reproducing the experimental data. The values of Y and K are shell charges and spring constants, respectively.

Crystal Structure of LiGe 2 (PO 4 ) 3
LiGe 2 (PO 4 ) 3 crystallizes in the rhombohedral space group R3C (a = b = 8.275 Å, c = 20.470 Å, α = β = 90 • and γ = 180 • ) [38] and consists of PO 4 tetrahedra and GeO 6 octahedra. These two polyhedrons share their corners, forming a three-dimensional network (see Figure 1a). Both classical and DFT simulations were employed to relax the ionic positions and lattice parameters of LiGe 2 (PO 4 ) 3 . The calculated lattice parameters were in good agreement with the values reported in the experiment, showing the suitability of the potential parameters, pseudopotentials, and basis sets used in this study (see Table 2). Sustain. Chem. 2022, 3, FOR PEER REVIEW 3 The electronic structure calculation performed by using DFT showed that LiGe2(PO4)3 is an insulator with a band gap of 3.50 eV (see Figure 1b). The value calculated in this study is in a good agreement with a value (3.38 eV) reported in a previous DFT simulation study [39].

Defect Properties
In this section, the formation of intrinsic defect processes such as Frenkel, Schottky, and anti-site defects, is discussed. Initially, the point defect energies were calculated, and then they were combined to calculate the intrinsic defect processes with the appropriate energies of the lattices. The following equations were written using Kröger-Vink notation [40] to explain the defect reactions. Anti-site defects were considered to be in an isolated or clustered form. Isolated anti-site defects were modeled by considering defects separately. In the clustered form, defects were considered in the same supercell.  The electronic structure calculation performed by using DFT showed that LiGe 2 (PO 4 ) 3 is an insulator with a band gap of 3.50 eV (see Figure 1b). The value calculated in this study is in a good agreement with a value (3.38 eV) reported in a previous DFT simulation study [39].

Defect Properties
In this section, the formation of intrinsic defect processes such as Frenkel, Schottky, and anti-site defects, is discussed. Initially, the point defect energies were calculated, and then they were combined to calculate the intrinsic defect processes with the appropriate energies of the lattices. The following equations were written using Kröger-Vink notation [40] to explain the defect reactions. Anti-site defects were considered to be in an isolated or clustered form. Isolated anti-site defects were modeled by considering defects separately. In the clustered form, defects were considered in the same supercell.
Li/Ge antisite (isolated) : Li/Ge antisite (cluster) : Ge/P antisite (isolated) : Ge X Ge + P X P → Ge P + P • Ge (10) Ge/P antisite (cluster) : Ge X Ge + P X P → Ge P : P • Ge X (11) Table 3 reports the calculated defect energies. The calculations show that the most prominent defect in this material was the Li Frenkel defect (0.75 eV/defect). The concentration of this defect was higher than those of the other defect processes. In particular, Li vacancies in the material facilitate vacancy-assisted Li-ion migration. The second most favorable defect was the Ge-P anti-site defect cluster (1.26 eV/defect). This defect cluster was formed by the aggregation of isolated defects with a binding energy of −0.75 eV. Anti-site defects have been observed in many as-prepared oxide materials [41][42][43]. Both experimental and theoretical studies have shown that this defect would change the properties of materials [44][45][46]. The Li-Ge anti-site defect energy was high due to the charge difference between Li + and Ge 4+ . The Li 2 O partial Schottky energy was relatively low (3.10 eV/defect). Other Frenkel and Schottky defect energies were high (>4 eV), meaning that the defects were not significant in this material at room temperature.

Diffusion of Li Ions
Materials that exhibit a high Li-ionic conductivity are important in the construction of promising Li-ion batteries. The experimental determination of long-range Li-ion diffusion pathways is generally challenging. Computational simulation studies have been useful in predicting the diffusion pathways or validating the pathways determined in experimental studies. For example, the Li-ion diffusion pathway and its migration barrier in Li 1+x Al x Ti 2−x (PO 4 ) 3 as calculated by Case et al. [47] were in excellent agreement with the corresponding experimental observations [19,48].
Here, a classical simulation was employed to examine the diffusion of Li ions. We identified a local Li hop with a jump distance of 5.89 Å and calculated its vacancy-assisted migration pathway and activation energy. The activation energy for this hop was 0.44 eV. This Li hop was extended through the lattice to construct a long-range Li-ion diffusion pathway. This resulted in a three-dimensional pathway, as shown in Figure 2a. The energy profile diagram (see Figure 2b) shows the activation energy calculated for this local Li hop. As this long-range diffusion pathway consisted of only one Li hop, the overall migration barrier was also 0.44 eV. The experimental activation energy of the Li ions in LiGe 2 (PO 4 ) 3 reported by Martínez-Juárez et al. [49] was 0.60 eV. In another experimental study by Nikodimos et al. [15], it was reported that the activation energy was 0.286 eV, and the DFT-based calculation by the same authors revealed an activation energy of 0.329 eV. The activation energy calculated in this study was between two experimental values and was closer to the value calculated in the DFT simulation. Nevertheless, this material exhibited high Li-ion conductivity. and the DFT-based calculation by the same authors revealed an activation energy of 0.329 eV. The activation energy calculated in this study was between two experimental values and was closer to the value calculated in the DFT simulation. Nevertheless, this material exhibited high Li-ion conductivity.

Solutions of Dopants
Ionic dopants are generally used to increase the ionic conductivity and chemical and mechanical stability of materials. Here, we considered a range of dopants to screen and predicted promising dopants that can be of interest for future experiments. For the dopant calculations, a classical simulation was employed. Charge-compensating defects, such as interstitials, vacancies, and lattice energies, were used to calculate solution energies. In the electronic supplementary information (ESI), the Buckingham potentials used for the oxides of dopants are listed (see Table S1).
First, alkali dopants (M = Na, K and Rb) on the Li site were considered. The following reaction equation explains the doping process.
The most favorable dopant was Na + , with an exothermic solution energy of −0.90 eV (Table 4). This was due to the ionic radius of Na + , which is close to that of Li + (0.76 Å). There was a gradual increase in the solution energy with the increase in the ionic radii of the dopants. Both K + and Rb + ions exhibited endoergic solution energies. A possible doped configuration that can be synthesized is Li1-xNaxGe2(PO4)3 (x = 0.0-1.0). The relaxed structure of Na-doped LiGe2(PO4)3 is shown Figure 3a. The doped Na atom forms a distorted octahedral structure (NaO6) with longer Na-O bond distances than Li-O bond distances in the LiO6 octahedral structure. As mentioned earlier, this is due to the fact that the ionic radius of Na + is larger than that of Li + . The DOS plots show that the doping had a negligible impact on the electronic structure and that there were no additional states present in the gap.

Solutions of Dopants
Ionic dopants are generally used to increase the ionic conductivity and chemical and mechanical stability of materials. Here, we considered a range of dopants to screen and predicted promising dopants that can be of interest for future experiments. For the dopant calculations, a classical simulation was employed. Charge-compensating defects, such as interstitials, vacancies, and lattice energies, were used to calculate solution energies. In the electronic supplementary information (ESI), the Buckingham potentials used for the oxides of dopants are listed (see Table S1).
First, alkali dopants (M = Na, K and Rb) on the Li site were considered. The following reaction equation explains the doping process.
The most favorable dopant was Na + , with an exothermic solution energy of −0.90 eV (Table 4). This was due to the ionic radius of Na + , which is close to that of Li + (0.76 Å). There was a gradual increase in the solution energy with the increase in the ionic radii of the dopants. Both K + and Rb + ions exhibited endoergic solution energies. A possible doped configuration that can be synthesized is Li 1-x Na x Ge 2 (PO 4 ) 3 (x = 0.0-1.0). The relaxed structure of Na-doped LiGe 2 (PO 4 ) 3 is shown Figure 3a. The doped Na atom forms a distorted octahedral structure (NaO 6 ) with longer Na-O bond distances than Li-O bond distances in the LiO 6 octahedral structure. As mentioned earlier, this is due to the fact that the ionic radius of Na + is larger than that of Li + . The DOS plots show that the doping had a negligible impact on the electronic structure and that there were no additional states present in the gap. Sustain. Chem. 2022, 3, FOR PEER REVIEW 6 The most promising dopant was Ga 3+ , and the solution energy for this dopant was 2.20 eV (see Table 5). The solution energy calculated for Al 3+ was higher than that calculated for Ga 3+ by only 0.14 eV. The favorability of both dopants (Ga 3+ and Al 3+ ) was due to the fact that their ionic radii matched closely with the ionic radius of Ge 4+ (0.53 Å). The solution energy gradually increased with the increase in the ionic radius. The most unfavorable dopant was La 3+ , and its solution energy was 4.28 eV. The same doping process can introduce oxygen vacancies, as explained in the following reaction equation.
The solution energies were higher for this process than those calculated for the process in which Li interstitials acted as charge-compensating defects (see Equation (13)). This indicates that the concentration of oxygen vacancies will not be high at room temperature. Although the solution energies were larger than 5 eV for all dopants, the candidate dopant was Ga 3+ . The trend in the solution energy was the same as that calculated for the previous process.
The relaxed structures of Ga-doped LiGe2(PO4)3 and electronic structures are shown in Figure 4. The Ga-O bond distances in the GaO6 octahedral unit were larger than the Ge-O bond distances. The band-gap value was not altered much upon doping. The atomic The most promising dopant was Ga 3+ , and the solution energy for this dopant was 2.20 eV (see Table 5). The solution energy calculated for Al 3+ was higher than that calculated for Ga 3+ by only 0.14 eV. The favorability of both dopants (Ga 3+ and Al 3+ ) was due to the fact that their ionic radii matched closely with the ionic radius of Ge 4+ (0.53 Å). The solution energy gradually increased with the increase in the ionic radius. The most unfavorable dopant was La 3+ , and its solution energy was 4.28 eV. The same doping process can introduce oxygen vacancies, as explained in the following reaction equation.
The solution energies were higher for this process than those calculated for the process in which Li interstitials acted as charge-compensating defects (see Equation (13)). This indicates that the concentration of oxygen vacancies will not be high at room temperature. Although the solution energies were larger than 5 eV for all dopants, the candidate dopant was Ga 3+ . The trend in the solution energy was the same as that calculated for the previous process.
The relaxed structures of Ga-doped LiGe 2 (PO 4 ) 3 and electronic structures are shown in Figure 4. The Ga-O bond distances in the GaO 6 octahedral unit were larger than the Ge-O bond distances. The band-gap value was not altered much upon doping. The atomic DOS plots show that states arising from Ga were mainly localized in the valence band (see Figure 4c,f).
Chem. 2022, 3, FOR PEER REVIEW 7 DOS plots show that states arising from Ga were mainly localized in the valence band (see Figure 4c,f). Finally, tetravalent dopants (Si, Ti, Sn, Zr, and Ce) were substitutionally doped on the Ge site. The solution energies were calculated using the following equation.

MO
Ge → M GeO (15) The solution energy calculated for Si 4+ was 0.67 eV (see Table 6). The other dopants exhibited higher solution energies. Thus, Si 4+ was the most promising. This was partly due to the fact that the ionic radius of Si 4+ closely matched with that calculated for Ge 4+ . The next most favorable dopant was Ti 4+ . Its solution energy was almost double the solution energy calculated for Si 4+ . Both Sn 4+ and Zr 4+ had solution energies that were closer to each other. The solution energy calculated for Ce 4+ was highly endoergic, and the doping of this cation required a high temperature. The relaxed structure of Si-doped LiGe2(PO4)3 and its electronic structure are shown in Figure 5. The Si-O bond distances in SiO6 unit are shorter than Ge-O bond distances in GeO6 unit. This is partly owing to the fact that the ionic radius of Si 4+ (0.40 Å) is smaller than that of Ge 4+ (0.53 Å). The total and atomic DOS plots show that the insulating nature of this material is unchanged. Finally, tetravalent dopants (Si, Ti, Sn, Zr, and Ce) were substitutionally doped on the Ge site. The solution energies were calculated using the following equation.
The solution energy calculated for Si 4+ was 0.67 eV (see Table 6). The other dopants exhibited higher solution energies. Thus, Si 4+ was the most promising. This was partly due to the fact that the ionic radius of Si 4+ closely matched with that calculated for Ge 4+ . The next most favorable dopant was Ti 4+ . Its solution energy was almost double the solution energy calculated for Si 4+ . Both Sn 4+ and Zr 4+ had solution energies that were closer to each other. The solution energy calculated for Ce 4+ was highly endoergic, and the doping of this cation required a high temperature. Table 6. Solution energies calculated for tetravalent dopants on the Ge site. The relaxed structure of Si-doped LiGe 2 (PO 4 ) 3 and its electronic structure are shown in Figure 5. The Si-O bond distances in SiO 6 unit are shorter than Ge-O bond distances in GeO 6 unit. This is partly owing to the fact that the ionic radius of Si 4+ (0.40 Å) is smaller than that of Ge 4+ (0.53 Å). The total and atomic DOS plots show that the insulating nature of this material is unchanged.

Diffusion of Li Ions in Ga-Doped LiGe2(PO4)3
The migration energy was calculated for a local Li hop in the presence of Ga on the Ge site. The migration energy was calculated to be 0.39 eV. The dopant concentration was 2.08%. It is expected that the migration energy can be decreased with further substitution.

Li Incorporation
The incorporation of extra lithium was considered by calculating the incorporation energies with respect to a Li atom and bulk Li. Table 7 reports the incorporation energies, Bader charges on the incorporated Li atoms, and volume changes upon Li incorporation. The incorporation energies were negative in all cases with respect to both a Li atom and bulk Li as reference states. The incorporation energies calculated with respect to the bulk Li were less exothermic. This was partly due to the additional energy needed to dissociate bulk Li to form Li atoms. The Bader charge analysis showed that all incorporated Li atoms became single positively charged ions (Li + ). The formation of positively charged ions was due to the strong Li-O bond formation with the nearest neighboring oxygen ions in the lattice (see Figure 6a-d). The incorporation became stronger with increase in the concentration of Li atoms. This enhancement in the incorporation was facilitated by the expansion of the volume. DOS plots were constructed for the Li-incorporated structures. The gap region consisted of states associated with the incorporated Li, as evidenced by the atomic DOS plots (see Figure 6i-l). Incorporation slightly reduced the band gap without changing the insulating character of this material. The migration energy was calculated for a local Li hop in the presence of Ga on the Ge site. The migration energy was calculated to be 0.39 eV. The dopant concentration was 2.08%. It is expected that the migration energy can be decreased with further substitution.

Li Incorporation
The incorporation of extra lithium was considered by calculating the incorporation energies with respect to a Li atom and bulk Li. Table 7 reports the incorporation energies, Bader charges on the incorporated Li atoms, and volume changes upon Li incorporation. The incorporation energies were negative in all cases with respect to both a Li atom and bulk Li as reference states. The incorporation energies calculated with respect to the bulk Li were less exothermic. This was partly due to the additional energy needed to dissociate bulk Li to form Li atoms. The Bader charge analysis showed that all incorporated Li atoms became single positively charged ions (Li + ). The formation of positively charged ions was due to the strong Li-O bond formation with the nearest neighboring oxygen ions in the lattice (see Figure 6a-d). The incorporation became stronger with increase in the concentration of Li atoms. This enhancement in the incorporation was facilitated by the expansion of the volume. DOS plots were constructed for the Li-incorporated structures. The gap region consisted of states associated with the incorporated Li, as evidenced by the atomic DOS plots (see Figure 6i-l). Incorporation slightly reduced the band gap without changing the insulating character of this material.

Conclusions
This work focused on the computational characterization of defects, diffusion, solutions of dopants, reaction energies for the formation of LiGe2(PO4)3, and the energetics of the incorporation of Li into LiGe2(PO4)3. The Li Frenkel defect was found to be the most favorable defect in this material, indicating that the concentration of Li vacancies would be higher than those of the other vacancies. The diffusion of Li-ions is three dimensional and fast, with an activation energy of 0.44 eV. The Li and Ge sites were favorably doped by Na and Si, respectively. Doping of Ga on the Ge site can produce Li interstitials or oxygen vacancies as charge-compensating defects. There is a slight improvement in the diffusion of Li + ions upon Ga doping on the Ge site. Incorporation of Li is exoergic and becomes more exothermic with the concentration of Li.

Conclusions
This work focused on the computational characterization of defects, diffusion, solutions of dopants, reaction energies for the formation of LiGe 2 (PO 4 ) 3 , and the energetics of the incorporation of Li into LiGe 2 (PO 4 ) 3 . The Li Frenkel defect was found to be the most favorable defect in this material, indicating that the concentration of Li vacancies would be higher than those of the other vacancies. The diffusion of Li-ions is three dimensional and fast, with an activation energy of 0.44 eV. The Li and Ge sites were favorably doped by Na and Si, respectively. Doping of Ga on the Ge site can produce Li interstitials or oxygen vacancies as charge-compensating defects. There is a slight improvement in the diffusion of Li + ions upon Ga doping on the Ge site. Incorporation of Li is exoergic and becomes more exothermic with the concentration of Li.