Short-Range Effects in the Special Pair of Photosystem II Reaction Centers: The Nonconservative Nature of Circular Dichroism

Photosystem II reaction centers extract electrons from water, providing the basis of oxygenic life on earth. Among the light-sensitive pigments of the reaction center, a central chlorophyll a dimer, known as the special pair, so far has escaped a complete theoretical characterization of its excited state properties. The close proximity of the special pair pigments gives rise to short-range effects that comprise a coupling between local and charge transfer (CT) excited states as well as other intermolecular quantum effects. Using a multiscale simulation and a diabatization technique, we show that the coupling to CT states is responsible for 45% of the excitonic coupling in the special pair. The other short-range effects cause a nonconservative nature of the circular dichroism spectrum of the reaction center by effectively rotating the electric transition dipole moments of the special pair pigments inverting and strongly enhancing their intrinsic rotational strength.

The coordinates of the D1D2cytb 559 complex were extracted from that of the PSII core complex (pdb file 3WU2.pdb). 2 After adding hydrogen atoms, the protonation probabilities of titratable residues were calculated and molecular dynamics (MD) simulations were performed.110 snapshots of the MD simulations were geometry optimized by QM/MM, including the special pair in the QM part.Based on an analysis of the snapshots, we define the snapshot with the lowest root mean square deviation (RMSD) to the average special pair coordinates as the representative special pair structure.Excited states of the special pair were calculated treating the remaining pigments and the protein environment by atomic partial charges (electrostatic embedding), for the representative structure and the remaining snapshots.These calculations are analyzed with the diabatization revealing local excited and charge transfer states as well as couplings.
From the distribution functions of these quantities, obtained from the snapshots, static disorder is estimated.
For the representative structure, in addition, quantum chemical calculations of excited states of the isolated monomers are performed, which are used to obtain the long-range excitonic coupling in the special pair.The energies and couplings obtained for the representative structure are post-processed, based on independent experimental data.The effects of charge transfer states are accounted for by perturbation theory, i.e., Q y energies and dipoles, as well as couplings comprising at least one Q y state are perturbed.The excitonic couplings between the special pair and the remaining RC pigments are obtained with the Poisson-TrEsp method for the representative structure.The parameters for the remaining RC pigments are taken from our earlier work. 1 Using these parameters of the Frenkel exciton Hamiltonian, the circular dichroism spectra are calculated in the full model, including higher excited states of the pigments, and a minimal model that takes into account the Q y transitions only, but includes an enhanced excitonic coupling in the special pair due to coupling to CT states and a rotation of electríc and magnetic transition dipole moments by non-CT short-range effects.

Isolated monomers
Poisson−TrEsp  The method for determining the protonation probabilities of titratable residues has been reviewed 3 and applied to the reaction center of PSII. 4 The overall process is given by i) preparing the structure for the PSII reaction center, ii) solving the Linearized Poisson-Boltzmann Equation (LPBE), and iii) running a Monte-Carlo titration.

Remaining reaction center (RC) pigments
The used structure is based on the crystal structure for the cyanobacterial PSII of T. vulcanus (PDB 3WU2) provided by Umena et al. (2011), 2 which shows a resolution of 1.9 Å.The structural elements contained in the full reaction center model are based on ref, 4 but omitting the J and X chain as well as the atoms of the water-oxidizing complex Mn 4 CaO 5 (WOC), since these structural elements are expected to be absent in the D1D2cytb 559 PSII-RC preparations on which the optical spectra are measured.Special attention has been paid to the protonation states of the ligands of the missing WOC, which were supposed to compensate for its charges.All ligands belonging to the reaction center have been taken into account.Neglecting the CP43-Glu-354 ligand can be justified by its absence in the PSII-RC preparations.
In the second step, the LPBE is solved using TAPBS. 5,6After the definition of the protein and membrane, as well as the assignment of dielectric constants (ϵ mem = 2, ϵ p = 4, ϵ solv = 80), the LPBE is solved on a grid using a finite difference method and a grid focusing approach.The model for the PSII reaction center is depicted in Figure S2.Standard parameters have been used for the calculation (e.g., temperature: 300 K, ion radius: 2.00 Å, etc.).
Finally, the calculation of the average protonation state for each titratable group has been conducted with the Karlsberg2 software, 7 where a range of pH-values between 4 and 11 has been considered.

Results
In the following, the most probable protonation states for the histidines, the titratable groups of the ligands, as well as other titratable groups are presented.Regarding its protonation state, histidine represents an exception to the other considered amino acids.At physiological pH-values, histidine may occur in its protonated (HSP), or two different deprotonated forms (HSD, HSE).The obtained results are given in Table S1.The protonation states of the WOC ligands are presented in Table S2.It can be readily seen, that the missing charges of the Mn 4 CaO 5 complex are compensated by the obtained protonation states of its ligands, as expected.Finally, only one other titratable group (Glu-D312) has been found to deviate from its standard Table S1: Protonation states of histidines protonation state (deprotonated) at a pH-value of 6.5.However, since Glu-D219 shows a probability of only 50% to be protonated, it is still modelled in its standard form.The results can be compared to current findings in the literature.In this regard, Narzi and Guidoni 8 conducted MD simulations of the Apo-PSII with different protonation states for the WOC ligands.The considered protonation pattern which is closest to the results presented above has been taken from Han et al. (2019). 9This pattern has been found to be the most stable one in terms of RMSD during the MD simulation of Narzi and Guidoni.

