1. Introduction
The discovery of p-type transparent conductivity of CuAlO
2 films with delafossite structure [
1] was a milestone in the race for finding p-type transparent conducting oxides (TCOs). These kinds of materials, among which electronic doped ZnO and SnO
2 are perhaps the most representative [
2,
3,
4,
5], have some very unusual physical properties; they are transparent, similar to glass, but conductive almost like a metal. They are widely used in a variety of electronic applications, from panel displays and touch screens to solar cell applications. In the design of heterostructure solar cells, they are mainly used as electron (ETL) and hole (HTL) transport layers. Problems arise with the fact that most TCOs are n-type conductors, and finding p-type conductors with wide bandgap and good conductivity has been elusive.
It has been discussed [
3,
4,
6] that for semiconductor oxides, the mobility of doping holes at the top of the valence band, typically filled by O-
bands, tends to be very low due to the strong localization of these levels. Possible solutions to achieve higher mobility of holes would be the use of materials for which the valence band edge is not O-
but cations with
d electrons at the top of the band or modifying the valence band edge by mixing these orbitals with appropriate cations that have energy-filled levels comparable to those of the O-
. Allegedly, this would be effective in reducing the strong Coulomb interaction of oxygens and therefore in delocalizing positive holes. The choice of a monovalent copper Cu
+ as a candidate, with filled
orbitals in a closed shell after giving its
electron, and energy of the Cu-
levels compared to those of the O-
leads to the idea of Cu oxides as good candidates to overcome the localization of holes via covalent bonding between the Cu
+ and the O
−2 ions. However, the simplest Cu oxide, CuO
2, has a rather small bandgap and it is not transparent in visible light, despite not having a considerably large conductivity.
Therefore, the discovery of p-type conductivity in CuAlO
2 triggered an increased interest in the experimental and theoretical characterization of this compound [
3,
7,
8,
9]. CuAlO
2 was followed by the family of CuM
III−AO
2 metals, in which the M
III−A cation is a metal from the group III-A in the periodic table; CuGaO
2 was found to also have p-type conductivity and CuInO
2 was found to exhibit the ability to be both n- and p-type doped [
2,
5,
10,
11,
12,
13,
14,
15]. Also, CuScO
2 [
16,
17] and CuYO
2 [
18,
19] have been studied, as well as CuCrO
2, which has been doped with Mg, exhibiting the highest conductivity but not good transparency [
20,
21,
22,
23,
24]. To date, numerous Cu ternary oxides with delafossite structure CuMO
2 (M = B, Al, Ga, In, Fe, Cr, Sc, Y, etc.), and even some with general formula AMO
2 (where A = Ag, Cu), have been reported with p-type conductivity as high as
S/cm and high optical transparency (50–85%), depending not only on their chemical compositions but also grain sizes or film deposition methods [
4,
25,
26,
27,
28,
29,
30,
31,
32]. Theoretically, special attention has been received by the Fe-delafossite for the difficulties in its calculation [
33,
34,
35,
36].
The effect of doping with divalent metals in the delafossite compounds has also been investigated. Other studies of mixed structures are common for different kinds of applications and materials [
37,
38] and are useful not only to predict properties in novel materials but also for the understanding of underlying mechanisms of the different properties of interest withing given materials. The race to achieve good transparency and conductivity has included the study of Mg-doped delafossites [
4,
21,
22,
24,
29], Ni- or Ca-doped delafossites [
20,
39], as well as doping with another trivalent metal in structures with high percentages of substitution [
27,
40].
In this paper, we aim to study a variety of Cu-based delafossite semiconductors as promising materials for photovoltaic applications in the design of heterostructure solar cells, as hole transport layers (HTLs). We focus on the CuMO
2 structure, with M = Al, Ga, In, Co, Sc, Y, Fe, Cr. Geometrical and electronic characterization were performed using the density functional theory for these materials, as described in
Section 2: Methodology. The main geometrical characteristics of the two phases in which these delafossites can be found are discussed in
Section 3: Structural characterization. A discussion about bandgap values, density of states, and band structure can be found in
Section 4: Electronic Structure. Then, the effects in the geometry and bandgap when one of the trivalent metals within the oxides is substituted are tested.
Section 5 shows results for CuM
xN
1−xO
2 structures, with N, M being the same species previously mentioned, and
. Different relative concentrations for N, M require larger supercells, making the more accurate calculations using hybrid exchange correlation functionals such as HSE prohibitive. Thus, a discussion on statistical models and the use of techniques based on machine learning to describe the correlation between the bandgap values described with the GGA and HSE functionals is presented in
Section 6. Finally, some concluding remarks highlighting the fundamental properties of these materials are given.
2. Methodology
Since theoretical simulations have proven to be useful for the understanding of the properties of systems at the atomic level, first-principles calculations were performed to gain an insight into the structural and electronic properties of these Cu-based ternary oxides with delafossite structures.
For this purpose, the atomic positions were optimized and the electronic properties were computed within the density functional theory (DFT) approach, based on the framework of the generalized Kohn–Sham scheme [
41,
42] in combination with the projector augmented wave (PAW) method [
43] for the description of the core electrons and the generalized gradient approximation (GGA) within the Perdew–Burke–Ernzerhof (PBE) formalism [
44] for the exchange and correlation term, as implemented in the Vienna ab initio simulation package (VASP) [
45,
46,
47]. It is well-known that GGA typically underestimates gap values, especially when localized orbitals such as metal
d or
f states are present. Hence, for a more accurate electronic description, the Heyd–Scuseria–Ernzerhof hybrid functional with the modified fraction of screened short-range Hartree–Fock exchange (HSE06) [
48,
49,
50] was also used to calculate the density of states and band energies. Hybrid functional calculation results in a more accurate description of electronic properties for some systems compared to simple DFT calculations with a generalized gradient approximation for the exchange and correlation term in the Kohn–Sham scheme. The higher accuracy, however, is reached at the cost of a higher computational time needed to achieve convergence.
In the search for an alternative method to hold acceptable accuracy with less demanding computational resources, calculations using DFT corrected for on-site Coulomb interactions (DFT + U) [
51] were also performed. The inclusion of the so-called Hubbard term (+U) provides an approximate correction for strongly correlated
d or
f electrons that are usually inaccurately described with standard GGA functionals. The choice of the U parameter is typically an ad hoc choice, sometimes with empirical information. In the present study, the U value has been chosen such that the valence band description in GGA + U reproduces the shape of the band in the density of state curve obtained with HSE06 calculations, i.e., the value of U for which the position of peaks in the valence band edge is the same as in the HSE simulations. Even when the HSE gap values are not reproduced, having a correctly described valence band allows further calculations with a reasonable computational time cost.
Table 1 shows the U values used in this work.
To model the tetragonal phase of CuMO2, a unit cell with 12 atoms (Cu3M3O6) was used, while for the hexagonal cell, a smaller cell with 8 atoms (Cu2M2O4) was needed. To model the CuMxN1−xO2 mixed structures when , the same 12-atom cells were used. However, for the machine learning techniques tested, a small number of calculations were run for which and , and for these structures, 1 × 2 × 1 supercells, containing 24 atoms, were built.
The electronic wave functions were expanded in a plane wave basis setup to a kinetic energy cutoff of 450 eV, and the atomic positions were optimized using the conjugate gradient method up until the forces on each atom were less than 0.01 eV
and the energy convergence less than
eV for the relaxation and GGA calculations, and less than
eV for the calculations with the hybrid functional. The structures with the heaviest atoms—the delafossites, including In and Y—were optimized, including spin–orbit coupling (SOC) [
52], but the
falls under
eV, leading to the conclusion that the effect of spin–orbit coupling does not need to be taken into account with these structures. For the Brillouin zone integration, a 12 × 12 × 6 Monkhorst–Pack scheme k-point mesh was used [
53] both for the optimization and electronic structure calculation.
As mentioned before, due to the high cost of running HSE calculations for each of these structures to obtain accurate bandgap values, we also explored the use of statistical models and more complex techniques based on machine learning to describe the correlation between the bandgap values described with the GGA functional as a first approximation, and those obtained with the HSE functional over already relaxed structures.
The applicability of this approach lies not only in the significant reduction in computational cost when calculating the bandgap values of these materials. It would also allow us to extend the study to new configurations with this delafossite structure which, due to their different relative concentration of elements occupying cationic positions, may require larger supercells for modeling. This would make the cost of HSE calculations practically prohibitive due to the greater number of atoms, but through the implementation of this approach, a simple GGA calculation would be enough to infer a meaningful value of the bandgap for a given configuration.
Several models and combinations of chemical and structural information to describe the compounds have been tested, searching for the best description of the structure–property correlation in this group of materials. In the present manuscript, only the main metrics of the best-performing pair of descriptors (set of features or parameters that better describe the system in an appropriate way for its digestion by the statistical model) and models are shown, together with some remarkable comments, in the corresponding section.
3. Structural Characterization
The delafossite structure of these Cu-based ternary oxides can be appreciated fairly in
Figure 1. It shows the atom positions and geometrical structure of the converged unit cell for the CuAlO
2 delafossite, as a prototype of the structure. The aforementioned delafossite structure consists of close-packed alternate layers of Cu atoms linearly coordinated with O atoms, forming O-Cu-O dumbbells parallel to the
z axis, and M-centered octahedra (M = Al, Ga, In, Co, Sc, Y, Cr, Fe) parallel to the
xy plane. The oxygens on the O-Cu-O dumbbells coordinate to three metal ions each, while metal ions coordinate to eight oxygens, forming the mentioned octahedra. The stacking can occur following an ABCABC pattern, thus forming a trigonal system with a rhombohedral Bravais lattice unit cell with space group
(symmetry number 166), as shown in
Figure 1a, or with an ABAB stacking pattern, thus showing a hexagonal unit cell with
(no. 194) space group, as shown in
Figure 1b.
3.1. Rhombohedral Structures
The structural parameters of the unit cells for the
phase were calculated and
Table 2 summarizes the optimized lattice parameters for all the studied structures, with experimentally measured values for some of the compounds and the deviation percentage that our simulated value represents with respect with the measured one.
Our simulations are in good agreement with experimentally measured values for the cells, with deviations that fall under 2% for almost all the compounds, except for the , which exhibits the highest deviation for both a and c parameters. Particularly good agreement between experimental values and our simulations was obtained for -, -, and - delafossites.
The role of the trivalent metal ion in the structural and electronic properties of these materials has been discussed [
16,
19,
28,
58]. The Cu-Cu distance (equal to the
a lattice constant, since Cu atoms are in the corners of the unit cell) has been pointed out as playing a relevant role for p-type mobility [
58] with holes hopping between Cu ions. On the other hand, the size of the metal ion seems decisive for the
a constant. Note, in
Table 2, the trends followed by the
a constant for CuAlO
2, CuGaO
2, CuInO
2, CuScO
2, and CuYO
2, for example. It is strictly true for CuAlO
2, CuGaO
2, and CuInO
2, as well as for CuScO
2 and CuYO
2 in separate groups. Put together, with CuInO
2 and CuScO
2, the trend to be expected is not clear; however, both delafossites show similar lattice constants, while In
3+ and Sc
3+ ions also exhibit similar sizes according to [
59]. The same can be said for CuCrO
2 and CuFeO
2 if we look for the Fe
3+ size in the high-spin configuration. An unexpectedly low value of
a compared to what could be intuitively expected is the one for CuCoO
2, which seems not to be an error of the simulation, since the experimental value [
55] is very similar. However, it is consistent with the fact that the electronic structure calculated for this delafossite shows no spin polarization, indicating the Co
3+ is in its low-spin configuration, and again the ion size of Co
3+ in the low-spin configuration, as reported in [
59], is very similar to the size of Al
3+, similar to the
a lattice constant.
Therefore, the size of the metal ion, which is decisive to the
a lattice constant, should play an important role in the conductivity of the delafossites, especially those with p-type conductivity. This might be the explanation behind the fact that undoped CuAlO
2 has higher conductivity than both undoped CuScO
2 and CuYO
2 [
16,
28]. More discussion on the effect of the M cation size and
a lattice constant on the conductivity can be found in [
4,
29] for some Mg-mixed Cu-based delafossites. It is of course not the only factor that contributes, with the covalency of the M-O bond being another important factor, as discussed in [
19].
3.2. Hexagonal Structures
As mentioned before, in the hexagonal cells, the same layers of O-Cu-O dumbbells and M-centered octahedral are observed, but with an ABAB stacking pattern. It is expected, then, to find the lattice constant
a to be practically identical between both phases, while the lattice constant
c is expected to have a significantly lower value. The structure parameters found in our calculations, with some experimental values for comparison, are shown in
Table 3.
The mentioned difference in
c lattice constant between both phases is caused by the fact that a smaller distance is enough to represent the structure, because of the stacking sequence. Hence, important distances within the unit cell, such as the distance between two Cu neighbors in the
plane (equal to the
a lattice constant), the distance between Cu-O atoms in the O-Cu-O dumbbells, or the distance from cornered O to centered M in the octahedra are practically the same, as can be appreciated in
Table 4.
Once more, the trend shown in the Cu-Cu distances (equal to the a lattice constant) is in agreement with the size of the trivalent metal ion. The same goes for the M-O distances, as expected. The O-Cu-O dumbbell distances, however, seem to be practically the same for all the structures, varying in a very short range of no more than .
The relative stability of the two different phases is shown in
Table 5. A negative
means the first structure has a lower energy, i.e., is the most stable one: the rhombohedral structure, in all our cases. However, the energy difference for all compounds falls under 0.3 eV, meaning both structures are relatively equally stable, and must be found indistinctly in nature [
62].
5. Mixed Structures
To explore the effects of the substitution of one of the trivalent metals in the delafossite compounds, we have studied the substitution of one of the trivalent metals (Al, Ga, In, Co, Sc, Y, Cr, Fe), having as possible substitutes the same set of metals. Changes in geometry, gap value, and electronic structure have been analyzed in CuMxN1−xO2 structures, where . In this section, as for the pure delafossites, structural results arise from GGA calculations while bandgaps were obtained in HSE simulations.
The inclusion of ions with a larger radius in a delafossite with a smaller M metal is expected to increase both
a and
c parameters, as well as the cell volume. However, in general, we have observed that while the parameter
a increases with the increase in a larger metal in the composition, it is not the case for parameter
c. This can be fairly appreciated in
Figure 11 for a selection of these mixed structures.
As the amount of Ga increases in
Figure 11a, the parameter
a also increases, suggesting a linear tendency. The same can be said for the amount of In in
Figure 11b, In in
Figure 11c, Sc in
Figure 11d, Fe in
Figure 11e, or Cr in
Figure 11f. For every mix, we have calculated only four structures, and a rigorous fit can not be performed with such a small amount of data points.
The case is completely different for the cell parameter
c. It can be appreciated in the inset plots in
Figure 11. In
Figure 11a–c, this parameter also seems to increase with the increase in the amount of larger metal. In these three examples, the substituted metal has a very similar electronic composition as the one to be substituted, since Al, Ga, In belong to the same group in the periodic table. But in all other cases, it is fairly appreciated that there is no clear dependency of
c on the amount of substituted metal, since the dispersion of the four points is notable. This behavior was also found previously [
40] for CuFe
xGa
1−xO
2. The referred study concluded that the
c parameter is almost independent of the concentration of the substituted cation. A possible explanation is that the strong repulsion of the M cations in the
plane affects the MO
6 octahedra, rearranging the oxygen atom positions. As the M cation becomes larger, the M–O distance increases, but the O–O contact distance varies only slightly, and therefore, there is little impact on the
c cell parameter.
The complete set of values for the geometry of these structures can be found in
Table 8.
Bandgap values have been calculated and the DoS structures analyzed. In general, it is impossible to find any clear trend in the behavior of the bandgap values, as was true for the geometry. An exception to this is the CuGa
xAl
1−xO
2 case. As can be appreciated in
Figure 12a, the gaps seem to follow a linear decrease with the increase in Ga concentration. This is not surprising, since in both CuAlO
2 and CuGaO
2 electronic structures, the states of the trivalent metal are very deep in the valence band, and very high in energy in the conduction band. All the interaction near the gap is due to Cu atoms and the tightness of the bond shall only be determined by its Cu-Cu distance. As mentioned before, the structure increases linearly in
a and
c for these delafossites, thus causing the electrons to feel a softer potential and be less tightly bonded to the Cu atoms. A second exception is found by inspecting every delafossite that includes Fe in the substitution. The Fe-states appear so energetically low in the conduction band that the inclusion of any Fe atom causes the bandgap to exhibit very similar values, despite which was the host delafossite (
Figure 12b).
However, apart from these two mentioned examples, there is no clear trend in the gap values that can be analyzed. The complete set of HSE values for these structures can be found in
Table 9.
6. Statistical Models and Machine Learning
The computational effort to calculate HSE bandgaps can be, as has been said, demanding in terms of time and memory needed. This is especially true for structures with relative compositions M/N different than 1/3 (or 2/3). Therefore, we discuss the use of statistical models and techniques based on machine learning to predict HSE bandgaps given the structural information of the material and its GGA calculated bandgap. A step of ab initio calculation is still needed, but the cost of a calculation with a GGA functional is far less expensive than those using HSE functionals.
To use machine learning techniques or other statistical models, a numerical way of fully describing structures is needed. These numerical representations are the descriptors, which must be invariant to symmetry operations, complete and unique, meaning able to distinguish between any two different structures [
67]. The descriptor for which the best performance has been obtained in our current work is the sine matrix [
68], consisting of an
ansatz that mimics the periodicity and the basic features of the elements in the Ewald sum matrix for a periodic system, using a sine function of the crystal coordinates of the atoms.
Apart from the sine matrix description, some other parameters of chemical nature have been added to the descriptor: the covalent radius, electronegativity, and ionization potential of the metal cations occupying M, N positions. Additionally, given the aim of this approach, the GGA results have been added as well as an extra feature of the descriptor. The featurization of the structures was carried out with the Python libraries matminer [
69] and smact [
70,
71,
72]. This featurization was performed for structures with 24 atoms (2 × 1 × 1 rhombohedral supercells), the ones that will be used to test results for new unbalanced configurations, but 12-atom structures (the cells used to model all results in previous sections) were processed as well, in order to check the consistency of results. We ensured that identical models were trained separately for both sets of structures of 12- and 24-atom unit cells, and the average deviation of bandgap results between both was 56 meV, a very small value in terms of bandgaps of semiconductors in the range between approximately 2 and 4 eV. This result suggests that scaling up to even larger unit cells should also yield consistent results.
Regarding the model, because of the limited number of structures considered relevant to describe this family of compounds with delafossite structure, high-complexity models like neural networks, which require large training sets, were automatically discarded even for initial evaluation. Beyond those, other smaller models have been tested, like simple regularized and non-regularized linear regressors, support vector regressors, and gradient-boosting regressors based on decision trees. For all of them, a proper search of optimal hyperparameters has been performed. The best results were achieved with a LASSO regressor (linear regressor with L1 regularization) with an alpha parameter for the regularization of 0.01.
For the evaluation of the models and descriptors, a standard train–test splitting was implemented within a leave-one-group-out cross-validation scheme [
73], using the 64 structures of 12 atoms described in
Section 5. The size of the test is equal to the square root of the total number of compounds in all the cases, with the training set made up of the remaining compounds. Groups for cross-validation were chosen so the balance between training and test sets was always maintained, meaning that the number of compounds in the set with any of the different metal cations studied was always the same. Following that scheme, average and best results between the series of folds in the cross-validation are here presented in
Table 10, which shows main performance evaluation metrics like the mean absolute error (MAE) and maximum error (
) between the predicted and DFT calculated results. According to common procedures of model evaluation, the errors presented in the table have been calculated only for the test sets. The results for the best fold are also presented in
Figure 13a. Given that all the Fe compounds present a GGA bandgap of 0, an analogous study was performed without including these compounds. Results for the sets without Fe compounds are also presented in the table, for comparison purposes, as it is clearly seen how the performance metrics improve filtering them.
Average results show that this method yields approximate results, which certainly should not be taken as definitive nor exact. However, the aim of this process is not to obtain the most accurate bandgap values but to obtain results within an acceptable range of accuracy that allow us to preselect which compositions would be the most suitable for our purposes.
With the choice of the descriptor and the model already defined, the same process has been followed to evaluate the performance of new structures, with the same crystal structure and same set of species in N, M positions but with relative concentration ratios of metal cations not used for the training of the model. A set of eight new structures (CuAl
1/6Ga
5/6O
2, CuAl
1/2Ga
1/2O
2, CuAl
5/6Ga
1/6O
2, CuAl
5/6In
1/6O
2), CuAl
5/6Cr
1/6O
2, CuAl
1/6Co
5/6O
2, CuGa
1/2In
1/2O
2, CuSc
5/6Y
1/6O
2) have been used as the test set, using the sixty-four original compounds as the training set in this case. These results, with the corresponding performance metrics (MAE,
and r2), are presented in
Figure 13b. Slightly higher errors are obtained in this case, compared with those obtained for the original compounds; however, they are very similar, considering that we are predicting concentration ratios not used to train the model. These results, with errors that fall within a tolerance margin of +/−0.2 eV, could still be of perfect use for preselecting compounds based on this property with high levels of confidence.
These results prove that this machine learning-based approach could be extremely useful to further explore this family of compounds efficiently and affordably, especially considering the titanic computational cost of running HSE calculations for several of these compounds that requires a description with large unit cells and the very inaccurate results obtained with GGA functionals. Indeed, a similar approach could be used as well to predict other properties of interest to sketch a map of properties to efficiently select the most promising compounds within the compositional landscape of this family of compounds, given those properties can be directly inferred from DFT calculations like the ones performed in this work.