Assessment of the Ab Initio Bethe-Salpeter Equation Approach for the Low-Lying Excitation Energies of Bacteriochlorophylls and Chlorophylls.

Bacteriochlorophyll and chlorophyll molecules are crucial building blocks of the photosynthetic apparatus in bacteria, algae, and plants. Embedded in transmembrane protein complexes, they are responsible for the primary processes of photosynthesis: excitation energy and charge transfer. Here, we use ab initio many-body perturbation theory within the GW approximation and Bethe-Salpeter equation (BSE) approach to calculate the electronic structure and optical excitations of bacteriochlorophylls a, b, c, d, and e and chlorophylls a and b. We systematically study the effects of the structure, basis set size, partial self-consistency in GW, and the underlying exchange-correlation approximation and compare our calculations with results from time-dependent density functional theory, multireference RASPT2, and experimental literature results. We find that optical excitations calculated with GW+BSE are in excellent agreement with experimental data, with an average deviation of less than 100 meV for the first three bright excitations of the entire family of (bacterio)chlorophylls. Contrary to state-of-the-art time-dependent density functional theory (TDDFT) with an optimally tuned range-separated hybrid functional, this accuracy is achieved in a parameter-free approach. Moreover, GW+BSE predicts the energy differences between the low-energy excitations correctly and eliminates spurious charge transfer states that TDDFT with (semi)local approximations is known to produce. Our study provides accurate reference results and highlights the potential of the GW+BSE approach for the simulation of larger pigment complexes.

Introduction

Electronic excitations form the foundation of some of the most fundamental natural processes. In photosynthesis, plants, algae, and bacteria convert solar energy into chemical energy, utilizing a cascade of coupled energy and charge transfer excitations that are performed by pigment–protein complexes with high quantum efficiency. Bacteriochlorophyll (BCL) and chlorophyll (CL) molecules are among the most important building blocks of these pigment–protein complexes.[1] They are responsible for the absorption and transfer of excitation energy and for the charge separation necessary for establishing a proton gradient that eventually drives the synthesis of chemical energy in plants and bacteria.[2] Accurately calculating the electronic structure and excitations of these molecules from first principles is the prerequisite for understanding their interactions with each other and with the surrounding proteins and, consequently, energy and charge transfer in natural photosynthesis.

BCL and CL molecules constitute a family of substituted tetrapyrroles with varying absorption properties depending on conjugation and the number and nature of substitutions. CL a and b are present in plants and green algae, whereas green bacteria mostly rely on BCL c, d, and e for excitation energy transfer and BCL a for concentrating excitations close to the reaction center of the photosynthetic unit.[3] BCL a is also the main pigment in purple bacteria, whose light harvesting apparatus and reaction center are among the most thoroughly studied natural light-harvesting systems.[4] The optical excitation spectrum of these pigments possesses two characteristic absorption bands: (1) the Q band in the visible part of the spectrum, composed of excitations Q and Q with high- and low-oscillator strength, respectively and (2) the B (or Soret) band in the near ultraviolet.

In the field of finite organic and biological molecular systems, neutral excitations and optical spectra are predominantly calculated using time-dependent density functional theory (TDDFT). In conjunction with model Hamiltonian approaches, TDDFT has been employed for the simulation of large photosynthetic pigment–protein complexes.[5,6] The accuracy of its approximations and implementations has been tested for a variety of biochromophores.[7,8,47] However, TDDFT’s standard approximations are inadequate for describing long-range charge transfer excitations[9] and high-energy Rydberg states[10] due to self-interaction errors and an incorrect asymptotic behavior. Exchange–correlation (xc) functionals that contain long-range exact exchange, such as optimally tuned range-separated hybrid functionals (OT-RSH), can be employed as a remedy in such cases,[11,12] but the use of such functionals requires a tedious per-system tuning procedure.

Multireference wavefunction-based methods have scarcely been used for molecules as large as BCL and CL. Vertical excitation energies of CL a, based on ADC(2) and different coupled cluster approaches show a spread of ∼0.4 eV, strongly depending on the method, basis set, and structural model used in these calculations.[13 −15] In 2016 and 2019, Anda et al. reported multistate RASPT2/RASSCF excitation energies of several BCL units within the light-harvesting system 2 of a purple bacterium.[16,17] A RASPT2 approach was also combined with electrostatic embedding of fixed-point charges to simulate the effect of the protein environment on excitation energies of the same system by Segatta et al.[18] While these reports constitute important advances in the use of wavefunction-based methods for complex biological molecules, they were performed with relatively small basis sets and show a dependence on the choice of the restricted active space (RAS).