Method
Molecular dynamics preparation After the most probable protonation states were determined, the system was embedded in a DOPC membrane with KCl salt concentration of 0.15 mol using CHARMM-GUI. 10The protein and cofactors were parameterized using the Amber ff14SB Force-Field (FF), 11 while the lipid17 FF used to describe the membrane and the TIP3P 12 model for water molecules.Missing parameters for metalligand complexes were calculated using the Seminario method 13 as implemented in the MCPB protocol. 14veral types of metal-ligand complexes are found in our PSII model.These include: six chlorophylls (Chls), the heme of Cytb 559 , and the non-heme iron at the interface of proteins D1 and D2.Two sets of parameters were calculated for the Chls.The first describes a Mg 2+ ion coordinated to five nitrogens, the four nitrogens of the chlorin ring and the ϵ nitrogen of a histidine residue.The second set of parameters only describes the coordination between the Mg atoms and the four chlorin nitrogens.The former was used to describe the Chls P D1 , P D2 , Chlz D1 and Chlz D2 , whereas the latter was used to describe the Chls Chl D1 and Chl D2 .The cytb 559 heme was parameterized as an octahedral complex of low-spin Fe +2 iron coordinated to six nitrogens: two ϵ nitrogens of histidine, and the four nitrogens of the porphyrin.The heme propionate functional groups were taken to be deprotonated.The non-heme iron was also parameterized as an octahedral complex, however this time the ligands are four ϵ nitrogens of histidine and two oxygen atoms of bicarbonate.The bicarbonate was protonated at the oxygen facing away from the iron, and the complex was modeled as a high-spin Fe +2 complex.The last two residues which require special attention are alanine with protonated C-terminus,and a methionine with formylated N-terminus.For these residues, we used the prepgen tool from the Amber toolbox. 15All missing parameters and RESP charges were calculated in the gas-phase at the B3LYP 16 /6-31G* 17 level of theory using Gaussian16. 18lecular dynamics simulation Starting with our reduced PSII model, initial geometry optimization was carried out with spatial restraints applied to all heavy atoms, excluding water oxygens.For backbone atoms, a restraint weight of 10.0 kcal/mol was used, while for all other atoms (side chain, ligands, lipids and ions), a restraint of 5.0 kcal/mol was applied.Following the optimization step, a heating simulation from 100 K to 300 K was carried out using a NVT ensemble for 800 ps while maintaining the spatial restraints.After the heating simulation was complete, the final geometry and velocities were used as a starting point for a ten step equilibration protocol, in which the spatial restraints were gradually removed (see table S3).Each equilibration step was carried out for 500 ps using a NPT ensemble with semi-isotropic pressure scaling.
After all restraints were removed, the simulation was continued for additional 5 ns.During these parts of the simulation, a time step of 1 fs was used.The last geometry of the equilibration simulation was used as the starting point of a 240 ns production run for the NPT ensemble.Here, new velocities were assigned based on a Boltzmann distribution, and the time step was increased to 2 fs.During all steps of the molecular dynamics, SHAKE constraints were applied to all hydrogen involving bonds, and a Langevin thermostat with collision rate of 5 ps −1 was used.S3: Restraint weights during the equilibration steps in kcal/mol .

Results
To measure the quality of the MD simulation, we investigated the backbone RMSD with respect to the crystal structure.In this way, we can see, if our reduced model is stable, and examine how far it drifts away from the original structure.Figure S3 shows the backbone RMSD of the initial production run.Overall, the RMSD shows stability over the course of the MD simulations and the backbone does not deviate much from the crystal structure.However, within the first 40 ns of the simulation, the fluctuations in the RMSD are still relatively rapid.Therefore, these first 40 ns were taken as equilibration phase.
Figure S4 shows the backbone RMSDs of the second stage production run.Since all MD simulations reveal a stable structure, these simulations are used for sampling geometries representing the static disorder in the system.New velocities were assigned and additional 100 ns production simulations were performed.From each trajectory, eleven snapshots with equal time spaces were sampled.All the sampled snapshots (including the representative structure) were optimized at the MM level, using the conjugate-gradient algorithm in AMBER.
QM/MM Optimization All MM optimized structures went through a further QM/MM optimization protocol using the ORCA 5.0.3 software. 19QM/MM partitioning was done by truncating the C α − C β bonds of the Chl-ligating histidines and the P 1 − P 2 bonds of the phytyl chains (atoms denoted as C 1 − C 2 in the pdb file).The truncated bonds were saturated with hydrogen link atoms (Fig. S5).The QM-MM interaction was calculated by electrostatic embedding.Since ORCA does not support periodic boundary conditions, all MM atoms were fixed to their positions and only the QM part was optimized, and non-bonded interactions were calculated explicitly (no cutoff).In order to reduce computation time, we initially optimized the structures at the XTB semi-empirical DFT level of theory, 20 using the ORCA L-BFGS optimizer.Next, a second QM/MM optimization step was carried out at the B3LYP/def2-SVP 21 level of theory with Grimme's D3 dispersion correction with Becke-Johnson damping. 22In this stage, the RIJCOSX approximation 23 was used for calculating the Coulomb and HF-exchange integrals.
Figure S5: QM/MM partitioning.The QM region is made of the special pair Chls and their their fifth-ligand histidine side chains.Phytyl tails were omitted for the QM region by truncating the P 1 − P 2 bond and saturating it with a hydrogen link atom.

Results
The representative structure of the special pair is provided as an additional xyz-file "SP_Repr_Structure.xyz".
5 Parameterization of special pair Hamiltonian

Method
Excited states calculation The first 30 adiabatic excited states of the special pair structures mentioned above have been calculated by linear response TD-DFT at the ωB97X-D3 24 /6-31G* level using the AMBER-Terachem interface 25 with a direct cutoff of 50 Å for the QM atoms.In this way, we can take advantage of Terachem 26,27 GPU acceleration and calculate excited states for a large number of snapshots within a feasible time frame.In addition, the representative structure was subject to excited state calculation at the ωB97X-D3(BJ) 28 /def2-SVP level of theory using ORCA. 19In order to save computation time, the RIJCOSX approximation was applied.The preference for the def2-SVP over the similar 6-31G* basis set is rooted in the need for an auxiliary basis set when using the RIJCOSX approximation i , which is explicitly given for the def2-SVP basis set (def2/J). 29Furthermore, initial comparative calculations revealed only very minor differences regarding relevant results and very similar computation times between the 6-31G* and def2-SVP basis sets.

Diabatization
In the following, the Multi-FED-FCD 30 method is outlined, which has been applied for diabatization.As its name suggests, this method builds on two previously proposed diabatization procedures, i.e., the Fragment Charge Difference (FCD), 31 and the Fragment Excitation Difference (FED) 32 method.
These methods take two adiabatic states, the diagonal Hamiltonian in the adiabatic representation, and a certain operator in the adiabatic basis as inputs and yield two diabatic states and the Hamiltonian in the diabatic representation.The underlying rationale is to find an operator which i) is diagonal in the diabatic representation with diagonal elements being at their maximum/minimum, and ii) can be used to sort diabatic states according to their eigenvalues with respect to this operator.Taking the FCD method with the charge difference operator ∆q mn = ∫︁ rϵA ρ mn (r)dr − ∫︁ rϵB ρ mn (r)dr (fragments: A, B; transition densities of states m, n: ρ mn ) as an illustrative example, the procedure of the fragment difference methods follows the subsequent steps: 1. Calculation of ∆q in the adiabatic basis and diagonalization to obtain the adiabatic-to-diabatic transformation U: with ∆q ′ being the charge difference operator in the basis of diabatic states a and b.

