Discovery of zirconium dioxides for the design of better oxygen-ion conductors using efficient algorithms beyond data mining

Abstract

It is important to find crystal structures with low formation (E_v) and migration-barrier (E_m) energies for oxygen vacancies for the development of fast oxygen-ion conductors. To identify crystal structures with lower E_v and E_m than those of ground-state ZrO₂, we first reoptimize the crystal structures of various oxides reported in the database, and then directly construct them using an evolutionary algorithm. For efficient searching, we employ the linearized ridge regression model for E_v using descriptors obtained from density functional theory calculations of the unit cells and apply the predicted E_v as a fitness value in the evolutionary algorithm. We also find a correlation between the E_v and E_m for the crystal structures of ZrO₂. On the basis of this correlation, we confirm that the newly constructed crystal structures, as well as certain reoptimized structures from the database, that possess low E_v also have E_m lower than that of ground-state ZrO₂. Our successful strategy consisting of a combination of the evolutionary algorithm, first-principles calculations, and machine-learning techniques may be applicable to other oxide systems in finding crystal structures with low E_v and E_m.

---

1. Introduction

A high oxygen ion conductivity (σ_O) is an important property for applications^1,2 such as the electrolytes of solid-oxide fuel cells (SOFC), oxygen separation membranes, and gas sensors. Low vacancy formation (E_v) and migration-barrier (E_m) energies of oxygen vacancies V_O are favorable for achieving a higher σ_O.³ Currently, Y-doped ZrO₂ (YSZ) is widely used because of its advantages such as abundance, chemical stability, non-toxicity, and low cost. This material shows the σ_O of ∼10⁻² S cm⁻¹ at the high temperature of 1000 K.⁴ For useful industrial applications, it is necessary to have similar σ_O values at lower temperatures or a higher σ_O at a similar temperature. Indeed, several oxides such as Gd-doped CeO₂ (GDC),⁵ pure or Er-doped δ-phase Bi₂O₃,^6,7 and Sr- and Mg-doped LaGaO₃ (LSGM)^8,9 have been reported to show higher σ_O than YSZ at the same temperature. However, the need for the development of new oxygen-ion conductors remains.

The type of crystal structure strongly affects the σ_O as well as E_v and E_m. The monoclinic ground-state structure of ZrO₂ (space group P2₁/c) is known to have a low σ_O. However, the σ_O can be significantly increased through phase transformation into either the tetragonal (space group P4₂/nmc) or cubic fluorite structure (space group Fm3̄m) upon doping with Y.⁵ The δ-phase Bi₂O₃ in the defective fluorite structure (its high-temperature phase) is known to have a high σ_O, whereas its α-phase counterpart is a low-temperature phase with a low σ_O.⁶ Reported crystal structures with high σ_O are mainly fluorites, perovskites,⁵ and their derivatives, such as the defective fluorite and melilite¹⁰ structures. However, we expect that other undiscovered crystal structures may have high σ_O values.

The evolutionary algorithm is useful for the construction and search of crystal structures.¹¹ This has been applied in combination with first-principles calculations for determining crystal structures with maximized or minimized functional properties (fitness values) such as enthalpy at high pressure,^12–14 hardness,¹⁵ and dielectric constant.¹⁶ This algorithm is especially powerful for easily obtained fitness values, because it requires many individual crystal structures. From this perspective, the direct usage of the E_v or E_m as a fitness value may be inefficient because these properties should be obtained from supercell calculations.¹⁷ In this case, one useful method is machine-learning techniques.

A regression analysis, one of the machine-learning techniques, is used to predict a target variable on the basis of already accumulated data, which can be used for descriptors (also called “predictor variables” or “representations”). If the descriptors are composed of “density functional theory (DFT)-unit-cell descriptors” that contain information regarding the unit-cells obtained only from DFT calculations,^18,19 the prediction model can be constructed with the advantage of computational efficiency. Some regression models have employed descriptors obtained from DFT calculations of the unit-cells to predict target variables such as the melting temperatures of binary metals,²⁰ GW-level band gaps for inorganic compounds,²¹ and interatomic potentials.²² In addition, Deml et al.²³ constructed a prediction equation for E_v using ordinary least-square regression (OLSR) with four unit-cell descriptors and 45 oxides. However, it is not yet guaranteed whether this prediction model maintains accuracy in predicting E_v for various types of crystal structures.

The purpose of the present study is to propose an efficient method to find various crystal structures with lower E_v and E_m values than those of ground-state ZrO₂. To this end, we reoptimized the reported crystal structures of various oxides in the Materials Project Database (MPD),²⁴ one of the inorganic materials databases based on first-principles calculations, and reconsidered their E_v and E_m. We constructed a prediction model for E_v based on linearized regression analyses using DFT-unit-cell descriptors and subsequently constructed crystal structures based on the predicted E_v using the evolutionary algorithm. We also investigated the relationship between E_v and E_m for the crystal structures of ZrO₂. Finally, we found crystal structures having E_v and E_m values lower than those of ground-state ZrO₂.

2. Methodology

2.1. Overview

First, we reoptimized the crystal structures of several oxides from the MPD²⁴ for ZrO₂. The computed E_v and other properties of the unit-cells are used as the target variable and descriptors, respectively, for the regression analyses.

Fig. 1 shows the workflow for the construction of the crystal structures and prediction of the E_v. New structures were constructed using the evolutionary algorithm. First-principles calculations were performed to optimize the constructed crystal structures and extract the DFT-unit-cell descriptors. The regression analyses were then applied to obtain the predicted E_v. In total, 462 crystal structures were constructed through the repetition of cycles comprising the construction of crystal structures, calculation of DFT-unit-cell descriptors from the constructed structures, and application of regression analysis to predict E_v.

Fig. 1 Flowchart of screening for crystal structures of ZrO₂ with low E_v constructed by evolutionary algorithm and the regression analysis. The procedures in the large rectangular box are repeated to construct 462 crystal structures of ZrO₂.

We identified crystal structures with predicted E_v in the target range and performed structural optimization again under tighter conditions. The accurate E_v of the remaining crystal structures were computed from their supercells. Finally, we obtained 14 newly constructed crystal structures that have both the predicted and computed E_v.

More details regarding each method are found in the next subsection.

2.2. Crystal structures for the training set in regression analysis