The ab initio Bethe–Salpeter equation (BSE) approach, when rigorously based on many-body Green’s function theory, is an alternative method for describing neutral excitations of correlated many-electron systems.[19] It is based on a framework of charged excitation energies that correspond to electron addition and removal energies and that are most frequently calculated within the GW approximation. The GW+BSE approach has been shown to be successful in predicting the optical spectra of bulk solids[20,21] and low-dimensional materials.[22] In recent years, it has also begun to be applied to finite systems, such as small molecules,[23−25] and larger molecular complexes,[26−28] for which its accuracy has been shown to be comparable to single-reference wavefunction methods for both localized and charge transfer excitations,[29] at a substantially reduced computational cost.

In this article, we assess the accuracy of the ab initioGW+BSE approach for the Q, Q, and the first bright B excitation of several members of the BCL and CL family and the chemically closely related bacteriochlorin (BC) molecule. We compare two different approaches for approximating the electronic self-energy Σ = iGW: (1) G0W0, a one-shot method, in which the zeroth-order Green’s function G0 and screened Coulomb interaction W0 are constructed from a DFT eigensystem and directly used to correct DFT eigenvalues perturbatively and (2) partially self-consistent GW (evGW), in which the corrected eigenvalues are used to iteratively recalculate G and/or W until self-consistency is reached. We compare our results to TDDFT calculations with the local density approximation (LDA), two global hybrids and an OT-RSH functional, with RASPT2 literature results,[16,17] and with experimental data.[30,31]

We find that the GW+BSE approach used in a partially self-consistent manner results in excitation energies in the visible and near-ultraviolet within less than 100 meV from experiment for the entire family of pigments studied here. Our results are almost completely independent of the DFT eigensystem used as an input for the GW+BSE calculations. In fact, even a simple and computationally inexpensive LDA starting point leads to excellent agreement with experiment and eliminates spurious charge transfer excitations between Q and Q that TDDFT with (semi)local functionals produces. Contrary to TDDFT, GW+BSE also correctly predicts the energy difference between the two Q-band excitations, a crucial prerequisite for understanding the coupling of excitations in systems consisting of more than one pigment. Finally, we show that differences between evGW+BSE and state-of-the-art TDDFT calculations using an OT-RSH functional can be explained almost entirely based on differences in how electron–hole interactions are described by the xc kernel of TDDFT and the BSE kernel, respectively. Eigenvalue differences as computed with evGW and DFT with an OT-RSH functional are almost identical.

The remainder of this article is structured as follows: we start by briefly reviewing the GW+BSE approach and report computational details and numerical convergence. We then show the effect of different DFT starting points and partial self-consistency on the excitation energies of the BC molecule. After this, we discuss our results for BCL a, b, c, d, and e and CL a and b, followed by a comparison with literature results and an in-depth discussion of differences between our GW+BSE and TDDFT results for BCL a.

Methods

GW+BSE Approach

In Green’s function-based many-body perturbation theory, the calculation of charged excitations, corresponding to electron removal and addition energies, is based on knowledge of the exact interacting single-particle Green’s function G, which can, in principle, be computed from a set of self-consistent integro-differential equations—Hedin’s equations—linking G to the electronic self-energy Σ, the screened Coulomb interaction W, the irreducible polarizability χ, and the vertex function Γ.[32] The lowest-order expansion of Σ with respect to W leads to the GW approximation, in which the electronic self-energy Σ = iGW.[33] Quasiparticle (QP) eigenvalues can be obtained by solvinghere, Vion is the ionic potential, VH is the Hartree potential, and εQP and φQP are QP energies and wavefunctions, respectively.

To avoid the high computational cost of a self-consistent solution of eq , the GW approach is commonly used within a one-shot scheme, in which G0 and W0 are constructed from a (generalized) Kohn–Sham (gKS) eigensystem obtained from a preceding DFT calculation. We use the notation G0W0@gKS to refer to G0W0 based on the gKS eigensystem (φgKS;εgKS) computed with the xc functional ExcgKS. In this approach, QP corrections are calculated to first order in Σ aswhere Vxc is the xc potential, and it is assumed that φQP ≈ φgKS.

While the G0W0 approach has been used with much success, in particular, for the calculation of band gaps and band structures of solids, a well-known and well-documented dependence on the gKS eigensystem used to construct G0 and W0 limits its predictive power.[34−36] Partial self-consistency in the QP eigenvalues can often mitigate this problem. In eigenvalue self-consistent GW, the gKS eigenvalues used to construct G and/or W are replaced with those from the output of a prior GW step; the self-energy corrections are then iterated until the QP eigenvalues converge. This approach, which we call evGW in the following (n refers to the number of iterations), has been shown to remove much of the starting point dependence for a range of different systems.[37,38]

