Al$_x$Ga$_{1-x}$As crystals with direct 2 eV band gaps from computational alchemy

We use alchemical first order derivatives for the rapid yet robust prediction of band structures. The power of the approach is demonstrated for the design challenge of finding Al$_x$Ga$_{1-x}$As semiconductor alloys with large direct band gap using computational alchemy within a genetic algorithm. Dozens of crystal polymorphs are identified for $x>\frac{2}{3}$ with direct band gaps larger than 2\:eV according to HSE approximated density functional theory. Based on a single generalized gradient approximated density functional theory band structure calculation of pure GaAs we observe convergence after visiting only $\sim$800 crystal candidates. The general applicability of alchemical gradients is demonstrated for band structure estimates in III-V and IV-IV crystals as well as for H$_2$ uptake in Sr and Ca-alanate crystals.


I. INTRODUCTION
Major advancements of many technologies hinge upon the performance of underlying materials. The efficiency of photovoltaic materials, for example, heavily depends on the electronic properties of the solid-state or organic semiconductors, or the stability of energy storage materials is governed by their energy landscapes and vibrational modes. With the help of modern compute hardware and ab initio theories, properties of materials can nowadays be predicted in silico before their assessment is carried out experimentally. 1-7 Subsequent synthesis and characterization is costly and should focus only on those materials predicted to exhibit the most promising prop- The inset exemplifies the mating of randomly picked parent A and B which produce a top-performing child C in the third generation. C subsequently enters the updated pool of parents. This process is repeated until convergence, based on rankings exclusively obtained through alchemical perturbation estimates. Initialization corresponds to a single unperturbed self-consistent DFT calculation and 20 random perturbations. (b) Self-consistent DFT validation is carried out for all converged parents, leading to a slightly altered ranking.
erties. Computational materials design procedures therefore hold much promise to lift some of the chemical challenges the world is facing today.
Finding materials with superior properties in chemical compound space (CCS), the space consisting of all possible materials, 8,9 can be seen as a numerical optimization problem. Due to the immense amount of possible compounds and atomic configurations in CCS naive screening is prohibitive. The combinatorial nature of CCS implies that fast yet accurate property estimates, as well as efficient algorithms for searching, are desirable. 10-13 Alchemical derivatives provide useful, chemically local, gradient information, suggesting the applicability of gradient-based algorithms to the design new materials. Alchemical derivatives were already employed to model the stability of binary solid mixtures within virtual crystal approximations 14,15 . They have also afforded accurate predictions of various properties and compound classes within less approximate perturbation theory [16][17][18][19][20][21][22][23] . Here, we present numerical evidence which suggests that alchemical first order derivatives can also be used for the rapid yet robust prediction of band structures. For molecules and ionic crystals, even chemical accuracy can be achieved in terms of alchemical first order based estimates of relative energies when fixing number of valence electrons and geometry. 24,25 Such constraints can also be met by several material classes of great interest, such as III-V semiconductors. Here, we rely on a materials design algorithm which combines alchemical gradients with stochastic sampling and which holds promise for general computational materials design campaigns. Using this algorithm, we have optimized III-V solid alloy solution mixtures consisting of Al x Ga 1−x As with respect to band structure, arguably one of the most important properties of semiconductors due to their many electronic applications. [26][27][28][29] The algorithm is illustrated in Fig. 1 where |Ψ ref is the wavefunction of the reference system. 9,18,24 To estimate changes in band structure, the same formula is applied to each eigenvalue ofĤ tar at any given wavevector k using the corresponding eigenfunction |φ ref k . In Fig. 2, the predictive performance for band structures is illustrated using GaAs (AlAs) as a reference in order to predict AlAs (GaAs): Note how most of the details in the band structure are faithfully reproduced, and how the estimates can even account for the change from direct band gap with conduction band minimum at Γ for GaAs to indirect band gap with minimum at X for AlAs. The prediction error ∆ε LUMO (Γ) = ε pred LUMO (Γ) − ε tar LUMO (Γ) amounts to 0.25 (-0.2) eV, and ∆ε LUMO (X) = 0.14 (-0.15) eV for AlAs (GaAs) using GaAs (AlAs) as a reference. This implies a low band gap prediction error of only ∆E g = 0.14 eV and ∆E g =-0.2 eV for AlAs and GaAs, respectively. A scatter plot for estimates corresponding to every k-point and every band is shown in Fig. 2(c), indicating a very good correlation. The mean absolute error (MAE) made when predicting eigenvalues for all occupied bands and k-points is as low as 0.11 eV.
Such predictive accuracy is no coincidence. Averaged MAEs for alchemical band structure predictions among all possible III-V and IV-IV semiconductors are provided in Fig. 3, Not surprisingly, predictions "close-by", e.g. among III elements for the same V element, are generally very accurate, and largest errors are found for combinations associated with dramatic changes in electronic states, e.g. estimates of III-Sb or Sn crystals using III-P or Si as respective reference. Integrated changes in electron density, shown in Fig. 3, are found to correspond to an upper error bound, suggesting the possibility to construct meaningful error measures by means of inexpensive approximate estimates of electron density changes.

