Optical Excitations of Chlorophyll $a$ and $b$ Monomers and Dimers

TL;DR

This study develops efficient computational methods using LCAO-based TDDFT to accurately predict optical absorption spectra of chlorophyll monomers and dimers, aiding the understanding of light harvesting in biological systems.

Contribution

It introduces a computationally efficient LCAO-TDDFT approach that reproduces plane wave results and semi-quantitatively matches experimental data for chlorophyll systems.

Findings

LCAO-TDDFT methods reduce computational effort compared to plane wave approaches.

Inclusion of GLLB-SC correction improves agreement with experimental photoinduced dissociation data.

Methods enable first-principles optical excitation calculations in large biomolecular systems.

Abstract

A necessary first step in the development of technologies such as artificial photosynthesis is understanding the photoexcitation process within the basic building blocks of naturally-occurring light harvesting complexes (LHCs). The most important of these building blocks in biological LHCs such as LHC II from green plants are the chlorophyll $a$ (Chl $a$) and chlorophyll $b$ (Chl $b$) chromophores dispersed throughout the protein matrix. However, efforts to describe such systems are still hampered by the lack of computationally efficient and accurate methods that are able to describe optical absorption in large biomolecules. In this work we employ a highly efficient linear combination of atomic orbitals (LCAOs) to represent the Kohn--Sham (KS) wave functions at the density functional theory (DFT) level and perform time dependent density functional theory (TDDFT) in either the reciprocal space and frequency domain (LCAO-TDDFT-$k$-$\omega$) or real space and time (LCAO-TDDFT-$r$-$t$) calculations of the optical absorption spectra of Chl $a$ and $b$ monomers and dimers. We find our LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ calculations reproduce results obtained with a plane wave (PW) representation of the KS wave functions (PW-TDDFT-$k$-$\omega$), but with a significant reduction in computational effort. Moreover, by applying the GLLB-SC derivative discontinuity correction $\Delta_x$ to the KS eigenenergies, with both LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ methods we are able to semi-quantitatively reproduce the experimentally measured photoinduced dissociation (PID) results. This work opens the path to first principles calculations of optical excitations in macromolecular systems.

Tables

Table 1. Table 1: Q band energies in eV and wavelengths in nm for Chl  a and b monomers and dimers with and without the 1 + charge tag.

method	structure	(eV)	(nm)		(eV)	(nm)
		Chl a			Chl b
PIDa	exp$\cdot$1+	1.97	629		1.99	623
$\frac{4\pi \omega }{c}\mathrm{Im}{\left[\alpha \left(\omega \right)\right]}^b$	full$\cdot$1+	2.08	596		2.15	577
	full	2.04	608		2.13	582
	full$\cdot$1+	1.96	633		1.87, 2.02	662, 615
$ℱ{\left\{d_m\left(t\right)\right\}}^c$	cut$\cdot$1+	1.95	636		1.89, 2.00	657, 621
	cut	2.01	618		1.96	631
	full$\cdot$1+	2.13, 2.30	582, 540		2.21, 2.31	562, 537
$\mathrm{Im}{\left[\epsilon \left(\omega \right)\right]}^d$	cut$\cdot$1+	2.13, 2.27	583, 545		2.21, 2.29	562, 542
	cut	2.21, 2.37	561, 524		2.32, 2.41	534, 515
$\mathrm{Im}{\left[\epsilon \left(\omega \right)\right]}^e$	cut	2.21, 2.35	561, 528		2.33, 2.39	533, 518
		(Chl a)2
PIDf	exp$\cdot$1+	1.90	652
$\frac{4\pi \omega }{c}\mathrm{Im}{\left[\alpha \left(\omega \right)\right]}^g$	full$\cdot$1+	2.05	605
$ℱ{\left\{d_m\left(t\right)\right\}}^c$	full$\cdot$1+	1.97	630
$\mathrm{Im}{\left[\epsilon \left(\omega \right)\right]}^d$	full$\cdot$1+	2.15, 2.32	578, 534
aPhotoinduced dissociation from Ref. 17. bTD-CAM-B3LYP/Def2-SVP from Ref. 17. cLCAO-TDDFT-$r$-$t$ from this work. dLCAO-TDDFT-$k$-$\omega$ from this work. ePW-TDDFT-$k$-$\omega$ from this work. fPhotoinduced dissociation from Ref. 19. gTD-CAM-B3LYP/Def2-SVP from Ref. 19.

Equations

$$\mathrm{Im}[\epsilon (\omega )]=\frac{1}{\Omega }\sum _{nm}\frac{4\pi \eta [f({\epsilon }_m)-f({\epsilon }_n)]}{(\omega -{\epsilon }_n+{\epsilon }_m-{\Delta }_x{)}^2+{\eta }^2}{\left|\frac{{\hat{\text{e}}}_{\text{q}}\cdot ⟨{\psi }_n|\nabla |{\psi }_m⟩}{{\epsilon }_n-{\epsilon }_m+{\Delta }_x}\right|}^2$$