As mentioned in the previous subsection, we collected various crystal structures of oxides with cation-to-O ratios of 1 : 2 in the MPD²⁴ for use as training data in regression analysis and for considering their E_v and E_m. We set the number of atoms in the unit-cells to ≤30. We performed first-principles calculations for structural optimization and computation of the E_v. To identify crystal structures with reliable E_v, we adopted several filters. Firstly, the crystal structures must satisfy the convergence of electronic and ionic relaxations. Secondly, they should be non-metallic. Thirdly, they should maintain a displacement of the center of the V_O site of 0.1 Å at most after structural relaxations; this is for avoiding the convergence of crystal structures into other crystal structures. On the basis of these filters, we collected 16 crystal structures of ZrO₂ with the computed E_v, as summarized in Table 1. Henceforth, this set of 16 crystal structures is called “reoptimized-ZrO₂” and each structure in the set is called by its space group type. Excluded crystal structures are summarized in Table S1 in ESI.† All the crystal structures in this study are displayed using the VESTA program.²⁵

Table 1 Computed properties, including E_v, for the 16 crystal structures of reoptimized-ZrO₂

Space group type	Oxide reported in MPD^24	E − E(ground-state) (eV per atom)	Eg^GGA+U (eV)	Computed Ev (eV)^a
a Values in the parentheses are the minimum and the maximum values in sequence, which depend on the kinds of O sites.
P21/c	ZrO2	0.000	3.53	6.15 (6.13, 6.16)
Pbca	ZrO2	0.007	3.52	6.15 (6.14, 6.15)
I41/amd	ZrO2	0.013	3.89	6.40
P21/m	TiO2	0.017	3.85	6.40 (6.38, 6.41)
C2/c	TiO2	0.018	3.85	6.38
Pca21	ZrO2	0.019	3.78	5.93 (5.93, 5.93)
P41212	SiO2	0.021	3.51	6.71
P42/mnm	ZrO2	0.026	3.55	6.68
Pnma	TiO2	0.028	4.11	6.06 (6.01, 6.10)
P42/nmc	ZrO2	0.029	3.89	5.81
Pbcn	ZrO2	0.032	3.85	5.82
Fm3̄m	ZrO2	0.037	3.38	5.84
P4/n	ZrO2	0.054	3.34	5.71 (5.35, 5.80)
Pna21	SiO2	0.116	2.80	5.14 (4.97, 5.30)
R3̄	SiO2	0.126	3.75	5.93 (5.87, 6.00)
P63mc	TiO2	0.179	3.75	6.51 (6.19, 6.76)

---

2.3. Regression analysis

The linearized ridge regression (RR)²⁶ model uses a minimization function of the OLSR with a L2-norm penalty term that is given bywhere β is an n-dimensional vector of the regression coefficients of the descriptors, X is an (n × p) descriptor matrix, y is an n-dimensional vector of the target property for the training set, and λ is a coefficient of the penalty term. Here, the penalty term is used for controlling the coefficients to avoid overfitting.

Table 2 shows the list of the DFT-unit-cell descriptors. To predict E_v of the crystal structures of ZrO₂, RR was employed with all 11 descriptors. The OLSR with three descriptors (x₀₁–x₀₃), suggested by Deml et al.,²³ was also used for comparison. Henceforth, the former and latter prediction models are referred to as RR-11-descriptors and OLSR-3-descriptors, respectively.

Table 2 DFT-unit-cell descriptors used for regression analyses

Name	Unit-cell-descriptor	Note
x01	Formation energy (ΔEf)	Relative energy from the Zr metal and the half-energy of O2 gas molecule
x02	Band-gap (Eg^GGA+U)	From the electronic DOS
x03	Center of O 2p-band (EO2p)	From the minimum energy of the electronic DOS to the VBM
x04	Volume per atom
x05	Average bond length (rZr–O)	Cutoff radius of 2.4 Å
x06	Coordination number of O (CNO)	Cutoff radius of 2.4 Å
x07	Mean of |x02 − x02(ground-state)|
x08	Mean of |x03 − x03(ground-state)|
x09	Mean of |x04 − x04(ground-state)|
x10	Mean of |x05 − x05(ground-state)|
x11	Mean of |x06 − x06(ground-state)|

---

The training and test data were randomly divided in the ratio of 75% to 25%. To avoid overfitting, three-fold cross-validation (CV) was also used. The prediction errors were defined as root-mean-square-errors (RMSE) that were averaged for 30 different random samplings of the training set.

2.4. Structure construction using an evolutionary algorithm

We used the USPEX^11,27,28 code to generate crystal structures based on an evolutionary algorithm. Here, the predicted E_v was used as a fitness value for minimizing the quantity |predicted E_v − 5 eV|. This form was applied to focus on the search for crystal structures with E_v values moderately lower than that (6.15 eV) of the ground-state P2₁/c structure.

Each newly constructed unit-cell of ZrO₂ was set to have four Zr and eight O atoms. Sixty first-generation crystal structures were produced using randomly selected space group types. The fitness value was obtained after structural relaxation based on the first-principles calculations. From the second generation onwards, new crystal structures were produced by genetic operators, namely, heredity (30%), random symmetric algorithm (30%), and mutation (40%). Each generation consisted of 40 crystal structures. The evolutionary algorithm was terminated if the best-ranked crystal structure was not changed over eight generations or after the number of generations reached 30.

2.5. First-principles calculations

All first-principles calculations were performed using the projector augmented-wave (PAW)^29,30 method implemented in the Vienna Ab-initio Simulation Package (VASP)^31,32 within the framework of the Generalized Gradient Approximation (GGA) of Perdew–Burke–Ernzerhof (PBE) form,³³ including on-site Coulomb interaction³⁴ with an effective U − J of 3 eV (ref. 35) (GGA+U) for the d-orbitals of Zr. Here, we employed the GGA+U rather than the GGA to decrease the number of metallic crystal structures constructed by the evolutionary algorithm (See Fig. S1 in ESI†). The valence electrons occupied the 4s, 4p, 5s, and 4d orbitals for Zr, and the 2s and 2p orbitals for O.