III. COMPUTATIONAL DETAILS
Alchemical derivative and PBE results are computed using the ABINIT package 30 with Goedecker normconserving pseudopotentials 31,32 and planewave cutoff of 100 Ry. Monkhorst-Park 33 6 × 6 × 6 mesh for fcc primitive cell band structure, and 3 × 3 × 3 mesh for band structure optimization (both for 3 × 3 × 3 and 4 × 4 × 4 fcc supercell) are used respectively. VASP 34 package with PAW 35 pseudopotentials with default cutoff are used for HSE 36 results for Al x Ga 1−x As 3 × 3 × 3 fcc supercell. 10 and 16 virtual orbitals are used for fcc primitive and supercells respectively. Experimental lattice constants (if available) are used for every crystal considered in this report. Otherwise calculated values are used from the literature (see table I). Alchemical derivatives are printed by restarting with the reference wavefunction in ABINIT, and setting the SCF iteration step to 0 to evaluate ∂ λ E. Three independent optimizations have been carried out for 3 × 3 × 3 fcc supercells, and one optimization for a 4 × 4 × 4 fcc supercell. Gradient sampling is carried out with population size of 20 and mutation rate of 0.5%. Pseudocode and convergence criteria have been included in the appendix. The 3 × 3 × 3 fcc supercell crystal is extended to 9 × 9 × 9 fcc supercell to compute the peri-odic Coulomb matrix in order to account for all possible periodic images.

A. Band structures of III-V and IV-IV Semiconductors
Alchemical prediction from the reference material of lattice constant a to any target material (tar) on the i th band at a given crystal momentum k is constructed as follows: , calculated from the i th orbital of the reference material of lattice constant a with crystal momentum k. The prediction of band gap is defined as: where the maximum of conduction band of materials investigated is located at Γ-point. A predicted band gap is direct if the predicted minimum of ε LUMO is located at Γ-point. Lattice constants for the semiconductors are scanned from 5.4Å to 6.5Å in step of 0.05Å including extra points correspond to the lattice constant reported in the literature. This range covers all equilibrium constants reported for all the semiconductors considered. We quantify the performance of our predictions using the mean absolute error (MAE) of a prediction to the entire i th band of a target material at lattice constant a, made from some reference material as where ε i,a (k) is the i th eigenvalue of Hamiltonian at lattice constant a and crystal momentum k. w(k) is the corresponding Monkhorst-Pack 33 weight for the sampled special k-points. Throughout this paper, the error reported corresponds to average over all bands and lattice constants. The subscripts i and a are omitted unless noted otherwise. The averaged MAE over all lattice constants and bands for all combinations of reference/target pair, or alchemical path, is on display in the upper panels in Fig. 3. Overall, decent agreements are found with most predictions being in agreement with the target by less than 0.5 eV. Alchemical estimates with averaged MAE less than 0.2 eV are highlighted by the white circles. A general trend can be observed: Predictions using semiconductors containing 3 rd -row elements in the reference give overall the largest averaged MAE, when the target material has elements from the 5 th row, as can be seen by the red/orange corner of the upper panels (left for III-V semiconductors and right for IV-IV semiconductors) in Fig. 3. A similar target-reference pattern has been observed for alchemical predictions of covalent bonding, and is due to lack of similarity of electron densities. 24,25 A scatter plot between averaged MAE and integrated absolute electron density difference, defined as |∆ρ(r)| = |ρ tar (r) − ρ ref (r)|, is shown in the lower left panel in Fig. 3. The results suggest that there is an upper error bound: Any alchemical path with small density changes will give a small predictive error to band structure, as all the points lie beneath the red dotted line. However, a small error does not necessarily imply small density changes, as many results with small predictive error correspond to large density changes. This is because the band structure is determined by the relative differences in Hamiltonian eigenvalues and the corresponding orbital structure. Cancellation of higher order (curvature) effects along the alchemical path can lead to a small predictive error for first order estimates yet large density change. Similar error cancellation has also been discussed in the context of covalent bonding Ref. 24.
The fact that small density changes imply accurate predictions can be used as a sufficient condition to detect good predictive power. This might be useful for future studies if decent approximations to inexpensively estimate electron density changes can be found. The region with dr|∆ρ(r)| < 0.5 a.u. electron/primitive cell is highlighted by the white background in the lower left panel. The band gap predictions of the corresponding alchemical paths highlighted as white crosses in the upper panels of Fig. 3. A scatter plot versus the true band gap is shown in the lower right panel, remarkable agreement is found with MAE of less than 0.05 eV. The linear fit gives E true g ≈ 0.93E pred + 0.019 eV, with R 2 = 0.93, RMSE = 0.047 eV and MAE = 0.036 eV. The negative E g in the lower right panel of Fig. 3 correspond to GaSb and InSb with small lattice constants. In such environment, the conduction band minimum is lower than the valence band maximum. Notice that decent predictive power can be achieved by alchemical predictions, satisfying the sufficient condition above, even at such extreme situation.
It is possible to use III-V semiconductor reference calculations to predict the band structure of other IV-IV semiconductors. However, such interpolations do not provide satisfactory predictive power when using only first order alchemical derivatives, due to the dissimilarity of the electronic density. Similar conclusions hold when attempting to predict band structures of II-VI semiconductors using III-V reference densities.

