Introduction

In 2000, during his noble prize lecture Prof. Herbert Kroemer stated “The interface is the device”.1 He was referring to the phenomenal success in design and application of semiconductor heterojunction devices in microelectronics. Arguably, Kroemer’s statement has an ever increasing relevance across a range of technologies far beyond transistors, where heterojunctions found their initial success. As electronic devices move towards smaller scales, the interface area, as a proportion of the device, increases, growing the importance of interfaces. Moreover, researchers are increasingly finding that interfaces between materials represent a rich space for the exploration of exotic properties that are not present in bulk materials, such as two-dimensional electron gases (or liquids) and quantum topological states. It is clear that the importance of the interface can only grow with the evolution of modern technological applications.

The fields of energy conversion and storage are no exception to this general principle. Interfaces are important both as sites of device operation, for example splitting charges in solar cells,2 or catalytic activity in photo/electrocatalysts3 and equally as important considerations in terms of the mechanical stability and lifetime of systems. While ab initio, or first-principles, methods are becoming an increasingly routine tool for designing and screening bulk properties of materials in both batteries and solar cells, the methods for modelling interfaces in these systems quickly and reliably are considerably less developed. However, the general field of theoretical surface science has thrived significantly in the last few decades, via close collaborations with experimental investigations, which have yielded an “atomic” resolution of relevant properties.4

In this review article we start by considering what useful and important interfacial properties are accessible from modern electronic structure methods, particularly density functional theory (DFT). Having provided the background to calculating these important properties, we consider how they serve as parameters for understanding and predicting the performance of battery, photovoltaic (PV), and photocatalytic (PC) devices. We survey the state-of-the-art in the use of DFT in the aforementioned applications.

We find that these fields have developed sophisticated frameworks to a number of different problems, while sharing only a few approaches. We consider what the fields can learn from one another and lay out the challenges that need to be addressed for the future of predictive interface modelling in energy materials. We reiterate Kromer’s original message and assert that the future of the device is the interface.

Theory

In this section we introduce the important properties for interface design, which are are accessible from first-principles calculations. We have separated the properties of interest into two broad categories—thermodynamic and electron energies.

Theoretical predictions, particularly using DFT calculations, depend significantly on the functional used to describe the electronic exchange-correlation (XC) interactions. For example, thermodynamic properties, such as electrochemical voltages and stability windows, and surface energies are described accurately5 either with a generalised gradient approximation (GGA) functional6 or a GGA corrected by a Hubbard U term7 (i.e., GGA + U). On the other hand, electronic properties, particularly band gap and band alignment, typically require post-DFT8,9 level of theory, such as hybrid functionals,10,11 including a portion of exact exchange, or the GW approximation.12 Additionally, the predicted adsorption properties of liquids on materials (see the following section) depend appreciably on the type of description of van der Waals (vdW) interactions within DFT.13,14 Detailing the precise parameters required for accurate thermodynamic, electronic or adsorption property prediction is beyond the scope of this review, and the reader is referred to other work on this topic.15

Thermodynamic properties

Thermodynamic (meta-)stability provides the most fundamental criterion for predicting the applicability of a system. With this in mind, we begin by considering the thermodynamics of interface and surface formation. The process to create a surface from the bulk is represented schematically on the left hand side of Fig. 1.

Fig. 1
figure 1

Important interface properties. Left: the thermodynamic properties that can be calculated from first-principles. Right: the electron energy properties that can be calculated from first-principles

The surface free energy γsurface in (J m−2) defines the penalty to cleave the chemical bonds of a bulk material to form a surface at a particular Miller index identified by its h, k, l plane.

$$\gamma _{{\mathrm{surface}}} = \frac{1}{{2A}}\left[ {G_{{\mathrm{surface}}} - G_{{\mathrm{bulk}}} - \mathop {\sum}\limits_i^{{\mathrm{species}}} {\kern 1pt} {\mathrm{\Delta }}n_i\mu _i} \right].$$
(1)

where A is the surface area (in m2) and Gbulk is Gibbs energy of the reference bulk material. In practice, the surface is modelled by a slab model consisting of a periodic 2D “film” formed by a few atomic layers parallel to the hkl crystalline plane of interest. Therefore, Gsurface of Eq. (1) is the Gibbs energy of the slab and has to be converged with the thickness of the slab along the non-periodic direction, z. Gsurface and Gbulk are generally approximated by the respective computed internal energies (i.e., the total energies provided by DFT—Esurface and Ebulk), thus neglecting entropic effects (configurational and vibrational), the pV term and the zero-point energy.

In the case of off-stoichiometric surfaces—surfaces whose stoichiometry is not commensurate with the bulk—the final surface energy depends on the environment set by the chemical potential μi for species i and amounting to an off-stoichiometry of Δni. Δni is negative (positive) if species i is removed (added) to the surface. The selection of appropriate references μi is crucial to identify regions of the chemical space where specific surface terminations dominate over others, in turn impacting the particle morphologies of solids (see discussion below).16,17

Although γsurface is sensitive to the thickness of the slab used in calculations, Sun and Ceder18 demonstrated that γsurface converges much faster to accurate values and minimal thickness is required (thus reducing the computational resources) if the reference bulk is oriented with the Miller index h, k, l direction of the desired surface. Another approach to reduce the number of layers required to calculate accurate γsurface is to saturate the dangling bonds with pseudo-hydrogens containing fractional core and electronic charges,19 especially on surfaces that are not charge-neutral,20 as recently demonstrated on materials for PV applications.21

The surface energies of a number of Miller index planes can be combined to predict equilibrium particle morphologies and aspect ratios. This is achieved by implementing the Wulff construction, which minimises the surface energy for a given enclosed volume. Using the Wulff construction one can study the equilibrium shape of particles of arbitrary sizes, as demonstrated for TiO2 by Barnard and Curtiss,22 and analogous studies of several minerals23,24 and energy materials.17,25 Recently, Ong and co-workers, in collaboration with the Materials Project,26 categorised the equilibrium morphologies of more than 72 elements in the periodic table using DFT calculations. However, the Wulff construction does not account for contributions from lines (or edges) that join two distinct surfaces (i.e., line energies), and can predict particle morphologies that considerably deviate from shapes observed experimentally.22

The characterisation of surfaces and their surface energies is crucial to build heterogenous interfaces between materials, such as heterogenous semiconductor,27 metal//ceramic,28 and thin-film interfaces, all of which are relevant for energy applications. The energy penalty to separate two materials forming a heterogeneous solid//solid interface or a grain boundary in the same material is defined as:

$$E_{{\mathrm{interface}}} = \left[ {V_{{\mathrm{particle - A}}} + V_{{\mathrm{particle - B}}}} \right] \ast \varepsilon _{{\mathrm{strain}}} + + A_{{\mathrm{interface}}}\gamma _{{\mathrm{interface}}}.$$
(2)

Ainterface is the interface area, γinterface is the interface energy per unit area (see Eq. (3)), Vparticle is the volume of each particle, and εstrain is the volume-averaged strain energy. The interfacial contribution scales with the area shared at the interface Ainterface, whereas the strain contribution scales with the volume of each particle. As a result, the interfacial energy will dominate at small particle sizes, while the strain energy contribution will prevail at large particle sizes.

γinterface (in J m−2) describes the energetics of the interaction between two materials A and B to form an A//B interface.

$$\gamma _{{\mathrm{interface}}} = \mathop {{{\mathrm{lim}}}}\limits_{N \to \infty } \frac{{G_{{\mathrm{interface}}}(N) - N_AG_{{\mathrm{bulk}},{\mathrm{A}}} - N_BG_{{\mathrm{bulk,B}}}}}{{2A_{{\mathrm{interface}}}}}.$$
(3)

Ginterface, Gbulk,A and Gbulk,B are the free energies (approximated as the total energies in DFT calculations) of the interface system, bulk A and bulk B, respectively. Ni is the number of formula units of i in the interface system. Similar to surface energies (see Eq. (1)), γinterface is converged within a finite number of formula units in typical calculations, with the main constraint arising from the computational cost of first-principles calculations. An off-stoichiometric interface can be mathematically treated similar to an off-stoichiometric surface by accounting for chemical potentials of added (removed) species (Eq. (1)).

In contrast to γsurface, which is always positive, γinterface can assume positive or negative values, signifying that a A//B heterogeneous interface may be thermodynamically more stable than the respective bulks in isolation. Particularly, γinterface is defined by the chemical bonding between the materials forming the interface, thus controlling the adhesion properties.

Another physical quantity of interest is the adsorption energy of a chemical species (ΔGAds) onto a solid surface, which is given in Eq. (4).

$${\mathrm{\Delta }}G_{{\mathrm{Ads}}} = G_{{\mathrm{Surf + Ads}}} - G_{{\mathrm{Surf}}} - n\mu _{{\mathrm{Ads}}}.$$
(4)

where GSurf+Ads and GSurf are the free energies of the species formed by the solid surface and the adsorbing species in the solid state, and the surface in the solid state, respectively. μAds is the chemical potential of the n-species adsorbed.

Computing meaningful adsorption energies requires an adequate choice of the thermodynamic reference for μAds. The value of μAds varies significantly depending on the environment, such as liquid or gas. While μAds(v, 0K) in vapour can be easily approximated by the DFT total energy (or computed enthalpy),3 determining an accurate reference for liquid species is not trivial. One empirical strategy to reference the adsorbate to the liquid state is to subtract the experimental heat of vaporisation ΔHvap, when available, to the gas reference as in Eq. (5).29,30

$$\mu _{{\mathrm{Ads}}}(l,T,P) = \mu _{{\mathrm{Ads}}}(v,{\mathrm{0K}}) - {\mathrm{\Delta }}H_{{\mathrm{vap}}}(0{\mathrm{K}}) - TS(l).$$
(5)

where S(l) is the tabulated experimental value of entropy of the adsorbate in the liquid state. A similarly effective empirical correction has been implemented to define μAds(l, T, P) from a computed reference.31,32 The model uses the experimental fugacity f at the saturated gas pressure of the liquid adsorbate at a given T. Note that the liquid and vapour phases of the adsorbate are in equilibrium at the saturated vapour pressure, i.e., μl = μv. Hence, μAds(l, T, P) can be defined as:

$$\begin{array}{*{20}{l}} {\mu _{{\mathrm{Ads}}}(l,T,P)} \hfill & = \hfill & {\mu _{{\mathrm{Ads}}}(v,{\mathrm{0K}}) + {\mathrm{\Delta }}H(v)(0 \to T) + } \hfill \\ {} \hfill & {} \hfill & { - TS(v,T) + {\mathrm{RT}}\,{\mathrm{ln}}f(T).} \hfill \end{array}$$
(6)

Notably, μAds(l, T, P) in Eq. (5) is identical to Eq. (6) given that μ is a thermodynamic state function. Hence, both the aforementioned approaches use the same experimental thermodynamic data and should yield identical results. However, the common limitation of using either Eqs. (5) or (6) is the need of a priori experimental data (i.e., ΔHvap(T) or f(T)), which may not be available for all substances of interest. A brief survey of the NIST tables suggests that experimental ΔHvap is more readily available than f(T), across different chemistries.

For liquid phase adsorbate reactions involving a proton transfer, especially relevant in catalytic processes, Nørskov and co-workers31,32 formulated a thermodynamic scheme called the “computational hydrogen electrode” (CHE) to obtain the liquid phase reference from a computed gas-phase of H2(g). The CHE is based on the reversible hydrogen electrode (RHE) reaction:

$${\mathrm{H}}^ + + {\mathrm{e}}^ - \leftrightarrow 0.5{\mathrm{H}}_2({\mathrm{g}})$$
(7)

