Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

The Bethe–Salpeter formalism with polarisable continuum embedding: reconciling linear-response and state-specific features

Ivan Duchemin *a, Ciro A. Guido bc, Denis Jacquemin b and Xavier Blase *d
aUniv. Grenobles Alpes, CEA, INAC-MEM, L_Sim, F-38000 Grenoble, France. E-mail:;
bLaboratoire CEISAM – UMR CNR 6230, Université de Nantes, 2 Rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France
cLaboratoire MOLTECH – UMR CNRS 6200, Université de Angers, 2 Bd Lavoisier, 49045 Angers Cedex, France
dUniv. Grenobles Alpes, CNRS, Institut Néel, F-38042 Grenoble, France

Received 1st February 2018 , Accepted 2nd April 2018

First published on 5th April 2018

The Bethe–Salpeter equation (BSE) formalism has been recently shown to be a valuable alternative to time-dependent density functional theory (TD-DFT) with the same computing time scaling with system size. In particular, problematic transitions for TD-DFT such as charge-transfer, Rydberg and cyanine-like excitations were shown to be accurately described with BSE. We demonstrate here that combining the BSE formalism with the polarisable continuum model (PCM) allows us to include simultaneously linear-response and state-specific contributions to solvatochromism. This is confirmed by exploring transitions of various natures (local, charge-transfer, etc.) in a series of solvated molecules (acrolein, indigo, p-nitro-aniline, donor–acceptor complexes, etc.) for which we compare BSE solvatochromic shifts to those obtained by linear-response and state-specific TD-DFT implementations. Such a remarkable and unique feature is particularly valuable for the study of solvent effects on excitations presenting a hybrid localised/charge-transfer character.

1 Introduction

The exploration of the excited-state (ES) properties of chemical systems certainly stands as a central question in theoretical chemistry. Indeed, ES phenomena govern many applications such as solar energy conversion, photocatalysis, light-emission or optical information storage. Further, while experimental characterisations can provide reference absorption and/or emission spectra, they are less suited to obtain some key information, e.g., ES geometries, nature of the excitation (localised, charge-transfer, etc.) or time evolution of hot electrons. Such a need for quantum mechanical formalisms allowing us to study realistic systems certainly explains the formidable popularity of time-dependent density functional theory (TD-DFT)1,2 that can be used to study the optical properties of systems comprising up to a few hundred atoms, thanks to a (formal) image file: c8sc00529j-t1.tif scaling with system size. Further, the availability of analytical TD-DFT derivatives3–7 together with the extension of efficient continuum models, such as the polarisable continuum model (PCM),8,9 to TD-DFT10–15 has dramatically helped in bridging the gap between quantum simulations and realistic systems, by respectively allowing us to explore ES potential energy surfaces and to take into account the impact of the surroundings. In TD-DFT, the coupling with the PCM was initially performed within a linear-response (LR) formalism,10,11 that is, using the electronic transition density for including solvent effects. While such a LR model is generally accurate for describing the local ES, it is less suited for the charge-transfer (CT) ES, in which a large reorganisation of the electron density occurs. To tackle such an ES, state-specific (SS) PCM-TD-DFT models, in which the solvatochromic effects depend on the total electronic density of the ES, have been developed.7,13–15 At this stage, let us point out that, in using these SS models, one can encounter some cases for which the exact details of the chosen SS approach as well as the selected exchange-correlation functional have a very large impact on the results, especially when self-consistent iterative methods are selected.7,16–18 In addition, as first pointed out by Corni et al.,19 who used a simple formal model explicitly including two states for the solute and two solvent macrostates, there is a need to simultaneously account for both LR and SS effects. However, to date, only the ad hoc sum of both LR and SS terms, determined in the context of a corrected linear response (cLR) approach,13 was proposed in a TD-DFT context.20 Alternatively one can turn towards single-reference electron-correlated wavefunction approaches, such as ADC(2), CC2, CCSD or SAC-CI, that have all been coupled to continuum models,21–30 but these models imply a significantly increased computational effort compared to TD-DFT. In this framework, we underline that the importance of the inclusion of both LR and SS effects was also clearly underlined by Lunkenheimer and Köhn in their work describing the coupling of the ADC(2) theory to a continuum approach of solvation effects.26

As another alternative to TD-DFT, the Bethe–Salpeter equation (BSE) formalism31–35 has been recently experiencing a growing interest in the study of molecular systems due to its ability to overcome some of the problems that TD-DFT is facing, including charge-transfer36–44 and cyanine-like45,46 excitations, while preserving the same image file: c8sc00529j-t2.tif scaling in its standard implementations. Extensive benchmark studies on diverse molecular families have been performed,47–53 demonstrating that excellent agreement with higher-level many-body wavefunction techniques, such as coupled-cluster (CC3) or CASPT2, could be obtained for all types of transitions, provided that they do not present a strong multiple-excitation character. Singlet–triplet transitions constitute the only notable exception as they may present the same instability problems with BSE and TD-DFT.51,52 We have recently reviewed the differences between BSE and TD-DFT formalisms in a chemical context, and we refer the interested reader to that original work for more details.35

As compared to TD-DFT, the BSE formalism relies on transition matrix elements between occupied and virtual energy levels calculated at the GW level, where G and W stand for the one-body Green's function and the screened-Coulomb potential. These GW energy levels, including HOMO and LUMO frontier orbital energies, were shown to be in much better agreement with reference wavefunction calculations as compared to standard DFT Kohn–Sham (KS) eigenvalues. Namely, GW HOMO and LUMO energies can be directly associated with the ionisation potential (IP) and electronic affinity (AE). Further, the TD-DFT exchange-correlation kernel matrix elements in the occupied-to-virtual transition space are replaced within the BSE formalism by matrix elements involving the screened Coulomb potential interaction between the hole and the electron.

The coupling of the GW and PCM formalisms was recently achieved,54 within an integral equation formalism9 (IEF-PCM) implementation including the so-called non-equilibrium (neq) effects related to the separation of the solvent response into its “fast” electronic and “slow” nuclear contributions. The electron and hole polarisation energies, namely the shifts of the ionisation potential and electronic affinity from the gas phase to solution, were very well reproduced with GW taking as a reference standard ΔSCF + PCM calculations where total energy calculations of the solvated ions and neutral species were performed at the DFT or CCSD levels. Similar studies were also performed adopting discrete polarisable models to study organic semiconductors and complexes of interest for optoelectronic applications.55–58 As central (GW + PCM) formalism features, we underline that the polarisation energies for all occupied (Pi+) and virtual (Pa) energy levels can be calculated, not only those associated with frontier orbitals, and that these polarisation energies are SS. Such polarisation energies, of the order of an eV in the case of water solvated nucleobases,54 dramatically reduce the HOMO–LUMO gap as compared to its gas phase value (see Fig. 1).

image file: c8sc00529j-f1.tif
Fig. 1 Schematic representation of the state-specific contribution to solvatochromic shifts within the (BSE/GW + PCM) formalism. The GW energy levels are first renormalised by the SS solvent-induced polarisation energies Pa(ε) and Pi+(ε) for electrons (occupied) and holes (virtual), respectively. Further, the strength of the screened electron–hole interaction Weh is also renormalised, namely reduced by a quantity that we label here Peh(ε), thanks to the additional screening provided by the PCM. Ground-state polarisation effects associated with the slow (ε0) degrees of freedom are incorporated at the initial (DFT + PCM) level (not represented here).

In the present study, we demonstrate that the BSE formalism combined with the PCM in a non-equilibrium formulation intrinsically combines the LR and SS solvatochromic contributions associated with the effect of the polarisable environment. Such a remarkable feature hinges in particular on the proper inclusion of dynamic polarisation energies not only for the occupied and virtual energy levels, calculated within the GW formalism, but also for the screened Coulomb electron–hole interaction. This renormalisation by polarisation of energy levels and electron–hole interactions leads to the SS contribution to solvatochromic shifts, in addition to the familiar transition matrix elements stemming from the LR contributions. This simultaneous account of both LR and SS effects was, to the best of our knowledge, only achieved using an a posteriori sum of the two terms in the context of the computationally more expensive ADC(2)/COSMO approach,26 and more recently with a TD-DFT/discrete polarizable scheme.20 Here, the obtained (BSE + PCM) solvatochromic shifts are computed for paradigmatic transitions in acrolein, indigo, p-nitro-aniline (PNA), a small donor–acceptor complex and a solvatochromic probe (see Fig. 2), and are compared to the sum of the shifts obtained at the TD-DFT level within the standard LR and SS implementations of the PCM, respectively.

image file: c8sc00529j-f2.tif
Fig. 2 Schematic representation of the molecules and complexes studied: (a) acrolein (A), (b) indigo (I), (c) para-nitroaniline (PNA), (d) “twisted” para-nitroaniline (PNAperp), (e) the donor–acceptor benzene–tetracyanoethylene (B-TCNE) complex and (f) the 4-nitropyridine N-oxide solvatochromic probe.

