Electrode and Electrolyte Materials From Atomistic Simulations: Properties of LixFEPO4 Electrode and Zircon-Based Ionic Conductors

LixFePO4 orthophosphates and fluorite- and pyrochlore-type zirconate materials are widely considered as functional compounds in energy storage devices, either as electrode or solid state electrolyte. These ceramic materials show enhanced cation exchange and anion conductivity properties that makes them attractive for various energy applications. In this contribution we discuss thermodynamic properties of LixFePO4 and yttria-stabilized zirconia compounds, including formation enthalpies, stability, and solubility limits. We found that at ambient conditions LixFePO4 has a large miscibility gap, which is consistent with existing experimental evidence. We show that cubic zirconia becomes stabilized with Y content of ~8%, which is in line with experimental observations. The computed activation energy of 0.92eV and ionic conductivity for oxygen diffusion in yttria-stabilized zirconia are also in line with the measured data, which shows that atomistic modeling can be applied for accurate prediction of key materials properties. We discuss these results with the existing simulation-based data on these materials produced by our group over the last decade. Last, but not least, we discuss similarities of the considered compounds in considering them as materials for energy storage and radiation damage resistant matrices for immobilization of radionuclides.


INTRODUCTION
MPO 4 orthophosphates, with M cations being transition metals, rare-earth elements or actinides, are ceramic materials of interest in various research fields, including geochronology (Williams et al., 2007), geothermometry (Andrehs and Heinrich, 1998;Mogilevsky, 2007), energy storage (Iyer et al., 2006;Yamada et al., 2006;Dunn et al., 2015;Dong et al., 2017;Cerdas et al., 2018;Phan et al., 2019), and nuclear waste management (Ewing and Wang, 2002;Neumeier et al., 2017a;Schlenz et al., 2018), to name but a few. Most of the potential applications come from high durability (e.g., radiation damage resistance) of these materials (Neumeier et al., 2017a;Phan et al., 2019). There exist large varieties of phosphate-based ceramics of different crystalline structures (e.g., cheralite, apatites, olivine, kosnarite, see Iyer et al., 2006;Neumeier et al., 2017a;Phan et al., 2019). The orthophosphates of interest for energy storage, FePO 4 and LiFePO 4 , have olivine-type structure of orthorombic space group symmetry of Pnma in which Fe exists in an octahedral environment (Iyer et al., 2006;Maxisch and Ceder, 2006, see Figure 1 for visualization of the structures). It consists of eight-fold coordinated Fe 3+ cations, which become reduced to Fe 2+ upon incorporation of Li. Besides, there exist less stable FePO 4 monoclinic, orthorombic, and trigonal phases with Fe in tetrahedraly coordinated crystallographic positions (Iyer et al., 2006;Yamada et al., 2006;Dong et al., 2017). LiFePO 4 materials are one of the best known candidates for energy storage electrodes (Dunn et al., 2015;Cerdas et al., 2018). These exhibit high theoretical energy density of 170 mAh/g and voltage of 3.5 V (Phan et al., 2019). Because of this, these materials are considered for large scale energy storage devices, including batteries for hybrid and electric automobiles (Prosini et al., 2002;Dunn et al., 2015;Cerdas et al., 2018). Intercalation of Li ions into FePO 4 results in formation of a solid solution between Li cations and vacant sites (Phan et al., 2019). Depending on the sizes of mixing species, a solid solution compounds could form a thermodynamically stable solid solution or a compound with a mixture of two phases, each rich in one of the cations Ji et al., 2019a). The second case indicates formation of a temperature dependent miscibility gap. The two endmember phases may be of the same type, like in the case of Li x FePO 4 , or different phases, like in the monazite-xenotime system (Mogilevsky, 2007;Ji et al., 2019a). The formation of miscibility gap is correlated with the maximum solubility (Ji et al., 2019a). Such a temperature-dependent maximum solubility offers opportunity to use such a system as geothermometer (Andrehs and Heinrich, 1998;Mogilevsky, 2007). It is known from experiment that Li x FePO 4 solid solution posses a wide miscibility gap at ambient temperature (Yamada et al., 2006;Meethong et al., 2007;Phan et al., 2019). In particular, studies of Yamada et al. (2006) indicate relatively low Li solubility limits in Li x FePO 4 system at x = 0.05 and 0.89. The thermodynamics of such a system, however, although has been modeled by CALPHAD method (Phan et al., 2019), has not been modeled in details using ab initio atomistic modeling methods and the results of various studies differ substantially (Phan et al., 2019).
In the last two decades, atomistic modeling became a widely used research technique in various research fields, including energy materials (Chroneos et al., 2013;Jahn and Kowalski, 2014;Wu et al., 2019). We used it intensively over the past decade for computation of, for instance, various physical and chemical properties of orthophosphate-and zircon-based ceramics (e.g., Kowalski et al., 2015;Ji et al., 2019a). This is because steady advancements in high performance computing and computational software, especially in ab initio methodsbased codes, allows nowadays for computation of complex systems containing hundreds of atoms from first principles (Jahn and Kowalski, 2014). Regarding ceramic compounds considered here, computational studies have been used to deliver information on: the structural (Rustad, 2012;Feng et al., 2013;Blanca-Romero et al., 2014;Beridze et al., 2016;Huittinen et al., 2017), the electronic structure (Tang and Holzwarth, 2003;Blanca-Romero et al., 2014;Kowalski et al., 2017a;Lee et al., 2017), the elastic (Wang et al., 2005;Feng et al., 2013;Ali et al., 2016;Ji et al., 2017a;Kowalski et al., 2017b), the thermodynamic (Mogilevsky, 2007;Feng et al., 2013;Li et al., 2014;Kowalski et al., 2015Ji et al., 2017b;Neumeier et al., 2017b;Eremin et al., 2019), the thermochemical (Rustad, 2012;Beridze et al., 2016;Kowalski, 2020), the electrochemical (Krishnamurthy et al., 2005b;Lee et al., 2017), and the radiation damage resistance Ji et al., 2017c;Jolley et al., 2017) parameters as well as materials at high-pressure (López-Solano et al., 2010;Stavrou et al., 2012;Ali et al., 2016;Shein and Shalaeva, 2016;Gomis et al., 2017). The relevant research activity increases steadily worldwide, with most of the papers published just recently. One important aspect is the correct calculations of compounds with d-and f -elements that contain strongly correlated electrons [e.g., Fe, Ni, lanthanides (Ln), actinides (An)]. In a series of papers we have shown that standard DFT approach often fails for such cases and these compounds must be carefully computed, including proper accounting for the correlation effects Blanca-Romero et al., 2014;Kowalski et al., 2015;Li and Kowalski, 2018). These simulations must be performed with methods beyond the standard DFT+U approach and include derivation of the Hubbard U parameter and careful choice of projectors for estimation of occupancy of d− and f − levels within the DFT+U scheme Kvashnina et al., 2018;Kick et al., 2019). In particular, we apply the linear response method (Cococcioni and de Gironcoli, 2005) with Wannier orbitals as representation of d or f states (Kvashnina et al., 2018) and here we will demonstrate impact of these procedures on the estimation of formation enthalpies and solubility limits of Li x FePO 4 compound.
In this contribution we provide an overview of the recent atomistic modeling activities on the orthophosphates and zirconates, focusing on the information that have been delivered by atomistic modeling activities at Forschungszentrum Jülich and that allowed on many occasions for better characterization of these materials, including long-term thermodynamic stability, thermochemical parameters, and thermal conductivity. Besides such overview, we present results of computation of thermochemical and thermodynamic parameters of Li x FePO 4 solid solution, with focus on the formation of miscibility gap in this system, as well as simulation of yttria-stabilized zirconia with focus on prediction the phase stability and ionic conductivity in this class of materials. We especially highlight a cross-linking, interdisciplinary character of our research, from which the general science community could highly benefit.

