Prediction of perovskite-related structures in ACuO_3-x (A $=$ Ca, Sr, Ba, Sc, Y, La) using density functional theory and Bayesian optimization

In general, it is impossible to estimate the configurational energy for a complete set of nonequivalent structures by DFT calculation owing to the considerable number of nonequivalent structures. Therefore, it has been a practical way to employ a machine learning model of the configurational energy estimated from a set of DFT calculations for sampled structures. The CE method Sanchez et al. (1984); de Fontaine (1994); Ducastelle (1994) has been a widely used approach, and gives a reliable model of configurational energy; hence, it has enabled us to accurately predict the ground-state structures and finite-temperature thermodynamics in many multicomponent systems.

In a binary system, the configurational energy $E$ is described in a linear form using the pseudospin configurational variable $\sigma_{i}$ for the respective lattice site $i$ as

where $\varphi_{\alpha}$ denotes the correlation function of cluster $\alpha$. Coefficients $V_{\alpha}$ are called the effective cluster interactions, which are estimated from a set of DFT calculations for sampled configurations using general linear regression methods such as conventional linear regression and the least absolute shrinkage and selection operator (LASSO) Nelson et al. (2013); Hastie et al. (2009).

We prepared a training dataset that contains the formation energies of all nonequivalent structures with up to $n=4$ estimated by DFT calculation performed using the plane-wave-basis projector augmented wave (PAW) method Blöchl (1994) within the Perdew–Burke–Ernzerhof generalized gradient approximation (GGA) exchange-correlation functional Perdew et al. (1996), as implemented in the vasp code Kresse and Hafner (1993); Kresse and Furthmüller (1996); Kresse and Joubert (1999). The cutoff energy was set to 400 eV. The total energies converged to less than 10^-3 meV/supercell. The atomic positions and lattice constants of the nonequivalent structure were optimized until the residual forces were less than 10^-2 eV/Å to consider the relaxation effect of the ideal cubic perovskite lattice. We assumed the ferromagnetic configuration for copper atoms in the DFT calculation for simplicity, although ACuO_3-x is expected to show an antiferromagnetic configuration as shown experimentally in Ca_1-xSr_xCuO₂ Vaknin et al. (1989). We adopted the $+U$ approach to improve the description of the ground state of correlated systems Dudarev et al. (1998). The value of $U$ for copper was set to 4 eV taken from the literature Wang et al. (2006).

Figure 2 shows the formation energies of nonequivalent oxygen-deficient structures with up to $n=4$ in ACuO_3-x systems. They are measured from the energy of the cubic perovskite ACuO₃ and the lowest energy among the structures in ACuO₂. Since the DFT calculation fails to finish successfully in some structures, we exclude the structures from the dataset. The sizes of the resultant training datasets are 485, 473, 480, 457, 468, and 480 for BaCuO_3-x, CaCuO_3-x, SrCuO_3-x, ScCuO_3-x, YCuO_3-x, and LaCuO_3-x, respectively.

Then, we construct various CE models for each system using the training dataset. We employ spin values of $+1$ and $-1$ for the oxygen atom and vacancy, respectively, which are required to define the correlation function. Table 2 shows the fitting errors of the four CE models, which indicates that all of the CE models have significant errors. They are much more significant than the order of the difference of the formation energies between structures. Figure 3 also shows the distribution of the formation energies of nonequivalent oxygen-deficient structures estimated using a CE model for SrCuO_3-x. The formation energies deviate from those estimated by DFT calculations. Therefore, it is impossible to find the true ground states using only the energy predicted using the CE model.

The significant error in the CE model arises from the local geometry optimization of the oxygen-deficient structures in the DFT calculation. Although the local geometry optimization is essential for predicting oxygen-deficient structures, the CE method describes the energy of the optimized structure using a linear polynomial model of the correlation functions defined for the ideal lattice of cubic perovskite. A general procedure of the CE method has been successfully applied to most of the substitutional alloys, while it does not work well in systems involving the significant atomic relaxation from the ideal lattice sites in a straightforward manner.

To find stable oxygen-deficient perovskite structures, we employ Bayesian optimization Jones (2001) based on GP Rasmussen and Williams (2006) specified by its mean function and covariance function instead of a linear CE model. Representing a structure by vector $\bm{d}$, a radial basis function covariance or kernel between structures $\bm{d}_{i}$ and $\bm{d}_{j}$ for noise-free observation is given by

