Assessment of interstitial potentials for rapid prediction of absolute band energies in crystals

Photoelectrochemical cells (PECs) can be applied to photogenerated charges to “split” water molecules into hydrogen and oxygen, which can then be subsequently burnt as a renewable, carbon free fuel.¹ While many device architectures exist, they all require an appropriate alignment of the electronic band edges of the photoactive material with the hydrogen and oxygen redox potentials.¹ To determine whether a material has such an alignment, the ionization potential (IP) and electron affinity (EA) must be known. In semiconductors, the IP and EA correspond to the absolute positions of the valence band maximum (VBM) and the conduction band minimum (CBM) with respect to an external reference.² In contrast, the work function is defined from the Fermi level, which typically falls between the VBM and CBM.

There exist several methods to extract values for the IP and EA experimentally, the most common of which require the energy of a transferred electron to be compared to a standard reference potential. This is normally the vacuum level for photoemission based techniques, such as photoelectron spectroscopy,^3,4 or the standard hydrogen electrode (SHE) for electrochemical techniques, such as cyclic voltammetry.² The measured IP and EA are influenced by the bulk and surface properties of the material, and this, in turn, affects the consistency of reported values.^5,6 The surface dependence can make band alignment measurements hard to compare, as it is difficult to source data that are both consistent and comprehensive. With this in mind, computational screening of the IP and EA would be desirable, given its high automatability and expediency. However, the first challenge lies in the choice and definition of the aforementioned reference potentials.

For crystalline systems, the total energy and electronic eigenstates are efficiently calculated with the use of periodic boundary conditions to model an infinite single crystal. However, the eigenstates alone are not calibrated with respect to an external potential, such as vacuum.⁷ In many implementations of density functional theory (DFT) for crystals, the average electrostatic potential of a given unit cell is used as an internal reference and set to 0 V. Consequentially, only relative quantities, including the bandgap, which are not dependent on an external reference potential, are valid for experimental comparison.

In practice, this limitation is often overcome in one of two ways. The first involves retaining the periodic boundary conditions but purposefully creating a supercell slab model containing the bulk, surface, and vacuum regions of the crystal. The electronic eigenvalues can then be properly referenced by taking the numerical value of the electrostatic potential in the local vacuum region away from the calculated eigenvalues.⁸ This approach calibrates with respect to a local vacuum level that depends on the crystal orientation/termination. The second method involves building a cluster model, embedded in the potential of the periodic crystal,⁹ which includes a vacuum region in its construction. In both cases, the computational overhead is large enough to make a high throughput materials screening procedure impractical. Other approaches developed to address this problem include comparison of charge neutrality,^10,11 impurity,¹² core,¹³ or defect transition¹⁴ levels between different materials.

In this study, we assess the use of an alternative reference potential based on the bulk electrostatic potential. Such an approach was recently applied in the context of porous solids³ but can be traced back to the work by Frensley and Kroemer.^15–18 For this method, the average electrostatic potential of the crystal interstices selected by the crystal symmetry is used as an effective local vacuum reference. We revisit the case of zinc blende (ZB) semiconductors, but the approach is more general and adaptable to other crystal classes, which we demonstrate for rock salt and perovskites.

The zinc blende structure is composed of two interpenetrating tetrahedral sublattices. Between these two sublattices (A and B), charge transfer is expected due to differences in electronegativity, resulting in a more positively charged sublattice (A^δ+) and negatively charged sublattice (B^δ−). In the complete zinc blende structure (Fig. 1), there are two symmetry unique interstices: the octahedral site (O_h) and the tetrahedral site (T_d). The electrostatic potential of the O_h and T_d sites are influenced in equal measure by each of the two sublattices. If O_h is most closely surrounded by ions from sublattice A, then T_d is surrounded by ions from sublattice B. This structure produces the two distinct maxima in the electrostatic potential for O_h and T_d. By symmetry, O_h in sublattice A is identical to T_d in sublattice B and vice versa. Thus, depending on how the unit cell basis is defined, O_h and T_d could correspond to either of the electrostatic potential maxima. However, their average (ϕ_ZB) is single valued and 0 V in the case of a lattice of point charges,

The concept described can be verified numerically. Treating the atoms as point charges and using a suitable Ewald summation technique, such as in the General Utility Lattice Program (GULP),²⁰ probe atoms can be placed at each of the interstices, and the corresponding Madelung potentials are calculated. The O_h and T_d sites in zinc blende possess antisymmetric Madelung potentials, making their average equal to 0 V, which is independent of the magnitude of the atomic charges in this binary system.

