Photosystem II (PSII) is the only biological system capable of splitting water to molecular oxygen. Its reaction center (RC) is responsible for the primary charge separation that drives the water oxidation reaction. In this work, we revisit the spectroscopic properties of the PSII RC using the complex time-dependent Redfield (ctR) theory for optical lineshapes [A. Gelzinis et al., J. Chem. Phys. 142, 154107 (2015)]. We obtain the PSII RC model parameters (site energies, disorder, and reorganization energies) from the fits of several spectra and then further validate the model by calculating additional independent spectra. We obtain good to excellent agreement between theory and calculations. We find that overall our model is similar to some of the previous asymmetric exciton models of the PSII RC. On the other hand, our model displays differences from previous work based on the modified Redfield theory. We extend the ctR theory to describe the Stark spectrum and use its fit to obtain the parameters of a single charge transfer state included in our model. Our results suggest that ChlD1+PheoD1 is most likely the primary charge transfer state, but that the Stark spectrum of the PSII RC is probably also influenced by other states.

Photosynthesis begins with the absorption of solar light by specific pigments, (bacterio)chlorophylls and carotenoids, that are bound to light-harvesting proteins. Subsequently the energy is transferred to reaction centers, where the excitation energy triggers a fast charge separation, through which it is transduced into chemical potential. Oxygenic photosynthesis in plants, cyanobacteria, and algae begins in photosystem II (PSII), which is the only known biological complex with strong enough redox potential to split water. PSII is a dimeric supercomplex, consisting of major and minor light harvesting complexes and the reaction center (RC). Excitation within the light-harvesting antennae of PSII is transferred to the RC, leading to the oxidation of water, the reduction of plastoquinone, and the formation of a proton gradient across the thylakoid membrane. The PSII is thus the main producer of oxygen on Earth.1,2

The PSII RC itself contains six chlorophyll (Chl) and two pheophytin (Pheo) molecules arranged into two almost symmetric branches, corresponding to the D1 and D2 proteins.3 The pigment arrangement in the PSII RC is shown in Fig. 1. Two closely situated Chls, PD1 and PD2, are known as the “special pair” (from the analogous attribution in the bacterial RC). Each branch also contains an accessory chlorophyll ChlD1/D2 and a pheophytin PheoD1/D2. Two peripheral chlorophylls ChlzD1/D2 are about 25 Å away from the other six RC pigments, and they do not participate in the primary charge separation. In vivo, the electron transfer proceeds through the D1 branch to PheoD1 and then further to plastoquinone QA and finally to QB.2 The most often studied PSII RC preparation (the so-called D1D2cytb559 complex), however, lacks the oxygen evolving complex, the redox active tyrosine TyrZ and the quinone acceptors,4,5 thus allowing the electrons to accumulate on PheoD1. Meanwhile, the hole remains mostly localized on PD1.

FIG. 1.

The structural arrangement of pigments and other cofactors in the PSII RC. Data taken from the 3WU2 PDB structure.3 The tails of chlorophylls, pheophytins, and plastoquinones are truncated for clarity. Figure made with VMD software.32 

FIG. 1.

The structural arrangement of pigments and other cofactors in the PSII RC. Data taken from the 3WU2 PDB structure.3 The tails of chlorophylls, pheophytins, and plastoquinones are truncated for clarity. Figure made with VMD software.32 

Close modal

The energy and charge transfer processes in the isolated PSII RC, the PSII core (which includes the CP43 and CP47 antenna complexes), or the whole PSII have been extensively studied over the last few decades, as reviewed in Refs. 1 and 6–12. Nonetheless, the pathways of charge separation and the identities of the participating states are still under debate. Many studies claim that the primary electron donor is the ChlD1 pigment; thus, the primary charge transfer (CT) state is ChlD1+PheoD1.13–17 In addition, a series of studies18–20 have highlighted the importance of the special pair CT state PD2+PD1, which has been proposed to be strongly mixed with the exciton states, providing an alternative charge separation pathway. These suggestions have to be tested against an increasing amount of experimental data using appropriate theoretical modeling.

Most of the theoretical work on the PSII RC is based on the exciton concept.21 The first exciton model that captured some features of the PSII RC spectroscopic data was the so-called multimer model.22 It was formulated by taking the structural data of the bacterial reaction center into account and additionally assuming equal transition energies for all the core pigments. Later, as the PSII structural data became available with increasing resolution,3,23–25 the theoretical models were refined,16,18,20,26–30 and many spectroscopic properties of the PSII RC were successfully described. Some models have included one or more CT states phenomenologically, treating them as additional exciton states.18,20,26,28 Only a few studies have accounted for the relevant properties of electrons and holes.29,31

Although there have been attempts to obtain the PSII RC site energies directly from the structural data,30,33,34 the majority of the theoretical models derive their parameters from fits to the optical spectra. Since exact calculations are almost always too computationally expensive, approximate methods have to be applied to calculate optical lineshapes and dynamics. We have recently analyzed various approaches to calculating absorption lineshapes35 and have found that the complex time dependent Redfield (ctR) theory gives results very close to the exact ones. On the other hand, the modified Redfield theory, used in Refs. 18, 20, 26, and 28 and in our previous work on the PSII RC29, was found to perform significantly worse. The Markovian Redfield approach, used in Refs. 16, 27, and 30, showed improvements but did not perform as well as the ctR theory. These findings raise some doubts about the optimization of the parameters of the previous models of the PSII RC16,18,20,26–29 and motivate our use of ctR theory for describing its spectroscopic properties.

The charge transfer states, being optically dark, are insensitive to most spectroscopic measurements like absorption, linear, or circular dichroism. Nonetheless, due to their strong static dipole moments, they are sensitive to the presence of static electric fields and thus can influence the results of Stark spectroscopy. Stark experiments have been performed on the PSII RC,36 but the results were not easy to understand without a detailed theoretical description. Previously, the PSII RC Stark spectrum was fitted using the modified Redfield theory.18 Even though significant progress in understanding the observed features was achieved, a check of the obtained conclusions, especially the parameters of the CT states, using a more accurate lineshape theory is definitely called for.

In this work, we revisit the spectroscopic properties of the PSII RC using the ctR theory.35 We construct a model that includes a single CT state together with molecular excitations. We obtain model parameters by fitting a number of linear optical spectra. We also extend the ctR theory to calculate the Stark spectrum, which is then used to tune the parameters of the CT state. Due to the more advanced lineshape theory, in our current model, we improve upon the parameters obtained in our previous work based on the modified Redfield theory.29 

To describe the molecular system, we use the tight-binding Hamiltonian,29,31 which is an extension of the Frenkel exciton model21 that enables the inclusion of CT states. This method is similar to other recent descriptions of Frenkel excitons and CT states in molecular systems.37–39 In Secs. III and IV, this model will be applied to the PSII RC, but its formulation is independent of any particular molecular system.

We consider the ground and the single-excited state manifolds that are sufficient to calculate linear spectra. The global ground state is denoted by g, and its energy is set to zero. In the single-excited state manifold, we have the molecular (pigment) excitation states n*, with excitation localized on a single site n, and the charge transfer states m+n, which are not accounted for in the Frenkel exciton model. To describe the charge transfer states, we introduce separate creation and annihilation operators for electrons and holes. Thus e^ne^n creates (annihilates) an electron on site n, and h^nh^n creates (annihilates) a hole on site n. Then molecular excitation (ME) states are described by e^nh^n|g=|n* and CT states by e^nh^m|g=|m+n. The electron and hole creation and annihilation operators satisfy the Fermi commutation relations.29,31

The total Hamiltonian is composed of Hamiltonians of the system, the bath, and their interaction,

Ĥ=ĤS+ĤB+ĤSB.
(1)

The electronic degrees of freedom compose the system that is described by the tight-binding Hamiltonian

ĤS=mNnNtmnee^me^n+mNnNtmnhh^mh^nmNnNVmnehe^mh^nh^ne^m+mNnmnJmne^mh^mh^ne^n.
(2)

Here tmnetmnh is the electron (hole) hopping rate between LUMO (HOMO) orbitals of different pigments. Diagonal elements tmme/h correspond to the orbital energies. N is the number of molecules, and the total number of states in the single excited state manifold is Ntot = N2. Jmn is the resonance coupling between the molecular excitations on sites m and n. The term Vmneh corresponds to the electron-hole Coulomb attraction energy between sites m and n. Note that electron-electron and hole-hole repulsion terms are not needed to describe the single excitation manifold. For this Hamiltonian, the energies of the excited states are

n*|ĤS|n*=εn*=tnne+tnnhVnneh,
(3)
m+n|ĤS|m+n=εm+n=tnne+tmmhVnmeh.
(4)

Here εn* is the energy of the ME state n* (site energy) and εm+n is the energy of the CT state m+n. Note that for calculations only the energies and not their constituent terms need to be specified. The couplings between different excited states are

n*|ĤS|m*=Jnm,
(5)
n*|ĤS|k+l=tnleδnk+tnkhδnl,
(6)
m+n|ĤS|k+l=tnleδmk+tmkhδnl.
(7)

Thus, couplings involving the CT states are related, which is not the case if the CT states are treated as additional molecular excitation states.18,20,26