The BSE is an equation for the two-particle electron–hole Green’s function and allows for the calculation of the polarizability including electron–hole interactions through the screened Coulomb interaction W. In practice, the BSE is usually solved by neglecting the frequency dependence of W. Within this static approximation, it can be written in a form equivalent to Casida’s equations of TDDFTwhere Ωs represents neutral excitation energies and (Xs, Ys) the corresponding eigenvectors.[19]A and −A represent resonant and antiresonant transitions that can be expressed asand that are coupled through B and −B, defined asfor singlet excitations. In these expressions, i and j are occupied and a and b are unoccupied states and (ia|bj) stands forNote that φQP = φgKS, whenever the G0W0 or evGW approaches are used to construct A and B.

Computational Details

Our calculations of charged and neutral excitations were performed using the GW+BSE and TDDFT implementation in the open-source MOLGW software package (version 2B), which relies on Gaussian basis functions.[39] We used the frozen-core approximation throughout, which changes excitation energies by less than 1 meV. We also employed the resolution-of-the-identity (RI) method, in order to reduce the calculation of four-center integrals to two- and three-center integrals. For BCL a, the RI changes the QP highest occupied molecular orbital (HOMO)–lowest unoccupied molecular orbital (LUMO) gap by less than 50 meV using a 6-31G basis set and BHLYP as a starting point, but we expect the effect of the RI to be even smaller for the larger basis sets used in the remainder of this article.[40] To further reduce the computational cost of the evaluation of the GW polarizability, we use the single-pole approximation (see the Supporting Information for details). The Tamm–Dancoff approximation, which corresponds to neglecting the B matrix elements in eq is not used as it consistently increases both GW+BSE and TDDFT results by ∼0.3 eV, in agreement with previous findings.[6,27]

We tested the influence of the Gaussian basis set size on HOMO–LUMO gaps and Q and Q excitations of BCL a (using a structure from ref (16)) with G0W0@BHLYP+BSE, considering seven different basis sets, namely, the Pople basis sets 6-31G, 6-311G, 6-311++G**, and 6-311++G(2d,2p), combined with the DeMon auxiliary basis set,[41] and the Karlsruhe basis sets def2-SVP, def2-TZVP, and def2-TZVPP and their corresponding auxiliary basis sets.[42]

Figures and S2 show the convergence of the HOMO–LUMO gap and Q and Q excitation energies as a function of the inverse number of basis functions, 1/Nbasis for GW+BSE and TDDFT, respectively (raw data are presented in Tables S1 and S2). We find that the HOMO–LUMO gap depends significantly more on 1/Nbasis than Q and Q excitation energies and that TDDFT results are less sensitive to the choice of the basis set than GW+BSE. Based on these tests, we use the 6-311++G(2d,2p) basis set for all calculations reported in the following. We estimate the error in the GW(+BSE) HOMO–LUMO gap and the Q and Q excitation energies by linearly extrapolating to an infinite basis set. By excluding the very small 6-31G and 6-311G basis sets from these fits, we obtain extrapolated values of 3.57 eV for the HOMO–LUMO gap, 1.11 eV for Q, and 1.81 eV for Q. We conclude that using the 6-311++G(2d,2p) basis set for all further calculations, we likely overestimate GW(+BSE) HOMO–LUMO gaps and Q and Q excitation energies by ∼0.1 eV with respect to the complete basis set limit. Conversely, the use of the single-pole approximation leads to a similar underestimation of the HOMO–LUMO gap and the Q and Q excitations, resulting in a fortuitous cancellation of errors.

Figure 1

Convergence as a function of the number of basis functions 1/Nbasis for (a) HOMO–LUMO gap and (b) Q and Q excitation energies, calculated with G0W0@BHLYP+BSE. Dashed lines represent a linear fit.

We test the effect of different xc functionals on our TDDFT and GW+BSE results. We use the LDA, two global hybrid functionals (B3LYP and BHLYP), and the range-separated hybrid (RSH) functional ωPBE. In RSH functionals, the Coulomb repulsion is separated into a short-range part and a long-range part, for numerical convenience expressed aswhere ω is called the range separation parameter. The ωPBE functional uses PBE exchange in the short range and the exact exchange energy in the long range, allowing for a self-interaction-free description at large electron–electron distances. We obtain the range separation parameter ω through the tuning procedure outlined in ref (43), in which ω is chosen such that the HOMO eigenvalue is as close as possible to the negative ionization potential both for the neutral and anionic system. Consequently, and by construction, the resulting HOMO–LUMO gap is a very good approximation to the fundamental gap of the neutral molecule. We use the Q-Chem code and a 6-31G(d,p) basis set for the tuning.[44] The tuned range separation parameters for all systems discussed in the following can be found in Table S3.