2 Formalism

We will not detail here the BSE formalism31–34,59–61 for which several reviews are available.35,62,63 We restrict ourselves to highlight the important points that allow understanding the merging with the PCM formalism. Since the BSE approach requires as an input the screened Coulomb potential W and accurate occupied/unoccupied {εGWi/a} energy levels calculated at the GW level, we start by a short introduction to the Green's function GW formalism62,64–67 from which the BSE approach derives.

2.1 GW formalism

The GW approach belongs to the family of many-body perturbation theories (MBPT) that takes as central quantity the Green's function (G) instead of, e.g., the charge density within DFT. In its time-ordered formulation, the Green's function G reads:
image file: c8sc00529j-t3.tif(1)
where {εi,εa} are the “true” electronic energy levels as defined experimentally in a direct/inverse photoemission experiment. The small positive infinitesimal η controls the proper analytic properties of G in the complex energy plane.67 More precisely, we define: εa = E(N + 1, a) − E(N, 0) for “virtual” energy levels and εi = E(N, 0) − E(N − 1, 0) for “occupied” energy levels, where E(N + 1, a) is the total energy of the (N + 1)-electron system in its ath quantum state, E(N − 1, i) the total energy of the (N − 1)-electron system in its ith quantum state, and E(N, 0) the N-electron system ground-state energy. {ϕi, ϕa} are called “Lehman amplitudes”. It can be formally shown that G verifies a simple Dyson equation:
image file: c8sc00529j-t4.tif(2)
with, e.g., 1 = (r1,t1) a space-time coordinate. G0 is the independent-electron Green's function and the “self-energy” operator ΣHXC contains all electron–electron interactions (Hartree plus exchange and correlation). Plugging the expression for G into the Dyson equation leads to a familiar eigenvalue equation for the (photoemission) excitation energies defined here above, namely:
image file: c8sc00529j-t5.tif(3)
where the h^0 Hamiltonian contains the kinetic, ionic and classical Hartree operators. Such an equation is formally equivalent to the DFT KS equation, but with an exchange-correlation “potential” that is both non-local and energy-dependent.

While eqn (3) is exact, an expression for ΣXC should be defined. The GW formalism provides an approximation for ΣXC to first order in the screened Coulomb potential W, with:

image file: c8sc00529j-t6.tif(4)
image file: c8sc00529j-t7.tif(5)
image file: c8sc00529j-t8.tif(6)
image file: c8sc00529j-t9.tif(7)
where v(r,r′) is the bare Coulomb potential, χ0 the independent-electron susceptibility and W the screened Coulomb potential. EF is the energy of the Fermi level. {fi/j} are occupation numbers and the input {ϕn, εn} eigenstates are typically KS eigenstates that will be corrected within the GW formalism.

The GW eigenvalues have been shown in many benchmark studies on gas phase molecular systems to be significantly more accurate than KS or Hartree–Fock eigenstates, providing, e.g., frontier orbital energies within a few tenths of an eV with respect to reference CCSD(T) calculations.68–72 For the sake of illustration, our gas phase ionization potential (IP) and electronic affinity (EA) are 10.35 eV (−εGWHOMO) and 0.68 eV (+εGWLUMO) respectively for acrolein within our GW approach, comparable to 10.01 eV and 0.70 eV within CCSD(T) for the same atomic basis (cc-pVTZ) and the same geometry (see the ESI). In practice, GW calculations proceed traditionally by correcting input KS eigenstates, substituting the GW self-energy contribution to the DFT exchange-correlation potential:

image file: c8sc00529j-t10.tif(8)

The obtained εGWa/i eigenvalues and the screened Coulomb potential W serve as input quantities for the BSE excitation energy calculation.

2.2 Bethe–Salpeter equation

While TD-DFT starts from the evaluation of the density–density susceptibility χ(12) = ∂n(1)/∂Uext(2) that measures the variation of the charge-density with respect to an external local perturbation, the excitation energies within BSE are obtained through the poles of a generalized susceptibility L(1234) = ∂G(12)/∂Uext(34), where Uext(34) is a non-local external perturbation. Deriving the Dyson equation for G (eqn (2)) leads in particular to the introduction of the (∂ΣXC/∂G) derivative, the analog to the exchange-correlation kernel within DFT. In the standard BSE/GW formalism, we consistently assume the GW approximation for ΣXC. Expressing then the 4-point susceptibility L in the transition space between occupied and unoccupied one-body eigenstates yields a standard linear algebra eigenvalue representation:
image file: c8sc00529j-t11.tif(9)
where ΩBSEλ represents the BSE excitation energies. Xλ represents the components of the two-body electron–hole ψλ(re,rh) eigenstates over the ϕi(rh)ϕa(re) transition basis, where (i,a) index (occupied, virtual) eigenstates, while Yλ gathers the ϕi(re) ϕa(rh) de-excitation components. Such a linear algebra representation in the transition space resembles the standard TD-DFT formalism within the so-called Casida's formulation1 but with modified matrix elements. For the resonant block, one obtains:
ABSEai,bj = δabδij(εGWaεGWi) − 〈ab|W|ij〉 + 〈ai|bj(10)
ab|W|ij〉 = 〈ϕa(r)ϕb(r)W(r,r′)ϕi(r′)ϕj(r′)〉(11)
ai|bj〉 = 〈ϕa(r)ϕi(r)v(r,r′)ϕb(r′)ϕj(r′)〉(12)
where εGWa/i are the GW unoccupied/occupied energy levels and W the screened Coulomb potential. The (εGWaεGWi) energy differences replace the TD-DFT KS (εKSaεKSi) energy differences between virtual and occupied states, while the W matrix elements can be interpreted as electron–hole interaction terms through the screened Coulomb potential. {ϕi/a} are the starting KS one-body eigenstates that are not corrected in the GW implementation selected here (see below).

2.3 Merging with the PCM formalism

The merging of the GW formalism with the PCM is described in detail in ref. 54. The seminal point is that the solvent reaction field can be straightforwardly incorporated into the screened Coulomb potential W that accounts for both the polarisability of the solute and that of the solvent when combined with the PCM. Under the assumption that the solute and solvent eigenstates do not spatially overlap, it can be shown that the screened Coulomb potential W is changed from its gas phase value:
image file: c8sc00529j-t12.tif(13)
to its condensed phase expression:
W = + ṽχQM0W(14)
= v + PCM0(15)
where we have dropped the frequency and space variables in the second equation for conciseness. χQM0 is the gas phase independent-electron susceptibility of the quantum subsystem to be solvated, while χPCM0 is that of the PCM solvent. The modified potential is thus the bare Coulomb potential renormalized by the solvent reaction field, equivalent to:
= v + PCMv(16)
where χPCM is the full (interacting) susceptibility of the solvent. The quantity [PCMv] (r1,r2) represents the reaction field (vreac) generated in r2 by the PCM surface charges developed in response to a unity point-charge added in r1, with (r1,r2) located in the cavity carved in the solvent to accommodate the solute. Details of our implementation using localised Gaussian bases and the Coulomb-fitting resolution-of-identity (RI-V) approach can be found in ref. 54 for the merging of GW with the PCM and ref. 56 and 58 for the merging with a discrete polarisable model.

Importantly, in the present non-equilibrium formulation of the (BSE/GW + PCM) implementation, the reaction field that modifies W is associated with the “fast” electronic excitations, namely with the ε dielectric constant (e.g., ε = 1.78 in water) that is the low-frequency optical dielectric response. The inclusion of the “slow” response of the solute (given by, e.g., ε0 = 78.35 for water) is accounted for in the preliminary ground-state DFT + PCM(ε0) run that serves as a starting point for GW and BSE calculations, namely in the construction of χ0 and G in eqn (5) and (7). The flow of calculations can then be summarised as follows:

(1) Calculation of input KS {εn, ϕn} eigenstates within a DFT + PCM(ε0) scheme. These eigenstates contain ground-state solvation effects,

(2) Calculation of the “fast” reaction field vreac(ε) in an auxiliary basis representation (see ref. 54, 56 and 58) that is incorporated inside the screened Coulomb potential W following eqn (5) and (7),

(3) Correction of the {εn} KS eigenstates, that include “slow” polarisation effects, with the GW self-energy operator that contains the “fast” vreac(ε) reaction field in order to yield the GW + PCM neq {εnGW+PCM} energy levels,

(4) Resolution of the BSE excitation energy eigenvalue problem with the {εnGW+PCM} eigenstates and the W and potentials that include vreac(ε) = PCM(ε)v.

2.4 LR and SS solvatochromic contributions