All vibrational degrees of freedom of the molecules themselves and their environment compose the bath and are represented by an infinite set of harmonic oscillators with the Hamiltonian,

ĤB=jωj2p^j2+x^j2.
(8)

Here ωj is the frequency, p^j is the dimensionless momentum, and x^j is the dimensionless coordinate of the jth bath mode. We set =1 throughout the paper. The system-bath coupling is assumed to be linear in the bath coordinates

ĤSB=nNjωjdnn,jehx^je^nh^nh^ne^nmNnmnjωjdmn,jehx^je^mh^nh^ne^m.
(9)

Here the two terms are separated for future convenience. The parameter dnn,jehdmn,jeh is the dimensionless displacement of the equilibrium configuration of the jth bath mode between the ground state and the excited state of the nth molecule (n+m CT state). This formulation differs from the previous applications of the tight-binding model, where the bath was coupled to electron and hole orbitals and not to the excited states themselves.29,31 The present approach is more natural when describing excitations in a molecular aggregate and not in a solid state system.

We introduce bath operators F^n=jωjdnn,jehx^j and F^mn=jωjdmn,jehx^j, enabling us to rewrite the system-bath coupling Hamiltonian as

ĤSB=nNF^ne^nh^nh^ne^n+mNnmnF^mne^mh^nh^ne^m.
(10)

For simplicity, we subsequently denote both the ME and the CT states by letters with tildes ã,b̃, Then summation over all single excited states is expressed as

ãNtotzã=nNzn+mNnmnzmn.
(11)

Here zã denotes any quantity related to the state ã. The reorganization energies of the excited states are defined as λã=jωjdã,jeh2/2.

For a harmonic bath, the key quantities that describe the bath effects are the correlation functions,40–42 

Cãb̃τ2τ1=TrBF^ãIτ2F^b̃Iτ1ρ^Beq.
(12)

Here

F^ãIt=eiĤBtF^ãeiĤBt
(13)

defines the operator in the interaction representation. The equilibrium bath density operator is

ρ^Beq=expβĤBTrBexpβĤB
(14)

with β=kBT1 being the inverse temperature in energy units, and kB is the Boltzmann constant. The spectral densities

Cãb̃ω=πjωj2dã,jehdb̃,jeh2δωωjδω+ωj
(15)

describe the bath in the frequency domain. They show how strongly the system states are coupled to the specific frequency bath modes. From the spectral densities, the bath correlation functions can be calculated as

Cãb̃t=1πdω11eβωeiωtCãb̃ω.
(16)

The reorganization energies are also related to the spectral densities

λã=1π0dωCããωω.
(17)

Diagonalization of the system Hamiltonian [Eq. (2)] leads to the introduction of the generalized exciton basis (denoted as such, because it includes not only Frenkel but also CT excitons),

|e=ANtotψe,ã|ã
(18)

with ψe,ã=e|Û|ã being the elements of transformation matrix Û, which satisfies

ĤSEB=ÛĤSSBÛ.
(19)

Here the superscript EB (SB) means the generalized exciton (site) basis. We assume that the matrix elements ψe,ã are real. In the generalized exciton basis, the bath correlation functions are

Ce1e2,e3e4τ2τ1=TrBF^e1e2Iτ2F^e3e4Iτ1ρ^Beq=ãNtotb̃Ntotψe1,ãψe2,ãψe3,b̃ψe4,b̃Cãb̃τ2τ1.
(20)

with F^e1e2=e1|ĤSB|e2. The lineshape functions, required to describe the optical lineshapes (see Sec. II B), are defined as

ge1e2,e3e4t=0tdτ20τ2dτ1Ce1e2,e3e4τ2τ1.
(21)

When static disorder is present in the system, the energies of the states are different in each particular case of the ensemble. Thus in the system Hamiltonian [Eq. (2)], one needs to change VnnehVnneh+δεn* and VmnehVmneh+δεn+m, with δεn* and δεn+m being random numbers from a Gaussian distribution with zero mean and some chosen standard deviations σn* and σn+m. All observable quantities then have to be averaged over disorder, i.e., X¯=Xdis.

1. Linear absorption

For calculations of the linear absorption, one has to include the interaction with a time-dependent electric field in the Hamiltonian, which induces the transitions between the ground and excited states,

ĤSF=μ^Et=nNμnoEte^nh^n+h^ne^n.
(22)

Here μ^ is the electric dipole operator, μn is the transition dipole vector for pigment n, and o is the polarization vector of the electric field. For now we disregarded the permanent dipole moments (the diagonal elements of the dipole operator). Note that the CT states are dark and carry no transition dipole strength. We also define the scalar dipole moment operator μ^=μ^o for later convenience.

Assuming linearly polarized incoming light, the linear absorption spectrum is calculated as a one-sided Fourier transform of the linear response function,35,40,42–44

Aω=ωRe0dteiωtR1tdis.
(23)

The linear response function is formally a tensor

R1t;o2,o1=TreiĤtμ^o2eiĤtμ^o1ρ^ggor.
(24)

Here ρ^gg=|gg|ρ^Beq is the ground state equilibrium density operator, o1 and o2 denote the polarizations of the incoming light and registered signal, and or denotes orientational averaging, which has to be performed for randomly oriented systems. In this expression, we have assumed the rotating wave approximation (RWA). For isotropic systems, the tensorial nature of the response function can be neglected and scalar response function can be defined as R1t=R1t;o,o, with o being an arbitrary polarization vector.

In the generalized exciton basis, the response function is expanded as

R1t=e1,e2Ntotμge2μe1gorTrBeiĤBte2|eiĤt|e1ρ^Beq.
(25)

Assuming the secular approximation (considering only terms e1 = e2, thus no coherence transfer) and using the cumulant expansion40 and the ctR theory for the off-diagonal fluctuations, we obtain35 

R1t=13e1Ntotμe1g2expiεe1tge1tξe1t.
(26)

Here μe1g=e1|μ^|g is the excitonic transition dipole moment, and ge1t=ge1e1,e1e1t. The ctR theory term is

ξe1t=e2e1e20tdτ20τ2dτ1eiωe1e2τ2τ1Ce1e2,e2e1τ2τ1.
(27)

The factor of 13 in the response function expression [Eq. (26)] comes from the orientational averaging.45 

2. Linear dichroism

Linear dichroism (LD) spectroscopy is performed on oriented samples (they are no longer isotropic). For simplicity we orient the molecular coordinate system in such a way that the membrane normal lies along the z axis. The linear dichroism spectra then corresponds to a difference between two absorption measurements: one where the incoming light is polarized parallel to the membrane plane, and the other, when the incoming light is polarized perpendicular to the membrane plane46 

LDω=AωAω.
(28)

The required absorption spectra are then defined by

A/ω=ωRe0dteiωtR1t;o/,o/dis.
(29)

The response functions are calculated using the same approximations as above. The relevant orientational average is

μe1go2μe1go2orμe1g,x2+μe1g,y22μe1g,z2.
(30)

Finally, the LD spectrum can be calculated as

LDω=ωRe0dteiωtRLD1tdis,
(31)

with

RLD1t=e1Ntotμe1g,x2+μe1g,y22μe1g,z2×expiεe1tge1tξe1t.
(32)

3. Circular dichroism

The circular dichroism spectrum is defined as a difference between the absorption of left- and right-circularly polarized light21,46,47

CDω=ALωARω.
(33)

To properly describe the CD spectrum, one has to take into account the coupling with the magnetic field. Therefore, in this case, the system-field coupling Hamiltonian is

ĤSF=μ^Etm^Bt=nNEtμnoEe^nh^n+h^ne^nnNBtmnoBe^nh^n+mn*oBh^ne^n.
(34)

Here m^ is the magnetic dipole moment operator and oB is the polarization vector of the magnetic field. Following Ref. 44, we use an effective molecular magnetic dipole moment defined as

mn=i2εn*Rn×μn.
(35)

Here Rn is the position vector corresponding to the center of the pigment n. The CD spectrum can also be calculated from the response function

CDω=ωIm0dteiωtRCD1tdis.
(36)

The relevant response function is found to be44 

RCD1t;o2,o1=TreiĤtμ^o2eiĤtm^o1ρ^ggor.
(37)

For an isotropic medium, we can again set RCD1t=RCD1t;o,o, though note that the same vector o is taken as both arguments. Using the same level of theory as above, the response function is

RCD1t=13e1Ntotμe1gme1g  ×expiεe1tge1tξe1t.
(38)

Here me1g=e1|m^|g is the excitonic magnetic dipole moment.

We note that this level of theory neglects the intrinsic CD of the pigments. In addition, the CD spectra can be affected by couplings to the higher energy states,21,48 which we ignore in the present treatment.

4. Stark spectroscopy

Stark spectroscopy relies on the effect of an applied static electric field on an absorption spectrum and is defined as the difference between the spectrum with the field on and the field off49 

Sω=AEωAω.
(39)

Theoretically, the key issue is to calculate the absorption spectrum with the presence of an electric field. For this, the system-field Hamiltonian is extended to also include the static field

ĤSF=μ^Et+E.
(40)

