Towards the description of charge transfer states in solubilised LHCII using subsystem DFT

Light harvesting complex II (LHCII) in plants and green algae have been shown to adapt their absorption properties, depending on the concentration of sunlight, switching between a light harvesting and a non-harvesting or quenched state. In a recent work, combining classical molecular dynamics (MD) simulations with quantum chemical calculations (Liguori et al. in Sci Rep 5:15661, 2015) on LHCII, it was shown that the Chl611–Chl612 cluster of the terminal emitter domain can play an important role in modifying the spectral properties of the complex. In that work the importance of charge transfer (CT) effects was highlighted, in re-shaping the absorption intensity of the chlorophyll dimer. Here in this work, we investigate the combined effect of the local excited (LE) and CT states in shaping the energy landscape of the chlorophyll dimer. Using subsystem Density Functional Theory over the classical \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}μs MD trajectory we look explicitly into the excitation energies of the LE and the CT states of the dimer and their corresponding couplings. Upon doing so, we observe a drop in the excitation energies of the CT states, accompanied by an increase in the couplings between the LE/LE and the LE/CT states facilitated by a shorter interchromophoric distance upon equilibration. Both these changes in conjunction, effectively produces a red-shift of the low-lying mixed exciton/CT states of the supramolecular chromophore pair. Supplementary Information The online version contains supplementary material available at 10.1007/s11120-022-00950-7.


Introduction
Light harvesting complex II (LHCII) or the major antenna complex of the Photosystem-II (PSII) in higher plants and green algae, is a complex structure accounting for almost 80% of the light absorbed by PSII. LHCII exists in the native trimeric form where one monomer consists of eight chlorophyll a (Chla), six chloropyll b (Chlb), two luteins (Lut), one neoxanthin (Nx) and one violaxanthin (Vx) chromophore as cofactors. The trimeric (major) LHCII and the monomeric (minor) LHCII complexes (CP24, CP26 and CP29) along with the PSII core form the PSII supercomplex, the crystal structure of which has been determined at high resolution (Wei et al. 2016;Su et al. 2017). Previous studies on LHCs, have highlighted the importance of the excitation energy transfer (EET) processes in fine tuning the spectral properties of the complexes (Segatta et al. 2019;Cupellini et al. 2020b;Balevičius et al. 2017;Maity et al. 2019;Duffy et al. 2013;Chmeliov et al. 2015;Lapillo et al. 2020;Khokhlov and Belov 2019). [We refer the reader to reference (Segatta et al. 2019) for a focused review on the different theoretical methods and to reference (Cupellini et al. 2020b) for the challenges therein]. In addition Electron Transfer (ET) effects, giving rise to Charge Transfer (CT) states, have also been shown to play a crucial role in modulating the overall spectral and energy landscape of the complexes especially for closely interacting pigments (Kell et al. 2014;Wahadoszamen et al. 2014;Novoderezhkin et al. 2007;Miloslavina et al. 2008;Cupellini et al. 2020a). Although CT states by themselves are optically dark states, previous theoretical studies have shown that these states can mix with the local states thereby gaining some dipole strength Nottoli et al. 2018). Previous studies have also Souloke Sen and Lucas Visscher have contributed equally to this work. elucidated a strong correlation between the LHC energetics and protein conformational changes (Ruban et al. 2012;Krüger 2014;Liguori et al. 2013), showing that the so-called 'quenched' states of LHCs are linked to particular conformations where the chromophore-chromophore interactions are strengthened (Pascal et al. 2005;Ruban et al. 2007;Ahn et al. 2008;Bode et al. 2009;Miloslavina et al. 2008). In this respect, interactions between a particular chlorophyll dimer (originally labelled Chla611-Chla612 in reference Liu et al. 2004) is particularly interesting because it is often assumed to be a critical site for the switch between light harvesting and quenched conformations of LHCs (Ruban et al. 2007;Mozzo et al. 2008;Mascoli et al. 2019).
The advent of Classical Molecular Dynamics (MD) simulations in the microsecond ( s ) time scale allows for probing the chromophore dynamics at atomistic resolution, also allowing the inclusion of lipid membrane and solvent effects (Curutchet and Mennucci 2017). Snapshots from these long simulations can be used as input for quantum chemical calculations to provide a computationally efficient strategy for studying spectral properties of LHCs over such time scales. We note that care has to be taken that the force field reproduces the structure of the chromophores well enough (Andreussi et al. 2017;Jurinovich et al. 2015) and have checked in earlier work ) that the difference of excitation energies between MM and QM optimized structures does not exceed 0.1 eV. In that work, we performed Time-Dependent Density Functional Theory (TDDFT) (Runge and Gross 1984;Casida 1995) calculations on the above chlorophyll dimer (i.e. Chla611-Chla612) on selected snapshots from a particular trajectory obtained from a s classical MD simulation which showed significant structural disorder from the crystal structure (Liguori et al. 2015). This particular trajectory [trajectory "A" in the original work of reference (Liguori et al. 2015) that was also used in reference (López-Tarifa et al. 2017)] was one of six independent simulations where the excitonic interactions between the chlorophyll dimer displayed marked differences compared to the others (Liguori et al. 2015). In our previous work, it was shown that a drop in the total oscillator strength of the Q bands ( Q y and Q x ) could be explained using a transition density decomposition analysis, wherein an increase in the CT component of the oscillator strength was observed upon equilibration. Furthermore a characterization of the excitonic states using charge transfer descriptors of Plasser and Lischka (2012) revealed an overall increase of the CT character of the low-lying states. Both of these analysis clearly showed an increase in the mixing of the CT states with the local excited (LE) Q states, making it difficult to make a quantitative calculation of the energies of the LE and the CT states towards the end of the trajectory.
In this previous work, the delocalized orbital picture resulting from the use of a supramolecular basis made the interpretation of the TDDFT states difficult when the distance between the chromophores was small. In the current work, we instead use a diabatic representation built on the framework of subsystem DFT in which the identification as LE or CT states becomes straightforward. A pictorial description of the involved electronic states for both these processes for two identical pigments is shown in Fig. 1. We construct a model Hamiltonian that can be used to predict the importance of the ET and EET processes by computing the subsystem DFT energies of the individual diabatic states (diagonal elements of the model Hamiltonian) as well as their couplings (the off-diagonal elements). We thereby focus on the pair Chla611-Chla612 and compare computed matrix elements from the beginning and the end of the MD trajectory. Upon doing so, we observe a decrease in the energies of the CT states which are accompanied by an increase in the coupling between the LE and CT states and among the LE states themselves, resulting from a shorter interchromophoric distance between Chla611 and Chla612 after equilibration. We show that this conformational change results in a red-shift of the low-lying mixed excitonic/CT state, induced by a drop in the CT state energy and an increased coupling between the local and CT state and between the local states themselves.
This paper is organized as follows: In the "Methods" section, we introduce the model Hamiltonian and give an overview of the subsystem DFT framework. Therein, the defining equations for the approximations used in the calculation Fig. 1 Excitation energy transfer and electron transfer are represented in the monomer coordinate system. The red arrows represents Förster (Förster 1948) energy transfer flow, whereas the green arrows represent the Marcus (Marcus and Sutin 1985) electron transfer between pigments A and B. Vertical excitation processes are represented with broken arrows. The ground, local excited and ionized potential surfaces of the corresponding pigments A and B are shown in black, red and green (denoted by Ψ A∕B i ) respectively in the increasing order of energy.
A∕B i denotes the reorganization energy associated with the electronic transition 1 3 of the matrix elements of the effective Hamiltonian (representing the local and CT energies and their couplings) are discussed. In the "Computational details" section, the technical details of the workflow and the calculations performed to obtain the matrix elements are specified. In the "Results and discussion" section, we then discuss the obtained values of the matrix elements of the model Hamiltonian and compare with the supramolecular picture. In the final concluding section we also present a brief outlook for further studies.

