Thermodynamic properties of LiNiO₂, LiCoO₂, and LiMnO₂ using density-functional theory

Abstract

The formation energies of LiCoO₂, LiNiO₂ and LiMnO₂ were calculated using a combination of adequately selected Hess cycles and DFT computations. Several exchange–correlation functionals were tested and PBE for solids (PBEsol) turned out to be the most accurate. The enthalpies of formation at 0 K are −168.0 kJ mol at⁻¹ for LiCoO₂, −173.2 kJ mol at⁻¹ for LiNiO₂, −209.9 kJ mol at⁻¹ for o-LiMnO₂ and −208.8 kJ mol at⁻¹ for r-LiMnO₂. In comparison to experimental formation energy data, a difference of 1.6 and 0.01 kJ mol at⁻¹ was obtained for LiCoO₂ and LiMnO₂, respectively. By contrast, a much larger discrepancy, around 24 kJ mol at⁻¹, was obtained for LiNiO₂ and confirmed by using an additional and independent Hess cycle. The influence of slight crystallographic distortions associated with magnetism and/or the Jahn–Teller effect on energy was carefully searched for and taken into account, as well as corrections arising from vibrational contributions. Hence, these results should motivate future measurements of the thermodynamic properties of LiNiO₂, which are currently scarce. Vibrational contributions to the structural and energetic properties were computed within the harmonic and the quasi-harmonic approximations. The LiCoO₂ heat capacity at constant pressure is in excellent agreement with experimental data, with a difference of only 3.3% at 300 K. In the case of LiNiO₂ the difference reaches 17% at 300 K, which could also motivate further investigation. The C_p(T) value for the orthorhombic phase o-LiMnO₂, for which no previous data were available, was computed. Structural properties such as specific mass, bulk modulus and coefficient of thermal expansion are presented.

---

Introduction

Most of the lithium-ion batteries (LIBs) found on the market use NMC (which stands for nickel, manganese and cobalt) as the cathode active material. The NMC crystals have the general formula LiNi_xMn_yCo_zO₂ (x + y + z = 1), with compositions that vary from LiNi_0.33Mn_0.33Co_0.33O₂ to LiNi_0.8Mn_0.1Co_0.1O₂,¹ systematically increasing in Ni content. The growing need for the storage of electrical energy in the world is a strong motivation to optimize NMC-containing batteries. The thermal stability of the cathode material is related to its enthalpy of formation and the voltage profile during delithiation is related to the chemical potential of lithium.² Hence, good knowledge of the thermodynamic properties of the NMC system is essential for the development of these materials.

The CALPHAD (CALculation PHAase Diagrams) method^3,4 is a proven approach for predicting the thermodynamic, kinetic, and other properties of multicomponent materials systems. It couples phase diagram data and thermodynamic properties. The key thermodynamic property is the Gibbs energy from which all other thermodynamic properties can be derived. For each phase, a model is selected and the Gibbs energy is expressed with a polynomial function versus temperature and composition. The parameters of the models are fitted from selected experimental data such as the melting point, phase transition temperature, enthalpy of formation, enthalpy of mixing, heat capacity, and others. The functions and parameters obtained are stored in a database. Once the database is finalized, it provides the thermodynamic properties of the system for any composition. A CALPHAD model is usually refined as new data pertaining to the system are published, which may take place over periods of several years. The CALPHAD approach to the NMC system may be expected to provide a continuous description of properties related to the material according to composition and during the battery cycle, as demonstrated in the studies by Chang et al. in their CALPHAD descriptions of the pseudo-binary equilibria LiCoO₂–CoO₂⁵ and LiNiO₂–NiO₂.⁶

The NMC crystals have a rhombohedral structure that belongs to the R3̄m space group.⁷ We can represent them as a distorted rocksalt structure made of an FCC network of oxygen anions and another FCC network occupied by cations, in which Li planes alternate with planes containing transition metal cations along the [111] direction. The CALPHAD modelling of the LiNi_xMn_yCo_zO₂ solid solution includes the thermodynamic quantities of its three end-members which by definition are the compounds obtained assuming x = 1, y = 1 or z = 1, i.e. LiCoO₂, LiNiO₂ and LiMnO₂. In this work, we will refer to these compounds generically as LiMO₂, given M = Co, Ni or Mn. At room temperature, while LiCoO₂ and LiNiO₂ are stable and share the same rhombohedral structure as NMC, LiMnO₂ is only stable in its orthorhombic structure in the space group Pmmn.⁸ Additionally, LiNiO₂ is particularly complex experimentally, as it is almost impossible to synthesize the perfect stoichiometric compound.⁹ The thermodynamic properties of these end-members are necessary input data for the CALPHAD model which will eventually allow the optimized design of new NMC compositions.

Several measurements of the enthalpies of formation and heat capacities for the end-member LiCoO₂ have been published.^10–16 By contrast the data regarding LiNiO₂ are very scarce and limited to the enthalpy of formation of LiNiO₂¹⁷ and heat capacity at low temperatures.¹⁸ The end-member LiMnO₂ is not stable as a crystal belonging to the space group R3̄m¹⁹ and, therefore, there are no experimental measurements of its thermodynamic properties.

In the last few years, the boost in computational power allowed the ab initio calculation of physical properties of materials that agree with experimental data and permit the prediction of properties hitherto not known, which provides an invaluable source of input data for the construction of CALPHAD databases. Density-functional theory (DFT)^20,21 has been used by many authors to calculate the physical properties of lithium-ion battery cathode materials.^22–30 The heat capacity of LiCoO₂ at constant volume (C_v) was obtained by means of DFT by Wu et al.²² In a similar work, Du et al.²³ calculated the heat capacity of LiMO₂ as a function of temperature using density-functional perturbation theory (DFPT)³¹ and a harmonic approximation (HA), although it is mentioned that a quasi-harmonic approximation (QHA) was employed. This approach may be acceptable at low temperatures, where the calculated C_v agrees well with the experimental C_p values of LiCoO₂.

Regarding the study of formation energies, Chang et al.²⁴ calculated the formation energies of LiCoO₂ and LiNiO₂ from pure elements with a combination of DFT and an empirical method, which consists in adding or subtracting a certain amount of energy of a M_yO_x species, depending on the oxidation state of the transition metal. Xu et al.²⁵ used the Generalized Gradient Approximation (GGA) to compute the formation energies of LiCoO₂, LiMnO₂, LiMn₂O₄ and LiFePO₄ using the PBE description of GGA. Kramer and Ceder²⁶ studied the morphology of LiCoO₂ by means of DFT+U. Longo et al.²⁷ focused on the study of layered and over-lithiated Mn oxides, Tuccillo et al.²⁸ studied the stability of the compounds LiMO₂ in three different structures and Tsebesebe et al.²⁹ assessed the thermodynamic, electronic, and mechanical properties of the same family of compounds. Finally, Lee et al.³⁰ computed the formation energy of LiCoO₂ with the support of molecular dynamics. A compilation of the results obtained by these authors is given in Table 1.

Table 1 Formation energies from elements obtained in previous works. The relative error is calculated based on standard enthalpy of formation experimental data^{10,15,17,32,33}

Ref.	Calculation method	Compound	Space group	Calculated formation energy (kJ mol at^−1)	Absolute error (kJ mol at^−1)	Relative error (%)	Transition metal 3d electrons U value (eV)
24	DFT + empirical method	LiNiO2	R3̄m	−150.5	2.2	1.5	—
		LiCoO2	R3̄m	−173.4	3.7	2.2
25	GGA	LiCoO2	R3̄m	−91.5	78.1	46.1	—
		LiMnO2	C2/m	−178.9	31.0	14.7
26	GGA+U	LiCoO2	R3̄m	−171.8	2.1	1.3	3.3
27	GGA+U	LiMnO2	R3̄m	−192.6	17.3	8.2	5.2
28	GGA+U	LiCoO2	R3̄m	−183.5	13.9	8.2	4
		LiNiO2	C2/m	−154.8	6.4	4.4	4
		LiMnO2	Pmmn	−209.0	0.9	0.4	4
29	GGA+U	LiCoO2	R3̄m	−41.7	127.9	71.9	6.1
		LiNiO2	R3̄m	−43.1	105.3	71.0	2.45
		LiMnO2	R3̄m	−32.0	178.0	78.5	4.25
30	Molecular dynamics	LiCoO2	R3̄m	−165.4	4.2	2.5	—

---

Table 1 displays a great variability of the results. A likely cause of this variability is the diversity of formation reactions used in the calculations. In some works, the reactants are the pure elements Li and M in their crystal forms and molecular O₂ while, in other works, they are combinations of various oxides where a redox reaction takes place. The involvement of molecular O₂ in the reaction requires a correction since the GGA underestimates the O₂ binding energy by about 150 kJ mol⁻¹.³⁴ Additionally, the choice of Hubbard-U values is not straightforward³⁵ and this is obviously yet another source of variability. The U value adds an energy correction to selected electron states in order to obtain a more localized electron wavefunction and is commonly used in calculating 3d transition metals oxide energies.^36–39 Some studies, however, find better agreement with experiments without using U values.⁴⁰ Additionally, the Hubbard U approach is especially problematic for electrochemical processes.³⁹

The main objective of this study is to compute the thermodynamic properties of the three end-members as a function of temperature using DFT, which is part of an effort aiming at building a CALPHAD model of the NMC system.

In order to estimate and minimize the uncertainties in our results, we studied in detail several aspects of the methodology: (1) the various possibilities for writing the DFT-computed chemical reaction that gives access to the formation energy; (2) the most appropriate theory level for computing the vibrational contributions; (3) the crystallographic structure used as input data to the DFT computations and (4) The specific influence of the exchange–correlation (x–c) functional on the quality of the results, which has been discussed above, in view of the data in Table 1. In this paper, we report on the calculated thermodynamic properties and also on the methodology, the latter being an essential element on which the CALPHAD input data assessments are based.