Transformation of the adiabatic Hamiltonian into the diabatic representation:
Following the same procedure but utilizing the excitation difference operator ∆x mn = ∫︁ rϵA ρ ex mn (r)dr − ∫︁ rϵB ρ ex mn (r)dr, the FED scheme allows to identify different LE states instead of CT states.Here, the excitation i https://sites.google.com/site/orcainputlibrary/basis-sets/ri-and-auxiliary-basis-setstransition density ρ ex mn = ρ elec mn + ρ hole mn is defined in the CIS-picture as with Φ denoting a molecular orbital, and c (m) ia being the CIS coefficient which describes the excitation from an occupied molecular orbital Φ i to a virtual molecular orbital Φ a of state m.At this point, the similarity between the transition density and excitation transition density, ρ ex mn = ρ elec mn + ρ hole mn and ρ mn = ρ elec mn − ρ hole mn , should be highlighted, suggesting that the combination of FCD and FED is efficient to compute.A nice illustration of the above idea of utilizing electron and hole densities to identify LE-and CT states can be found in Fig. 1 of the recent work by Wang et al. 33 The above methods show some shortcomings like, e.g., the open question of how to best select the two adiabatic states (or the one state per fragment), which should span a similar space like the diabatic ones.Furthermore, a multi state treatment is not possible.Motivated by these shortcomings, the Multi-FED-FCD 30 has been developed.The Multi-FED-FCD method relies on previous works which address more than two states in the calculation and automatically assign states based on their eigenvalues, namely the 3-/4 state FCD approach 34 and the Multi state FCD approach. 35Building on this knowledge, the Multi-FED-FCD scheme i) takes many adiabatic states as input to avoid manual adiabatic state selection, and ii) automatically assigns LE-and CT states according to their eigenvalues by an efficient combination of the FCD and FED operators.In more detail, the Multi-FED-FCD process follows the subsequent steps: 1. Definition and calculation of a new operator in the adiabatic basis, which serves to distinguish between LE-and CT states: D = (∆q) 2 − (∆x) 2 with eigenvalues +1 (CT state) and -1 (LE state), 2. Diagonalization of D and sorting of eigenvectors according to eigenvalues to separate the subspaces: 3. Transformation of ∆q and ∆x into this new basis and diagonalization within subspaces to obtain the block-diagonal transformation matrix U 2 : ∆q 4. At this stage, the transformation represented by U 1 U 2 already transforms into a diabatic basis.However, after a second separation into LE1-, LE2-, CT1-, and CT2-subspaces according to the eigenvalues of the charge and excitation difference operators, a third transformation U 3 is applied which introduces the assumption that diabatic states of the same type should have zero electronic coupling.This transformation is indicated in Figure S1 of the supplementary information of ref 30 by the resulting ϵ-blocks which represent diagonal submatrices.Therefore, the couplings between LE states on the same pigment, and between CT states where the hole and particle are localized on the same fragments, are zero.This assumption, although frequently used (e.g., [33][34][35] ), has been questioned in the literature 36 for the case of multimers though.After transforming the Hamiltonian to the intermediary diabatic representation with U 1 and U 2 , the block-diagonal transformation matrix U 3 is constructed by rediagonalizing the Hamiltonian in the respective four diabatic subspaces.The resulting overall scheme is visualized nicely in Figure S1 of the supplementary information of ref, 30 noted already above, and the resulting adiabatic-to-diabatic transformation matrix is represented by Although the Multi-FED-FCD method removes the necessity of selecting adiabatic input states manually, there remains the question of how many adiabatic states are required for a sufficient description of the obtained diabatic states and elements of the Hamiltonian in the diabatic representation.A possible answer to this quality-related question is already provided in 35 with a convergence criterion of the quantity of interest (e.g., a certain CT-LE coupling value) with regard to the number of adiabatic input states.
To further assess the quality of the obtained diabatic states, we checked i) the eigenvalues resulting from the Multi-FED-FCD method, and ii) the character of the diabatic states.The latter is calculated by summing up the contributions from the transition density matrix in the atomic orbital basis, which lie on the respective fragments that characterize the kind of diabatic state of interest. 37Specifically, the quality of diabatic states is computed as with T being the transition density matrix of diabatic state i in the atomic orbital basis with orbitals µ and ν.Index A indicates that the sum runs over atomic orbitals belonging to atoms of fragment A. Thus, X AA represents the quality of a locally excited state of fragment A. Accordingly, computes the quality of a A + B − CT state.
The Multi-FED-FCD diabatization procedure 30 has been applied to obtain an initial set of diabatic states.
The adiabatic-to-diabatic transformation obtained from the diabatization has been also applied to compute the diabatic electronic and magnetic transition dipole moments.

Perturbation theory for inclusion of CT states: resonant case
We consider the special pair with local excited Q y states Because of the degeneracy, the excited state wavefunction will be delocalized and the adiabatic states read with energy E (ad) 1 with energy is the direct excitonic coupling between local excited Q y states obtained from the diabatization.Now, we switch on the coupling to CT states where, e.g., V 1k = V P D1 Q y ,kCT is the coupling between the Q y state of P D1 and the kth CT state, and h.c.denotes the hermitian conjugate.
In second-order perturbation theory, the following energies are obtained for the two adiabatic states and Taking into account that and that the coupling elements are real (V 1k = V k1 and V 2k = V k2 ), the difference between the energies of adiabatic states becomes The splitting of eigenenergies in a homodimer is two times the coupling and, hence, we obtain the following where is the direct coupling between local excited states (containing long-range Coulomb as well as non-CT-SR contributions) and the second term on the r.h.s.contains the superexchange-type coupling mediated by the coupling between local excited and CT states.

Perturbation theory for inclusion of CT states: non-resonant case
We consider the special pair with local excited Q y state of P D1 , |Ψ 1 ⟩ = |P D1 Q y ⟩, with energy E 1 and a high-energy local excited state of P D2 , |Ψ 2 ⟩ = |P D2 a⟩, with energy E 2 .The coupling between these local excited states and the CT states (eq 7) is treated in perturbation theory up to second-order, giving for the adiabatic states and Taking into account that the coupling matrix elements are real, that is, V 1k = V k1 and V 2k = V k2 we may express the above equations as and where we have introduced the effective couplings and Obviously, it holds that V eff 12 ̸ = V eff 21 .In order to restore the required symmetry, we determine an energy difference ∆E such that its inverse has minimal deviations from both 1/(E 1 − E k ) and 1/(E 2 − E k ), occurring in V eff 12 and V eff 21 , respectively.Hence, we have The solution for 1/∆E reads and the effective coupling is obtained as If one of the two states, say |Ψ 2 ⟩ has an energy E 2 that is not off-resonant enough from the CT state, we can only give a lower bound to V eff 12 by setting E 2 = E 1 in the above equation, resulting in The term on the r.h.s. is used as an estimate of the CT contribution to the excitonic coupling V (Q y ,a) P D1 P D2 between local excited Q y state of P D1 and the ath high-energy excited state of P D2 in eq 3 of the main text.