COMPUTATIONAL APPROACH
In all ab initio 1 calculations discussed here we used a density functional theory (DFT)-based quantum chemistry approach and calculations were performed with Quantum-ESPRESSO simulation package (Giannozzi et al., 2009). We applied the PBEsol exchange-correlation functional (Perdew et al., 2008), the ultrasoft pseudopotentials to represent the core electrons of the atoms (Vanderbilt, 1990) and the plane-wave energy cutoff of 50 Ryd. The PBEsol functional is specifically selected because it correctly reproduces slowly varying electron density limit and results in good structural parameters of solids (Perdew et al., 2008). This is important for consideration of, for instance, thermodynamics of solid solutions (Li et al., 2014;Ji et al., 2019a). Following our broad experience on computation of lanthanide orthophosphates and zirconates (e.g., Blanca-Romero et al., 2014;Li et al., 2015;Beridze et al., 2016) we applied the self-consistent DFT+U approach. The DFT+U calculations were performed with the Hubbard U parameter values computed from first principles using the linear response method of Cococcioni and de Gironcoli (2005). In order to apply realistic projectors for occupations of d orbitals of Fe we used Poor Man Wannier scheme implemented in Quantum-ESPRESSO. This computational setup was extensively tested by us in several studies and, among others, proved to give very good results for orthophosphates and zirconates Beridze et al., 2016;Finkeldei et al., 2017). The activation barriers were computed using Nudged Elastic Band (NEB) method implemented in Quantum-ESPRESSO, with 10 images and climbing image procedure to compute the transition state. The exponential pre-factors were computed with the aid of transition state theory (Moynihan et al., 1982) and probability of vacancy migration derived by Bukaemskiy et al. (2021), so that the ionic conductivity, σ , is estimated using the modified here formula of Moynihan et al. (1982): where Ze is the charge of the carrier, k b is the Boltzmann constant, N a is Avogadro number, d is the distance of the jump, E a is the activation barrier, T is the temperature and ν 0 is the attempt frequency. The first part in the square bracket reflects the vacancy migration probability contribution to the pre-exponential factor derived by Bukaemskiy et al. (2021). The attempt frequency ν 0 was estimated from the computation of phonon spectra of the initial (IS) and transition states (TS), as The computation of Li x FePO 4 phases were performed with 2x2x1 supercells (96 atoms for FePO 4 and 112 atoms for LiFePO 4 phases) using the 2x2x2 k-point grid. The oxides were computed as: Fe 2 O 3 (Pnmm symmetry, with supercell containing 30 atoms and 4x6x2 k-point grid), hexagonal P 2 O 5 as is in Blanca-Romero et al. (2014) and Y 2 O 3 as cubic oxide. The magnetic arrangements in iron phases were computed using the models of Whittingham et al. (2005) and Lee et al. (2017). The Hubbard U parameters computed with the linear response method are listed in Table 1. As in our previous studies Beridze et al., 2016;Kvashnina et al., 2018;Sun et al., 2020) we see strong dependence on the Fe redox state, with the U parameter for Fe(III) being ∼0.5 eV larger than for Fe(II). This is well-consistent with previous studies of FePO 4 and LiFePO 4 phases by Maxisch and Ceder (2006), who obtained 4.9 and 3.7 eV for both phases (taken as effective value, U eff = U − J), respectively.