Here the static electric field is denoted by E, and we assume that the electric dipole operator μ^ stays constant even in the presence of the electric field—therefore we are describing only the linear Stark effect. The term μ^E can be added to the system Hamiltonian because it is constant. Here we can no longer neglect the diagonal elements of the dipole moment operator. In fact, they cause a shift in the excited state energies,

ã|ĤSμ^E|ãg|μ^E|g=εãΔμãE.
(41)

Here,

Δμã=ã|μ^|ãg|μ^|g
(42)

is the difference between the permanent dipole moments of the excited and the ground states and will be called the static dipole moment for convenience. We have to subtract the ground state energy because we have set it to zero. Thus, for simplicity, we redefine the dipole moment operator

μ^μ^Îg|μ^|g
(43)

with Î being the identity operator.

In this work, we follow Ref. 50 and use perturbation theory to calculate the Stark spectrum. To this end, we first write the expression for the linear response function in the presence of the static electric field

RE1t;o2,o1;oE=TreiĤμ^oEEtμ^o2eiĤμ^oEEtμ^o1ρ^ggor=Trexpi0tdτμ^IτoEEeiĤtμ^o2eiĤtexp+i0tdτμ^IτoEEμ^o1ρ^ggor.
(44)

Here oE is the polarization vector of the static electric field, and E is its magnitude. Note that here the interaction picture is assumed with respect to Ĥ [see Eq. (1)], contrary to Eq. (13). If we expand this expression in powers of E (this is done by expanding the time-ordered exponentials and collecting the terms), then the zeroth order term gives the usual absorption spectrum, the first order term vanishes for isotropic samples, and the second order term is responsible for the Stark signal.

Thus, the Stark spectrum response function is

RS3t;oE;oE=E20tdτ20τ2dτ1Trμ^Iτ1oEμ^Iτ2oEeiĤtμ^oEeiĤtμ^oEρ^ggorE20tdτ0tdτTrμ^IτoEeiĤtμ^oEeiĤtμ^IτoEμ^oEρ^ggorE20tdτ20τ2dτ1TreiĤtμ^oEeiĤtμ^Iτ2oEμ^Iτ1oEμ^oEρ^ggor.
(45)

Here we have set o1=o2=oE. Expansion of these three terms in the eigenstate basis of the system Hamiltonian produces many contributions to the Stark spectrum. Nonetheless, the majority of them are nonresonant, i.e., oscillating with factors other than eiω¯et, ω¯e being the band gap between the ground and the single excited state manifolds. In fact, only a single contribution survives the RWA (this limit was denoted as pure RWA in Ref. 18),

RS3t;oE;oE=E20tdτ20τ2dτ1e3Ntote2Ntote1Ntotμge3Eμe3e2Eμe2e1Eμe1gEor×TrBeiĤBte3|eiĤtτ2|e3e2|eiĤτ2τ1|e2e1|eiĤτ1|e1ρ^gg.
(46)

Here we have denoted μabE/E=a|μ^oE/E|b. Then, using the secular approximation and performing the second order cumulant expansion, we obtain the following expression, which includes the ctR theory terms to describe the off-diagonal fluctuations:

RS3t;oE;oE=E20tdτ20τ2dτ1e3Ntote2Ntote1Ntotμge3Eμe3e2Eμe2e1Eμe1gEorexpiωe3gt+ωe2e3τ2+ωe1e2τ1×expge1e1,e1e1τ1ξe1τ1+ge2e2,e1e1τ1+ge2e2,e1e1τ2τ1ge2e2,e1e1τ2ge2e2,e2e2τ2τ1ξe2τ2τ1ge3e3,e1e1τ2τ1+ge3e3,e1e1τ2+ge3e3,e1e1tτ1ge3e3,e1e1t+ge3e3,e2e2τ2τ1+ge3e3,e2e2tτ2ge3e3,e2e2tτ1ge3e3,e3e3tτ2ξe3tτ2.
(47)

If we assume magic angle cosθm.a.=1/3 between the static and probing electric fields, then the orientational averaging results in18,45 

μge3Eμe3e2Eμe2e1Eμe1gEorm.a.=19μe1gμge3μe2e1μe3e2.
(48)

The Stark spectrum is finally calculated as

Sω=ωRe0dteiωtRS3t;oE;oEdis.
(49)

Expression (47) is one of the major results of the present work.

Since the PSII RC has 8 pigments, in addition to 8 ME states, there could be 56 possible CT states. If we limit the CT states to the particular pigments (ChlD1, PheoD1, PD1, and PD2) that have been proposed to be involved in charge separation and take into consideration the directionality of the charge transfer, this number is significantly reduced to 6. Restricting the charge transfer states to neighboring pigment pairs further reduces this number to 3. Given that the majority of the spectra that we fit (absorption, LD, CD) are not significantly affected by the dark CT states, we include a single CT state in our model. We note that only the Stark spectrum is sensitive to CT states, and previous attempts to fit the Stark spectrum of the PSII RC included a single CT state.18 We performed calculations that considered either the PD2+PD1 or ChlD1+PheoD1 CT states and from the fits of the Stark spectrum (see Sec. IV for details) we have chosen ChlD1+PheoD1 as the primary CT state. Thus, our model includes 9 states in the single excitation manifold.

We obtain the site energies of the pigments from the fitting procedure. The resonance couplings between the pigments were taken from Ref. 51, where they were calculated using the Poisson TrEsp method.52–54 The CT state energy and the coupling between CT state and molecular states were also obtained from the fits. Formally, there are two couplings involving a single CT state,

PheD1*|ĤS|ChlD1+PheoD1=tPheoD1ChlD1h,
(50)
ChlD1*|ĤS|ChlD1+PheoD1=tPheoD1ChlD1e,
(51)

but we have chosen tmne=tmnh for simplicity. The entire single excitation Hamiltonian matrix is given in Table I.

TABLE I.

Single excitation Hamiltonian matrix elements (in cm1), energies of the states in nanometers, and standard deviations of energy distributions (in cm1). See the text for details.

PD1PD2ChlD1ChlD2PheoD1PheoD2ChlzD1ChlzD2ChlD1+PheoD1Eλσ
PD1 15 080 158 −27.3 −41.8 −4 12.6 0.5 0.6 663 80 
PD2 158 15 015 −46.8 −22 15.1 −3 0.6 0.6 666 80 
ChlD1 −27.3 −46.8 14 800 3.5 43.5 −2.2 1.7 −0.1 125 675.5 80 
ChlD2 −41.8 −22 3.5 15 010 −2.4 41.7 −0.1 1.8 666 80 
PheoD1 −4 15.1 43.5 −2.4 14 950 1.5 −2.5 −0.2 125 669 80 
PheoD2 12.6 −3 −2.2 41.7 1.5 14 850 −0.2 −2.6 673.5 80 
ChlzD1 0.5 0.6 1.7 −0.1 −2.5 −0.2 14 985 0.1 667.5 80 
ChlzD2 0.6 0.6 −0.1 1.8 −0.2 −2.6 0.1 15 065 664 80 
ChlD1+PheoD1 125 125 14 270 701 550 
PD1PD2ChlD1ChlD2PheoD1PheoD2ChlzD1ChlzD2ChlD1+PheoD1Eλσ
PD1 15 080 158 −27.3 −41.8 −4 12.6 0.5 0.6 663 80 
PD2 158 15 015 −46.8 −22 15.1 −3 0.6 0.6 666 80 
ChlD1 −27.3 −46.8 14 800 3.5 43.5 −2.2 1.7 −0.1 125 675.5 80 
ChlD2 −41.8 −22 3.5 15 010 −2.4 41.7 −0.1 1.8 666 80 
PheoD1 −4 15.1 43.5 −2.4 14 950 1.5 −2.5 −0.2 125 669 80 
PheoD2 12.6 −3 −2.2 41.7 1.5 14 850 −0.2 −2.6 673.5 80 
ChlzD1 0.5 0.6 1.7 −0.1 −2.5 −0.2 14 985 0.1 667.5 80 
ChlzD2 0.6 0.6 −0.1 1.8 −0.2 −2.6 0.1 15 065 664 80 
ChlD1+PheoD1 125 125 14 270 701 550 

Our model takes into account the disorder value of the energies of the pigments and CT states. We assume the disorder value for all pigments to be the same, and the CT state to have different disorder values. In both cases, the disorder value was described by a Gaussian distribution with standard deviations obtained from the fits (given in Table I).

The strengths of the transition dipole moments of Chl and Pheo pigments are taken to be 4.4 and 3.4 D, respectively.27,29 The transition dipole vectors are assumed to be directed from NB to ND atoms (for nomenclature see, e.g., Ref. 61) as obtained from the 3WU2 crystal structure.3 The CT state was assumed to have no transition dipole strength. The static dipole moments are assumed to be rotated by 15° in the direction of the NC atom from the NB to ND axis (this was based on the calculations of absorption difference spectra corresponding to the PD1+PheoD1 CT state, see Sec. IV for details). Their strengths are taken to be 1.55 and 1.3 D for Chls and Pheos, respectively. The static dipole moment of the CT state was calculated by assuming charges fully localized at the geometric centers of the four N atoms of the involved pigments. This gives a static dipole strength of 49.3 D for the ChlD1+PheoD1 state. The transition dipole moments, the static dipole moments, and the geometric centers of the pigments of the PSII RC are given in Table II. Here the z axis is perpendicular to the membrane plane, with data taken from the 3WU2 crystal structure3 obtained from the OPM database.55 