During the construction of crystal structures by the evolutionary algorithm, structural relaxations of the unit-cells (with the associated changes in lattice constants and atomic coordinates) were modeled until the interatomic force on each atom was reduced to within 0.01 eV Å⁻¹. To avoid bad convergence, we constrained the number of ionic iteration cycles to the maximum of 120. The cutoff energy was set to 400 eV. The Brillouin zone was sampled by Γ-centered meshes with the density of 2π × 0.12 Å⁻¹ for structural optimization such as fast screening. The calculations for the total energy and electronic density of states (DOS) of the optimized unit-cells were performed with a finer sampling of 2π × 0.08 Å⁻¹ of the Brillouin zone. Subsequently, structural optimizations were performed for the survivors of the evolutionary algorithm again, but in tighter computational conditions: the cutoff energy was increased to 500 eV, and the Brillouin zone was sampled using Γ-centered 8 × 8 × 8 meshes.

The computed E_v were obtained using the supercell method¹⁷ as follows,

where is the energy of a supercell including a V_O, E(ZrO₂) is the energy of a supercell of pure ZrO₂, 1/2E(O₂) is the half-energy of an O₂ gas molecule (used as the chemical potential of O for the O-rich condition), q is the unit charge, E_VBM is the energy of the valence band maximum (VBM), and E_Fermi is the Fermi level, that is a variable changes between the VBM and the conduction band minimum (CBM). The chemical potential of O depends on the processing condition, but it provides a constant shift to the E_v of all of the crystal structures of ZrO₂ and the general consequence related to the ranks of the E_v are unchanged. Therefore, we only considered the O-rich condition in this study (See other O-rich conditions in Table S2 in ESI†). We calculated both of the E_v for a neutral oxygen vacancy and a doubly charged oxygen vacancy . For a prediction model based on the evolutionary algorithm and regression analysis, was used as the computed E_v. The calculations for the supercell with a V_O were performed until the interatomic forces on each atom were reduced to within 0.02 eV Å⁻¹ with fixed lattice constants. The sizes of the supercells were set such that the lengths of the three lattice constants were ∼10 Å. Γ-centered 2 × 2 × 2 meshes were used to describe the k-space. The types of O sites were confirmed by SPGLIB implemented in PHONOPY.^36,37 For a crystal structure containing more than two types of O sites, the computed E_v values were averaged.

The minimum energy path (MEP) for a V_O migration was obtained by the climbing-images nudged-elastic-band (CI-NEB) method.^38,39 The E_m of a V_O can be defined as the energy difference between the transition state and the initial or final state from the obtained MEP, as shown in Fig. 2. The CI-NEB calculations were performed with three intermediate images (states) until the forces decreased below 0.03 eV Å⁻¹ with a spring constant of 5 eV Å⁻² between images. Here, the climbing images with even numbers were used to increase the probability of finding the transition state. We employed a for E_m, assuming that V_O were formed by the reduction of cation valency upon doping. The supercells with a were described by decreasing two electrons of the background charge for compensation. The internal coordinates for the supercells with a were relaxed in the same computational conditions as those for .

Fig. 2 Energy differences of the intermediate states (images) of the P4₂/nmc structure of ZrO₂, which were obtained by the CI-NEB method. The atomic coordinates in the figure were shown from the normal direction to the (001) plane. Green and light gray spheres denote Zr and O atoms, respectively. A brown sphere is used to emphasize the migrating O atom along the 〈110〉 path. A pink circle denotes a V_O site. The V_O at the initial or final state is formed by removing an O atom and breaking four O–Zr bonds. The E_m is defined as the energy difference between the transition state and the initial or final state. The migrating O atom has four O–Zr bonds at the initial state, while it maintains only two O–Zr bonds at the transition state and breaks the other two O–Zr bonds.

3. Results and discussion

3.1. Construction of a prediction model for E_v using DFT-unit-cell descriptors

As mentioned in the Introduction, Deml et al.²³ constructed a linear prediction model for the E_v with four unit-cell descriptors, namely, the formation energy of oxides (ΔE_f) obtained as the fitted elemental-phase reference energy (FERE),³⁵ the computed band-gaps obtained by GGA+U (E_g^GGA+U) or GW⁴⁰ calculations, the O 2p band center (E_O2p), and the differences in electronegativity between the cations and O, for 45 different oxides in six crystal structures. Their prediction error was reported as ∼0.20 eV as the mean absolute error (MAE). To evaluate whether this model is well-fitted for the different crystal structures upon fixing the substituted chemical elements and composition ratio as listed in Table 1, we employed the OLSR-3-descriptors model with the descriptors x₁–x₃ listed in Table 2. For the RR-11 descriptor, the average difference in electronegativity between the cations and O is not used because it does not vary among the different types of crystal structures. The 16 crystal structures of the reoptimized-ZrO₂ in Table 1 were used to construct the prediction model.

Fig. 3(a) shows the relationship between the predicted and computed E_v obtained by using the OLSR-3-descriptors model. The distribution of data is largely scattered from the diagonal line that indicates equality between the predicted and computed E_v. As arranged in Table 3, the CV-score and RMSE of the test set are 0.76 and 0.51 eV, respectively, which are much larger errors than the MAE of 0.20 eV for the 45 oxides.²³ The large difference between the CV-score and RMSE of the test set also indicates that the prediction model has a large dependence on the type of training set. Therefore, it implies that the prediction model using the OLSR-3-descriptors is not good enough to predict the E_v of the various crystal structures of ZrO₂.

Fig. 3 Relationship between the computed and predicted E_v by using regression analyses with the 16 crystal structures of reoptimized-ZrO₂ with (a) 3-OLSR-descriptors and (b) 11-RR-descriptors. One snapshot is drawn from thirty repeated random samples for training and test data. The RMSE of test data averaged from thirty samplings from the prediction models using (a) the 3-OLSR-descriptors and (b) the 11-RR-descriptors are 0.51 and 0.16 eV, respectively.

Table 3 Errors between the predicted and computed E_v of two regression analyses using the 16 crystal structures of reoptimized-ZrO₂. The prediction error values are obtained as RMSE in eV. The parentheses denote the standard deviation of values averaged over thirty readings

Prediction Model	Data divided into training (75%) and test (25%) set			Whole data
	RMSE of training set	CV-score	RMSE of test set	RMSE of training set	CV-score
OLSR-3-descriptors	0.27 (0.05)	0.76 (0.36)	0.51 (0.22)	0.31	0.60 (0.19)
RR-11-descriptors	0.04 (0.04)	0.17 (0.06)	0.16 (0.16)	0.02 (0.01)	0.14 (0.07)

---