The renormalisation of v into v + vreac(ε) and the resulting change ΔW, accounting for switching on the fast ε solvent response, leads to modifying eqn (10) as follows:
image file: c8sc00529j-t13.tif(17a)
+δabδijΔ(εGW/PCMaεGW/PCMi) − 〈abW|ij(17b)
where ABSE/PCM0 is the BSE Hamiltonian that includes the ground-state charges only, namely built with ε = 1. Correspondingly, ΔεGW/PCMa/i corresponds to the shift of the GW energy levels when switching on the fast solvent response, that is, image file: c8sc00529j-t37.tif. Such a decomposition allows a direct correspondence with the solvatochromic shifts calculated within TD-DFT in both the LR and SS formulations.

As the easiest identification, the 3rd line contribution (eqn (17c)) straightforwardly corresponds to LR reaction field matrix elements, i.e., to transition density polarisation effects. More explicitly, 〈ai|PCMv|bj〉 describes the action on ϕb(r′)ϕj(r′) of the reaction field image file: c8sc00529j-t14.tif, where image file: c8sc00529j-t15.tif is the surface charge generated by the PCM susceptibility χPCM(ε) in response to the field generated by the transition density ϕa(r)ϕi(r). In the original notations of ref. 10, such matrix elements are strictly equivalent to the Bfai,bj = 〈ai|Kf|bj〉 linear response terms (eqn (31) in ref. 10) where Kf is the “fast” reaction potential integral operator. As such, the BSE + PCM formalism includes the LR solvent contributions.

Let us now demonstrate that the second line (eqn (17b)) recovers the SS contribution. Part of the demonstration relies on the specificity of the GW + PCM formalism that captures accurately SS dynamic polarisation energies, namely:

image file: c8sc00529j-t16.tif(18a)
image file: c8sc00529j-t17.tif(18b)
with positive (negative) shifts for occupied (virtual) levels since ΔW corresponds to a reduced interaction. This was demonstrated numerically in ref. 54 and justified in Appendix A. As such, the term of eqn (17b), that we label ΔASSai,bj, reads:
image file: c8sc00529j-t18.tif(19)

In the case of a “pure” transition between levels (i) and (a), one straightforwardly obtains:

image file: c8sc00529j-t19.tif(20)
where (aaii) represents the unrelaxed variation Δρia of the charge density between the excited and the ground states. Identifying ΔW to the reaction field PCM(ε)v as justified in Appendix B, we recognise the action on Δρia of the reaction field image file: c8sc00529j-t20.tif, where Qia (r) is the surface charge generated by the PCM susceptibility χPCM(ε) in response to the field generated by Δρia itself, with the (1/2) factor indicating that it is a self-interaction term. This corresponds to the SS expression within TD-DFT for the solvatochromic shift, namely with the notations of ref. 12:
image file: c8sc00529j-t21.tif(21)
where here i represents the ith ES. Physically, the terms in eqn (17b) represent the quantum mechanical version of the classical solvation energy associated with the interaction of a charge redistribution (dipolar, quadrupolar, etc.) with the polarisable medium, where this charge redistribution corresponds to the charge variation between the ground and the excited states.

We now turn to the case of general BSE electron–hole eigenstates image file: c8sc00529j-t22.tif assuming for simplicity the TDA approximation. Hole and electron densities are now described by a correlated codensity ρλ(re,rh) = |ψλ(re,rh)|2 that cannot be expressed in terms of individual eigenstate densities as in the simple two-level model. This however does not affect the central result that the screened Coulomb potential W, and thus the electron–hole interaction:

image file: c8sc00529j-t23.tif(22)
are properly renormalized by the reaction field. Further, in the presence of the PCM, the BSE electron–hole Hamiltonian of eqn (10) is first dressed with the reaction field operator vreac and then fully rediagonalised. As such, beyond first-order perturbation theory, the (BSE + PCM) eigenstates and codensities are relaxed with respect to the presence of the polarisable medium and the variation from gas to solvent of the electron–hole interaction results from the variation of both ρλ and W. It should be emphasized however that in the present implementation, the molecular orbitals and corresponding codensities are not relaxed with respect to the change in QM charge density upon excitation. Namely, the so-called Z-matrix relaxation effects are here neglected. From this point of view, the present BSE/PCM methodology is at that stage in analogy with the unrelaxed density (UD) approximation used in the Vertical Excitation Method (VEM) approach (VEM-UD).7,15

Concerning now the contribution from the GW electron and hole quasiparticle energies, we obtain, considering e.g. the unoccupied (electron) energy levels εGW/PCMa:

image file: c8sc00529j-t24.tif(23)
with image file: c8sc00529j-t25.tif. The variation of the electron and hole energies from the gas to the solvent thus stems from the variation of the Xia coefficients (relaxation effects) and from that of the individual εGW/PCMi/a energy levels, namely their related polarization energy ΔεGWi/a. In the correlated electron–hole BSE description, the variation of the electronic energy levels contributes to the excitation energy solvatochromic shift through a weighted average of the individual level polarization energy.

3 Computational details

The geometries of all compounds have been obtained at the MP2/6-311G(2d,2p) level in the gas-phase, but for the benzene/TCNE complex the geometry was taken from ref. 73 (see the ESI). Our GW and BSE calculations are performed with the Fiesta package44,47,53 that implements these formalisms with Gaussian basis sets and the Coulomb-fitting RI-V approach.74 Input KS eigenstates are generated with the NWChem package75 at the PBE0/cc-pVTZ level.76–78 The corresponding cc-pVTZ-RI auxiliary basis set79 is used in the many-body RI-V calculations. For improved accuracy and removal of most of the dependency on the selected DFT functional, our GW calculations are performed at the partially self-consistent evGW level68,70,72,80,81 with GW-updated eigenvalues and frozen eigenvectors.

Within the present neq approach, the starting DFT calculations are performed in combination with the COSMO formalism82 as implemented in NWChem, using the (ε0) dielectric constant associated with the (slow) nuclear and electronic degrees of freedom (e.g., ε0 = 78.39 for water). As such, the KS eigenstates used to build our electronic excitations carry information regarding the ground-state polarisation effects. In a second-step, the screened Coulomb potential used in the GW and BSE calculations, that describes fast electronic excitations out of the ground-state, is “dressed” with the reaction field generated with the (fast) electronic dielectric response in the low-frequency limit (e.g., ε = 1.78 for water). In GW and BSE, the reaction field is described at the full IEF-PCM level, following the implementation detailed in ref. 54.

Our BSE calculations are compared to TD-DFT calculations performed with the Gaussian16 package83 using the IEF-PCM non-equilibrium implementation and the same cc-pVTZ atomic basis set. The LR shifts have been obtained with the default Gaussian16 implementation,10,11 whereas the SS shifts have been determined with the so-called corrected LR (cLR) formalism,13 that is a perturbative approach. For the TD-DFT calculations we selected the same PBE0 functional, complemented in the case of CT transitions with calculations performed with the range-separated hybrid CAM-B3LYP functional,84 known to more accurately describe CT transitions in the TD-DFT context. The character of the electronic transitions in TD-DFT has been determined by estimating the effective electronic displacement induced by the excitation using the Γ index.85,86

In order to differentiate static and dynamic contributions within our approach, we make use of a ground state frozen polarization excitation energy (Ω0),15,19 obtained by considering that the polarisable medium does not respond to the fast electronic excitation, namely by setting ε = 1 while keeping the correct static dielectric constant of the solvent ε0. As such, labelling Ω the final BSE excitation energies, accounting for both static and dynamic PCM responses, the quantity (ΩΩ0) quantifies the impact of switching on the fast PCM(ε) response. We underline that our BSE calculations are performed beyond the Tamm–Dancoff approximation (TDA), that is, include the full BSE matrix. However, the inclusion of contributions from de-excitation processes complicates the simple analysis provided above concerning the LR and SS contributions with a clear distinction between 〈ab|W|ij〉 terms, namely the coupling between density terms, and 〈ai|vreac|bj〉 contributions, that is the coupling between transition dipoles. Beyond TDA, we can decompose the expectation value 〈ΨBSE|HBSE|ΨBSE〉 into

image file: c8sc00529j-t26.tif(24a)
image file: c8sc00529j-t27.tif(24b)
image file: c8sc00529j-t28.tif(24c)
where we have dropped the λ excitation index for simplicity and where the factor 2 in front of the bare Coulomb potential indicates singlet transition energies. In such a decomposition, all terms involving occupied-virtual orbital products, including the (boxed) 〈aj|W|bi〉 term, contribute to the LR response, whereas the remaining matrix elements contribute to the SS response.

4 Results