Results and Discussion

Bacteriochlorin

To validate our methodological setup and investigate the starting point dependence of the G0W0 approach and the effect of eigenvalue self-consistency on our calculated HOMO–LUMO gaps and excitation energies, we start by examining the BC molecule, for which GW+BSE results have been reported in ref (27). We use a BC structure from ref (27) and denote the lowest energy excitations Q and Q, according to the direction of their transition dipole moments. Table contains our calculated HOMO–LUMO gaps, Q and Q excitation energies, and oscillator strengths using TDDFT and several flavors of the GW+BSE approach. We find that, as expected, (generalized) Kohn–Sham HOMO–LUMO gaps show a large dependence on the xc functional, with LDA, B3LYP, and BHLYP leading to a significantly lower HOMO–LUMO gap and ωPBE leading to a HOMO–LUMO gap similar to the HOMO–LUMO gap calculated with G0W0 and eigenvalue-self-consistent evGW. In turn, Q and Q excitation energies from TDDFT are considerably less dependent on the xc functional than HOMO–LUMO gaps. In agreement with previous studies, we find that TDDFT overestimates the experimental values for Q and Q by up to ∼0.4 eV, depending on the xc functional.[27] TDDFT with the OT-RSH ωPBE is in best agreement with experiment, overestimating it by ∼0.2 eV for both excitations. We further find that G0W0@LDA+BSE underestimates Q by 0.4 eV and Q by 0.6 eV, whereas the use of a BHLYP and ωPBE starting point results in excitations within 0.1 eV of the experimental results. In accordance with prior studies, we observe that most of the starting point dependence of the G0W0+BSE results is inherited from the starting point dependence of the HOMO–LUMO gaps.[25]

Table 1

HOMO–LUMO Gaps, Q and Q Excitation Energies (in eV), and the Corresponding Oscillator Strengths, Γ and Γ, for BC Calculated with the 6-311++G(2d,2p) Basis Set

method	xc functional	H-L gap	Qx	Γx	Qy	Γy
TDDFT	LDA	1.38	2.04	0.18	2.39	0.03
	B3LYP	2.17	2.06	0.23	2.51	0.04
	BHLYP	3.27	1.93	0.28	2.55	0.04
	ωPBE	4.38	1.87	0.23	2.42	0.05
G0W0+BSE	LDA	4.15	1.21	0.09	1.67	0.04
	B3LYP	4.36	1.44	0.14	1.97	0.04
	BHLYP	4.56	1.67	0.19	2.23	0.05
	ωPBE	4.59	1.64	0.19	2.26	0.05
evGnW0+BSE	LDA	4.31	1.41	0.13	1.97	0.04
	B3LYP	4.43	1.54	0.16	2.12	0.05
	BHLYP	4.56	1.66	0.19	2.24	0.05
	ωPBE	4.57	1.62	0.18	2.27	0.04
evGnWn+BSE	LDA	4.42	1.51	0.17	2.21	0.05
	B3LYP	4.55	1.63	0.19	2.24	0.05
	BHLYP	4.60	1.69	0.20	2.27	0.05
	ωPBE	4.56	1.61	0.18	2.26	0.04
Expa			1.60		2.30

Data for bacteriopheophorbide from refs.[31,27]

In order to investigate the effect of eigenvalue self-consistency in the GW+BSE approach, we tested the effect of updating the eigenvalues in the construction of G only (evGW0) and of both G and W (evGW). Eigenvalue self-consistency in G alone only slightly changes the results as compared to G0W0. In contrast, full eigenvalue self-consistency largely eliminates the starting point dependence. In particular, using an LDA starting point results in excitation energies within 0.1 eV from experiment—similar to the ωPBE starting point, but at considerably reduced computational cost. In Table S4, we report similar results for the more complex pigment BCL a. In the remainder of this article, we therefore focus primarily on eigenvalue self-consistent results based on LDA and ωPBE starting points.

Excitation Energies of Bacteriochlorophylls and Chlorophylls