Diabatic Hamiltonian
The full Hamiltonian of the special pair as obtained from the diabatization is provided as an additional file "SP_Hamiltonian_diab.xlsx".The values for the CT states and the Q y states are reported in Table S4.The energies of the lowest CT states are found to be around 10000 cm −1 higher than the Q y states.However, relatively large couplings of up to 1127 cm −1 suggest a potential impact of the CT states on the Q y states.Indeed, the site energies after perturbation by CT states (eq 2 of the main text) have been computed to be 15919 cm −1 and 15922 cm −1 for P D1 and P D2 , respectively, as compared to 16079 cm −1 and 16078 cm −1 , respectively, without the perturbation (Table S5).The coupling between the Q y states of P D1 and P D2 is increased from 137 cm −1 to 235 cm −1 by the coupling to CT states (eq 3 of the main text with a = Q y ).Note, that the full Hamiltonian after perturbation by CT states is provided as an additional file ("SP_Hamiltonian_pert.xlsx").Furthermore, an analysis of the character of the first two adiabatic states according to eq 3-4 yields a composition of CT1: 2.6%, CT2: 2.5%, LE1: 47.9%, LE2: 46.9% for the first adiabatic state and a composition of CT1: 1.7%, CT2: 1.7%, LE1: 46.5%, LE2: 50.0% for the second adiabatic state.
In order to assess the quality of the obtained CT states (cf.Table S4), their character has been computed according to eq 4 as 95%, 93%, 82%, and 89% for the CT1 states, and 93%, 92%, 91%, and 88% for the CT2 states.Thus, although yielding diabatic states of high quality, the diabatization procedure cannot provide strictly separated diabatic states.The latter result justifies to include the effects of the quasi CT states perturbatively for the respective energies, couplings, and dipole moments of the Q y states despite strict CT states would be optically dark and, thus, should not affect the electric transition dipoles.However, as stated in the main text, we do not find the effects of the CT states to be the main cause for the observed non-conservativity in the CD-spectra.

Excitation energies
The Q y transition energies (site energies) obtained from quantum chemistry calculations on isolated monomers and (after diabatization) from dimers are further analyzed regarding the impact of the environment.The results are summarized in Table S5.Vacuum calculations are compared with calculations were the protein environment was included in the quantum chemical calculations by point charges (electrostatic embedding).In vacuum, the site energies obtained after diabatization are blue-shifted by about 60 cm −1 with respect to the values obtained for the isolated monomers.Including the perturbation by CT states (eq 2 of the main text) leads to a 150 cm −1 red shift.Hence, in vacuum the non-CT-SR effects and the CT coupling effects partially compensate each other.Including the electrostatic coupling to the protein environment adds a 50 cm −1 red shift to P D1 and a 15 cm −1 blue shift to P D2 , such that the site energies of the two pigments are in resonance.The overall effect of CT state coupling, non-CT-SR and long-range electrostatic effects on the site energies of the special pair are small.

Poisson-TrEsp couplings
The long-range couplings between special pair states have been computed using the Poisson-TrEsp method. 38,39This method takes into account the polarizability of the protein environment by that of a homogeneous dielectric with an average optical dielectric constant of ϵ opt = 2, estimated from an analysis of the oscillator strength of Chl a bound to photosystem I. 40 The required atomic partial transition charges for the special pair are obtained from quantum chemical calculations on the isolated P D1 and P D2 monomers.The transition densities between the ground state and the first 5 LE states were calculated using TDDFT with the ωB97X-D3 XC functional and a def2-SVP basis set, using Q-Chem. 41From a fit of the ESP of the transition densities, the respective atomic transition charges were obtained.The MEAD software tool ii has been used for the solution of the Poisson equation for the electrostatic potential of the transition charges used ii https://rtullmann.de/index.php?name=extended-mead P D1 P D2 (Q y , a): post-processed long-range excitonic couplings considering rescaling of electronic transition dipole moments as well as polarization effects in the environment;V diab P D1 Q y ,P D2 a : excitonic couplings obtained from diabatization (i.e., including non-CT-SR effects); V (Q y ,a) P D1 P D2 (qc): overall excitonic coupling after perturbation by CT-states (eq 3 of the main text); V (Q y ,a) P D1 P D2 (pp): overall post-processed excitonic coupling values obtained as V (Q y ,a) P D1 P D2 (qc)-V (LR, qc) (LR, qc) for the coupling calculations.The transition charges of the LE states of the special pair have been rescaled to correct for uncertainties in the quantum chemical transition dipole moments, based on experimental oscillator strengths of Chl a, as it has been conducted in ref. 42 On the one hand, the scaling factor for the Q y states is obtained by a comparison of the vacuum transition dipole moments of the monomers with the value of for this transition, based on an estimate by Reimers et al. 44 For the reaction field factor, introduced recently into the Poisson-TrEsp method, 45 we use the lower limit 1.0, since it provides a better fit of the experimental CD spectrum (Fig. S10).The long-range contributions to overall excitonic couplings in Eq. 3 of the main text, obtained with unscaled transition charges from the diabatization and subsequent perturbation theory, were then corrected by subtracting the vacuum coupling obtained from the unscaled transition charges and adding the above Poisson-TrEsp values obtained for rescaled transition charges and in the presence of a dielectric medium.This procedure yields coupling values between Q y and higher excited states which are summarized in Tables S6 and S7, where long-range contributions are provided together with the respective couplings after diabatization and perturbation by CT states (eq 3 of the main text).Whereas the excitonic coupling between Q y transitions is strongly enhanced by the coupling to CT states, the excitonic coupling between the 3rd LE state of one special pair pigment and the Q y transition of the other is barely affected by SR effects.In two cases, V (Q y ,4) P D1 P D2 and we find a sign inversion of the overall coupling by SR effects.(eq of the main text), middle left : same as top but for P D2 , bottom left:

Fluctuations of Hamiltonian parameters
(eq of the main text).
Table S8: Q y transiton energies of the special pair pigments as well as excitonic couplings between Q y transitions, obtained after diabatization and including a perturbation by CT states (eqs 2 and of the main text for site energies and couplings, respectively).The values obtained for the representative structure are compared with the mean values of the 110 MD snapshots.couplings are very close to the transition energies obtained for the representative structure (Table S8).In the calculation of spectra, the distribution functions are approximated by Gaussians with standard deviations obtained from the snapshots.Note, that the snapshot analysis has been performed based on 20 adiabatic input states, whereas 30 adiabatic input states have been used for the representative structure calculations.
Although the histograms in Figure S6 roughly indicate a Gaussian distribution, the possibility of having multiple possible conformations for the special pair cannot be eliminated completely with this relatively small sample size.The standard deviation in units of cm −1 of the energy fluctuations of the diabatic transition energies E diab P D1 a and E diab P D2 a of the first a = 1...5 LE states are obtained as.where the FWHM values (in cm −1 ) are given in brackets.Taking into account the coupling to CT states (main text, eq 2) leads to an increase of the variation to 92 cm −1 (217 cm −1 FWHM) and 94 cm −1 (221 cm −1 FWHM) for the Q y state of the P D1 and P D2 pigments, respectively.
The standard deviations in diabatic couplings V diab P D1 Q y ,P D2 a and V diab P D1 a,P D2 Q y between Q y state and the ath LE state of the other special pair pigment (FWHM in brackets, all in units of cm −1 ) are obtained as Taking into account the coupling via CT states (eq 3 of the main text) results in the couplings V Additionally, the explicit fluctuations of CT state energies have been computed to be ≈1550 cm −1 and the fluctuations of the coupling between CT states and LE states are in the range of 4-230 cm −1 .The latter fluctuations are about 5 times higher than the fluctuations of the direct LE couplings.In particular, variations of LE-LE couplings between Q y transitions are about 4 times higher when taking into account their indirect coupling via CT states, as seen also in the right panels of Fig. S6.
6 Computation of CD-spectra 6.1 Method