B. Band-gap maximization in AlxGa1−xAs
Numerical results with such predictive power suggest that first order alchemical estimates are sufficiently accurate for robust yet efficient optimizations of Al x Ga 1−x As based band structures. Unfortunately, due to the large dimensionality of the problem, it is still challenging to op- timize the band structure using reliable gradients alone. For example, modeling a 3x3x3 supercell of pure GaAs, corresponding to 27 Ga atoms, implies a total compositional space of 2 27 ∼ 134 million possible Al x Ga 1−x As combinations (not accounting for symmetry). Therefore, we have used the genetic optimization algorithm, described in Fig. 1 and based on alchemical perturbations towards the pseudopotential of Al at all Ga sites. Starting with a single unperturbed reference PBE DFT calculation, the corresponding direct band gap maximization history, on display in Fig. 4, indicates rapid convergence towards predicted PBE band gaps of ∼1.6 eV after less than 800 generations. The optimization history for each generation (Fig. 4(a)) clearly illustrates how the average trend of direct E pred g is moving upward as the gradient sampling iterations proceed These results imply that the mating procedure during hybrid optimization Fig. 1(a) systematically steers the population towards larger direct E g . Since the band structure is determined by the structure of occupied and unoccupied orbitals, this also indicates that crystal truncation and catenation roughly preserve the local structure of orbitals around each atom. As a result, the algorithm identifies alloys with E pred g of around 1.6 eV after only several hundred iterations. Among the best ten alloys out of 1444 identified (table II) within hybrid optimization, a crystal Al 0.67 Ga 0.33 As (structure shown in Fig. 5), with direct E tar g = 1.4 eV has been selected, and its unfolded 41,42 target band structure is plotted in Fig. 4(c), where the spectral weight at Γ is 70.7%, and the band gap prediction error |∆E g | = 0.18 eV. We also note that the average value of the top 20 candidates converges quickly within ∼1000 generations.
For the select examples (best, good, bad, worst) a decent correlation between E pred g and E tar g is shown in Fig. 4(b) with a linear fit yielding MAE = 0.014 eV. Since the trends are preserved between E tar g and E pred g , one genetic optimization is sufficient to identify the global optima. Optimization performance obtained for other starting structures or larger super cells corroborates this observation, indicating significant robustness and generality of the design approach. Consecutive optimizations confirm that the best alloys identified during each hybrid optimization are structurally similar see Fig. 5(b) . HSE band gaps are known to be in better agreement with experiments than PBE. 43,44 . We have therefore calculated the corresponding HSE band gaps 45   the PBE weights for unfolding, assuming that they will negligibly affect the HSE gap. Throughout the optimization history, E pred g values cover a range of 0.9 eV to 1.6 eV, as plotted as a sorted sequence in Fig. 5(a). The cross-over from direct to indirect gap occurs, in line with previous calculations. [46][47][48] , for mole fractions exceeding x ≈ 0.7, substantially larger than what has been realized in experiments so far, i.e. x ≈ 0.4.
To better understand what and how structural features impact changes in direct band gap, we have analyzed the crystal structures visited throughout the alloy optimization history. We compare alloys with large and small band-gap using a periodic variant of the sorted Coulomb matrix representation C 49,50 , where with nuclear charges Z and minimal distances min(R IJ ) between periodic images of atom positions R I and R J . C represents any crystal in a unique fashion and allows us to quantify the similarity between crystal C A and C B by the normalized matrix norm ||C A − C B ||/||C B ||. Rank- (a) Normalized sorted coulomb matrix based distance distribution between best crystal (shown as inset with Al, Ga, As in blue, green, and red, respectively) and crystals from three colored bands highlighted in (a) (200 crystals from each band).
ing as well as distributions of similarities to the top alloy (crystal structure shown as inset) with the widest alchemically predicted direct band gap are presented in Fig. 5 for all alloys falling into the 2-200 (blue), 500-700 (green), and 1000-1200 (red) windows of candidate rank. Overall, and not surprisingly, the top alloys clearly indicate higher probability of being more similar to the best candidate (blue bars larger than green than red for decreasing dissimilarity values). The histogram even suggests that for similarity values smaller than 0.45, the sorted periodic Coulomb matrix alone can be used as a descriptor to identify large direct band gap crystals with high probability. Encouragingly, the relative mutual similarity distributions for intermediate (green) and distant (red)  crystal ranks also coincide with the corresponding trends in band-gaps.