By utilizing ϕ as a bulk electrostatic potential reference, the Frensley–Kroemer methodology has a potential advantage over the more conventional methods by only requiring a calculation of the bulk electrostatic potential and corresponding electronic eigenvalues. Once ϕ is known, the bulk VBM can be adjusted to the absolute value according to

In principle, this would require far less setup and computational time per calculation, a seemingly ideal candidate for materials screening. The two key assumptions being made are that the bulk and surface components of the IP are separable and that ϕ is a reasonable approximation to the vacuum potential.

Indeed, the validity of this approach has been a subject of debate, with Van de Walle and Martin finding poor agreement with their results using model solid theory.²¹ They highlighted pitfalls in assuming a negligible interstitial charge density and flat interstitial electrostatic potential. However, a comparison of band offset models by Yu et al. found that the performance of the Frensley–Kroemer method remained on par with other procedures of the time.²² Later, this was re-iterated by Tung who compiled the data from other charge density based IP calculation methods.²³ These results are plotted in Fig. 2, alongside a list of photoemission measurements compiled by Pelatt et al.⁴ 

Two recent studies employed similar methodologies to the Frensley–Kroemer method based on DFT calculations. Choe et al. implemented an internal reference method of band alignment but with a focus on metal–metal and metal semiconductor interfaces,²⁴ whereas Butler et al. used the centers of metal–organic framework (MOF) pores as the vacuum potential reference.³ We note that with modern day increases in computational power, higher resolution fast Fourier transform (FFT) grids enable a more precise interrogation of the electrostatic potential in and around the interstitial region.

Due to the crystal symmetry of zinc blende, the potentials at the O_h and T_d sites are equal and opposite. Hence, the magnitude of their difference may give insight into the effective charges of A^δ+ and B^δ−. Frensley and Kroemer attempted to exploit the interstitial potential difference ΔV = V(O_h) − V(T_d) to calculate the effective charges.¹⁷ This offers an alternative perspective to electron density and wavefunction based partial charge analysis for atoms in compounds.²⁵ 

If one considers two zinc blende crystals—one more ionic and another more covalent—the ionic crystal will have a greater value of ΔV as the potential difference is greater in the more ionic case. Indeed, for a covalent crystal with the diamond structure, such as Si and Ge, ΔV = 0 by virtue of symmetry.²⁶ For a binary ZB compound, the partial atomic charge (δ) can be calculated as follows:¹⁷ 

This adjustment normalizes ΔV with respect to the lattice parameter a, where 8.827 is a constant derived from the crystal topology. Since absolute charges are limited in practical use,²⁷ further normalization with respect to the formal ionicity has been performed using feature scaling for each compound,

where I′ is the normalized ionicity, I, while I_max and I_min are the maximal and minimal oxidation states. For all compounds studied here, I_min = 0; thus, the above expression reduces to

For example, in ZnS, the formal oxidation state is Zn(II). If the partial charge was found to be +1, then I′ = 0.5.

A set of 27 zinc blende type compounds was investigated for the local-density approximation (LDA),²⁸ solid state Perdew-Burke-Ernzerhof (PBEsol),²⁹ strongly constrained and appropriately normed semilocal (SCAN),³⁰ and Heyd-Scuseria-Ernzerhof (HSE06)³¹ functionals. The diamond structured compounds Si and Ge were also tested. The O_h and T_d eigenvalues and electrostatic interstitial potentials were calculated using Density Functional Theory (DFT) within the Vienna Ab Initio Software Package (VASP)^32–34 with Projector Augmented Wave (PAW) frozen-core potentials.^35,36 For all calculations, a 13 × 13 × 13 k-point gamma-centered mesh was utilized with a plane wave cutoff energy of 500 eV. The tetrahedron method with Blöchl corrections³⁷ was the chosen method for Brillouin zone sampling.

The ground-state lattice parameters were calculated via a primitive unit cell volume scan fitted to the Murnaghan equation of state.³⁸ Using the optimized lattice parameters, the electrostatic Hartree potential, charge density, and eigenvalues were extracted using an enhanced 170 × 170 × 170 real-space FFT grid. All calculations were performed using the zinc blende primitive cell, unless otherwise stated.

The Macrodensity code³⁹ was utilized to extract the Hartree potential at the interstices. The calculation of the absolute VBM was performed by taking the average interstitial Hartree potential from the highest occupied Kohn–Sham eigenvalue. The calculation of the absolute CBM was performed in the same fashion but with the lowest unoccupied eigenvalue. The Bader effective atomic charges were obtained from the VASP charge density output using the Henkelman code.⁴⁰ 