To improve the prediction accuracy, we constructed RR-11-descriptors by adding more DFT-unit-cell descriptors on top of the three descriptors; these DFT descriptors are intuitively thought to be related to the forming and breaking of chemical bonds. We employed the RR model, which can control the coefficients of various descriptors through a penalty term and thereby minimize overfitting. Fig. 3(b) shows the relationship between the predicted and computed E_v obtained by using the RR-11-descriptors model. The scattering of points from the diagonal line is significantly decreased. The CV-score and RMSE of the test sets are 0.17 and 0.16 eV, respectively, which are much smaller than those obtained from the OLSR-3-descriptors model. The similarity in the two values suggests the absence of any serious overfitting.

Furthermore, we applied the RR-11-descriptor model with a training set consisting of all 16 crystal structures of reoptimized-ZrO₂. In the results, the CV-score slightly decreased to 0.14 eV. This improvement can be attributed to the increased number of training data points. This value is also only 0.02 eV different from that of the RMSE of the test set of the prediction model with a 75% training set, which suggests the absence of any serious overfitting. Finally, we identified the coefficients for a linear fitting function applied as the fitness value in the evolutionary algorithm by using the RR-11-descriptors model with all 16 crystal structures of the reoptimized-ZrO₂.

3.2. Crystal structures with low E_v

As Fig. 1 shows, we generated 462 crystal structures through combinations of the regression analysis for the E_v and evolutionary algorithm, and screened out 71 crystal structures based on the criterion |predicted E_v − 5| < 0.5 eV. After crystal relaxation once again under tighter conditions, we obtained 14 remaining crystal structures with the predicted and computed E_v, excluding those that were metallic or could not satisfy the criteria for structural optimization. The detailed properties of the 14 crystal structures are listed in Table 4. Henceforth, this set of crystal structures is referred to as “generated-ZrO₂”, and the name of each such structure is prefixed with “Gen-” (“generated or constructed”); these structures are numbered in the order of the computed E_v. These crystal structures have very low symmetry; most of them are P1. To identify similar crystal structures with higher symmetry, a looser tolerance for finding the space group was also applied.

Table 4 Computed properties of the 14 crystal structures of generated-ZrO₂

Name of constructed crystal structure	Space group type^a	E − E(ground-state) (eV per atom)	Eg^GGA+U (eV)	Predicted Ev (eV)	Computed Ev (eV)^b
a The space group in parentheses is identified using looser tolerance in SPGLIB in PHONOPY code^36,37 to determine similar higher-symmetry crystal structures.b Values in the parentheses are the minimum and the maximum values in sequence, which depend on the kinds of O sites.
Gen-01	P1 (Pnma)	0.113	2.77	5.11	5.11 (4.94, 5.30)
Gen-02	Pnma (Pnma)	0.110	2.82	5.18	5.21 (5.11, 5.30)
Gen-03	P1 (Pnma)	0.110	2.83	5.18	5.21 (5.11, 5.30)
Gen-04	P1̄ (Pnma)	0.110	2.83	5.18	5.21 (5.11, 5.30)
Gen-05	P1 (Pnma)	0.110	2.83	5.18	5.21 (5.11, 5.31)
Gen-06	P212121 (Pnma)	0.110	2.83	5.18	5.21 (5.11, 5.30)
Gen-07	Pca21 (P63mc)	0.223	2.86	5.80	5.69 (5.56, 5.83)
Gen-08	P1 (P42/nmc)	0.024	3.79	5.80	5.79 (5.79, 5.79)
Gen-09	P1 (P42/nmc)	0.023	3.97	5.75	5.85 (5.85, 5.85)
Gen-10	P1 (P21/m)	0.021	4.20	5.85	6.03 (5.90, 6.16)
Gen-11	P1 (P21/m)	0.021	4.20	5.86	6.03 (5.90, 6.17)
Gen-12	P1̄ (P21/m)	0.021	4.18	5.86	6.03 (5.90, 6.17)
Gen-13	P1 (Fdd2)	0.030	3.44	6.00	6.20 (6.07, 6.33)
Gen-14	P1 (P21/m)	0.113	3.70	6.53	6.39 (6.33, 6.53)

---

Fig. 4 shows the relationship between the computed and predicted E_v of the crystal structures of reoptimized-ZrO₂ and generated-ZrO₂, which are used as the training and test data, respectively. The E_v of the ground-state P2₁/c structure is emphasized with a star-shaped mark for easier viewing. Note that the distribution of the predicted E_v for the 14 test-data crystal structures is in the range 5.11–6.53 eV, which is wider than the screening criteria. This is because some crystal structures are relaxed into more stable structures, with changes in the predicted E_v.

Fig. 4 Relationship between the computed and predicted E_v of the training- and test-data crystal structures. The crystal structures for the training and test data are from reoptimized-ZrO₂ and generated-ZrO₂, respectively. The predicted E_v is predicted by the 11-RR-descriptors model. The RMSE of the 14 test-data crystal structures is 0.11 eV. Considering that the three-fold CV-score for the 16 training data is almost the same as 0.14 eV, overfitting is avoided.

The distribution of the test data is close to the diagonal line. The prediction error is only 0.11 eV, which is almost the same as the CV-score of the RR-11-descriptors model with the 16 training data (0.14 eV). This means that the E_v of the 14 crystal structures of generated-ZrO₂ have been successfully predicted using only DFT-unit-cell descriptors. Note that the distribution of the 16 training data points is very close to the diagonal line in Fig. 3(b), because their predicted E_v are obtained as the training data, and not as the test data.

Fig. 5 shows Pearson's linear correlation coefficient (r_p) for each DFT-unit-cell descriptor with respect to the computed E_v. The absolute value of r_p is larger than 0.7, which indicates a strong correlation between the two properties for a given data point.⁴¹ For the 16 crystal structures of reoptimized-ZrO₂, the |r_p| for the volume per atom, average bond lengths (r_Zr–O), and O coordination number (CN_O) are 0.82, 0.85, and 0.73 (high), respectively, while the |r_p| for the formation energy (ΔE_f) and E_g^GGA+U are 0.22 and 0.42 (low), respectively. This observation differs from that of Deml's model²³ for 45 various oxides, which showed strong correlations between the FERE and E_v, and between the E_g^GGA+U and E_v, with |r_p| values of 0.84 and 0.89, respectively. The different accuracies of the prediction models between this study and ref. 23 can arise from several factors. Firstly, the 45 training-data oxides used in ref. 23 included six types of crystal structures with high symmetries, namely, anti-fluorite, fluorite, rocksalt, rutile, perovskite, and spinel. Most of these oxides are energetically stable for various combinations of constituent elements. However, the crystal structures of ZrO₂ in this study all have different types of metastable lower-symmetry crystal structures. The contribution of each descriptor in the regression analysis may differ. Secondly, the range of training data can also affect the accuracy of prediction model. The E_v of the 45 training-data structures in ref. 23 are widely distributed between 2.6 and 7.0 eV, whereas the E_v of ZrO₂ in this study are distributed between only 5.1 and 6.7 eV.