Next, we turn to reporting the vertical excitation energies of several members of the BCL and CL family of pigments. All structures were obtained from ref (45) and geometry-optimized using DFT as implemented in the Turbomole code with a def2-TZVP basis set and the B3LYP xc functional.[46] Atomic coordinates of all relaxed structures can be found in the Supporting Information. We used both LDA and ωPBE starting points for our evGW+BSE and ωPBE for our TDDFT calculations. Unlike the Q excitation of BCL a and b, which has significant oscillator strength, the Q excitation of BCL c–e is dark. Following ref (8), we therefore also compare our calculations with experimental results for the higher-energy B band.[31] We report the vertical excitation energies and corresponding oscillator strengths of the first six excitations of all pigments in Tables S6 and S7. In these calculations, we included a total of 20 excitations, in order to ensure that the higher lying excitations are well converged.

Table demonstrates that evGW+BSE is in excellent agreement with experiment for the entire family of BCL and CL molecules. The MAE is about 50 meV for the Q and Q and between 100 and 200 meV for the B excitation. Our evGW+BSE results also accurately reflect the spectral shifts of the Q excitation when comparing different BCL pigments with each other. For example, the BCL b molecule differs from BCL a through an ethyliden side group, which shifts the Q excitation by 40 meV to the red. This red shift is perfectly reproduced in our GW+BSE calculations. This is the first main result of this study. The second one is that our results are essentially independent of the DFT eigensystem used as an input for the GW+BSE approach: a computationally inexpensive LDA starting point results in the same level of agreement with experiment as the more tedious ωPBE calculation that involves a system-dependent tuning procedure for the range separation parameter ω. This is in stark contrast to TDDFT. TD-LDA leads to spurious excitations with charge transfer character in between Q and Q, as well as slightly above Q, depending on the structure, as discussed below and in the literature.[47] TDDFT with the optimally tuned ωPBE results in good agreement with experiment for all three excitations, albeit with slightly higher MAEs of 170, 50, and 250 meV for Q, Q, and B, respectively.

Table 2

Q, Q and First B Band Excitation Energy of BCLs and CLs Calculated Using a 6-311++G(2d,2p) Basis Seta

	GW@LDA+BSE			GW@ωPBE+BSE			TD-ωPBE			expb
molecule	Qy	Qx	B	Qy	Qx	B	Qy	Qx	B	Qy	Qx	B
BCL a	1.52	2.08	3.25	1.50	2.10	3.16	1.75	2.16	3.33	1.60	2.15	3.46
BCL b	1.48	2.07	2.95	1.45	2.09	3.05	1.69	2.15	3.19	1.56	2.14	3.37
BCL c	1.85	2.05	2.94	1.84	2.11	3.02	2.05	2.21	3.21	1.88		2.89
BCL d	1.90	2.15	2.94	1.89	2.21	3.05	2.08	2.29	3.19	1.90		2.93
BCL e	2.01	2.04	2.78	1.96	2.13	2.88	2.10	2.23	3.02	1.92		2.72
CL a	1.85	2.13	2.91	1.86	2.19	3.02	2.06	2.29	3.16	1.87	2.14	2.88
CL b	1.95	2.17	2.79	1.93	2.20	2.85	2.10	2.29	2.97	1.92	2.26	2.72
MAE	0.05	0.06	0.12	0.04	0.05	0.17	0.17	0.05	0.25

GW+BSE results are based on eigenvalue self-consistent evGW.

Experimental results in diethyl ether from ref (8).

In Figure , we plot the difference between our calculated results and experiment, averaged over all three excitations, to further highlight qualitative differences between evGW+BSE and TDDFT. For BCL a and BCL b, evGW+BSE, on average, underestimates experiment by ∼100 meV, whereas the average TDDFT deviation is close to zero, because TDDFT slightly overestimates the Q and Q excitations but underestimates the B excitation of these pigments. For all other BCL and the two CL molecules studied here, we consistently find that the average deviation of evGW+BSE is significantly smaller than that of TDDFT. Similar to our results for the BC molecule and to other benchmark studies of complex organic molecules,[6] TDDFT tends to overestimate all three excitations by between 200 and 300 meV. evGW+BSE is in much closer agreement with experiment for these pigments, on average, overestimating their excitation energies by less than 100 meV. We stress again that these results are independent of the DFT starting point, whereas our TDDFT results rely on a per-system tuning procedure.

Figure 2

Colored bars denote the average difference between calculated and experimental excitation energies for evGW@LDA+BSE (blue), evGW@ωPBE+BSE (red), and TDDFT (green) with ωPBE. The dots represent the maximum deviation in each case.