$${\Delta }_x=\frac{8\sqrt{2}}{3{\pi }^2}\sum _{n=1}^N\left(\sqrt{{\epsilon }_{N+1}-{\epsilon }_n}-\sqrt{{\epsilon }_N-{\epsilon }_n}\right)⟨{\psi }_{N+1}|\frac{{\psi }_n^{\ast }{\psi }_n}{\rho }|{\psi }_{N+1}⟩,$$

$$⟨{\psi }_n|\nabla |{\psi }_m⟩=\sum _{\mu \nu }{c}_{\nu n}^{\ast }{c}_{\mu m}⟨{\tilde{\varphi }}_{\nu }|{\mathcal{T}}^{†}\nabla \mathcal{T}|{\tilde{\varphi }}_{\mu }⟩,$$

$$\mathcal{T}=1+\sum _{ai}\left(|{\phi }_i^a⟩-|{\tilde{\phi }}_i^a⟩\right)⟨{\tilde{p}}_i^a|,$$

$$\mathcal{F}\{d_m(t)\}=\frac{1}{\pi \mathcal{K}}\int \sin (\omega t)e^{-\frac{t^2}{2{\sigma }^2}}[d_m(t)-d_m(0)]dt$$

$$d_m(t)=-\iiint \rho (\text{r},t)(\text{r}\cdot {\hat{\text{e}}}_{\text{q}})dV.$$

$$i\frac{\partial }{\partial t}{\psi }_n(\text{r},t)={\mathcal{H}}_{\mathrm{KS}}(t){\psi }_n(\text{r},t),$$

$$\sum _{\mu }i⟨{\tilde{\varphi }}_{\nu }|{\mathcal{T}}^{†}\mathcal{T}|{\tilde{\varphi }}_{\mu }⟩\frac{\partial }{\partial t}{c}_{\mu n}(t)=\sum _{\mu ,m}⟨{\tilde{\varphi }}_{\nu }|{\mathcal{T}}^{†}{\mathcal{H}}_{\mathrm{KS}}^{nm}(t){\mathcal{T}}^{†}|{\tilde{\varphi }}_{\mu }⟩c_{\mu m}(t).$$

Full text

I Introduction

Chlorophyll a (Chl a or C55H72MgN4O5) and chlorophyll b (Chl b or C55H70MgN4O6) Moss (1988); Scheer (2006) are the fundamental functional units of the light harvesting complex (LHC II) Siefermann-Harms (1985) present in green plants. For this reason, understanding the photoexcitation process within Chl a and b is of great importance in the development of technologies such as those involved in the optimization of food crop production Long (2014); Ray et al. (2013) and conversion of solar radiation into a usable form of energy directly through methods such as conventional solar cells Ruban, Johnson, and Duffy (2011); Fleming et al. (2012) or via subsidiary technologies such as photosynthetically-driven (bio)reactor systems for hydrogen combustion Allakhverdiev et al. (2010); Krassen et al. (2009); Berardi et al. (2014); Searle and Malins (2014). Moreover, such information is more generally applicable to the in silico design and optimization of dye-sensitized solar cells Gong et al. (2017), organic photovoltaic cells Agnihotri (2014), photocatalytic systems Ahmad et al. (2016), optoelectronic devices Li (2009), and plasmonics Geddes (2016).

Although much progress has been recently made in both the experimental measurement of individual monomer and dimer Chl a and b spectra Milne et al. (2015); Stockett et al. (2015); Milne et al. (2016), and their theoretical description at the time-dependent density functional theory (TDDFT) Runge and Gross (1984) level Hedayatifar et al. (2016); Linnanto and Korppi-Tommola (2006); Dreuw and Head-Gordon (2004); Parusel and Grimme (2000), the lack of reasonably accurate yet highly efficient computational methods has hampered efforts to describe the optical absorption of large Chl-containing biomolecules. Initial attempts to investigate the optical absorption characteristics of biomacromolecules such as the LHC II using first-principles electronic structure methods have helped to clarify several aspects of the functioning of these systemsJornet-Somoza et al. (2015); Konig and Neugebauer (2011). However, these calculations neglected the protein matrix, and thus failed to fully explain the role played by the surrounding proteins in photosynthesis and why they are produced. Moreover, the computational resources required for a complete treatment of systems of this size lie considerably beyond what is generally available to most researchers.

