Accurate Computation of the Absorption Spectrum of Chlorophyll a with Pair Natural Orbital Coupled Cluster Methods.

The ability to accurately compute low-energy excited states of chlorophylls is critically important for understanding the vital roles they play in light harvesting, energy transfer, and photosynthetic charge separation. The challenge for quantum chemical methods arises both from the intrinsic complexity of the electronic structure problem and, in the case of biological models, from the need to account for protein-pigment interactions. In this work, we report electronic structure calculations of unprecedented accuracy for the low-energy excited states in the Q and B bands of chlorophyll a. This is achieved by using the newly developed domain-based local pair natural orbital (DLPNO) implementation of the similarity transformed equation of motion coupled cluster theory with single and double excitations (STEOM-CCSD) in combination with sufficiently large and flexible basis sets. The results of our DLPNO-STEOM-CCSD calculations are compared with more approximate approaches. The results demonstrate that, in contrast to time-dependent density functional theory, the DLPNO-STEOM-CCSD method provides a balanced performance for both absorption bands. In addition to vertical excitation energies, we have calculated the vibronic spectrum for the Q and B bands through a combination of DLPNO-STEOM-CCSD and ground-state density functional theory frequency calculations. These results serve as a basis for comparison with gas-phase experiments.

Introduction

Chlorophylls are photosynthetic pigments that play crucial roles in the absorption of sunlight, excitation energy transfer, and primary charge separation in photosynthetic organisms.[1−4] Their low-energy electronic transitions comprise the Q band (including the Q and Q transitions) and the B band, also known as the Soret band. Vibrational overtones accompany each major peak in the absorption spectra. The Q and Q transitions (subscripts denote the idealized polarization direction within the macrocycle plane) are the most important for determining the biophysical properties of chlorophyll species in particular and chlorophyll-containing photosystems in general.[5−12] The transitions are conventionally described by the four-orbital Gouterman model,[13,14] which, however, has known limitations. These arise from the fact that the electronic nature of actual excitations is more complex than that suggested by the Gouterman model, particularly for the B band, and because of the possible mixing of transitions due to vibronic coupling.[9,11,12,15−19]

The accurate calculation of the low-energy spectrum of the chlorophyll excited states has been a long-standing target of quantum chemistry.[20−22] A wide variety of computational investigations have been used for computing excitation energies and absorption profiles, ranging from time-dependent density functional theory (TD-DFT)[23−30] and DFT/multireference configuration interaction (MRCI)[31] to various wave function methods such as SAC-CI,[32] CC2,[33,34] and ADC(2).[34] Advanced multireference approaches have so far been applied only to bacteriochlorophylls.[35−37] TD-DFT remains a popular choice for photosynthetic chromophores because of its low computational cost and sometimes acceptable accuracy. However, questions have been raised about the performance of the approach for chlorophylls and similar systems.[26] In addition, the inadequacy of TD-DFT to properly predict charge-transfer states, and hence be reliably applied in the context of photosynthetic charge separation, is well known.[25,27,38] Wave function-based approaches are able to deliver systematically better results than DFT, albeit at increased computational cost; therefore, it is important to seek new ways for their practical usage.

In recent years, the domain-based local pair natural orbital (DLPNO) approach has been successfully applied to accelerate the methods of the coupled cluster (CC) hierarchy for ground-state calculations. Significant progress has also been made in the use of localization techniques for excited states. The first breakthrough in this regard was the implementation of the local CC variant of EOM-CCSD by Korona and Werner,[39] which was soon followed by a local CC2 method making use of the density fitting[40] and Laplace transform[41] techniques. The concept of state-specific pair natural orbitals was proposed by Helmich and Hättig[42] and was soon implemented for CC2 and ADC(2)[43] as well as for EOM-CCSD.[44] Excited-state approaches based on state-specific natural orbitals have also been proposed at the CC2[45] and ADC(2)[46] levels by Mester et al. In contrast, Peng et al. rely on state-averaged PNOs in their EOM-CCSD implementation.[47,48] Further details on this and other methods are provided in a recent review.[24] Among these methods, our group[49−51] has recently proposed a combination of ground-state DLPNOs[52,53] with the similarity transformed equation of motion (STEOM) approach of Nooijen and Bartlett,[54] which has the advantage that single excitations are decoupled from the doubles. This property allows for reducing the final diagonalization step to the size of the space of single (hole/particle) excitations. Thus, the final excitation space is similar to that of TD-DFT and other particle/hole theories. However, unlike TD-DFT, the effect of the doubles is retained, and a further implicit triple correction improves the description of charge-transfer states. This is the reason why frequently the results of STEOM-CCSD calculations surpass the accuracy of the EOM-CCSD method itself.[54] We emphasize that the STEOM-CCSD method and its DLPNO variant do not rely on a perturbative approximation to EOM-CCSD, as the popular and more approximate CC2 and ADC(2) methods do. One caveat of the method is that STEOM-CCSD introduces an active space of occupied and virtual orbitals. However, this active space is not a complete active space in the sense of multiconfigurational approaches, but it serves to identify the orbitals that dominantly participate in the electron depopulation and repopulation processes that accompany electronic excitations. This active space is not practically restricted to about 14 orbitals (as would be the case for a complete active space approach). Fortunately, recent work has demonstrated that the selection of this active space can be automatized.[55,56]