To illustrate the methodology and to confirm that the PCM-BSE formalism provides both LR and SS shifts, we perform calculations on a series of standard test molecules (Fig. 2), providing comparison with TD-DFT/PCM calculations conducted with the standard LR and SS (cLR) response formalisms. We start with transitions that have a rather local character and hence are expected to be better described within TD-DFT by the LR formalism (Subsection 4.1), before focusing in a second Subsection 4.2 on transitions having a strong CT character that requires the SS formalism. Eventually, in Subsection 4.3, we evidence that some systems may require both LR and SS contributions, highlighting the merits of the (BSE + PCM) formalism that can treat all systems on an equal footing.

4.1 Local ES: dominating LR contributions

We start our analysis with the well-known cases of acrolein and indigo (Fig. 2a and b). Acrolein is characterised by a negative solvatochromism for the lowest n–π* (A′′) transition but a positive solvatochromism of the lowest π–π* (A′) excitation. Experimental data are available in both the gas phase and water,87,88 and acrolein has been used as a benchmark system for studying continuum solvation models in conjunction with several levels of theory, e.g., TD-DFT,12,13,15 ADC(2),26 SAC-CI,21,23 and CCSD.22,25,27 Hybrid two-layer QM/MM88,89 or three-layer QM/MM/PCM90 calculations, combining TD-DFT with an atomistic force field description of the explicit solvent, have also been used to understand the pros and cons of the PCM model for water, a protic medium well-known to be challenging for continuum approaches.

Our data are compiled in Table 1. The gas phase (Ωg) and solvated (Ω) theoretical and experimental excitation energies are provided, together with the overall shift (ΩΩg). In our non-equilibrium formalism, the shift is decomposed into a ground-state static contribution (Ω0Ωg) and a dynamic one (ΩΩ0) that is itself partitioned into LR and SS contributions. Within TD-DFT, the LR and SS (cLR) shifts are obtained as two separate calculations, whereas within BSE, this decomposition is obtained by partitioning (see above) the non-equilibrium BSE/GW + PCM shift. Wavefunction approaches, such as ADC(2) or CCSD, can also be used to provide both contributions simultaneously26 or separately25 and we give some literature examples in Table 1.

Table 1 Lowest singlet excitation energies for acrolein and indigo in the gas phase (Ωg) and in water (Ω), combining TD-DFT or BSE with neq PCM (ε0 = 78.39, ε = 1.78 for water). The Ω0 energies are obtained by setting ε = 1, namely accounting only for ground-state polarisation effects with no additional surface charges induced by the excitation. The TD-PBE0 (ΩΩ0) shifts are provided using either the LR or SS (cLR) formalisms. The (ΩΩ0) BSE shifts are also partitioned into LR and SS contributions following the analysis in Section 3. All values are in eV
Ω g Ω(ΩΩg) Ω 0Ωg ΩΩ0 Ref.
a In the breakdown approach used in that work, the SS contribution is +0.17 eV and the LR contribution is negligible. b In the breakdown approach used in that work, the SS contribution is −0.21 eV and the LR contribution is −0.18 eV. c In ethanol, the most polar protic solvent in which indigo is soluble experimentally.
Acrolein n–π* in water
TD-PBE0 (LR) 3.599 3.785(+0.186) +0.189 −0.003 This work
TD-PBE0 (cLR) 3.599 3.736(+0.137) +0.189 −0.052 This work
BSE 3.736 3.988(+0.252) +0.232 −0.011 +0.031 This work
CC3 3.74 This work
CCSDR(3)/MM 3.81 4.08(+0.27) 91
SAC-CI 3.85 3.95(+0.10) 21
CCSD 3.94 4.14(+0.20) 22
ADC(2) 3.69 3.86(+0.17)a 26
CCSD (LR) 3.88 4.10(+0.22) 25
CCSD (SS) 3.88 4.05(+0.17) 25
Exp. 3.69 3.94(+0.25) 87 and 88
[thin space (1/6-em)]
Acrolein π–π* in water
TD-PBE0 (LR) 6.383 6.174(−0.209) −0.073 −0.136 This work
TD-PBE0 (cLR) 6.383 6.281(−0.102) −0.073 −0.029 This work
BSE 6.498 6.214(−0.284) −0.112 −0.163 −0.004 This work
CC3 6.82 This work
CCSDR(3)/MM 6.73 6.22(−0.51) 91
SAC-CI 6.97 6.75(−0.22) 21
CCSD 6.89 6.54(−0.35) 22
ADC(2) 6.79 6.40(−0.39)b 26
CCSD (LR) 6.80 6.39(−0.41) 25
CCSD (SS) 6.80 6.54(−0.26) 25
Exp. 6.42 5.89(−0.53)a 87 and 88
[thin space (1/6-em)]
Indigo in water
TD-PBE0 (LR) 2.304 2.160(−0.144) −0.068 −0.076 This work
TD-PBE0 (cLR) 2.304 2.229(−0.075) −0.068 −0.007 This work
BSE 2.259 2.047(−0.212) −0.122 −0.082 −0.008 This work
Exp. 2.32 2.04(−0.19)c 92, 93, 94, 95, 96, 97 and 98

The acrolein n–π* transition in the gas phase (Ωg) is found to be located at 3.60 eV and 3.74 eV within TD-PBE0 and BSE respectively, in good agreement with the CC3 value of 3.74 eV, as well as with previous wavefunction estimates and experiment. The analysis of the intermediate (Ω0Ωg) and total (ΩΩg) solvatochromic shift indicates that the positive solvatochromism for this transition is entirely dominated by ground-state effects and that the additional shift associated with the fast optical excitation is negligible. TD-PBE0 and BSE calculations performed on top of the DFT + PCM(ε0) ground-state yield similar Ω0Ωg shifts, namely +0.19 eV and +0.23 eV, respectively. Concerning the effect of switching ε (1.78 in water), both TD-PBE0+PCM, within LR or cLR, and BSE + PCM yield very small additional shifts, ranging from −0.05 eV (cLR) to +0.02 eV (BSE). The BSE computed shift of +0.25 eV turns out to be in close agreement with the experimental values (+0.25 eV), as well as with previous wavefunction approaches (see Table 1).

While the n–π* transition does not allow us to clearly discriminate between LR and SS responses, a more interesting test of the effect of the fast response (ε) polarisable environment comes with the higher-lying π–π* transition. For this transition, the reaction field associated with the optical excitation, namely vreac(ε), leads to a shift that is larger than the one associated with ground-state solvation charges. In the gas phase, the BSE transition energy (6.50 eV) is reasonably close to the CC3 (6.82 eV) and experimental (6.42 eV) values. The most salient feature is that within TD-DFT, only the LR scheme can significantly contribute to the redshift, while the cLR approach fails to deliver any sizeable solvation effect, with (ΩΩ0) being equal to −0.14 eV for the former model, and −0.03 eV for the latter. This effect was expected for a local π–π* transition not involving a strong density reorganisation between the two states: the LR-PCM-TD-DFT is more suited as it captures the dominating contributions originating from the transition densities, that can be viewed as “dispersion-like” terms.19 At the CCSD level, the results of Caricato25 also demonstrated that the LR contribution is dominant; in that work the decomposition of the total response into various contributions was not performed, but rather two different models have been applied as in TD-DFT. We also underline that in their ADC(2) study, Lunkenheimer and Köhn also found that the LR term, negligible for the n–π* case,26 becomes large for the π–π* transition, though the approach used to compute the relative contributions is not straightforwardly comparable to ours. In any case, the BSE + PCM formalism, with a −0.17 eV (ΩΩ0) shift, captures the correct bathochromic effect. Further, consistent with the TD-DFT calculations, we observe that the BSE shift is dominated by the LR contribution, while the SS term provides a negligible shift. This shows, consistent with our analysis of eqn (17c), that the BSE formalism correctly captures the LR response. When compared to experiment, the BSE + PCM shift remains too small, but we recall that we neglect here, as in any continuum approach, the explicit solvent–solute interactions that are known to be significant in the present case.26,91

To confirm the present observations, we consider the case of the lowest transition in indigo, a hallmark centro-symmetric dye presenting a low-lying dipole-allowed π–π* transition. This compound was studied previously at the TD-DFT level with the LR PCM model, and it was shown that this approach nicely reproduces the experimental solvatochromic shifts.99 With TD-DFT and the LR formalism, we found that the ground-state and dynamic polarisation effects, as measured respectively by (Ω0Ωg) and (ΩΩ0), have the same order of magnitude, whereas the cLR correction does not lead to any significant (ΩΩ0) effect, as expected for a dye in which both the ground-state and excited-state total dipoles are strictly null. As in the case of the π–π* transition in acrolein, the BSE + PCM scheme leads to a clear redshift, with in this case a good agreement with the experimental trends as well. As a matter of fact, the (ΩΩ0) BSE LR (SS) shift is very close to the corresponding TD-DFT LR (SS) shift, demonstrating the relevance of the analysis and partitioning of the BSE overall (ΩΩ0) difference discussed in Section 3.