Embedding the Special Pair Hamiltonian into the Reaction Center Hamiltonian
The long-range excitonic couplings between special pair states and LE states of the other reaction center pigments have been computed using the Poisson-TrEsp method, 38,39 described above.The required atomic partial transition charges for the special pair are calculated based on the diabatic transition densities and are determined with Mulliken partitioning, whereas the transition charges of other reaction center pigments are taken from prior work. 1 The coupling values between as well as transition energies of excited states which do not involve the special pair pigments have been taken from ref. 1 Note that we have used different transition charges for the intra special pair couplings, described above, since in the latter case we wanted to avoid artefacts due to the fact that the diabatic transition densities are not perfectly localized.For the couplings between special pair pigments and the other RC pigments, these effects are negligible.
The system-field interaction is characterized by the electric transition dipole moments µ ma of LE states.
For the non-special pair pigments, these transition dipole moments refer to those of monomeric Chl a, Pheo a and Car, as described previously. 1 For the special pair pigments, we use the transition dipole moments of the diabatic states instead.For the Q y transition of the special pair pigments, in addition, the coupling to CT states is included (main text, eq 9).The diabatic electronic transition dipole moments have been scaled for the Q y -region and the Soret-band, respectively, as it has been conducted in ref. 42 The atomic transition charges for the Poisson-TrEsp calculations for couplings between special pair states and other RC pigments have been scaled analogously.The scaling factors for the transition dipole moments have been calculated to be 0.840 (Q y of P D1 ), 0.857 (Q y of P D2 ), 0.905 (higher excited LE states of P D1 )), and 0.912 (higher excited LE states of P D2 )).In order to keep the scaling procedure consistent, we keep only the first 5 LE states.This truncation has no significant effect on the results based on rotational strength calculations of the special pair.
Regarding the other reaction center pigments, the Q y -, B y -, B x , and N x+xy are considered for the 4 chlorophyll a and pheophytin a pigments.Additionally, four excited states of the β-carotene (Car) are taken into account, corresponding to the 0-0.0-1, 0-2, and 0-3 vibronic transitions between the electronic ground state and the first dipole-allowed excited state (S 2 ).
The coupling between Q y transitions in the special pair is found to be 193 cm −1 (see above and eq 4 of the main text).In the calculation of the CD spectra, we nevertheless use a value of 150 cm −1 that has been inferred from a fit of optical spectra. 46Indeed, we find that the non-conservativity of the CD spectrum obtained with this coupling is slightly closer to the experimental value.An assessment regarding the influence of this coupling value on the CD spectrum is provided in Fig. S8.Prior to the computation of optical spectra, the special pair structure has been embedded into the crystal structure using the Kabsch algorithm. 47 quantum chemical methods are known to yield too high excitation energies, they have been shifted a posteriori.The Q y transition energies of the special pair pigments P D1 and P D2 have been shifted to resemble the values from our earlier studies, 1,4 and the transition energies of the higher LE states are shifted by the same constant (950 cm −1 ).Finally, the site energy values for the reaction center pigments obtained from the literature 1,4 come with a tolerance of ±60 cm −1 .Therefore, a final manual fit of Q y transition energies (site energies) has been conducted for the final CD-spectrum using the full model.Please find a comparison between the refined site energies and those of the literature (that correspond to that of the present minimal model) in Table S9.The special pair excitation energies used for the calculation of spectra are for P D1 : 14975 cm −1 , 18934 cm −1 , 26690 cm −1 , 28242 cm −1 , 29251 cm −1 ; and for P D2 : 14975 cm −1 , 19215 cm −1 , 26357 cm −1 , 28402 cm −1 , 28954 cm −1 .Note, that the full reaction center Hamiltonian is provided as an additional excel file "RC_Hamiltonian.xlsx".