where $l$ and $\sigma_{f}^{2}$ are tuning parameters controlling the length scales for structure $\bm{d}$ and the observation, respectively. Here, we represent oxygen-deficient structures by the correlation functions of point clusters and 11 pair clusters used in the CE models shown in Fig. 6 (c). As described above, the discrepancy between the DFT calculation and the CE model is too large to predict the stable oxygen-deficient perovskite structures using the CE model itself. On the other hand, the CE models have a strong correlation with the DFT calculation, as shown in Fig. 3. This implies that a set of correlation functions is useful for representing a structure in GP-based Bayesian optimization.

Describing the formation energies of structures in the training dataset as $\bm{y}$, the mean function $\mu$ at structure $\bm{d}_{*}$ and the variance function $\sigma_{*}^{2}$ are given as

and

respectively, where ${\bm{k}}_{*}=\left[k(\bm{d}_{*},\bm{d}_{1}),\cdots,k(\bm{d}_{*},\bm{d}_{N})\right]^{\top}$ is the vector of kernel functions between structure $\bm{d}_{*}$ and $N$ structures in the training dataset. A symmetric kernel matrix $\bm{K}$ is composed of kernel functions for all pair arrangements of the training data, expressed as

Our procedure of Bayesian optimization is as follows. First, a GP model is developed from an initial set of formation energies. The initial set is composed of the formation energies of all nonequivalent structures with up to $n=4$ and 100 nonequivalent structures that are randomly selected from the pool of 773,613 nonequivalent structures with up to $n=8$. The model is then iteratively updated by (i) sampling the structure for which the formation energy is expected to be the lowest among the remaining candidate structures for each composition and (ii) updating the GP model using the training dataset updated by the structures observed in step (i). These steps are repeated until the convex hull of the formation energy converges. We exclude some structures from the dataset, i.e., structures for which DFT calculations fail to converge.

There are several well-known procedures to sample the structure that is expected to show the lowest formation energy in step (i). Here, we adopt a simple procedure of the probability of improvement (PI) strategy Jones (2001). In the minimization problem of formation energy, the PI involves sampling the structure at which the probability that the formation energy is lower than $y_{\rm best}$ is maximized, where $y_{\rm best}$ denotes the current best formation energy among the training data. This means that structure $\bm{d}_{i^{\prime}}$ is selected by maximizing the probability as

where $\Phi\left(y-\mu(\bm{d}_{*})/\sigma_{*}\right)$ denotes the cumulative distribution function of normal distribution $N(\mu,\sigma^{2})$.

Figure 4 shows the convergence behavior of the convex hull of formation energy in Bayesian optimization. We measure the convergence degree of the convex hull using the area of the convex hull below the line indicating zero formation energy. As can be seen in Fig. 4, the convex hull converges well in all systems. The numbers of DFT calculations required for the convergence are 1076, 1085, 1032, 701, 945, and 886 in CaCuO_3-x, SrCuO_3-x, BaCuO_3-x, ScCuO_3-x, YCuO_3-x, and LaCuO_3-x, respectively.

Figure 5 shows the formation energies of the oxygen-deficient structures with up to $n=8$ selected by Bayesian optimization including the structures in the initial training dataset of Bayesian optimization. The broken line in Fig. 5 indicates the convex hull of formation energy representing the set of oxygen-deficient perovskite structures that are more stable than the other oxygen-deficient perovskite structures. We have found three, five, eight, five, four and three oxygen-deficient perovskite structures on the convex hull in CaCuO_3-x, SrCuO_3-x, BaCuO_3-x, ScCuO_3-x, YCuO_3-x, and LaCuO_3-x, respectively. Four of the oxygen-deficient perovskite structures on the convex hull, i.e., CaCuO₂, SrCuO₂, SrCuO_2.5, and LaCuO_2.5, are identical to the experimental structures in the literature Chen et al. (1996); Bringley et al. (1990); Karpinski et al. (1994); Takano et al. (1989). They are denoted by p-1 and p-2 in Fig. 5.