Structural Data
The first test of a computational method is its ability to reproduce the measured lattice parameters of the computed crystalline solid. The lattice parameters of considered materials were measured by different studies and are well-known. These are collected in Table 2 and compared to the computed data.
Our previous studies showed that the structural parameters of lanthanide-orthophosphates are very sensitive to the applied computational method, especially to the exchange-correlation functional . A correct treatment of strongly correlated 4f electrons also plays an important role in those cases. Here, with our computational setup we got much better fit to the measured data than Maxisch and Ceder (2006), who applied the PBE exchange-correlation functional that tends to overestimate lattice parameters and volumes, which is also evident in the data collected in Table 2. Blanca-Romero et al. (2014) performed extensive tests of the capability of different DFT-based approaches to reproduce the measured lattice parameters and bond-distances of monazitetype lanthanide-orthophosphates. In that paper, we found that the standard DFT method with explicitly computed 4f electrons overestimates the lattice parameters and bond-lengths by up to 3%, which is consistent with previous studies (Wang et al., 2005;Rustad, 2012). A very good match to the experimental values of structural parameters was achieved applying the parameter free DFT+U method, with the PBEsol exchangecorrelation functional (Perdew et al., 2008) and the Hubbard U parameters derived ab initio. An improved description of structures by the PBEsol exchange-correlation functional is an important property of that functional. It recovers the known solution for slowly varying electron densities and with that it improves the description of structural parameters over widely used GGA functionals such as PBE (see discussion by Perdew et al., 2008). It is evident also for zirconate-based materials.
In our past studies of Nd 2−x Zr 2+x O 7+x/2 compound, with the PBEsol exchange-correlation functional we got perfect match to the lattice parameter in case of the pyrochlore phase (Finkeldei et al., 2017). In the most recent contribution we computed the MUO 4 compounds with M = Ni, Fe, Co, Cd, and found that only by using the correct projectors for estimation of occupation of d orbitals, e.g., Wannier functions, we could reproduce experimentally seen structural distortions (Murphy et al., 2021). The problem arises from the fact that with the standard DFT+U approach, when using atomic orbitals as projectors, the total occupancy of the d or f states of interest is much higher than the actual one (Kick et al., 2019;Murphy et al., 2021). As illustrated in Table 3, in our case it gives ∼1.3 excess electrons for Fe(III) and ∼0.5 for Fe(II). This is cured when using Wannier orbitals as projectors. Unfortunately, the forces and cell optimizations with this scheme are not yet implemented in Quantum-ESPRESSO or any equivalent codes, so we could not perform geometry optimization with such a more realistic approach.