To apply this alignment procedure, there must be a plateau in the electrostatic potential at the interstitial site of interest. The flatness of the potential at each interstice was investigated by increasing the sampling volume, which is averaged to give ϕ. For clarity, we outline two terms used in this work: “FFT grid” and “sampling volume.” The FFT grid is the grid density used in the underlying DFT calculation and sets the resolution of the electrostatic potential across the unit cell. The sampling volume is a sample cube of points taken from the FFT grid about the interstitial region (see Fig. S1 in the supplementary material for further details). Thus, in the case of ZB, two sampling volumes are extracted from each FFT grid. The VBM calculated for each sampling volume is given in Fig. 3. The PBEsol functional and an enhanced 300 × 300 × 300 FFT grid were used here; however, the same qualitative conclusions were reached for testing LDA and PBE despite some small quantitative differences. 100 × 100 × 100 and 50 × 50 × 50 grids were tested with the same sampling volume sizes.

As the FFT grid becomes more coarse, adjacent points are spaced further apart, and sampling of the electrostatic potential will extend over a greater real space volume. The variance of V(T_d) was significantly greater than V(O_h), often by 2 orders of magnitude for the 10 × 10 × 10 sampling volume size. Geometrically, the O_h volume is expected to be larger than the T_d volume. This observation supports the notion that the interstitial volume, for any given FFT grid, is a significant factor. For convergence testing before a production run, we recommend using single point sampling (a 1 × 1 × 1 sampling volume) at each interstice and converging with respect to the FFT grid density. However, this does not overcome inaccuracies inherent to the method itself, such as the interstitial volume being too small. The use of progressively larger sampling volumes allows for a quick evaluation of the local curvature about the interstitial region. A deviation from pointwise sampling implies a greater curvature and possibly less accuracy, which may be what we see for the case of BN.

An FFT sampling volume of 100 × 100 × 100 was found to be sufficiently converged for our band alignment case study. However, the Bader charge analyses required finer resolutions of the electron density, which is why 170 × 170 × 170 was used as standard in the rest of this work.

The convergence tests suggest classes of materials that may be ill-suited for this method. In Fig. 3, the greatest variation across different sampling volume sizes was for the boron compounds, which also had the smallest interstitial volumes, by virtue of their small lattice parameters. On the other hand, compounds with comparably small lattice parameters—such as SiC—varied much less between sampling volume sizes. This implies that the observed variability of the boron compounds may be due to something other than a small interstitial volume. This is reflected in the top half of Fig. 4, which gives the VBM difference between the 10 × 10 × 10 and single point sampling volume sizes plotted against the lattice parameter of each compound.

While the compounds with the smallest lattice parameters tend to produce the greatest divergence from single point sampling, there was a large degree of variation alongside this trend. This phenomenon is not expected to be a result of the DFT functional, as similar behavior is found in each case, but is related to the nature of the chemical bonding including the extent of wavefunction delocalization. Further evidence that the accuracy of the Frensley–Kroemer method is limited by the interstitial volume lies in the comparison of the calculated IP with reference data. The difference in the calculated VBM and IP reference data used in Fig. 2 is plotted against the lattice parameter of the compound in the lower half of Fig. 4. Here, we observe a similar trend as in the upper half, indicating that the small lattice parameter compounds are less accurate with regard to the IP themselves, even when a single point sampling volume size is used with a reasonably fine FFT grid.

In general, the choice of the exchange–correlation functional has little impact on the qualitative predictions. There is a systematic shift whereby the hybrid HSE06 functional predicts a deeper valence band energy, which is consistent with established behavior.^41,42 Figure 5 demonstrates that modern DFT calculations are able to reproduce Frensley and Kroemer’s original dataset. The most notable deviation from expected behavior is BN, which has the smallest lattice parameter in the dataset.

We note that BN is the only compound to exhibit a negative electron affinity (energy above vacuum), which is qualitatively consistent with reported UPS measurements. However, it was not conclusively determined whether this result was due to the surface or bulk properties, defects, or the mixing of cubic and hexagonal phases.⁴³ In contrast, other sources^4,44 report a positive electron affinity. Hence, there is more ambiguity concerning BN compared to the rest of our dataset.

The final band alignments, with an optimized FFT grid and sampling volume size, are displayed in Fig. 6. Upon observation of Figs. 2 and 5, it can be deduced that the HSE06 functional appears to perform better than the other functionals (excluding AlAs). Considering that the Frensley–Kroemer method overestimates the IP (see Fig. 2), this results in a slightly better quantitative agreement with the reference data compiled by Pelatt.⁴ As such, we expect that this effect is likely coincidental rather than physical. A comparison of each functional with the reference data from Fig. 2 is given in Fig. S3.