Fig. 5 Pearson's linear correlation coefficients (r_p) between the computed E_v and the 11 DFT-unit-cell descriptors for the 16 training data points (reoptimized-ZrO₂) and the 30 total data points (the crystal structures of reoptimized-ZrO₂ and generated ZrO₂).

The r_p changes for all 30 crystal structures of the reoptimized-ZrO₂ and generated-ZrO₂. The values of |r_p| for the volume per atom and r_Zr–O are similarly large as 0.88 and 0.82 respectively, while |r_p| for CN_O decreases to 0.49. The |r_p| for ΔE_f and E_g^GGA+U changed to 0.51 and 0.69, respectively. When the RR-11-descriptors model is applied to all 30 crystal structures, the prediction accuracy is improved, with the RMSE of the test data of 0.08 eV compared with that of the same prediction model (0.16 eV) with the 16 training data crystal structures. Because the descriptors also maintain their correlations with each other in the regression analysis, the important unit-cell descriptors for the prediction of E_v can be changed according to the target materials and training data. Therefore, we suggest employing more types of unit-cell descriptors but with regression analysis using a penalty term; this is a useful way to construct the prediction model of the E_v.

We identify many crystal structures having E_v lower than that of the ground-state P2₁/c structure. As shown in Table 1 and Fig. 4, eight crystal structures, including the P4₂/nmc and Fm3̄m structures that are well-known high-temperature structures of ZrO₂ and that can be stabilized by extrinsic doping,⁵ reveal lower E_v than that of the ground-state P2₁/c structure. In particular, the Pna2₁ structure, which is obtained from the structural data of SiO₂, shows the lowest E_v of 5.14 eV, which is almost 1 eV less than that of the ground-state P2₁/c structure. This implies that exchanging the substituting elements in the crystal structures of other oxides can be an efficient way to reconsider a functional property. Additionally, we obtain 14 crystal structures based on the evolutionary algorithm that have lower E_v than that of the ground-state P2₁/c structure. This suggests that the discovery of a crystal structure with a special functional property can be accelerated by a combination of the machine-learning technique and the direct construction of crystal structures based on the evolutionary algorithm.

Meanwhile, we reveal why the neutral oxygen vacancy is mainly employed in the calculations of the E_v for the prediction model. The V_O in the crystal structures of ZrO₂ has been widely investigated. The P2₁/c and P4₂/nmc structures are known to form a deep level inside the E_g based on the GGA^42–44 and HSE06 (ref. 45 and 46) calculations. Additionally, the E_v of the in the dilute limit is physically well defined regardless of the doping level of the cation sites, and quickly converged with respect to the supercell size.

However, it is also worth investigating the relationship between and considering the cases of extrinsic doping by aliovalent ions. Here, for , we fixed the E_Fermi at the center of the E_g because the crystal structures of reoptimized-ZrO₂ are all insulators (E_g^GGA+U > ∼3 eV). Among the 16 crystal structures of the reoptimized-ZrO₂ in Table 1, the Pnma and the P6₃mc structures were excluded from the computation of because crystal structures with a are considered too unstable for convergence. Fig. 6 shows the relationship between the computed and of the 14 crystal structures of the reoptimized-ZrO₂. It is noticeable that the and have a very strong correlation with the r_p of 0.86. This suggests that the crystal structures with low also have low , which implies that extrinsic doping with aliovalent ions may also more easily generate V_O for these crystal structures (numerical values are summarized in Table S3 in ESI†).

Fig. 6 Relationship of the computed and of the 14 crystal structures of reoptimized-ZrO₂. The r_p between the and is 0.86, which implies a very strong correlation between the two properties.

Before proceeding, it is worthwhile to discuss the comparison of the computed values (E_g and E_v shown in Fig. S1 in ESI†) between the GGA+U and GGA methods. It is noticeable that the E_g obtained by the GGA+U are larger than those by the GGA, which implies that the GGA+U is useful to reduce the underestimation of the E_g. For the computed E_v, it is found that the E_v obtained by the GGA are slightly smaller or almost the same as those obtained by the GGA+U. This suggests that the GGA+U does not lose its reliability compared with the GGA, despite including the empirical U parameter.

3.3. Crystal structures with low E_m

We summarize the E_m values of the Fm3̄m and the P4₂/nmc structures obtained herein and in former computational studies^42,47–50 in Table 5. Former computational studies have mainly used the GGA, rather than the GGA+U. To compare the results, we also employ the GGA with the PBE form.³³ The E_m values obtained herein using the GGA are almost the same as those mentioned in previous computational reports. The E_m values obtained using the GGA+U are larger than those obtained using the GGA for both the Fm3̄m and P4₂/nmc structures.

Table 5 E_m values of the Fm3̄m and the P4₂/nmc structures

Crystal structure	Em (eV)^a		Reference	Method	GGA parameterization
	From GGA	From GGA+U
a Values outside and in the parentheses are obtained using and , respectively.b (110) direction in a tetragonal cell.c (001) direction in a tetragonal cell.d Ultrasoft pseudopotentials.^66
Fm3̄m	0.26 (0.54)	0.39 (1.80)	This study	PAW (VASP code)	Perdew et al. (PBE)^33
	0.25 (0.45)		47	Atomic orbital (SIESTA code)^63	Perdew et al. (PBE)^33
	0.24		48	PAW (VASP code)	Perdew et al. (PBE)^33
	0.26 (0.54)		49	PAW (VASP code)	Perdew et al. (PBE)^33
	0.28		50	Plane-wave with US-PP^d (Quantum-espresso code)^64	Perdew et al. (PBE)^33