In order to account for vibronic transitions, we followed the ideas outlined by Baiardi et al.[57] in developing a time-dependent approach[58−60] that relies on DLPNO-STEOM vertical excitation energies (VEEs), transition dipole moments (TDMs) (Franck–Condon approximation[61−63]), and their derivatives (Herzberg–Teller effects[64]), whereas the ground- and excited-state geometries as well as harmonic vibrational frequencies are calculated at the DFT/TD-DFT level. This approach essentially amounts to calculating Fermi’s golden rule by Fourier transforming it into the time domain and applying the path integral formula for the harmonic oscillator to evaluate the resulting time-correlation function. The resulting composite method has been used with success to characterize and interpret a variety of spectroscopic measurements.[58−60,65−67] It is also worthwhile to mention a similar approach used by Barone and co-workers on a variety of chemical compounds.[68−74]

In this study, we apply the above DLPNO–STEOM-CCSD approach to all states that comprise the absorption spectrum of the chlorophyll a (Chl a) model (Figure ). The results, which can be considered as reference values for lower-level coupled-cluster approximations and for DFT approaches, are compared with the results obtained by a variety of other methods and with the experimental gas-phase results[75] on Chl a. Moreover, we demonstrate that treating vibronic coupling, including Herzberg–Teller effects,[74,76] is, in fact, required in order to obtain an accurate representation of the linear absorption spectrum of Chl a.

Figure 1

Structure of the gas-phase optimized Chl a model. The phytyl-carrying group at position 17 in the present calculations is replaced by a methyl group.

Methodology

All quantum chemical computations in this work were performed using the development version of the ORCA 4.2 program package.[77] The initial molecular geometry of a Chl a molecule was extracted from the crystal structure of photosystem II (PDB ID: 3WU2).[78] As a first step, we built different variants of Chl a with respect to the length of the group at position 17 of the chlorin ring, either with the phytyl chain terminated after the ester group (capped with a methyl group) or with the group replaced by a methyl group. Both models were fully optimized with the B3LYP hybrid functional.[79] No variations in macrocyclic ring curvature were observed. VEEs (obtained with DLPNO–STEOM-CCSD, see Table S1) are invariant with respect to the phytyl chain length of the model. Based on these results, we chose the model with a minimal representation of the chain at position 17 of the chlorin ring that is capped with a methyl group (Figure ). Final gas-phase geometry optimizations were performed using CAM-B3LYP[80] along with def2-TZVP[81−83] basis sets, with the inclusion of dispersion corrections via the D3(BJ) scheme.[84,85] The choice of the functional was motivated by the observation that it provides better balance in the description of low- and high-energy states compared to other functionals, as well as in order to ensure the consistency of methods between geometry optimization and calculation of ground- and excited-state frequencies. In addition, a recent study[86] demonstrated an improved performance of CAM-B3LYP over B3LYP for the prediction of the vibronic spectrum for both Chl a and Chl f. The RIJCOSX approximation[87,88] was used to speed up the calculations. Tight thresholds were employed for the self-consistent field (SCF) and optimization (verytightscf and verytightopt keywords in ORCA nomenclature). The resulting geometries had no imaginary frequencies. TD-DFT VEEs were computed using range-separated CAM-B3LYP[80] and ωB97X-D3(BJ),[89−92] as well as hybrid functionals B3LYP,[79] PBE0,[93] and BHandHLYP,[94] and the double-hybrid B2PLYP[95,96] and ωB2PLYP.[97] A total of 10 roots were computed. This number of roots was selected so as to ensure that all excitations that could possibly correspond to the B band are fully recovered. All results are reported in detail in the Supporting Information.