C. H2 uptake in Ca and Sr alanates
While large band-gap semi-conductors are of widespread interest, one might wonder if the alchemical approach is also applicable to other materials design problem. In order to explore the general applicability of alchemical derivatives for ranking crystals, we now also consider the alchemically predicted ranking among Ca and Sr alanates with respect to their potential for H 2 uptake. Ca(AlH 4 ) 2 and Sr(AlH 4 ) 2 release hydrogen when heated, and have been proposed as reversible hydrogen storage materials. 51,52 We have investigated 42 Ca-and 44 Sr-alanate crystal polymorph structures from Ref. 53. We comply with the aforementioned constraint of fixed structure by generating all Sr-alanate (Ca-alanate) crystals with identical crystal structures for each Ca-alanate (Sr-alanate) crystal. The H 2 absorption energy of a Ca-alanateĤ tar is alchemically estimated using the Sr-alanateĤ ref in the same crystal structure. Results in Fig. 6 indicate that the alchemical estimates correctly predict the energetically strongest Ca-alanate, and the three strongest Sr-alanates. Decent predictive power regarding the H 2 absorption energy ordering is also illustrated in Fig. 6(c) and (d) for the 60 most stable alanate crystals with a linear fit corresponding to a MAE of only 0.02 eV. The ranking error distribution Fig. 6(d) inset and a Spearman's rank correlation coefficient of 0.96 also indicate substantial predictive power, sufficient for most materials design purposes.