P42/nmc	0.26 (1.44)^b, 0.63 (1.57)^c	0.34 (1.74)^b, 0.64 (1.66)^c	This study	PAW (VASP code)	Perdew et al. (PBE)^33
	0.22 (1.35)^b, 0.61 (1.43)^c		42	PAW (VASP code) with cutoff energy of 250 eV	Perdew et al.^65

---

We also compare the E_m values of and . Table 5 shows that the E_m values of are much larger than those of , independent of the crystal structure and the exchange-correlation functional. Larger E_m values for than for have been reported in other oxides like MgO⁵¹ and ZnO.⁵² When a V_O jumps into its nearest-neighboring O site, the O atom at the nearest-neighboring site migrates in the opposite direction. For , no occupied levels are present inside E_g. On the other hand, for , the two electrons occupying the deep level have transition into the conduction band when the moving O atom is in the transition state (center of the hopping path),⁴⁷ which may need more energy for V_O hopping. Considering the abovementioned reasons, we employed for the E_m in this study.

Hereafter, we will refer to the E_m obtained by the GGA+U calculations.

In the case of the initial-state structure for the CI-NEB method, we employed the relaxed supercell possessing a having the lowest energy among the different O sites. For the 14 crystal structures of reoptimized-ZrO₂ that have well-converged structures with a , we considered various final-state structures that have the V_O at the nearest-neighboring O site with the cutoff distance of 3.5 Å for the CI-NEB method (see Fig. S2 in ESI†). Among the E_m values of different migration paths, the lowest value was extracted; we then checked whether the V_O could migrate to the original site in the nearest-neighboring supercell via the one type of migration path with the lowest E_m, because the migration of the V_O through the bulk system is practically meaningful. In 12 of the 14 crystal structures, V_O could migrate to the original site in the nearest-neighboring supercell via the migration path with the lowest E_m. However, for the P4/n and R3̄ structures, a V_O needs different paths with different E_m values to migrate to the original site in the nearest-neighboring supercell. We define the E_m for these structures as the minimum barrier energy required for migration from the original site to that in the nearest-neighboring supercell with several kinds of CI-NEB calculations. The energy diagrams for the migration of a V_O are shown in Fig. S3 and Table S4 in ESI.†

Fig. 7 and 8 show the relationships between the computed E_v and E_m values, and between the relative ΔE_f (compared with that of the ground-state P2₁/c structure) and E_m values of various ZrO₂ crystal structures, respectively. Firstly, we compare the E_m of the 14 crystal structures of the reoptimized-ZrO₂. The E_m values (0.39 and 0.34 eV, respectively) of the Fm3̄m and P4₂/nmc structures are lower than that (0.61 eV) of the ground-state P2₁/c structure. Ten crystal structures, including the Fm3̄m and the P4₂/nmc structures, show lower E_m than that of the ground-state P2₁/c structure.

Fig. 7 Relationship between the computed E_v [(a) and (b) , respectively] and E_m of various crystal structures of ZrO₂. The r_p between the E_m and and between the E_m and of the 14 crystal structures from reoptimized-ZrO₂ are 0.64 and 0.62, respectively, which implies that the correlation in a reasonable level is found. The crystal structures with asterisks (*) need several different types of migration for a V_O to move to the same site in the nearest-neighboring computational cell (see Fig. S3 in ESI†). On the basis of this information, the E_m of three crystal structures from generated-ZrO₂ are obtained (red × marks).

Fig. 8 Relationship between the E_m and relative ΔE_f (compared with that of the ground-state P2₁/c structure) of the 14 crystal structures of reoptimized-ZrO₂ and three crystal structures of generated-ZrO₂ which have E_m values.

Among these ten, eight crystal structures have only 0.05 eV per atom higher ΔE_f than the ground-state P2₁/c structure (Table 1 and Fig. 8). This suggests that high σ_O can be realized for these structures with stabilizing methods such as extrinsic doping and the formation of epitaxial heterostructures. The P4₂/nmc and Fm3̄m structures are confirmed to have both relatively low E_m and ΔE_f. The Pbcn structure, which also shows relatively low ΔE_f, is reportedly produced by biaxial tensile strain⁵³ on fluorite-structured CeO₂, which forms the (100) interface. This structure, obtained by 7% tensile strain, is also reported to show a higher O diffusivity than that of YSZ based on first-principles molecular dynamics calculations.⁵⁴ The I4₁/amd, Pbca, Pca2₁, C2/c, and P2₁/m structures also have the potential for stabilization, although they have somewhat larger E_m than the P4₂/nmc and Fm3̄m structures do. These crystal structures have not been investigated as oxygen-ion conductors to the best of the authors' knowledge; however, some of them have been reported as metastable structures. The Pbca structure was experimentally synthesized at 600 °C and 6 GPa.⁵⁵ This structure was also confirmed by first-principles calculations to be stabilized by hydrostatic pressure over 4.1 GPa.⁵⁶ Mg-doped tetragonal ZrO₂ cooled to ∼30 K reportedly transforms to the orthorhombic structure with the space group of Pbc2₁ (Pca2₁).⁵⁷

The E_m values for the 14 crystal structures of the reoptimized-ZrO₂ have similarly large r_p of 0.64 and 0.62 with the computed and , respectively. The E_v is the energy required for breaking all the chemical bonds to separate an O atom from an oxide, and the E_m is the energy required for breaking some of the chemical bonds locally in an oxide, as shown in Fig. 2. Therefore, this correlation may be reasonable in terms of chemical bond cleavage. Correlations between E_v and E_m have also been reported for other oxides, like perovskite oxides⁵⁸ with various kinds of constituent elements and a multinary oxide of Ba_1−xSr_xCo_1−yFe_yO_3−δ.⁵⁹ On the other hand, the relative ΔE_f and E_m do not show any correlation, having an |r_p| of only 0.10, which implies that the relative ΔE_f is not a good single descriptor for the E_m. To realize crystal structures with low values for both relative ΔE_f and E_m, special challenges in processing may be necessary.

The above information suggests that we are likely to identify a low-E_m crystal structure when we initially investigate crystal structures with low E_v. Therefore, for E_m calculations, we choose the four crystal structures (Gen-01, Gen-02, Gen-07, and Gen-08) of the generated-ZrO₂ that possess the lowest E_v values in ascending order. The crystal structures between Gen-03 and Gen-06 are not chosen because we expect these structures to have properties similar to Gen-02 (see Fig. S4 in ESI†). The crystal structures from Gen-09 to Gen-14 are also not considered for the E_m because they may have relatively high E_m values.