Methods
In order to obtain approximate solutions of the Schrödingers equation of a many-body system involving electrons and nuclei one usually employs the Born-Oppenheimer approximation (1927). In this approximation the total wavefunction of the system, Ψ i ( , ) , where and represent the set of coordinates of the nuclei and the electrons, respectively, is written as a product of the electron and nuclear wavefunctions, i.e.
where Φ i is ith electronic wave function obtained for a fixed nuclear configuration , and the i are nuclear wave functions for vibrational states . The product of these two wave functions, Φ i ( , ) i ( ) , describes so-called vibronic states i and their energy, E i . Vibronic transitions are associated with a change in both the electronic and the vibrational states, the corresponding rate of which is given by the Fermi's golden rule, where k i j is the rate, Ĥ 1 a perturbation operator and E i , E j are the energies of the initial and final vibronic states i and j respectively. In the Franck-Condon approximation, the matrix element, H 1 i j is separated into two factors, i.e., The first term on the rhs of Eq. 4 represents the coupling matrix element between the electronic states i and j that is our main interest. The nuclear wavefunction overlap matrix elements, given by the second term on the rhs of Eq. 4, determine the homogenous (vibrational) broadening of the transition between the two electronic states i and j and will not be considered explicitly. We aim to calculate the matrix elements of the electronic Hamiltonian Ĥ el in a diabatic basis, consisting of the electronic  Fig. 2 The ground state, Φ GS , two locally excited states, Φ A * B and Φ AB * , and two CT states, Φ CT1 , Φ CT2 along with the distribution of electrons, where each fragment is considered as a two-level two electron system, are shown. a † and a are the corresponding creation and annihilation operators, and i, j and a, b represent occupied and unoccupied orbitals of the fragments A and B, respectively 1 3 and in the above definitions for convenience. In this diabatic basis the Hamiltonian matrix assumes a blocked structure with diagonal blocks representing the electronic ground state ( el GS ), interactions between LE states, ( el LE ), as well as interactions between CT states ( el CT ). Off-diagonal blocks ( el X−Y , where X, Y ∈ (GS, LE, CT) and X ≠ Y ) couple these physically distinct states and lead to states with mixed LE-CT character.
In order to describe the final states we need to compute the values of the matrix elements in each of the blocks. We can set the energy of the ground state to zero and define the diagonal matrix elements of el LE and el CT block, i.e. ⟨Φ i �Ĥ el �Φ i ⟩ , where, i ∈ (A * B, AB * , CT1, CT2) as LE and CT excitation energies. The calculation of off-diagonal matrix elements inside the el LE block, i.e. ⟨Φ i �Ĥ el �Φ j ⟩, where i, j ∈ (A * B, AB * ) and i ≠ j is facilitated by the fact that the two transition densities are (almost) non-overlapping. We furthermore note that off-diagonal elements in el CT block, i.e. the coupling between CT states in which the transfer goes in opposite directions can be assumed negligible as they correspond to a two-electron process. The couplings between the GS and CT state and between the GS and LE states are accounted for in the calculation of the CT and LE energies and need not be considered explicitly. The final significant coupling is between LE and CT states and will be calculated in a generalization of Fragment Orbital DFT (Senthilkumar et al. 2003(Senthilkumar et al. , 2005) that we describe below.
As the diabatic states will be constructed independently from each other, they are in general non-orthogonal and adiabatic eigen states Ψ ad and energies E ad need therefore be determined by solving a generalized eigenvalue problem where dia denotes the diabatic matrix representation of el . The overlap matrix dia is constructed in addition to evaluation of the matrix elements of dia and ad are the coefficients of the adiabatic eigenstates in the diabatic basis. As the coupling between the ground state and the locally excited (10) dia ad = E ad dia ad , states is zero by construction, and the coupling between the ground state and the charge transfer states is very weak for the systems we are considering, we thereby focus only on the set of diabatic states { Φ A * B , Φ AB * , Φ CT1 Φ CT2 }. The structure of the matrices dia and dia used in this work is shown in Fig. 3. We note that dia shares a similar structure with other ab-initio exciton models used to describe such processes Nottoli et al. 2018). In the next sections we will introduce a composite approach to obtain the non-zero matrix elements.