We now briefly introduce the first three methodological aspects. The choice of the DFT-computed reaction is a central feature. Any chemical reaction involving the LiMO₂ compound of interest and other compounds for which reference enthalpy data are available can be used for this purpose. Extracting the desired standard formation energy of the LiMO₂ compound from the DFT results is then performed classically using Hess law.⁴¹ Although the availability of the required enthalpy reference data somewhat constrains the choice of the DFT-computed reaction, many reactions can be considered. In the following, the possible options are sorted according to the fact that this DFT-computed reaction involves oxidation and reduction, i.e. a redox reaction, and more precisely according to the change in the oxidation state of the transition metal. Since the scheme with no-redox reactions comes out as the best, we shall present it first.

As for the vibrational contributions, we compute the corrections to the energetic properties, namely zero-point energy and enthalpy increments with temperature, both within the harmonic approximation (HA) and the quasi-harmonic approximation (QHA). At no additional computational cost, this gives us access first to the heat capacities at constant volume (C_v) and constant pressure (C_p), within the QHA, and second to structural properties such as specific mass, bulk modulus and thermal expansion coefficient.

The choice for the crystallographic description of the DFT input structure sometimes has a significant influence on the energy results so that it is also a major aspect of our methodology. Ideally, the only input data for a DFT study should be the chemical composition and the total number of electrons. However, most studies of crystalline compounds use the experimentally determined structure as input. This makes the DFT result sensitive to possible errors or uncertainties in the experimental data and, consequently, less “ab initio”. In our work, we depart from this standard procedure to some extent by retaining a structure that leads to a lower energy, in case we find one, while remaining compatible with the available structural data.

In the next section, we present the first step of our methodology, which consists of the selection of suitable chemical no-redox reactions according to available reference data, the design of Hess cycles and the choice of options regarding the DFT levels. The following section describes the computational details. The calculation results are reported and discussed subsequently, where we compare and assess the various aspects of the methodology. Finally, we draw conclusions.

Methodology

We calculate the formation energy of each end-member using a Hess cycle that mixes a DFT-computed chemical reaction with several auxiliary reactions for which reference reaction energy data are available in the literature. We classify the DFT-computed reactions into three categories. Type 1 reactions involve only oxide compounds and no oxidation or reduction of the various chemical elements, referred to as “no-redox reactions” in this work. Type 2 reactions are all-solid-state reactions involving charge transfer. Type 3 reactions are standard formation reactions involving all elements in their standard state at 300 K (gaseous molecular oxygen, metallic Li and transition elements). Type 2 and 3 reactions involve electron transfer and are referred to as “redox reactions”. This section is restricted to Type 1 reactions, which are the cornerstone of our methodology, while Type 2 and Type 3 reactions will only be used in the following sections for comparison purposes. The energy of each compound taking part in the reaction, reactant or product, is computed using DFT to derive the reaction energy.

Table 2 shows the DFT-computed formation reactions evaluated for each end member.

Table 2 Type 1 DFT-computed reactions used to calculate the formation energy of NMC end-members

End-member	Reactions assessed	Reaction #
LiCoO2	Li2O + Co3O4 → 2LiCoO2 + CoO	1
LiMnO2	Li2O + Mn2O3 → 2LiMnO2	2
LiNiO2	Li2O + 2NaNiO2 → 2LiNiO2 + Na2O	3
	Li2O + 2NiTiO3 → 2LiNiO2 + Ti2O3	4

---

Experimental energies of formation are usually expressed using the pure elements in their standard state. Thus, we choose an appropriate Hess cycle in order to compare experimental and calculated values. Table 3 lists all Hess cycles used for each reaction and all reactions making up each cycle, as well as the associated reaction energies.

Table 3 Thermodynamic cycles employed to obtain formation energies of LiMO₂ (M = Co, Ni, Mn) from pure elements in their standard state. Each cycle includes a Type 1 DFT-computed formation reaction (i.e. without any redox process)

Compound	Reaction #	Hess’ cycle	Energy
LiCoO2	1	Li2O + Co3O4 → 2LiCoO2 + CoO	ΔER1, calculated with DFT
		2Li + 0.5O2 → Li2O	ΔELi2O, experimental
		3Co + 2O2 → Co3O4	ΔECo3O4, experimental
		Co + 0.5O2 → CoO	ΔECoO, experimental
		Li + Co + O2 → LiCoO2	ΔfE,R1,LiCoO2 = 0.5(ΔER1 + ΔELi2O + ΔECo3O4 − ΔECoO)
LiMnO2	2	Li2O + Mn2O3 → 2LiMnO2	ΔER2, calculated with DFT
		2Li + 0.5O2 → Li2O	ΔELi2O, experimental
		2Mn + 1.5O2 → Mn2O3	ΔEMn2O3, experimental
		Li + Mn + O2 → LiMnO2	ΔfE,R2,LiMnO2 = 0.5(ΔER2 + ΔELi2O + ΔEMn2O3)
LiNiO2	3	Li2O + 2NaNiO2 → 2LiNiO2 + Na2O	ΔER3, calculated with DFT
		2Li + 0.5O2 → Li2O	ΔELi2O, experimental
		Na + Ni + O2 → NaNiO2	ΔENaNiO2, experimental
		2Na + 0.5O2 → Na2O	ΔENa2O, experimental
		Li + Ni + O2 → LiNiO2	ΔfE,R3,LiNiO2 = 0.5(ΔER3 + ΔELi2O + 2ΔENaNiO2 – ΔENa2O)
	4	Li2O + 2NiTiO3 → 2LiNiO2 + Ti2O3	ΔER4, calculated with DFT
		2Li + 0.5O2 → Li2O	ΔELi2O, experimental
		Ni + Ti + 1.5O2 → NiTiO3	ΔENiTiO3, experimental
		2Ti + 1.5O2 → Ti2O3	ΔETi2O3, experimental
		Li + Ni + O2 → LiNiO2	ΔfE,R4,LiNiO2 = 0.5(ΔER4 + ΔELi2O + 2ΔENiTiO3 − ΔETi2O3)

---

We computed the energies of the reactions listed in Table 2 using various choices of DFT parameters. We shall first report results obtained with our reference parameter set which consists of the GGA-PBEsol⁴² x–c function at 0 K and does not include vibrational contributions, a condition denoted “static”. From then on, in an effort to minimize the uncertainty, we shall evaluate various modeling options by comparison with this reference set of results. As the reaction energies of the auxiliary reactions listed in Table 3 are actually experimental enthalpies of formation at 300 K, in a second step, we shall report the DFT computed reaction energies also at 300 K, which were obtained by complementing static data with phonon calculations within the harmonic (HA) or quasi-harmonic (QHA) approximations. The fundamental mathematical formulae used for obtaining the thermodynamic quantities calculated in this work are listed in the ESI.†

Due to the significant computational cost of the latter correction, however, we compute it only for LiCoO₂ and show that the error introduced by ignoring this correction is small. In a later section, we shall also compare these reference results with static formation energies obtained with other functionals, namely PBE,⁴³ PW91⁴⁴ and rSCAN,⁴⁵ as well as with reactions involving a redox process. Then, we also compare them with GGA+U calculations using different U values in ranges that have been recommended in previous works. For Co-3d, we set the Hubbard-U value as 2.6, 3.3 and 4.0 eV. For Ni-3d, we chose 4.0, 6.4 and 8.0 eV. The DFT+U calculations involving species with Mn could not be performed because the calculations did not converge using the set of parameters listed in the Computational details section. The second value of each set of Hubbard-U values corresponds to a value found in the studies by Wang et al.,³⁶ where the U values were fitted according to the energy of the oxidation reactions 6 CoO + O₂ → 2 Co₃O₄ and 2 NiO + O₂ → 2 NiO₂. We did not have enough reliable data to reasonably assign different U values to different oxidation states of a given chemical element, leading to inaccuracies when computing the energy of Co₃O₄, for instance, a mixed-valence compound that contains both Co⁺² and Co⁺³. Table 1 illustrates the variety of U values among different studies that sometimes use an average of values found in the literature, and sometimes set a value that best fits a specific property. Thus, we arbitrarily choose a larger and a lower value for each transition metal in order to quantify its influence on the results.

We make use of the fact that for all considered solid-state compounds at a pressure P of 1 atm, the difference between enthalpy and energy is negligible. Thus for data calculated statically at P = 0 or with the QHA at P = 1 atm, we use the terms “formation energies” or “enthalpies of formation” interchangeably. We stick to the term “energy” when computed with the HA at 300 K, since the system is then under much higher pressures or in situations where the pressure is not specified. For the gaseous O₂ phase, we take into account the difference between molar energy and enthalpy by assuming a constant value C_p − C_v = R.

As for the vibrational contributions, the heat capacity of LiMO₂ at constant volume (C_v) is obtained using the HA and the heat capacity at constant pressure (C_p, P = 1 atm) is obtained within the QHA.⁴⁶ We calculated the C_p curves in a temperature range from 0 to 1200 K, which corresponds to static pressures varying from −2 to −9.5 GPa, depending on the LiMO₂ compound. At each temperature, the volume dependence of the Helmholtz free energy (F) was obtained with phonon calculations over at least four different volumes and it was least-squares fitted to a parabolic equation, as detailed in the ESI.† All phonon calculations are performed with the PBEsol functional. In cases where the full QHA model was unusable due to the presence of imaginary modes of vibration at large negative static pressures, we replaced it with the simplified QHA model of Fleche.⁴⁷ This model is based on DFT data consisting of the static pressure–volume relationship, which we fitted using the Birch–Murnaghan equation of state, and the phonon spectrum computed at zero static pressure. The Fleche model then uses the Debye model to process the acoustic phonons and the Mie–Grüneisen theory to approximate the optical Grüneisen coefficient. In cases where both can be used, the full QHA model and the Fleche's model only show small differences. We noted that for all the compounds studied in this work, we always achieved a complete absence of unstable vibration modes at zero static pressure. To a large extent, this is due to the fact that we have always searched for the structure with the lowest static energy at zero static pressure, starting from the experimentally known structure and allowing for small distortions by considering a variety of supercells. Neumann–Kopp approximations from binary oxides are also presented along with our QHA results. With the QHA, we also obtain the cell specific mass, the bulk modulus and the thermal expansion coefficient of the LiMO₂ compounds.