The E_m values of these crystal structures from the generated-ZrO₂ are also shown in Fig. 7 and 8. The E_m values of Gen-01, Gen-02, and Gen-08 are 0.42, 0.39, and 0.28 eV, respectively. Therefore, we succeed in constructing crystal structures with lower E_m as well as lower E_v than those of the ground-state P2₁/c structure.

It is worth analyzing the Gen-08 structure further because it is very similar in crystal structure and energetic properties to P4₂/nmc, despite its space group being P1. Fig. 9 shows the crystal structures of the P4₂/nmc and Gen-08 structures. We confirm that Gen-08 and P4₂/nmc do not become the same by structural optimizations with denser k-space sampling. However, the space group defined by using looser tolerance and structural properties, such as the radial distribution function and the bond lengths (see Fig. S4 in ESI†), imply that these crystal structures are quite similar; the Gen-08 structure is slightly distorted in terms of its bond lengths and lattice parameters (see Table S5 in ESI†). In the macroscopic view, these crystal structures may be difficult to distinguish from each other. For the construction of the Gen-08 structure, we only prepared the prediction model of E_v by the regression analysis and evolutionary algorithm that requires the initial structure-related input of only the number of constituting atoms of Zr and O in empty space. Therefore, the construction of Gen-08, which has a similar crystal structure to P4₂/nmc, implies a successful result of a combination of evolutionary algorithm and regression analysis.

Fig. 9 Unit-cells of (a) P4₂/nmc and (b) Gen-08 structures, and supercells of (c) P4₂/nmc and (d) Gen-08 structures of ZrO₂. Numerical values in (c) and (d) are O–Zr bond lengths (in Å). Green and light gray spheres denote Zr and O atoms, respectively.

The Gen-07 structure, possessing an energy of 0.2 eV per atom higher than that of the ground-state P2₁/c structure, is too unstable to contain a . This suggests that the energetic stability and the low E_v of the crystal structures should be considered in choosing candidates for low E_m. The Gen-01–Gen-07 crystal structures have energies more than 0.1 eV per atom higher than that of ground-state P2₁/c structure, despite having E_v lower than 5.7 eV. On the other hand, the crystal structures Gen-08–Gen-13 have less than 0.03 eV per atom higher energy than the ground-state P2₁/c structure, although they have higher E_v values. The latter group may have higher stabilization potential than the former.

Before concluding, we note our exclusion of dynamical stabilities based on the phonon calculations with harmonic approximations^36,37 from this study. The P4₂/nmc and Fm3̄m structures are known to be stabilized at temperatures exceeding ∼1173 °C and ∼2370 °C, respectively. However, as Fig. S5 in ESI† shows, the P4₂/nmc structure includes no imaginary phonon frequencies, whereas the Fm3̄m structure includes imaginary phonon frequencies despite its capacity for high-temperature stabilization, which is consistent with the phonon calculations by Wang et al.⁶⁰ Therefore, anharmonic phonon effects,^61,62 which are particularly significant at high temperatures, must be considered for these materials.

4. Conclusion

We have constructed a prediction model for the E_v values of various crystal structures of ZrO₂ using a linearized RR with 11 unit-cell descriptors. The predicted E_v was used as the fitness value for crystal structure constructions based on the evolutionary algorithm in order to mainly construct crystal structures having E_v lower than that of the ground-state P2₁/c structure. Our prediction model guarantees prediction errors of less than ∼0.15 eV for various crystal structures of ZrO₂.

We also obtained the E_m values of various crystal structures of ZrO₂, and found that this property is correlated with E_v. On the basis of this correlation, we calculated the E_m values of several newly constructed crystal structures having low E_v, and confirmed that these E_m values are also lower than those of the ground-state P2₁/c structure and lower than or similar to those of the P4₂/nmc and Fm3̄m structures.

We have successfully discovered various crystal structures with low E_v and E_m values for ZrO₂ by the direct construction of crystal structures based on a combination of evolutionary algorithm and regression analysis and by reconsidering the crystal structures present in other oxides. Our DFT-unit-cell descriptors are constructed on the basis of the properties related to the chemical bonds between cations and O. Therefore, the good accuracy of the prediction model for E_v and the correlation between the E_v and E_m may be reasonable because these properties are all related to the chemical bond cleavage. We expect that the method employed in this study can be applied to determine unknown low-E_v and -E_m crystal structures of other material systems, which have not been investigated closely. In addition, the strategy described herein can be also applied to search for crystal structures with other functional properties.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Dr N. Nagasako in TCRDL for fruitful discussions.

References

