Evolutionary Algorithm-Based Crystal Structure Prediction of CuxZnyOz Ternary Oxides

Binary zinc(II) oxide (ZnO) and copper(II) oxide (CuO) are used in a number of applications, including optoelectronic and semiconductor applications. However, no crystal structures have been reported for ternary Cu-Zn-O oxides. In that context, we investigated the structural characteristics and thermodynamics of CuxZnyOz ternary oxides to map their experimental feasibility. We combined evolutionary crystal structure prediction and quantum chemical methods to investigate potential CuxZnyOz ternary oxides. The USPEX algorithm and density functional theory were used to screen over 4000 crystal structures with different stoichiometries. When comparing compositions with non-magnetic CuI ions, magnetic CuII ions, and mixed CuI-CuII compositions, the magnetic Cu2Zn2O4 system is thermodynamically the most favorable. At ambient pressures, the thermodynamically most favorable ternary crystal structure is still 2.8 kJ/mol per atom higher in Gibbs free energy compared to experimentally known binary phases. The results suggest that thermodynamics of the hypothetical CuxZnyOz ternary oxides should also be evaluated at high pressures. The predicted ternary materials are indirect band gap semiconductors.


Introduction
d-metal oxides are used in a variety of technological applications such as catalysis [1], optoelectronic applications [2,3], and thermoelectric energy conversion [4]. Recently, many copper-based and iron-based ternary, Earth-abundant oxides have been investigated as photocathode materials in photoelectrochemical cells for water splitting [5]. Considering the widely useful optoelectronic and semiconducting properties of binary ZnO and CuO oxides, one could consider that ternary Cu-Zn-O oxides could possess interesting properties for photoelectrochemical applications. However, ternary Cu x Zn y O z oxides reported to date are typically Cu-doped ZnO or Zn-doped CuO and true ternary phases like CuZnO 2 are not known. For example, Jin et al. studied the hexagonal-to-cubic phase transition of Zn 0.854 Cu 0.146 O under high pressure, where the material had a ZnO-like wurtzite or rock salt crystal structure [6]. Prabhakaran and Boothroyd studied the single-crystal growth of Zn-doped CuO, Cu 1−x Zn x O, where x was 0.05 or 0.1 and the crystal structure was the monoclinic CuO (tenorite) crystal structure [7]. Analogously, Bououdina et al. studied nanostructured wurtzite-ZnO doped with up to 10% Cu using neutron diffraction and ab initio calculations [8]. The optical properties of such Cu-doped ZnO were also studied with density functional theory (DFT) by Volnianska and Bogusławski [9]. Wattoo et al. investigated the effect of zinc concentration on the physical properties of CuO, where x was varied between 0 and 0.08 in Cu 1−x Zn x O [10]. When Amaral et al. studied the Zn doping effect on the structural properties of CuO, their Cu 1−x Zn x O samples showed small amounts of spurious phases above x = 0.05 [11]. Overall, it appears that Cu x Zn y O z ternary phases have been studied only for compositions where one metal corresponds to about 90% of the total metal content.