4.2 Charge-transfer ES: dominating SS contributions

We now turn to the study of charge-transfer (CT) excitations for which the SS contribution should increase with increasing CT character. We start with the p-nitro-aniline (PNA) molecule (Fig. 2c), a typical push–pull system studied both theoretically5,100 and experimentally,101 and characterised by a low-lying CT excitation from the amine (–NH2) to the nitro (–NO2) group. Inspired by ref. 5, we first consider an artificial configuration obtained by rotating the amino group perpendicularly to the conjugated plane (PNAperp in Fig. 2d) so as to produce a transition having a pure dark CT character. We focus on the low-lying dark CT excitation showing a dominant HOMO–LUMO character. The “standard” planar PNA molecule is studied below. Further, we consider an intermolecular CT excitation in the donor–acceptor benzene/TCNE complex (Fig. 2e),38,73 for which gas phase experimental data are available.102 The lowest energy excitation is a clear HOMO–LUMO transition with a strong CT character even though a small delocalization of the HOMO (LUMO) on TCNE (Benzene) leads to a non-zero oscillator strength (see e.g.ref. 38). Since CT transitions are poorly described by global hybrid functionals with a moderate amount of exact exchange such as PBE0, TD-DFT calculations are also performed with the CAM-B3LYP range-separated hybrid functional. Our data are compiled in Table 2.
Table 2 Data as in Table 1 but for the lowest CT transitions in PNAperp and in the benzene–TCNE complex. The experimental value for the latter case is taken from ref. 102
Ω g Ω(ΩΩg) Ω 0Ωg ΩΩ0
PNA perp in water (“dark” CT excitation)
TD-PBE0 (LR) 3.686 3.316(−0.370) −0.370 0.000
TD-PBE0 (cLR) 3.686 2.861(−0.795) −0.370 −0.425
TD-CAM-B3LYP (LR) 4.621 4.308(−0.313) −0.312 −0.001
TD-CAM-B3LYP (cLR) 4.621 4.038(−0.583) −0.312 −0.271
BSE 5.112 4.399(−0.713) −0.423 +0.015 −0.304
[thin space (1/6-em)]
Benzene–TCNE in water (“bright” CT excitation)
TD-PBE0 (LR) 2.157 2.081(−0.076) −0.065 −0.011
TD-PBE0 (cLR) 2.157 1.747(−0.410) −0.065 −0.345
TD-CAM-B3LYP (LR) 2.944 2.876(−0.068) −0.061 −0.007
TD-CAM-B3LYP (cLR) 2.944 2.492(−0.452) −0.061 −0.291
BSE 3.503 3.121(−0.382) −0.086 −0.009 −0.287
Exp. 3.59

While the contribution of the fast optical dielectric response (ε) to the solvatochromic shifts, as measured by (ΩΩ0), mainly originates from the LR contribution in both acrolein and indigo, the solvent-induced dynamic shift associated with CT transitions can only be described by adopting a SS (cLR here) formalism in the TD-DFT context. This is clearly illustrated by the PNAperp system where the TD-DFT LR shift is trifling, a logical consequence of the dark nature of the considered transition, whereas the SS contribution is very large, as a result of the large density reorganisation associated with excitation. Consistent with previous studies,17 the magnitude of the TD-DFT SS shift strongly depends on the chosen functional, and it goes from −0.42 eV to −0.27 eV upon replacing the PBE0 functional by CAM-B3LYP, which is more suited for such an excited-state. The BSE (ΩΩ0) shift originates mainly from its SS component as well and lies in between the PBE0 and CAM-B3LYP values, though much closer to the latter, as expected.

The same conclusions are reached when considering the inter-molecular CT transition in the benzene–TCNE complex. First, we observe that the BSE formalism very nicely reproduces the gas phase excitation energy that is available experimentally. The TD-CAM-B3LYP calculation also provide a reasonable value, while the TD-PBE0 approach yields a much too small Ωg, a logical consequence of its lack of long-range corrections. As expected for TD-DFT, the LR approach again provides a negligible (ΩΩ0) dynamic shift, while the SS formalism yields a large redshift. Again, the BSE + PCM formalism captures the correct physics, with a negligible LR contribution and a large SS contribution. Such calculation clearly demonstrates that the proposed BSE + PCM approach can also, following eqn (17b), capture SS dynamic shifts.

4.3 Hybrid ES

We now turn to difficult cases where neither the LR nor SS contributions can be neglected, and our results are collected in Table 3. As such, TD-DFT calculations adopting one or the other response formalisms cannot capture the entire solvatochromic shifts, though, as discussed in ref. 26 both contributions should in principle be accounted for. This is first evidenced by considering the planar (standard) PNA molecule. The much smaller disagreement between TD-PBE0 and TD-CAM-B3LYP, as compared to the PNAperp case, is indeed a first indication that the CT character is here significantly reduced. As a result (see Table 3), both SS and LR contributions to the (ΩΩ0) shift are large: they are of equal magnitude within TD-PBE0 while the SS contribution is larger with TD-CAM-B3LYP. Clearly, the BSE formalism simultaneously accounts for both contributions, the SS contribution being slightly larger than the LR one, consistent with the TD-DFT results. The total (ΩΩ0) dynamic shift amounts to −0.18 eV and −0.27 eV when adding the LR and cLR contribution determined at TD-PBE0 and TD-CAM-B3LYP, respectively, comparable to the −0.19 eV solvatochromic shift obtained at the BSE level.
Table 3 Data as in Table 1 but for the transitions of mixed character
Ω g Ω(ΩΩg) Ω 0Ωg ΩΩ0
PNA in water (partial CT excitation)
TD-PBE0 (LR) 4.202 3.802(−0.400) −0.314 −0.086
TD-PBE0 (cLR) 4.202 3.796(−0.406) −0.314 −0.092
TD-CAM-B3LYP (LR) 4.513 4.105(−0.408) −0.321 −0.087
TD-CAM-B3LYP (cLR) 4.513 4.005(−0.508) −0.321 −0.187
BSE 4.527 3.864(−0.663) −0.470 −0.090 −0.102
[thin space (1/6-em)]
4-Nitropyridine N-oxide in benzene (mixed excitation)
TD-PBE0 (LR) 3.989 3.815(−0.174) −0.048 −0.126
TD-PBE0 (cLR) 3.989 3.797(−0.192) −0.048 −0.144
TD-CAM-B3LYP (LR) 4.196 4.023(−0.173) −0.033 −0.140
TD-CAM-B3LYP (cLR) 4.196 4.109(−0.087) −0.033 −0.054
BSE 3.966 3.658(−0.267) −0.001 −0.193 −0.073
[thin space (1/6-em)]
4-Nitropyridine N-oxide in water (mixed excitation)
TD-PBE0 (LR) 3.989 3.797(−0.192) −0.097 −0.095
TD-PBE0 (cLR) 3.989 3.827(−0.162) −0.097 −0.065
TD-CAM-B3LYP (LR) 4.196 4.028(−0.168) −0.065 −0.103
TD-CAM-B3LYP (cLR) 4.196 4.067(−0.129) −0.065 −0.064
BSE 3.966 3.687(−0.279) −0.070 −0.141 −0.068

4-Nitropyridine N-oxide is an organic probe used to assess the nature of solvents following a Kamel–Taft type of analysis.103 A previous throughout theoretical analysis of the solvatochromism of this probe is available,104 and demonstrates that, none of the available LR or SS PCM model is able to describe the solvent effects in a TD-DFT context, as the observed spectral changes come from a fine interplay of several effects. For this compound, the obtained conclusions are similar to the PNA case: both the LR and SS contributions are non-negligible even though the SS dynamic shift is somehow larger. Adding the LR and SS TD-DFT contributions, the overall (ΩΩ0) shift amounts to −0.16 eV and −0.17 eV with PBE0 and CAM-B3LYP in water, respectively, a shift that can be compared to the −0.21 eV value obtained within BSE. The experimental gas-phase excitation energy obtained through extrapolation in ref. 104 is 3.80 eV, and one notices that the BSE value is reasonably close to that estimate. In benzene, the measurement gives 3.52 eV,103 corresponding to a solvatochromic shift of −0.28 eV, a value that BSE can reproduce (−0.27 eV), whereas TD-DFT in unable to do so. In water, the hydrogen bonds with the negatively charged oxygen atom of the probe play a crucial role, and a blueshift is observed (the excitation maximum takes place at 3.94 eV),103 an effect that all continuum approaches logically fail to capture.

5 Further discussion