1. P. Knauth and H. L. Tuller, J. Am. Ceram. Soc., 2002, 85, 1654–1680.
2. S. J. Skinner and J. A. Kilner, Mater. Today, 2003, 6, 30–37.
3. S. C. Singhal and K. Kendall, High-temperature solid oxide fuel cells: fundamentals, design and applications, Elsevier, 2003.
4. P. J. Gellings and H. Bouwmeester, Handbook of solid state electrochemistry, CRC Press, 1997.
5. M. Mogensen, N. M. Sammes and G. A. Tompsett, Solid State Ionics, 2000, 129, 63–94.
6. N. Sammes, G. Tompsett, H. Näfe and F. Aldinger, J. Eur. Ceram. Soc., 1999, 19, 1801–1826.
7. K. Shitara, T. Moriasa, A. Sumitani, A. Seko, H. Hayashi, Y. Koyama, R. Huang, D. Han, H. Moriwake and I. Tanaka, Chem. Mater., 2017, 29, 3763–3768.
8. K. Huang and J. B. Goodenough, J. Alloys Compd., 2000, 303, 454–464.
9. T. Ishihara, H. Matsuda and Y. Takita, J. Am. Chem. Soc., 1994, 116, 3801–3803.
10. X. Kuang, M. A. Green, H. Niu, P. Zajdel, C. Dickinson, J. B. Claridge, L. Jantsky and M. J. Rosseinsky, Nat. Mater., 2008, 7, 498–504.
11. A. R. Oganov and C. W. Glass, J. Chem. Phys., 2006, 124, 244704.
12. A. R. Oganov, C. W. Glass and S. Ono, Earth Planet. Sci. Lett., 2006, 241, 95–103.
13. A. R. Oganov, J. Chen, C. Gatti, Y. Ma, Y. Ma, C. W. Glass, Z. Liu, T. Yu, O. O. Kurakevych and V. L. Solozhenko, Nature, 2009, 457, 863–867.
14. Y. Shen, A. R. Oganov, G. Qian, J. Zhang, H. Dong, Q. Zhu and Z. Zhou, Sci. Rep., 2015, 5, 14204.
15. A. R. Oganov and A. O. Lyakhov, J. Superhard Mater., 2010, 32, 143–147.
16. Q. Zeng, A. R. Oganov, A. O. Lyakhov, C. Xie, X. Zhang, J. Zhang, Q. Zhu, B. Wei, I. Grigorenko and L. Zhang, Acta Crystallogr., Sect. C: Struct. Chem., 2014, 70, 76–84.
17. C. G. Van de Walle and J. Neugebauer, J. Appl. Phys., 2004, 95, 3851–3879.
18. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
19. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
20. A. Seko, T. Maekawa, K. Tsuda and I. Tanaka, Phys. Rev. B, 2014, 89, 054303.
21. J. Lee, A. Seko, K. Shitara, K. Nakayama and I. Tanaka, Phys. Rev. B, 2016, 93, 115104.
22. A. Seko, A. Takahashi and I. Tanaka, Phys. Rev. B, 2014, 90, 024101.
23. A. M. Deml, A. M. Holder, R. P. O'Hayre, C. B. Musgrave and V. Stevanović, J. Phys. Chem. Lett., 2015, 6, 1948–1953.
24. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, Appl. Phys. Lett. Mater., 2013, 1, 011002.
25. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
26. A. E. Hoerl and R. W. Kennard, Technometrics, 1970, 12, 55–67.
27. A. R. Oganov, A. O. Lyakhov and M. Valle, Acc. Chem. Res., 2011, 44, 227–237.
28. A. O. Lyakhov, A. R. Oganov, H. T. Stokes and Q. Zhu, Comput. Phys. Commun., 2013, 184, 1172–1182.
29. P. E. Blöchl, Phys. Rev. B, 1994, 50, 17953–17979.
30. G. Kresse and D. Joubert, Phys. Rev. B, 1999, 59, 1758–1775.
31. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
32. G. Kresse and J. Furthmüller, Phys. Rev. B, 1996, 54, 11169–11186.
33. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
34. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B, 1998, 57, 1505–1509.
35. V. Stevanović, S. Lany, X. Zhang and A. Zunger, Phys. Rev. B, 2012, 85, 115104.
36. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B, 2008, 78, 134106.
37. A. Togo, L. Chaput, I. Tanaka and G. Hug, Phys. Rev. B, 2010, 81, 174301.
38. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
39. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
40. M. Shishkin and G. Kresse, Phys. Rev. B, 2007, 75, 235102.
41. J. D. Evans, Straightforward statistics for the behavioral sciences, Brooks/Cole, 1996.
42. A. Eichler, Phys. Rev. B, 2001, 64, 174103.
43. J. Zheng, G. Ceder, T. Maxisch, W. Chim and W. Choi, Phys. Rev. B, 2007, 75, 104112.
44. A. S. Foster, V. Sulimov, F. L. Gejo, A. Shluger and R. M. Nieminen, Phys. Rev. B, 2001, 64, 224108.
45. J.-H. Hur, S. Park and U.-I. Chung, J. Appl. Phys., 2012, 112, 113719.
46. W.-J. Yin, S.-H. Wei, M. M. Al-Jassim and Y. Yan, Appl. Phys. Lett., 2011, 99, 142109.
47. S. Kasamatsu, T. Tada and S. Watanabe, Appl. Phys. Express, 2009, 2, 061402.
48. J. Hirschfeld and H. Lustfeld, Phys. Rev. B, 2011, 84, 224308.
49. B. Liu, H. Xiao, Y. Zhang, D. S. Aidhy and W. J. Weber, Comput. Mater. Sci., 2014, 92, 22–27.
50. O. I. Malyi, P. Wu, V. V. Kulish, K. Bai and Z. Chen, Solid State Ionics, 2012, 212, 117–122.
51. J. Mulroue and D. Duffy, Proc. R. Soc. A, 2011, 467, 2054–2065.
52. P. Erhart and K. Albe, Phys. Rev. B, 2006, 73, 115207.
53. D. S. Aidhy, B. Liu, Y. Zhang and W. J. Weber, J. Phys. Chem. C, 2014, 118, 30139–30144.
54. M. Oka, H. Kamisaka, T. Fukumura and T. Hasegawa, Phys. Chem. Chem. Phys., 2015, 17, 29057–29063.
55. O. Ohtaka, T. Yamanaka, S. Kume, N. Hara, H. Asano and F. Izumi, Proc. Jpn. Acad., Ser. B, 1990, 66, 193–196.
56. J. Dewhurst and J. Lowther, Phys. Rev. B, 1998, 57, 741.
57. E. H. Kisi, C. J. Howard and R. J. Hill, J. Am. Ceram. Soc., 1989, 72, 1757–1760.
58. T. T. Mayeshiba and D. D. Morgan, Solid State Ionics, 2016, 296, 71–77.
59. R. Merkle, Y. A. Mastrikov, E. A. Kotomin, M. M. Kuklja and J. Maier, J. Electrochem. Soc., 2011, 159, B219–B226.
60. G. Wang, G. Luo, Y. L. Soo, R. F. Sabirianov, H.-J. Lin, W.-N. Mei, F. Namavar and C. L. Cheung, Phys. Chem. Chem. Phys., 2011, 13, 19517–19525.
61. O. Hellman, I. A. Abrikosov and S. I. Simak, Phys. Rev. B, 2011, 84, 180301.
62. P. Souvatzis, O. Eriksson, M. I. Katsnelson and S. P. Rudin, Phys. Rev. Lett., 2008, 100, 095901.
63. E. Artacho, E. Anglada, O. Diéguez, J. D. Gale, A. García, J. Junquera, R. M. Martin, P. Ordejón, J. M. Pruneda and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2008, 20, 064208.
64. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni and I. Dabo, J. Phys.: Condens. Matter, 2009, 21, 395502.
65. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B, 1992, 46, 6671–6687.
66. D. Vanderbilt, Phys. Rev. B, 1990, 41, 7892–7895.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8ra02958j

---