The DLPNO–STEOM-CCSD calculations were performed with the def2-TZVP basis set (total of 1464 basis functions). Calculations with larger basis sets and basis sets with diffuse functions could not be successfully completed for the present system. A total of 10 roots were requested depending on the requirement for a particular given case. No change in excitation energies was observed when a higher number of roots were computed; however, the cost of the calculation increases steeply beyond this point. The calculations were performed with the “TightPNO” setting, the TCutPNOsingles keyword set to 6.6 × 10–10 and the active space selection keywords “Othresh” and “Vthresh” set to 1.0 × 10–3. As an indication of computational cost, one such calculation for 10 roots took 5 days 7 h (wall clock time) on 8 cores of an Intel Xeon E5-2640 v3 CPU. Slightly lower thresholds in active space selection (“Othresh” and “Vthresh” set to 5.0 × 10–3) reduce the total computational time to 3 days 20 h with no change in the computed VEEs. Additionally, if the basis set is reduced to def2-SVP, the required time is reduced to less than 12 h with no noticeable loss in accuracy (see the Supporting Information). At this level, the cost is about 2 to 3 times that of TD-DFT calculations with double-hybrid functionals. Nevertheless, all the reported values in the present work utilize the more conservative thresholds and large basis sets stated above.

The vibronic spectrum of the excited states associated with both x- and y-polarized components of the Q and B states was computed using the excited-state dynamics (ESD) module of the ORCA package. One of the basic requirements for the computation of the vibronic spectrum using this methodology is the ground- and excited-state Hessians, which were computed using the CAM-B3LYP/def2-TZVP level of theory. The influence of the vibronic coupling on the absorption spectra was computed using the Herzberg–Teller effects. The temperature was set to 300 K during the calculations to match with the experimental conditions. We have used a Gaussian line shape to model the absorption spectra together with a computed linewidth equal to the experimental full width at half-maximum. The vibronic progression (sticks format) was computed using an intrinsic linewidth of 0.1 cm–1 (using the INLINEW keyword of the ORCA_ESD module).

A direct comparison of the computed vertical excitation or the 0–0 energies with the experimental band maximum is conceptually incorrect,[98] as the experimental bands contain contributions from the vibrational levels and geometric relaxation in the excited states. In this study, we report an effort to derive “quasi-experimental” vertical excitation and 0–0 energies to allow for a direct comparison of the experimental results with the quantum chemical calculations (Table S2). Our approach is based on band shape calculations using the vertical gradient (VG) approach for the x- and y-polarized components of the Q and B states.

Results and Discussion

Experimental Gas-Phase Spectra and the Gouterman Model

Experimental data on the absorption spectra of chlorophylls and bacteriochlorophylls refer almost exclusively to the species in solution. Depending on the solvent, magnesium may coordinate either one or two solvent molecules axially, with a concomitant geometric distortion of the chlorin ring. Additionally, solvent molecules may directly interact with the ring substituents. Both effects result in the alterations of the electronic properties and hence of the absorption spectrum of the molecule. The treatment of the solvation effects represents a challenge to theoretical approaches as significant as the intrinsic electronic structure problem. Indeed, most theoretical studies of chlorophylls have dealt with solvated, or in general axially coordinated, systems.[16,99−108] However, to disentangle the two problems and be able to focus on the electronic structure itself, it is useful to have experimental gas-phase values that can serve as reference for quantum chemical calculations. Such gas-phase spectra have been reported only recently and only for Chl a and Chl b using action spectroscopy (Figure ).[28,29,75,109] However, gas-phase experiments with neutral chlorophyll molecules are challenging to perform as these molecules decompose easily.[28] Therefore, the experiments are generally performed using a charge tag approach, where a cationic species is used alongside the chlorophyll molecules. Variable sizes of charge tags are used to observe the influence on the chlorophyll absorption profiles. Interestingly, a study by Wellman and Jockusch[109] demonstrated the influence of the alkali metal-tagged Chl a on the Soret region. The experimental data[75] mentioned in the present study were obtained using tetramethylammonium cations.

Figure 2

Gas-phase action spectra (at 300 K) of charge-tagged Chl a after one-photon absorption. The original data are from the study of Gruber et al.(75) Band maxima in the Q band are observed at 1.94, 2.07, 2.23, and 2.38 eV (the first and third are assigned to Q and Q), whereas B-band maxima are observed at 3.08 and 3.38 eV (assigned to B and B).