Results
To span the configuration space of an unknown crystal structure such as Cu x Zn y O z , different compositions must be considered. For non-magnetic Cu x Zn y O z with Cu(I) and Zn(II) ions, we investigated the following compositions for the primitive cell (number of screened crystal structures in parentheses): Cu 2 ZnO 2 (~600), Cu 4 Zn 2 O 4 (~1600), Cu 6 Zn 3 O 6 (~400), Cu 8 Zn 4 O 8 (~400). The compositions were chosen by increasing the number of Cu(I) ions stepwise from 2 to 8.
Magnetic compositions for the crystal structure predictions were divided into two groups. The first group represents compositions with only Cu(II) ions and Zn(II) ions: CuZnO 2 (~200 crystal structures) and Cu 2 Zn 2 O 4 (~400 crystal structures). In other words, we first studied the simplest 1:1:2 composition and then a doubled composition with two Cu(II) ions in the primitive cell. A similar strategy was used for the second group of magnetic crystal structures with both Cu (I) and Cu(II) ions: simplest composition  Cu II Cu I  2 Zn 2 O 4 (~400 crystal structures), double composition Cu II  2 Cu I  4 Zn 4 O 8 (~200 crystal  structures), and finally, the composition Cu II  2 Cu I  2 Zn 4 O 7 (~400 crystal structures). The Cu II 2 Cu I 2 Zn 4 O 7 composition was included in the screening due to the equal number of Cu and Zn in the crystal structure. In general, the number of formula units used in the USPEX simulations is limited by the available computational resources, and evolutionary searches with two magnetic Cu II ions already proved to be rather demanding with the hybrid DFT-PBE0 method.
Because the studied hypothetical Cu x Zn y O z ternary materials are completely new and there is no information on their magnetic ground states, we had to investigate both ferromagnetic and antiferromagnetic spin settings to find the magnetic ground state of each crystal structure. However, we noted that many of the Cu 2 Zn 2 O 4 crystal structures that we tried to predict an antiferromagnetic ground state ended up as analogs of ferromagnetic CuZnO 2 crystal structures with the same magnetic moments, space groups, and relative energies. The magnetic moments of the spin-up and spin-down Cu(II) ions in the predicted antiferromagnetic Cu 2 Zn 2 O 4 structures were found to be identical to each other. To simplify further simulations with larger unit cells in the case of mixed Cu I -Cu II compositions, we carried out the crystal structure predictions only for ferromagnetic configurations.
Relative energies (∆E) of the predicted Cu x Zn y O z crystal structures were estimated by comparing their total energy with the related binary oxides using the following equation: In particular, in case of non-magnetic Cu 2 ZnO 2 , relative energies per atom were estimated in the following way: In case of magnetic CuZnO 2 , relative energies were estimated in the following way: In case of mixed Cu II Cu I 2 Zn 2 O 4 , relative energies were estimated in the following way: The relative energies were evaluated with Equations (1)-(4) for all the lowest-energy crystal structures produced by the USPEX simulations. In the detailed discussion here, we focus on crystal structures where ∆E was found to be smaller than 10 kJ mol −1 per atom. For these crystal structures, we also estimated the Gibbs free energy ∆G, where the phonon contributions were evaluated in two different ways: (1) at the Γ point only, and (2) based on phonon supercells. Relative Gibbs free energies ∆G were obtained analogously to ∆E, using Equations (1)-(4). Figure 1 illustrates the relative total energies and Gibbs free energies of the predicted lowest-energy Cu x Zn y O z crystal structures and Table 1 lists detailed information for them. Optimized lattice parameters are presented in Table 2. The crystal structures are labeled using a scheme that is explained in the caption of Figure 1. The unit cell parameters and atomic coordinates of the predicted lowest-energy structures are given as Supplementary Materials in CIF format.
It seems clear that magnetic Cu x Zn y O z crystal structures possess lower energies than non-magnetic crystal structures. The same is true for Gibbs free energies. None of the mixed crystal structures with both Cu II and Cu I ions were found to have ∆E smaller than 10 kJ mol −1 per atom. The predicted non-magnetic crystal structures possess a smaller band gap (2.2-2.4 eV) than the magnetic compounds (3.1-3.4 eV).   Table 1.  The relative energies are listed in Table 1. Table 1. Space group, relative energy ∆E, relative Gibbs free energy ∆G, band gap, Cu II magnetic moment, and coordination numbers for the lowest-energy Cu x Zn y O z crystal structures predicted in this study. The crystal structures are labeled using a scheme that is explained in the caption of Figure 1. Coordination numbers of the metal atoms also vary depending on the composition. In the reference binary phases, the coordination numbers of the metal atoms are as follows: Cu in Cu 2 O is two-coordinated (linear), Cu in CuO is four-coordinated (square planar), and Zn in ZnO is four-coordinated (tetrahedral). Non-magnetic crystal structures NM1-NM4 contain linearly coordinated Cu I ions with a coordination number of two, similar to Cu 2 O. In the case of NM3, the linear Cu shows clearly bent coordination, with O-Cu-O angles of about 158 • . In crystal structures NM1-NM3, the Zn II ions are four-coordinated with tetrahedral coordination, but the tetrahedrals are somewhat distorted from ideal. In NM4, the Zn II atoms are somewhat unexpectedly five-coordinated with square pyramidal coordination. This structure has the highest relative energy of the non-magnetic crystal structures, but its Gibbs free energy is lower compared to NM3. Table 2. Optimized lattice parameters of the lowest-energy Cu x Zn y O z crystal structures predicted in this study. The crystal structures are labeled using a scheme that is explained in the caption of Figure 1. The lattice parameters of the experimentally known binary oxides are provided, as well a . The lowest-energy magnetic crystal structure M1 has Cu II ions with a coordination number of four and square planar coordination, similar to CuO. Zn II ions in the structure have a coordination number of six, with slightly distorted octahedral coordination. The coordination number of both Cu II and Zn II ions in the magnetic crystal structures M2-M4 is six, with octahedral coordination. Finally, in the magnetic crystal structure M5, both Cu II and Zn II ions show a somewhat unexpected five-coordinated square pyramidal coordination sphere.