On the one hand, methods based on the Kohn–Sham (KS) density of states Kohn and Sham (1965), while being quite efficient, often underestimate energy gaps by more than half. This is because an independent-particle picture fails to describe electronic screening of the excited states Migani et al. (2013). On the other hand, quasiparticle-based calculations of spectra from the Bethe–Salpeter equation (BSE) Salpeter and Bethe (1951), while often achieving quantitative accuracy, are extremely heavy computationally Onida, Reining, and Rubio (2002). As a result, only recently have even the smallest dye-sensitized solar cells (DSSC) been described at the BSE level Mowbray and Migani (2016). Moreover, such methods are intrinsically ill-suited to the description of isolated and/or non-periodic systems.

Although TDDFT Runge and Gross (1984); Marques et al. (2006) real time Yabana and Bertsch (1996, 1999); Ghosh et al. (2017); Koh, Nguyen-Beck, and Parkhill (2017) and frequency domain Casida (1995, 2009); Yan, Thygesen, and Jacobsen (2011); Yan et al. (2011); Mowbray (2014) calculations provide an attractive alternative, implementations based on real-space (RS) or plane-wave (PW) representations of the KS orbitals Kohn and Sham (1965) are both computationally expensive, and exhibit a strong exchange and correlation (xc) functional dependence for their accuracy. Time propagation RS (RS-TDDFT-$r$-$t$) calculations require time steps much shorter than what is needed to resolve the features of the spectra in order to ensure the stability of the calculation Jornet-Somoza et al. (2015). Such instabilities of RS calculations may be related to their freedom in representing the KS wavefunctions, which are only constrained by the grid spacing.

In this work we employ linear combinations of atomic orbitals (LCAOs) to provide a more efficient representation of the KS orbitals, while retaining the accuracy of PW-based TDDFT reciprocal space and frequency domain PW-TDDFT-$k$-$\omega$ calculations of the optical absorption. Moreover, we use the recently developed derivative discontinuity correction based on the exchange part of the GLLB-SC Kuisma et al. (2010) functional $\Delta_{x}$ to correct the KS eigenenergies and provide a semi-quantitative agreement with experimental measurements Castelli et al. (2012); Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019). Furthermore, the constraints imposed by an LCAO representation of the KS wavefunctions may be expected to improve the stability of time-propagation TDDFT calculations (LCAO-TDDFT-$r$-$t$), allowing one to use larger time steps. However, the reliability of LCAO-TDDFT is inherently basis set dependent Rossi et al. (2015). This makes an assessment and benchmarking for the fundamental functional units with PW-TDDFT-$k$-$\omega$ calculations and photoinduced dissociation (PID) experiments essential before applying LCAO-TDDFT-$k$-$\omega$ or LCAO-TDDFT-$r$-$t$ to the complete macromolecular system.

By applying LCAO-TDDFT methods Glanzmann et al. (2015); Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019); Kuisma et al. (2015) to the Chl a and b monomer and dimer systems, we may clearly explain their advantages and disadvantages for describing light-harvesting systems, with the aim of applying these methods to macromolecules such as the complete LHC II. Note that since experimental measurements of the Chl b dimer are currently unavailable, we have restricted consideration to the Chl a dimer herein.

This paper is organized as follows. In Sec. II we describe the computational parameters employed to model the Chl a and b systems and provide a theoretical background to LCAO-TDDFT-$k$-$\omega$ Glanzmann et al. (2015); Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019) and LCAO-TDDFT-$r$-$t$ Kuisma et al. (2015) calculations in the projector augmented wave (PAW) method. In Sec. III, we perform a basis set convergence test and direct comparison with PW-TDDFT-$k$-$\omega$ calculations and PID measurements for our LCAO-TDDFT-$k$-$\omega$ calculations of Chl a and b, model the excitonic density of their first four bright excitations $\omega_{i}$ using the electron hole density difference $\Delta\rho(\textbf{r},\omega_{i})$ from LCAO-TDDFT-$k$-$\omega$, compare our LCAO-TDDFT-$r$-$t$ calculations with PID measurements for Chl a and b, compare both LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ calculations with PID measurements for the Chl a dimer (Chl a)2, and finally tabulate and compare our results for the Q band with those available in the literature Milne et al. (2015, 2016). We then provide concluding remarks in Sec. IV. Atomic units ($\hslash=e=m_{e}=a_{0}=1$) are employed throughout unless stated otherwise.

II METHODOLOGY