Although the convex hull of oxygen-deficient perovskite structures has many vertices or oxygen-deficient perovskite structures, experimental structures that are not classified into the perovskite family have also been known in ACuO_3-x systems. In the six ACuO_3-x systems, seven structure prototypes with no partial occupancy of Sc₂Cu₂O₅, Y₂Cu₂O₅, Nd₂Cu₂O₅, SrCuO₂, BaCuO₂, delafosite FeCuO₂, and CaCuO₂ types are known in addition to the oxygen-deficient perovskite Ca₂Mn₂O₅ type (p-1) and LaNiO₂ type (p-2). Therefore, we estimate the formation energy for the seven prototypes in each system by DFT calculation.

Figure 5 shows the formation energy distribution of the experimental prototype structures. The convex hull obtained from the oxygen-deficient perovskite structures and the experimental prototype structures is also shown in Fig. 5. In all compositions of ACuO_2.5, the Sc₂Cu₂O₅ type structure (e-1) is stable although the Y₂Cu₂O₅-type structure is competitive with the Sc₂Cu₂O₅-type structure in LaCuO_2.5. Also, the SrCuO₂ type (e-3) is stable in SrCuO₂ and BaCuO₂, whereas the delafosite FeCuO₂ type (e-2) is stable in ScCuO₂, YCuO₂, and LaCuO₂. The stable structure of CaCuO₂ is of the LaNiO₂ type (p-2), which is consistent with the experimental structure. On the other hand, most of the oxygen-deficient perovskite structures are not vertices of the convex hull, which indicates that they are not thermodynamically stable at 0 GPa.

To understand the common features of the oxygen-deficient perovskite structures on the convex hull, we carry out a machine learning clustering of the oxygen-deficient perovskite structures. Here, we introduce the Warren–Cowley short-range order (WC-SRO) parameter to represent the oxygen-deficient perovskite structures. The WC-SRO parameters can be derived from the correlation functions used in both the CE method and Bayesian optimization. In ACuO_3-x, the WC-SRO parameter of pair $\alpha$ for structure $\bm{d}$ is defined as

where $\varphi_{\alpha}(\bm{d})$ denotes the correlation function of pair $\alpha$ and $q=(3-2x)/3$.

We employ a hierarchical clustering approach using Pearson correlation coefficients between the oxygen-deficient perovskite structures. For the distance metric to define the similarity between the structures, we adopt a complete-linkage clustering, which uses the maximum correlation-based distances between all observations of the two sets. Figure 6 (a) shows Pearson correlation coefficients between the oxygen-deficient perovskite structures on the convex hull, represented by WC-SRO parameters of the eleven pairs shown in Fig. 6 (b). As can be seen in Fig. 6 (a), the hierarchical clustering classifies the oxygen-deficient perovskite structures into four large groups.

Figure 6 (c) shows the WC-SRO parameters averaged over the oxygen-deficient perovskite structures belonging to each group. Group 1 shows positive large values of the WC-SRO parameters for the second, tenth, and eleventh pairs; hence, most of those pairs are composed of the same species in the structures of group 1. This indicates that the structures in group 1 consist mainly of planar fourfold copper units. An example is SrCuO₂ with the p-2 structure illustrated in Fig. 5, where all copper atoms are planar fourfold units. Also, group 4 is characterized by a positive value of the WC-SRO parameter of the nearest pair along the $\langle 111\rangle$ direction (the eighth pair). This means that $\langle 111\rangle$ nearest pairs of vacancies are components of the structures in group 4. An example is SrCuO_2.5 with the p-1 structure illustrated in Fig. 5.

Note that we evaluate WC-SRO parameters of the oxygen-deficient perovskite structures with the ideal cubic perovskite lattice. Some of the structures optimized by the DFT calculation significantly differ from their ideal structures.

As shown in Sec. IV.3, most of the oxygen-deficient perovskite structures are not stable at 0 GPa. On the other hand, oxygen-deficient perovskite compounds such as CaCuO₂, SrCuO₂, and SrCuO_2.5 have been synthesized under high-pressure conditions Chen et al. (1996); Karpinski et al. (1994); Takano et al. (1989). Therefore, we investigate the pressure dependence of the phase stability of the oxygen-deficient perovskite structures. We restrict the DFT calculation to the experimental prototype structures and the oxygen-deficient perovskite structures satisfying the condition that their energies measured from the convex hull of the oxygen-deficient perovskite structures are less than 0.1 eV/ACuO_3-x.