Before proceeding to the presentation and discussion of the computational results, it is useful to recall that absorption spectroscopy of porphyrin-like macrocyclic compounds is conventionally discussed in the context of the Gouterman model, which emphasizes excitations among the four frontier molecular orbitals HOMO – 1, HOMO, LUMO, and LUMO + 1, creating the Q and B bands (Figure ). The Gouterman model originates from the description of the porphyrin ring, where these four orbitals transform as the accidentally degenerate pair a2u and a1u, plus the unoccupied e degenerate pair (eg and eg) in the idealized D4 point group. As Chl a loses all of the formal symmetry elements of the parent porphyrin system, all degeneracies in the frontier MOs are lifted, leading to the splitting of transitions into distinct components. These are conventionally labeled as x and y according to the transition dipole orientation in the macrocyclic plane, although the situation is not as clear-cut in chlorophylls and bacteriochlorophylls as implied by the formal analogy of the labels with the porphyrin case. The y-labeled features are principally associated with HOMO → LUMO and HOMO – 1 → LUMO + 1 transitions, whereas the x-labeled features are mostly associated with HOMO – 1 → LUMO and HOMO → LUMO + 1 transitions.

Figure 3

Chl a frontier molecular (Hartree–Fock) orbitals associated with the Q- and B-band transitions according to the Gouterman model. Although no symmetry elements or degeneracies are present, the orbitals resemble at least qualitatively those of the parent porphyrin system. In parentheses, the symmetry labels of the parent porphyrin (D4 point group) and chlorin (C2 point group) systems are provided for reference. The arrows show which orbital pairs contribute mostly to the corresponding excitations, with the labels adhering to the nomenclature originally defined for the parent porphyrin system to indicate the approximate polarization of individual transitions.

As we will discuss in more detail in a later section presenting the band shape calculations, we have attempted to back-correct the experimental data shown in Figure . This gives us a set of VEEs, 0–0 energies, and band maxima that we will refer to as “quasi-experimental” in the following. These energies will make it easy to compare the calculated VEEs from various quantum chemical approaches with the quasi-experimental energies. Table presents the results of this analysis that are used as reference in the following.

Table 1

Quasi-Experimental 0–0 Energies and VEEs Associated with the Individual Excited States (All Values in eV)

	Qy	Qx	Bx	By
band max	1.94	2.23	3.08	3.38
0–0 energy	1.93	2.22	3.07	3.37
VEE	1.99	2.30	3.12	3.38

Results of DLPNO–STEOM-CCSD Calculations

DLPNO–STEOM-CCSD calculations predict two excited states in the Q region, at 1.75 and 2.24 eV, which can be identified with Q and Q, respectively. Two higher energy transitions at 3.17 and 3.39 eV fall within the Soret region and are assigned to the B and B components, respectively. These assignments are straightforwardly based on the nature of the molecular orbitals involved (Figure ), the computed configuration interaction coefficients for the excitations (Table ), the corresponding natural transition orbitals (NTOs) (Figure ), and the orientation of the TDM vector relative to the macrocyclic plane (Figure ). The Q transition mainly involves the HOMO → LUMO excitation, which largely conforms to the Gouterman model, whereas Q is composed of a smaller contribution from the HOMO – 1 → LUMO excitation combined with the HOMO → LUMO + 1 transition. It is noted that both LUMO and LUMO + 1 have some contribution from the vinyl group at position 3 of the chlorin macrocycle. The calculated oscillator strengths are consistent with the above assignment, with Q higher than Q (fosc = 0.22 vs 0.04) and the Soret bands having significantly higher oscillator strengths of fosc = 1.30 and 1.23 for B and B, respectively. The relevant information of the higher lying excited states is provided in the Supporting Information (Table S3 and Figure S1).

Table 2

VEEs (E, in eV) of Gas-Phase Optimized Chl a Using DLPNO–STEOM-CCSDa

excited state	E (eV)	frontier molecular orbitals	% contribution	|μ| (Debye)
S1 (Qy)	1.746 (0.219)	HOMO → LUMO, HOMO – 1 → LUMO + 1	0.69, 0.16	5.76
S2 (Qx)	2.240 (0.040)	HOMO – 1 → LUMO, HOMO → LUMO + 1	0.59, 0.28	2.15
S3 (Bx)	3.172 (1.304)	HOMO → LUMO + 1, HOMO – 1 → LUMO	0.49, 0.22	10.34
S4 (By)	3.397 (1.231)	HOMO – 1 → LUMO + 1, HOMO → LUMO	0.61, 0.15	9.79