All DFT Hohenberg and Kohn (1964); Kohn and Sham (1965) calculations were performed using the PAW Blöchl (1994); Blöchl, Först, and Schimpl (2003); Rostgaard (2006) implementation within the GPAW code Mortensen, Hansen, and Jacobsen (2005); Enkovaara et al. (2010); Larsen et al. (2009). We employed the generalized gradient approximation (PBE Perdew, Burke, and Ernzerhof (1996)) for the exchange and correlation functional. To represent the electron density and wave functions we employed either LCAOs Larsen et al. (2009) with a double-$\zeta$-polarized (DZP) basis set, after performing convergence tests with basis sets of varying quality up to a quadruple-$\zeta$-polarized (QZP) basis set, or PWs with a plane-wave cutoff of 340 eV. Gas phase structures for the full Chl a and b monomers (see Figure 1(a) full$\cdot$1+) and the Chl a dimer (see Figure 2) were based on those of Refs. 17 and 19, respectively. To determine the influence of the hydrocarbon chain and the tetramethylammonium charge tag (N(CH3)${}_{4}^{+}$ or 1+) on the optical absorbance of Chl a and b monomers we have also considered structures with the hydrocarbon chain cut from C20H39 to C5H9 with and without the N(CH3)${}_{4}^{+}$ charge tag, as shown in Figure 1 (a) cut$\cdot$1+ and (b) cut, respectively.

Each structure was optimized using the ASE code Bahn and Jacobsen (2002); Larsen et al. (2017) until a maximum force less than 0.03 eV/Å was obtained. All structures were modelled in supercells with more than 6 Å of vacuum using non-periodic boundary conditions, that is, the electron density and wave functions were set to zero at the cell boundary. This prevents long-range interactions between repeated images that would affect a periodic calculation. More importantly, the use of non-periodic boundary conditions has been previously shown to be essential for modelling charged structures in gas phase Mowbray et al. (2013).

Optical absorbance spectra were modelled using the imaginary part of the dielectric function, $\operatorname{Im}[\varepsilon(\omega)]$, from LCAO-TDDFT calculations in the frequency domain (LCAO-TDDFT-$k$-$\omega$) Glanzmann et al. (2015); Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019); LCA , and using the Fourier transform of the dipole moment, $\mathcal{F}\{d_{m}(t)\}$, from LCAO-TDDFT in the real time domain (LCAO-TDDFT-$r$-$t$) Kuisma et al. (2015).

The optical absorbance from LCAO-TDDFT-$k$-$\omega$ in the $\hat{\textbf{e}}_{\textbf{q}}$ direction, neglecting local field effects, is given by Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019)

[TABLE]

where $\Omega$ is the supercell volume, $\eta$ is the electronic broadening, $f$ is the Fermi–Dirac function, $\varepsilon_{n}$ is the eigenenergy and $\psi_{n}$ is the KS wave function of the $n$th orbital, and $\Delta_{x}$ is the derivative discontinuity correction to the eigenenergies Kuisma et al. (2010) from the exchange part of the GLLB-SC functional given by Preciado-Rivas, Torres-Sánchez, and Mowbray (2019)

[TABLE]

where $N$ is the number of electrons. Incorporating this derivative discontinuity correction has been shown to provide better agreement with experimental band gaps Castelli et al. (2012) and optical absorption spectra Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019).

The matrix elements in Eq. (1) may be expressed as Glanzmann et al. (2015)

[TABLE]

where the sum is over $\mu\equiv\{i,a\}$ for the $i$th state centered on atom $a$, $\tilde{\phi}_{\mu}$ are the localized basis functions, so that $|\tilde{\psi}_{m}\rangle=\sum_{\mu}c_{\mu m}|\tilde{\phi}_{\mu}\rangle$ with coefficients $c_{\mu m}$ for the $m$th KS wave function, and $\mathcal{T}$ is the PAW transformation operator Blöchl (1994); Blöchl, Först, and Schimpl (2003); Rostgaard (2006)

[TABLE]

where $\tilde{\varphi}_{i}^{a}$ and $\varphi_{i}^{a}$ are the pseudo and all-electron partial waves for state $i$ on atom $a$ within the PAW formalism, and $|\tilde{p}_{i}^{a}\rangle$ are the smooth PAW projector functions.

Since the matrix elements of Eq. (3) must already be calculated to obtain the forces at the DFT level, the calculation of $\operatorname{Im}[\varepsilon(\omega)]$ using Eq. (1) simply involves the multiplication of previously calculated matrices. For this reason, the calculation of $\operatorname{Im}[\varepsilon(\omega)]$ using LCAO-TDDFT-$k$-$\omega$ is very efficient, with scaling of $\mathcal{O}(NM^{2})$ or better Blackford et al. (1997), where $N$ is the number of KS wavefunctions and $M\geq N$ is the total number of basis functions used in the LCAO calculation. This provides a significant speed-up compared to the $\mathcal{O}(N^{5})$ scaling of the Casida LCAO-TDDFT-$r$-$\omega$ formalism Kuisma et al. (2015); Casida (1995, 2009). We found ten unoccupied KS wave functions per chlorophyll molecule were already sufficient to converge the spectra up to 3.5 eV. This is because local field effects may be safely neglected in the optical limit for gas phase structures Mowbray (2014); Glanzmann et al. (2015).