Figure 7 (a) shows the formation enthalpy distribution of the experimental structures and the oxygen-deficient perovskite structures at 10 GPa. Contrary to the phase stability at 0 GPa, all of the experimental oxygen-deficient perovskite structures, namely, CaCuO₂, SrCuO₂, SrCuO_2.5, LaCuO_2.5, and LaCuO_2.6, are stable at 10 GPa. In addition, the convex hull indicates that some of the oxygen-deficient perovskite structures are stable at 10 GPa in SrCuO_3-x, BaCuO_3-x, and LaCuO_3-x. Figure 7 (b) shows the crystal structures of the oxygen-deficient perovskites with the composition $x\geq 0.4$ on or almost on the convex hull, where only LaCuO_2.6 (p-7) is experimentally reported Bringley et al. (1990). The crystal structures with the other compositions are also shown in Appendix. In CaCuO_3-x, ScCuO_3-x, and YCuO_3-x, the structures on the convex hull at 10 GPa are almost the same as those at 0 GPa. However, some oxygen-deficient perovskite structures such as CaCuO_2.5 and YCuO_2.571 are close to the convex hull. They are expected to be stable at higher pressures.

As can be seen in Fig. 7, the oxygen configuration in the predicted structure of SrCuO_2.375 (p-3) is similar to that in the experimental structure of La_0.5Ca_0.5CuO_2.5 Guloy et al. (1994). The difference between them is the Cu–O polyhedron located at the center of the unit cell, i.e., they are planar fourfold and octahedral sixfold structures, respectively. The predicted structure of SrCuO_2.25 (p-4) is also similar to the structures of SrCuO_2.375 (p-3) and La_0.5Ca_0.5CuO_2.5. At the same time, since La_0.5Ca_0.5CuO_2.25 Guloy et al. (1994) and La_0.25Sr_0.75CuO_2.25 Fu et al. (1990) are reported as existing compounds, it seems to be natural to compare SrCuO_2.25 (p-4) with them. However, their crystal structures remain unclarified. The results in Ref. Fu et al., 1990 only reveal that La_0.25Sr_0.75CuO_2.25 has the space group of $P$42₁2 or $P\bar{4}$2${}_{1}m$. Nevertheless, they are subgroups of the space group $P$4/$mbm$ of SrCuO_2.25 (p-4) with index 2. Hence, the oxygen configuration in SrCuO_2.25 (p-4) may be identical to those in La_0.5Ca_0.5CuO_2.25 and La_0.25Sr_0.75CuO_2.25.

Besides, we compare the oxygen-deficient structures of BaCuO_2.333 with that of YBa₂Cu₃O_7-d (YBCO), because the present structure dataset contains the oxygen configuration of YBCO. As shown in Appendix (Fig. 9), although the YBCO-like structure cannot be found as an oxygen-deficient structure with the lowest enthalpy at $x=2/3$ for BaCuO_3-x, the formation enthalpy of the YBCO-like structure is only 52 meV/BaCuO_3-x larger than that of the lowest formation enthalpy. This implies that the formation of the YBCO-type structure consisting of the alternate stacking of Cu–O chains and Cu–O double layers is dominated not by the cation ordering but by the oxygen vacancy ordering.

Finally, we demonstrate the pressure dependence of phase stability between the oxygen-deficient perovskite structures and the experimental structures. For each structure, the pressure–enthalpy curve is obtained from the energy–volume curve fitted to the Vinet equation of states Vinet et al. (1987) using the DFT energies of ten structures constructed by expansions/contractions of the equilibrium structure. Figure 8 shows the pressure–enthalpy curves of the oxygen-deficient perovskite structure (p-1) and Sc₂Cu₂O₅-type structure (e-1) in SrCuO_2.5. As can be seen in Fig. 8, the oxygen-deficient perovskite structure becomes stable as the pressure increases. It is found that the transition pressure from the Sc₂Cu₂O₅-type structure to the oxygen-deficient perovskite p-1 structure is 5.6 GPa in SrCuO_2.5. Similarly, LaCuO_2.5 involves the phase transition from the Sc₂Cu₂O₅-type structure to the Nd₂Cu₂O₅-type structure at 1.9 GPa and the phase transition from the Nd₂Cu₂O₅-type structure to the oxygen-deficient perovskite p-1 structure at 4.4 GPa.