Oscillator strengths are shown in parenthesis. The major contributions in the canonical basis associated with the first four electronic states are described. The magnitude of the TDM (μ, in Debye) associated with the excited state is also provided.

Figure 4

NTOs associated with the first four electronic excited states computed by DLPNO–STEOM-CCSD. Isosurface maps are drawn with a cutoff value of ±0.3 au.

Figure 5

Orientation of TDMs associated with the Q- and B-band transitions computed using DLPNO–STEOM-CCSD. The TDM vectors of Q, Q, B, and B deviate (with reference to the y-axis) by 10.1, 45.3, 71.3, and 6.4°, respectively.

Compared to the quasi-experimental VEEs, the computed excitation energy for the lowest excited state (Q) appears to be red-shifted by 0.24 eV. This represents the largest deviation from the reference among the computed values. On the other hand, the second excited state (Q) matches well with the quasi-experimental value of 2.30 eV. In the B-band region, the computed values for both the B and B bands agree very well with the quasi-experimental values. It is particularly encouraging that the new method delivers a balanced description of all relevant transitions. Only the apparent slight underestimation of the lowest energy excitation (i.e., Q) seems to require further analysis given the importance of the lowest excited state in photosynthetic processes. A similar red shift has been observed in a STEOM-CCSD study[110] of free base porphyrin. In addition to a possible yet-to-be determined methodological origin, one possible reason is the reference geometric structure used for the calculations. As discussed in greater detail below, the results for certain transitions can be highly sensitive to this choice. Another possible reason for the apparent divergence of the computed Q from the quasi-experimental reference value could be related to the charge-tagged approach employed in the experiments: depending on the specificity of the binding of cationic species around Chl a, electrochromic shifts can be produced in the absorption spectra,[109,111] and hence the low-energy part of the action spectrum may deviate from the net gas-phase spectrum of Chl a. However, such deviations are predicted to be small.[30]

Accurate determination of the TDM and its orientation associated with the low-energy excited states of the biological chromophores[112] hold the key to model important energy-transfer processes. Based on our DLPNO–STEOM-CCSD calculations, the TDM associated with Q slightly deviates from the xy plane and is not perfectly aligned to the y axis (Figure ). Interestingly, similar results were obtained in the experimental study (polarization-resolved femtosecond visible pump–infrared probe spectroscopy) performed with d8-toluene as the solvent.[113] In addition, we find the TDMs associated with Q and Q to be nonorthogonal,[9,12] with an angle of 35° between them. Similarly, the TDMs associated with B and B were found to be nonorthogonal, with a mutual angle of 65°.

Comparison of DLPNO–STEOM-CCSD to the Literature Results

In comparison with other wave function methods, SAC-CI[27,32] predicts two excitations each in the Q-band and B-band regions, whereas non-Gouterman-type excitations (N bands) were found nearly degenerate with the B component (0.01 eV above B), consistent with the present DLPNO–STEOM-CCSD results. The lowest excited state (Q) obtained from the present DLPNO–STEOM-CCSD approach is within the range previously reported by SAC-CI,[27,32] that is, from 1.75 to 1.81 eV, depending on the geometry used and the computational details.

Literature values for the first excited state obtained with other approximations to the CC theory, such as CC2 and ADC(2), range from 1.99 to 2.26 eV,[34] but the performance of such approximate methods for higher excited states that comprise the complete spectrum remains to be evaluated. A collection of literature values is summarized in the Supporting Information (Table S4). It is important to note that most past studies that applied wave function methods to this problem have used diverse geometric models, in some cases obtained directly from a protein crystal structure without further optimization[32] (a practice that is strongly discouraged[114,115]) and in other cases optimized at the levels of theory that perhaps would not have been used today. A direct comparison of correlated wavefunction-based methods may not be reliable in the absence of results obtained with a common high-quality geometric reference, converged technical settings, and full coverage of all excited states of interest.