TABLE II.

The transition dipole moments in D, the static dipole moments in D, and geometric center of the pigments in Å. Data taken from the 3WU2 crystal structure3 obtained from OPM database.55 

μxμyμzΔμxΔμyΔμzRxRyRz
PD1 1.08 −3.74 2.05 0.32 −1.07 1.07 37.87 4.05 11.99 
PD2 −0.66 3.73 2.24 −0.13 1.03 1.15 32.57 9.55 11.12 
ChlD1 −0.33 4.06 1.67 0.29 1.34 0.73 30.65 −2.13 9.30 
ChlD2 0.11 −4.04 1.74 −0.36 −1.31 0.75 40.19 15.45 8.24 
PheoD1 1.65 −0.55 −2.92 0.84 0.10 −0.99 33.73 −4.13 −0.31 
PheoD2 −1.77 0.19 −2.90 −0.89 −0.19 −0.93 36.67 16.71 −1.41 
ChlzD1 −1.91 −0.64 3.91 −1.01 −0.32 1.13 31.92 −27.24 8.07 
ChlzD2 2.07 0.98 3.76 1.06 0.41 1.05 38.91 39.35 4.76 
μxμyμzΔμxΔμyΔμzRxRyRz
PD1 1.08 −3.74 2.05 0.32 −1.07 1.07 37.87 4.05 11.99 
PD2 −0.66 3.73 2.24 −0.13 1.03 1.15 32.57 9.55 11.12 
ChlD1 −0.33 4.06 1.67 0.29 1.34 0.73 30.65 −2.13 9.30 
ChlD2 0.11 −4.04 1.74 −0.36 −1.31 0.75 40.19 15.45 8.24 
PheoD1 1.65 −0.55 −2.92 0.84 0.10 −0.99 33.73 −4.13 −0.31 
PheoD2 −1.77 0.19 −2.90 −0.89 −0.19 −0.93 36.67 16.71 −1.41 
ChlzD1 −1.91 −0.64 3.91 −1.01 −0.32 1.13 31.92 −27.24 8.07 
ChlzD2 2.07 0.98 3.76 1.06 0.41 1.05 38.91 39.35 4.76 

We assume that fluctuations of the energies are uncorrelated, thus CABω=CAδAB. All pigments (both Chls and Pheos) are described by the same spectral density consisting of two parts, Cm*ω=CLFω+CHFω. In our previous work,29 we assumed that the low frequency part was given by a Debye spectral density. However, this form is unsuitable for low temperature spectra calculations due to the divergence of the total Huang-Rhys factor.62 Instead we adopt the superohmic spectral density that is free from this deficiency, consisting of two terms,

CLFω=j=1,2λjπ2ωωjc3expω/ωjc.
(52)

The two cutoff frequencies were chosen to be ω1c=20cm1 and ω2c=80cm1. The values for reorganization energy were obtained from the fits of spectra and found to be λ1=15cm1 and λ2=35cm1. The high frequency part is represented by 8 underdamped vibrational modes,

CHFω=j=184sjWjω3Γω2Wj2Γ22+4ω2Γ2.
(53)

The damping parameter was chosen to be the same for all modes, Γ1=1ps. The frequencies Wj and the Huang-Rhys factors sj were taken from ΔFLN measurements of Chl a in solution56 and are given in Table III. We included only the modes with center frequencies less than 750 cm1, as higher frequency modes are not expected to influence system dynamics, being completely nonresonant. For the CT state, we keep the high frequency part of the spectral density the same but introduce an additional multiplicative factor for the low frequency part, as CT states are expected to couple more strongly to the environment. Thus, CCTω=ηCLFω+CHFω. From the fits of the Stark spectrum, we obtained η=12 for the ChlD1+PheoD1 state.

TABLE III.

Parameters of the underdamped vibrations included in the spectral density (from ΔFLN measurements of Chl a in solution56).

Mode12345678
Wj,cm1 93 189 262 352 389 519 573 745 
sj 0.012 0.010 0.013 0.023 0.015 0.016 0.017 0.034 
Mode12345678
Wj,cm1 93 189 262 352 389 519 573 745 
sj 0.012 0.010 0.013 0.023 0.015 0.016 0.017 0.034 

Most of the model parameters were obtained from a fit of nine linear spectra shown in Fig. 2. In particular, we fitted absorption at 6, 77, and 298 K, linear and circular dichroism at 6 and 298 K, and absorption difference spectra with exchanged pheophytins at 6 K. The latter spectra correspond to the so-called RC1x and RC2x complexes,57,58 where in the RC1x complexes, 80% of the PheoD2 is estimated to be replaced by 131-OH-Pheo, and in the RC2x complexes, all PheoD2 and 50% of PheoD1 are replaced. When calculating RC1x/RC2x spectra, we took this into account accordingly. The energies of the exchanged Pheos were also free parameters of the fit. In the fitting procedure, we minimized the sum of the deviations from the experimental data at each wavelength for each spectra and added them with different weights (1.0 for Abs 6 K and 298 K, 2.0 for Abs 77 K, 5.0 for Abs RC-RC1x/RC2x 6 K, 0.05 for LD 6 K, 0.005 for CD 6 K and 298 K, and 0.01 for LD 298 K). Smaller weights were chosen for some spectra due to the poorer quality of the fit. Before calculating the differences, the calculated spectra were normalized in the following way: absorption spectra were normalized by setting the maximum amplitudes to unity, CD/LD spectra were normalized in such a way that the difference between the experimental and calculated spectra would be the smallest, and for absorption difference RC and RC1x/RC2x spectra were normalized to unity separately (this was done also to the experimental data) and then subtracted. Fitting was performed using R software63 with a differential evolution algorithm64 as implemented in Ref. 65.

FIG. 2.

Fits of the PSII RC spectroscopic data. Open circles denote experimental values, while red lines denote theoretical calculations. (a) Absorption at 6 K. (b) Linear dichroism at 6 K. (c) Circular dichroism at 6 K. (d) Absorption at 77 K. (e) Absorption difference RC − RC1x preparations (with pheophytin exchange) at 6 K. (f) Absorption difference RC − RC2x preparations (with pheophytin exchange) at 6 K. (g) Absorption at 298 K. (h) Linear dichroism at 298 K. (i) Circular dichroism at 298 K. Experimental data for (a) are taken from Ref. 57, for (b), (c), (e), (f), and (i) taken from Ref. 58, for (d) from Ref. 59, and for (g) from Ref. 60.

FIG. 2.

Fits of the PSII RC spectroscopic data. Open circles denote experimental values, while red lines denote theoretical calculations. (a) Absorption at 6 K. (b) Linear dichroism at 6 K. (c) Circular dichroism at 6 K. (d) Absorption at 77 K. (e) Absorption difference RC − RC1x preparations (with pheophytin exchange) at 6 K. (f) Absorption difference RC − RC2x preparations (with pheophytin exchange) at 6 K. (g) Absorption at 298 K. (h) Linear dichroism at 298 K. (i) Circular dichroism at 298 K. Experimental data for (a) are taken from Ref. 57, for (b), (c), (e), (f), and (i) taken from Ref. 58, for (d) from Ref. 59, and for (g) from Ref. 60.

Close modal

We obtain excellent fits for all the absorption experiments, and good qualitative fits for circular dichroism (the CD spectrum of PSII RC is nonconservative, while our theoretical approach can give only conservative spectra in the Qy region). There is some quantitative difference for the calculated linear dichroism spectra, with positive feature 680nm having less amplitude than in the experimental data. Overall, our results are of similar quality to previous studies on the PSII RC.16,18

We have adjusted the CT state parameters (its energy, coupling with molecular states, disorder values and multiplicative factor describing coupling with low frequency vibrations) from the fits of the Stark spectrum measured at 77 K,36 with results given in Fig. 3. During the fitting procedure, both the experimental and the simulated Stark spectra were normalized with respect to the absorption spectrum. We tried including either PD2+PD1 or ChlD1+PheoD1 CT states in our model. We found that a better fit was obtained with the ChlD1+PheoD1 state. Results with the PD2+PD1 state are also reasonable but are somewhat worse—the amplitude of the 680nm is underestimated, and a negative feature 700nm appears, which is absent in the experimental data. Therefore, we hold that ChlD1+PheoD1 should be the primary CT state, contrary to the conclusions obtained previously.18 We also note that calculations without any CT states roughly reproduce the amplitude of the 670nm feature but severely underestimate the amplitude of the 680nm feature. This is a clear indication that CT states significantly influence the Stark spectrum of the PSII RC, much more than in previous calculations.18 

FIG. 3.

Fit of Stark spectrum of the PSII RC at 77 K. Open circles denote experimental values (data taken from Ref. 36), red line corresponds to calculations with (a) ChlD1+PheoD1 CT state, (b) PD2+PD1 CT state, (c) both CT states. The parameters for CT states were fitted separately in (a) and (b). For (c), calculations were made by taking parameters corresponding to the best fit from (a) and (b) and were not further adjusted. Blue dashed line denotes calculations without any CT states.