Our results are in excellent agreement with correlated excited-state methods for those systems for which such studies have been reported, primarily CL a and BCL a. ADC(2) excitation energies of the first three excitations of CL a reported by Suomivuori et al. are 1.97, 2.11, and 2.95 eV, within ∼0.1 eV of our evGW@LDA+BSE results.[13] In another study by the same authors, the ADC(2) Q excitation energy of histidin-ligated BCL a was reported to be 1.46 eV, again within ∼0.1 eV of our results, although it should be noted that the structures of ligated and free-standing BCL a slightly differ, leading to excitation energy differences of 10–30 meV at the ADC(2) level.[14] Furthermore, Sirohiwal et al. used a pair-natural orbital coupled cluster approach to study CL a and reported Q and Q excitation energies of 1.75 and 2.24 eV, respectively, for CL a, also within ∼0.1 eV of our GW+BSE results for these excitations.[15]

Not only the absolute energies of Q and Q excitations are important for understanding and predicting excitation energy and charge transfer in photosynthetic systems but also their relative energy difference, Δ, plays a role, in particular, in coupled systems of several pigment units. It is therefore reassuring that evGW+BSE predicts Δ in very good agreement with experiment, with a deviation of only 10 meV for BCL a, BCL b, and CL a and 120 meV for CL b for the LDA starting point and a slightly larger deviation of, on average, 60 meV for the ωPBE starting point. TDDFT based on ωPBE tends to underestimate Δ, by, on average, 110 meV for these four pigments. For BCL a, we also show in Table that Δ strongly depends on the xc functional used in the TDDFT calculations, primarily because of the strong dependence of the Q excitation on the amount of exact exchange, which can be seen by comparing the results based on the LDA (0% of exact exchange), B3LYP (∼23%), and BHLYP (50%). As before, evGW+BSE is in excellent agreement with experiment and almost independent of the underlying xc functional.

Table 3

Difference between Q and Q Excitation Energies (in eV) Using TDDFT and evGW+BSE for BCL a

method	xc functional	ΔQx–Qy
evGnWn+BSE	LDA	0.57
	B3LYP	0.54
	BHLYP	0.54
	ωPBE	0.59
TDDFT	LDA	0.25
	B3LYP	0.40
	BHLYP	0.64
	ωPBE	0.43
exp[8]		0.55

The experimental results reported in Tables and 3 are based on measurements in diethyl ether, whereas our calculations are for gas-phase molecules. To approximately account for the effect of the solvent, we extracted experimental reference values for Q and Q excitations from a study by Limantara et al.,[30] in which electronic absorption spectroscopy was used to obtain Q and Q for a large number of nonpolar and polar solvents at room temperature. This study reports regression lines for Q and Q excitations of BCL a as a function of R(n) = n2 – 1/n2 + 2, where n is the refractive index of the solvent. The extrapolated values for n = 1 (vacuum) are 1.68 eV (nonpolar) and 1.67 eV (polar) for the Q and 2.25 eV (nonpolar) and 2.21 eV (polar) for the Q excitation. Based on these regression parameters, we estimate that the experimental reference values in Table lie ∼50–70 meV below the gas-phase excitation energies. We also calculated the Q and Q excitation energies of BCL a with TDDFT (using ωPBE), approximating solvent effects with the COSMO approach as implemented in Turbomole. We used a dielectric constant of 4.33ε0 corresponding to the value in diethyl ether. COSMO red-shifts the Q and Q excitation energies by 70 and 50 meV, respectively, supporting our estimate. We conclude that solvent effects are small—within the numerical accuracy of our GW+BSE calculations—and do not change our main conclusions. Note that we also neglect the effects of temperature and the 0–0 vibrational energy contribution in our comparison with experimental results. Exact agreement of our calculated results with experiment is therefore not expected.

Bacteriochlorophyll a

In the remainder of this paper, we will use the BCL a molecule as a case study to compare to available computational literature results for this pigment, discuss the origin of differences between our evGW+BSE and TDDFT results, and comment on the effects of the choice of structure on excitation energies.

Comparison with RASPT2

For the Q and Q excitations of BCL a, we compare our GW+BSE and TDDFT calculations to multistate, second-order perturbation theory (RASPT2) calculations by Anda et al.[16,17] For this comparison, we use the molecular geometry reported in ref (16), which is a BCL a unit from the light-harvesting system LH2 of Rhodoblastus acidophilus. This structure was extracted from an experimental X-ray crystallographic structure of the LH2 complex (unit 302 within the structure 1NKZ in the RCSB Protein Data Bank).[48] The phytyl tail was truncated and replaced by a hydrogen atom, and no further geometry optimization was carried out. In the following, we will call this structure “A”. Our geometry-optimized version of “A”, which we relaxed using DFT as implemented in the Turbomole code with a def2-TZVP basis set and B3LYP,[46] will be called “R”. A visual comparison between “A” and “R” is shown in Figure . The large differences that we observe between these two structures are unsurprising, given that we perform our geometry optimizations without taking into account the protein environment in which BCL a “A” is embedded in vivo. Table shows our GW+BSE and TDDFT results for “A” in comparison with the RASPT2 excitation energies from refs.[16,17] We find, as before, that when eigenvalue self-consistency is used in GW, HOMO–LUMO gaps and Q and Q excitation energies differ by a maximum of 0.1 eV. Most notably, however, our GW+BSE excitation energies substantially differ from those calculated with RASPT2, with Q 0.4 eV and Q 0.5 eV lower than the RASPT2 result.