Formation Enthalpies
The formation enthalpies from oxides for series of lanthanide orthophosphates have been measured by Ushakov et al. (2001) and of FePO 4 and LiFePO 4 by Iyer et al. (2006). Rustad (2012) noticed that there is a systematic offset between the computed and measured values for LnPO 4 of ∼40 kJ/mol, with the computed enthalpies being less exothermic. Blanca-Romero et al. (2014) have shown that this offset is present also in the DFT+U calculations and is to a large extent Ln-cation independent, which rules out the 4f electrons correlations as a contributing factor. They attributed this to the overestimation of P-O bond lengths, and thus volumes, of the LnPO 4 and P 2 O 5 compounds. Beridze et al. (2016) have found an identical offset for xenotime phase. When a constant shift of ∼30 kJ/mol is applied to the computed formation enthalpies, the measured values are nicely reproduced. The computed formation enthalpies for FePO 4 and LiFePO 4 compounds are reported in Table 4. It is evident that DFT heavily underestimates the values (taking absolute values) by 40-100 kJ/mol. This is significantly improved with the DFT+U approach with an error of 20 kJ/mol. Best result, however, is obtained with the DFT+U method when Wanner functions are used as projectors of Fe d states occupations. In this scheme, the computed formation enthalpies are within 10 kJ/mol. Accurate prediction of formation enthalpies are crucial for correct estimate of thermodynamic parameters of solid solution, including solubilities (section 3.4).
The most stable phase of ZrO 2 is monoclinic. However, upon doping with tri-valent elements it undergoes phase transition to cubic phase, with possible triclinic phase as an intermediate. The experimental evidence shows transition to that phase at The energies are reported in kJ/mol.
∼8% of YO 1.5 (Lee et al., 2003;Götsch et al., 2016;Ahamer et al., 2017). In Figure 2, we show the results of computation of formation enthalpies for the three phases of ZrO 2 (monoclinic, tetragonal, and cubic). The results are plotted together with the experimental data of Lee et al. (2003). The computed values show clearly that at ∼ 8% content of Y the enthalpy of cubic phase becomes the lowest and that phase most stable, which is well consistent with the aforementioned experimental data. Moreover, the computed formation enthalpy as a function of Y content is well consistent with the measured values. Also the computed enthalpy difference between monoclinic and cubic phases of ∼14 kJ/mol is well consistent with the previous measurements and estimates (ranging from 6 to 22 kJ/mol, with the best measured value of 10 kJ/mol, Lee et al., 2003).