Crystal Structure
As mentioned above, we found that the ferromagnetic and antiferromagnetic configurations of magnetic Cu x Zn y O z crystal structures are practically isoenergetic. In fact, ferromagnetic lowest-energy CuZnO 2 was found to have the same space group and structural parameters as antiferromagnetic M1. Energetically, the FM and AFM (M1) crystal structures are also very close to each other (for AFM M1, ∆E = 3.5, ∆G Γ = 2.7, ∆G = 3.0 kJ mol -1 per atom). Thus, to save computational resources, the harmonic frequency calculations of some crystal structures were only carried out with ferromagnetic spin settings. In particular, the antiferromagnetic M3 structure has space group P-1, whereas the ferromagnetic structure possesses space group C2/c. Antiferromagnetic M5 has space group Pm and its ferromagnetic configuration has space group Pmn2 1 . In the case of other M crystal structures, no changes in space groups were identified if the magnetic configuration was changed from antiferromagnetic to ferromagnetic.
∆E and ∆G estimated at the Γ-point and based on the phonon supercell calculations are in good agreement for the NM crystal structures, the maximum difference being 0.8 kJ mol −1 per atom. Magnetic M crystal structures have similar trends for ∆E and ∆G. Only in the case of M3 was the difference between ∆E and ∆G Γ found to be slightly larger (1.4 kJ mol −1 per atom). ∆G evaluated at the Γ-point and based on phonon supercell calculations are typically in good agreement with each other. The largest difference is found in the case of crystal structure M3, where it is 1.0 kJ mol −1 per atom.
We also evaluated the Gibbs free energy for the studied crystal structures at 300, 400, 500, 600, and 700 K (Figure 2). The goal was to see whether the ternary oxides would become more favorable at higher temperatures. ∆G evaluated at the Γ-point only decreases for all structures with increasing temperature (top panel in Figure 2). However, when more accurate phonon supercell calculations are considered, the changes are smaller (bottom panel in Figure 2). In fact, for non-magnetic structures, increase in the temperature leads to a small increase in ∆G. The largest effect on ∆G Γ can be seen for M3, where ∆G Γ decreases from 3.9 (298 K) to 1.7 kJ mol −1 per atom (700 K). However, ∆G evaluated from more accurate phonon supercell calculations shows that ∆G remains at about~4 kJ mol −1 per atom even at 700 K. Overall, increasing the temperature does not seem to favor the ternary Cu x Zn y O z phases over the experimentally known binary phases Cu 2 O, CuO, and ZnO. However, for the lowest-energy system M1, ∆G with respect to the binary phases is only 2.8 kJ/mol per atom. This is almost identical to the Gibbs free energy difference between graphite and diamond [14]. Therefore, high-pressure studies on the hypothetical ternary phases should be carried out to see whether higher pressures could favor them over the binary phases.
decreases for all structures with increasing temperature (top panel in Figure 2). How when more accurate phonon supercell calculations are considered, the changes smaller (bottom panel in Figure 2). In fact, for non-magnetic structures, increase i temperature leads to a small increase in ΔG. The largest effect on ΔG Г can be seen for where ΔG Г decreases from 3.9 (298 K) to 1.7 kJ mol −1 per atom (700 K). However evaluated from more accurate phonon supercell calculations shows that ΔG remai about ~4 kJ mol −1 per atom even at 700 K. Overall, increasing the temperature doe seem to favor the ternary CuxZnyOz phases over the experimentally known binary ph Cu2O, CuO, and ZnO. However, for the lowest-energy system M1, ΔG with respect t binary phases is only 2.8 kJ/mol per atom. This is almost identical to the Gibbs free en difference between graphite and diamond [14]. Therefore, high-pressure studies o hypothetical ternary phases should be carried out to see whether higher pressures c favor them over the binary phases.  Å and two Zn-O distances are 2.09 and 2.13 Å. M1 also has the largest band gap among the studied structures (3.4 eV). Meanwhile, M2 crystallizes in the monoclinic crystal system with the same C2/c space group as M1. However, the crystal structure is more complicated compared to M1, showing three different Cu-O distances (1.96 Å, 2.05 Å, and 2.38 Å) and two Zn-O distances (2.05 Å and 2.18 Å). The band gap of the M2 structure is 3.1 eV, which is slightly smaller than that of M1. Increasing the temperature does not change the thermodynamic stability order and M1 remains more favorable than M2.  rather similar down to −3 eV, while the contribution from Cu increases and surpasses th of O at around −3 eV. Similarly, in the case of M2, Cu contributes more to the topmo valence bands than Zn and the contribution from Zn stays rather similar down to −3 e while the contribution from Cu increases and surpasses that of O at around −3 eV. F both materials, the conduction bands are dominated by Cu, with contributions from and even smaller contributions from Zn (practically negligible in the case of M2).  Regarding the magnitude of the predicted band gaps, they are likely to be affected by the known behavior of DFT-PBE0 [15]. For systems where the experimental band gap is smaller than 1 eV, DFT-PBE0 typically significantly overestimates the band gap. For band gaps between 2 and 5 eV, DFT-PBE0 produces more reasonable estimates. For binary Cu 2 O, DFT-PBE0/TZVP yields a band gap of 2.39 eV, while an often-cited experimental result is 2.17 eV [16]. For CuO, DFT-PBE0/TZVP yields a band gap of 3.8 eV, which is clearly overestimated in comparison to experimental 1.7 eV [13]. This behavior should be taken into account when considering the band gaps predicted for the hypothetical ternary oxides here.