Theory of Circular Dichroism Spectra
The transition of the system from the ground state to a one-exciton state is caused by couplings to the fields of the electromagnetic wave with frequency ω and is included as H ex−rad into the overall Hamiltonian where the exciton-radiation interaction Hamiltonian is defined in the rotating wave approximation 48 as with where rm j is the coordinate of the j-th electron of pigment m, pm j = −ih∇ j is its momentum, m e is its mass, ) are the amplitudes of the incoming electromagnetic field, and ϵ is its polarization unit vector.The light propagates in the direction of n = k k .The rate constant for the transition to the Mth exciton state is given by Fermi's golden rule as with eigenvalues of H ex given by E M = hω M0 .The δ-distribution is replaced by a line shape function D M (ω) as a result of the exciton-vibrational couplings.For the calculation of optical spectra, the exciton-vibrational coupling as well as the coupling of the RC to external light fields has to be taken into account.As described previously, 1 the exciton-vibrational coupling gives rise to a line shape function of optical transitions that contains the excitation of vibrational sidebands and a lifetime broadening due to exciton relaxation.As before, we use a line shape function obtained with a partial ordering prescription (POP) cumulant expansion approach and a spectral density J(ω) that has been extracted from fluorescence line narrowing spectra of the B777 pigment-protein complex. 49The line shape function reads 50 where τ M describes the lifetime broadening of the exciton states, and ∆ω M a frequency shift as defined in ref. 50e function F MM (t) depends on the generalized spectral density The latter is commonly approximated to be of the form 49 with the Huang-Rhys factor S and the normalized function ( ∫︁ ∞ 0 dω J 0 (ω) = 1) with parameters s 1 = 0.8, s 2 = 0.5, hω 1 = 0.069 meV, and hω 2 = 0.24 meV. 49The shape of the spectral density has been found to be very similar for different pigment-protein complexes, whereas the Huang-Rhys factor may be different. 3The latter can be obtained from the temperature dependence of the linear absorbance, where S = 0.65 was determined for the PSII RC. 51 Note that the fluctuations of excitonic couplings and the correlations in site energy fluctuations can be neglected, based on an earlier normal mode analysis of the spectral density. 52,53 order to calculate the rate constants k 0→M , the matrix elements h m0 are determined using a multipole expansion 54 in the wave vector k.This approximation yields up to first order h m0 ≈ h m0 with the zeroth-order contribution h and the first-order term h (1) Using the relation the zeroth-order contribution becomes with i |0⟩ being the electric transition dipole moment of pigment m in so-called length gauge.The first-order contribution splits into two parts which are symmetric and antisymmetric with respect to the interchange of k and ϵ as h (1) Using the vector identity and the definition of the magnetic transition dipole moment the first contribution on the r.h.s. of (33) becomes The second contribution h results from the electric quadrupole moment and reads m , (37)   where in the last step the definition of the quadrupole operator Q ˆ(αβ) where R m denotes the center position of pigment m and r ′(m) i the relative coordinate with respect to this center.Using this transformation, the magnetic transition dipole becomes with the intrinsic magnetic transition dipole m The transition quadrupole moment is transformed as with the intrinsic transition quadrupole moment Q Inserting the expressions for the magnetic transition dipole as well as the electric transition quadrupole into (36) and (37) yields, using ( 34) It has been shown in ref 48 that the contribution of the intrinsic transition quadrupole vanishes for CD spectra of isotropic samples.As we only investigate spectra of the latter kind, that contribution is neglected already at the present stage, and the overall first-order contribution is obtained from eqs (33), (41), and (42) as Using the multipole extension, mentioned above, the rate constant in eq (25) is approximated as with the zeroth-order part The first-order contribution to the rate constant is given by Inserting ( 32) and ( 43) yields Isotropic circular dichroism The rate constant k 0→M is related to the absorption coefficient α(ω) by 55 as where n mol denotes the density of molecular aggregates in a homogeneous sample.The zeroth-order contribution α (0) (ω) and first-order contribution α (1) (ω) result from ( 45) and ( 47), respectively.The CD signal is defined as the difference between left and right circularly polarized light.Without loss of generality, the light is assumed to propagate in z-direction (n = (0, 0, 1)), resulting in polarization vectors ϵ r = 1 √ 2 (1, −i, 0) and (1, i, 0).The zeroth-order terms α (0) l/r (ω) in ( 48) are identical because ϵ r = ϵ * l .Therefore, their contribution vanishes for the CD signal.The latter is solely obtained from the first-order contributions by where again the vector identity (34) has been used.It can be seen from (49), that the CD signal splits into an excitonic contribution and an intrinsic contribution CD(ω) = CD(ω) (exc) + CD(ω) (intr) .For an isotropic sample, a rotational average over randomly oriented complexes gives the isotropic signal 48 CD(ω) where (ϵ * l × ϵ l ) = in, as well as ∑ m,n ... = 1 2 (∑ m,n ... + ∑ m,n ...), interchanging of variable names in the second sum, and R mn = R m − R n have been used.The orientational average has been carried out using Euler angles 56 as conducted in. 48Similarly, the intrinsic isotropic contribution reads CD(ω) Please note that eqs 50 and 51 reduce to expressions derived earlier 48 if the local excitation energies E m and E n are approximated by the photon energy hω.This approximation becomes questionable, when the coupling between transitions of very different excitation energies are included in the theory.We have, therefore, repeated the whole derivation without using this approximation.Note also, that magnetic transition dipole moments are purely imaginary, that is, }.The CD signal can also be obtained using an alternative theory which is commonly used in the literature and known as the Rosenfeld equation. 57,58ith (51) represent the most accurate result, approximating the diabatic energies by the photon energy results in well-known expressions for the excitonic part.Although the commonly used Rosenfeld equation, ( 53) and (54), includes contributions originating from the intrinsic magnetic transition dipole moment, the resulting CD signal is not origin-independent due to the bespoke approximation.Moreover, if also higher energy transitions are taken into account in the calculation of CD signals, the respective energies can be substantially different from the photon energy related to the investigated region of the spectrum.Please note, that in our calculations the high-energy LE states are only included to study their influence on the low-energy (Q y ) region of the spectra.An open question was whether the non-conservativity of the low-energy region of the experimental CD spectrum is mainly caused by the excitonic coupling between high-energy and lowenergy excited states of the pigments that redistributed rotational strength, or by the intrinsic (non-excitonic) rotational strength of the special pair.The present calculations show that the latter effect dominates for the present system.

Theory of OD and LD spectra
The expressions for the linear absorption (OD) and linear dichroism (LD) spectra read and respectively, with the electric transition dipole moment of the Mth exciton state and the angle θ M between this transition dipole moment and the membrane normal.
Note that we assume an isotropic sample in the case of the CD and the OD spectrum.In case of LD, the complexes have been oriented in a plane with isotropic orientation with respect to the normal on that plane.
This plane was assumed to be parallel to the membrane, in which PSII is embedded.

CD spectra for alternative parameters and theories
In the following, a comparison of CD spectra calculated for different parameters and models is presented.
First, the spectrum of the full model in Fig. 2 of the main text is compared to the spectrum obtained with the same parameters but using the Rosenfeld equation (eqs 57 and 58) for the calculation of the signal (Figure S7, right).Whereas the long-wavelength positive peak obtained in the two theories agrees well, the relative intensity of the two negative peaks is somewhat different.The non-conservativity ratio is very similar (3.6 ± 0.1).
In the left part of Figure S7, the spectra of the full model obtained with the refined site energies (Fig. 2 of the main text) are compared to those obtained with the same model but a site energy set from the literature. 4e former have been fitted based on a tolerance of 60 cm −1 with respect to the latter.Numerical values of these site energy sets are given in Table S9.The CD spectra obtained with the original site energy set contain a slightly narrower main positive band and a somewhat red-shifted main negative band, as compared to the calculations with the refined site energies.
If instead of the 150 cm −1 coupling between the Q y transitions of the special pair pigments the directly calculated value of 193 cm −1 (eq 4 of the main text) is used, a larger intensity and a somewhat red-shifted position of the positive peak in the CD spectrum are obtained (Figure S8).Note, that the ratio of non-  S9.
Furthermore, the effects of the short-range contributions on the magnitude of the electric transition dipole moment were investigated in Fig. S8.The SR contributions are found to slightly enhance the dipole strengths of the Q y transition of P D1 and P D2 by about 4%.However, this enhancement has practically no effect on the CD spectrum.Note that the rotation of the transition dipole moments by the SR effects is taken into account.
The model by Lindorfer et al. 1 took into account the coupling between the Q y and high-energy transitions of the pigments and postulated increased excitonic couplings mediated by CT states in the special pair.Here, we find that this mediation is much weaker than assumed (Tables S6 and S7).In order to further connect to our previous calculations, 1 we replaced the present diabatic transition energies and dipole moments with those obtained for the isolated monomers.Thereby, the effects of higher excited states are included, but short-range effects are omitted, except for the increased special pair coupling, which was kept at 150 cm −1 .The ratio of non-conservativity drops from 3.4 to 1.0, with the respective spectra shown in Figure S9.In addition to assuming monomer quantities, the model proposed in ref 1 (i) does not consider intrinsic CD contributions, and (ii) approximates transition energies by the photon energy of the incoming light hω.Introducing those approximations on top of the considered monomer values yields ratios of non-conservativity of 1.8 (i) and 1.4 (i+ii), respectively.Based on the original model and data of ref, 1 we computed a non-conservativity ratio of 2.0 (without increase of coupling of Q y states to high-energy LE states via CT states).Remaining discrepancies which could explain the prevailing difference between the non-conservativity ratio of 1.4 and the ratio of 2.0 concern (i) different quantum chemical methods for the computation of monomer excited states, and (ii) fixing the direction of electric transition dipole moments to predefined directions rather than working with quantities obtained from quantum chemistry computations.Note, that the spectra are particularly sensitive to the direction of transition dipole moments (cf.main text).
Finally, we discuss the influence of the reaction field corrections, introduced recently into the Poisson-TrEsp method for the calculation of excitonic couplings. 45So far, we have used the lower limit f rf = 1.0.We find that for the present system, a larger reaction field factor leads to less agreement between calculated and measured spectra.The CD spectrum obtained for the upper limit f rf = 1.44 are compared in Figure S10 to the experimental spectrum and the spectrum obtained for f rf = 1.0 (same as in Fig. 2, of the main text, full model).The spectra obtained for f rf = 1.44 agree significantly less with the experiment than those   obtained for f rf = 1.0.Whereas the redshift of the spectrum calculated for f rf = 1.44 could be corrected by a blue-shift of all site energies, there are remaining discrepancies concerning the shape.We note that in a recent application of Poisson-TrEsp to Fenna-Matthews-Olson-protein-reaction center supercomplexes, also a rather low reaction field correction factor f rf = 1.15 has been inferred.A discussion of the uncertainties related to an independent estimation of this factor is given in ref. 45