whose potential is set to 0 V. Following Eq. (7), the chemical potential of the pair, μ(H+) + μ(e), is equal to \(0.5\mu _{{\mathrm{H}}_2({\mathrm{g}})}\). An estimate of \(\mu _{{\mathrm{H}}_2({\mathrm{g}})}({\mathrm{0K}})\) is derived from DFT calculations in combination with appropriate enthalpic and entropic corrections (obtained either from advanced theoretical calculations or experiments). Subsequently, μ(H+) + μ(e) can be shifted to the desired potential by applying the Nernst’s equation ΔG = −zFV (see later Eq. (10)), where V is the applied bias. A pH correction to Eq. (7) is not required since CHE is set at 0 V at all pHs. The CHE method has been successfully applied to a number of studies in electrocatalysis.31,32 With suitable modifications, the CHE method can be extended to other solvated ion-electron pairs (see discussion in the section ‘Challenges and opportunities').

Electron energies

In this section, we consider the important electron energy levels used to predict and design materials for energy applications, as depicted in the right panel of Fig. 1 for a typical semiconductor.33 These physical quantities define which way charges will flow, how much energy will be lost by carriers (electrons and holes) and what electrical bias is required for operation across materials interfaces. Pedagogical discussions of these quantities are addressed by Peressi et al.34 and works from Kahn and Bredas.35

Vacuum level (Evac): The potential experienced by an electron ‘a few nanometres’ outside the material surface, and should be differentiated from the vacuum level at infinite distance from a materials surface, Evac(∞). Valence band maximum (VBM): The value of the highest occupied eigenvalue, when referenced against Evac gives the ionisation potential (IP). Conduction band minimum (CBM): The value of the lowest unoccupied eigenvalue, when referenced against Evac gives the electron affinity (EA). Band offset: The difference in VBM or CBM in two materials forming a junction. Dipole (Ds): The difference in potential far from the interface, either between two materials or material and vacuum; illustrated in Fig. 1. In electrochemistry, this is known as the Galvani potential. Fermi level (EF): the electrochemical potential of electrons (set by the Fermi Energy) in the system; EF depends on carrier concentrations. Workfunction (Φ), EF referenced to Evac.

The lineup of bands in materials and (un)occupied electron energy levels in liquids and gases is one of enormous importance in semiconductor heterojunction devices.36,37,38,39 For gas molecules the calculation of the electron energies is relatively straightforward. Calculating the electronic structure of a molecule in gas phase provides an inherent reference potential and the electron eigenvalues are on an absolute scale. In solids the situation is more complex, as the absolute value of the potential has no meaning in a periodic lattice.15 Thus, it is necessary to introduce a reference potential.40 There have been several attempts to introduce reference potentials, such as model solid theory and the hydrogen defect level.41,42

Anderson’s rule43 states that when two semiconductors are brought into contact, their vacuum levels must coincide at the interface. By linking this to the definitions of IP and EA, Anderson’s rule provides an intuitive route to estimating the offsets of the valence and conduction bands of materials. Since IP = Evac − VBM, we can write the valence band offset (VBO) between the two materials as VBO = IPA − IPB. Similar arguments give the conduction band offset, CBO = EAA − EAB. Estimating band offsets via Anderson’s rule is similar to that of the Mott and Schottky procedure to determine energy level alignments at metal-semiconductor interfaces, where the offset is ΔE = EA − Φ, with Φ being the workfunction of the metal. When a junction is formed between two semiconductors the Fermi levels align through charge transfer, this results in a bending of bands and vacuum levels in the region of the interface, the space-charge region.44 The bending of bands and vacuum levels means that while the vacuum level aligns at the interface, they can be quite different on different sides of the junction far from the interface, this difference is referred to built-in potential. Often there are additional electronic states at a surface or interface, filling or emptying of these states can result in microscopic dipoles, that cause deviations from Anderson’s rule predicted offsets.44 The deviation of an interface band offset from the Anderson’s rule predicted offset has been characterised in terms of a density of induced states.45 Similar deviations can be induced via a mismatch in lattice parameters of materials (see Eq. (2)) that form an interface, typically quantified by a “deformation potential.”41,46

The theoretical band alignment approach most closely related to experimental measurements is the slab model approach,47,48 where the electrostatic potential in the vacuum (Evac, see Fig. 2) provides the reference. Thus, the IP can be calculated as

$${\mathrm{IP}} = E_{{\mathrm{vac}}} - {\mathrm{VBM}}_{\mathrm{s}}.$$
(8)

where VBMs is the VBM of the slab. Generally, slab models introduce surface states from dangling bonds, which means that the IP from a single slab calculation does not represent the IP to remove an electron from the bulk of the material. However, dangling bonds are an artefact of slab models that only consider ideal surface cuts. Typically, surfaces undergo reconstructions to avoid surface electronic states (there are exceptions),20 and satisfy the “electron counting rules”49 as a result. Nevertheless, the effect of dangling bonds can be removed by using the VBM of the bulk, VBMb. The IP is then

$${\mathrm{IP}} = D_s - {\mathrm{VBM}}_{\mathrm{b}}.$$
(9)

where Ds is the surface dipole of the slab calculation (see Fig. 2). To obtain the potential in the bulk we need to account for the effect of the oscillating potential—this is achieved by calculating the macroscopic average of the potential, as seen in Fig. 2. The macroscopic average is commonly defined as the running average over the repeating lattice unit normal to the surface.34

Fig. 2
figure 2

Important electron potentials for calculating band offsets. The planar and macroscopic averages of the Hartree potential derived from a DFT slab calculation are plotted in blue and orange. The surface dipole (Ds) and vacuum potential (Evac) are indicated

The branch point energy EBP is another approach used to align bands across bulk semiconductors. EBP is defined as a level at which defect states in the band gap change from donor-like to acceptor-like states and is said to be universal across semicondutors.45,50 The most popular EBP procedure is to calculate the band gap centre (BGC),51 which is the average of the VBM and CBM across the Brillouin zone. The primary advantage of calculating the BGC is its insensitivity to the choice of exchange-correlation functional in the DFT framework,51 which allows the use of computationally inexpensive calculations. However, the BGC approach is sensitive to the choice of computational parameters.52,53

The natural band alignment approach starts by first calculating explicitly the heterojuctions of a series of semiconductors and the band offsets are inherently obtained. By assuming transitivity of offsets, e.g., Δ(A/C) = Δ(A/B) + Δ(B/C), one derives offsets for a wide range of materials.54,55 The strength of this approach is that the methodology is closely linked to spectroscopic approaches for band alignment, and thus provides a strong link to experimental studies. A drawback of the natural band alignment is the assumption of transitivity, which is only reasonable between materials of similar structure and chemistry. Natural band alignment has found success within families of materials, e.g., III–V, II–IV, and chalcopyrite semiconductors, and recently, hybrid halide perovskites.56

The presence of impurities at a surface can lead to a change in the electron energies of materials. Addition of thin layers at the surface, either self-assembled monolayers of molecules or films of hetero-structure crystals have been shown to be a promising route to tuning IPs and EAs for particular applications57,58 (see section ‘Interfaces in photovoltaics’).

In the case of a semiconductor-liquid (or a semiconductor-metal) interface, the band alignment at the surface differs significantly from the bulk, due to band bending (depicted in the PC cell in Fig. 2). This has been attributed to the difference in Ds at the semiconductor-liquid interface compared to semiconductor–vacuum and liquid–vacuum interfaces.59 In the case of photocatalysts that aim at splitting water, band bending can enable (or disable) the catalytic activity of a semiconductor60 (see section ‘Interfaces in photocatalysts’).

Energy applications

Interfaces in batteries

This section describes the properties regulating interfaces (and interphases) in battery devices that can be accessed via computational materials science methods.

In rechargeable batteries, such as Li-ion cells, the active species (Li) is reversibly intercalated between the low potential anode—the negative electrode that undergoes oxidation—and the high potential cathode—the positive electrode that undergoes reduction, across a liquid/solid/polymer electrolyte medium during successive discharging and charging cycles. Thus, the potential difference, also known as the open circuit voltage (VOC), between the cathode and the anode drives the battery discharge, and is defined by Nernst’s equation (below).

$$V_{{\mathrm{OC}}} = - \frac{{{\mathrm{\Delta }}G}}{{xz{\mathrm{F}}}} = - \frac{{\mu _{{\mathrm{Li}}}^{{\mathrm{cathode}}}(x) - \mu _{{\mathrm{Li}}}^{{\mathrm{anode}}}}}{{xz{\mathrm{F}}}}$$
(10)

where ΔG is the free energy gain driving the battery discharge of x active species (Li) atoms. μLi is the chemical potential of Li (or the equivalent active species), F is the Faraday constant, and z is the charge of the active specie (=1 for Li). Equation (10) illustrates that the VOC is directly proportional to the chemical potential of the active species. Thus, electrodes can be ranked according to their potentials relative to the standard reduction potential (SRP) of the active species (e.g., Li/Li+ ~ −3.04 V vs. SHE), and a combination of a high potential cathode with a low potential anode can store significantly more energy than the vice-versa.61 Often, in batteries’ research the terms “potential” and “voltage” are used interchangeably (e.g., a high-voltage (potential) cathode assembled with a low voltage (potential) anode), since the electrode potentials are always referenced to the SRP (or another suitable electrode) for all practical measurements and the numbers reported in studies are, actually, voltages.

First-principles calculations, especially DFT, utilising the concepts exposed in Eqs. (10) and (1) have been applied systematically to calculate the bulk intercalation VOC of electrode combinations,25,30,61 the particle morphologies of electrode materials,25 and solid electrolytes,17 and understand the electrode–electrolyte interactions to form interphases.62,63 Specifically in electrolytes, the maximum oxidising (against the cathode) and the minimum reducing (against the anode) potentials define the “voltage window” of the electrolyte, as shown in Fig. 3.

Fig. 3
figure 3

Comparison of battery (left), PV (centre) and PC (right) devices. Left shows a schematic of a battery, the open circuit voltage (VOC) across the cell (Vcat → Van) is depicted by a solid black line, this determines the chemical potential of the moving species (e.g., μLi), materials interfaces are formed at the electrode//electrolyte interfaces. The anodic (against the cathode) and cathodic (against the anode) potentials of the electrolyte VAc and VCv are also shown. Centre shows a typical pn solar cell, the solid black lines show the valence and conduction band edges, which track the voltage across the cell, offsets between valence and conduction bands (ΔEv and ΔEc) and Schottky barrier (Φ) are shown. Right indicates a PC cell, with the valence and conduction band edges displaying the bending behaviour at the bulk semiconductor//electrolyte interface. The conduction (valence) band edge at the interface has to be above (below) the H2 (O2) evolution potentials for spontaneous water-splitting to occur upon photoexcitation of electrons and holes

Electrolytes are inherently “stable” against battery electrodes if their voltage windows span across the potentials of the anode and the cathode. As an example, if the anode in a Li-ion battery is assumed to be Li-metal and the cathode LiCoO2, the average voltage of operation is ~3.2 V vs. Li+/Li, which is equivalent to a ΔμLi of ~3.2 eV or ~308.8 kJ mol−1. Hence, an electrolyte will be stable against both Li-metal and LiCoO2 if its voltage window spans at least 0–3.2 V. However, an energy scale of ~3.2 eV is normally sufficient to drive chemical reactions and most electrolytes, as a result, are unstable against both electrodes. Additionally, there is an abrupt change of the chemical potential at the interface between electrolyte and electrodes (Fig. 3) that needs to be accommodated by the electrolyte. Therefore, most electrolytes display a thermodynamic driving force towards decomposition, and the resulting new phases form the solid-electrolyte-interphase (SEI) along the electrode//electrolyte interface in batteries.

If the thermodynamic decomposition of the electrolyte is not kinetically limited, three situations are likely to occur:64,65 (i) Formation of new phases blocking electrons and transparent to the migration of the active ion (e.g., Li+), which is ideal as the new interphase “protects” the electrolyte from further decomposition. (ii) A mixed-conduction interphase, where new phases can transport both electrons and the active ion are formed. Mixed-conducting interphases are undesired; they consume the active ions in the electrolyte, contributing to the phenomenon of cell self-discharge, and eventual cell failure. (iii) A non-conducting interface that blocks both electrons and ions (e.g., MgO formation in Mg batteries66), which “passivates” the electrode surface causing electrochemical reactions in the battery to stop.

The magnitude of the electrolyte stability windows can be used to inform whether an electrolyte will decompose at the potential of the electrodes. To a first approximation, decomposition requires that the highest occupied molecular orbital (HOMO) of the liquid electrolyte is lower than the CBM of the cathode, and/or the lowest unoccupied molecular orbital (LUMO) of the electrolyte is higher than the VBM (or Fermi energy for a metal) of the anode.

A method has been proposed to chart the stability windows of a number of liquid electrolytes via Eq. (11),67 where VCv, VAv are the cathodic and the anodic potentials, and VCv − VAv is the electrolyte stability window (ΔV).

$${\mathrm{\Delta }}V \equiv V_{{\mathrm{Cv}}} - V_{{\mathrm{Av}}} = - \frac{{\varepsilon _{{\mathrm{LUMO}}}}}{{e^ - }} - \frac{{\varepsilon _{{\mathrm{HOMO}}}}}{{e^ - }}.$$
(11)

The authors suggest a computationally viable but accurate approximation of both the HOMO and LUMO levels (ε) via the Δ-SCF method,15 which derives ε as the difference between the computed total energy of an ion and the same ion with an electron added (LUMO) or removed (HOMO). The authors also considered the effect of a solvent in theoretically determining VAv, and VCv, and found the use of explicit solvent molecules to significantly improve the agreement with experiments.67 Furthermore, the vacuum electron potential (Fig. 2) was used to rescale VCv and VAv to a common reference, e.g., Li-metal.

First-principles calculations of electrode//liquid electrolytes

The decomposition of liquid electrolytes can result in the formation of several intermediate products, requiring extensive computations or experiments to characterise accurately. For example, the exact composition of the SEI that forms due to the decomposition of LiPF6 in carbonates at the reducing potentials of the Li-graphite anode in Li-ion batteries has been debated for the past 20 years.68 Recently, theoretical advancements have been made in the description of the SEI69 via extensive ab initio molecular dynamics (AIMD) simulations of liquid carbonates in the presence of a graphite surface. For example, Leung et al.69 identified the preliminary stages of electrolyte decomposition, providing a useful fingerprint to interpret spectroscopy experiments. Suo et al.70 studied the SEI combining a number of experimental techniques with ab initio calculations of the solid//liquid interface to estimate the electronic tunnelling barriers and the thickness of the SEI. Additionally, a combination of static adsorption calculations (see Eq. (4)) and AIMD simulations of the liquid electrolyte//solid anode interface was recently used to study the initial stages of Mg deposition/stripping in Mg batteries.29

The primary challenge in describing solid//liquid interfaces in batteries is the accurate treatment of an electrified surface by introducing an appropriate potential at the electrode, which mimics the formation of a polarized double layer (discussed later in the context of polarisation).71,72 Also, ab initio methods can only access the initial stages of electrolyte decomposition in Li-ion batteries since they cannot capture the length and time scales relevant to the mesoscopic nature of the SEI. Nevertheless, the ongoing development of reliable many-body potentials derived from first-principles data, especially via machine-learning strategies,73,74 may mitigate this constraint soon.

Mapping electrolyte decomposition of solid electrolytes

A solid electrolyte (SE), which exhibits significant ionic conductivity similar to liquid electrolytes, can also form decomposition products upon contact with the electrodes. Notably, Zhu et al.64 and Richards et al.65 developed a methodology based on first-principles DFT calculations to evaluate the grand-potential (Φ[c, μLi], see Eq. (12)) and chart the stability windows of SEs for Li-batteries, while Aykol et al. developed a similar approach to identify new coating materials to protect high-voltage cathodes.75 Specifically, Φ can be written as:

$${\mathrm{\Phi }}[c,\mu _{{\mathrm{Li}}}] = G[c] - n_{{\mathrm{Li}}}[c]\mu _{{\mathrm{Li}}}.$$
(12)

where G[c] is the Gibbs energy for a given Li-containing phase, c. nLi and μLi are the number of Li atoms in c and the Li chemical potential, respectively. A phase c can be thermodynamically stable across a range of Φ, i.e., across a range of μLi, which in turn is equivalent to a range of potentials via Eq. (10), or the stability window of the SE (ΔV, see Eq. (11)). Analogous to liquid electrolytes, the SE will form other thermodynamically stable phases beyond its stability window, which can be compared directly with (or provide guidance to) experimental observations.76,77 For example, Fig. 4 shows the computed stability window (~1.71–2.14 V) of Li10GeP2S12 (LGPS), as well as other secondary phases (e.g., Li3P, Li2S and Li-Ge inter-metallics) that can form as a result of LGPS decomposition.77 The primary drawback in the thermodynamic approach of Eq. (12) is that it neglects interphases that can be stabilised by kinetic effects.78

Fig. 4
figure 4

Electrolyte stability window, voltage profile, and phase equilibria of Li10GeP2S12 upon Li-insertion and extraction. The new thermodynamically stable phases formed outside the stability windows of Li10GeP2S12 are shown. Reprinted with permission from ref. 77. Copyright 2016 Wiley-VCH

Recently, Tang et al.79 employed AIMD simulations to assess the existence of kinetically stabilised interphases by creating an explicit heterogenous SE//electrode interface, following Eq. (3), and minimising the lattice mismatch (see the section ‘Electron energies') of the SE and the electrode surfaces upon their alignment. Subsequently, the authors performed short AIMD simulations and matched the radial distribution functions (RDF) of various atomic pairs against the RDFs of decomposition compounds predicted to form via the formalism of Eq. (12). Through a process of matching and elimination the authors were able to identify specific interphases that may be more kinetically stabilised than others. Although promising, a major challenge with AIMD simulations of solid//solid interfaces is the limited time propagation, analogous to solid//liquid interfaces.

Polarisation of interfaces

In circumstances where interphases (if any) between electrode materials and the electrolyte are stable, polarisation effects at the interface should be taken into account.63,80 Polarisation typically causes a sharp potential variation at the interface, as shown in the left panel of Fig. 3. For liquid electrolytes, the potential drop and the formation of double layer at each electrode is driven by the simultaneous alignment of the dipoles of solvent molecules (e.g., carbonates) at the electrode surfaces and the concentration of the ionic species in the electrolyte. In the case of SEs, the double layer is formed by a concentration gradient of the mobile active species (e.g., Li+ and Li-vacancies) since the anion framework (and other cations, if any) remains “frozen”. The width of the double layer, in both liquid and solid electrolytes, is determined by the corresponding dielectric constants.

Using DFT, Haruyama et al. were the first to attempt the description of the space-charge layer of various oxide cathode//sulphide SE interfaces in Li-systems.63 However, the authors entirely neglected the occurrence of electrochemical reactions driving the formation of new phases at the cathode//SE interface,64,65 which clearly dominate interfacial reactions over polarisation. Stegmaier et al.81 adapted some of the concepts used to describe semiconductor interfaces to address the polarisation of a SE//cathode interface. The authors reported that significant energies (~1 eV) are required to stabilise the segregation of Li-vacancies in the interfacial region. In addition, the authors estimated a significant potential drop (~3 V) across the interface (see also the following section), which needs further verification.

Interfaces in photovoltaics

This section describes the properties regulating interfaces in different PV cells that can be accessed with first-principles methods. First, we consider the general operating principles of a solar cell and how interfaces affect efficiency, then look at specific types of PV technology and the important interfaces therein. Also note that a number of factors determine the optimal alignment of energy levels across the interfaces in a functioning solar cells, which depend on the specific device architecture.

Interface effects on photovoltaic efficiency

The efficiency of an operational solar cell can be described by an equivalent electrical circuit, where resistors in series represent the loss of efficiency. The sources of series resistance are the electrical resistances of the absorber, contact layers, and the interfaces between those layers. The resistance from an interface is labelled the contact resistance (Rc).

Rc depends on a number of properties, including carrier effective masses and space-charge-region (center panel of Fig. 3) widths. However, the most important factor, which is also the most readily tractable from first-principles calculations is the energy barrier to charge flow across the interface.

pn and pin solar cells

In pn junction solar cells, the photo-generated charge is separated at an interface between p- and n-type semiconductors. In this case we require a band offset at the interface between the two semiconductors, as depicted in Fig. 3. Examples of pn junction solar cells are crystalline Si cells and thin-film devices, such as CdTe, Cu(In, Ge)Se2 (CIGS) and Cu2ZnSnS4 (CZTS). In pin solar cells, such as those based on amorphous Si and hybrid perovskites, the absorber layer is not heavily doped (i.e., intrinsically doped, i). In this architecture a field persists across the absorber layer and carriers are separated within the absorber.

The offset between the absorber and partner layer is critical for the operation of p (−i)−n type devices. The offset should be such that electrons have a driving force to move in one direction, while holes move in the opposite direction. One follows the principle that electrons will move towards lower potentials (sink) and holes towards higher potentials (float). In the case of contact layers in a pn PV junction, the band alignment should be such that the extraction of the majority carriers from each layer (i.e., e from n and h+ from p) is optimised.

If the absorber layer is p-type, then electrons (the minority carriers) need to be extracted to the heterojunction partner (and subsequently to the contact layer) and the CBO is the important property (see Fig. 3). It has been found that for p-type absorbers with a small barrier (<0.2 eV) to electron transport into the partner layer results in the most efficient devices.82 The reason that this barrier is favourable relates to the peculiar alignment of Fermi levels to enhance a depletion layer in the absorber near the interface, which results in a space-charge region that repels holes limiting the recombination of electrons and holes at the interface. Although there is a small barrier to transport at the interface the electrons can still tunnel through.

If the absorber is n-type, holes need to be extracted to the heterojunction partner. In general, holes have higher effective masses than electrons and therefore, do not tunnel as readily. As a result, offsets with no barrier to hole transport are generally required. Recently, multi-scale approaches have been developed to model the space-charge region in PV systems,83 with a similar model adopted in battery research.81

In thin-film solar cells, the combination of X-ray photoemission spectroscopy (XPS) and electronic structure calculations has proved to be highly successful in identifying efficiency-limiting factors.84,85 For example, XPS and theory results have shown that the architecture of CdTe solar cells with CdS buffer layers is not appropriate for CuSbS2 or SnS absorber layers, due to poor band alignment.84,85 Similarly, the important effects of interface stoichiometry have been demonstrated.86 However, the existence of different reference values and limited instrument resolution has often hampered direct comparisons between XPS and DFT data, with such comparisons requiring a series of tedious and ad hoc post-processing steps. The free software package, galore,87 addresses this issue by applying systematic corrections to calculated data for direct comparison with experiments. Another important consideration is the orientation of interfaces/surfaces used, see Fig. 5.88,89

Fig. 5
figure 5

The interplay of surface energy (top panel) and electron energy (bottom panel). The equilibrium morphology and band alignments of SnS for various surface termination of SnS have been calculated88 and highlight the importance of termination for measured IPs/EAs. Reproduced with permission from Appl. Phys. Lett. 104, 211603 (2014). Copyright 2014 American Institute of Physics

In halide perovskite solar cells there has been extensive work on modelling band alignments which has shown how varying the halide anion can affect the band matching.56 There have been many studies of explicit perovskite//oxide interfaces, where properties, such as interface adhesion (see Eq. (3)), band alignment, charge transfer and defect states are calculated directly.90 Size effects can also play a crucial role, often very thin layers of interface are used for models, which means that true bulk alignments are not obtained. The presence of dipoles across the system (a 2D slab in periodic DFT) can create artefacts—systems should be carefully checked to ensure that no electric dipoles are introduced by the asymmetry of the slab.

Bulk heterojunction and dye-sensitised solar cells

Bulk heterojunctions are typically organic solar cells and have a fundamentally different architecture to cells with inorganic absorber layers. In BHJ cells the absorber layer is a blend of either small molecules or polymers of p- and n-type. This blend results in a morphology where there are interfaces between the p- and n-layers extended throughout the absorber layer.

In BHJ cells the absorption of light form excitons (bound electron–hole pair) in the absorber molecule. Excitons must be split apart to extract current. The field across the cell is insufficient for splitting an exciton, therefore a second molecule with a different IP and EA is required.

In a BHJ cell there is typically a spacer layer between the BHJ blend and the electrodes. These spacer layers and interfaces perform a number of critical functions including, modulating the offset between absorber and electrode, selectively blocking minority carriers at the respective electrodes, and stabilising the interface of absorber//electrode.91

How the bands align in BHJ solar cells between molecules is generally relatively inexpensive to compute, as molecules are treated rigorously without periodic boundary conditions and there is an absolute reference potential level available. Therefore, offsets for molecule types can be calculated in a high-throughput fashion and optimal donor–acceptor pairs chosen.92 While these approaches are useful for screening, accurate understanding of interfaces in BHJs requires consideration of excited states, using, for example, costly time-dependent DFT (TD-DFT) or many-body perturbation methods. Effects introduced by structural disorder are also crucial.93

The degree to which alignment at molecule//solid interfaces is controlled by the levels of the constituent materials depends on the strength of chemical bonding, which can be assessed via Eq. (5). With weak bonding of physisorbed interfaces the alignment can be predicted by vacuum level alignment.94 Strong bonding (chemisorption) promotes charge transfers at the interface, resulting in an additional dipole layer that shifts the vacuum levels and affects the band alignment.95,96 This can further depend on molecular orientation at the interface and the presence of water layers.96 It was recently shown that there is a surprisingly homogeneous band alignment between organic molecules and metal oxides (often used as interlayers).97 This observation has been explained and critically assessed by numerical models of the organic//metal oxide interface.98

Computational methods have also been extensively applied to understand the importance of surface modifications to band alignments in BHJ solar cells. Polymers at a metal surface can be used to reduce contact workfunctions and improve device performance.99 A combination of spectroscopy and quantum chemistry calculations was recently used to elucidate a set of design rules for high VOC organic solar cells.99

In a dye-sensitised solar cell (DSSC), the excitons are formed in organic light-absorbing molecules, adsorbed to a porous metal oxide (typically TiO2). The carriers are generated by the dye then immediately transferred to the oxide, the original charge state of the dye is restored by charge transfer from an electrolyte solution in contact with the dye. In DSSCs the interfaces between the dye and the oxide and the dye and the electrolyte are of the utmost importance and interfacial adsorption regulates the kinetics of charge transfer. Thus, describing adsorption energies of the liquid dye (or redox mediators) on substrates is important.

Many of the modelling methods and concerns for DSSC are related to those for BHJ solar cells. Notably, the interface between dyes and oxides are very similar to an organic//interlayer interface in a BHJ solar cell. Excitonic effects of electron–hole coupling can be particularly important in DSSC interfaces. Notably, there has been a great effort to build more efficient many-body perturbation methods to accurately treat these systems.100

Interfaces in photocatalysts

Similar to BHJ and DSSC solar cells, photocatalysis relies on semiconductors that generate electron–hole pairs due to photoexcitation, which is then subsequently used to drive reduction (using generated electrons) and oxidation (holes) reactions. A common example of photocatalysis is water-splitting,101 where the hydrogen evolution reaction (HER) and the oxygen evolution reaction (OER) form the reduction and oxidation components. Importantly, the CBM (VBM) of a candidate photocatalyst, especially at the interface with the aqueous solution (or electrolyte), has to be positioned above (below) the potential of the HER (OER) reaction. Thus, the band gap of a PC semiconductor must be at least 1.23 eV and must span the potentials for HER and OER reactions. Analogous constraints can be extended to PC for CO2 reduction that eventually yields to fuels, such as CH3OH,101 and interfacial reactions in a fuel cell electrode with recombination of H2 and O2 to yield H2O.31

As outlined in the section ‘Electron energies', several approaches exist to calculate the band-edge positions within the bulk semiconductor and liquid, the alignment of which is important for a material to be able to split water.102 However, evaluating the band-edge changes at the interface is more challenging especially across a semiconductor-liquid due to band bending (see the section ‘Electron energies'). Theoretically, the extent of band bending can be directly evaluated via AIMD simulations of a large supercell that contains both the semiconductor surface and the bulk liquid in contact with each other.103 Nevertheless, such calculations are computationally demanding.