Methods
Crystal structure predictions were carried out by using USPEX 9.4.4 code for evolutionary crystal structure prediction [17][18][19]. A typical workflow for a USPEX simulation is presented in Figure 5. The only chemical input for the crystal structure prediction is the composition of the studied system, such as Cu 2 Zn 2 O 4 . In addition, technical parameters related to the USPEX crystal structure prediction algorithm are provided, though these are not highly system-dependent. Similar input files were used for different chemical compositions in this study. An example of the USPEX input file used is given in the Supplementary Materials (Table S1). The ternary CuxZnyOz compositions investigated here are explained in detail in Section 2 "Results". At the beginning of the simulation, completely random structures are used as a starting point. The space group is chosen randomly, and the Wyckoff positions are filled randomly. Next, local optimizations with quantum chemical methods are carried out for the crystal structures in the starting population and the energies (enthalpies) of the optimized structures are compared with each other. The fittest (lowest-energy) structures are chosen as parents for the new generation, after which new structures are generated by applying heredity and mutations to the parent structures. A typical USPEX run involves thousands of local optimizations, producing hundreds of crystal structures (detailed numbers reported in Section 2 "Results"). The USPEX procedure can be considered converged when the fittest structure is no longer changing after several (e.g., 10) generations. After convergence, the fittest crystal structures still need to be re-optimized at a higher level of theory, as the local optimizations in the USPEX simulations are typically carried out with lower accuracy due to the vast number of calculations.
Quantum chemical calculations within the USPEX simulations were performed using density functional theory (DFT). The CRYSTAL17 and Quantum Espresso (QE) version 6.0 program packages were utilized for the DFT calculations [20,21]. The PBE exchangecorrelation functional with GBRV ultrasoft pseudopotentials was used for all DFT calculations carried out with QE [22,23]. It is known that standard generalized gradient approximation (GGA) functionals such as PBE can fail to treat the magnetic moments and electronic structure of systems such as strongly correlated d-metal oxides, sometimes even favoring the wrong magnetic ground state [24][25][26][27][28][29][30][31]. Therefore, QE was used only for nonmagnetic CuxZnyOz compositions, and magnetic CuxZnyOz compositions were studied only with CRYSTAL and hybrid density functional methods combined with all-electron basis sets. For CRYSTAL, we used the CRYSTAL interface developed for USPEX, enabling The ternary Cu x Zn y O z compositions investigated here are explained in detail in Section 2 "Results". At the beginning of the simulation, completely random structures are used as a starting point. The space group is chosen randomly, and the Wyckoff positions are filled randomly. Next, local optimizations with quantum chemical methods are carried out for the crystal structures in the starting population and the energies (enthalpies) of the optimized structures are compared with each other. The fittest (lowest-energy) structures are chosen as parents for the new generation, after which new structures are generated by applying heredity and mutations to the parent structures. A typical USPEX run involves thousands of local optimizations, producing hundreds of crystal structures (detailed numbers reported in Section 2 "Results"). The USPEX procedure can be considered converged when the fittest structure is no longer changing after several (e.g., 10) generations. After convergence, the fittest crystal structures still need to be re-optimized at a higher level of theory, as the local optimizations in the USPEX simulations are typically carried out with lower accuracy due to the vast number of calculations.
Quantum chemical calculations within the USPEX simulations were performed using density functional theory (DFT). The CRYSTAL17 and Quantum Espresso (QE) version 6.0 program packages were utilized for the DFT calculations [20,21]. The PBE exchangecorrelation functional with GBRV ultrasoft pseudopotentials was used for all DFT calculations carried out with QE [22,23]. It is known that standard generalized gradient approximation (GGA) functionals such as PBE can fail to treat the magnetic moments and electronic structure of systems such as strongly correlated d-metal oxides, sometimes even favoring the wrong magnetic ground state [24][25][26][27][28][29][30][31]. Therefore, QE was used only for non-magnetic Cu x Zn y O z compositions, and magnetic Cu x Zn y O z compositions were studied only with CRYSTAL and hybrid density functional methods combined with allelectron basis sets. For CRYSTAL, we used the CRYSTAL interface developed for USPEX, enabling the prediction of magnetic ordering in addition to the crystal structure [32]. The hybrid DFT-PBE0 functional with 25% Hartree-Fock exchange was utilized in CRYSTAL calculations [22,33]. Overall, for a 3d metal such as Cu, the use of hybrid DFT over GGA or GGA + U is expected to increase the accuracy of the predictions [16,24,[34][35][36][37][38]. Therefore, we also carried out USPEX/CRYSTAL simulations for non-magnetic structures to ensure the reliability of the USPEX/QE results. Successful use of the USPEX/CRYSTAL for nonmagnetic and magnetic d-metal oxides and fluorides has been reported earlier [32,39,40]. All-electron, Gaussian-type split-valence + polarization (SVP) basis sets based on Karlsruhe def2 basis sets were used within the USPEX/CRYSTAL calculations [13,41]. To accelerate the evolutionary searches, the local structure optimizations within USPEX were carried out by using relatively loose convergence criteria. Space group P1 was used for the local structure optimizations. The CRYSTAL and QE input files used within the USPEX simulations are given as Supplementary Materials (Tables S2 and S3).
The lowest-energy structures produced by USPEX/CRYSTAL were re-optimized at the DFT-PBE0 level of theory using tighter convergence criteria ("accurate quantum chemical calculations" in Figure 5). Structural optimizations were carried out in the space groups found with the FINDSYM program package [42]. All-electron, Gaussian-type triple-ζ-valence + polarization (TZVP) basis sets based on the Karlsruhe def2 basis set were used within the crystal structure predictions [13,41]. Cu(I) ions are expected to show weak "cuprophilic" d 10 -d 10 interactions [43][44][45][46], which were taken into account using Grimme's D3 dispersion correction with zero-damping (ZD) [47,48]. For the reoptimization, the reciprocal space k-point meshes were chosen depending on the magnitude of the corresponding direct space lattice parameter d: d < 4 Å → 12 k-points along d; 4 Å < d < 6 Å → 8 k-points, 6 Å < d < 8 Å → 6 k-points; 8 Å < d < 12 Å → 4 k-points; d > 12 Å → 2 k-points. Tightened tolerance factors (TOLINTEG) of 8,8,8,8, and 16 were used for the evaluation of the Coulomb and exchange integrals. All reported structures were confirmed to be true local minima by means of a harmonic frequency calculation [49,50]. Calculations of phonon dispersions were carried out for supercells where a, b, and c lattice parameters were about 10 Å. The optimized structures of the lowest-energy ternary Cu x Zn y O z crystal structures are included as Supplementary Materials in CIF format.
For estimating the energetics and thermodynamics of the predicted ternary Cu x Zn y O z crystal structures, we also optimized the structures of binary Cu 2 O, CuO, and ZnO at the DFT-PBE0-D3(ZD)/TZVP level of theory. All relative energies reported in this paper have been obtained at this level of theory. CuO was studied using an antiferromagnetic ground state [51][52][53][54][55]. Monkhorst-Pack-type 8 × 8 × 8, 4 × 8 × 4, and 12 × 12 × 8 k-point meshes were used for Cu 2 O, CuO, and ZnO, respectively. Including the D3(ZD) dispersion correction does not change the lattice parameters of the binary oxides significantly in comparison to DFT-PBE0 without dispersion correction [13]. For Cu 2 O, the optimized lattice parameter a is 4.28 Å with D3 and 4.32 Å without D3. For ZnO, the optimized lattice parameters are 3.25 Å and 5.19 Å with D3 and 3.27 Å and 5.21 Å without D3.

Conclusions
We combined evolutionary crystal structure prediction and quantum chemical methods to discover potential Cu x Zn y O z ternary oxides. Over 4000 crystal structures with different stoichiometries were screened in the ternary system. When comparing compositions with non-magnetic Cu I ions, magnetic Cu II ions, and mixed Cu I -Cu II compositions, the magnetic Cu 2 Zn 2 O 4 system is thermodynamically the most favorable. At ambient pressures, the thermodynamically most favorable crystal structure (M1) is still 2.8 kJ/mol per atom higher in Gibbs free energy compared to experimentally known binary phases. However, this energy difference is almost identical to the Gibbs free energy difference between graphite and diamond [14], suggesting that thermodynamics of the hypothetical ternary systems studied here should also be evaluated at high pressures.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules28165986/s1, Table S1: USPEX input file for Cu x Zn y O z crystal structure prediction; Table S2: Quantum Espresso input files used in USPEX simulations; Table S3: CRYSTAL input files used in USPEX simulations. Crystal structure data are included in CIF format.