Construction of a minimal model
In order to extract the essential parts of our full model, which includes 5 LE states per RC pigment, we investigated several approximations.In case of the special pair, these 5 LE states correspond to the 5 lowest diabatic states of P D1 and P D2 , where the one with the lowest energy (Q y ) is perturbed by superexchange-type coupling to CT states at higher energies.In a first attempt, we treated all high-energy LE states, that is, all LE states except for the two lowest (Q y and Q x ) for the Chls and Pheos and all excited LE states of the Cars in perturbation theory.Similar to the perturbation theory with respect to the coupling to CT states, as provided in the main text and in Section 6.1.1,we obtain perturbed energies , where a = Q y for all Chls and Pheos, except the special pair, and a = Q x , Q y for P D1 and P D2 .Here, V (0) ma,kc is the excitonic coupling between the ath excited state of site m and the cth high-energy LE state of site k.
In the case of the special pair, the energies (E (0) P D1 Q y and E (0) P D2 Q y ), couplings (V (0) P D1 a,P D2 c ) and transition dipole moments involving the Q y states contain already a perturbation by CT states, as described in the main text (eqs 2, 3, and 9 and 11, respectively).Note, that in case of the energies, these effects enter only indirectly since the final energies of LE states of the full model are shifted by a constant obtained from a fit of the experiment, as described in Section 6.1.1.coordinates, (ii) using long-range coupling values except for the intra special pair coupling, (iii) omitting the Q x -states, and (iv) omitting the effects of higher LE-states.All these approximations result in a negligible change of the spectra and ratios of non-conservativity.The latter are smaller than 0.1 and, therefore, below the tolerance limit.Using the unfitted site energies as originally proposed in ref 4 even leads to an increase of the non-conservativity from 3.1 to 3.3 compared to the use of the site energies from the full model calculations (Table S9).Therefore, we take the originally proposed energies in our minimal model.Omitting the influence of CT-states but still accounting for electron-electron exchange, results in a reduction of the non-conservativity ratio from 3.3 to 2.8, whereas the use of transition dipoles from monomer calculations, i.e., neglecting all short-range contributions, leads to a ratio of 0.9 as described in the main text.
In Figure S13, the OD, LD, and CD spectra are shown for both sets of site energies using the full model and in Figure S14 these spectra are shown again for the two sets of site energies, but using the minimal model, derived above.  and fitted site energies for the minimal model calculations compared to experimental spectra from 61 Table S10: Hamiltonian parameters for minimal model, excerpt of full reaction center Hamiltonian as well as implicit high-energy (h-e) model (minimal/full/implicit h-e); note that the Q x -states are not included in the minimal model.The original site energies (Table S9) are taken here in the minimal and full model and are used as the umperturbed energies E (0) ma in the h-e effective model (eq 63).The excitonic coupling 150 cm −1 between Q y transitions of P D1 and P D2 in the minimal and full model has been adopted from the directly calculated value (eq 4 of the main text) to fit the experimental CD spectrum (see also Fig. S8).Table S11: Electric transition dipole moments of the five lowest excited states (a = 1 . . .5) of the special pair pigments (m = P D1 , P D2 ) obtained by calculations on isolated monomers (µ mon ma ) or dimer calculations and diabatization (µ m , eq 8 of the main text).For the Q y states (a = 1), we also included the perturbation by CT states (numbers in brackets, eq 9 of the main text).Angles ϕ and θ are defined in Fig. 4 of the main text.The quality of LE states after diabatization is defined in section.5  originates from intra-dimer non-CT-SR effects, which are only present in the full dimer calculation.BSSE and polarization effects by the protein environment have only very minor effects on the transition dipole moments.
To further demonstrate this effect, we performed a geometry scan.Starting with our representative structure, we increased the distance between the two pigments by moving the P D2 pigment along the normal on the plane of the P D1 pigment.In this way, we maintain the pigment's relative orientation while increasing their distance.The scan results (Fig. S18) clearly show that as the inter-pigment distance is increased, the diabatic transition dipole moments are rotating toward their orientations obtained for the isolated monomers (Fig. S18 A and B, and Table S13).In Fig. S18 C, the magnitude of the two lowest (M = 1, 2) adiabatic transition dipole moments are shown as a function of distance increase ∆r.For ∆r = 0 (native distance), the transition dipole moment of the lowest adiabatic state (M = 1) is larger than that of M = 2, reflecting the "in-line" type geometry of transition dipole moments.At ∆r around 5 Å, the magnitudes of the adiabatic transition dipole moments become equal to that of the adiabatic transition dipole moments.In this case, the excitonic coupling becomes zero.For larger separations ∆r, the dipole geometry changes towards "sandwich"-type and the higher-energy adiabatic state (M = 2) gets the largest dipole moment magnitude.With increasing distance, the excitonic coupling decreases and the dipole strengths of adiabatic states slowly converge towards that of the diabatic states.
Table S13: Q y electric and intrinsic magnetic transition dipole moments of the special pair pigments, with ("Protein") and without ("Gas Phase") the protein environment.Angles Θ and Φ are defined in Fig. 4 of the main text.The values for the protein environment are identical to those of 1 µ mQ y as defined in Eq. 9 of the main text 2 m mQ y as defined in Eq. 11 of the main text 3 defined in Fig. 4 in the main text.
Quantum chemical calculations of excited states of special pair, electrostatic embedding of pigment/protein environment Post−processing QM (special pair)/MM optimization of MD snapshots Energies, couplings, transition dipole moments,