Additionally, the bandgap underestimation found for local and semi-local functionals (LDA, PBEsol, and SCAN) is corrected for HSE06 through a systematic lowering of the VBM and raising of the CBM. Since we calculate the electron affinity via EA = bandgap + IP, the value is dependent on the accuracy of the bandgap and the accuracy of the IP. For reference, the associated bandgap values are given in Fig. S2 for each functional.

The dataset includes species with formal oxidation states of 0 (e.g., Ge), +1 (e.g., Ag), +2 (e.g., Zn), +3 (e.g., Ga), and +4 (e.g., Si in SiC). A comparison between the derived Frensley–Kroemer ionicity and Bader effective ionicity is made in Fig. 7. There is little quantitative agreement between the metrics, but there is qualitative agreement within several of the periodic groups and significant differences between these groups. For example, both methods predict a similar ionicity trend across the In, Ga, Cd, Zn, and Al compounds. However, the two metrics disagree on where these groups are positioned on the ionicity scale. Both metrics are in agreement that Si and Ge are covalent.

Aside from the case of BAs, the ionicity derived from Bader partial charge analysis is greater for each compound compared to their Frensley–Kroemer counterpart, albeit to a varying degree. There is little agreement between the two metrics for the B, Cu, and Be compounds. When we compare the two effective ionic charge metrics with the measure of Phillips,²⁶ we find that there is better quantitative agreement between Frensley–Kroemer and Phillips for the Be and Al compounds and for SiC and BN (see the lower half of Fig. 7). However, it is possible that this agreement is coincidental as the overall agreement with the Phillips ionicity is poor for both the ionicity metrics used here. This is not surprising, given that the Bader and Frensley–Kroemer metrics have more in common with their interpretation of the ionicity—via effective ionic charges—than Phillips, which is based on a spectroscopic theory.

While the dataset considered here was based on zinc blende compounds, a similar analysis can be extended to other structure types. For example, the $(14,14,14)$ position in a cubic perovskite structure formally has an electrostatic potential of 0 V owing to its centrosymmetric environment. This relation has been confirmed numerically through an Ewald summation for a monopolar crystal. We have tested the approach for CsPbI₃, and it predicts a VBM of −5.13 eV (HSE06) below the vacuum level. This value is in the range typically found for lead iodide perovskites.^45–47 

In rock salt compounds, the $(14,14,14)$ interstitial position could be used as a valid 0 V reference based on the crystal symmetry. Again this was confirmed through a numerical simulation. We have tested the approach for MgO, and it predicts a VBM of −4.85 eV (HSE06). The value is reasonable but shallower than literature values (e.g., 4.8–7 eV with various treatments⁵). These preliminary values were obtained using a 100 × 100 × 100 FFT grid, with a plane wave cutoff energy of 500 eV and Γ-centered k-point grids (CsPbI₃: 4 × 4 × 4 and MgO: 9 × 9 × 9). The contour plots of the Hartree potential in these crystals are given in Fig. 8, which confirm that electrostatic plateaus are formed around the interstitial sites of interest. We also note that the ionicity metrics cannot be applied in these examples, as there is only one interstitial site that is sampled in each case.

We have assessed an approach to rapidly screen the absolute band energies of crystals. We found qualitative agreement with the original reports of Frensley and Kroemer and the available literature data. A notable exception in this regard was BN. This is likely due to the small interstitial volume of the interstices. Indeed, we observed an overall worsening of the precision and accuracy of the method for compounds with lattice parameters below 5 Å.

The relatively low computational overhead and simplicity of the alignment approach make it an attractive basis to build and compare low-cost band alignment models. Further exploration of interface dipole corrections for the IP,⁴⁸ the orientation dependence of the IP, and alternative pseudopotentials may be helpful in mitigating the inaccuracies and bridging the gap to what can be measured in experiments. The original model has its limitations, as we identify, but we expect that the interstitial potential may become a useful feature for constructing non-linear machine learning models of absolute band energies.

See the supplementary material for additional supporting data (pdf).

We acknowledge the contributions of Daniel W. Davies and Young-Kwang Jung for the compilation of the initial zinc blende structure set. We also acknowledge useful discussions with Keith T. Butler, Matthew Okenyi, and Matthias J. Golomb. This project was supported by the Plastic Electronics Centre for Doctoral Training, and calculations were facilitated by the Imperial College London Research Computing Service. We are grateful to the UK Materials and Molecular Modeling Hub for computational resources, which is partially funded by the EPSRC (Grant No. EP/T022213/1).

The workflows that support the findings of this study are integrated in the Macrodensity package available at http://doi.org/10.5281/zenodo.884521.