FIG. 3.

Fit of Stark spectrum of the PSII RC at 77 K. Open circles denote experimental values (data taken from Ref. 36), red line corresponds to calculations with (a) ChlD1+PheoD1 CT state, (b) PD2+PD1 CT state, (c) both CT states. The parameters for CT states were fitted separately in (a) and (b). For (c), calculations were made by taking parameters corresponding to the best fit from (a) and (b) and were not further adjusted. Blue dashed line denotes calculations without any CT states.

Close modal

Having performed the fits and obtained the optimal values of the parameters, we checked our model by calculating other linear spectra, for which experimental results are available in the literature.

We calculated the CD spectra of different PSII RC preparations (normal RC, RC1x, and RC2x) at 6 K. The comparison of calculations with experimental data (taken from Ref. 58) is presented in Fig. 4. Our simulations qualitatively reproduce the changes observed in the CD spectra upon the increase of exchanged Pheo (going from RC to RC1x to RC2x). First, the main positive feature at 680nm loses a significant amount of amplitude. Second, a positive feature at 675nm appears. Third, the negative feature at 670nm loses its amplitude. Fourth, a negative feature at 655nm appears.

FIG. 4.

Comparison of 6 K CD spectra of different PSII RC preparations (with partial exchange of Pheos). Experimental data are (taken from Ref. 58) denoted by symbols on the left, and simulations are denoted by lines on the right.

FIG. 4.

Comparison of 6 K CD spectra of different PSII RC preparations (with partial exchange of Pheos). Experimental data are (taken from Ref. 58) denoted by symbols on the left, and simulations are denoted by lines on the right.

Close modal

We also calculated LD spectra at 6 K of RC, RC1x, and RC2x preparations. Comparison with experimental data is presented in Fig. 5. Again, our simulations qualitatively reproduce the main changes upon the increase of exchanged Pheo—appearance of negative feature at 655nm, lack of changes in negative feature at 670nm, and increase of amplitude of positive feature at 680nm. Nonetheless, we have to note that the overall amplitudes of both the negative feature at 655nm and the positive feature at 680nm are considerably smaller in our simulations.

FIG. 5.

Comparison of 6 K LD spectra of different PSII RC preparations (with partial exchange of Pheos). Experimental data (taken from Ref. 58) are denoted by symbols on the left, simulations are denoted by lines on the right.

FIG. 5.

Comparison of 6 K LD spectra of different PSII RC preparations (with partial exchange of Pheos). Experimental data (taken from Ref. 58) are denoted by symbols on the left, simulations are denoted by lines on the right.

Close modal

In Fig. 6, we compare the calculated triplet minus singlet absorption difference spectrum of the PSII RC with the experimental data.58 Calculations were performed assuming triplet localized at ChlD1 and ignoring its couplings and dipole strength. We obtain good agreement with experimental data, similar to other studies.16,18,30

FIG. 6.

Triplet minus singlet absorption difference spectrum at 10 K. Open circles denote experimental values (data taken from Ref. 58), red line corresponds to calculations, which were performed assuming triplet localized at ChlD1.

FIG. 6.

Triplet minus singlet absorption difference spectrum at 10 K. Open circles denote experimental values (data taken from Ref. 58), red line corresponds to calculations, which were performed assuming triplet localized at ChlD1.

Close modal

We calculated absorption difference spectra of normal RC preparations (here called RC6) and the so-called RC5 preparations, which lack one of the peripheral chlorophylls.66,67 Comparison with the experimental data is presented in Fig. 7. Experimental spectra were calculated by normalizing both RC6 and RC5 absorption at the 678nm peak and subtracting the spectra, then normalizing with respect to the RC6 spectra. Simulated spectra were obtained by calculating normal RC6 spectra, then calculating the spectra without the ChlzD1/D2 pigment, subtracting, and normalizing with respect to the RC6 spectra. Our calculations are consistent with ChlzD2 being lost in the RC5 preparations, with both peak position at 666nm and its half-width agreeing with the experimental data. The amplitude of the 666nm peak is higher in the calculations, which could be related to differences in the low temperature experimental absorption spectra, as the amplitude of 672nm peak varies in different reports.57–59,66

FIG. 7.

Comparison of experimental and calculated RC6-RC5 absorption difference spectra at 4 K (left, experimental data taken from Ref. 66) and 77 K (right, experimental data taken from Ref. 67). Experimental data are denoted by open circles, red lines denote calculations done by neglecting ChlzD1, and blue lines—neglecting ChlzD2.

FIG. 7.

Comparison of experimental and calculated RC6-RC5 absorption difference spectra at 4 K (left, experimental data taken from Ref. 66) and 77 K (right, experimental data taken from Ref. 67). Experimental data are denoted by open circles, red lines denote calculations done by neglecting ChlzD1, and blue lines—neglecting ChlzD2.

Close modal

In Fig. 8, we show a comparison of experimental and calculated RC5-RC5PheoD1 absorption difference spectra (here RC5PheoD1 means RC5 with reduced PheoD1), with experimental data taken from Ref. 68. Calculations were done assuming that in RC5 preparations the ChlzD2 is lost, consistent with our RC6-RC5 absorption difference calculations. When calculating spectra with reduced PheoD1, we neglected the couplings and dipole strength of this pigment. In addition, as suggested in Ref. 16, we took the electrochromic shifts into account. This was done by assuming a negative charge localized at the geometric center of the four nitrogen atoms of PheoD1 and then calculating the interaction of that charge with the static dipole moments of the other pigments, that is,

Δεn*=14π𝜖0𝜖reRPheoD1RnΔμnRPheoD1Rn3.
(54)

Here 𝜖r is the relative dielectric constant, which was assumed to be equal to 2.5 in our calculations. The CT state was neglected in the calculations for the reduced PheoD1 since we expect that due to its huge energy shift, it would effectively decouple from the ME states. As seen from Fig. 8, we obtain excellent agreement for 77 K data, with both the main negative feature at 682nm and the positive feature at 673nm reproduced. In calculations for 277 K, we still reproduce the main negative feature but obtain a larger amplitude for the positive one compared to the experimental data.

FIG. 8.

Comparison of experimental and calculated RC5-RC5 with PheoD1 reduced absorption difference spectra at 77 K (left) and 277 K (right). Experimental data (taken from Ref. 68) are denoted by open circles, and red lines denote calculations.

FIG. 8.

Comparison of experimental and calculated RC5-RC5 with PheoD1 reduced absorption difference spectra at 77 K (left) and 277 K (right). Experimental data (taken from Ref. 68) are denoted by open circles, and red lines denote calculations.

Close modal

We also calculated the absorption difference spectrum corresponding to the CT state PD1+PheoD1 taken from pump–probe measurements.19,60,69 This is presented in Fig. 9 for 10, 77, and 298 K. In calculations of the absorption with the CT state, we neglected the couplings and dipole strength of the involved pigments and accounted for the electrochromic shifts. The effective dielectric constant 𝜖r was assumed to be equal to 2.5, as before. We obtain good agreement between experiment and calculations for all temperatures. It is important to note that agreement was obtained only by assuming an angle of 15° between NB and ND axes towards the NC atom for the static dipole moments of both Chls and Pheos.

FIG. 9.

Comparison of experimental and calculated absorption difference spectra corresponding to the CT state PD1+PheoD1. Experimental data denoted by open circles, red lines denote calculations. (a) 10 K (curve 2 from Fig. 2 of Ref. 69). (b) 77 K (EADS corresponding to 20 ns component of 660 nm excitation from Ref. 19). (c) 298 K (DADS corresponding to non-decaying component of 400 nm excitation from Ref. 60).

FIG. 9.

Comparison of experimental and calculated absorption difference spectra corresponding to the CT state PD1+PheoD1. Experimental data denoted by open circles, red lines denote calculations. (a) 10 K (curve 2 from Fig. 2 of Ref. 69). (b) 77 K (EADS corresponding to 20 ns component of 660 nm excitation from Ref. 19). (c) 298 K (DADS corresponding to non-decaying component of 400 nm excitation from Ref. 60).

Close modal

We obtained our site energies from the fit of nine linear spectra and then verified them by calculating more independent spectra. In this subsection, we will compare our results to other work. We have to note, however, that direct comparison of the site energies can sometimes be misleading, due to the type of spectral density employed and the amount of reorganization energy included in different models. Therefore, comparing the site energies relatively, e.g., between the corresponding pigments of different branches of the RC, is a better approach.

Our model finds that the lowest site energy is of the ChlD1 pigment, in agreement with many previous studies.16,18,29,30 This is consistent with the triplet minus singlet absorption difference spectrum calculations, which can be well described assuming that at low temperatures, the triplet is localized at the ChlD1 pigment. On the other hand, some studies assign a considerably higher energy for this pigment.70,71 Our fitting shows that the latter assignment is not consistent with the experimental data.