Non-negligible variations in the VEEs were noticed with respect to the DFT functional employed for the ground-state geometry optimization. Two important structural differences can be observed: PBE-[116] and B3LYP-optimized geometries produced more in-plane vinyl conformations; however, CAM-B3LYP preferentially rotates the vinyl group more out of the plane vinyl rotation (Figure S2). On the other hand, CAM-B3LYP geometries have slightly shorter Mg–N distances compared to PBE- and B3LYP-optimized geometries. Single-point DLPNO–STEOM-CCSD calculations on the CAM-B3LYP-optimized geometry produce more blue-shifted VEEs for both the Q and B bands (by 0.1–0.2 eV) compared to the B3LYP and PBE ground-state geometries. As noted earlier, the first non-Gouterman-type excitation lies 0.13 eV above B using the CAM-B3LYP geometry. However, this excitation is found to be nearly degenerate with the B excitation with the B3LYP-optimized geometry and falls in between B and B in case of the geometry optimized with the PBE functional (Table S5). Therefore, small variations in geometry lead to non-negligible changes in the energetics of the low-energy excitations and lowering of non-Gouterman-type excitations, which further complicate the B-band description; however, these changes do not alter the fundamental description of the spectrum.

To probe the effect of the vinyl group rotation on the structure–property correlations, we computed VEEs using DLPNO–STEOM-CCSD on structures with exactly the same macrocyclic geometry, where the dihedral angle associated with the vinyl group rotation is scanned from 0 to 90° in 15° steps (Table S6). We found that vinyl rotation does not impact the x- and y-polarized Q and B components but only affects the stabilization of a non-Gouterman state. Therefore, it is the variation produced by the various DFT functionals in the bond distances of the macrocyclic ring that produces the spectral shift in the Q and B bands. Beyond the implications for understanding the biophysical role of specific chlorophylls, this has consequences for the evaluation of theoretical methods used for excitation energies, as the choice of geometry becomes a critical parameter in assessing the quality of theoretical methods, affecting the perceived error of a given method for the calculation of excited states.

Comparison of DLPNO–STEOM-CCSD Calculations with DFT

Table compares the DLPNO–STEOM-CCSD values with the excited-state energies predicted by seven diverse functionals using the same geometry for Chl a. All functionals predict two excitations in the Q region, and the overarching impression from the TD-DFT results is the tendency to produce blue-shifted values for Q, and even more so for Q, for which range-separated functionals perform surprisingly poorly.

Table 3

VEEs (in eV) Associated with the x- and y-Polarized Components of the Q and B Bands of Chl a Computed with DLPNO–STEOM-CCSD and TD-DFT, Compared to Quasi-Experimental Valuesa

method	Qy	Qx	Bx	By
DLPNO–STEOM-CCSD	1.75 (0.219) [1]	2.24 (0.040) [2]	3.17 (1.304) [3]	3.40 (1.231) [4]
PBE0	2.19 (0.248) [1]	2.40 (0.020) [2]	3.26 (0.509) [4]	3.45 (0.786) [7]
B3LYP	2.16 (0.237) [1]	2.35 (0.019) [2]	3.17 (0.481) [4]	3.35 (0.731) [7]
BHandHLYP	2.18 (0.269) [1]	2.58 (0.021) [2]	3.50 (0.966) [3]	3.80 (0.899) [5]
CAM-B3LYP	2.14 (0.239) [1]	2.53 (0.029) [2]	3.44 (0.885) [3]	3.71 (0.895) [5]
ωB97X-D3(BJ)	2.09 (0.222) [1]	2.71 (0.047) [2]	3.56 (1.058) [3]	3.93 (0.730) [5]
B2PLYP	2.12 (0.262) [1]	2.23 (0.019) [2]	3.17 (0.905) [3]	3.27 (0.748) [4]
ωB2PLYP	2.04 (0.233) [1]	2.49 (0.048) [2]	3.45 (1.158) [3]	3.78 (0.704) [4]
quasi-exp. VEE	1.99	2.30	3.12	3.38

Values in parenthesis indicate oscillator strengths. Numbers in brackets represent the rank of the root based on increasing excitation energies; DFT methods often predict additional spurious low-energy states (see the Supporting Information for complete data).