Wu et al.104 devised a three-step procedure to capture interfacial band bending that involved calculating bulk properties: (i) the CBM \(\left( {E_{{\mathrm{bulk}}}^c} \right)\) referenced to the bulk Hartree potential (Hbulk) of the semiconductor, (ii) the acceptor level of the liquid (Abulk) with respect to its bulk Hartree potential (Hsol–bulk), and (iii) the difference in the Hartree potentials at the semiconductor-liquid interface (Hedge − Hsol–edge). While step (iii) can be completed with AIMD simulations, steps (i) and (ii) are performed with suitable DFT approximations104 or GW approaches105—GW describes the band edges more accurately. Subsequently, the bending of the CBM at the semiconductor surface \(\left( {E_{{\mathrm{edge}}}^c} \right)\) is related to \(\left( {E_{{\mathrm{bulk}}}^c} \right)\) via the difference in the Hartree potentials of the interface and the bulk, as Hbulk − Hedge. By performing an analogous estimate of band bending (i.e., the acceptor level) within the liquid, the authors104 were able to accurately capture the band alignment of several semiconductor–water interfaces. Although this approach was further extended to a high-throughput search for photocatalysts, it still suffers from poor predictions of band-edge positions.

Notably, the kinetic barriers to separate the generated electron–hole pairs (reflected as overpotentials in experiments) and thermodynamic losses (e.g., losses incurred in the external circuit or other supporting components) typically lead to net losses of ~0.5 V.52 To account for both band bending and kinetic barriers, several computational searches of feasible photocatalysts have involved screening of semiconductors with a band gap of ~2 eV, significantly higher than the thermodynamic 1.23 eV value, which covers the OER and HER voltages.52,101,106,107 As an example, Fig. 6 illustrates the results from a computational screening of 3d-metal oxides as PC candidates, which identified MnO and NiO as promising PCs.101