We summarise in Fig. 3 the principal findings of the present study, namely that the BSE/GW formalism accounts for both LR and SS shifts, a fact demonstrated for a large variety of transitions, including excitations where one of the two terms is dominant and excitations for which both are needed. While the partitioning between BSE LR and SS terms is provided in eqn (17b)–(17c), it is interesting to analyse now the contributions to the BSE dynamic shift originating from the solvent-induced renormalisation of the ionisation potential and electronic affinity, namely the variation of the GW HOMO–LUMO gap in the solvent, and from the reduction of the electron–hole interaction contained in the 〈ab|W|ij〉 terms.
image file: c8sc00529j-f3.tif
Fig. 3 Schematic representation of the (ΩΩ0) dynamic solvent induced shifts for the π–π* excitation in acrolein, the CT excitation in the benzene–TCNE complex, and the (planar) PNA mixed excitation. By dynamic shift, we mean the effect of switching the fast vreac(ε) PCM reaction field on top of the (slow) ground-state vreac(ε0) PCM response. The represented data correspond to the (ΩΩ0) shifts given in Tables 1–3.

For the π–π* transition in acrolein, with a shift captured by the LR scheme only at the TD-DFT level, the decomposition of the (ΩΩ0) energy difference confirms that the variation of the BSE SS-like 〈ΔεW〉 component (eqn (17b)) does not contribute to the shift. Such a result may seem surprising since, as expected, the 〈ΔεGW〉 HOMO–LUMO gap becomes smaller by 2.34 eV upon “switching” vreac(ε), consistent with an (absolute) polarisation energy of ca. 1.1–1.2 eV that affects the IP and AE of hydrated compounds compared to the gas phase.54 However, very remarkably, the 〈W〉 electron–hole binding energy is also decreased by 2.35 eV. Namely, the screening by the fast reaction field reduces similarly the occupied-to-virtual (εaGWεiGW) GW gap and the electron–hole binding energy (see Fig. 1). As such, only the 〈vreac(ε)〉 LR terms (eqn (17c)) explain the solvatochromic effect, consistent with the TD-DFT results.

We now turn to the opposite situation of a CT excitation for which the solvatochromic shift can only be described by SS implementations, such as cLR, within TD-DFT. In the case of the benzene–TCNE complex, the LR 〈ia|vreac(ε)|bj〉 contributing matrix elements become very small, as a logical consequence of the small overlap between the donor and acceptor wavefunctions. As such the only contribution to the (ΩΩ0) shift originates from the variation of the SS-like 〈ΔεW〉 term. The impact of vreac(ε) leads to the reduction of the 〈Δε〉 average gap by 1.60 eV, while the electron–hole binding energy 〈W〉 reduces by a smaller 1.31 eV amount, accounting for most of the (ΩΩ0) = −0.29 eV redshift reported in Table 2.

An important consequence of the present analysis, showing that, in BSE/GW theory, SS shifts result from the competition between the reduction of the GW occupied-to-virtual (εaGWεiGW) energy gaps and the electron–hole 〈ab|W|ij〉 binding energies, is that both electron–electron and electron–hole interactions must be treated on the same footing, namely here through the screened Coulomb potential W. In particular, BSE + PCM calculations starting from KS eigenstates generated with exchange-correlation functionals optimally-tuned, so as to generate the correct gas phase HOMO–LUMO gap, might not deliver the correct SS contribution to optical excitation energy shifts.

Finally, while Table 2 indicates a large variability of the SS (cLR) TD-DFT (ΩΩ0) dynamic shift as a function of the chosen XC functional, it is interesting to emphasize the stability of the BSE/evGW data with respect to the input KS eigenstates used to build the evGW electronic energy levels and the screened Coulomb potential W. This is illustrated in Fig. 4a where we plot the PNAperp excitation energy in water with the cLR-TD-PBE(α) (open triangles) and BSE/evGW@PBE(α) (open circles) methods. Here α indicates the percentage of exact exchange (EEX) used in the hybrid functional and going from 0% (PBE) to 100% (HFPBE) in Fig. 4. While the TD-DFT excitation energy is shown to increase steeply as a function of α, the BSE data are much more stable, with a variation of ca. 0.12 eV from (α = 0%) to (α = 90%). The stability of gas phase BSE/evGW excitation energies with respect to the starting KS starting point has been documented in previous benchmark studies,47,49 and is therefore shown here to pertain for condensed-phase calculations. In fact, plotting now in Fig. 4b the dynamic (ΩΩ0) solvatochromic shifts, we observe again a large variability of the TD-PBE(α) cLR shifts (the LR contribution is vanishingly small), while again the BSE shifts are extremely stable, evolving from −0.29 eV to −0.25 eV with α. To rationalize this result, we recall that BSE calculations are performed on top of partially self-consistent evGW calculations where the corrected electronic energy levels are self-consistently reinjected during the construction of G and W. This self-consistent treatment of the quasiparticle energies {εnGW} and screened Coulomb potential W leads to a large stability of the BSE Hamiltonian and reaction field, with a small residual dependency on the starting KS eigenstates due to the fact that the {ϕn} eigenstates are unchanged. Such a stability of the BSE/evGW excitation energies and solvent induced shifts, that allows us to alleviate the standard problem of the proper choice of the XC functional central in TD-DFT, is one of the interesting features of the present scheme. As a fair tribute to TD-DFT, we observe however that the cLR shift matches closely that of BSE for a large exact exchange ratio. It is indeed known that large amounts of exact exchange are required in TD-DFT to obtain an accurate description of the ground-to-excited density variation for CT excitations.17

image file: c8sc00529j-f4.tif
Fig. 4 (a) Evolution of the PNAperp excitation energy Ω in solution as a function of the percentage of exact exchange α in the PBE(α) functional. Triangles correspond to the cLR-TD-DFT results while circles indicate BSE/evGW@PBE(α) data. (b) Dynamic shift (ΩΩ0) obtained with cLR-TD-DFT (triangles) and BSE/evGW@PBE(α) (circles). The inconsistent value appearing with TD-DFT for α = 40% is due to the fact that two excited states of mixed character are nearly degenerate for that EEX ratio. Energies are in eV.

6 Conclusions

We have demonstrated that the Bethe–Salpeter equation (BSE) formalism combined with the PCM approach straightforwardly includes, on the same footing, both linear-response (LR) and state-specific (SS) contributions to the solvatochromic shifts. This result was established on analytic grounds as well as by comparing the numerical results obtained for several molecules or complexes using BSE + PCM or TD-DFT + PCM approaches, the latter being performed with both the LR and cLR formalisms. To the best of our knowledge, such a consistent theoretical approach simultaneously providing both LR and SS solvent shifts remains unavailable in the TD-DFT world, and has been developed previously, only for theories presenting a less favourable scaling with system size, e.g., ADC(2).26

Beyond the standard LR 〈ab|vreac(ε)|ij〉 matrix element contributions, where vreac(ε) is the PCM reaction field to the electronic excitation, the proper inclusion of SS shifts hinges on the incorporation of vreac(ε) in the screened Coulomb potential W. This dressed screened Coulomb potential renormalises both electron–electron and electron–hole interactions, namely both the GW quasiparticle energies for occupied/virtual electronic levels and the BSE electron–hole (excitonic) binding energy. Ground-state polarisation effects are accounted for, in the present non-equilibrium scheme, by starting our BSE/GW calculations with KS eigenstates generated at the DFT + PCM(ε0) level.

Following previous studies related to merging the GW formalism with discrete polarisable models,56–58 the same analysis can be straightforwardly applied to BSE calculations combined with other polarisation models.55,57,105 Specifically, the use of an explicit (molecular) description of the environment, combined with, e.g., empirical force fields, should allow us to account for ground-state electrostatic field effects induced by the solvent molecule static multipoles, an important contribution in the case of polar solvents that cannot be captured by the PCM model. As shown in the case of the π–π* transition in aqueous acrolein,89,91 such contributions can explain the difference between the experimental shift (−0.53 eV) and the ∼−0.25 eV redshifts obtained with the present BSE or TD-DFT-LR formalisms combined with the PCM.

A From ΔεGW to the Born solvation energy model

In this Appendix, we justify eqn (18a). As shown in ref. 54, the variation of the εnGW quasi-particle energies can be simply expressed within the so-called static Coulomb-hole (COH) plus screened-exchange (SEX) approximation to the self-energy ΣGW, namely
with the following matrix elements:
image file: c8sc00529j-t29.tif(26)
image file: c8sc00529j-t30.tif(27)
where W is taken in the low-frequency optical limit and [W with combining tilde] = Wv. We therefore get for the variation ΔεnGW:
image file: c8sc00529j-t31.tif(28)
where {ϕn} represents the frozen ground-state eigenstates generated at the DFT + PCM(ε0) level. Here, one can expect the principal contribution to eqn (28) to come from the reaction field induced by the monopolar terms associated with the unit charge density ϕ*n(r)ϕn(r) = |ϕn(r)|2, i.e.,
image file: c8sc00529j-t32.tif(29)
with the + or − sign depending on whether n describes an empty (a) or an occupied (i) state. All other contributions are expected to be smaller, implying most dipole–dipole interactions between neutral codensities since, by orthogonalization, image file: c8sc00529j-t33.tif and image file: c8sc00529j-t34.tif. A similar argument was provided and validated by Neaton and co-workers in their study of the screening of molecular energy levels by metallic electrodes,106 in which they assume that ΔW is slowly varying over the solute cavity, that is, the monopolar term is dominant.