Crystal structures

All the structures were built based on low-temperature experimental information, which includes cell parameters, Wyckoff positions and magnetic structure. Magnetism is taken into account in DFT by spin-polarization. When the antiferromagnetic structure is the most stable, the set of symmetry operations must be changed accordingly compared to the symmetry of the mass distribution. In such cases, or wherever our calculated structure is a slightly distorted supercell of the experimentally known structure, we present the crystallographic details including the spin direction for each of the transition metal orbits in the ESI.† This enables a complete description of these structures without ambiguity. The procedure to obtain such distorted structures is a simplified version of the procedure of Togo and Tanaka.⁴⁸ We build a set of supercells with the original structure taken from the literature for the same compound. In each supercell, we introduce a slight atomic displacement or cell strain and we submit it to a geometry optimization task with CASTEP. The result we reported is the most stable structure according to CASTEP. Compared to the original structure from the literature, it usually contains small distortions of typically 0.05 Å. The geometry-optimized structure is symmetrized using a symmetry finder software tool with a tight tolerance such as 0.005 Å. The Wyckoff positions are then easily read either from the software tool output itself or from e.g. the International Tables of Crystallography Vol A or the Bilbao Crystallographic Server.⁴⁹ We stop the screening as soon as a significant energy lowering is achieved. We do not claim that our reported structure is the minimum energy state because there is always the possibility that some other supercell, not included in our screening procedure, would have a lower energy.

Table 4 summarizes the experimental structures and enthalpies of formation of all the compounds involved, as well as the structures used for the DFT calculations. The list contains additional species not mentioned in Table 3 that will be employed in a later section when evaluating reactions involving a redox process.