Unlike DLPNO–STEOM-CCSD, all DFT functionals that do not incorporate perturbational correlation contribution (double-hybrids) predict more than two excitations in the Soret region. Besides the predicted B and B, these additional non-Gouterman-type[27] excitations have a charge-transfer character with diminishing oscillator strengths. Similar observations were reported earlier, and some of these excitations were labeled as N transitions.[27] They might be considered as artifacts or “ghost states”, at least in the sense that they appear unphysically low in energy. We observe three categories of functionals giving distinct results. PBE0 and B3LYP predict three non-Gouterman-type excitations, one between Q and B and two between B and B (Table S7). For CAM-B3LYP, ωB97X-D3(BJ), and BHandHLYP, the results are better as they predict the first non-Gouterman-type excitation situated right in between B and B (Table S8 and Figures S3–S5). In contrast, DLPNO–STEOM-CCSD predicts such excitations only above the B band (see Table S3). Double-hybrid functionals B2PLYP and ωB2PLYP provide a similar description of the ordering of excited states as that obtained by the DLPNO–STEOM-CCSD computations (Tables S9 and S10). However, both double-hybrid functionals predict the first N transitions to be nearly isoenergetic with the B band (a difference of 0.003 and 0.045 eV in the case of B2PLYP and ωB2PLYP, respectively). Overall, B2PLYP performs bests among the DFT functionals and provides results of practically equal quality to the DLPNO–STEOM-CCSD results. In contrast, CAM-B3LYP, ωB97X-D3(BJ), BHandHLYP, and ωB2PLYP functionals overestimate the B-band excitation energies and suffer from the presence of ghost states. Therefore, with the important exception of B2PLYP, no other DFT functional can provide a one-stop-solution, in the framework of TD-DFT, for all low-energy absorption features of Chl a. It is stressed that the lack of balance and the occurrence of ghost states are characteristic of TD-DFT and are by no means unique to the chlorophyll system.[38,117,118]

Band Shape Calculations

Understanding the features of the Q band of biochromophores is central for understanding various photophysical phenomena related to the protein–pigment interaction.[119,120] In general, the Q-band regions of important chromophores such as bacteriochlorophyll a and pheophytin a are easier to interpret (in terms of Q and Q assignment) because of their large Q–Q energy gaps.[9] The situation is more complicated for Chl a because of the small Q–Q energy gap and the difficulty in assigning the Q origin. The most recent recommendations regarding the Q origin by Reimers et al.(9) are based on a model which includes the vibronic picture to obtain satisfactory results. In addition, they report that Q and Q mix, influencing exciton transport and coherence. Q and Q mixing was reported to depend on the solvent, as solvent molecules can ligate axially,[103] form peripheral hydrogen bonds, influence the macrocyclic ring curvature, and so forth. These structural perturbations collectively can influence the Q-band properties and the overall vibronic structure. Here, we present an intrinsic vibronic picture of the Q band of Chl a in vacuo for comparison against the experimental spectrum.

As part of this investigation, we have first used the comparison of the computed and experimental band shapes to estimate the shifts required to adjust the experimental band maxima to VEEs. On the basis of these data, we have determined the energy difference between the band maxima, the vertical excitation, and the corresponding 0–0 energy. We found that most band maxima are slightly red-shifted compared to the VEE for each excited state under consideration; however, the extent of the red shift is slightly different for each excited state (Table S2). Subsequently, we used these obtained shifts for each individual excited state and applied them to the band maxima of the experimental spectrum as a correction, which allows us to extract “quasi-experimental” VEEs and 0–0 energies that can be compared to the computed values (see Table ).

In the following, we describe how we have employed band shape calculations in combination with the DLPNO–STEOM-CCSD energies to analyze the spectrum of Chl a. We have used two distinct approximations: the VG approach and the adiabatic Hessian after-step (AHAS) approach. Both methods require an optimized ground-state geometry, whereas the excited-state structure is approximated by a single augmented Hessian optimization step. The VG approach assumes that the ground-state Hessian adequately approximates the excited-state one, whereas AHAS relies on an explicitly computed Hessian for the excited state. For the present calculations, all Hessians were computed with CAM-B3LYP, and the VEEs and TDMs were those of DLPNO–STEOM-CCSD.

Both approaches show significant overlap between the vibronic progressions associated with the Q and Q excitations, and the overall spectrum is in very good agreement with the experimental spectrum. The main difference between the VG and AHAS results is that the AHAS approach produces an amplified shoulder in the region next to Q (Figures S8 and S9). On the other hand, both approaches produce equivalent results for Q and the respective vibronic progression. The shift results in an increased overlap, but even though there are vibrational modes that appear at the same region in the vibronic progression of Q and Q, there are no modes with high intensities among them. Overall, the VG result approximates better the experimental envelope and is to be preferred over the much more expensive AHAS scheme (Figures S6 and S7). It is also important to note here that previous studies[121] using the VG approach in the context of the Chl a absorption spectra yielded satisfactory results.