The optical absorbance from LCAO-TDDFT-$r$-$t$ is given by the Fourier transform of the dipole moment, $d_{m}$, so that

[TABLE]

where $\mathcal{K}=1\times 10^{-5}$ is the strength of the “kick” in the $\hat{\textbf{e}}_{\textbf{q}}$ direction, and $\sigma$ is the peak broadening. The dipole moment is then obtained from the time dependent charge density $\rho(\textbf{r},t)$ using

[TABLE]

The time dependent charge density is simply $\rho(\textbf{r},t)=-e\sum_{n}f(\varepsilon_{n})|\psi_{n}(\textbf{r},t)|^{2}$, where $\psi_{n}(\textbf{r},t)$ are the time-dependent KS wave functions which are obtained by propagating

[TABLE]

where $\mathcal{H}_{\mathrm{KS}}(t)=-\frac{1}{2}\nabla^{2}+v_{\mathrm{KS}}[\rho(\textbf{r},t)](\textbf{r},t)$ is the KS Hamiltonian for independent electrons, and $v_{\mathrm{KS}}$ is the KS potential functional.

Employing LCAO to represent the time-dependent KS wave functions, we may express Eq. (7) as Kuisma et al. (2015)

[TABLE]

After applying the initial kick, the time propagation is performed using the reliable and numerically stable semi-implicit Crank–Nicolson method Crank and Nicolson (1947), as described in Ref. 48.

As was the case for LCAO-TDDFT-$k$-$\omega$, each step of the LCAO-TDDFT-$r$-$t$ is very efficient and involves solving the system of linear equations, Eq. (8), with ScaLAPACK Blackford et al. (1997). For this reason, the calculation of $\mathcal{F}\{d_{m}(t)\}$ using LCAO-TDDFT-$r$-$t$ scales as $\mathcal{O}(NM^{2})$. Although this is worse than the $\mathcal{O}(NG)$ scaling of RS-TDDFT-$r$-$t$, where $G$ is the number of grid points, the constant prefactor for the grid propagation is so large that even for systems with thousands of electrons it is outperformed by several orders of magnitude by LCAO-TDDFT-$r$-$t$ Kuisma et al. (2015).

All spectra have been calculated with intrinsic widths of 50 meV. In the frequency domain we employed an electronic broadening of $\eta=50$ meV to the individual Lorentzian peaks, while in the real time domain we convoluted the dipole moment with a Gaussian of width $\sigma\approx 13$ fs.

The LCAO-TDDFT-$r$-$t$ spectra remained stable with a 0.01 fs time step, which was sufficient to resolve the shape of the dipole spectra. Such a time step is an order of magnitude greater than the 0.002 fs time step required for real space time propagation calculations to maintain stability Jornet-Somoza et al. (2015). We were thus easily able to propagate up to 80 fs, twice that recorded with real space methods Jornet-Somoza et al. (2015).

To provide insight into the spatial distribution of the exciton, we employ the electron hole density difference, $\Delta\rho(\textbf{r},\omega)=\rho_{h}(\textbf{r},\omega)+\rho_{e}(\textbf{r},\omega)$, as described in Refs. 31; 44; 45. Here, the electron/hole densities $\rho_{e/h}(\textbf{r})$ are obtained by averaging the two-particle exciton density $\rho_{ex}(\textbf{r}_{e},\textbf{r}_{h})$, with respect to the hole/electron coordinate $\textbf{r}_{h/e}$. We calculate the electron hole density difference $\Delta\rho(\textbf{r},\omega)$ for the first four bright excitations using the transitions from LCAO-TDDFT-$k$-$\omega$. Since these molecular transitions are both well separated and composed of single transitions $m\to n$ at the LCAO-TDDFT-$k$-$\omega$ level, we find $\Delta\rho(\textbf{r},\omega_{mn})\approx e(|\psi_{n}|^{2}-|\psi_{m}|^{2})$, where $\psi_{m}$ and $\psi_{n}$ are the KS orbitals corresponding to the hole and electron, respectively, and $\omega_{mn}\approx\varepsilon_{n}-\varepsilon_{m}+\Delta_{x}$ is the energy of the transition $m\to n$, including the derivative discontinuity correction $\Delta_{x}$ of Eq. (2).

III RESULTS AND DISCUSSION

In Figure 3 we show the optical absorption spectrum of the neutral cut structure of Chl a calculated using either multiple-$\zeta$ basis sets (SZ, DZ, TZ and QZ) or the polarized versions of the basis sets (SZP, DZP, TZP and QZP) within the LCAO mode of the gpaw code. We systematically increase the number of functions to assess the sensitivity of the optical absorbance and observe that the spectrum converges differently for basis sets with and without polarization functions.