Fig. 6
figure 6

Evaluation of the valence and conduction band-edge positions (solid black lines) of various candidate 3d transition metal oxides for PC applications, using a combination of DFT and G0W0 calculations. The dashed black lines for each oxide signifies its band gap centre position, while the coloured dotted lines correspond to various reduction/oxidation reactions, as indicated in the text legend. Two sets of black lines for Fe2O3 correspond to two different surface facets. Reproduced from ref. 101 with permission from The Royal Society of Chemistry

Experimental and computational studies, have strived to improve PC (and in general, catalytic) processes by identifying, and hopefully reducing, kinetic barriers (i.e., overpotentials) to a given reaction. Often, such barriers exist because of (un)stable adsorbed species on the catalyst surface and/or the poor kinetics of charge transfer across the catalyst//adsorbate interface. For example, nanocomposites of rutile- and anatase-TiO2 have shown enhanced photo-activity than the corresponding bulk forms due to better charge transfer.108 Thus, surface science, specifically the process of adsorption/desorption of catalyst//electrolyte interface plays a significant role in determining the overall rate and the efficiency of a given reaction. In this context, several studies have identified key aspects that influence kinetic barriers in general catalytic reactions, such as the surface Fermi level,109 d-band centre of metallic catalysts,4 lattice strain,110 hot electrons due to plasmonic excitations,111 electronic properties of the adsorbate,112 and other descriptors that correlate with reaction kinetics.113