Another somewhat controversial issue is the assignment of the PheoD1/D2 energies. We find that PheoD2 has the second lowest energy, while that of PheoD1 is considerably higher. This is in general agreement with the findings of Ref. 30, although contrary to other reports in the literature.72 Since we can accurately reproduce not only absorption difference spectra RC-RC1x/RC2x but also absorption difference corresponding to reduced PheoD1 and the CT state PD1+PheoD1, we hold that our assignment accurately reflects the system.

In previous work, the energies of the special pair pigments have mostly been found to be the highest of the pigments in the RC.16,18,20,26,29 Our findings suggest high, but not necessarily the highest energies for these pigments, with ChlD2 having a similar energy to PD1 and ChlzD2 to PD2. Due to a relatively large resonance coupling between the special pair pigments, it is rather difficult to pinpoint their energies from a specific experiment. Another issue to consider is that due to very close proximity between the special pair pigments,3 their orbitals are likely overlapping. Thus, even the molecular excited states of PD1 and PD2 pigments might actually have partially delocalized or even CT character. Their excited state energies are likely to be only effective. Nonetheless, having simulated the absorption difference spectra corresponding to the PD1+PheoD1 CT state and obtained good agreement with experimental data, our energy for the PD1 pigment should be accurate within the limits of our model.

Regarding the peripheral chlorophylls, our results are in general agreement with other studies.16,18,29,30 Nonetheless, previous models often suggested very similar energies between the two Chlz molecules, while in the present work, we found a somewhat larger energy for ChlzD2, which suggests that it is the pigment lost in the RC5 preparations.66,67 This is in agreement with the assignment from Ref. 30.

A convenient way to compare the site energies is to present the excited state energy distributions that show the energies of the eigenstates in which each pigment or CT state participates. These distributions are defined as

Dãω=eNtotψe,ã2δωωegλee,eedis.
(55)

Here λee,ee=ãψe,ã4λã is the eigenstate reorganization energy. For our model, these distributions are presented in Fig. 10. In addition, we also present distributions from our previous model29 and from Ref. 30, where the energies of the sites were obtained by taking into account electrostatic interactions as calculated from the structural data.

FIG. 10.

The excited state energy distributions showing the energies of the states that each pigment or CT state participate in. See the text for more details. Black lines denote the present model, blue lines denote our previous work,29 and red lines denote recent work from Müh et al.30 

FIG. 10.

The excited state energy distributions showing the energies of the states that each pigment or CT state participate in. See the text for more details. Black lines denote the present model, blue lines denote our previous work,29 and red lines denote recent work from Müh et al.30 

Close modal

It can be seen that there is good agreement between our current model and that of Ref. 30. We note that this agreement is because Eq. (55) accounts for differences in reorganization energy, thus even though our model has higher state energies, it also has larger reorganization energies and the combined result is very similar to Ref. 30. The only noticeable difference is in the special pigment energy distributions that are somewhat blueshifted in our model.

On the other hand, the differences between the excited state distributions of the present model and our previous work29 are larger, though the overall pattern is somewhat similar. The largest differences are for ChlD2, PheoD1, and both Chlzs. One main difference is that the energy distributions for our previous work show more structure. This is because of the differences in disorder values between both models. Our previous work,29 following previous models based on the modified Redfield theory,18,20,26 assumed about three times smaller disorder than in the current model. Taking into account that presently we use a more advanced lineshape theory and a spectral density suitable for low temperature calculations and that the disorder value that we obtain is very similar to Refs. 16, 27, and 30, we believe that our current values are more accurate and that disorder in the previous studies18,20,26,28,29 might have been severely underestimated. These differences in disorder are expected to have significant influence on the time-resolved spectroscopic data and will have to be investigated more closely in the future.

In this work, we extended the ctR theory to describe the Stark spectrum. This was done using the perturbative expansion in terms of the strength of the static electric field. Our expressions thus significantly improve upon Ref. 18, where the Markovian approximation and the modified Redfield lineshapes were used. Even though here we presented only the term corresponding to the RWA, other terms can be easily formulated. Nonetheless, our approach, though straightforward to compute, is computationally expensive, involving a double time integral over a triple sum. This means that it scales rather badly for larger systems. Additionally, the Stark spectrum is sensitive to the CT states and our results indicate that they possess very large disorder. Thus, Stark spectrum calculations need a large number of realizations to converge. We performed our fits with 25 000 realizations, and the obtained spectra were still somewhat noisy. The Stark spectra presented in this paper (Fig. 3) were calculated using 600 000 realizations. For comparison, 5000 realizations are sufficient to achieve good convergence for the absorption spectrum at 77 K. Although this suggests that our determined CT state parameters might possess some uncertainty, they still should be reasonably accurate.

Several conclusions can be made from our calculations of the Stark spectrum of the PSII RC. First, our results show that CT states influence the Stark spectrum more than was obtained previously.18 Without the CT states, we cannot obtain a reasonable agreement with experimental data. This implies that the coupling between the ME states and the CT states is somewhat larger than was estimated before.18 If so, this should affect calculations of the time-resolved spectroscopic data.

Second, our results indicate that if the Stark spectrum is mostly influenced by a single CT state, it is probably ChlD1+PheoD1 rather than PD2+PD1, contrary to what was previously concluded.18 Although the agreement between theory and experiment is poorer with the PD2+PD1 state than with ChlD1+PheoD1, the inclusion of the PD2+PD1 state does provide better agreement than when CT states are neglected entirely. The agreement that we obtain with the single ChlD1+PheoD1 state [Fig. 3(a)] is good but not excellent. It is likely that instead of a single state, more CT states are coupled to the ME states strongly enough to influence the Stark spectrum. We did not try to fit the experimental Stark spectrum with 2 CT states because the fit would have been very poorly constrained (fitting a single spectrum with eight free parameters). Instead we chose to make separate optimizations with the individual ChlD1+PheoD1 and PD2+PD1 states. Using the parameters obtained from the individual optimizations, if we include both CT states [see Fig. 3(c)] we obtain improved agreement of the 680nm feature [if compared to Fig. 3(a)] but at a cost of an appearance of a spurious negative feature at 700nm.

Third, we have obtained very large disorder of the CT state (σChlD1+PheoD1=550cm1), which is clearly visible in its energy distribution in Fig. 10. The obtained value of σChlD1+PheoD1 is about seven times larger than that of the ME states. This should be related to the fact that CT states have large static dipole moments, thus they are more susceptible to the changes in the environment. Also we note that our value of disorder is considerably larger than previous estimates.18,20,26,29 This could be due to the shortcomings of the modified Redfield theory, which was used in those studies. This result is in line with previous suggestions that disorder might control different charge separation pathways in the PSII RC.

Overall, we would like to suggest that calculations of the Stark spectrum of the PSII RC should be made not to extract the values of the system parameters (due to the computational expense) but to constrain them. Thus, models which include a number of CT states needed to calculate the charge separation dynamics20,29 should be tested against the experimental Stark spectrum using the level of theory developed in this paper. Since each CT state necessarily adds a number of parameters to the model, their values could be severely constrained by their sensitivity to the Stark spectrum.

We would also like to draw attention to our calculations of the absorption difference spectra corresponding to the CT state PD1+PheoD1 at different temperatures (Fig. 9). The good agreement between theory and experiment implies that even at room temperature the relaxed pump-probe spectra correspond to this specific state. Therefore, the energy gap between this state and all the other CT states should be relatively large, a few times more than the thermal energy at room temperature. This information should be useful in simulating time-resolved spectroscopic data.

In previous theoretical studies on the PSII RC, fitting of relaxed fluorescence was employed to obtain the model parameters.16,18,26 In our view, this approach is problematic and could lead to erroneous results. Since a normal RC performs charge separation with high efficiency, it should barely fluoresce. Thus, the majority of the fluorescence may originate from those reaction centers that cannot perform charge separation, due to their particular realization of static disorder that provides an unfavorable energy landscape. Fits of the relaxed fluorescence would then depend not on the normal RC, but only on a specific subensemble. A better approach to calculations of the fluorescence of the RC may be to simulate the time-resolved fluorescence and then integrate over all times. Unfortunately this requires including all the relevant CT states in the model and thus is out of scope of the present study.

In the present work, we used the ctR theory to describe the optical lineshapes of the PSII RC. We previously tested the ctR theory on a molecular dimer35 and found that it provided excellent agreement with an exact approach. In addition, it was found to be considerably more accurate than the modified Redfield approach, which has often been used for calculations of the PSII RC optical spectra.18,20,26,28,29 The ctR theory was also somewhat more accurate than the Markovian Redfield approach (denoted by us as cR theory in Ref. 35), which was also applied to the PSII RC.16,27,30,73 As the cR and the ctR theories give reasonably close results, it is not surprising that our present work is in general agreement with the latest work where the cR theory was used.30 

We note that the ctR theory is formulated using the secular approximation. In usual circumstances, non-secular terms are expected to make small contributions to spectral lineshapes.35,44 On the other hand, when reorganization energies are very large74 or differing between the states, like between ME and CT states,75 the effects of non-secular terms might increase considerably. On the other hand, they should be strongly reduced, when mixing between the ME and CT states is relatively small. This is the case for our model, due to the large energy gap between the ME and CT states and huge disorder of the latter. Thus we expect that non-secular terms would play only a minor role in our model. Another issue is that the inclusion of non-secular terms is non-trivial and would increase the computational costs substantially, especially for the Stark spectrum.