It can be seen in Figure 3(a) that the SZ basis set yields a red-shifted spectrum compared to the other spectra, the TZ and QZ basis sets spectra are almost the same, and the DZ basis set yields the same energies of the Q and Soret bands compared to those of the TZ and QZ but the peaks within the band are more separated than in the TZ or QZ spectra. The polarized multiple-$\zeta$ basis sets are shown in Figure 3(b) and we observe that the SZP basis set also yields a red-shifted spectra compared to other spectra and there are more peaks above $3$ eV. The spectra obtained using TZP and QZP are nearly identical and the DZP basis set yields the same energy of the Q and the Soret bands although the peaks within the bands are a little bit separated when compared to TZP and QZP.

Even though it seems that the spectrum of Chl a will not change if we increase the number of functions of the basis set, we cannot rely solely on these results to claim that an LCAO basis set will guarantee an accurate description of the optical absorption spectra. This is because LCAO basis sets cannot be systematically converged to the complete basis set limit, for molecules of this size. This is clear from comparing the result of the inclusion of polarization functions with the LCAO basis sets in Figure 3. Although double-$\zeta$ basis sets are converged with respect to the number of radial functions, the inclusion of polarization systematically alters the spectra. For this reason, we must compare the spectra of Chl a and Chl b obtained using LCAO and a PW representations of the KS wavefunctions.

In Figure 4, we show the optical absorption spectra of the cut structures of both Chl a and Chl b calculated using either the DZP basis set within the LCAO mode, or using a PW representation of the KS wavefunctions that ensures converged optical absorption spectra. We find that the DZP basis set yields an optical absorption spectrum with transition energies very close to those of the spectrum calculated using PW. We observe that the LCAO spectra is only slightly (much smaller than the $0.1$ eV DFT accuracy) red-shifted compared to the PW spectra, but underestimates the intensities of the transitions. Nevertheless, we can affirm that the atomic basis set with double-$\zeta$ and polarization functions, DZP, is sufficient to ensure semi-quantitative agreement of the optical absorption spectra for the Chl a and Chl b with the PW calculations.

The DZP basis sets have been shown to be sufficient (and often necessary) to converge to the results calculated with PW representations of the KS wavefunctions in other works as well Strange et al. (2008); Larsen et al. (2009). From hereon we shall restrict consideration to the DZP basis set and simply refer to our calculations as LCAO.

In Figure 5 we plot the spectra obtained from LCAO-TDDFT-$k$-$\omega$ including the derivative discontinuity correction $\Delta_{x}$ for the Chl a and b monomers shown in Figure 1. The onset of the calculated spectra are in semi-quantitative agreement with the experimental spectra obtained with either the monocationic tetramethylammonium 1+ or acetylcholine 3+ tag. The difference between the maximum of the Q band and the calculated first excitation energy is less than $0.2$ eV and the difference between the maximum of the Soret band and the calculated fourth excitation is less than $0.1$ eV. Also, the spectra of the *cut$\cdot$*1+ structure shows the same qualitative behavior as the *full$\cdot$*1+ structure, suggesting that none to negligible changes in the optical absorption spectra are caused by the carbon chain. However, when the charge tag is removed, that is, regarding the cut structure, the Q band peaks are blue shifted and the band gap is widened, whereas the Soret band peaks are red shifted.

For Chl b, Figure 5(b), the calculated spectra for the full$\cdot$1+ structure is also in semi-quantitative agreement with the experimental spectra. Specifically, the peaks of the Q band and the Soret band of the *full$\cdot$*1+ structure are blue and red shifted, respectively, when compared to the experimental data, with a difference between the Q band maximum and the calculated first excitation of $0.23$ eV and a difference between the Soret band maximum and the calculated fourth excitation of $0.26$ eV. Again, the spectrum of the *cut$\cdot$*1+ structure is qualitatively the same as the spectrum of the *full$\cdot$*1+, reinforcing the idea that the carbon chain has no impact in the optical absorbance and that it should be centered on the Mg atom and the chlorin ring. When removing the charge tag (cut structure) the Q band peaks are again blue shifted, and the band gap is widened, whereas the Soret band peaks are red shifted.

The intensities of the excitations of the Soret band are being underestimated about $75\%$ for Chl a and $52\%$ for Chl b, comparing the *full$\cdot$*1+ to the experimental data. Also, the intensity of the Soret band peaks is less compared to that of the Q band peaks. This is not the case for the experimental data. Although the relative intensities are adjusted, that is, the Soret band peaks intensities are increased by applying the GLLB-SC correction to the spectra following Eq. 1, they are still underestimated.