Challenges and opportunities

The fields of PV, PC and battery research have developed many useful approaches to modelling interfaces that we have reviewed in this work. As suggested by the qualitative overview in Fig. 7, modelling approaches have been typically applied to a specific application domain. For example, studies of interfaces (and interphases) in the field of batteries and electrochemical systems (blue area in Fig. 7) have predominantly focussed on thermodynamic properties (e.g., electrochemical stability and surface energies), while PV and PC applications have employed models to accurately estimate electronic properties (e.g., band gaps and band alignment, red area in Fig. 7). In the current scenario, as highlighted by the dashed line in Fig. 7, the thermodynamic and electronic regimes are separated by a fictitious “barrier” that reduces the exchange of ideas between scientific communities. Hence, we hope that our work, in conjunction with scientific studies that follow, will lead to a universal theoretical understanding and design of interfaces, across all energy applications (white area above the dashed line in Fig. 7).

Fig. 7
figure 7

Qualitative overview of the interfacial properties either computed or measured by scientists working in materials design for battery, photovoltaic and photocatalytic applications. The dashed line is an indicator of the current scenario where thermodynamic properties (e.g., surface energies, chemical potentials/voltages and electrochemical stabilities) are engineered for battery materials (blue area), while electronic properties (e.g., band gaps, band alignments and Fermi level) are the primary focus in the development of photovoltaic and photocatalytic devices (red area). These two seemingly distinct realms are separated by a “barrier” that mitigates the exchange of ideas between scientific communities. The white area above the dashed line identifies the ideal situation where a universal theoretical understanding of interfaces is identified