Figure S3 :
Figure S3: Root mean square deviation (RMSD) of the protein backbone, with respect to crystal structure, during the MD simulation.For further exploration of configuration space, ten snapshots with equal time spaces (marked by dashed lines) were taken from the last 200 ns of the trajectory.

Figure S4 :
Figure S4: Root mean square deviation (RMSD) of the protein backbone, with respect to the crystal structure, of ten additional MD simulations.Eleven snapshots from each trajectory (marked by dashed lines) were taken for calculating the static disorder.
structure and static disorder In order to choose a representative structure of our PSII model, we selected the snapshot showing the smallest RMSD deviation from the average backbone geometry.For generating static disorder, we further explored the phase space in the MD simulations.Ten snapshots, with equal time spaces between them, were taken out of the last 200 ns of the production run.

4. 58 D
for the Q y transition of chlorophyll a obtained by Knox & Spring43 from an analysis of experimental oscillator strengths measured in different solvents.The scaling factor for the higher excited states on the other hand results from a comparison of calculated transition dipole strengths relative to that calculated for the Q y transition and corresponding areas of the absorption spectrum.42For the present quantum chemical calculations we obtain scaling factors that lead to calibrated vacuum transition dipole moments of 4.58 D (1.LE), 6.99 D (3.LE), 6.93 D 4.LE) and 4.38 D (5.LE).In the present work we include the Q x -transition (2.LE) of the Chls in the special pair.An experimental vacuum transition dipole moment of 1.47 D has been assumed As outlined in the main text, we monitored the change in the parameters of the special pair Hamiltonian along 110 MD snapshots, that were each geometry-optimized and analyzed by our quantum chemical/diabatization procedure.Examples of distribution functions obtained from these analyzes are shown in Fig.S6, where the Q y transition energies of P D1 and P D2 (eq 2 of the main text) as well as the coupling V (Q y ,Q y ) P D1 P D2 between Q y transitions, obtained after diabatization with and without including the perturbation by CT states are investigated.The mean values of the transition energies and the main text), the standard deviations of which (FWHM in brackets, all in units of cm−1 β has been used and r ˆ(m) i,α denotes the α-th Cartesian coordinate of r(m)

Figure S7 :
Figure S7: Comparison of spectra calculated in the full model using Q y energies refined in the present work (same as in Fig. 2 of the main text) with those obtained by using Q y energies from our earlier work 4 left panel) or with spectra calculated in the Rosenfeld model (right panel), using the refined site energies.Numeric values of site energies are given in TableS9.

Figure S8 :
Figure S8: Comparison of full model spectrum (Fig. 2 of the main text) to spectra obtained (i) with altered Q y − Q y coupling (150 cm −1 to 193 cm −1 ), and (ii) with monomer electric transition dipole strengths where, in contrast to the full model, short-range effects on the magnitude of dipoles are omitted.Note, that short-range effects regarding the direction of dipoles are still considered.

Figure S10 :
Figure S10: Comparison of final CD-spectrum with additional reaction field corrections as proposed by Friedl et al.; 45 experimental spectrum is from Hall et al.60

Figure S11 :
Figure S11: CD spectra obtained in the minimal model with (green line) and without (blue line) including the Q x transition of the P D1 and P D2 pigments.

Figure S12 :
Figure S12: CD spectra obtained in the full model, using individual Hamiltonian parameter fluctuations (blue) for the special pair and effective site energy fluctuations of the Q y energies of P D1 and P D2 (green).

Figure
FigureS13: CD (left), OD (middle) and LD (right) spectra obtained with original site energies from4 and fitted site energies for the full model calculations compared to experimental spectra from61

Figure
FigureS14: CD (left), OD (middle) and LD (right) spectra obtained with original site energies from4 and fitted site energies for the minimal model calculations compared to experimental spectra from61

Figure S15 :
Figure S15: Electric transition dipole moments of 2.-5.LE-state of P D1 (blue) relative to the reference vector from the pigment center to the N B -atom (red); top: µ mon ma from monomer calculations, bottom: µ (a) m from dimer calculations after diabatization.

Figure S16 :
Figure S16: Intrinsic magnetic transition dipole moments of 2.-5.LE-state of P D1 (blue) relative to the reference pigment plane normal vector (red); top: m mon ma from monomer calculations, bottom: m (a) m from dimer calculations after diabatization.

Figure S17 :
Figure S17: Upper part: Gas-phase calculated Q y monomeric electric (red) and magnetic (orange) transition dipole moments are compared with the respective transition dipole moments resulting from dimer calculations and subsequent diabatization and perturbation by CT states (eqs 9 and 11, respectively, of the main text) shown in blue and cyan, respectively.The adiabatic electric transition dipole moments for the M = 1 (yellow) and M = 2 (gold) states are shown in yellow and gold, respectively, and the respective adiabatic magnetic transition dipole moments are shown in purple and pink.Lower part: Same as upper part, but including the protein environment as point charges (electrostatic embedding).

Figure S18 :
Figure S18: Orientation and magnitude of electric transition dipole moments of the special pair, as a function of distance increase ∆r.Starting with the native structure (∆r = 0), the distance was increased by moving the P D2 pigment along the normal on the plane of the P D1 pigment.A) Angle Φ of the diabatic electric transition dipole moments of P D1 and P D2 .B) Angle Θ of the diabatic electric transition dipole moments of P D1 and P D2 .(For a definition of angles see Fig. 4 of the main text.C) The magnitudes of electric transition dipole moments of the adiabatic states M = 1, 2 are compared to that of the diabatic states of P D1 and P D2 .

Table S4 :
Excerpt of diabatic Hamiltonian comprising the four lowest CT states P −

Table S5 :
Q y transition energies (in cm −1 ) as obtained from quantum chemistry calculations on (i) isolated monomers, (ii) dimers, after diabatization, (iii) dimers, after diabatization and perturbation by CT-states (eq 2 of the main text) in vacuum and by including the protein environment as point charges (for (ii) and (iii)).

Table S6 :
Excitonic couplings between Q y state of P D1 and five lowest excited states of P D2 (a = 1, ..., 5) ; D2 (Q y , a): long-range excitonic coupling based on transition charges of quantum chemical calculations; V (LR, pp)

Table S7 :
same as TableS6but between Q y state of P D2 and higher excited states of P D1

Table S9 :
4 y transition energies (in units of cm −1 ) of the present full model are compared to those of our earlier work.4Pleasenote that the latter values are used also in the minimal model, introduced in the present work.

Table S12 :
.1.Pigment m LE state a Quality [%]|µ| [A.U.] Same as in TableS11but for magnetic transition dipole moments m mon ma (isolated monomers) and m (eqs 10 and 11 of the main text). m Table 2 of the main text.In "Gas Phase with Ghost Orb." and in "Gas Phase with Point Charges" the BSSE and the Coulomb coupling, respectively, were investigated (see text).