Subsystem DFT
While the theory has thus far been general and applicable with any diabatic wave function ansatz, we will now specialize to the wave functions provided by subsystem DFT (Jacob and Neugebauer 2014;Wesolowski et al. 2015). Subsystem DFT is typically formulated such that a monomer basis set is used to represent a subsystem and will thereby by construction yield local states. We, therefore construct each of the aforementioned diabatic states from separate sets of (spin) orbitals of A and B resulting from an underlying subsystem DFT calculation. The ground state determinant, Φ GS in Eq. 5 ( Fig. 2) of the combined system A + B in the subsystem framework is thereby written as combination of two sets of subsystem spinorbitals In the following section we briefly outline the basic formulation of subsystem DFT and thereafter introduce the different flavours that are used in approximating the matrix elements of dia . The basic idea of subsystem DFT is to partition the electron density of a total system into subsystem contributions (Wesolowski and Warshel 1993) (11) where n is the number of subsystems. The total energy of the system is then minimized with respect to all K , where each of the subsystem densities is usually represented by a Kohn-Sham like system of independent particles. In case of two subsystems A and B, the ground state energy of A + B is then given by, where Therefore in subsystem DFT, in addition to the usual E xc functional, we also need to approximate   where and s are the spatial and spin coordinates, respectively. For simplicity, we assume that both subsystems have an even number of electrons so that we can integrate out the spin coordinate. The one-electron Kohn-Sham potentials are then and (12) in which a non-additive exchange-correlation energy is defined for convenience: The above procedure can be made fully self-consistent by including the mutual polarization between the subsystems by the effect of the embedding potential v emb (r) in so-called freeze-and-thaw cycles Gritsenko 2013;Jacob and Neugebauer 2014). Compared to other embedding schemes one may in this way account for polarization while still allowing for the Pauli repulsion that is missing in classical embedding (Jacob et al. 2006). In applying these Kohn-Sham wave functions, there are two sources of non-orthogonality (Pavanello et al. 2013), the first one stemming from the overlap between the full electron diabatic states as in Eq. 10, and the second one arising at the subsystem level from the fact that the two set of orbitals } of subsystem A and B are non-orthogonal with respect to each other.