Several examples exist where an exchange of concepts and methodologies can be useful in improving the design of interfaces. For example, the thermodynamic stability (see the section ‘Mapping electrolyte decomposition of solid electrolytes') of a solid//solid interface (or interphase) is explicitly examined while designing an electrode//SE combination in batteries. Nevertheless, theoretical (and experimental) studies of PV (e.g., CdTe/CdS-based devices114) and PC devices largely ignore such thermodynamic contributions. Note that the precise nature of the semiconductor interface that forms the space-charge region in a PV (and PC) device can influence its overall efficiency. Similarly, polarisation of interfaces (right panel of Fig. 2) is a recurring theme in studies on PV and PC interfaces, utilising decades of knowledge built upon engineering diodes and transistors, while battery scientists have only recently begun explicitly including the impact of polarisation on the overall electrochemical potential of electrode//SE interfaces.115

Notably, the formalism of the computational hydrogen electrode (see Eq. (7)) in photo/electrocatalytic interfaces31,32 can be extended (i) to include solvent effects that are typically ignored116 at the light absorber//solid substrate interface in a DSSC, (ii) to model the redox-mediator-interface in Li-air batteries, and (iii) to study the (potential) decomposition of liquid electrolytes against an electrode under an external potential in common battery systems. The computational hydrogen electrode formalism, specifically developed to calculate an accurate reference for μ(H+) + μ(e) in an aqueous solution, can be further generalised to calculate μ(Li+) + μ(e), μ(Na+) + μ(e), or μ(Mg2+) + 2μ(e), while accounting for variations in the solvent environment by either calculating or experimentally measuring the corresponding ion solvation energies.

On a positive note, there are studies where concepts and methods originally developed for PV (or PC) applications have been successfully exchanged to the battery community (and vice-versa). For example, the band alignment strategy (see the section ‘Electron energies') has been successfully applied to battery systems, including the determination of stability windows of liquid electrolytes,67 the study of electrolyte decomposition at the anode electrode,69,70 and the elucidation of the charge transport and tunnelling mechanisms in Li–O2 batteries.117 Such instances shed hope on eventually arriving at models that enable a universal understanding and design of interfaces, covering all types of energy applications (white region in Fig. 7).

Apart from the thermodynamic and ground-state electronic properties, kinetic, excited-state electronic, and mechanical properties of interfaces are also highly relevant in determining the overall efficiency of energy devices and need to be precisely determined. For example, ion transport across electrode//electrolyte interfaces (that determine the power density) in batteries,118 charge separation, transport, and excited-state dynamics across the space-charge region (carrier density and short circuit current) in PV,119,120,121 electron/hole transport at a photocatalyst//electrolyte interface (overpotential)101 are crucial kinetic and excited-state properties that need to be described accurately.

Similarly, the mechanical properties of interfaces, including epitaxial strain and elastic or plastic deformation in solid//solid interfaces and the extent of wetting in solid//liquid interfaces, can influence device efficiencies as well. Reviewing kinetic and mechanical properties that are typically computed using ab initio methods is beyond the scope of this work, and we refer the reader to relevant literature on these topics.110,122,123 Also, we note that the contribution of defects at an interface (e.g., space-charge region), requires more detailed investigations that add to a few recent studies.63,81,124

Nano-scale interface effects are increasingly becoming technologically important. In quantum dot solar cells (PV) and PC, nanoparticles are of particular interest. Nano-sized materials present a new set of modelling challenges, as quantum confinement can affect electronic and structural properties. Often it is necessary to simulate the entire nano-sized object for accurate results. Efficient DFT methods and multi-scale approaches show promise for accurate calculations of these important effects.125 Detailed review of this field is beyond our remit here, but there is extensive literature on this topic.126,127

In general, caution is warranted as the selection of computational parameters, in particular the treatment of the exchange and correlation functional in DFT, can provide significantly different results while simulating interfaces. The overall accuracy and predicting capabilities of GGA, GGA + U and hybrid functionals do benefit from significant error cancellation. Therefore, the development of strategies, in the spirit of the Bayesian error estimation functional by Wellendorff et al.,128 represents a pathway to improve the reliability of predicted results.

For modelling practical interfaces, there is the pressing need to develop the ability to simulate “large” scale models of interfaces that can also simulate “long” range behaviour with high accuracy. In this regard the development of efficient linear-scaling DFT and GPU accelerated approaches are going to be important for future studies.129,130 Another approach to tackling large systems is to use multi-scale models, where first-principles calculations are used to parameterise larger scale simulations (e.g., interatomic potentials, mesoscale, continuum modelling and phase-field).131 Although, effort is required in estimating the propagation of uncertainties across space-time scales in multi-scale models.

Materials informatics is quickly changing the face of modelling studies, signified by the growth in material-property databases132,133 and machine-learning driven studies.74,134 For bulk systems, reliable descriptors of crystal structures are now available to represent the structures and the corresponding calculated (or measured) properties in a suitable form for machine-learning models, which can subsequently be used to predict properties of novel materials. Particularly, descriptors that can accurately predict macroscopic interfacial properties, such as the extent of band bending at interfaces, will be extremely useful.135

However, such descriptors are lacking for interfacial systems, which is probably justified by the limited number (if any) of freely available datasets of interface structures and properties. Hence, it is desirable to create databases that can compare heterogenous interfaces of thousands of materials, where calculations are performed in a high-throughput64,75,132,133,136 fashion with a systematic calibration of the computational parameters. Thus, a community effort is required to develop resources to take advantage of the “fourth paradigm” of scientific research137 for modelling interfaces.

Conclusions

We have reviewed the state-of-the-art in materials modelling of interfaces in battery, PV and PC technologies. We started by providing an overview of the “pencil-and-paper” theory behind the properties that can be calculated from electronic structure simulations. Subsequently, we focussed on the applications of these concepts in calculating the material properties for practical energy applications. Over the course of reviewing the various theoretical approaches, we find that all the aforementioned energy fields have developed a range of unique techniques for calculating different properties, many of which can be extremely useful for allied energy applications—we have outlined a few of these prospects. Finally, we identify some common challenges in modelling interfaces across energy disciplines, where progress will be crucial for the continued success of computational materials science in understanding, designing and improving interfaces. It is our hope that this review will help to encourage the cross-fertilisation of ideas across a range of modelling communities to the general benefit of computational materials design.