Elastic, Thermodynamic, and Thermal Conductivity Parameters
Besides formation enthalpies, our previous studies show good ability of atomistic modeling to predict the elastic, thermodynamic, and thermal parameters of considered systems. Li et al. (2014), , and Kowalski et al. (2017b) computed the elastic parameters of series of lanthanide orthophosphates and obtained good agreements with the available experimental data. With these computed data they provided good estimates of parameters that are key for modeling of solid solutions within these compounds. Maxisch and Ceder (2006) computed elastic parameters of FePO 4 and LiFePO 4 compounds, which we will use for estimates of thermodynamic parameters of solid solutions in Li x FePO 4 system (section 3.4). In the follow-up studies of lanthanide phosphate we computed heat capacities for series of these compounds and explained the quasirandom-like behavior of heat capacity along lanthanide series by the lanthanide cation-dependent contribution from thermal excitation of 4f electrons (Schottky effect, Kowalski et al., 2015Kowalski et al., , 2017bJi et al., 2017a). Ji et al. (2019b) computed thermal conductivity for series of lanthanide phosphates and were able to derive accurate information on phonon mean free path in these systems. These studies demonstrate that with appropriate computational setup, accurate information on various physical parameters can be delivered by atomistic simulations and materials effectively screened for desired physical characteristics. This has been used by us to deliver crucial information on orthophosphate-based ceramics for immobilization of actinides (Huittinen et al., 2017(Huittinen et al., , 2018Ji et al., 2019a). Besides computing orthophosphate-type ceramics, we derived various parameters for fluorite-and pyrochlore-type compounds. Among those, with atomistic modeling we computed defect FIGURE 2 | The computed formation enthalpies of YO 1.5 − ZrO 2 system. The lines represent results for monoclinic (dashed green), trigonal (dot-dashed red), and cubic (solid black) phases. The data comes from Lee et al. (2003).

Solid Solutions and Solubility
Formation and thermodynamic stability of solid solutions is a topic of studies in various research fields, including battery research (Meethong et al., 2007;Phan et al., 2019). In the past we published a series of studies on lanthanide phosphate ceramic-type solid solutions in the context of geothermometry and nuclear waste management (Li et al., 2014;Hirsch et al., 2017;Kowalski et al., 2017b;Neumeier et al., 2017b;Ji et al., 2019a). This is because the information gained allows for the assessment of long term stability of ceramic nuclear waste forms against phase separation. It was shown experimentally (Popa et al., 2007;Phan et al., 2019) and by ab initio simulations (Li et al., 2014;) that single phase orthophosphate-based solid solutions are highly regular. The excess enthalpy of mixing, H E , of a A 1−x B x PO 4 mixed cation compound could be described by a simple equation (Popa et al., 2007): where W is a Margules interaction parameters (Prieto, 2009). A solid solution is stable against formation of a miscibility gap if W < 2RT, where R is the gas constant. It has been demonstrated experimentally (Yamada et al., 2006;Phan et al., 2019) that Li x FePO 4 solid solution has a wide miscibility gap, for x between 0.05 and 0.89 at room temperature. However, there is a problem with derivation of consistent model for thermodynamic parameters of mixing for this system, including mixing enthalpy or Margules interaction parameters (Phan et al., 2019). It is thus of great interest to characterize W parameter for Li x FePO 4 system using the computational methods that were proven by us for lanthanide phosphates Ji et al., 2019a). The first systematic ab initio calculations of W parameters for monazite-type, orthophosphate solid solutions were performed by Li et al. (2014). In follow-up studies , we explained the quadratic dependence of W parameter on the difference in molar volumes of endmembers ( V) by a strain energy-based model, in which where E is the Young's modulus and V is the volume. These studies show that V is an important parameter that determines the value of W parameter. Considering the elastic properties of both endmember phases, this equation can be also written as : where G A is the shear modulus of the doped phase and B B is the bulk modulus of the dopant phase. The Young's, shear and bulk moduli of FePO 4 and LiFePO 4 phases as computed by Maxisch and Ceder (2006) are given in Table 5. The values of Margules interaction parameter computed with the Equations (4) and (5) are also reported. These values are too low and would result in thermodynamically stable solid solution at all Li content (such forms at room temperature for W < 5 kJ/mol, Li et al., 2014), which is inconsistent with experimental data of Yamada et al. (2006).
they obtained a good match to the measured values. However, in our estimate we used the experimental volumes ( Table 2) and any such a correction would require significant error in one of the reported volumes. In order to obtain the Margules interaction parameters that are consistent with the measured solubilities (Yamada et al., 2006, see below), instead of the reported V = 20 Å 3 between the two endmember phases it 1,2 , and W s 1,2 are the Margules interaction parameters estimated from Equations (4), (5), and solubility data of Yamada et al. (2006) (with Equations 7 and 8), respectively. should be V = 26 Å 3 . Our computed volumes, however, give even smaller difference of V = 14 Å 3 ( Table 2).
The experimental maximum solubilities of Yamada et al. (2006) show slightly asymmetric solid solution with maximum solubility of Li in FePO 4 of x 1 = 0.05 and content/depletion of Li in LiFePO 4 of x 2 = 0.89. Such an offset between both solubilities is consistent with the prediction of Equation (5). With these measured solubilities, the Margules interaction parameters can be derived. The maximum solubilities x 1 and x 2 at the two ends of slightly asymmetric solid solution or the relevant Margules interaction parameters W 1 and W 2 can be derived by solving self-consistently two equations (Mogilevsky, 2007;Ji et al., 2019a): These equations realize chemical equilibrium between Li (equilibrium concentration of x 1 in FePO 4 ) and Li-vacancy (equilibrium concentration of x 2 in LiFePO 4 ) in FePO 4 (first equation) and Li x FePO 4 (second equation) phases, respectively, and are derived by equality of the respected chemical potentials (Mogilevsky, 2007). The self-consistent solution of Equations (7) and (8) can be easily done in a numerical way. The resulted Margules interaction parameters for Li x FePO 4 solid solution obtained taking the measured solubilities as an input (Yamada et al., 2006) and assuming room temperature are given in  Li et al., 2014], which is in perfect agreement with the experimental measurements (∼2.5±1.0 kJ/mol, Stevens et al., 2006;Phan et al., 2019). Interestingly, the same values can be derived from Equation (5) and elastic moduli reported in Table 5, assuming V = 26 Å 3 . We note that such an offset is plausible because we consider a solid solution between Li cation and a vacant site, while the considered models were designed strictly for mixing of two cation species. We also attempted direct computation of Margules interaction parameter by ab initio methods. We obtained values between 13 and 28 kJ/mol, depending on the computational approach. However, we notice that such a small value of energy (enthalpy) is very sensitive to the computational setup, and can not be derived here precisely because of the inability of computing forces with the Wannier projectors scheme (see section 3.2). Other studies also result in large spread of predicted values (Figure 2 in Phan et al., 2019) and indicate sensitivity of this parameter to other effects, including the electronic entropy of Fe in different redox state in the Li x FePO 4 solid solution (Zhou et al., 2006). In Figure 3, we plot the free energy of mixing at room temperature, using the Margules interaction parameters extracted form the maximum solubilities measurements of Yamada et al. (2006). With the horizontal line we indicate the widths of wide miscibility and spinodal gaps. This information is important for understanding of charging relationships in the Li x FePO 4 system (Phan et al., 2019).  (Li et al., 2014), assuming configurational entropy of ideal mixing and neglecting other entropy contributions. R is the gas constant. W 1 and W 2 Margules interaction parameters are these from last column of Table 5. The miscibility and spinodal gaps are indicated by minima and inflection points marked with blue circles.