B Relation between ΔW and vreac

Working out eqn (14) and (16) to lower order in vreac = PCMv leads to:
image file: c8sc00529j-t35.tif

image file: c8sc00529j-t36.tif
with χ0, ε = 1 − 0, and W being the solute free-electron susceptibility, dielectric matrix and screened Coulomb potential in the absence of the PCM. Neglecting the internal response of the QM solute in front of the PCM “reservoir”, namely taking χ0 → 0, one obtains that indeed ΔW is equal to the reaction field.

Conflicts of interest

There are no conflicts to declare.


This project has received funding from the European Union Horizon 2020 research and innovation programme under grant agreement no. 646176 (EXTMOS). Calculations were performed using resources from the GENCI-CINES/IDRIS French national supercomputing centers, the CCIPL (Centre de Calcul Intensif des Pays de la Loire), and Grenoble CIMENT (Calcul Intensif, Modélisation, Expérimentation Numérique et Technologique) computational mesocenters. C. A. G and D. J. are indebted to the LumoMat consortium for support through the FCPol-Resp project.


  1. M. E. Casida, J. Mol. Struct.: THEOCHEM, 2009, 914, 3–18 CrossRef CAS.
  2. C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications, Oxford University Press, New York, 2012 Search PubMed.
  3. C. van Caillie and R. D. Amos, Chem. Phys. Lett., 1999, 308, 249–255 CrossRef CAS.
  4. F. Furche and R. Ahlrichs, J. Chem. Phys., 2002, 117, 7433–7447 CrossRef CAS.
  5. G. Scalmani, M. J. Frisch, B. Mennucci, J. Tomasi, R. Cammi and V. Barone, J. Chem. Phys., 2006, 124, 094107 CrossRef PubMed.
  6. J. Liu and W. Z. Liang, J. Chem. Phys., 2011, 135, 184111 CrossRef PubMed.
  7. C. A. Guido, G. Scalmani, B. Mennucci and D. Jacquemin, J. Chem. Phys., 2017, 146, 204106 CrossRef PubMed.
  8. S. Miertǔs, E. Scrocco and J. Tomasi, Chem. Phys., 1981, 55, 117–129 CrossRef.
  9. E. Cancès, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032–3041 CrossRef.
  10. R. Cammi and B. Mennucci, J. Chem. Phys., 1999, 110, 9877–9886 CrossRef CAS.
  11. M. Cossi and V. Barone, J. Chem. Phys., 2001, 115, 4708–4717 CrossRef CAS.
  12. R. Cammi, S. Corni, B. Mennucci and J. Tomasi, J. Chem. Phys., 2005, 122, 104513 CrossRef CAS PubMed.
  13. M. Caricato, B. Mennucci, J. Tomasi, F. Ingrosso, R. Cammi, S. Corni and G. Scalmani, J. Chem. Phys., 2006, 124, 124520 CrossRef PubMed.
  14. R. Improta, V. Barone, G. Scalmani and M. J. Frisch, J. Chem. Phys., 2006, 125, 054103 CrossRef PubMed.
  15. A. V. Marenich, C. J. Cramer, D. G. Truhlar, C. A. Guido, B. Mennucci, G. Scalmani and M. J. Frisch, Chem. Sci., 2011, 2, 2143–2161 RSC.
  16. A. Pedone, J. Chem. Theory Comput., 2013, 9, 4087–4096 CrossRef CAS PubMed.
  17. C. A. Guido, D. Jacquemin, C. Adamo and B. Mennucci, J. Chem. Theory Comput., 2015, 11, 5782–5790 CrossRef CAS PubMed.
  18. C. A. Guido, B. Mennucci, G. Scalmani and D. Jacquemin, J. Chem. Theory Comput., 2018, 14, 1544–1553 CrossRef CAS PubMed.
  19. S. Corni, R. Cammi, B. Mennucci and J. Tomasi, J. Chem. Phys., 2005, 123, 134512 CrossRef CAS PubMed.
  20. R. Guareschi, O. Valsson, C. Curutchet, B. Mennucci and C. Filippi, J. Phys. Chem. Lett., 2016, 7, 4547–4553 CrossRef CAS PubMed.
  21. R. Cammi, R. Fukuda, M. Ehara and H. Nakatsuji, J. Chem. Phys., 2010, 133, 024104 CrossRef PubMed.
  22. M. Caricato, B. Mennucci, G. Scalmani, G. W. Trucks and M. J. Frisch, J. Chem. Phys., 2010, 132, 084102 CrossRef PubMed.
  23. R. Fukuda, M. Ehara, H. Nakatsuji and R. Cammi, J. Chem. Phys., 2011, 134, 104109 CrossRef PubMed.
  24. M. Caricato, J. Chem. Theory Comput., 2012, 8, 4494–4502 CrossRef CAS PubMed.
  25. M. Caricato, J. Chem. Phys., 2013, 139, 044116 CrossRef PubMed.
  26. B. Lunkenheimer and A. Köhn, J. Chem. Theory Comput., 2013, 9, 977–994 CrossRef CAS PubMed.
  27. M. Caricato, Comput. Theor. Chem., 2014, 1040, 99–105 CrossRef.
  28. R. Fukuda and M. Ehara, J. Chem. Phys., 2014, 141, 154104 CrossRef PubMed.
  29. J.-M. Mewes, Z.-Q. You, M. Wormit, T. Kriesche, J. M. Herbert and A. Dreuw, J. Phys. Chem. A, 2015, 119, 5446–5464 CrossRef CAS PubMed.
  30. J.-M. Mewes, J. M. Herbert and A. Dreuw, Phys. Chem. Chem. Phys., 2017, 19, 1644–1654 RSC.
  31. E. E. Salpeter and H. A. Bethe, Phys. Rev., 1951, 84, 1232–1242 CrossRef.
  32. L. J. Sham and T. M. Rice, Phys. Rev., 1966, 144, 708–714 CrossRef CAS.
  33. W. Hanke and L. J. Sham, Phys. Rev. Lett., 1979, 43, 387–390 CrossRef CAS.
  34. G. Strinati, Phys. Rev. Lett., 1982, 49, 1519–1522 CrossRef CAS.
  35. X. Blase, I. Duchemin and D. Jacquemin, Chem. Soc. Rev., 2018, 47, 1022–1043 RSC.
  36. D. Rocca, D. Lu and G. Galli, J. Chem. Phys., 2010, 133, 164109 CrossRef PubMed.
  37. J. M. Garcia-Lastra and K. S. Thygesen, Phys. Rev. Lett., 2011, 106, 187402 CrossRef CAS PubMed.
  38. X. Blase and C. Attaccalite, Appl. Phys. Lett., 2011, 99, 171909 CrossRef.
  39. I. Duchemin, T. Deutsch and X. Blase, Phys. Rev. Lett., 2012, 109, 167801 CrossRef CAS PubMed.
  40. B. Baumeier, D. Andrienko and M. Rohlfing, J. Chem. Theory Comput., 2012, 8, 2790–2795 CrossRef CAS PubMed.
  41. S. Sharifzadeh, P. Darancet, L. Kronik and J. B. Neaton, J. Phys. Chem. Lett., 2013, 4, 2197–2201 CrossRef CAS.
  42. P. Cudazzo, M. Gatti, A. Rubio and F. Sottile, Phys. Rev. B, 2013, 88, 195152 CrossRef.
  43. V. Ziaei and T. Bredow, J. Chem. Phys., 2016, 145, 174305 CrossRef PubMed.
  44. D. Escudero, I. Duchemin, X. Blase and D. Jacquemin, J. Phys. Chem. Lett., 2017, 8, 936–940 CrossRef CAS PubMed.
  45. P. Boulanger, D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2014, 10, 1212–1218 CrossRef CAS PubMed.
  46. C. Azarias, I. Duchemin, X. Blase and D. Jacquemin, J. Chem. Phys., 2017, 146, 034301 CrossRef PubMed.
  47. D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2015, 11, 3290–3304 CrossRef CAS PubMed.
  48. F. Bruneval, S. M. Hamed and J. B. Neaton, J. Chem. Phys., 2015, 142, 244101 CrossRef PubMed.
  49. D. Jacquemin, I. Duchemin and X. Blase, J. Chem. Theory Comput., 2015, 11, 5340–5359 CrossRef CAS PubMed.
  50. D. Jacquemin, I. Duchemin, A. Blondel and X. Blase, J. Chem. Theory Comput., 2016, 12, 3969–3981 CrossRef CAS PubMed.
  51. D. Jacquemin, I. Duchemin, A. Blondel and X. Blase, J. Chem. Theory Comput., 2017, 13, 767–783 CrossRef CAS PubMed.
  52. T. Rangel, S. M. Hamed, F. Bruneval and J. B. Neaton, J. Chem. Phys., 2017, 146, 194108 CrossRef PubMed.
  53. D. Jacquemin, I. Duchemin and X. Blase, J. Phys. Chem. Lett., 2017, 8, 1524–1529 CrossRef CAS PubMed.
  54. I. Duchemin, D. Jacquemin and X. Blase, J. Chem. Phys., 2016, 144, 164106 CrossRef PubMed.
  55. B. Baumeier, M. Rohlfing and D. Andrienko, J. Chem. Theory Comput., 2014, 10, 3104–3110 CrossRef CAS PubMed.
  56. J. Li, G. D'Avino, I. Duchemin, D. Beljonne and X. Blase, J. Phys. Chem. Lett., 2016, 7, 2814–2820 CrossRef CAS PubMed.
  57. J. Li, G. D'Avino, A. Pershin, D. Jacquemin, I. Duchemin, D. Beljonne and X. Blase, Physical Review Materials, 2017, 1, 025602 CrossRef.
  58. J. Li, G. D'Avino, I. Duchemin, D. Beljonne and X. Blase, Phys. Rev. B, 2018, 97, 035108 CrossRef.
  59. M. Rohlfing and S. G. Louie, Phys. Rev. Lett., 1998, 80, 3320–3323 CrossRef CAS.
  60. L. X. Benedict, E. L. Shirley and R. B. Bohn, Phys. Rev. Lett., 1998, 80, 4514–4517 CrossRef CAS.
  61. S. Albrecht, L. Reining, R. Del Sole and G. Onida, Phys. Rev. Lett., 1998, 80, 4510–4513 CrossRef CAS.
  62. G. Onida, L. Reining and A. Rubio, Rev. Mod. Phys., 2002, 74, 601–659 CrossRef CAS.
  63. G. Bussi, Phys. Scr., 2004, 2004, 141 CrossRef.
  64. L. Hedin, Phys. Rev., 1965, 139, A796–A823 CrossRef.
  65. M. S. Hybertsen and S. G. Louie, Phys. Rev. B, 1986, 34, 5390–5413 CrossRef CAS.
  66. R. W. Godby, M. Schlüter and L. J. Sham, Phys. Rev. B, 1988, 37, 10159–10175 CrossRef.
  67. B. Farid, R. Daling, D. Lenstra and W. van Haeringen, Phys. Rev. B, 1988, 38, 7530–7534 CrossRef.
  68. C. Faber, C. Attaccalite, V. Olevano, E. Runge and X. Blase, Phys. Rev. B, 2011, 83, 115123 CrossRef.
  69. K. Krause, M. E. Harding and W. Klopper, Mol. Phys., 2015, 113, 1952–1960 CrossRef CAS.
  70. F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers and M. J. van Setten, J. Chem. Theory Comput., 2016, 12, 2528–2541 CrossRef CAS PubMed.
  71. J. W. Knight, X. Wang, L. Gallandi, O. Dolgounitcheva, X. Ren, J. V. Ortiz, P. Rinke, T. Körzdörfer and N. Marom, J. Chem. Theory Comput., 2016, 12, 615–626 CrossRef CAS PubMed.
  72. T. Rangel, S. M. Hamed, F. Bruneval and J. B. Neaton, J. Chem. Theory Comput., 2016, 12, 2834–2842 CrossRef CAS PubMed.
  73. T. Stein, L. Kronik and R. Baer, J. Am. Chem. Soc., 2009, 131, 2818–2820 CrossRef CAS PubMed.
  74. I. Duchemin, J. Li and X. Blase, J. Chem. Theory Comput., 2017, 13, 1199–1208 CrossRef CAS PubMed.
  75. M. Valiev, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, H. V. Dam, D. Wang, J. Nieplocha, E. Apra, T. Windus and W. de Jong, Comput. Phys. Commun., 2010, 181, 1477–1489 CrossRef CAS.
  76. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023 CrossRef.
  77. J. P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., 1996, 105, 9982–9985 CrossRef CAS.
  78. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170 CrossRef CAS.
  79. F. Weigend, J. Chem. Phys., 2002, 116, 3175–3183 CrossRef CAS.
  80. X. Blase, C. Attaccalite and V. Olevano, Phys. Rev. B, 2011, 83, 115103 CrossRef.
  81. C. Faber, J. L. Janssen, M. Côté, E. Runge and X. Blase, Phys. Rev. B, 2011, 84, 155104 CrossRef.
  82. A. Klamt, C. Moya and J. Palomar, J. Chem. Theory Comput., 2015, 11, 4220–4225 CrossRef CAS PubMed.
  83. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian-16 Revision A.03, Gaussian Inc., Wallingford CT, 2016 Search PubMed.
  84. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57 CrossRef CAS.
  85. C. A. Guido, P. Cortona and C. Adamo, J. Chem. Phys., 2014, 140, 104101 CrossRef PubMed.
  86. C. A. Guido, P. Cortona, B. Mennucci and C. Adamo, J. Chem. Theory Comput., 2013, 9, 3118–3126 CrossRef CAS PubMed.
  87. A. Moskvin, O. Yablonskii and L. Bondar, Theor. Exp. Chem., 1966, 2, 469 CrossRef.
  88. K. Aidas, A. Møgelhøj, E. J. K. Nilsson, M. S. Johnson, K. V. Mikkelsen, O. Christiansen, P. Söderhjelm and J. Kongsted, J. Chem. Phys., 2008, 128, 194503 CrossRef PubMed.
  89. J. M. Olsen, K. Aidas and J. Kongsted, J. Chem. Theory Comput., 2010, 6, 3721–3734 CrossRef CAS.
  90. A. H. Steindal, K. Ruud, L. Frediani, K. Aidas and J. Kongsted, J. Phys. Chem. B, 2011, 115, 3027–3037 CrossRef CAS PubMed.
  91. K. Sneskov, E. Matito, J. Kongsted and O. Christiansen, J. Chem. Theory Comput., 2010, 6, 839–850 CrossRef CAS PubMed.
  92. S. E. Sheppard and P. T. Newsome, J. Am. Chem. Soc., 1942, 64, 2937–3946 CrossRef CAS.
  93. P. W. Saddler, J. Org. Chem., 1956, 21, 316–318 CrossRef.
  94. M. Klessinger and W. Lüttke, Chem. Ber., 1966, 99, 2136–2145 CrossRef CAS.
  95. E. Wille and W. Lüttke, Angew. Chem., Int. Ed., 1971, 10, 803–804 CrossRef CAS.
  96. A. R. Monahan and J. E. Kuder, J. Org. Chem., 1972, 37, 4182–4184 CrossRef CAS.
  97. G. Haucke and G. Graness, Angew. Chem., Int. Ed., 1995, 34, 67–68 CrossRef CAS.
  98. R. Gerken, L. Fitjer, P. Müller, I. Usón, K. Kowski and P. Rademacher, Tetrahedron, 1999, 55, 14429–14434 CrossRef.
  99. D. Jacquemin, J. Preat, V. Wathelet and E. A. Perpète, J. Chem. Phys., 2006, 124, 074104 CrossRef PubMed.
  100. R. Cammi, L. Frediani, B. Mennucci and K. Ruud, J. Chem. Phys., 2003, 119, 5818–5827 CrossRef CAS.
  101. S. Millefiori, G. Favini, A. Millefiori and D. Grasso, Spectrochim. Acta, Part A, 1977, 33, 21–27 CrossRef.
  102. I. Hanazaki, J. Phys. Chem., 1972, 76, 1982–1989 CrossRef CAS.
  103. A. F. Lagalante, R. J. Jacobson and T. J. Bruno, J. Org. Chem., 1996, 61, 6404–6406 CrossRef CAS PubMed.
  104. S. Budzak, A. D. Laurent, C. Laurence, M. Miroslav and D. Jacquemin, J. Chem. Theory Comput., 2016, 12, 1919–1929 CrossRef CAS PubMed.
  105. D. Varsano, S. Caprasecca and E. Coccia, J. Phys.: Condens. Matter, 2017, 29, 013002 CrossRef PubMed.
  106. J. B. Neaton, M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett., 2006, 97, 216405 CrossRef CAS PubMed.


Electronic supplementary information (ESI) available: Cartesian coordinates of the compounds. See DOI: 10.1039/c8sc00529j

This journal is © The Royal Society of Chemistry 2018