V. CONCLUSIONS
To conclude and summarize, we have presented numerical evidence for the usefulness of alchemical first order derivatives for the prediction of band structures. To illustrate the power of this method, we have optimized band structures in III-V semiconductors, namely to maximize the direct band gap in Al x Ga 1−x As semiconductor alloys. Thanks to the computational efficiency of a gradient based genetic optimization algorithm, multiple Al x Ga 1−x As alloys with direct E g ≈ 2.1 eV and mole fraction x ≈ 0.7 have been identified. We note that after having identified such promising materials candidates, the synthetic procedure to realize them in an experiment remains a largely outstanding challenge. The qualitative identification of alloys with the largest E g within single optimization runs implies that the E g surface is fairly flat in the crystal space of Al x Ga 1−x As systems. Alchemical derivatives can also tackle other challenges in the realm of computational materials design, such as estimating H 2 absorption energy ordering across Ca and Sr-alanates of varying crystal structure. Future work will deal with (i) extensions to higher orders, and (ii) systematic comparisons to alternative materials design approaches such as special quasirandom structures 54 , cluster expansion 55 , or machine learning approaches 56 . The objective of the gradient-based genetic algorithm (GA) is to maximize the band gap while the bottom of the conduction band must be located at Γ-point in the Brillouin zone. To this end, we vary the mixing ratio of Al and Ga, as well as their location. Instead of performing a DFT evaluation for each candidate crystal, the band structure of a candidate crystal is estimated using alchemical derivatives. Each alchemical prediction takes less than 1% of the computational cost for a full DFT  band structure calculation. This enables the rapid exploration of many configurations and mixture ratios with decent accuracy. The Brillouin zone of the supercell is unfolded and corresponds to the one from a primitive cell using spectral weights. 41,42 Three independent gradient-based GA optimizations have been carried out, start with three starting crystals: Pure GaAs and two random crystals, where either Ga or Al is randomly chosen for each of 27 III-sites. All three optimizations converged within two optimization steps, which give total six GA sampling histories. 3137 crystals have been evaluated in the first GA sampling history. Out of these, 1444 are predicted to be direct band gap crystal. Roughly half of the searched crystals are of indirect band gap. E true g are calculated for the following four groups of crystals from the first GA sampling step history, depending on the sorted direct band gap predictions: the best 50 E best g , the top 491 to 504 E good g , the top 1001 to 1010 E bad g , and the worst 15 E worst g . The analysis of the first GA sampling history is shown in Fig. 4, while the analysis of all three optimization histories is presented in Fig. 5. Pseudocode for the gradient-based optimization is presented in the pseudocode below, where each of the routines corresponds to the following.
• optimization: The main routine consist of: A stochastic sampling procedure; and a gradient step that updates wavefunction by full DFT evaluation, as schematically shown in the left panel of Fig. 1.
.001 eV at h th optimization step.
• full DFT: It performs full DFT evaluations on the requested crystals and returns the best corresponding true band gap, E true g , and orbitals, {φ i }.
• GA sampling: The routine stochastically samples the alchemical estimates, using reference orbitals {φ i }. The initial population of twenty crystals is randomly generated. The best five crystals are used for full DFT evaluation. The sampling criterion is at least 1400 crystals with direct band gap are found. And the difference between the average of the top-20 population,Ē pred g , and the best-predicted crystal, E best g , is less than 0.02 eV.
• get parents: At each iteration during GA sampling, two parents, ParentA and ParentB, are randomly drawn from the top-20 of all searched crystals. In other words, the population size of GA is twenty.
• mate: The routine to generate a child from two parent crystals. The occupancy at each of the IIIsites in the child crystal is randomly inherited from either parent with 50/50 possibility, as illustrated in the right panel of Fig. 1. A mutation rate of 0.05% is used at each of the III-sites, where the atom type would flip from Al to Ga or vise versa as indicated the yellow atom in the right panel of Fig. 1.
• alchemical evaluation estimate the band structure of the requested crystal using reference orbital {φ i }. If the prediction of the band gap is indirect, the value of the direct band gap is set to 0. * Electronic address: anatole.vonlilienfeld@unibas.ch 1 F. Capasso, "Band-gap engineering: From physics and materials to new semiconductor devices," Science, vol. 235, p. 172, 1987.  Table. II. The HSE ordering rank and the crystal formula with mole fraction is shown for each crystal. The 3 × 3 × 3 supercell is extended to 6×6×6 to illustrate 2-fold symmetry, where Al, Ga, and As are represented by blue, green, and red spheres.   Table I in main text). Each column corresponds to a crystal, each row to the atomic site. Upper and lower block correspond to III and V elements, respectively Al Al Ga Ga Ga Al Al Ga Al Al (16.0, 9.2, 6.5) 28 As As As As As As As As As As ( 2.0, 1.2, 0.8) 29 As As As As As As As As As As ( 6.0, 1.2, 0.8) 30 As As As As As As As As As As (10.0, 1.2, 0.8) 31 As As As As As As As As As As ( 4.0, 4.6, 0.8) 32 As As As As As As As As As As ( 8.0, 4.6, 0.8) 33 As As As As As As As As As As (12.0, 4.6, 0.8) 34 As As As As As As As As As As ( 6.0, 8.1, 0.8) 35 As As As As As As As As As As (10.0, 8.1, 0.8) 36 As As As As As As As As As As (14.0, 8.1, 0.8) 37 As As As As As As As As As As ( 4.0, 2.3, 4.1) 38 As As As As As As As As As As ( 8.0, 2.3, 4.1) 39 As As As As As As As As As As (12.0, 2.3, 4.1) 40 As As As As As As As As As As ( 6.0, 5.8, 4.1) 41 As As As As As As As As As As (10.0, 5.8, 4.1) 42 As As As As As As As As As As (14.0, 5.8, 4.1) 43 As As As As As As As As As As ( 8.0, 9.2, 4.1) 44 As As As As As As As As As As (12.0, 9.2, 4.1) 45 As As As As As As As As As As (16.0, 9.2, 4.1) 46 As As As As As As As As As As ( 6.0, 3.5, 7.3) 47 As As As As As As As As As As (10.0, 3.5, 7.3) 48 As As As As As As As As As As (14.0, 3.5, 7.3) 49 As As As As As As As As As As ( 8.0, 6.9, 7.3) 50 As As As As As As As As As As (12.0, 6.9, 7.3) 51 As As As As As As As As As As (16.0, 6.9, 7.3) 52 As As As As As As As As As As (10.0, 10.4, 7.3) 53 As As As As As As As As As As (14.0, 10.4, 7.3) 54 As As As As As As As As As As (18.0, 10.4, 7.3)