Ionic Conductivity
Because of high ionic conductivity, yttria-stabilized zirconia (see Figure 4 for the structure) is commonly used as solid state electrolyte. Its ionic conduction properties have been investigated in many studies, including experimental Carlson, 1964, 1965;Ioffe et al., 1978;Li et al., 1994;Lee et al., 2001Lee et al., , 2003Kilo et al., 2003;Zhang et al., 2007) as well as theoretical and atomistic modeling approaches (Kilo et al., 2003;Krishnamurthy et al., 2005b;Sizov et al., 2014;Bukaemskiy et al., 2021). The conductivity is usually described using Arrhenius-type equation (Ioffe et al., 1978;Ahamer et al., 2017) and the ionic conductivity can be written in the following simple form: where σ 0 is the pre-exponential factor and E a is the activation energy. In Table 6, we report the available experimental and computed data on E a together with the results of our simulations. We note that molecular dynamics simulations could be also used for simulation of ionic diffusion, but for systems with activation barriers close to ∼1 eV such a method requires long simulation times (e.g., 10 ns as applied by Sizov et al., 2014). Such simulation times are beyond the capability of ab initio molecular dynamics methods (capable of simulations at ps time scales) and could be performed only with less accurate description of interatomic interactions by simple interatomic potentials, like in studies of Sizov et al. (2014).
The computed activation energy of 0.92 eV is well-consistent with the measured values as well as with some previous theoretical predictions. Here, following the studies of Ahamer et al. (2017) and Guan et al. (2020) we assume diffusion along Y-Zr edge as a diffusion rate determining step. Interestingly, the activation barrier for transition along Zr-Zr edge is twice smaller as the one for transition along Y-Zr edge, which is in line with previous findings (Krishnamurthy et al., 2005b). The diffusion paths along these two edges are depicted in Figure 4. In the performed simulations with Eu as dopant, the activation barrier is comparable (0.97 eV). This is consistent with previous studies showing similar effects of other trivalent cations on stabilization of cubic zirconia and its ionic FIGURE 5 | The computed (lines) ionic conductivity in YO 1.5 − ZrO 2 system. The data (black filled circles) come from Strickler and Carlson (1964) and Liu et al. (2016). Frontiers in Energy Research | www.frontiersin.org conductivity (Krishnamurthy et al., 2005a). With the derived activation barrier we also computed the attempt frequency, ν 0 , using the transition state theory (Equation 2). We obtained value of 2.37 · 10 12 s −1 . With the computed activation energy and the attempt frequency, using Equation (1) we derived temperature and Y-content dependent ionic conductivity. The results are plotted in Figure 5 against the experimental data. Our model matches the measured values reasonably well, including temperature and Y content dependence. This is in great part due to inclusion of vacancy migration probability as estimated by Bukaemskiy et al. (2021) (see Equation 1). This shows that atomistic modeling and appropriate theoretical consideration can deliver accurate prediction for ionic conductivity in solid state electrolyte candidate materials. Nevertheless, as discussed by Ahamer et al. (2017) and Guan et al. (2020), the oxygen diffusion in yttria-stabilized zirconia is a complex process, which details, however, could be revealed by combination of computed and measured data.