Such a difference can be explained by the fact that we are neglecting charge transfer excitations at the linear density response level. The higher intensity of these peaks is often attributed in the literature to this type of charge transfer excitation Milne et al. (2015), which we are unable to model at the LCAO-TDDFT-$k$-$\omega$ level. In this way, our underestimation of the Soret band intensity provides indirect insight into the nature of the experimentally observed peaks.

In Figure 6 we plot the spatial distribution of the electron and hole densities for the first four bright transitions of the Chl a and b cut structures. As shown in Figure 5, the first and fourth excitations are induced by optical absorption along the N–Mg–N bond in the $x$ polarization direction Moss (1988); Scheer (2006) (blue arrows in Figure 6), while the second and third excitations are induced by optical absorption along the $y$ polarization direction Moss (1988); Scheer (2006) (red arrows in Figure 6), as expected. In each case, the excitations are $\pi\rightarrow\pi$ transitions involving the two highest occupied molecular orbitals (HOMOs) and the two lowest unoccupied molecular orbitals (LUMOs).

For both species the electron density of the first and third excitations has weight on the $2p_{z}$ levels of the N atoms parallel to the direction of excitation towards the Mg atom. Conversely, the hole density of the second and third excitations on both Chl a and b has significant weight on the $2p_{z}$ levels of both the N and Mg atoms. Overall, the first and third transitions tend to move charge density from the chlorin ring to either its edges or towards the central Mg atom, whereas the second and fourth transitions tend to move charge from the central Mg atom to the edge of the chlorin ring along the direction of excitation.

Although the excitations for Chl a and b are qualitatively the same, we do notice differences in the region of the methyl/aldehyde. Specifically, the electron of the third and fourth transitions has significantly more weight around the aldehyde group of Chl b compared to the methyl group of Chl a. Also, the hole density of the second and fourth transitions has weight on all four N atoms for Chl b, whereas the electron density has some weight on one of the N atoms perpendicular to the direction of excitation.

Finally, the excitations of the charge tagged species with a cut (cut$\cdot$1+) and full (full$\cdot$1+) hydrocarbon chain, shown in Figure 7, are basically the same as those for the neutral species (cut shown in Figure 6). This provides additional evidence justifying the experimental use of charge tagged Chl a and b molecules to describe the optical absorption of the neutral isolated species Milne et al. (2015); Stockett et al. (2015).

In Figure 8 we directly compare the spectra obtained from LCAO-TDDFT-$r$-$t$ for the Chl a and b monomers shown in Figure 1 with those measured experimentally from Refs. 17 and 18. The Q-band energy or onset of the LCAO-TDDFT-$r$-$t$ spectra agrees quantitatively with that measured experimentally for both Chl a and b in Refs. 17 and 18. However, whereas we found distinct Q band transitions with polarizations in the $x$ and $y$ directions for LCAO-TDDFT-$k$-$\omega$ (see Figure 5), with LCAO-TDDFT-$r$-$t$ we find the Q band transitions coincide in energy. Moreover, for Chl b cut$\cdot$1+, we find that whereas there are two Q-band peaks, they both have polarization in the $x$ direction, from ring B to ring D along the N–Mg–N bond. This splitting of the Q band peak in the $x$ direction and the weak Q band peak in the $y$ direction are both suppressed once the charge tag is removed (see Figure 8(b) cut).

While we find removing the charge tag induces a blue shift of the Q band for Chl a and b, as was the case for LCAO-TDDFT-$k$-$\omega$, we find the Soret band is less affected by the charge tag with LCAO-TDDFT-$r$-$t$. In the region of the Soret band LCAO-TDDFT-$r$-$t$ yields many peaks from 2.5–4 eV. This may be related to the difficulty in converging higher energy peaks with TDDFT-$t$, although our time propagation of $\sim 80$ fs should be sufficient to converge these peaks.

Overall, we find the spectra from LCAO-TDDFT-$r$-$t$ provide excellent agreement with the experimentally measured spectra with tetramethylammonium (exp$\cdot$1+) and acetylcholine (exp$\cdot$3+) charge tags. This provides further evidence that the exact structure of the charge tag has little to no impact on the onset of the spectra Milne et al. (2015); Stockett et al. (2015). Moreover, the spectra exhibit rather minor changes upon removing the charge tag, providing further justification for using charge tagged Chl a and b to describe the optical absorption of the neutral isolated species Milne et al. (2015); Stockett et al. (2015).