The complete vibronic progression associated with the Q and Q excitations using the VG approach compared to the experimental spectrum is shown in Figure . For simplicity, the spectral features observed in the experimental spectra are broken down into four regions (P, i = 1–4). The computed vibronic progression associated with the Q excitation typically covers all four regions, with a higher contribution in the P1 and P2, whereas its tail spans regions P3 and P4. P1 corresponds to the 0–0 transition, with a minor contribution from the low-frequency out-of-plane chlorin deformation modes and substituent rotations, as well as low-energy in-plane bending modes. The P2 region corresponds to the vibronic progression, with the contribution from various in-plane stretching–bending modes, whereas the accompanying shoulder region corresponds to the in-plane chlorin deformations combining the C–C and C–N stretching and C–H bending modes. No fundamentals of the first excited state occur in the P3 region, which is associated with the Q excitation. The P4 region encompasses modes of the second excited state ascribed to the P2 shoulder from the first excited state. Based on the present results, the Herzberg–Teller contribution has a noticeable effect on the overall form of the P3 and P4 regions (see Figures S10 and S11), whereas P1 and P2 can be adequately described in the Frank–Condon picture.

Figure 6

Computed vibronic spectrum of the Q band using the VG approximation. The spectrum associated with the Q excitation (i.e., S0 → S1) is shifted by 0.18 eV to align with the experimental band maximum. The “raw” band shape is depicted in Figure S5. Various absorption maxima observed in the experimental spectra are labeled accordingly.

Our calculations do not allow us to directly discuss the mixing of states, as we are operating in an adiabatic picture that excludes the effects of nonadiabatic coupling matrix elements. Although this is a limitation to be addressed in future developments, we observe that the experimental spectrum can be reproduced sufficiently well without explicitly accounting for the mixing of states. Therefore, we may suggest that mixing is not a dominant feature in the gas-phase absorption spectra of Chl a.

The vibronic spectrum associated with the B and B excited states of the B band is presented in Figure . The B band is significantly less well resolved experimentally compared to the Q band; therefore, we do not consider it feasible to delve into detailed assignments. However, it is important to stress that the overall shape of the spectrum is reproduced quite accurately by the present simulations. This includes not only the two principal features, whose separation is only slightly underestimated, but also the shoulders that are attributed to the vibrational progressions of B and B.

Figure 7

Computed vibronic spectrum of the B band using the VG approximation.

The present study has focused strictly on gas-phase Chl a models. As a final note, we would therefore like to point out that the separation between the Q and Q band maxima, which Figure indicates to be ca. 2500 cm–1 (ca. 0.3 eV), is suggested to decrease upon solvent coordination. For example, the vibronic model developed by Reimers et al.(9) predicted this energy gap to be 0.2 eV in ether, 0.12 eV in pyridine, 0.10 eV in MeOH/EtOH mixture. A recent study by Ogilvie and co-workers[12] showed a significant overlap between the vibrational modes associated with the Q side-band and Q band for an axially isopropanol-coordinated Chl a (Q–Q gap of 0.14 eV). Based on their study, this overlap was evident as the excitation energy of S2 decreases as a result of the axial ligation due to HOMO – 1 destabilization. On the same principle, the heterogeneous protein environment can also modulate the Q–Q energy gap and other properties of the chromophores[122] through axial ligation,[78] hydrogen bonding,[123,124] macrocyclic ring distortion,[125,126] and electrostatic effects.[119] The vibrational modes of Chl a are also fine-tuned in a protein matrix for efficient energy transfer and charge separation. The vibronic picture shown above can thus be considered as a reference point for gas-phase Chl a and as a guide for comparison with environmentally or chemically modulated forms of the system.

Conclusions

In conclusion, the present study demonstrates that the DLPNO–STEOM-CCSD approach is eminently useable for systems as large and complex as chlorophylls and provides accurate results in reasonable turnaround times and in a black box fashion. It is noted that in terms of cost, the method in the current implementation is only a few times more expensive than double-hybrid DFT. Importantly, DLPNO–STEOM-CCSD offers a balanced treatment of both the low-energy and high-energy regions of the spectrum, in contrast to any of the investigated DFT methods with the possible exception of B2PLYP. The study serves as a springboard for further large-scale applications that will use DLPNO–STEOM-CCSD to calculate the absorption features of photosynthetic pigments in solvents or within proteins. The treatment of the environment, axial ligation, ring deformation, and peripheral hydrogen bonding, among others, will all affect the excited-state profile of the chromophore. A reliable electronic structure method such as the one described here will be essential for ensuring that the fundamental physics of the system is correctly reproduced.