Figure 3

Overlay of structures “A” (red) and “R” (blue) in (a) top view and (b) side view.

Table 4

HOMO–LUMO Gaps and Q and Q Excitation Energies (in eV) for BCL a Structure “A”

method/basis set	xc functional	H-L gap	Qy	Qx
evGnWn+BSE/6-311++G(2d,2p)	LDA	3.62	1.17	1.90
	B3LYP	3.67	1.19	1.90
	BHLYP	3.72	1.23	1.92
	ωPBE	3.68	1.16	1.91
evGnWn+BSE/ANO-RCC-vDZP	LDA	3.76	1.38	2.18
TDDFT/6-311++G(2d,2p)	LDA	0.92	1.59	1.99
	B3LYP	1.60	1.64	2.17
	BHLYP	2.61	1.57	2.34
	ωPBE	3.70	1.48	2.02
RASPT2/ANO-RCC-vDZP			1.61	2.40

We find that about half of this difference can be traced back to the use of a smaller basis set (ANO-RCC-vDZP) in ref (16). Repeating our evGW@LDA+BSE calculation with the same basis, we obtain excitation energies of 1.38 eV for Q and 2.18 eV for Q. In line with previous studies, we also find that TDDFT with global hybrid functionals (B3LYP and BHLYP) results in similar excitation energies to RASPT2 for the Q excitation.[17,49] We hypothesize that this agreement is fortuitous. The optimally tuned RSH functional ωPBE has been shown to better describe singlet excitation energies of a wide variety of organic compounds as compared to global hybrid functionals[50−52] and is more than 0.1 eV lower in energy than the RASPT2 Q excitation energy. Similar trends have also been shown for CL a, where DFT-based multireference CI, just as TDDFT with global hybrid functionals, tends to overestimate experiment by ∼0.2 eV for the Q and Q excitation.[53] All in all, given that comparisons with experimental data are complicated for an in vivo structure such as “A”, we consider it most likely that our GW+BSE calculations underestimate the excitation energies of structure “A” by ∼0.1 eV, similar to our results for gas-phase BCL a (Table ). The remaining deviations could be attributed to the multireference character of the Q excitation[16] and the choice of the RAS.

We also note that our GW+BSE results reproduce the energetic order and relative energy differences of the Q excitation of other BCL units within the LH2 ring that RASPT2 predicts, when using the ANO-RCC-vDZP basis set. However, the use of the significantly larger 6-311++G(2d,2p) basis leads to substantially larger excitation energy differences between these units (Table S7). Finally, it is worth mentioning that our GW+BSE calculations reproduce the relatively large energy difference Δ ≈ 0.8 eV that RASPT2 predicts, whereas TDDFT excitation energy differences are much less sensitive to details of the structure, with Δ ≈ 0.5 eV (using ωPBE) similar to the gas-phase structure of BCL a. We speculate that a geometry optimization of structure “A” within its protein environment would result in a smaller Δ for both RASPT2 and GW+BSE.

Role of the Electron–Hole Kernel

We find that the difference between GW+BSE and TDDFT excitation energies can be traced back almost entirely to differences in how electron–hole interactions are described in both schemes. The Q excitation is primarily (∼90%) a HOMO → LUMO transition, and the HOMO–LUMO gaps, as calculated with DFT-ωPBE and evGW@ωPBE, differ by only 0.02 meV (Table ). In fact, the density of states (DOS) in the energy range relevant for both the Q and the Q excitations based on evGW@ωPBE and DFT-ωPBE eigenvalues are almost identical (see Figure ). To further test our hypothesis, we construct the statically screened Coulomb interaction W (see eq ) and solve the BSE based on a DFT-ωPBE eigensystem (instead of first computing QP eigenvalues using eq ). We obtain values for the Q and Q excitation that are only 20 meV higher and 40 meV lower than the full GW+BSE solution, respectively, for structure “A”. Similarly, for structure “R”, the results are within less than 10 and 50 meV for the Q and Q excitation, respectively. This observation confirms that differences between the GW+BSE and TDDFT excitation energies are primarily due to differences in the xc and the BSE kernel. Generally, the overestimation of excitation energies that we observe with TDDFT is in line with results for other organic π chromophores such as rhodamine and rosamine[54] and of phenothiazine dyes,[55] for which it has been linked to an insufficient treatment of differential electron correlation between the ground and excited states by most TDDFT xc kernels.[56]