Going beyond the monomer, in Figure 9 we directly compare the optical absorbance spectra for the Chl a dimer measured experimentally from Ref. 19, obtained from LCAO-TDDFT-$k$-$\omega$, $\operatorname{Im}[\varepsilon(\omega)]$, and LCAO-TDDFT-$r$-$t$, $\mathcal{F}\{d_{m}(t)\}$. In both cases, the experimental onset of the peak is reproduced semi-quantitatively, with the derivative discontinuity corrected LCAO-TDDFT-$k$-$\omega$ results somewhat blue shifted by about $0.2$ eV. However, as was the case for the monomer, the experimentally observed splitting of the Q band is clearly visible with LCAO-TDDFT-$k$-$\omega$ but absent from the LCAO-TDDFT-$r$-$t$ Chl a spectrum.

Overall, these results show that increasing the system size, as we go from the monomer to the dimer, LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ provide the same level of accuracy. For LCAO-TDDFT-$k$-$\omega$, this means that after applying the derivative discontinuity correction $\Delta_{x}\approx 0.5$ eV, the spectra agree semi-quantitatively with the measured spectra for both the Chl a monomer and dimer. For LCAO-TDDFT-$r$-$t$, this correction is not required to obtain a quantitative description of the experimental spectra for both the Chl a monomer and dimer. Altogether, this suggests that for even larger systems, such as the LHC II, we may reasonably expect both LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ to continue to perform well.

In Table 1 we provide a direct comparison of the Q-band energy (first and second transitions shown in Figures 6 and 7) for Chl a, Chl b, and (Chl a)2 obtained from PID measurements, TD-CAM-B3LYP/Def2-SVP calculations of the Casida LCAO-TDDFT-$r$-$\omega$ type Casida (1995, 2009) using the imaginary part of the dynamic polarizability $\alpha(\omega)$ ($\frac{4\pi\omega}{c}\operatorname{Im}[\alpha(\omega)]$), LCAO-TDDFT-$r$-$t$ ($\mathcal{F}\{d_{m}(t)\}$), PW-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$k$-$\omega$ ($\operatorname{Im}[\varepsilon(\omega)]$). Overall, we find excellent agreement between PID, LCAO-TDDFT-$r$-$t$, and TD-CAM-B3LYP/Def2-SVP. In fact, our LCAO-TDDFT-$r$-$t$ calculations perform somewhat better than TD-CAM-B3LYP/Def2-SVP for describing the Q band energies of Chl a, Chl b, and (Chl a)2. These LCAO-TDDFT results, along with its computational speed-up, provide strong motivation for the future use of both LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ for describing larger systems such as the LHC II.

IV CONCLUSIONS

We have shown that LCAO-TDDFT-$k$-$\omega$ Glanzmann et al. (2015); Preciado-Rivas, Torres-Sánchez, and Mowbray (2019); Lyon et al. (2019) and LCAO-TDDFT-$r$-$t$ Kuisma et al. (2015) calculations are both significantly faster than their PW-TDDFT-$k$-$\omega$, TD-CAM-B3LYP/Def2-SVP, and RS-TDDFT-$r$-$t$ counterparts and also provide a similar level of accuracy to alternative TDDFT methods in the description of the optical absorption spectrum of the key photosynthetic light-harvesting molecules Chl a and b and the dimeric system (Chl a)2. Specifically, our LCAO-TDDFT-$k$-$\omega$ and LCAO-TDDFT-$r$-$t$ calculations provide both a semi-quantitative and qualitative description of the spectra of chlorophyll monomers and dimers, with the former yielding a significant reduction in memory requirements compared to PW-TDDFT-$k$-$\omega$ and the latter yielding an order of magnitude speed-up compared to RS-TDDFT-$r$-$t$ calculations due to the increased numerical stability of LCAO-TDDFT-$r$-$t$. These results pave the way for future studies applying these highly efficient and accurate LCAO-TDDFT methods to the characterization of biomacromolecules, such as the LHC II, and the full understanding of the physical and chemical processes involved in photosynthesis.

Acknowledgements.

This work employed the Imbabura cluster of Yachay Tech University, which was purchased under contract No. 2017-024 (SIE-UITEY-007-2017), and the Quinde I supercomputer of Public Company Yachay E. P., which was implemented under contract No. 0051-2015, corresponding to Component No. 7 of Group No. 2, Re-YACHAY-018-2015. D.J.M. thanks the management and operation team of the Quinde I Supercomputer for their help and assistance in this project. A.H.L. acknowledges funding from the European Union’s Horizon 2020 research and innovation program under grant agreement no. 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence; the European Research Council (ERC-2010-AdG-267374); Spanish grant (FIS2013-46159-C3-1-P); and Grupos Consolidados (IT578-13). The Coimbra Chemistry Centre is supported by FCT, through the Project PEst-OE/QUI/UI0313/2014 and POCI-01-0145-FEDER-007630. B.F.M. thanks the Portuguese Foundation for Science and Technology (projects CENTRO-01-0145-FEDER-000014 and POCI-01-0145-FEDER-032229) and the Donostia International Physics Centre for financial support.