CONCLUSIONS
In this contribution we presented an overview of our decadelong atomistic modeling contribution to the research on orthophosphates and zirconates type ceramics. We discussed the atomistic modeling derivation of structural, thermodynamic, and diffusion properties that are of importance for application of these materials as compounds in energy storage devices. In particular, we discussed the importance of the application of parameter free DFT+U approach and selection of realistic projector functions for counting the occupations of strongly correlated d and f orbitals for the Hubbard model-based DFT+U scheme. Only with this approach we were able to correctly reproduce the measured formation enthalpies of FePO 4 and LiFePO 4 phases. The consideration of the thermodynamic properties of Li x FePO 4 solid solutions indicates a system with a wide miscibility gap, which values and slight asymmetry are qualitatively consistent with the existing experimental data. Based on the measured solubility data we derived set of Margules interaction parameters that describe this solid solution. The resulting excess free energy of mixing shows wide miscibility and spinodal gaps at room temperature. We also discussed our studies of zirconium-based ceramics. In particular, we derived the stability diagram of yttrium-doped zirconia, showing that it stabilizes in cubic phase at Y content of ∼8%, well in line with the experimental measurements. The computed formation enthalpies along YO 1.5 − ZrO 2 solid solution are also well consistent with the measured data. With the computation of transition state and application of vacancy distribution model we were able to derive activation energies and temperaturedependent ionic conductivities for oxygen diffusion in this material, that are well-consistent with the measured data. This shows the power of carefully set up atomistic modeling for computation of various properties of ceramic materials as compounds for energy storage devices.
We discussed various successfully studies of application of atomistic modeling to prediction of a set of physical and chemical properties of orthophosphates and zirconates. In most cases, the best results have been obtained by a joint computational and experimental approach, or at least by extensive testing and comparison to the available experimental data. Application of a reliable, state-of-the-art ab initio approach is also a crucial factor contributing to this success. With the steady increase in the availability of computational power we expect that atomistic modeling research will be applied to tackle even more complex problems, such as kinetically driven dissolution or corrosion processes, and for effective screening of materials for energy applications.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
PK, ZH, and OC contributed equally to computing the data. PK performed collective analysis of the data and editing of the manuscript. All authors contributed to the article and approved the submitted version.