Many-body Green's function perturbation theories, such as the GW and Bethe–Salpeter formalisms, are starting to be routinely applied to study charged and neutral electronic excitations in molecular organic systems relevant to applications in photovoltaics, photochemistry or biology. In parallel, density functional theory and its time-dependent extensions significantly progressed along the line of range-separated hybrid functionals within the generalized Kohn–Sham formalism designed to provide correct excitation energies. We give an overview and compare these approaches with examples drawn from the study of gas phase organic systems such as fullerenes, porphyrins, bacteriochlorophylls or nucleobases molecules. The perspectives and challenges that many-body perturbation theory is facing, such as the role of self-consistency, the calculation of forces and potential energy surfaces in the excited states, or the development of embedding techniques specific to the GW and Bethe–Salpeter equation formalisms, are outlined.
Central to a large variety of fields and applications, such as photovoltaics , printed ‘flexible’ electronics  and light-emitting diodes [3,4], photochemistry  or photosynthesis , the study of the electronic and optical properties of organic systems are crucial issues that stand as a challenge to both theory and experiment. The typical timescale associated with electronic excitation processes and the localization of the optical excitations on a few molecules in an usually disordered environment (solvant, donor–acceptor interface, etc.) are challenging difficulties even for modern experimental characterization tools, despite the tremendous progresses in, e.g. pump-probe time-resolved optical experiments. As a result, important questions such as the nature of the transport mechanisms of free carriers in organic systems (‘band’-like versus polaronic)  and the mechanisms of dissociation of the strongly bound photoexcited electron–hole pairs (excitons) in organic photovoltaic cells1 remain very controversial.
In such situations, which require information at the nanoscale, parameter-free (ab initio) simulations may stand as a valuable tool to guide the interpretation of experiments. Further, the idea that in silico computer simulations may guide experiments in designing novel molecules to improve, e.g. organic solar cell efficiency, namely computer-aided screening of novel systems, is an appealing direction that requires, however, to find a theoretical framework computationally manageable and showing ‘sufficient’ accuracy to be reliable. First-principles simulations may also be valuable to set up parametrized Hamiltonians capable of tackling the complexity of systems where several competing mechanisms are involved. Indeed, the various energy scales associated with (i) the inter-molecular hopping terms that govern the dispersion of the ‘bands’ in organic crystals, (ii) the strength of the electron–vibration interactions that distinguishes adiabatic versus vertical absorption and limits the mobility of carriers, or (iii) the effect of disorder on both ‘onsite’ and ‘offsite’ energy terms are very similar in organic crystals, so that establishing a hierarchy between the related phenomena in order to set up a perturbative treatment may prove difficult.
In this review, we present and compare mean-field approaches, such as density functional theory (DFT), Hartree–Fock and their hybrids, together with their time-dependent extensions, to quasi-particle Green's function approaches such as the GW formalism, aiming at predicting accurate electronic energy levels (ionization energy, electronic affinity, ‘band’ gaps, etc.) and the Bethe–Salpeter formalism designed to provide accurate optical (neutral) excitations. Concerning density functional theory, we focus more specifically on the electronic excitation properties as provided by the Kohn–Sham (KS) formalism. We will not be concerned in particular with ground state properties (cohesive energies, geometry, etc.). Even though frontier orbital energies in finite size systems can, in principle, be calculated by differences of ground state total energies between the neutral and charged systems (the so-called ΔSCF technique), such an approach does not apply to extended systems and does not allow to obtain the entire excitation spectrum and related eigenstates.
We illustrate our comparisons between DFT/time-dependent density functional theory (TDDFT) and GW/BSE on the basis of calculations performed on gas phase molecules, including fullerenes, acenes, porphyrins, chlorophylls and nucleobases. In the case of optical properties, much emphasis will be given to charge-transfer (CT) excitations, where the photoexcited electron is displaced from one site to another, within one molecule or between separated molecules. These special excitations are crucial to photovoltaic, photosynthetic or photocatalytic processes, but stand as a central limitation to standard TDDFT.
2. Electronic properties: generalized DFT–Hartree–Fock Kohn–Sham formalisms versus GW calculations
The quasi-particle formalism, namely the mapping of the true many-body problem onto a single (quasi-)particle framework, allows to draw fruitful correspondence between the KS approach within DFT, the Hartree–Fock Roothan's equation, and the self-energy formulation of the exchange and correlation contribution to the single-particle eigenvalue problem within many-body perturbation theory (MBPT). In all cases, one attempts to solve a one-body eigenvalue problem: 2.1where we introduce a general Σ(r,r′;E) self-energy operator for the exchange and correlation contribution. The self-energy operator can be shown to be, in general, non-local, energy-dependent and non-Hermitian, so that the corresponding eigenstates present an imaginary part interpreted as the lifetime of the quasi-particles with respect to electron–electron scattering.
(a) Generalized DFT–Hartree–Fock Kohn–Sham formalisms
In the standard KS framework with semilocal exchange-correlation functionals, Σ becomes the adiabatic (energy-independent) and local V XC([n];r) exchange-correlation potential, a functional of the charge density [n]. Within Hartree–Fock, the non-locality is recovered, thanks to the true exchange Fock operator with explicit summation over the Hartree–Fock self-consistent occupied eigenstates: ΣHF(r,r′)=−ρHF(r,r′)v(r,r′), where v(r,r′) is the bare Coulomb potential and the one-particle density matrix, but the energy-dependence is absent, as a signature of the instantaneous character of the exchange process. Energy-electron loss spectroscopy, namely the loss of energy of an electron going through a piece of matter with various impinging kinetic energy, clearly reveals that the effective interaction of an electron with its surrounding should depend on its energy.
It is now well documented that the KS band gap of semiconductors or insulators is significantly underestimated when using standard semilocal functionals, such as the popular local density approximation (LDA) [8–10] or the 1996 Perdew, Becke and Ernzerhof generalized gradient approximation (PBE) , as a result of the lack of discontinuity of the density-dependent exchange correlation potential upon addition (or removal) of an electron to the system [12,13]. A compilation of results for extended inorganic semiconductors and insulators can be found in reference  evidencing that indeed DFT–LDA KS band gaps are too small. On the contrary, by lack of correlation, Hartree–Fock produces much too large band gaps [15,16]. For sake of illustration, in the DFT–LDA and Hartree–Fock worlds, the band gap of bulk silicon is about 0.6 and 6.0 eV, respectively, to be compared with the 1.2 eV experimental value.
A similar behaviour can be observed for gas phase organic molecules. We provide in figure 1 the calculated and experimental energy gaps between the highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals for three paradigmatic molecules in organic electronics or photovoltaics, namely the C60 fullerene, the tetraphenylporphyrin (H2TPP) and pentacene. The DFT–LDA KS approach dramatically underestimates the energy gap, whereas the Hartree–Fock overestimates it. This is consistent with what is known of semiconductors, even though in the present case the DFT–LDA failure is much larger in percentage, whereas the Hartree–Fock performs somehow better, as a sign of reduced correlations in ‘few electrons’ systems when compared with bulk extended materials.
An important step was taken in 1993 by Axel Becke , who emphasized the importance of exact exchange for improving ‘thermochemical accuracy’ in DFT calculations, yielding several schemes to mix standard (semi)local exchange-correlation functionals with some percentage of exact exchange. As such, the exchange-correlation potential is no longer a pure density functional, but contains as well an orbital-dependent contribution leading to a generalized KS (GKS) framework [22,23]. Such a proposition developed quickly and became extremely popular in the study of organic systems with the development of hybrid functionals such as the famous Becke 3 parameters B3LYP functional . Even though designed to improve ground state properties (binding energy, bond lengths, etc.), such a mixing somehow also improves the description of the electronic properties. From a pragmatic point of view, this was to be expected because standard DFT with semilocal functionals underestimates band gaps, whereas Hartree–Fock overestimates them. More fundamental arguments can be invoked relying on two very desirable properties verified by the Hartree–Fock formalism, namely (i) the absence of self-interaction, thanks to its cancellation between the Hartree and Fock exchange terms, and (ii) the correct long-range (LR) (−1/r) behaviour of the effective potential seen by an electron far away from a molecule or a surface [24–26].
Despite such desirable features, the popular B3LYP functional, with the mixing of 20% of exact exchange, does not allow to obtain accurate electronic properties in the case of molecular systems as evidenced in figure 1. The HOMO–LUMO gap remains too small by several electron volts. Increasing the amount of exact exchange would certainly open the gap, but would degrade the ground state properties. It is clearly a difficult task to improve with a single-parameter both the ground state and the excited-state properties. Further, even for structural properties in the excited states, the choice of the percentage of exact exchange remains a difficult issue, and may depend not only on the system of interest, but also on the property to be described. As a recent important example in the field of organic chromophores, it was shown that a proper description of the geometry, and related optical absorption spectra, of various isomers of the famous retinal protonated Schiff base (RPSB) chromophore  requires to include as much as 54% of exact exchange, as incorporated in the M06-2X hybrid functional , that is a much larger percentage of exact exchange than the original B3LYP functional.
To proceed further, we note that the correct LR behaviour in vacuum of the effective exchange-correlation potential would require 100% of exact exchange, a dramatically too large quantity at short distance that would spoil the delicate cancellation of errors between local correlation and local exchange within the DFT framework. As a means to provide more flexibility in designing improved functionals, Savin and co-workers [29,30]2 introduced range-separated hybrids (RSHs), thanks to a partitioning of the true (1/r) Coulomb potential into an short-range (SR) erfc(μr)/r and an LR erf(μr)/r contribution, with (erf) the error function and (erfc) its complementary. This decomposition allows to properly introduce an SR and LR contribution to the exchange energy, namely for the LR Fock operator: . The local exchange (and correlation) contributions are treated with standard (semi)local functionals of the density. Tuning the (μ) parameter, with typical values of 0.2–0.5 bohr−1, allows to adjust the typical decay length of the local part.
The range-separation scheme was further reformulated , and an interesting solution to tune the range-separation (μ) parameter for optimal electronic properties was proposed by Baer, Neuhauser and Livshits (BNL) [32–34]. The starting observation is that ionization energies and electronic affinities calculated, thanks to a ΔSCF approach, namely a difference of ground state total energies for the neutral and charged molecule, are usually in reasonable agreement with experiment. As such, the range-separation parameter can be tuned, e.g. to optimize the ionization potential by enforcing the relation: εμHOMO=E(N;μ)−E(N−1;μ), where, for example, E(N;μ) are the DFT RSH ground state total energies of the system with N electrons obtained with a fixed (μ) parameter. A crucial observation is that one can find an optimal μ value that provides a good description of both the ionization energy and electronic affinity, namely a good HOMO–LUMO gap. This was illustrated in the case of a large set of organic molecules , and the results are reported in figure 1 in the case of C60 and H2TPP (‘OT-BNL’ values). The HOMO–LUMO gaps are starting to reach a few tenth of an electron volt accuracy, with the caveats that the range parameter needs to be tuned for each molecule, and that the underlying ΔSCF technique cannot be applied to extended systems. However, as shown below, such a functional, ‘educated’ on the HOMO–LUMO gap, has been also shown to provide remarkably accurate optical properties, as a signature of the well groundedness of this pragmatic scheme.
Besides the question of optimally tuning the range-separation parameter, several variations have been proposed around the long-range corrected (LC) scheme, either by changing the (erf) function with another decaying function, such as the exponential function, mimicking the screened (Yukawa) Coulomb potential in metals (e.g. ), or by providing more flexibility in the range decomposition such as in the three-parameters CAM-B3LYP functional introduced by Yanai et al.  with an LR contribution to the bare Coulomb potential reading . This generalized form allows to better tune the amount of LR exchange and to keep some percentage of exact exchange at SR. The original parametrization retained 19% of exact exchange at SR, close to the B3LYP percentage, and 65% at LR with an attenuation parameter μ of 0.33, offering a good compromise to describe both atomization energies at the B3LYP level of accuracy and the proper LR behaviour of the exchange-correlation potential crucial for excited-state properties such as Rydberg or CT optical excitations as described below in §a.
Before turning to the many-body perturbation theory treatment, we briefly introduce another rationale to design functionals with accurate electronic properties. A generalization of the KS scheme to fractional occupancy of the electronic energy levels quickly leads to the observation that the total energy derivative (δE/δN) should be constant for any fractional value of the number N of electrons between two successive integer values, displaying further a discontinuity at integer values . In great contrast, the standard density-dependent semilocal functionals within DFT are convex functionals of N with no discontinuity of (δE/δN) at integer values. On the contrary, the Hartree–Fock total energy is concave. Such opposite behaviour clearly stand again as a justification for generalized KS approach mixing both local density-based functionals with the bare exchange operator. These observations have been recently used as a guideline to prepare hybrid functionals imposing the linearity of the total energy with respect to fractional occupancy, yielding much improved KS single-particle energies . Alternatively, ad hoc orbital-dependent corrections to existing functionals can be imposed in order to enforce the piecewise linear behaviour of the total energy in between integer values of N. Such a scheme was for instance demonstrated by Dabo et al.  to lead, for a large set of organic molecules, to ionization potentials and HOMO–LUMO gaps in good agreement with experiment as evidenced in figure 1 (α-NKC0 results).
(b) The many-body Green's function GW formalism
We now turn to the expression of the self-energy within many-body perturbation theory and the GW approximation in particular. We will not attempt here to re-derive the expression for the self-energy within the so-called GW technique [38–41]. We just schematically sketch the principal ideas associated with the underlying Green's function formalism  and the functional derivative approach  to many-body perturbation theory, namely
— the (charged) excitation energies associated with adding or removing an electron from the system can be properly defined as the poles in the energy representation of the one-particle time-ordered G(1,2) Green's function which measures the amplitude of probability to find an electron (a hole) in position and time (2=r2,t2) after its introduction in the system at (1=r1,t1).
— the equation of motion-, or time-dependent Schrödinger equation for operators, that may be used to determine G(1,2), involves a two-body Green's function G2(1,2,3,4), initiating a hierarchical chain of increasing order Green's function with increasing complexity.
— fortunately, the needed two-body G2 operator can be expressed as a function of the derivative of the one-body G(1,2) Green's function with respect to an external perturbation, allowing to introduce quantities familiar within linear response theory, namely the susceptibility P(1,2) of the system, that relates the change of the charge density in (1) with respect to an external perturbation acting in (2), the energy-dependent dielectric function ϵ(1,2) and the related screened Coulomb potential W(1,2).
The quantities defined here above are however not given in closed form, but are related by a set of coupled (self-consistent) equations, traditionally labelled Hedin's equations  in condensed matter physics. These relations read where v(12)=v(r1,r2)δ(t1−t2) is the bare Coulomb potential, and Γ(34;2) the so-called three-body vertex correction. Such a set of equations can in principle be solved iteratively, starting from a zeroth-order system where the self-energy is zero, namely the Hartree mean-field solution, yielding to first-order: Γ(12;3)=δ(12)δ(13). This simple approximation for the vertex correction yields the famous GW approximation for the self-energy, written here in the energy representation3 where we have introduced the zeroth-order one-body (εn,ϕn) mean-field eigenstates. P0(r,r′;ω) is the non-interacting polarizability (the fi/j are the occupation factors). The summations over occupied and empty states lead to an O(N4) scaling for GW calculations with respect to system size, a scaling larger than the standard O(N3) scaling for DFT calculations with (semi)local functionals.
In practice, the mean-field starting point is never the Hartree solution, but more traditionally the DFT KS eigenstates which represent in general the ‘best available’ mean-field starting point. This leads to the standard ‘single-shot perturbative’ G0W0 treatment where the exchange-correlation contribution to the DFT KS eigenvalues is replaced by the GW self-energy operator expectation value onto the ‘frozen’ KS DFT eigenstates, namely Taking, for example, the LDA to the exchange-correlation potential vXC,DFT leads to the so-called G0W0@LDA scheme, the most common approach for GW calculations in solids.
While GW calculations for periodic bulk organic systems [44–46] or polymers with short repeat unit [47–50] appeared more than a decade ago, the case of large gas phase organic molecules such as fullerene or porphyrins could only be tackled in the past few years [19,51–61] owing to the cost of performing GW calculations with the traditional implementation for solids relying mostly on planewaves with periodic boundary conditions. Besides the large number of atoms (a few dozens) contained in such systems, it is to be noted that large unit cells are required to avoid cell–cell interactions in periodic boundary conditions, because GW calculations are concerned with charged excitations, leading to very slow convergence rate with respect to unit cell size. Several developments were recently introduced to significantly reduce the cost of GW calculations for organic systems, from the development of localized-basis implementations [19,47,56,62–65] to techniques allowing to bypass the summation over empty states [46,64,66–69].
As indicated in figure 1 for our three test molecules, the standard G0W0@LDA calculation significantly improves the quasi-particle properties, leading to errors on the HOMO–LUMO gap reduced from 3.4 eV within DFT–LDA (mean absolute error) to 0.7 eV within G0W0@LDA. However, the G0W0@LDA approaches still underestimates the HOMO–LUMO gap, by as much as 0.9–1.0 eV in the case of pentacene. These results point to the quality of the starting point KS states when performing single-shot G0W0 calculations. In the case of organic molecules, the dramatically too small value of the DFT–LDA HOMO–LUMO gap leads to overscreening when building the polarizability and the resulting screened Coulomb potential, explaining the too small G0W0@LDA HOMO–LUMO gaps.
Several solutions to that problem have been proposed. A first cure consists of changing the starting point by using hybrid functionals instead of DFT–LDA [59,60,70]. Another long-standing line of work consists of self-consistent GW calculations where the screened Coulomb potential W and/or the Green function G are improved by re-injecting solutions from the previous GW iteration , a scheme that has been recently explored for isolated atoms [71,72] and small molecules [59,71,73]. As a simple approach to self-consistency, we present in figure 1 the results obtained by self-consistently updating the eigenvalues in G and W, whereas ‘freezing’ the starting DFT–LDA KS wave functions, a scheme we label the ev-GW@LDA approach. This leads to a significant improvement of the ionization potential and energy gap, as recently shown for a larger set of organic molecules [19,45]. In particular, the HOMO–LUMO gap of C60 and H2TPP are now within one-tenth of an electron volt from the experimental value.
We conclude this section by briefly emphasizing two important issues that have recently emerged in the GW study of organic molecular systems. A first important observation is that the GW approach leads not only to a much better description of the ionization energy and electronic affinity, but also corrects as well the ordering of levels which, for organic molecules presenting π and σ states of different nature and localization, can be wrong within LDA or PBE. This was, for example, demonstrated in the case of DNA and RNA nucleobases [53,54] where the ev-GW@LDA approach was shown to lead to an excellent agreement with high-level coupled-cluster quantum chemical approaches . This is illustrated in figure 2 where the absolute position and ordering of the states at the top of the occupied manifold of gas phase cytosine is represented at several theoretical levels. While DFT–LDA predicts the HOMO level to be a σO state strongly localized on the oxygen atom, the ev-GW@LDA and high-level multi-detrimental quantum chemistry techniques correctly describe the HOMO level as a more standard π-state lying ∼ 2.5 eV below the KS LDA value. As a result, the ev-GW@LDA ionization potential lies within 0.2 eV of the experimental value. The B3LYP functional corrects the ordering of the levels, but hardly improves the ionization energy. However, a good agreement with the ev-GW@LDA values can also be obtained within the framework of tuned ranged-separated hybrids .
Another important issue is the observation that the calculation of electron–phonon coupling constants, namely the variation <δεn/δRν> of the electronic energy levels with respect to an ionic motion along a ‘Rν’ vibrational eigenmode, can change significantly from DFT–LDA to GW. It was shown in particular in the case of C60 that the electron–phonon coupling potential to the LUMO level, a quantity that used to be quite studied in the framework of fullerides superconductivity, was 50% larger within GW, leading to a much better agreement with experimental values . The importance of the electron–phonon coupling strength in understanding the ‘adiabatic’ band structure  and mobility of carriers  in organic systems certainly call for a careful re-evaluation of previous studies at the DFT–LDA level. The same conclusions apply to the field of phonon mediated superconductivity.
3. Time-dependent density functional theory versus Bethe–Salpeter equations
While GW calculations aim at obtaining accurate quasi-particle energies, such as the electronic affinity, ionization potential, and more generally the entire ‘band structure’ of a given system, to be compared, e.g. to direct and inverse photoemission experiments, the Bethe–Salpeter formalism tackles the problem of the neutral optical excitations, namely excitations where the electron does not leave the system and interacts through the (screened) Coulomb potential with the hole left in the occupied manifold. As such, the Bethe–Salpeter equation (BSE), originally derived in the 50s  and adapted in the mid-1960s by Sham and Rice in the context of condensed matter physics [78–80], is a two-body electron–hole eigenvalue problem, very much as the proton–electron hydrogenoid textbook exercise, but in an environment that provides screening through Coulomb and Pauli repulsion between electrons. As another familiar reference for solid-state physicists, it can be considered as a generalization of the Elliott  (or Mott–Wannier) approach to delocalized excitons in semiconductors. Three independent ab initio implementations were developed simultaneously in 1998 [82–84]. Because BSE calculations require as an input the GW quasi-particle energy levels,4 ‘GW-BSE’ calculations on molecular systems such as fullerenes or porphyrins have also emerged rather recently [88–98], after pioneering studies on small molecules [99,100] or bulk organic semiconductors [101–103].
The BSE formalism can be straightforwardly compared with TDDFT [104–106] in the so-called Casida's formulation , which recast the TDDFT problem as a similar eigenvalue problem in the electron–hole two-body basis, namely 3.1where the indexes (i,j) and (a,b) indicate the occupied and virtual orbitals, and (re,rh) the electron and hole positions, respectively. In this block notation, the vector [ϕa(re)ϕi(rh)] represents all excitations (note, for example, that ϕa(re) means that an electron is put into a virtual orbital), whereas the vector [ϕi(re)ϕa(rh)] represents all de-excitations. As such, R (R*) describes the resonant coupling between electron–hole excitations (de-excitations), while the off-diagonal blocks C and C* account for non-resonant coupling between excitations and disexcitations. The TDDFT and BSE resonant parts can be directly compared, taking for TDDFT its expression for a general hybrid range-separated functional: where is the LR part of the Coulomb potential in the general CAM-B3LYP range-separated formulation. KDFT=(v+fXC) is the sum of the bare Coulomb potential and the exchange-correlation (semi)local kernel fXC. We use the notation 〈⋯ 〉 for the double integral. The middle (last) line gathers terms with occupied and virtual orbitals taken at different (identical) integration variables. For isolated systems, the wave functions can be taken to be real. The resonant terms have been given explicitly here above for singlet excitations and we have considered spin-unpolarized systems.
A clear difference between TDDFT and GW-BSE stems from the diagonal part (first line) of the resonant contributions. The BSE formulation starts from the quasi-particle energy differences as given within the GW formalism, whereas the TDDFT formulation uses the KS eigenvalues. As emphasized here above, the GW HOMO–LUMO gaps generally agree well with the experimental value, whereas their KS analogues are usually significantly underestimated. Consequently, the BSE equations lend support to the intuitive physical picture that the attractive electron–hole interaction stabilize excitons, leading to an optical gap smaller than the quasi-particle (photoemission) HOMO–LUMO one. On the contrary, the interaction terms in the TDDFT formalism contribute to open the usually too small KS HOMO–LUMO gap in order to obtain optical excitation energies. A second important remark is that the long-range vLRαβμ contribution within TDDFT plays the role of the screened Coulomb potential W that appears in the BSE framework. We will return to that point here below.
4. Optical properties and charge-transfer excitations
Despite the differences highlighted here above, the BSE and TDDFT formalisms usually agree very well for excitations where the hole and the electron strongly overlap spatially. Such transitions are often labelled ‘Frenkel excitations’. This concerns in particular intramolecular excitations in molecules such as the one studied here above (fullerenes, acenes, porphyrins, etc.) where the HOMO and LUMO eigenstates are uniformly distributed. The literature emphasizing the accuracy of TDDFT for predicting correctly the optical spectra of such systems, even in cases where standard local and adiabatic semilocal kernels are used, is too large to be reported here, and we refer the reader to a recent benchmark for a large number of functionals on a large set of reference molecular systems . One example is provided below.
We focus more specifically on the so-called CT excitations for which the photoexcited electron and the hole left behind are spatially separated. Such excitations are central to many applications in photovoltaics, photocatalysis, photosynthesis and biology. Further, it is now well-known that standard TDDFT dramatically fails to describe such excitations as explained here below. As such, these non-local excitations, as well as the so-called Rydberg transitions to high-energy delocalized weakly bound states, have been central to the motivations for developing range-separated hybrids (see e.g. the CAM-B3LYP seminal paper ).
(a) A historical problem: the zincbacteriochlorin–bacteriochlorin complex
One of the first systems revealing an important difficulty associated with TDDFT is a donor–acceptor organic complex mimicking photosynthetic processes in bacteria, namely the zincbacteriochlorin–bacteriochlorin (ZnBC–BC) complex, initiating several theoretical studies at various TDDFT or quantum-chemistry post-Hartree–Fock levels [108–112]. In such molecular dimers, the optical excitations can be classified into two categories, namely intramolecular excitations, where the promoted electron remains on the same molecule, and CT excitations for which the photoexcited electron ‘jumps’ from the donor to the acceptor. The ZnBC–BC atomic structure is provided as an inset in figure 3a.
Comparing now the various formalisms, it was shown that TD-LDA, TD-B3LYP and GW-BSE excitation energies associated with intramolecular transitions, such as the lowest ‘Q’ transitions, all agree with a maximum discrepancy of 0.2 eV (see table I in ). As emphasized here above, TDDFT with standard adiabatic semilocal kernels reproduces usually very well Frenkel excitations. On the contrary, CT excitations within TD-LDA were underestimated by about 1.6 eV when compared with GW-BSE. Performing TDDFT calculations with the global B3LYP hybrid and the range-separated CAM-B3LYP functional reduced the error to about 1 and 0.1 eV, respectively. Clearly, introducing a sufficient amount of exact exchange dramatically improves the description of CT excitations at the TDDFT level, without degrading the description of the Frenkel transitions.
Charge-transfer excitations verify a simple asymptotic behaviour in the so-called ‘Mulliken limit’  of a large distance D between the donor and the acceptor. Namely, the lowest lying CT excitation should converge to: ECT0(D)=EA(acceptor)−IP(donor)−1/D, where EA (acceptor) and IP(donor) are the acceptor electronic affinity and the donor ionization energy. The (−1/D) scaling with distance is just the electrostatic interaction (in atomic units here) between the hole left on the donor and the electron that has jumped onto the acceptor. It was shown very early  that within TDDFT using standard (semi)local exchange correlation kernels, the CT excitation energies were indeed underestimated by as much as several electronvolts, showing further a dispersion-less character, failing to follow the correct (−1/D) scaling behaviour.
To understand this behaviour, it is instructive to consider the projection of the resonant part R of the TDDFT or BSE two-body Hamiltonian into the subspace spanned by the HOMO(D) and LUMO(A) levels of the donor (D) and the acceptor (A) in the limit of non-overlapping donor and acceptor wave functions. We note Rhl,hl such a scalar energy. In the long distance limit, the previous equations for Rai,bj yield and noting that in this limit . We introduced the macroscopic (LR) dielectric constant (ϵM).
The standard TDDFT approach, namely a TDDFT based on a density-dependent local XC functional, can be recovered by setting (α) and (β) to zero. Inspection of the TDDFT asymptotic limit allows then to understand the lack of dispersion with donor/acceptor distance of the CT excitation, and the fact that the CT excitation energies reduced to the (too small) HOMO–LUMO KS gap. In the case of hybrid functionals without any range separation (namely β=0), then the scaling reduces to (−α/D), with e.g. (α=0.2) in the B3LYP formulation, allowing to partially recover the correct dispersion. In the original CAM-B3LYP formulation , (α+β) was set to 0.65, starting to approach the ideal asymptotic behaviour. The original LR-corrected (LC) range-separated hybrid (α=0,β=1) allows to recover exactly the correct asymptotic (−1/D) behaviour. Combined further with, e.g. the Baer–Neuhauser–Livshits (BNL) idea to tune the range-separation parameter in order to obtain the correct HOMO–LUMO gap (see §2a), both the correct (−1/D) scaling and asymptotic quasi-particle HOMO–LUMO gap can be recovered.
Clearly, in the Bethe–Salpeter formalism, the dielectric constant ϵM for donor and acceptor fragments separated by vacuum reduces to one, so that the correct (−1/D) limit is recovered, converging towards the correct HOMO–LUMO gap described at the GW level. This is illustrated in figure 3a where the dispersive character of the CT excitations is clearly reproduced. Further, in the limit of isolated BC or ZnBC molecules, the GW-BSE intramolecular (Frenkel) ‘Q’ transitions are found to agree within 0.1 eV when compared with experiment . As a result, both intramolecular and CT excitations are correctly reproduced. Similar results were obtained in the case of donor molecules (the coumarins) of interest for dye-sensitized solar cells (DSSCs), with an excellent description of excitations showing mixed Frenkel and CT character , providing good agreement with tuned range-separated hybrid calculations . This ability to describe correctly excitations with various spatial character is a crucial feature when attempting to understand the evolution of Frenkel excitations into resonant ‘hot’ CT states which are believed to be the initial step for exciton dissociation in organic solar cells . As another important remark, the comparison between formalisms shows that the CAM-B3LYP TDDFT (α,β) parameters would have to be re-adjusted in another screening environment (solvant, bulk, etc.). In the parameter-free Bethe–Salpeter formalism, the dielectric function and the screened Coulomb potential automatically account for the proper screening properties.
(b) Donor–acceptor complexes and comparison with experiment
An interesting family of systems allowing to further illustrate the above-mentioned considerations is a set of donor–acceptor gas phase complexes composed of the tetracyanoethylene (TCNE) acceptor molecule with acenes (benzene, naphthalene, etc.) and their derivatives. The interest of such systems is that the low-lying optical excitations are clear CT states well characterized by available gas phase experimental data , triggering several global and range-separated hybrid TDDFT  and GW-BSE studies [94,97,117]. The HOMO and LUMO states for the TCNE–anthracene complex are represented in figure 3(b, inset) revealing clearly the donor–acceptor character of such molecular dimers.
As shown in figure 3b, the TD-B3LYP calculations (blue up triangles), despite the 20% of exact exchange, dramatically underestimate the excitation energies when compared with the experimental data (red dots). The G0W0-BSE calculations (yellow diamonds) provide better results, but still underestimate the excitation energies, owing to the too small G0W0@LDA HOMO–LUMO gaps (see §2b). On the contrary, both the range-separated ‘optimally tuned’ hybrid TD-BNL calculations (green squares) , and the GW-BSE calculations (black diamond), with self-consistency on the eigenvalues , provide an excellent agreement with experiment. As emphasized above, the range-separation (μ) parameter involved in the range-separated hybrid TDDFT calculation has been tuned to reproduce the correct HOMO–LUMO gap, following the BNL prescription. On the contrary, the GW-BSE approach is parameter-free.
After decades of expertise in applying the GW-BSE formalism to inorganic semiconducting or insulating systems, there is an emerging line of work devoted to extending such techniques to organic molecular systems for applications in electronics, photovoltaics, photocatalysis and biology. As emphasized here above, excellent results can be obtained in describing the electronic charged and neutral (optical) excitation energies of a large variety of systems. While the GW-BSE formalism cannot be claimed to provide more accurate results than DFT or TDDFT calculations with ‘well-parametrized’ range-separated hybrids, a clear advantage of the many-body perturbative Green's function approach is that it is parameter-free and has been shown to provide reliable results in a large variety of inorganic or organic systems, with a metallic or strongly insulating character, in finite size or extended systems, demonstrating similar accuracy (a very few tenth of an eV) for local or non-local excitations. This is an important feature for upcoming studies of hybrid systems mixing, e.g. organic molecules with extended metallic or insulating inorganic substrates (TiO2 in DSSCs, nanotubes, graphene, metallic electrodes, etc.).
Even though GW-BSE calculations on systems comprising a few hundred atoms are nowadays feasible , offering as a matter of fact the same O(N4) scaling with size than DFT calculations with hybrid functionals, many-body perturbation theory calculations will certainly remain computationally more expensive that (TD)DFT implementations. However, the recent developments in the field (localized basis formulations, resolution of the identity techniques, removal of the summation over empty states, etc.) clearly allow to provide reference calculations on close to realistic systems. As such, using GW-BSE calculations to ‘educate’ improved hybrid functionals is certainly a valuable direction of work: clearly, the development of ‘locally screened’ range-separated hybrids, namely with a space-dependent Coulomb attenuation parameter , already points in the direction of mimicking within DFT the behaviour of the screened Coulomb potential W at the heart of the GW and Bethe–Salpeter formalisms.
Despite such recent progress in performing many-body perturbation theory for organic systems, there is still much to be done to be able to provide insights into complex mechanisms such as exciton dissociation at disordered bulk donor–acceptor interfaces in photovoltaic cells, or bond breaking in the excited state in photocatalytic applications. A first important goal along the line of applying MBPT to more realistic systems is the development of specific ‘embedding’ techniques in order to be able to treat in an accurate but affordable way the effect of the environment (solvant, disordered crystalline phase, etc.). Adapting the discrete or continuous  (and  for a recent review) polarizable models already available to TDDFT is certainly an important direction in the field.
A further important challenge in the GW-BSE framework is the calculation of gradients in the excited states. A specific case is that of the forces acting on the ions upon photoexcitation, a key ingredient in a field such as photochemistry. Even in less extreme applications not involving bond breaking, such as the understanding of the absorption and fluorescence spectra in organic molecules, the knowledge of the relaxation energy upon photoexcitation, and of the zero-point motion energy contribution, allowing to differentiate vertical versus adiabatic transition energies, has proved to be crucial in order to compare theory with experiments . Calculations of analytic gradients in the excited states within TDDFT are now routinely implemented in most codes. Even though the quality of the obtained potential energy surface and of the related relaxed molecular geometry in a given excited state may drastically depend on the chosen functional , the formalism for calculating the gradients, without using expensive finite difference techniques, is now well documented [122,123]. Such analytic gradients do not exist at the GW/BSE level and work along that line is required .
C.F. acknowledges a Ph.D. grant from the French CNRS and CEA agencies and X.B. partial support from the French Research Agency under contract no. ANR 2012-BS04 PANELS. Computer time for part of the work described in this work has been provided by the French national GENCI supercomputing center CEA-Curie through national (grant no. i2012096655) and European PRACE (grant no. 2012071258) projects. Part of the work presented in this review has been performed in collaboration with Valério Olevano, Prs. Erich Runge (Ilmenau, Germany) and Michel Côté (Montreal, Quebec). They are warmly acknowledged.
One contribution of 8 to a Theme Issue ‘Density functional theory across chemistry, physics and biology’.
↵1 See for example the Correspondence section of Nature Mater.12 (2013) 593–595 and the controversy about the role of ‘hot’ versus ‘cold’ CT excitons in the dissociation mechanism(s) of photogenerated bound Frenkel excitations.
↵2 This partitioning was originally introduced to combine post-Hartree–Fock many-body wave functions calculations, such as CI for LR exchange and correlation effects, with density functional theory for short-range (SR) contributions.
↵3 Time homogeneity allows to conclude that the defined operators depend only on the difference (t2−t1) leading via Fourier transform to the energy dependence. A proper mathematical treatment of such Fourier transforms leads to introduce the small infinitesimals (0+) present in some of the formulae.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.