Figure 4

(Generalized) Kohn–Sham (green), G0W0 (red), and evGW (blue) DOS calculated using the LDA and ωPBE. The HOMO energies are aligned to zero.

Charge Transfer Excitations with TD-LDA and GW@LDA+BSE

Finally, motivated by the excellent performance of evGW@LDA+BSE, we compare G0W0@LDA+BSE, evGWLDA+BSE, and TD-LDA results for structures “A” and “R” of BCL a. Figure shows the excitation spectrum calculated at these levels of theory. TD-LDA’s severe underestimation of charge transfer excitations is well known[9] and leads to spurious excitations with charge transfer character at energies between Q and Q for BCL a.[47] Our comparison of structures “A” and “R” shows that while the energy of Q and Q is changing only slightly when TD-LDA is used, the relative position of these spurious low-oscillator strength excitations depends strongly on the structure. G0W0@LDA+BSE results in a very different, albeit no more reassuring, picture. For both structures, the first excitation already appears at energies below or around 1 eV and its oscillator strength is considerably lower than with TD-LDA; for structure “A”, the oscillator strength of Q is even lower than that of Q. For structure “R”, excitations 2, 3, and 4 have similar, very low, oscillator strength. However, already at the G0W0@LDA+BSE level, no charge transfer excitations are found between Q and Q—a consequence of the inherent nonlocality of the BSE kernel. Finally, for both structures, eigenvalue self-consistency pushes all excitations to significantly higher energies and results in a quantitatively correct description of Q and Q.

Figure 5

First excitations for structures “A” (left) and “R” (right) as calculated with TD-LDA (top), G0W0@LDA+BSE (center), and evGW@LDA+BSE (bottom). Arrows indicate excitations with very low oscillator strength.

Inspection of the DOS calculated with DFT-xc, G0W0@xc, and GW@xc (xc = LDA, ωPBE) shown in Figure is instructive for understanding the contribution of eigenvalue differences to the TDDFT and GW+BSE excitation energies. The G0W0@LDA DOS underestimates the HOMO–LUMO gap and the energy difference between the HOMO and HOMO – 1. In contrast, there is virtually no difference between the HOMO, HOMO – 1, and LUMO energies as calculated with DFT-ωPBE, G0W0@ωPBE, evGW@ωPBE, and evGW@LDA. As expected, the DFT-LDA DOS is markedly different, not only underestimating the HOMO–LUMO gap but also significantly underestimating the energy differences between the HOMO – 1, HOMO – 2, and HOMO – 3. Notably, the spurious dark states between Q and Q that TD-LDA predicts have significant contributions from transitions involving these lower occupied states.

Conclusions

In this article, we performed a systematic first-principles study of the electronic structure and excitations of seven members of the (bacterio)chlorophyll family, which we validated through comparison with calculated and experimental literature results. The GW+BSE approach, when used in a partially self-consistent fashion, is in excellent agreement with experiment for excitations in the visible and near-ultraviolet part of the spectrum. GW+BSE also correctly predicts the energy difference between the low-energy Q and Q excitations of these pigments, relevant for the description of the coupling between pigment complexes, present in the light harvesting units and reaction centers of plants and bacteria and crucial for excitation energy and charge transfer. Most importantly, our results are almost entirely independent of the DFT eigensystem used as an input for the GW+BSE calculations. A computationally inexpensive LDA starting point leads to similar results as a more involved optimally tuned ωPBE starting point.

It should be noted that the GW approach, despite its implementation using Gaussian basis functions and the use of the RI approximation in MOLGW and other codes, remains a major bottleneck of these calculations due to its O(N4) scaling with system size. Furthermore, our results highlight that the GW approach, more so than DFT, requires careful convergence with respect to the basis set size. This limits its applicability to systems with a few (B)CL pigments at most, until algorithms with better scaling become more widely available.[57−59] Our study joins a growing number of results, demonstrating that the GW+BSE approach can accurately predict neutral excitations of complex molecules without empirical parameters.[29] With new approaches for combining GW+BSE with large-scale molecular mechanics simulations[28] and polarizable continuum embedding[60] emerging, an accurate prediction of excitation energy and charge transfer in complex molecular environments is within reach.