Calculation of dia LE block
Having discussed the calculation of the ground state, Φ GS , by subsystem DFT in the previous section, we now introduce the formulas used in obtaining the matrix elements of the blocks dia LE , dia CT and dia LE−CT that we need in this work. We will compute these blocks independently (each block is computed without the knowledge of the other one) in separate sets of calculations. We assume that all of the LE states Φ A * B and Φ AB * in Eqs. 6 and 7 ( Fig. 2) can formally be constructed from single-particle excitations from either A 0 or B 0 : where a † and a are creation and annihiliation operators and indices are chosen such that The coefficients c ai and c bj describe the coupling of these single orbital transitions for state localized on A and state localized on B. The diagonal matrix elements of the dia LE block are taken from two separate FDEu (where u denotes uncoupled) TDDFT (Casida and Wesołowski 2004;Neugebauer 2007Neugebauer , 2008Neugebauer , 2010König et al. 2013) calculations by selecting the states of interest. In this work, we are interested in the lowest two states Q y and Q x for each of the chromophores. These energies are also called site energies in the context of interacting chromophores and are obtained as solutions of Here, sub and sub is the subsystem partitioned response matrices (Neugebauer 2007;König et al. 2013) assuming real spatial orbitals, { } the set of eigenvalues (site energies) and { sub , sub } the set of corresponding excitation and de-excitation vectors for each of the Φ A * B and Φ AB * states, respectively [see for example Eq. 23 in Ref. König et al. (2013)]. Note that by choosing the FDEu starting ansatz the dia LE matrix becomes diagonal inside the Φ A * B and Φ AB * subblocks. The coupling between the LE states Φ A * B and Φ AB * , is subsequently approximated by the coupling between the transition densities of excitations localized on fragments A and B and yields the remaining non-zero elements of the dia LE block. In the Tamm-Dancoff approximation (denoted as FDEc-TDA) it is given by König et al. (2013), where A, ( ) and B, ( � ) are the transition densities of excitations and corresponding to Φ A * B and Φ AB * , respectively, and the kernel f ( , � ) is given by where E XC [ ] and T s [ ] denote the usual exchange-correlation and kinetic energy functional and tot is the summed density of subsystems A and B. The first term of Eq. 23 represents the Förster type Coulomb coupling between two localized excitations whereas additional short-range interactions are taken into account via the approximate exchangecorrelation and kinetic energy functionals in the second and third term. The remaining elements to be considered are the overlap matrix elements, these are taken as elements of the unit matrix so that diagonalization of the LE submatrix will reproduce the FDEc-TDA results. We would like to note here that the above formulation of subsystem TDDFT lacks an explicit description of charge transfer excitations (which we calculate in the next section using explicitly constructed charge localized states) that become significant when two chromophores are sufficiently close to each other . Recent works by Neugebauer and co-workers using projection based embedding techniques (as an alternative to the conventional non-additive functionals used here) in combination with a supramolecular basis set in the context of (21) subsystem TDDFT have been shown to be able to describe such inter-subsystem charge transfer excitations and reproduce the supramolecular results exactly (Tölle et al. 2019a(Tölle et al. , 2019bScholz et al. 2020;Tölle and Neugebauer 2022). We have not considered this extension of FDE here.

Calculation of dia CT block
A detailed discussion of the computation of CT energies with subsystem DFT and other approaches can be found in the literature (Pavanello and Neugebauer 2011;Pavanello et al. 2013;Solovyeva et al. 2014;Ramos et al. 2015). Here we will give the essential points, largely following Marcus Theory. We start by considering the charge transfer process in terms of (quasi-)diabatic charge localized states and first consider the particular case of Φ CT1 defined in Eq. 8 and describing the process where an electron is transferred from A to B. In order to calculate the excitation energy E CT1 of state CT1 ( A + + B − ) with respect to the GS ( A + B ) one then solves a 2 × 2 generalized eigenvalue equation (Pavanello and Neugebauer 2011;Pavanello et al. 2013;Migliore 2009Migliore , 2011. This yields two adiabatic state energies with their difference given by, where and we abbreviated CT1 by 1 and the GS by 0. H dia 00 , H dia 11 , H dia ab and S 01 represents the diabatic energies of the GS and CT1, the coupling between them and their overlap, respectively. When the energy difference between the diabatic states is much larger than the coupling ( H dia 00 − H dia 11 ≫ V 01 ), the diabatic and adiabatic energies coincide, i.e. Eq. 25 at small overlap and coupling becomes, Similar considerations apply for the reverse CT state Φ CT2 defined in Eq. (9).
In the subsystem DFT framework, these CT states are constructed from broken-symmetry KS determinants. Orbitals are optimized following Eq. 15 in a series of freeze-andthaw cycles, for two neutral fragments A and B comprising the ground state ( Φ GS ) and two charge-separated states A + and B − (or A − and B + ) comprising the product state Φ CT1 (or Φ CT2 ). These in total six sets of subsystem KS orbitals are combined to generate the three diabatic states Φ GS , Φ CT1 and Φ CT2 . In contrast to the procedure followed for the LE states, where it is easy to consider many excited states as they all are generated from the same GS orbital set, the procedure for CT states will produce only one CT1 state and only one CT2 state, each with it own set of optimized orbitals that include the effect of polarization via the freeze-thaw optimization step. For example, Φ CT1 in Eq. 8 is constructed as, where i A+ , j B− are the set of orbitals for the corresponding charge localized cation A + and anion B − . These states are then used to calculate the matrix elements in Eqs. 25 and 26 as, Here (01∕2) ( ) is the transition density between states 0 and 1 or 2 (i.e. between Φ GS and Φ CT1 or Φ CT2 ) and S 01∕2 the corresponding overlap. (ii) ( ) is the density calculated using orbitals of the appropriate diabatic states 0, 1 or 2 ( Φ GS or Φ CT1 /Φ CT2 ). We need to account for two kinds of non-orthogonality in these calculations. Since the subsystem orbitals are the results of an FDE calculation, their product wave function is not necessarily normalized so that S ii ≠ 1 . Furthermore, the off-diagonal elements S 01 and S 02 are in general non-zero ( S 12 is zero, as mentioned before). This means that rather than focusing on only the blue part of Fig. 3, as is our intention, we also need to check the magnitude of the GS-CT coupling, or equivalently, the difference between the diabatic and adiabatic CT energies. As we will discuss later, we are in the current work in the weak coupling regime, which means that the desired dia CT elements can be taken as the approximate energy difference Eq. 27 and we need not consider the GS explicitly.
We would like to point out here that the CT energies so computed using the subsystem formalism have the correct long-range behaviour, which the approximate functionals employed in conventional supramolecular TDDFT often fail to describe properly. The energy of a charge separated donor and acceptor in the subsystem formalism can be written as, where the E A+ , E B− are the diabatic energies of the two ionized states and E int is the interaction energy between the two charged species given by, At very long range, where there is zero overlap between A + and B − , the non-additive terms drop out, while the other terms yield the correct limit: E int ≈ − 1 R , where R is the distance between A and B (Solovyeva et al. 2014).
For the short-range or strong-coupling limit, we could in principle use Eq. 26 and include the GS explicitly in the calculation, but such an approach would require careful testing of the validity of the employed density functional approximation in this limit.

Calculation of dia LE−CT block
Having discussed the diagonal blocks of dia , we now turn our attention to the coupling of the CT and LE states given by the dia LE−CT block. The most rigorous approach would be to consider the broken-symmetry KS determinants of the preceeding section, as well as the full expansion of the TDA-TDDFT wave function, thereby accounting for the fact that these are constructed from different subsystem orbital sets and are not necessarily orthogonal. Here, we will introduce a simple approach using second quantization for non-orthogonal orbitals in order to calculate the necessary couplings. The approach described here represents a generalized form of the Fragment Orbital DFT (FODFT) approach previously proposed by Senthilkumar et al. (2003Senthilkumar et al. ( , 2005 for calculating charge transfer integrals for non-orthogonal orbitals and used later by Hernández-Fernández et al. (2016) for calculating electronic couplings for hole transfer in stacked porphyrin dyads.
In this approach we approximate the full interaction Hamiltonian by the fixed Kohn-Sham (KS) Hamiltonian matrix representation obtained in a supermolecular calculation. This Hamiltonian is diagonal in its eigenbasis of supermolecular orbitals and can be expressed in second quantization as: where { p } are the supermolecular orbital energies.
We then define a normal ordered Hamiltonian with the supermolecular ground state determinant Φ GS as the reference vaccum, where we have redefined the usual creation/annihilation operators ( a † /a) in Eq. 34 in terms of hole/particle operators ( b † /b) in Eq. 35 where the indices i, a denote the supermolecular occupied and virtual orbitals. We then construct a set of unpolarized reference fragment orbitals obtained from two independent fragment calculations of A and B which constitutes the above defined reference (supermolecular) vacuum (valid for weakly interacting fragments). Assuming that the excited states can be effectively described by singleparticle excitation in this basis and that these underlying set of reference orbitals do not change, the LE state Φ A * B and CT state Φ CT1 in the hole/particle formalism can be written as, Similarly, the coupling between the states Φ A * B and Φ CT2 is given as, Given that CT2 occurs only from the highest occupied molecular orbital of a fragment B ( H B ) to the lowest unoccupied orbital of another fragment A ( L A ), the coupling with the two LE states on fragment A becomes, Equations 39, 40 and 42, 43 along with the corresponding coupling terms for the LE states on fragment B (see SI, not shown here) form the eight unique couplings in which we recognize the contributions of hole as well as electron transfer, which are, respectively, governed by the overlap between the occupied and the virtual orbitals of the individual fragments. Note that, in the above calculation of the couplings, we assumed the LE states ( Q y and Q x ) to be pure states (i.e. composed of a single single orbital transition), thereby making their calculation fairly straightforward, obtained readily from a ground state Kohn-Sham calculation of A + B . In practice however, the LE states are composed of multiple single orbital transitions (see Eq. 19), in which case a more complete description of the couplings can be obtained either by a subsystem approach (Difley and Van Voorhis 2011) or by diabatization techniques (Voityuk and Rösch 2002;Yang and Hsu 2013;Hsu et al. 2008;Nottoli et al. 2018;Tölle et al. 2020). For the present purpose of the work, the above approximation is considered to be sufficient for a qualitative assessment.
A summary of the above approximations introduced in the context of the different methods and the structure of the matrix dia and dia used in Eq. 10 are shown in Table 1 and Fig. 3, respectively.

Effect of a polarizable environment on
where the first term v el is the electrostatic operator and describes the Coulomb interaction between the QM system and the permanent charge distribution of the MM environment. The second term, v pol , describes the induced polarization at the MM atoms due to interaction with the QM part and other MM atoms. These two terms are given by with where the interaction tensor to a given order, n, is given by T (n) pq, 1 ..., n . R pq is the distance between the interacting entities, R si, is a component of the distance vector and ind s is the induced dipole at site s. The induced dipoles depend on the QM density and are updated during the self-consistent field cycles used to solve the QM part.
The DRF model has also been extended to the calculation of excited states in the linear response formalism of TDDFT where, an additional explicit response of the MM region DRF potential v DRF is accounted for in the kernel resulting from the response in the induced dipole upon change in the QM charge distribution. Similar to the calculation of the induced dipoles using the ground state density, in the excited state formalism it utilises the transition density for the calculation of the induced dipoles. It is defined by Jensen et al. (2003b) as where ind s, [ ] is the induced dipole at site s due to the perturbation in the density and T (1) si, is the corresponding interaction tensor. In all our subsequent calculations we restrict ourselves to the adiabatic local density approximation (ALDA) of the kernel, so the dependence on is dropped from the above equations. In the presence of a DRF environment, the total effective Hamiltonian eff can be written as, where el is our block diagonal electronic Hamiltonian in the absence of the environment defined above, and DRF is the Hamiltonian describing the perturbation on each of the blocks of el due to the DRF environment. In the following section, we note the nature of perturbation of DRF on the blocks of the dia matrix Effect on dia LE block In the calculation of the matrix elements of the dia LE block , the effect of the environment is accounted for in two steps-(1) As an additional potential term, v DRF , given by Eq. 44 in the calculation of the ground state subsystem orbitals via Eq. 15 and (2) In the calculation of the diagonal (LE energies) and off-diagonal (LE couplings) matrix elements in Eq. 22 via the linear response of v DRF in Eq. 47.
Effect on dia CT block The effect of the environment is explicitly accounted for in the calculation of the set of subsystem orbitals for the individual fragments, i.e.
in the relevant diabatic states belonging to either the reactant or the product state via v DRF given by Eq. 44 and entering the ground state Kohn-Sham calculation in Eq. 15. It therefore affects the CT energy (i.e. the diagonal sub ∈ A * B, AB * ; FDEu-TDDFT; Refs. Neugebauer (2007), König et al. (2013) König et al. (2013) dia CT (Eq. 27) X ∈ CT1, CT2 ; FDE-ET; Refs. Pavanello and Neugebauer (2011), Pavanello et al. (2013), Solovyeva et al. (2014), Ramos et al. (2015) dia Generalized FODFT; Refs. Senthilkumar et al. (2003Senthilkumar et al. ( , 2005; see SI for remaining couplings 1 3 elements of this block) and the electronic coupling with the ground state through Eqs. 25 and 26. Effect on dia LE−CT block The effect of the environment is taken into account in the ground state Kohn-Sham Hamiltonian , Ĥ KS , in Eq. 33 of the combined system of A + B as well in the calculation of the reference fragment molecular orbitals of A and B via v DRF , given in Eq. 44.

Computational details
All calculations have been carried out with the Amsterdam Density Functional (ADF) program (Baerends et al. 2018;Te Velde et al. 2001) using the double-zeta polarised (DZP) basis set and the range-separated CAMY-B3LYP exchange and correlation (XC) functional. The CAMY-B3LYP XC functional is implemented with the XCFun library (Ekström et al. 2010, https:// dftli bs. org/ xcfun/) and is a modified version of the original CAM-B3LYP (Yanai et al. 2004) with a different switching function. The performance of this functional in conjunction with the basis set for the calculation of the Q band energies has shown reasonable agreements before (López-Tarifa et al. 2017) with Milne et al. (2015) and hence was used consistently throughout this work. For the Frozen Density Embedding (FDE) calculations, the GGA functional BLYP (Becke 1988;Lee et al. 1988) and PW91k (Lembarki and Chermette 1994) were used for the non-additive exchange-correlation and kinetic energy part of the embedding potential through out this work. For the calculation of the elements of dia LE block, the Time-Dependent DFT (TDDFT) extension of the Frozen Density Embedding (FDE) scheme in the linear response regime (Neugebauer 2007;König et al. 2013) as implemented in ADF was used to calculate two lowest excitation ( Q y and Q x ) for each of the chromophores, Chla611 and Chla612. In order to include the polarization effects, 3 freeze-and-thaw cycles were performed for each of the chromophores. The coupling between the local excitations was then subsequently calculated using the Tamm-Dancoff approximation in the FDEc formalism (König et al. 2013) as mentioned before. For the calculation of the diagonal elements of the dia CT block, i.e. the CT energies, separate unrestricted calculations in the FDE framework were performed for each of the chromophores in the neutral ground state and the charge-separated state to generate a total of 4 unrestricted fragments. Three freeze-and-thaw cycles were performed for each of the fragment in order to introduce polarization effects. For the inter-subsystem contributions,we used the BLYP (Becke 1988;Lee et al. 1988) XC functional to evaluate Eqs. 29 and 30. A schematic overview of the workflow is shown in Fig. 4. The elements of the dia LE−CT block were calculated using a recently developed stand-alone code Reduction of Orbital Space Extent (ROSE), specifically designed for localization of molecular orbitals. Senjean et al. (2021) The environment was modelled in the DRF framework as mentioned earlier, with Mulliken charges obtained from a Self-Consistent Charge Density Functional Tight Binding approach (Elstner et al. 1998) (SCC-DFTB) with third order corrections, using the parameter set 3ob-3-1 (Lu et al. 2015) as implemented within the AMS ) engine of ADF. The atomic polarizabilities were taken from the standard Thole's set of atomic polarizabilities (Thole 1981; Van Duijnen and Swart 1998) as given inside ADF. [In addition, for the Magnesium and Phosporous atom, polarizability values were taken from Ref. Stout and Dykstra (1995) and Lupinetti and Thakkar (2005), respectively]. All the elements of the dia matrix were calculated on a total of 103 frames with 51 in the first set and 52 in the last set separated by 400 ps in each window, from the classical trajectory [labelled as trajectory  (Liguori et al. 2015)] computed using GROMACS Molecular Dynamics package (Van Der Spoel et al. 2005). Such a separation of the frames was motivated by the conformational changes observed in the particular selected trajectory between Chla612 and Chla611 towards the end of the simulation in reference (Liguori et al. 2015). The input geometric structure from each frame of the trajectory for the QM calculations was prepared in two steps. First, we use the visualization program VMD (Humphrey et al. 1996) to visualize and select chromophores Chla612 and Chla611 and the neighbouring residues (protein, chromophores and lipids) within a distance of ≈ 7 Å from each of the chromophore porphyrin rings in the first and last set of snapshots of the trajectory and save the coordinates for the resulting snapshots. In the second step, we process each of the frames using the RDKit (http:// www. rdkit. org) library, wherein the DRF environment is separated from the QM region consisting of Chla612 and Chla611, followed by the addition of all the missing hydrogens, correction of the bond order of the porphyrin ring of the two chlorophylls and removal of the phytol chains of the chlorophylls in order to reduce the computational cost [removal of the phytol chain have been shown previously to not significantly affect the site energy calculations (López-Tarifa et al. 2017)]. Thereafter, we calculate all the elements of the dia in three different sets of calculations for each frame, for the dia LE , dia CT and dia LE−CT blocks, respectively, as mentioned above. The generation of the input geometries and the corresponding preparation of inputs for the subsequent ADF calculation were done in a fully automatised way using the Python Library PLAMS (Handzlik et al. 2018) as interfaced with ADF. As a final note, we would to mention that the numerical costs are dominated by the (not well optimized) calculation of the dia CT block which takes 17 h of the 20 h in total needed per snapshot on one node (with 128 cores) of the Dutch supercomputer Snellius.

Results and discussion
The supramolecular picture Förster resonance energy transfer (FRET) (Förster 1948) and Marcus electron transfer (ET) (Marcus and Sutin 1985) processes are both usually described in terms of local and non-local diabatic states. This is straightforward in subsystem DFT, but also in supermolecular TDDFT with approximate exchange-correlation functionals, the resulting adiabatic states can sometimes still be interpreted in a diabatic picture. This is especially so if the local Frenkel (LE) type excitons and the charge-separated (CT) states are energetically well separated. This identification is difficult, however, when low-lying CT states, occurring either as an artefact of the method itself or otherwise, start to mix with these local excitations. In a previous work ) on the Chla611-Chla612 dimer, we studied such a mixing of the CT states upon decreasing the distance between the two chlorophylls from 9 to 5 Å. In that work we used a supermolecular approach but quantified the amount of mixing in terms of a charge transfer number (Plasser and Lischka 2012) which indicated a sharp increase of LE-CT mixing upon decreasing the interchromophoric distance. However, due to the delocalized nature of the Highest Occupied and Lowest Unoccupied Molecular Orbitals (HOMO and LUMO) for that conformer, a clear interpretation in terms of diabatic LE and CT states remained elusive. In this work we provide such an interpretation and construct the individual matrix elements of dia to get a better picture of the couplings and energetics of these diabatic states. In the next section, we focus on the individual matrix elements of dia for selected snapshots of the MD trajectory spanning the beginning (sampled by first set of snapshots) of the simulation where the Chla611 and Chla612 chromophores are still far apart as well as the end of the simulation (sampled by last set of snapshots) where they are relatively close.

Local excitation (LE) and charge transfer (CT) energies
In Fig. 5 we show the four relevant ( Q x1 , Q y1 , Q x2 , and Q y2 ) local excitation energies of the two chlorophylls (Chla612 and Chla611) along with two lowest CT state energies, CT1 ( Chla611 + ∕Chla612 − ) and CT2 ( Chla612 + ∕Chla611 − ), which represent the diagonal elements of the blocks dia LE and dia CT of the diabatic matrix respresentation (Fig. 3). To quantify these values we also list the averaged LE and CT energies in Table 2. As can be seen from Fig. 5 and Table 2, the LE energies are minimally effected (the lowest Q y2 is red-shifted by ≈ 0.02 eV) by the changed orientation and reduction of dimer distance in the chromophore pair. These energies are also relatively robust with respect to the smaller geometry fluctuations occuring between different snapshots with a variance of only 0.001 eV. The CT (CT1 and CT2) energies on the other hand show both a much larger fluctuation between snapshots and a clear decrease in energy going from the more distant geometries in the first set of snapshots to the close arrangement in the last set. The fluctuations of the CT energies are not unexpected as CT states have substantially larger reorganization energies than the local excitations as the Franck-Condon point is further away from the minimum than is the case for LE states. Looking closer at the overall trend for these CT energies we note the drop in energies of both of the CT states that we anticipated from the earlier supermolecular analysis. We 1 3 see in Fig. 5 both the CT energies coming down and even crossing the Q x line at some points. In addition we observe that the relative gap between the two CT states decreases from an average of |ΔE CT | of 0.51 to only 0.07 eV. The latter observation can be explained by the asymmetric orientation and thereby polarization (as represented by the DRF embedding potentials) for CT1 and CT2 in the beginning of the trajectory. In the final part of the studied trajectory the arrangement of the chromophores is more symmetrical and the polarization in the CT states is more similar, thereby bringing the energies of these two CT states closer together. Still, with the absence of an exact symmetry, the two CT states are not fully equivalent also in this stage.
From Eq. 31, the CT1 and CT2 energies, |ΔE CT | not only depend on the absolute energies of the charge localized fragments (i.e. E A+ + E B− and E A− + E B+ ) but also on their interaction energies ΔE int . The latter, negative, contribution to the CT energy is the driving force for the lowering of the CT energy towards the end of the trajectory. This enhances the mixing of the CT states with the local states and makes switching to an adiabatic picture by either diagonalizing the diabatic matrix representation or by using supermolecular TDDFT necessary. In the first part of the trajectory, the adiabatic (black) curves align well with the diabatic LE states, with only for the CT states that are hard to describe by density functional approximations significant deviations visible. For the structures from the final part of the trajectory, also the lowest supermolecular state is clearly below the lowest diabatic LE state, corroborating the influence of CT states in lowering this energy.
Effect of the environment on LE and CT energies In order to assess the effect of the protein environment on the LE and CT energies we repeated the calculation of LE and CT energies of the dimer in vacuum. Figure 6 shows the effect of the environment on the diabatic states in the beginning and at the end of the trajectory. The average site energies and the CT energies with and without the environment are shown in Table S1 in the supplementary information. We note that: (i) all the averaged site energies are red-shifted by the environment, both in the beginning and at the end of the trajectory, with the exception of the Q y2 and Q x2 which remain unaltered in the beginning of the trajectory. (ii) The CT energies are most strongly influenced in the beginning of the trajectory, with both the CT energies shifted by the environment in opposite directions (CT1 red-shifted and CT2 blueshifted) thereby increasing the average gap between them. This effect can be explained by the asymmetric distribution of partial charges in the immediate proximity of both the chlorophylls (in particular the DPPG ligand near Chla611) in the beginning of the trajectory that we already mentioned in the preceding section. Towards the end of the trajectory, on the other hand, the stacked and symmetric arrangement of the chlorophylls with a reduced interchromophore distance,  renders the two charge-separated states less sensitive to external perturbations. Looking at the magnitude of the environment shifts we conclude that overall the environment plays a minor role in determining the energies of the LE and CT states.

Coupling between the diabatic states
We now consider the couplings between the diabatic states (off-diagonal elements of dia , see Fig. 3) that determine the final adiabatic picture. The off-diagonal elements in the LE block represent the predominantly Coulomb-like coupling between the local excitations which are responsible for Förster resonant energy transfer (Förster 1948), whereas the elements in the LE/CT block represent the overlap-dependent couplings between the local and CT states, which play an equally important role in Marcus electron transfer theory (Marcus and Sutin 1985). We note that these LE/CT couplings can describe both electron transfer between the fragment LUMOs (one coupling element), as well as hole transfer involving the two sets of HOMO and HOMO-1 located on each fragment (three possibilities), excluding transfer between the two HOMO-1 orbitals that is of minor importance). Table 3 and Figure 7 show the values of all the relevant LE/LE and LE/CT coupling elements: between the two LE Q y states ( V Q y ), the two LE Q x states ( V Q x ), and between these and the two possible CT states CT1 and CT2, as discussed in the previous sections: ), where A and B are Chla611 and Chla612, respectively. Note that the coupling V Q A x ,CT1 (and V Q B x ,CT2 ) is negligible owing to the mutual orthogonality of the HOMO and HOMO-1 for each of the fragments A and B (see Eq. 40) and is not shown. From Fig. 7 the much smaller variation for the LE/LE coupling compared to the LE/CT coupling stands out. The variation of the LE/LE coupling with the interchromophore distance R can be qualitatively described by an interacting dipole model (Liguori et al. 2015;López-Tarifa et al. 2017) and exhibits a ∼ R −3 (dipole-dipole) dependence. This can be contrasted to the exponential dependence ∼ e − R∕2 on this distance that is applicable for the LE/CT coupling (Cave and Newton 1996). Both couplings are enhanced in the final part of the studied trajectory in which interchromophoric distances are much smaller than that in the initial frames.
Effect of environment on the couplings The environmental effect on the LE/LE and the LE/CT couplings are shown in Fig. 8 as differences in the absolute couplings. Table S2 of the supplementary information also lists the averaged couplings in vacuum and in DRF environment. In case of non-orthogonal orbitals, the LE/CT couplings are not only determined by the different overlap terms, but also by the supramolecular KS orbital energies (see section on calculation of dia LE−CT block). The variations of these LE/ CT couplings upon inclusion of the environment are a direct consequence of the effect of the environment on both of these quantities. A further detailed analysis of this effect is beyond the scope of this current work. For the LE/LE couplings we see both for the Q y1 , Q y2 and Q x1 , Q x2 a more significant effect than seen previously for their energies. The increase in coupling upon including the environment effects (so relative to a vacuum treatment) can be explained by the  S1 in supplementary information). This explicit treatment of the environment makes it possible to avoid the often applied scaling of transition dipole moment to achieve better agreement with experimental data. As such scaling does not provide additional insights we do not apply this in the present work. We note that our results here are in reasonable agreement with the recent works of Sláma et al. (2020) with minor differences arising from the different levels of theory and MD simulations used in both works. For comparison of the LE/LE couplings using the FDEc-TDA approach for the lowest Q y bands to other methods, we refer to López-Tarifa et al.

The adiabatic picture
With the energy gap between the LE and the CT states decreased and their coupling increased at shorter distances (the last part of the trajectory), we expect to see an increased mixing between the CT states with the LE states. It had been shown previously from a disordered averaged exciton model, that such coupling of CT states to LE states can strongly modulate the optical spectra . Such a pronounced effect of the mixing of the CT states with the LE states at the end of the trajectory can indeed be seen when switching to the adiabatic picture. Solving for the generalized eigen value problem (Eq. 9) in the basis of the six diabatic states ( Q y1 , Q x1 , Q y2 , Q x2 , CT1 and CT2) for the first and last set frames of the trajectory, shows an average red-shift of about 0.02 eV of the lowest state. This state acquires clear mixed LE/CT character towards the end of the trajectory which is in qualitative agreement with the earlier supermolecular analysis. This indicates that the composite manner of constructing the matrix that is to diagonalized is adequate. To compare better with the supermolecular picture we show both the new set of adiabatic states as well as the six lowest supermolecular excitations in Fig. 9. Like already seen in the diabatic picture (Fig. 5), where there is little interaction between the CT states and the local states, the lowest lying adiabatic states follow closely to those of the supramolecular states and retain their local character. The two highest lying adiabatic and supramolecular states on the other hand, represent states which are predominantly CT in character, and get significantly different energies in the two approaches. The supermolecular ones are hereby likely to be to low in energy due to flaws in the treatment of CT states with the available density functional approximations where the adiabatized states of the current approach could be too high due to missing couplings with states at higher energies. Nonetheless, this study provides a qualitative assessment of the importance of accounting for CT states and electron delocalization when studying the energy landscape of pigment-protein complexes, particularly in the case of strongly interacting pigments. Furthermore, it shows how pigment reorganizations following protein conformational changes can greatly affect the spectral properties. This spectral tunability is fundamental for LHCs to be able to switch between different functional states. Moreover, an improved knowledge of the nature of these dimer excitations in different orientations can provide hints to understand the effects of artificial modifications of the Chla611-Chla612 binding pocket on the overall absorption spectrum of the complex.

Conclusion
In this work, we looked into changes in energies of the lowlying states of the chlorophyll dimer Chla611-Chla612 induced by a conformational change of the dimer in the highly disordered terminal emitter domain reported previously (Liguori et al. 2015). In order to get a better insight into these states, we elucidate the effect of this conformational change on the LE and CT state energies of the dimer along with their corresponding couplings using subsystem DFT calculations on snapshots from a s MD simulation generated previously. Subsystem DFT provides an alternate Fig. 9 The six adiabatic states labelled E1-E6 are shown in bold in increasing order of their energy from the beginning and the end (left and right of the split on x-axis) of the trajectory. Also shown in dotted lines are the lowest six supramolecular states. All units are in eV localized diabatic approach as compared to the delocalized picture arising from supramolecular DFT, providing an unambiguous identification of the local and CT states. We find that the re-orientation of Chla611-Chla612 dimer upon equilibration brings down the energies of the Chl-Chl CT states, causing an increase in the mixing of the CT and LE states (which remain almost unaffected), accompanied by an increased coupling between them. This, in combination with an increase in the coupling between the local states themselves, facilitated by an reduced interchromophoric distance, gives rise to red-shifted low-lying mixed excitonic/CT states. Previous studies have shown that such mixing of low-lying charge transfer states with the Q y band can be operative in the the Chla band (Romero et al. 2009;Ramanan et al. 2017;Novoderezhkin et al. 2016;Chmeliov et al. 2016) and are indicative of the so-called 'red-emitting states' in the fluorescence spectra as observed in aggregates, supported by femtosecond Transient Absorption(TA) kinetics and Stark spectroscopy studies (Kell et al. 2014;Krüger et al. 2010;Müller et al. 2010). Moreover, Chl-Chl CT states have also been shown to play an important role in regulating energy flow through antenna complexes and create energy sinks in the reaction centers (RC) of photosynthetic complexes (Madjet et al. 2009;Raszewski et al. 2008;Novoderezhkin et al. 2005Novoderezhkin et al. , 2007. Our study therefore provides further evidence that low-lying CT states, prompted by a favourable conformational change of the chlorophyll dimer, can play an pivotal role in regulating light harvesting and can create energy sinks facilitated by an increased excitonic interaction.