Table 4 List of compounds employed in this study. For each compound, the crystal lattice structure (prototype, Pearson's symbol and space group number), experimental standard formation energy and magnetic structure are given. NSP, FM and AFM stand for “non-spin-polarized”, “ferromagnetic” and “antiferromagnetic”, respectively. LiNiO₂, CoO, Mn₂O₃ and Mn are special cases discussed in the main text

Compound	Crystal lattice structure	Enthalpy of formation at 300 K (kJ mol at^−1)	Input structure for the DFT calculation (spin polarization, space group number) and reference for magnetic structure
LiCoO2	NaFeO2, hR12, 166	−169.9 ± 0.6^15	NSP, 166^50
		−168.3 ± 0.6^10
		−170.6 ± 0.5^16
LiMnO2	NaMnO2, oP8, 59	−209.9 ± 0.5^33	AFM, 13^51
LiNiO2	NaFeO2, hR12, 166	−148.3 ± 0.6^17	FM, 14
Co3O4	Fe3O4, cF56, 227	−130.1 ± 0.2^52	AFM, 216^53
CoO	NaCl, cF8, 225	−119.4 ± 0.6^54	AFM, 122
MnO	NaCl, cF8, 225	−191.3^55	AFM, 166^56
Mn2O3	Mn2O3, oP80, 61	−192.3 ± 0.2^57	AFM, 19
NiO	NaCl, cF8, 225	−119.7 ± 0.2^54	AFM, 166^56
NaNiO2	NaNiO2, mS8, 12	−157.3^58	AFM, 12^59
NiTiO3	TiFeO3, hR30, 148	−240.4^60	AFM, 146^61
Ti2O3	Al2O3, hR30, 167	−303.4^62	FM, 167^63
Na2O	CaF2, cF12, 225	−138.3^58	NSP, 225
Li2O	CaF2, cF12, 225	−199.1 ± 0.1^64	NSP, 225
Li2O2	Li2[O2]-a, hP8, 194	−161.1^65	NSP, 194
Li	W, cI2, 229	0	NSP, 229
Co	Mg, hP2, 194	0	FM, 194
Mn	Mn, cI58, 217	0	AFM, 215^66
Ni	Cu, cF4, 225	0	FM, 225
O2		0	FM

---

Among the species considered in this work, the CoO crystal shows a particular inconsistency between experimental and calculated data. The DFT calculations show that CoO in the sphalerite structure with space group I4̄2d (Table S4, ESI†) is about 5 kJ mol at⁻¹ more stable than its known rock-salt configuration. This phase is reported in the work of Redman and Steward,⁶⁷ in which the authors consider it metastable. Sui et al.⁶⁸ observed that the sphalerite phase is predominant in the early stages of CoO crystal growth and does not disappear in the later stages. It must be noted that the XRD patterns of the rock salt and sphalerite phases are fairly similar, so data with good signal-to-noise ratio and high resolution may be required to differentiate them. This arises the question of which structure was actually associated with the formation enthalpy measurements. The experimental enthalpy of formation of CoO used in this work was obtained in the study by Boyle et al.⁵⁴ where the X-ray diffraction spectrum data contained additional lines which did not belong to the rock-salt CoO structure. Ghosh et al.⁶⁹ studied in detail the DFT description of the energy ordering of these phases and their dependence on the x–c functional but did not include PBEsol in their study. The known antiferromagnetic (AFM) structure of the rock-salt CoO is also questionable. We find an AFM arrangement in the P2₁/m space group (Table S5, ESI†) that is 0.3 kJ mol at⁻¹ more stable than the one in the C2/m configuration proposed by Jauch et al.⁷⁰

The compound LiMnO₂ is stable in an orthorhombic configuration (o-LiMnO₂, space group Pmmn).^71,72 At low temperatures, o-LiMnO₂ is antiferromagnetic and a magnetic structure was proposed by Kellerman et al.⁵¹ This structure can belong to either P2/c or C2/c space group. We performed DFT calculations for both structures and found similar energies for both. Atom coordinates and cell parameters are given in Table 5 and Table 6. We proceeded with the P2/c structure in our calculations.

Table 5 Cell parameters of the AFM structure of o-LiMnO₂ with two different configurations of identical energy. Cell dimensions a, b and c in Å

Cell parameter	P2/c	C2/c
a	5.341	9.132
b	5.686	11.367
c	5.341	5.34
α and γ	90°	90°
β	117.515°	148.745°

---

Table 6 Atom coordinates of the AFM structure of o-LiMnO₂ with two different configurations of identical energy

Structure	Elements	Wyckoff position	x	y	z
P2/c	Li	2e	0	0.120	0.25
	Li	2f	0.5	0.880	0.25
	Mn	2e	0	0.634	0.25
	Mn	2f	0.5	0.366	0.25
	O	4g	0.25	0.398	0
	O	4g	0.25	0.142	0
C2/c	Li	4e	0	0.310	0.25
	Li	4e	0	0.810	0.25
	Mn	4e	0	0.067	0.25
	Mn	4e	0.5	0.062	0.25
	O	8f	0.25	0.179	0.75
	O	8f	0.25	0.449	0.75

---

The LiNiO₂ compound, as LiCoO₂, crystallizes in the α-NaFeO₂ rhombohedral system.^73,74 However, a local Jahn–Teller distortion is observed by EXAFS spectroscopy.⁷⁵ Neutron diffraction studies showed that the distorted structure is best fitted to the C2/m group at low temperatures.⁷⁶ Several ab initio studies were carried out to better understand the possible cooperative Jahn–Teller distortion at low temperatures^23,77–81 and all conclude that a distorted structure is more stable than the non-distorted rhombohedral structure. In this work, we evaluated some configurations of Jahn–Teller distortions of the rhombohedral LiNiO₂, and found that the structure belonging to the P2₁/m space group (Table S3, ESI†) to be the most stable, with a maximal deviation of 0.09 Å from the original non-distorted structure, in agreement with the results obtained by Chen et al.⁸⁰ and Sicolo et al..⁸¹

The stable configuration obtained for Mn₂O₃ belongs to the P2₁2₁2₁ space group (Table S1, ESI†), with a maximal deviation of 0.04 Å from its experimental Pbca structure.

The magnetic structure of Mn metal has not yet been completely solved, in spite of many experimental and theoretical efforts over the past 70 years (see e.g. Manago et al.⁸² and references therein). Since we are only interested in the energy, the finer details of the magnetic structure may be neglected. We followed the description given by Hobbs et al.⁶⁶ for collinear antiferromagnetism which can be computed in a space group P4̄3m. rSCAN is a special case as it leads to strong overmagnetization.⁸³ Since the magnetic moments in α-Mn are strongly dependent on the geometry of the crystal, our rSCAN energy result may be off by up to 3 kJ mol at⁻¹. As a matter of fact, rSCAN finds that a ferrimagnetic arrangement in I4̄3m more stable than the AFM arrangement in P4̄3m, by about 7 kJ mol at⁻¹.

The O₂ molecule was built in a cell with a, b and c parameters equal to 20, 21 and 22 Å, respectively and α, β and γ equal to 90° and a spin polarization corresponding to its triplet electronic ground state.

Computational details

Most DFT computations were performed using the CASTEP^84–87 code implementing the plane-wave pseudopotential method. The BIOVIA Materials Studio software⁸⁸ was employed to build all the structures and to analyze the results. Energies of formation were computed by performing geometry optimizations of all structures involved. The generalized gradient approximation (GGA) was used with and without the GGA+U correction for the strong correlations of d electrons. Several x–c functionals were used. Formation energies and phonons were performed with the PBEsol⁴² functional. The formation energy was also computed with PBE,⁴³ PW91⁴⁴ and rSCAN⁴⁵ functionals. The “2017R2” ultrasoft pseudopotentials library was used. For the computation of the static formation energies, the values of the numerical parameters were chosen after a convergence study so as to reach a precision of 0.1 kJ mol at⁻¹ on reaction energies. The energy cutoff was set to 630 eV and the self-consistent field (SCF) convergence threshold to 1 × 10⁻⁶ eV per atom. The k-point sampling was generated by Monkhorst and Pack's scheme⁸⁹ with a k-point interval of 0.04 Å⁻¹ on all axes. All calculations are spin-polarized unless stated otherwise. For formation energy calculations, the geometry of all crystals was fully relaxed up to a maximal force of 0.02 eV Å⁻¹, maximal stress of 0.02 GPa, maximal atom displacement of 0.001 Å and a relaxation energy convergence of 1 × 10⁻⁵ eV atom⁻¹. We assumed the experimental configuration as the starting point, with the exception of a few compounds found to have more stable configurations at 0 K, in which case their crystallography information is given in the ESI.†

Phonon calculations were performed using CASTEP with the finite displacement method on Li_xM_xO_2x (M = Co, Ni, Mn) supercells relaxed with an energy convergence of 5 × 10⁻⁶ eV atom⁻¹, a maximal force of 0.003 eV Å⁻¹, a maximal atom displacement of 5 × 10⁻⁴ Å and a maximal stress of 0.01 GPa. We use a SCF convergence of 1 × 10⁻¹¹ eV standard and fine grids of 1.75 and 2.1, respectively, and supercell sizes consistent with a cutoff of 4.5 Å for the force constants. The supercells are described in the ESI.† PBEsol was used to obtain both the Birch–Murnaghan equations of state and the lists of phonon frequencies. In cases where the phonon spectrum contained imaginary frequencies, we checked whether these modes actually reflected the instability of the structure, as opposed to numerical uncertainties associated with the calculation of force constants by recomputing the phonon spectra using the ABINIT code^90–94 with PAW pseudopotentials and linear response theory (DPFT).

Results

Static formation energies

The enthalpies of formation obtained with Type 1 no-redox reactions and the PBEsol functional for the three LiMO₂ compounds are outlined in Table 7:

Table 7 Enthalpies of formation in kJ mol at⁻¹ of LiMO₂ compounds computed with the static method. The two values assigned for LiNiO₂ refer to the two reactions reported in Table 2

End-member	Calculated enthalpy of formation	Experimental value	δ th-exp	Relative error (%)
LiCoO2	−167.97	−169.9 ± 0.6^15	1.93	1.1
		−168.3 ± 0.6^10	0.33	0.2
		−170.6 ± 0.5^16	2.63	1.5
o-LiMnO2	−209.86	−209.87 ± 0.5^33	0.01	0.003
r-LiMnO2	−208.82	—
LiNiO2	−173.05	−24.75	−24.75	16.7
	−173.32	−25.02	−25.02	16.9

---

It should be noted that the calculated enthalpy of formation of LiCoO₂ was obtained with the sphalerite crystal structure of CoO. The result obtained with the monoclinic structure reported by Jauch et al.⁷⁰ is −166.79 kJ mol at⁻¹, 1.18 kJ mol at⁻¹ above the value reported in Table 7.

We observe very good agreement with experimental data for LiCoO₂ and o-LiMnO₂. However, we find a striking discrepancy between the calculated and experimental values of LiNiO₂, which is an order of magnitude larger than what we obtained for the other end-members. The fact that the calculated values for LiNiO₂ agree within the two proposed reactions (#3 and #4 in Table 3) motivated some further investigation about this discrepancy with the experiment. To this end, in the first step, we searched for a possible issue with the DFT computation of Ni(III), as opposed to Co(III) and Mn(III). In the second step, we studied the influence of the x–c functional, and, in the third step, we recomputed the LiMO₂ formation energies using other formation reactions (i.e. Type 2 and Type 3 as defined above). All of these results are presented in the following subsections of this paper. First, we applied our method to compute the formation energy of NiTiO₃ according to the Hess’ cycle in Table 8, which involves the same compounds already employed in Table 2 except the Li-containing ones.

Table 8 Hess’ cycle used for obtaining the formation energy of NiTiO₃

Hess’ cycle	Energy	Value (kJ mol^−1)
Ti2O3 + 2NaNiO2 → 2NiTiO3 + Na2O	ΔER, calculated with DFT	−41.23
2Ti + 1.5O2 → Ti2O3	ΔETi2O3, experimental	−1517.05^62
Na + Ni + O2 → NaNiO2	ΔENaNiO2, experimental	−629.28^58
2Na + 0.5O2 → Na2O	ΔENa2O, experimental	−414.82^58
Ni + Ti + O2 → NiTiO3	ΔEf,NiTiO3 = 0.5(ΔER + ΔETi2O3 + 2ΔENaNiO2 − ΔENa2O)	−1201.03

---

The calculated NiTiO₃ enthalpy of formation is −240.21 kJ mol at⁻¹. A difference of 0.19 kJ mol at⁻¹ is found in comparison with the experimental value of −240.4 kJ mol at⁻¹ by Kubaschewski et al.⁶⁰ This excellent agreement supports our attribution of the origin of the discrepancy between the calculated and experimental values for LiNiO₂ to the Li containing compounds, and more precisely to LiNiO₂ itself since Li₂O was validated by the study of LiCoO₂ and LiMnO₂. This issue is analyzed further in the Discussion section.

To wrap up the static formation energies section, we now go back to LiMnO₂. As mentioned in the Methodology section, the stable LiMnO₂ adopts an orthorhombic structure (o-LiMnO₂) and this is used in our calculations. However, the LiMnO₂ compound of interest as an end-member of the NMC system has a rhombohedral configuration (r-LiMnO₂), which is metastable. Among the different rhombohedral structures that we tested, including FM and some AFM configurations, a distorted configuration belonging to the P2/c space group was the most stable. The calculated enthalpy of formation of this structure is about 1.0 kJ mol at⁻¹ above that of the stable orthorhombic o-LiMnO₂. The crystallographic description of the r-LiMnO₂ structure is included in the ESI† in Table S2.

Vibrational contribution to the formation energy

The corrections to the zero-point energy (ZPE) and vibrational energy at 300 K are often overlooked in computations of formation energies because they are usually marginal in comparison to the DFT precision with respect to the experimental data and they require calculations that are computationally heavy in contrast to static calculations. In the present study, however, where the discrepancy between static DFT and experiment was found to reach up to a few kJ mol at⁻¹, it becomes fully justified to study this correction. Thus, we compute this contribution to the formation energy of LiCoO₂.

The energy corrections obtained with the HA for all the compounds involved in the DFT-computed reaction (reaction 1 of Table 2) are given in Table 9. We choose energy references for the three chemical elements Li, Co and O in such a way that the static total energies of Li₂O, CoO and Co₃O₄ are all equal to zero. The fourth column (internal energy at 300 K) is the sum of the three other columns (static total energy + zero point energy + energy gain from 0 to 300 K).

Table 9 Corrections to internal energies of species involved in Reaction 1 of Table 2, and (last line) corrections to the reaction energy, arising from phonon contributions computed within the HA. Values in kJ mol⁻¹

Compound	Static total energy	Zero point energy	Energy gain from 0 to 300 K	Internal energy at 300 K
Li2O	0	23.28	7.55	30.83
CoO	0	10.15	5.11	15.26
Co3O4	0	51.71	19.32	71.03
LiCoO2	−36.77	29.71	10.56	3.49
Reaction 1	−73.55	−5.42	−0.64	−79.61

---

While the total HA contribution to the internal energy of each compound turns out to be of the order of several tens of kJ mol at⁻¹, the reaction energy of reaction 1 changes by only about 8% or, in absolute terms, −6.1 kJ mol⁻¹. This is in line with the common tendency of the heat of solid-state reactions to be little dependent on temperature. From the point of view of the differences reported in Table 7 however, this correction is significant. Thus, we shall take it into account below in the assessment of the precision of our calculation method.

We also computed the vibrational contributions to the LiCoO₂ energy within the QHA. At a pressure of 1 bar, the energy gain from 0 to 300 K is equal to 2.8 kJ mol at⁻¹, which is only 0.2 kJ mol at⁻¹ greater than that obtained with HA. Assuming that the difference between HA and QHA for Li₂O, CoO and Co₃O₄ is of the same order of magnitude and that some compensation between reactants and products also occurs, we may presume that the QHA does not bring any significant improvement for computing formation energies, compared to the HA.

As reported in the thermodynamic assessment of Chen et al.,⁹⁵ CoO undergoes an AFM to paramagnetic phase transition at a temperature close to 300 K. This transition affects its C_p as it shows an abrupt increase starting at around 200 K with a peak at 290 K. Integration of this experimental C_p gives an enthalpy increment from 0 to 300 K of 4.7 kJ mol at⁻¹. As a result of the difference with the value in Table 9, LiCoO₂ has a formation energy of 0.05 kJ mol at⁻¹ which is lower than that obtained without this correction. All these corrections to the formation energy of LiCoO₂ are reported in Table 10.

Table 10 Formation energy of LiCoO₂ in kJ mol at⁻¹ according to the static method and the corrections of ZPE and energy gain at 300 K. The energy difference and relative error are calculated in relation to the average of available experimental data, −169.6 ± 1.7 kJ mol at⁻¹.^10,15,16

Method	Formation energy	δ th-exp	Relative error (%)
Static	−167.97	1.63	1.0
Including ZPE + T = 300 K, HA	−168.73	0.87	0.5
Including ZPE + T = 300 K, HA + real CoO Cp	−168.77	0.82	0.5

---

Influence of x–c functional and Hubbard U on the formation energy

As mentioned in the Methodology section, we tested not only reactions with no redox process (Type 1), listed in Table 3, but also other reactions with oxides involving redox (Type 2) and the reactions with pure elements in their standard states at 300 K (Type 3). These additional reactions are listed in Table 11. Reactions 5, 7 and 9 belong to Type 2 whereas reactions 6, 8 and 10 are of Type 3 (with gaseous molecular oxygen, metallic Li and transition element).

Table 11 Thermodynamic cycles employed to obtain formation energies of LiMO₂ from reactions with redox processes. The experimental enthalpies of formation are listed inTable 4

Compound	Reaction #	Hess’ cycle and reactions	Energy	Transition metal electronic transfer
LiCoO2	5	Li2O2 + 2CoO → 2LiCoO2	ΔER5, calculated with DFT	2Co^+2 → 2Co^+3 + 2e^−
		2Li + O2 → Li2O2	ΔELi2O2, experimental
		Co + 0.5O2 → CoO	ΔECoO, experimental
		Li + Co + O2 → LiCoO2	ΔfE,R5,LiCoO2 = 0.5(ΔER5 + ΔELi2O2 + 2ΔECoO)
	6	Li + Co + O2 → LiCoO2	ΔfE,R6,LiMnO2, calculated with DFT	Co → Co^+3 + 3e^−
LiMnO2	7	Li2O2 + 2MnO → 2LiMnO2	ΔER7, calculated with DFT	2Mn^+2 → 2Mn^+3 + 2e^−
		2Li + O2 → Li2O2	ΔELi2O2, experimental
		Mn + 0.5O2 → MnO	ΔEMnO, experimental
		Li + Mn + O2 → LiMnO2	ΔfE,R7,LiMnO2 = 0.5(ΔER8 + ΔELi2O2 + 2ΔEMnO)
	8	Li + Mn + O2 → LiMnO2	ΔfE,R8,LiMnO2, calculated with DFT	Mn → Mn^+3 + 3e^−
LiNiO2	9	Li2O2 + 2NiO → 2LiNiO2	ΔER9, calculated with DFT	2Ni^+2 → 2Ni^+3 + 2e^−
		2Li + O2 → Li2O2	ΔELi2O2, experimental
		Ni + 0.5O2 → NiO	ΔENiO, experimental
		Li + Ni + O2 → LiNiO2	ΔfE,R9,LiNiO2 = 0.5(ΔER + ΔELi2O2 + 2ΔENiO)
	10	Li + Ni + O2 → LiNiO2	ΔfE,R10,LiNiO2 calculated with DFT	Ni → Ni^+3 + 3e^−

---

We also computed the same reaction energies with other x–c functionals, including PBE + U with different U values.

Fig. 1 shows all our results of computed enthalpies of formation with the various x–c functionals and several Hubbard U values. For CALPHAD applications, we consider a precision of ±4 kJ mol at⁻¹ to be acceptable with respect to the experimental value. This choice is of course somewhat arbitrary, but realistic. The calculations using PBE + U = 4.0 eV did not converge for CoO and, therefore, the enthalpy of formation of LiCoO₂ with reactions that use this compound could not be computed.

Fig. 1 Formation energy of (a) LiCoO₂ (b) o-LiMnO₂ and (c) LiNiO₂ according to the DFT-computed reaction and the x–c functional. The dashed line is the experimental value. The numbered label on each point is the absolute value of the energy difference between the result and the experimental value.^{10,15–17,33} The dotted lines represent the precision threshold level of ± 4 kJ mol at⁻¹ for CALPHAD applications.

Type 1 data are represented by triangles, Type 2 by squares and Type 3 by circles.

Type 1 data obtained with PBEsol, PW91, rSCAN, PBE or PBE + U with the lowest value of U (i.e. 2.6 eV for Co and 4.8 eV for Ni) show a very high degree of internal consistency as they differ by at most 1.2 kJ mol at⁻¹.

In contrast, each of the other two sets of data points, Type 2 and Type 3, clearly shows more scattering with respect to the x–c functional. While in the Type 2 data set, this scattering remains limited between 5 and 20 kJ mol at⁻¹ depending on the compound, and in the Type 3 data set, it becomes very large, ranging from 22 to 120 kJ mol at⁻¹.

As for the comparison to experimental data is concerned, for each of the three studied compounds, the experimental value falls more or less in the middle of the whole set of data points which covers an interval of width ranging from about 62 kJ mol at⁻¹ in the case of LiNiO₂ to 120 kJ mol at⁻¹ for LiMnO₂. In other words, an uninformed choice of the DFT-computed reaction, the x–c functional and the Hubbard U correction might result in an estimate of the LMO formation energy with an error that is most likely around ±30 kJ mol at⁻¹ and possibly up to ±60 kJ mol at⁻¹. This is indeed the picture that is shown in Table 1.

If however we focus on the Type 1 data points, we are led to two very different conclusions. First, in the case of LiCoO₂ and o-LiMnO₂, the agreement with experimental data is extremely good, with discrepancies ranging from 0.01 to 2.8 kJ mol at⁻¹. A small superior prediction capacity of PBEsol can be observed with a maximum difference of 1.6 kJ mol at⁻¹ w.r.t experiment. We note that the good quality of PBEsol could be expected since this functional is designed for crystals, whereas PBE and PW91 have more generic applications. Second, the discrepancy of about 24 kJ mol at⁻¹ in the case of LiNiO₂ already noted from the examination of Table 7 becomes very striking here, considering once again the consistency of the DFT data across the various x–c functionals.

The Type 2 data set displays a systematic discrepancy w.r.t. experiment which is much larger than the scattering associated to the x–c functional. For the three LiMO₂ compounds, this error is of the order of −20 to −30 kJ mol at⁻¹.

No systematic trend can be drawn for the Type 3 data set as the scattering across the x–c functional is large and the shift w.r.t. the experiment does not have a definite sign. The formation energies estimates however seem to order almost systematically with Type 2 being the lowest and Type 3 the highest.

We shall come back to this in the Discussion section. In conclusion, the results for LiCoO₂ and LiMnO₂ very clearly reveal the different precision levels obtained with the no-redox reactions compared to reactions involving electron transfer. While the latter lead to overestimates and underestimates with offsets from the experimental value in the order of several tens of kJ mol at⁻¹, the no-redox reactions consistently yield results within the ±4 kJ mol at⁻¹ acceptable uncertainty margin.

The influence of the Hubbard U over the enthalpy of formation of LiCoO₂ and LiNiO₂ is quite large among all tested reactions. We address this subject in the Discussion section.

Heat capacity

The heat capacity of LiMO₂ compounds was obtained within the HA and QHA. The HA only provides access to the heat capacity at constant volume (C_v). We calculate it using the static volume obtained at 0 K (C_v (T, V = V_static)). An important consequence is that pressure varies with temperature. Its value is out of reach of the HA but is an output of the QHA. For example, in the case of LiCoO₂, it reaches 2.5 GPa at 300 K and 7.5 GPa at 1000 K. Thus, C_v is arguably not the most appropriate quantity to compare with the experimental heat capacity, which, for solid state, is always obtained at constant pressure and varying volume. By employing the QHA, we obtain a heat capacity at constant pressure (C_p (T, p = 1 atm)). Both C_v and C_p are plotted along with experimental data and Neumann–Kopp approximation. The supercells used for phonon calculations are described in Tables S6, S7 and S8 (ESI†).

The calculated C_p for LiCoO₂ is in very good agreement with experimental data, with a maximal error of 10% around 150 K. The results are given in Fig. 2.

Fig. 2 Heat capacity of LiCoO₂ and comparison with experimental measurements.^10–14 A fit of experimental data was used to calculate the relative error of the calculated C_p. Two Neumann–Kopp approximations were derived: from Li₂O₂ + CoO and from Li₂O + CoO + Co₃O₄.⁹⁶

Regarding LiNiO₂ (Fig. 3), we are limited to experimental measurements at low temperatures and we find an error of about 17% at 300 K, which also motivates further investigations.

Fig. 3 Heat capacity of LiNiO₂ and comparison with experimental data¹⁸ and with a Neumann–Kopp approximation from oxides Li₂O₂ and NiO.⁹⁶

We have not found any report of experimental heat capacity measurements of LiMnO₂ in the literature. For the phonon calculations, we chose the P2/c structure. Some imaginary vibration frequencies were present for values of static pressure starting between −3 and −4 GPa, which corresponds to a temperature in the range 250–300 K. Imaginary frequencies in general can be due to either an imprecision in the calculations or a crystal instability. There is no definitive solution on how to deal with them.^97,98 In order to shed light on the origin of the problem we recomputed the phonon spectrum using the ABINIT code^90–92 using a linear response. The imaginary modes were still present, which lead us to the conclusion that the structure is unstable. To bypass this problem, we calculate the C_p of o-LiMnO₂ with a QHA proposed by Fleche,⁴⁷ in which a single phonon calculation at null pressure is necessary. The results are displayed in Fig. 4.

Fig. 4 Heat capacity of o-LiMnO₂ in P2/c symmetry calculated with a QHA from Fleche⁴⁷ and comparison with a two Neumann–Kopp approximation: from Li₂O + Mn₂O₃ and from Li₂O₂ + MnO.^57,96,99

Structural properties

Analogously to the enthalpy of formation, the cell parameters are properties measured at room temperature whereas the static DFT calculations correspond to 0 K and do not include zero-point vibrational (ZPV) corrections. We have corrected the ZPV and the volume expansion between 0 and 300 K as obtained with the QHA and we derived a specific mass. The results for LiCoO₂ are presented in Table 12.

Table 12 LiCoO₂ specific mass computed with the static and QHA methods at 300 K and difference with respect to the experimental value of Hertz et al.⁵⁰

Method	Specific mass (g cm^−3)	Difference (g cm^−3)	Relative error (%)
Static, 0 K and no ZPV	5.137	0.091	1.8
QHA, 300 K	5.048	0.002	0.04
Experimental value^50	5.046

---

The good agreement between the experimental and the QHA values shows that static DFT computations of the cell volume must be compared to a corrected experimental volume of 96.7–1.7 = 95.0 ų rather than to the exact experimental volume at 300 K. We use this corrected value to compare with the static volume obtained using other functionals and the results are presented in Fig. 5.

Fig. 5 Calculated cell volume of LiCoO₂ for different x–c and comparison with an estimated experimental static volume derived from ref. 50.

The calculated bulk modulus, the coefficient of volumetric thermal expansion (CTE) and the specific mass at 300 K are listed in Table 13.

Table 13 Structural properties of LiMO₂ end-members at 300 K obtained with the QHA

Compound	Bulk modulus (GPa)	CTE (K^−1)	Specific mass (g cm^−3)
LiCoO2	146.4	3.39 × 10^−5	5.048
o-LiMnO2	83.9	6.25 × 10^−5	4.213
LiNiO2	111.0	4.76 × 10^−5	4.774

---

The experimental bulk modulus determined by Hu et al.¹⁰⁰ for LiCoO₂ is 159.5 ± 2.2 GPa.

Discussion

Static formation energy

As illustrated by Fig. 1, the performance of DFT in computing the formation energies of LiCoO₂ and LiMnO₂ strongly depends on the type of DFT-computed reaction. Type 1 no-redox reactions give significantly better agreement between DFT and experiment than Type 2 and Type 3 reactions. Another striking feature displayed in Fig. 1 is that formation energies obtained with Type 1 reactions are much less dependent on the x–c functional than Type 3 data, while Type 2 data show intermediate sensitivity. A third point of relevance in this discussion is that, as shown by Table 14, the reaction energy of Type 1 reactions is much closer to zero than that of Type 2 or Type 3 reactions.

Table 14 Comparison of the reaction energy ΔE_R (in kJ per mol of reaction) of the DFT-computed reaction and the error δ_th-exp (in kJ mol at⁻¹) on the formation energy of LiMO₂

DFT-computed reaction	x–c functional	ΔER	δ th-exp
Li2O + Co3O4 → 2LiCoO2 + CoO (no-redox)	PBEsol	−73.6	1.6
	PW91	−67.8	2.4
	rSCAN	−64.2	2.8
	PBE	−67.4	2.4
Li2O2 + 2CoO → 2LiCoO2	PBEsol	−452.8	31.2
	PW91	−412.2	26.1
	rSCAN	−402.9	10.2
	PBE	−409.0	25.7
Li + Co + O2 → LiCoO2	PBEsol	−607.8	21.6
	PW91	−629.0	16.4
	rSCAN	−647.6	26.4
	PBE	−576.2	29.6
Li2O + Mn2O3 → 2LiMnO2 (no-redox)	PBEsol	−122.9	0.01
	PW91	−115.3	1.0
	rSCAN	−114.5	1.1
	PBE	−114.7	1.0
Li2O2 + 2MnO → 2LiMnO2	PBEsol	−503.0	29.7
	PW91	−464.6	24.9
	rSCAN	−435.2	21.3
	PBE	−463.0	24.7
Li + Mn + O2 → LiMnO2	PBEsol	−527.8	77.9
	PW91	−803.4	9.0
	rSCAN	−990.9	37.8
	PBE	−728.1	27.8

---

Lastly, Table 14 illustrates that the discrepancy δ_th-exp between calculation and experiment on the LiMO₂ formation energy Δ_rE_LiMO2 is much larger for Type 2 and Type 3 reactions than for Type 1 reactions. In summary, we point out a striking correlation between four properties of the DFT-computed reaction: (P1) constancy of the formal charges, (P2) insensitivity to the x–c functional, (P3) smallness of the energy of the DFT-computed reaction, and (P4) good agreement between calculation and experiment on Δ_rE_LiMO2. Type 1 and Type 3 reactions best illustrate the extreme cases, while Type 2 reactions are intermediate. Although we have not made any additional theoretical calculations to support this, this correlation lends itself to a simple qualitative interpretation based on the approximate character of the x–c functionals (w.r.t. the universal—and unknown—exact functional) and on a partition of the electron density among ions reflecting more or less their formal charge. In a Type 1 reaction, if we suppose that the sets of ions do not change much between reactants and products, the errors intrinsic to the x–c functionals on each ion separately compensate between reactants and products. The differences between the various functionals are hidden behind this compensation effect (P2) and the formation energy based on such a reaction is in good agreement with the experiment (P4). As for property P3, it is intuitively consistent with the similarities in ionic charges and coordination. In a Type 3 reaction, equivalent but opposite arguments may be given: large changes in formal charges reflect large electron density changes within each ion, so that the associated energy values differ according to the x–c functionals which implies that not all of these values can agree with the experimental value. What remains puzzling in this qualitative interpretation is that it is based on formal charges, which are difficult to reconcile with the details of the electronic structure, e.g. they hide the fact that bonding always includes a significant covalent contribution. In fact, the population analysis output by CASTEP in our calculations reveals changes in Mulliken (both charges and bond orders) and Hirshfeld populations during Type 1 reactions that are not necessarily smaller than in Type 2 reactions. Thus, the constancy of formal charges may just reflect a global compensation of various microscopic evolutions in ionic charge, covalency, bond lengths, etc, and this remains to be investigated more closely. As to the partitioning of the electron density among ions, it could be investigated using differences in charge density maps and Bader charges and this will be considered for future work.

In summary, we propose that the success of our methodology based on no-redox DFT-computed reactions relies on a compensation mechanism of the error due to the approximate character of the exchange–correlation (x–c) functional.

If we now go back to Fig. 1 and focus on Type 1 (triangles) data points, the question arises why agreement with experimental data is excellent for LiCoO₂ and LiMnO₂, with a discrepancy δ_th-exp ≈ 1 kJ mol at⁻¹, while in the case of LiNiO₂, the disagreement is very clear. Since stable oxides of nickel are limited to NiO with Ni(II), we first chose to base the DFT contribution on a reaction with Li₂O, NaNiO₂ and Na₂O, which does not involve any redox process. With this reaction, we have obtained an enthalpy of formation that is about 24 kJ mol at⁻¹ lower than the experimental value. We then studied another no-redox reaction with Li₂O, NiTiO₃ and Ti₂O₃ to help locate the source of this large discrepancy between DFT and experiment. Since this reaction leads to a very similar result as the previous one, we tested the accuracy of the proposed method in computing the formation energy of NiTiO₃, revealing a satisfactory result (Table 8) with a discrepancy of 0.19 kJ mol at⁻¹. This very good agreement, together with that obtained previously in cycles involving Li₂O, suggests LiNiO₂ as being the only source of the disagreement observed in Table 2. It is also worth noting that all the studies reported in Table 1 which calculated the LiNiO₂ formation energy using DFT and obtained satisfactory results used either an empirical method²⁴ or an adjustment of the Ni-3d Hubbard U value.²⁸ Thus, we conclude that our results motivate a reassessment of the experimental data for the LiNiO₂ heat of formation.

In our methodology, the choice of the no-redox DFT-computed LiMO₂ formation reaction is largely determined by the availability of experimental data for the formation energy of the other compounds taking part in the Hess’ cycle. While this might, at first sight, be considered as a weakness hindering a wide use of the methodology, in practice, this has so far not been a limitation. For instance, a straightforward no-redox reaction from binary oxides and LiCoO₂ is not possible, since Co₂O₃ is metastable, but we could easily bypass this issue by employing other stable binary oxides (CoO and Co₃O₄). Further, within this no-redox constraint, we have even been able to check that the choice of Hess’ cycle does not influence the final result. There remains to uncover the source of the discrepancy regarding the formation enthalpy of LiNiO₂, but overall we are confident that our methodology based on no-redox reactions deserved to be pursued further.

DFT+U results

Regarding the DFT+U results, we originally expected the influence of the U value on the Type 1 results to be small. This is true to some extent, given that the Hubbard correction shifts the energy of each of the compounds by several hundreds of kJ mol at⁻¹ while the resulting LiMO₂ formation energies only change by at most a few tens of kJ mol at⁻¹. But this remainder appears quite large on the scale of Fig. 1, which also displays the better results obtained without the Hubbard correction. We observe that the residual sensitivity to U is rather large compared to the precision required by our target application of CALPHAD modeling. Due to the uncertainty affecting the correct value of the input parameter U, this example also illustrates how difficult it is to improve the precision of DFT+U formation enthalpies below the average 5 kJ mol at⁻¹ level obtained by semi-empirical fitting methods such as those of Wang et al.³⁶ or Stevanovic et al.¹⁰¹ Further, the scattering of the LiCoO₂ Type 1 results in Fig. 1a, spread over about 25 kJ mol at⁻¹, may be due to some additional reason. We think that there may be an issue with Co₃O₄ since it contains both Co(II) and Co(III) whereas we assigned one and the same U value to all Co ions. Metastable electronic states arising with DFT+U are a known and complex problem.¹⁰² This possible issue with the computation of Co₃O₄ with GGA+U is out of the scope of the present work. Overall, the data in Fig. 1 illustrate a sensitivity of ΔE to U which, combined with the uncertainty about the value of U and the potential pitfall of metastable electronic states, may explain the rather poor performance of DFT+U in the present case. In comparison, the straightforward and parameter-free combination of standard DFT, PBEsol and no-redox reactions leads to a good agreement on thermodynamic and structural properties.

Vibrational contribution to the formation energy

The corrections due to zero point energy and temperature increments result in a formation energy of LiCoO₂ 0.76 kJ mol at⁻¹ closer to the average of the experimental data, while the discrepancy remains in the range of 1 to 4 kJ mol at⁻¹ depending on the source of experimental data. The smallness of this value is consistent with the widespread point of view that vibrational corrections to the formation energy are in generally below the so-called 1 kcal mol at⁻¹ chemical accuracy, but we think that here it is mostly due to the compensation effect behind the use of a no-redox reaction. Nevertheless, considering that the uncertainty in the experimental measurement of the LiCoO₂ enthalpy of formation was estimated by their authors to be 0.7 kJ mol at⁻¹,¹⁰ we conclude that the ZPE and the vibrational contribution at 300 K obtained within the HA make up a relevant and valuable complement to static energy estimation. The additional precision brought by the QHA seems unnecessary for the Type 1 DFT-computed reaction.

Structural properties

Even though the QHA does not contribute significantly to the formation energies obtained with Type 1 Hess’ cycles, it remains essential to reach a proper description of the structural data. For example, the LiCoO₂ specific mass output by the QHA model at 300 K has a discrepancy of only 0.002 g cm⁻³, as compared to 0.091 g cm⁻³ when only the static method is used (Table 12). Among the studied x–c functionals, PBEsol provided the most accurate specific mass. We recall that the differences in cell volumes obtained with various x–c functionals are of the same order of magnitude as the difference between the static volume and the QHA estimate at 300 K, so that a meaningful comparison with experimental data can only be done if the latter correction is taken into account, as we did in Fig. 5. Within the QHA, we calculate the bulk modulus B(300 K) of LiCoO₂ and an underestimation of 8% in relation to the experimental value was found.

Heat capacity

The QHA is also essential for heat capacity data at higher temperatures, as shown by the differences between C_p and C_v (Fig. 2, Fig. 3 and Fig. 4). Our computation method for the heat capacity could be validated with LiCoO₂ for which numerous measurements are available from 0 to 1200 K. As for LiNiO₂, the heat capacity agrees very well with experimental data in a temperature range of approximately 0–200 K. An increasing deviation is observed between 200 K and 300 K, the maximum temperature at which experimental C_p data are available. The maximal divergence occurs at 300 K, with a relative error of 17.2%. Our phonon data allowed us to compute C_p up to 1200 K. The discrepancy at 300 K requires a validation that will require new measurements to be made, in order to give more confidence in our prediction for the higher temperatures. As for o-LiMnO₂ for which no experimental data are available, we compute a C_p curve that might be employed as an entry in a thermodynamic database. We note however that the presence of unstable vibration modes beyond a specific volume corresponding to a temperature of about 300 K also needs experimental validation. Note that a significant structural cationic disorder has been observed in LiMnO₂ by Kellerman et al.⁵¹ and the thermodynamic contributions arising from this disorder remain to be assessed. For the CALPHAD application, as a first approximation, it is reasonable to assume that the C_p value of r-LiMnO₂ is the same as that obtained for o-LiMnO₂. Such an assumption is in fact a kind of Neumann–Kopp approximation, more precise than that obtained from binary oxides or pure elements. For LiMnO₂, the calculated C_p agrees with the two different Neumann–Kopp approximations proposed. For LiCoO₂, however, the calculated C_p value is in clear disagreement with two Neumann–Kopp approximations based on binary oxides, even with the one based on the Type 1 no-redox reaction, and these two N–K approximations also disagree with each other. Thus, more investigations are needed on the heat capacity of these compounds. This is also a reminder that the N–K approximation should always be used with caution.

Conclusion

In the present study, a new methodology for computing the structural and thermodynamic properties of LiMO₂ (M = Co, Mn and Ni) was proposed and assessed using experimental results available for comparison. This paves the way for the study of the NMC system of which they are end-members. The properties were computed by means of DFT and using three levels of approximation: static, harmonic and quasi-harmonic.

For each of these compounds, we compared several Hess’ cycles involving different DFT-computed formation reactions. The best performance is clearly obtained using a Hess’ cycle in which no redox process occurs in the DFT-computed reaction involving the target compound LiMO₂ (M = Co, Ni, Mn). These no-redox reactions consistently give the best accuracy for the formation energy, between 0.01 and 1.6 kJ mol at⁻¹ for LiCoO₂ and o-LiMnO₂. This scheme naturally avoids the use of DFT estimates of molecular O₂. This very high precision and the fact that the same results are reproduced with four different x–c functionals enabled us to uncover an unexplained and large discrepancy of 24 kJ mol at⁻¹ for the LiNiO₂ formation energy, confirmed by reproducing the calculation with two different Hess cycles. This motivates further investigation on the experimental value of the LiNiO₂ standard enthalpy of formation. The enthalpy of formation of the metastable rhombohedral LiMnO₂ phase was computed and it is about 1.0 kJ mol at⁻¹ more positive than that of o-LiMnO₂.

The harmonic vibrational contribution to the formation energy brings the result 0.76 kJ mol at⁻¹ closer to the LiCoO₂ experimental value, resulting in a prediction in line with the level of accuracy required for thermodynamic CALPHAD calculations. Our results also show that, even though the contributions obtained with the QHA are not significant to formation energy computations when using no-redox DFT-computed reactions, they are relevant for the best agreement with structural experimental data at 300 K and the heat capacity obtained at temperatures above 300 K. We found that PBEsol performed best for all the properties computed in this study. A very good agreement between calculated and experimental heat capacities at constant pressure (P = 1 atm) was obtained for LiCoO₂, enabling us to propose estimates for the heat capacity of LiNiO₂ at high temperatures as well as the heat capacity of o-LiMnO₂, hitherto unknown experimentally. As for the methodology to formation energy calculations using no-redox DFT-computed reactions, which gives very good results for LiCoO₂ and LiMnO₂, the question arises whether it might be generalized. It may be considered in a wide array of problems provided formal charges may be defined and a reaction without changes of formal charges can be found, involving compounds with known experimental enthalpies of formation. The definition of formal charges is not universally possible, but it is possible at least for the large family of oxides. While we have no theoretical support for this methodology, we trust that the literature will eventually provide more data on which to base an assessment of this methodology.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We are grateful to Dr Jean-Louis Flèche for providing us with his QHA code. This work was supported by the Cross-cutting basic research Program (RTA Program) of the CEA Energy Division.

References

1. N. Zhang, J. Li, H. Li, A. Liu, Q. Huang, L. Ma, Y. Li and J. R. Dahn, Chem. Mater., 2018, 30, 8852–8860.
2. W. Zhang, D. M. Cupid, P. Gotcu, K. Chang, D. Li, Y. Du and H. J. Seifert, Chem. Mater., 2018, 30, 2287–2298.
3. N. Li, D. Li, W. Zhang, K. Chang, F. Dang, Y. Du and H. J. Seifert, Prog. Nat. Sci.: Mater. Int., 2019, 29, 265–276.
4. A. A. Luo, Calphad, 2015, 50, 6–22.
5. K. Chang, B. Hallstedt, D. Music, J. Fischer, C. Ziebert, S. Ulrich and H. J. Seifert, Calphad, 2013, 41, 6–15.
6. K. Chang, B. Hallstedt and D. Music, Calphad, 2012, 37, 100–107.
7. R. Robert, C. Villevieille and P. Novák, J. Mater. Chem. A, 2014, 2, 8589.
8. B. Ammundsen and J. Paulsen, Adv. Mater., 2001, 13, 943–956.
9. M. Bianchini, M. Roca-Ayats, P. Hartmann, T. Brezesinski and J. Janek, Angew. Chem., Int. Ed., 2019, 58, 10434–10458.
10. P. Gotcu-Freis, D. M. Cupid, M. Rohde and H. J. Seifert, J. Chem. Thermodyn., 2015, 84, 118–127.
11. H. Kawaji, M. Takematsu, T. Tojo, T. Atake, A. Hirano and R. Kanno, J. Therm. Anal. Calorim., 2002, 68, 833–839.
12. A. L. Emelina, M. A. Bykov, M. L. Kovba, B. M. Senyavin and E. V. Golubina, Russ. J. Phys. Chem., 2011, 85, 357–363.
13. M. Ménétrier, D. Carlier, M. Blangero and C. Delmas, Electrochem. Solid-State Lett., 2008, 11, A179–A182.
14. O. Jankovský, J. Kovařík, J. Leitner, K. Ruzicka and D. Sedmidubský, Thermochim. Acta, 2016, 634, 26–30.
15. M. Wang, A. Navrotsky, S. Venkatraman and A. Manthiram, J. Electrochem. Soc., 2005, 152, J82–J84.
16. M. Masoumi, D. M. Cupid, T. L. Reichmann, K. Chang, D. Music, J. M. Schneider and H. J. Seifert, Int. J. Mater. Res., 2017, 108, 869–878.
17. M. Wang and A. Navrotsky, Solid State Ionics, 2004, 166, 167–173.
18. H. Kawaji, Solid State Ionics, 2002, 152–153, 195–198.
19. P. Strobel, J.-P. Levy and J.-C. Joubert, J. Cryst. Grow., 1984, 66, 257–261.
20. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
21. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
22. L. Wu, W. H. Lee and J. Zhang, Mater. Today: Proc., 2014, 1, 82–93.
23. T. Du, B. Xu, M. Wu, G. Liu and C. Ouyang, J. Phys. Chem. C, 2016, 120, 5876–5882.
24. K. Chang, B. Hallstedt and D. Music, Chem. Mater., 2012, 24, 97–105.
25. H. Xu, W. Xiao, Z. Wang, J. Hu and G. Shao, J. Energy Chem., 2021, 59, 229–241.
26. D. Kramer and G. Ceder, Chem. Mater., 2009, 21, 3799–3809.
27. R. C. Longo, F. T. Kong, S. Kc, M. S. Park, J. Yoon, D.-H. Yeon, J.-H. Park, S.-G. Doo and K. Cho, Phys. Chem. Chem. Phys., 2014, 16, 11233–11242.
28. M. Tuccillo, O. Palumbo, M. Pavone, A. B. Muñoz-García, A. Paolone and S. Brutti, Crystals, 2020, 10, 526.
29. N. T. Tsebesebe, K. Kgatwane, R. Ledwaba and P. Ngoepe, J. Phys.: Conf. Ser., 2022, 2298, 012010.
30. E. Lee, K.-R. Lee and B.-J. Lee, Comput. Mater. Sci., 2018, 142, 47–58.
31. K. Refson, P. R. Tulip and S. J. Clark, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 155114.
32. M. Masoumi, Dissertation, KIT, 2019.
33. M. Wang and A. Navrotsky, J. Solid State Chem., 2005, 178, 1230–1240.
34. D. C. Patton, D. V. Porezag and M. R. Pederson, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 7454–7459.
35. L. Wang, T. Maxisch and G. Ceder, Chem. Mater., 2007, 19, 543–552.
36. L. Wang, T. Maxisch and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 195107.
37. A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer, K. A. Persson and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 045115.
38. C. Franchini, R. Podloucky, J. Paier, M. Marsman and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 195128.
39. F. Zhou, M. Cococcioni, C. A. Marianetti, D. Morgan and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 235121.
40. A. Chakraborty, S. Kunnikuruvan, S. Kumar, B. Markovsky, D. Aurbach, M. Dixit and D. T. Major, Chem. Mater., 2020, 32, 915–952.
41. D. K. Chakrabarty, An Introduction to Physical Chemistry, Alpha Science International, Limited, 2001.
42. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
43. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
44. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
45. A. P. Bartók and J. R. Yates, J. Chem. Phys., 2019, 150, 161101.
46. M. Born and K. Huang, Am. J. Phys., 1955, 23, 474.
47. J. L. Fleche, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 245116.
48. A. Togo and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 184104.
49. M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov and H. Wondratschek, Z. Kristallogr. - Cryst. Mater., 2006, 221, 15–27.
50. J. T. Hertz, Q. Huang, T. McQueen, T. Klimczuk, J. W. G. Bos, L. Viciu and R. J. Cava, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 075119.
51. D. G. Kellerman, J. E. Medvedeva, V. S. Gorshkov, A. I. Kurbakov, V. G. Zubkov, A. P. Tyutyunnik and V. A. Trunov, Solid State Sci., 2007, 9, 196–204.
52. G. M. Kale, S. S. Pandit and K. T. Jacob, Trans. JIM, 1988, 29, 125–132.
53. W. L. Roth, J. Phys. Chem. Solids, 1964, 25, 1–10.
54. B. J. Boyle, E. G. King and K. C. Conway, J. Am. Chem. Soc., 1954, 76, 3835–3837.
55. I. Barin, O. Knacke and O. Kubaschewski, Thermochemical properties of inorganic substances, Springer, Berlin, Heidelberg, 1st edn, 1991.
56. W. L. Roth, Phys. Rev., 1958, 110, 1333–1341.
57. K. T. Jacob, A. Kumar, G. Rajitha and Y. Waseda, High Temp. Mater. Processes, 2011, 30, 459–472.
58. R. Shivaramaiah, S. Tallapragada, G. P. Nagabhushana and A. Navrotsky, J. Solid State Chem., 2019, 280, 121011.
59. C. Darie, P. Bordet, S. de Brion, M. Holzapfel, O. Isnard, A. Lecchi, J. E. Lorenzo and E. Suard, Eur. Phys. J. B, 2005, 43, 159–162.
60. O. Kubaschewski, C. B. Alcock and P. J. Spencer, Materials Thermochemistry, Pergamon, New York, 6th edn, 1993.
61. G. Shirane, S. J. Pickart and Y. Ishikawa, J. Phys. Soc. Jpn., 1959, 14, 1352–1360.
62. G. L. Humphrey, J. Am. Chem. Soc., 1951, 73, 1587–1590.
63. P. H. Carr and S. Foner, J. Appl. Phys., 1960, 31, S344–S345.
64. G. K. Johnson, R. T. Grow and W. N. Hubbard, J. Chem. Thermodyn., 1975, 7, 781–786.
65. K. Chang and B. Hallstedt, Calphad, 2011, 35, 160–164.
66. D. Hobbs, J. Hafner and D. Spišák, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 014407.
67. M. J. Redman and E. G. Steward, Nature, 1962, 193, 867.
68. Y. C. Sui, Y. Zhao, J. Zhang, S. Jaswal, X. Z. Li and D. J. Sellmyer, IEEE Trans. Magn., 2007, 43, 3115–3117.
69. A. Ghosh, S. Jana, M. K. Niranjan, F. Tran, D. Wimberger, P. Blaha, L. A. Constantin and P. Samal, J. Phys. Chem. C, 2022, 126, 14650–14660.
70. W. Jauch, M. Reehuis, H. J. Bleif, F. Kubanek and P. Pattison, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 052102.
71. R. Hoppe, G. Brachtel and M. Jansen, Z. Anorg. Allg. Chem., 1975, 417, 1–10.
72. G. Dittrich and R. Hoppe, Z. Anorg. Allg. Chem., 1969, 368, 262–270.
73. C. Pouillerie, E. Suard and C. Delmas, J. Solid State Chem., 2001, 158, 187–197.
74. J. Cao, H. Zou, C. Guo, Z. Chen and S. Pu, Solid State Ionics, 2009, 180, 1209–1214.
75. A. Rougier, C. Delmas and A. V. Chadwick, Solid State Commun., 1995, 94, 123–127.
76. J.-H. Chung, Th Proffen, S. Shamoto, A. M. Ghorayeb, L. Croguennec, W. Tian, B. C. Sales, R. Jin, D. Mandrus and T. Egami, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 064410.
77. H. Das, A. Urban, W. Huang and G. Ceder, Chem. Mater., 2017, 29, 7840–7851.
78. C. A. Marianetti, D. Morgan and G. Ceder, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 224304.
79. H. Chen, C. L. Freeman and J. H. Harding, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 085108.
80. Z. Chen, H. Zou, X. Zhu, J. Zou and J. Cao, J. Solid State Chem., 2011, 184, 1784–1790.
81. S. Sicolo, M. Mock, M. Bianchini and K. Albe, Chem. Mater., 2020, 32, 10096–10103.
82. M. Manago, G. Motoyama, S. Nishigori, K. Fujiwara, K. Kinjo, S. Kitagawa, K. Ishida, K. Akiba, S. Araki, T. C. Kobayashi and H. Harima, J. Phys. Soc. Jpn., 2022, 91, 113701.
83. A. Pulkkinen, B. Barbiellini, J. Nokelainen, V. Sokolovskiy, D. Baigutlin, O. Miroshkina, M. Zagrebin, V. Buchelnikov, C. Lane, R. S. Markiewicz, A. Bansil, J. Sun, K. Pussi and E. Lähderanta, Phys. Rev. B, 2020, 101, 075115.
84. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, Z. Kristallogr. - Cryst. Mater., 2005, 220, 567–570.
85. Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett., 1997, 78, 4063–4066.
86. J. H. Lloyd-Williams and B. Monserrat, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 184301.
87. M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Rev. Mod. Phys., 1992, 64, 1045–1097.
88. Materials Studio, BIOVIA, Dassault Systèmes (version 20.1.0.2728) 2020.
89. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
90. A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applencourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, F. Bruneval, G. Brunin, D. Caliste, M. Côté, J. Denier, C. Dreyer, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, F. Jollet, G. Jomard, A. Martin, H. P. C. Miranda, F. Naccarato, G. Petretto, N. A. Pike, V. Planes, S. Prokhorenko, T. Rangel, F. Ricci, G.-M. Rignanese, M. Royo, M. Stengel, M. Torrent, M. J. van Setten, B. Van Troeye, M. J. Verstraete, J. Wiktor, J. W. Zwanziger and X. Gonze, J. Chem. Phys., 2020, 152, 124102.
91. M. Torrent, F. Jollet, F. Bottin, G. Zérah and X. Gonze, Comput. Mater. Sci., 2008, 42, 337–351.
92. X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, N. Brouwer, F. Bruneval, G. Brunin, T. Cavignac, J.-B. Charraud, W. Chen, M. Côté, S. Cottenier, J. Denier, G. Geneste, P. Ghosez, M. Giantomassi, Y. Gillet, O. Gingras, D. R. Hamann, G. Hautier, X. He, N. Helbig, N. Holzwarth, Y. Jia, F. Jollet, W. Lafargue-Dit-Hauret, K. Lejaeghere, M. A. L. Marques, A. Martin, C. Martins, H. P. C. Miranda, F. Naccarato, K. Persson, G. Petretto, V. Planes, Y. Pouillon, S. Prokhorenko, F. Ricci, G.-M. Rignanese, A. H. Romero, M. M. Schmitt, M. Torrent, M. J. van Setten, B. Van Troeye, M. J. Verstraete, G. Zérah and J. W. Zwanziger, Comput. Phys. Commun., 2020, 248, 107042.
93. F. Bottin, S. Leroux, A. Knyazev and G. Zérah, Comput. Mater. Sci., 2008, 42, 329–336.
94. C. Lee and X. Gonze, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 8610–8613.
95. M. Chen, B. Hallstedt and L. J. Gauckler, J. Phase Equilib., 2003, 24, 212–227.
96. P. J. Linstrom and W. G. Mallard, J. Chem. Eng. Data, 2001, 46, 1059–1063.
97. G. J. Ackland, M. C. Warren and S. J. Clark, J. Phys.: Condens. Matter, 1997, 9, 7861.
98. C.-C. Lee, C.-E. Hsu and H.-C. Hsueh, J. Phys.: Condens. Matter, 2020, 33, 055902.
99. J. L. Shapiro, B. F. Woodfield, R. Stevens, J. Boerio-Goates and M. L. Wilson, J. Chem. Thermodyn., 1999, 31, 725–739.
100. Y.-Q. Hu, L. Xiong, X.-Q. Liu, H.-Y. Zhao, G.-T. Liu, L.-G. Bai, W.-R. Cui, Y. Gong and X.-D. Li, Chinese Phys. B, 2019, 28, 016402.
101. V. Stevanović, S. Lany, X. Zhang and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 115104.
102. A. B. Shick, W. E. Pickett and A. I. Liechtenstein, J. Electron Spectrosc. Relat. Phenom., 2001, 114–116, 753–758.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp01771k

---