Here we assumed the same form of the spectral density for CT as for ME states, with the only difference being a multiplicative factor for the low frequency term. Even though the actual shape of the spectral density of the CT state might be considerably different, there is no way to estimate it. Thus, our assumption is the simplest one currently possible. Note that we used the same spectral density for simulations of data ranging from 4 K to room temperature. This is based on the harmonic bath approximation. We believe it is valid, as our simulations match experimental data at all temperatures. Additionally, we note that the degree of anharmonicity was found to be small for a photosynthetic Fenna-Matthews-Olson complex.76 

We also assumed independent fluctuations for all states. Naturally, fluctuations of ME and CT states sharing the same molecules may be correlated. Nonetheless, currently there seems to be no reliable way of estimating these correlations for the PSII RC. On the other hand, parametrization in terms of correlation coefficients would increase the number of free parameters in our model. Thus, we chose not to do this, for simplicity.

Recently a focus has shifted to the question and relevance of quantum coherence, mostly related to the vibrational degrees of freedom, in the PSII RC. This interest was stimulated by simultaneous publications demonstrating coherent oscillations in two-dimensional electronic spectra.77,78 The possibility of coherence enhanced charge separation was either inferred78 or demonstrated in a proof of concept manner.77 To describe the influence of vibrations, one must either include them explicitly in the Hamiltonian79 or in the spectral density. These descriptions should be equivalent if accurate theories are used. The former approach, however, better captures the influence of vibrations if perturbative approaches are employed. However, the computational complexity then increases significantly. In the present work, we include several vibrational modes in the spectral density. The ctR approach allows us to qualitatively capture their influence,35 which is not the case for the cR or modified Redfield approaches. Of course, an exact approach, like the HEOM theory80 or stochastic calculations of path integrals,74 should be used if possible. Nonetheless, we believe that the ctR theory provides a very good balance between accuracy and computational cost. The latter is very important to enable extensive fitting of experimental spectra, as was done in this work.

We have revisited the spectroscopic properties of the PSII RC using the ctR theory. We constructed a model based on the tight-binding Hamiltonian, including a single CT state. The parameters of the model were obtained by fitting nine linear spectra and then further tested by calculations of other independent spectra. Our simulations are in good to excellent agreement with all the calculated spectra. The parameter values of our model differ from previous models based on the modified Redfield theory and are mostly in agreement with recent structure-based estimations that accounted for electrostatic effects. Additionally, we extended the ctR theory to calculate the Stark spectrum and applied it to the PSII RC. Our results showed that no agreement with experimental data could be obtained without an inclusion of CT states. We concluded that ChlD1+PheoD1 is the primary CT state, but the Stark spectrum is probably also influenced by other states. Finally, we suggest that calculations of the Stark spectrum should be used to constrain the CT state parameters obtained by fitting the time-resolved data. Our future work will be directed towards simulating the 2D spectrum of the PSII RC using the ctR theory.

This work was supported by the Research Council of Lithuania (LMT Grant No. MIP-080/2015). J.P.O. acknowledges support from the Office of Basic Energy Sciences, the US Department of Energy under Grant No. DE-SC0016384. Computations were performed using the resources of the High Performance Computing Center “HPC Sauletekis” at Faculty of Physics, Vilnius University.

1.
H.
van Amerongen
and
R.
Croce
,
Photosynth. Res.
116
,
251
(
2013
).
2.
R. E.
Blankenship
,
Molecular Mechanisms of Photosynthesis
, 2nd ed. (
Wiley Blackwell
,
Chichester
,
2014
).
3.
Y.
Umena
,
K.
Kawakami
,
J.-R.
Shen
, and
N.
Kamiya
,
Nature
473
,
55
(
2011
).
4.
O.
Nanba
and
K.
Satoh
,
Proc. Natl. Acad. Sci. U. S. A.
84
,
109
(
1987
).
5.
P. J.
van Leeuwen
,
M. C.
Nieveen
,
E. J.
van de Meent
,
J. P.
Dekker
, and
H. J.
van Gorkom
,
Photosynth. Res.
28
,
149
(
1991
).
6.
S. R.
Greenfield
and
M. R.
Wasielewski
,
Photosynth. Res.
48
,
83
(
1996
).
7.
J.
Dekker
and
R.
Van Grondelle
,
Photosynth. Res.
63
,
195
(
2000
).
8.
B. A.
Diner
and
F.
Rappaport
,
Annu. Rev. Plant Biol.
53
,
551
(
2002
).
9.
F.
Rappaport
and
B. A.
Diner
,
Coord. Chem. Rev.
252
,
259
(
2008
).
10.
G.
Renger
and
T.
Renger
,
Photosynth. Res.
98
,
53
(
2008
).
11.
S.
Vassiliev
and
D.
Bruce
,
Photosynth. Res.
97
,
75
(
2008
).
12.
T.
Renger
and
E.
Schlodder
,
J. Photochem. Photobiol., B
104
,
126
(
2011
).
13.
M. E.
van Brederode
and
R.
van Grondelle
,
FEBS Lett.
455
,
1
(
1999
).
14.
V. I.
Prokhorenko
and
A. R.
Holzwarth
,
J. Phys. Chem. B
104
,
11563
(
2000
).
15.
B. A.
Diner
,
E.
Schlodder
,
P. J.
Nixon
,
W. J.
Coleman
,
F.
Rappaport
,
J.
Lavergne
,
W. F. J.
Vermaas
, and
D. A.
Chisholm
,
Biochemistry
40
,
9265
(
2001
).
16.
G.
Raszweski
,
W.
Saenger
, and
T.
Renger
,
Biophys. J.
88
,
986
(
2005
).
17.
A. R.
Holzwarth
,
M. G.
Müller
,
M.
Reus
,
M.
Nowaczyk
,
J.
Sander
, and
M.
Rögner
,
Proc. Natl. Acad. Sci. U. S. A.
103
,
6895
(
2006
).
18.
V. I.
Novoderezhkin
,
J. P.
Dekker
,
H.
van Amerongen
, and
R.
van Grondelle
,
Biophys. J.
93
,
1293
(
2007
).
19.
E.
Romero
,
I. H. M.
van Stokkum
,
V. I.
Novoderezhkin
,
J. P.
Dekker
, and
R.
van Grondelle
,
Biochemistry
49
,
4300
(
2010
).
20.
V. I.
Novoderezhkin
,
E.
Romero
,
J. P.
Dekker
, and
R.
van Grondelle
,
ChemPhysChem
12
,
681
(
2011
).
21.
H.
van Amerongen
,
L.
Valkunas
, and
R.
van Grondelle
,
Photosynthetic Excitons
(
World Scientific
,
Singapore
,
2000
).
22.
J. R.
Durrant
,
D. R.
Klug
,
S. L.
Kwa
,
R.
van Grondelle
,
G.
Porter
, and
J. P.
Dekker
,
Proc. Natl. Acad. Sci. U. S. A.
92
,
4798
(
1995
).
23.
N.
Kamiya
and
J.-R.
Shen
,
Proc. Natl. Acad. Sci. U. S. A.
100
,
98
(
2003
).
24.
J.
Biesiadka
,
B.
Loll
,
J.
Kern
,
K.-D.
Irrgang
, and
A.
Zouni
,
Phys. Chem. Chem. Phys.
6
,
4733
(
2004
).
25.
B.
Loll
,
J.
Kern
,
W.
Saenger
,
A.
Zouni
, and
J.
Biesiadka
,
Nature
438
,
1040
(
2005
).
26.
V. I.
Novoderezhkin
,
E. G.
Andrizhiyevskaya
,
J. P.
Dekker
, and
R.
van Grondelle
,
Biophys. J.
89
,
1464
(
2005
).
27.
G.
Raszewski
,
B. A.
Diner
,
E.
Schlodder
, and
T.
Renger
,
Biophys. J.
95
,
105
(
2008
).
28.
K. L. M.
Lewis
,
F. D.
Fuller
,
J. A.
Myers
,
C. F.
Yocum
,
S.
Mukamel
,
D.
Abramavicius
, and
J. P.
Ogilvie
,
J. Phys. Chem. A
117
,
34
(
2013
).
29.
A.
Gelzinis
,
L.
Valkunas
,
F. D.
Fuller
,
J. P.
Ogilvie
,
S.
Mukamel
, and
D.
Abramavicius
,
New J. Phys.
15
,
075013
(
2013
).
30.
F.
Müh
,
M.
Plöckinger
, and
T.
Renger
,
J. Phys. Chem. Lett.
8
,
850
(
2017
).
31.
D.
Abramavicius
and
S.
Mukamel
,
J. Chem. Phys.
133
,
184501
(
2010
).
32.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics
14
,
33
(
1996
).
33.
L.
Zhang
,
D.-A.
Silva
,
H.
Zhang
,
A.
Yue
,
Y.
Yan
, and
X.
Huang
,
Nat. Commun.
5
,
4170
(
2014
).
34.
T. J.
Frankcombe
,
Phys. Chem. Chem. Phys.
17
,
3295
(
2015
).
35.
A.
Gelzinis
,
D.
Abramavicius
, and
L.
Valkunas
,
J. Chem. Phys.
142
,
154107
(
2015
).
36.
R. N.
Frese
,
M.
Germano
,
F. L.
de Weerd
,
I. H. M.
van Stokkum
,
A. Y.
Shkuropatov
,
V. A.
Shuvalov
,
H. J.
van Gorkom
,
R.
van Grondelle
, and
J. P.
Dekker
,
Biochemistry
42
,
9205
(
2003
).
37.
M.
Hoffmann
,
K.
Schmidt
,
T.
Fritz
,
T.
Hasche
,
V.
Agranovich
, and
K.
Leo
,
Chem. Phys.
258
,
73
(
2000
).
38.
N. J.
Hestand
,
R. V.
Kazantsev
,
A. S.
Weingarten
,
L. C.
Palmer
,
S. I.
Stupp
, and
F. C.
Spano
,
J. Am. Chem. Soc.
138
,
11762
(
2016
).
39.
R.
Tempelaar
and
D. R.
Reichman
,
J. Chem. Phys.
146
,
174703
(
2017
).
40.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
New York
,
1995
).
41.
J.
Rammer
,
Quantum Field Theory of Non-Eqilibrium States
(
Cambridge University Press
,
New York
,
2007
).
42.
V.
May
and
O.
Kühn
,
Charge and Energy Transfer Dynamics in Molecular Systems
, 2nd ed. (
Wiley-VCH
,
Weinheim
,
2011
).
43.
M.
Cho
,
Two-Dimensional Optical Spectroscopy
(
CRC Press
,
Boca Raton
,
2009
).
44.
T.-C.
Dinh
and
T.
Renger
,
J. Chem. Phys.
142
,
034104
(
2015
).
45.
S. S.
Andrews
,
J. Chem. Educ.
81
,
877
(
2004
).
46.
G.
Garab
and
H.
van Amerongen
,
Photosynth. Res.
101
,
135
(
2009
).
47.
D.
Abramavicius
and
S.
Mukamel
,
J. Chem. Phys.
122
,
134305
(
2005
).
48.
D.
Lindorfer
,
F.
Müh
, and
T.
Renger
,
Phys. Chem. Chem. Phys.
19
,
7524
(
2017
).
49.
G. U.
Bublitz
and
S. G.
Boxer
,
Annu. Rev. Phys. Chem.
48
,
213
(
1997
).
50.
O. J. G.
Somsen
,
V.
Chernyak
,
R. N.
Frese
,
R.
van Grondelle
, and
S.
Mukamel
,
J. Phys. Chem. B
102
,
8893
(
1998
).
51.
Y.
Shibata
,
S.
Nishi
,
K.
Kawakami
,
J.-R.
Shen
, and
T.
Renger
,
J. Am. Chem. Soc.
135
,
6903
(
2013
).
52.
J.
Adolphs
,
F.
Müh
,
M. E.-A.
Madjet
, and
T.
Renger
,
Photosynth. Res.
95
,
197
(
2008
).
53.
T.
Renger
and
F.
Müh
,
Photosynth. Res.
111
,
47
(
2012
).
54.
T.
Renger
,
M. E.-A.
Madjet
,
M.
Schmidt am Busch
,
J.
Adolphs
, and
F.
Müh
,
Photosynth. Res.
116
,
367
(
2013
).
55.
M. A.
Lomize
,
A. L.
Lomize
,
I. D.
Pogozheva
, and
H. I.
Mosberg
,
Bioinformatics
22
,
623
(
2006
).
56.
M.
Rätsep
,
J.
Linnanto
, and
A.
Freiberg
,
J. Chem. Phys.
130
,
194501
(
2009
).
57.
M.
Germano
,
A.
Shkuropatov
,
H.
Permentier
,
R.
Khatypov
,
V.
Shuvalov
,
A.
Hoff
, and
H.
Van Gorkom
,
Photosynth. Res.
64
,
189
(
2000
).
58.
M.
Germano
,
A. Y.
Shkuropatov
,
H.
Permentier
,
R.
de Wijn
,
A. J.
Hoff
,
V. A.
Shuvalov
, and
H. J.
van Gorkom
,
Biochemistry
40
,
11472
(
2001
).
59.
L.
Konermann
and
A. R.
Holzwarth
,
Biochemistry
35
,
829
(
1996
).
60.
E. G.
Andrizhiyevskaya
,
D.
Frolov
,
R.
van Grondelle
, and
J. P.
Dekker
,
Phys. Chem. Chem. Phys.
6
,
4810
(
2004
).
61.
F.
Müh
and
T.
Renger
, “
Structure-based calculation of pigment–protein and excitonic pigment–pigment coupling in photosynthetic light-harvesting complexes
,” in
The Biophysics of Photosynthesis
, edited by
J.
Golbeck
and
A.
van der Est
(
Springer
,
New York
,
2014
), pp.
3
44
.
62.
A.
Kell
,
X.
Feng
,
M.
Reppert
, and
R.
Jankowiak
,
J. Phys. Chem. B
117
,
7317
(
2013
).
63.
R Development Core Team
,
R: A Language and Environment for Statistical Computing
(
R Foundation for Statistical Computing
,
2015
).
64.
K. V.
Price
,
R. M.
Storn
, and
J. A.
Lampinen
,
Differential Evolution. A Practical Approach to Global Optimization
, Natural Computing Series (
Springer-Verlag
,
Berlin
,
2005
).
65.
D.
Ardia
,
K. M.
Mullen
,
B. G.
Peterson
, and
J.
Ulrich
, DEoptim: Differential Evolution in R (2015), version 2.2–3.
66.
C.
Eijckelhoff
,
F.
Vácha
,
R.
van Grondelle
,
J. P.
Dekker
, and
J.
Barber
,
Biochim. Biophys. Acta
1318
,
266
(
1997
).
67.
F.
Vácha
,
D. M.
Joseph
,
J. R.
Durrant
,
A.
Telfer
,
D. R.
Klug
, and
J.
Barber
,
Proc. Natl. Acad. Sci. U. S. A.
92
,
2929
(
1995
).
68.
F.
Vácha
,
M.
Durchan
, and
P.
Šiffel
,
Biochim. Biophys. Acta
1554
,
147
(
2002
).
69.
P. J.
van Kan
,
S. C.
Otte
,
F. A.
Kleinherenbrink
,
M. C.
Nieveen
,
T. J.
Aartsma
, and
H. J.
van Gorkom
,
Biochim. Biophys. Acta
1020
,
146
(
1990
).
70.
I. V.
Shelaev
,
F. E.
Gostev
,
V. A.
Nadtochenko
,
A. Y.
Shkuropatov
,
A. A.
Zabelin
,
M. D.
Mamedov
,
A. Y.
Semenov
,
O. M.
Sarkisov
, and
V. A.
Shuvalov
,
Photosynth. Res.
98
,
95
(
2008
).
71.
I.
Shelaev
,
F.
Gostev
,
M.
Vishnev
,
A.
Shkuropatov
,
V.
Ptushenko
,
M.
Mamedov
,
O.
Sarkisov
,
V.
Nadtochenko
,
A.
Semenov
, and
V.
Shuvalov
,
J. Photochem. Photobiol., B
104
,
44
(
2011
).
72.
K.
Acharya
,
B.
Neupane
,
V.
Zazubovich
,
R. T.
Sayre
,
R.
Picorel
,
M.
Seibert
, and
R.
Jankowiak
,
J. Phys. Chem. B
116
,
3890
(
2012
).
73.
T.
Renger
and
R. A.
Marcus
,
J. Chem. Phys.
116
,
9997
(
2002
).
74.
J.
Ma
and
J.
Cao
,
J. Chem. Phys.
142
,
094106
(
2015
).
75.
T.
Mančal
,
L.
Valkunas
, and
G. R.
Fleming
,
Chem. Phys. Lett.
432
,
301
(
2006
).
76.
S.
Valleau
,
A.
Eisfeld
, and
A.
Aspuru-Guzik
,
J. Chem. Phys.
137
,
224103
(
2012
).
77.
F. D.
Fuller
,
J.
Pan
,
A.
Gelzinis
,
V.
Butkus
,
S. S.
Senlik
,
D. E.
Wilcox
,
C. F.
Yocum
,
L.
Valkunas
,
D.
Abramavicius
, and
J. P.
Ogilvie
,
Nat. Chem.
6
,
706
(
2014
).
78.
E.
Romero
,
R.
Augulis
,
V. I.
Novoderezhkin
,
M.
Ferretti
,
J.
Thieme
,
D.
Zigmantas
, and
R. V.
Grondelle
,
Nat. Phys.
10
,
676
(
2014
).
79.
V.
Butkus
,
L.
Valkunas
, and
D.
Abramavicius
,
J. Chem. Phys.
140
,
034306
(
2014
).
80.
L.
Chen
,
R.
Zheng
,
Q.
Shi
, and
Y. J.
Yan
,
J. Chem. Phys.
131
,
094502
(
2009
).