Photophysical Properties of Natural Light-Harvesting Complexes

Jan 31, 2008 - Anjali Pandit , Francesco Buda , Adriaan J. van Gammeren , Swapna Ganapathy and Huub J. M. de Groot. The Journal of Physical Chemistry ...
0 downloads 0 Views 345KB Size
Download PDF
J. Phys. Chem. B 2008, 112, 2207-2217

2207

Photophysical Properties of Natural Light-Harvesting Complexes Studied by Subsystem Density Functional Theory Johannes Neugebauer* Laboratorium fu¨r Physikalische Chemie, ETH Zurich, Wolfgang-Pauli-Strasse 10, 8093 Zurich, Switzerland ReceiVed: October 12, 2007; In Final Form: NoVember 26, 2007

In this study, we investigate the excited states and absorption spectra of a natural light-harvesting system by means of subsystem density functional theory. In systems of this type, both specific interactions of the pigments with surrounding protein side chains as well as excitation energy transfer (EET) couplings resulting from the aggregation behavior of the chromophores modify the photophysical properties of the individual pigment molecules. It is shown that the recently proposed approximate scheme (J. Chem. Phys. 2007, 126, 134116) for coupled excitations within a subsystem approach to time-dependent DFT is capable of describing both effects in a consistent manner, and is efficient enough to study even the large assemblies of chromophores occurring in the light-harvesting complex 2 (LH2) of the purple bacterium Rhodopseudomonas acidophila. A way to extract phenomenological coupling constants as used in model calculations on EET rates is outlined. The resulting EET coupling constants and spectral properties are in reasonable agreement with the available reference data. Possible problems related to the effective exchange-correlation kernel are discussed.

1. Introduction Light-harvesting complexes accomplish the first two of the essential steps in photosynthesis: They are responsible for the absorption of light and the transport of the excitation energy to the photosynthetic reaction center, where the subsequent steps, i.e., charge separation and chemical storage of the energy, take place.1-3 A key step toward the understanding of the processes in light-harvesting complexes on a molecular level was the crystal structure determination of the integral membrane lightharvesting complexes from purple bacteria (for reviews, see refs 4 and 5), e.g., from Rhodopseudomonas molischianum6 and from Rhodopseudomonas acidophila.7-9 In the present study, we will concentrate on the lightharvesting complex 2 (LH2) of the latter bacteria, for which extensive studies of pigment properties and excitation energy transfer processes have been conducted.10-19 To investigate the energy transfer pathways in light-harvesting systems, it is necessary to calculate the excitation energy couplings between the individual chromophores.20-22 In principle, this information should be accessible from a quantum-chemical calculation of the excited states of the entire light-harvesting system, but there are two basic problems. The first is the size of the antenna systems, which issdepending on the particular model chosens in the order of several thousands to ten thousands of atoms. The second is the fact that this type of analysis necessarily requires a picture in which the excited states of the pigment complex are expressed in terms of local excited states of the individual chromophores. A quantum-chemical calculation, however, would typically result in a delocalized picture, which prohibits (or at least severely complicates) the extraction of such coupling constants. For these two reasons, the investigation of excitonic couplings in antenna systems is usually carried out by calculating excited states of individual, isolated chromophores. Subsequently, the * E-mail: [email protected].

coupling constants are determined by calculating the interaction matrix elements between the excited states of the monomers, which effectively corresponds to an interaction of the transition densities of the isolated monomers. With these coupling matrix elements, an effective Hamiltonian matrix can be set up from which energy levels of the light-harvesting system can be extracted upon diagonalization.23 The coupling matrix elements are often approximated by considering only Coulomb coupling, for which further approximations can be introduced. Among the methods which are most often used are multipole expansions, e.g., the Fo¨rster-type dipole-dipole coupling,24 transition monopole interactions in terms of atomic partial transition charges,25,26 or a real-space partitioning of the transition densities and a subsequent numerical integration of the Coulombic interactions as in the transition density cube method.14 For a comparison of some of these approaches, see ref 27. More recently, density functional and configuration-interaction-type methods based on an exciton-like interaction between two local electronic transitions on isolated pigments have been suggested, which may also incorporate the effect of a dielectric medium on the coupling.28-31 In view of the structure of the light-harvesting complex, which clearly contains specific interactions between the chromophores and their environment, e.g., the apoproteins surrounding the pigments,32 it appears desirable to be able to study the electronic excitations of the antenna system including such interactions by an atomistic, quantum-chemical model. As was demonstrated in recent work, the frozen-density embedding (FDE) method,33 which is a subsystem approach to density functional theory (DFT), and its time-dependent response generalization34 are efficient approaches to study chromophores in condensed phases35-37 or in complexes with specific interactions.38,39 It is also possible to include excitation energy transfer (EET) couplings between different pigment molecules in a related subsystem-based time-dependent density functional theory approach40 that is based on the extension of FDE to time-dependent density functional theory (TDDFT).34

10.1021/jp709956k CCC: $40.75 © 2008 American Chemical Society Published on Web 01/31/2008

2208 J. Phys. Chem. B, Vol. 112, No. 7, 2008

Figure 1. Structure of LH2 of Rhodopseudomonas acidophila.

Therefore, we are going to test the ability of this subsystem approach for the description of specific pigment-environment and excitation energy transfer interactions in the following. First, a brief summary of the structural features of LH2 is given in section 2 in order to make this work self-contained and to provide a basis for the discussion of specific effects presented in subsequent sections. In section 3, the key steps for the determination of excitation energies within the subsystemTDDFT approach will be reviewed and efficiency considerations will be presented. Moreover, that section explains how the coupling matrix elements obtained in the TDDFT treatment can be related to the phenomenological coupling constants used in model theories of EET. Technical details of the calculations are presented in section 4. Section 5 presents results for specific pigment-environment interactions on local excitations; couplings arising from an interaction of excitations on different subsystems are subsequently investigated in section 6 before the calculated absorption spectra are presented in section 7. Finally, conclusions are given in section 8. 2. The Structure of LH2 The structure of LH2 from Rhodopseudomonas acidophila has been described in detail before,7-9 so that only a brief summary of the most important features relevant for the present study will be given here. A crystal structure obtained at 2.0 Å resolution which was reported in ref 8 showed that this antenna system has, within the experimental accuracy, C9 symmetry. The structure is shown in Figure 1. It contains nine bacteriochlorophyll a (Bchl a) molecules forming a ring in such a way that the planes of the macrocyles are aligned parallel to the membrane. This set of Bchl a molecules is called the B800 system, since its absorption maximum is at a wavelength of 800 nm. Additionally, there is a second ring of 18 Bchl a molecules forming the so-called B850 system, in which the macrocycle planes are perpendicular to the membrane. These Bchl a molecules occur in pairs of so-called R- and βB850 molecules. The Mg-Mg distance in B800 is 21.2 Å, and the distances in B850 are 9.0 Å within an (R,β)-dimer and 9.5 Å between different dimers. There is an inner ring of so-called R-apoproteins, and an outer ring of β-apoproteins. The Mg atoms

Neugebauer in B800 interact with an oxygen atom from the COO--modified R-Met1 residue, while the Mg atoms in B850 ligate nitrogen atoms from R-His31 or β-His30 residues, respectively. These interactions may lead to shifts in the absorption properties of the Bchl a molecules, although this effect is typically small.17 The configuration-interaction singles (CIS) calculations in ref 15 show that a red-shift of about 0.02 eV results from the His residues, although this calculation was carried out with a very small basis set (3-21G*). Additionally, the residues R-Trp45 and R-Tyr44 can form hydrogen bonds to the ring I acetyl groups of RB850 and βB850, respectively. In ref 17, it was reported that this causes a measurable red-shift in the absorption spectra of the light-harvesting complex. The B800 Bchl a molecules can form hydrogen bonds to the β-Arg20 residue, for which a significant red-shift was calculated in ref 12. The asymmetric unit of the antenna system is called a protomer and consists of an R- and a β-apoprotein together with one B800 molecule and a B850 dimer as well as a carotenoid (rhodopin glucoside, RG1). In ref 8, part of the electron density was assigned to a second carotenoid, which was assumed to exhibit cis-double bonds in contrast to the all-trans-structure of RG1. In later work, this assignment was reversed.7 In sections 5 and 6, it will be demonstrated that a detailed understanding of the spectral features of the LH2 antenna system is possible on the basis of FDE calculations. We will concentrate on the interactions of the Bchl a pigments with their direct environment, as well as on EET coupling interactions between these chromophores, whereas the RG1 pigments are not explicitly taken into account in this study (see, e.g., ref 41 for an investigation on their role in excitation energy transfer processes). 3. Methodology In ref 40, it was shown that excitation energies within a subsystem approach to TDDFT can be obtained from the eigenvalue problem

[Ω ˜ - ωk2]F ˜k ) 0

(1)

which is similar to conventional TDDFT,42-44 but the matrix Ω ˜ has a subsystem structure in this case and contains contributions arising from the nonadditive kinetic-energy functional employed in FDE and related subsystem DFT approaches (see, e.g., refs 33 and 45-47). In that subsystem structure, Ω ˜ effectively couples local excitations of the interacting fragments rather than individual orbital transitions as in conventional TDDFT. These couplings are often small or even negligible when the excitations under study have local character. In such cases, good results can be obtained without including intersystem couplings.35-38 Critical cases are systems in which identical chromophores are present, or when (near-)degenerate states occur. Because of this, good approximations for coupled excitation energies can be obtained by diagonalizing subblocks of Ω ˜ which only include the couplings between these critical transitions, as was demonstrated in ref 40. The most important approximations introduced in that work are thus (i) the restriction to subsets of excitations for the subsystems, which can be obtained as in conventional TDDFT calculations by means of a Davidson-type diagonalization, and additionally (ii) the use of density fitting techniques to calculate the potential induced by a transition density of a local excitation. However, to make the approach work for larger aggregates, some further improvements are introduced here, which make the calculation of the elements of Ω ˜ much more efficient without significant effect

Properties of Natural Light-Harvesting Complexes

J. Phys. Chem. B, Vol. 112, No. 7, 2008 2209

Figure 2. Illustration of the subsystem grid technique for the calculation of coupling constants within the subsystem-TDDFT approach.

on the numerical accuracy. These improvements are described in detail below. The scheme presented in the following is implemented in a locally modified version of the ADF program package.49,50 Most of the actual computational effort is spent on the integration of the matrix elements

Ω ˜ µAνB )

∫dr1 ∑

(iaσ)A

x

2U(iaσ)AµA ω(iaσ)AφiAσ(r1)VνelB(r1)φaAσ(r1) (2)

where we adopted the notation from ref 40. In particular, U(iaσ)AµA are elements of the eigenvector for transition µA on subsystem A obtained in the (uncoupled) subsystem-TDDFT calculations, φ are occupied (index i) or virtual (index a) orbitals of a given subsystem, and VνelB is the induced potential due to transition ν located on subsystem B, which is calculated analytically from fitted transition densities. These matrix elements are calculated in ADF by numerical integration, i.e.,

Ω ˜ µAνB ≈

∑k w(rk) ∑

(iaσ)A

x

2U(iaσ)AµA ω(iaσ)AφiAσ(rk)VνelB(rk)φaAσ(rk) (3)

where rk and w(rk) are the coordinates and weights of the grid points, respectively. From this, it becomes apparent that the coupling matrix elements for a certain column µA in the matrix Ω ˜ will only require an accurate integration grid in the region where the molecular orbitals (MOs) of system A are nonvanishing. Since the MOs are expanded in basis functions of the subsystem, the numerical grid must be accurate in regions where the basis functions of the subsystem are nonvanishing. However, we usually have to calculate coupling matrix elements between electronic transitions of all the subsystems, and the numerical grid must be suited for all of them. This would correspond to the default grid construction used in ADF and would cover the entire supersystem for which the calculation is carried out. For the subsystem-TDDFT approach used here, we can significantly reduce the computational cost by using different subsystem grids for the coupling matrix elements involving different subsystems, so that the number of integration points per coupling matrix element is indeed small and independent of the size of the surrounding system (see Figure 2 for a schematic representation). Hence, the size of the integration grid can greatly be reduced without a significant effect on the numerical integration accuracy. Since the computer time needed scales linearly with the number of integration points, this results in a substantial speed-up of the calculations. This was imple-

mented in the subsystem-TDDFT module of a locally modified version of ADF. A similar effect could probably be achieved by applying distance-based cut-offs in the numerical integration as used in linear-scaling techniques. However, by constructing subsystem grids, we can avoid even the effort of constructing parts of the grid that are not important for the matrix elements and of performing distance checks in order to determine blocks of irrelevant integration points. Note that also the recently presented ground-state FDE implementation in ADF, in which the densities of several subsystems can be optimized in one run, allows the use of subsystem grids optimized for each particular subsystem.47 It should be noted that the numerical integration is very well suited for parallelization, and full advantage of this is taken in the calculation of the coupling constants. Additionally, a prescreening for the elements U(iaσ)Aµ of the subsystem eigenvectors is carried out in order to reduce the number of occupied-virtual pairs that have to be taken into account in the summation, which also increases the speed of the calculation while allowing for full control over the numerical error introduced. In this way, an efficient calculation of the excitation energies and transition dipole moments is possible for systems that are composed of many subunits, even if the individual subunits are already comparatively large (>100 atoms). Nevertheless, the limiting factor for the calculations is in general the size of the largest subsystem in the calculation. Another feature that is used in ground-state FDE calculations is the superposition of subsystem densities in order to approximate the density of the full system.35 In particular, it is possible to copy the density of a particular fragment to several positions in space in case of identical subsystems.36 This can also be applied for the present subsystem-TDDFT approach and has been extended so that also transition densities and (electric and magnetic) transition dipole moments can be used for multiple fragments. The cost for the preparation steps of the subsystem calculation is thus significantly reduced for systems in which several identical chromophores occur. Coupling constants as used in phenomenological excitonic coupling schemes are usually calculated as matrix elements between local excited states. With these coupling matrix elements, a configuration-interaction-like matrix can be set up, and excited-state energies of the total system are found by diagonalization of that matrix (see, e.g., ref 51). The matrix eigenvalue problem solved in the context of TDDFT is, however, somewhat different, as it effectively involves coupling matrix elements between (energy-weighted) transition densities, and it yields squared excitation energies instead of excited-state energies. For better comparison with data from the literature, it is thus desirable to establish a relationship between the two different sets of coupling constants or to calculate CI-like

2210 J. Phys. Chem. B, Vol. 112, No. 7, 2008

Neugebauer

coupling constants from the TDDFT data. This is possible as follows: Consider a pair of uncoupled excited states with energies Ea,b. In a CI-like case, the coupling between the two states, Vab, will lead to new energy levels

E+,- )

Ea + Eb ( 2

x(

)

Ea - Eb 2 + Vab2 2

(4)

The energy levels will, of course, change if couplings to other local excited states are included, although for (near-)degenerate states typically one single coupling is dominant. Nevertheless, this means that Vab2 can be obtained from the energies of the two-state problem as

Vab2 )

(

) (

)

E+ - E - 2 E a - Eb 2 2

2

(5)

The energy differences can also be expressed in terms of excitation energy differences ω, which are the quantities occurring in the TDDFT context

Vab2 )

(

) (

)

ω + - ω- 2 ω a - ωb 2 2

2

(6)

The excitation energies ωa,b of the local (subsystem) excitations are known in the calculation, and the coupled excitation energies ω+,- can be obtained from the eigenvalues of the corresponding 2 × 2 subblock of the coupling matrix. They can also be expressed directly in terms of the coupling matrix elements Ω ˜ ab, leading to (cf. the structure of Ω ˜ in ref 40)

ω + - ω- )

x x

ωa2 + ωb2 + 2

x( x(

ωa2 + ωb2 2

)

ωa2 - ωb2 2 +Ω ˜ ab2 2

)

ωa2 - ωb2 2 +Ω ˜ ab2 (7) 2

In combination with eq 6, this can be used to convert the TDDFT coupling matrix elements Ω ˜ ab into CI-like excitonic coupling integrals Vab. The information about the sign is taken from Ω ˜ ab, since only the relative signs of the matrix elements play a role, which are the same for matrix elements between excited states or their corresponding transition densities. It should be noted that the relationship above only guarantees the same energy gaps between the two excited states but not exactly the same excited-state energies. This is because the coupling in a CI context creates a symmetrical splitting of the energy around the mean value of the locally excited states, (Ea + Eb)/ 2, while the coupling in a TDDFT context induces a splitting that is symmetrical for the square of the excitation energy, ω2, around the mean value of the squared local excitation energies, (ωa2 + ωb2)/2. By diagonalizing the CI-like matrix constructed from such phenomenological coupling constants and the excitation energies, it turns out, however, that the differences between the original TDDFT results and those from a CI-like treatment are usually very small (in the order of 0.001 eV). It should be noted that the definition of the couplings on the basis of energy differences, i.e., following eq 4, is consistent with other DFT-based approaches,30,31 in particular for the case of identical monomer excitation energies, which is stated explicitly in ref 31. Equations 6 and 7 implement this definition in the context of the FDE scheme used here. A difference exists, however, in how the interactions between the monomers are taken into account, as has been explained in ref 40 in more

detail and will briefly be analyzed again in section 5: The approaches proposed in refs 30 and 31 evaluate the couplings starting from isolated monomers, whereas the method employed here is integrated into the FDE formalism, so that already the monomer properties and excitation energies are calculated in the presence of all other monomers included in the calculation. Thus, the starting point for the evaluation of the coupling constants are the (hypothetical) excitation energies of the monomers in the aggregate with excitonic coupling switched off. In a recent study based on the approaches in refs 30 and 31, this effect was approximately taken into account by means of a set of point charges that were fitted to the other monomers’ charge distributions48 and was found to be rather small in that particular case. In contrast to such an approximate treatment, FDE can be regarded as a method that is in principle exact for the description of the ground-state interactions between different molecules (in the limit of exact functionals, see ref 45), and its local response variant has proven to yield good results in many cases where the effects of such interactions on the absorption properties are dominant.35-38 Furthermore, the calculation of the interaction between the local excitations is carried out in a slightly different way compared to the other DFT-based schemes. Whereas Coulomb and exchange-correlation couplings are treated in formally the same way (although this may not hold for the approximation employed for the exchangecorrelation kernel), overlap effects are treated differently. Within the FDE framework, the nonadditive kinetic-energy functional and its functional derivatives that enter the embedding potential and effective kernel should, in principle, take care of nonorthogonality effects, whereas an explicitly overlap-dependent term is included in the perturbative treatment of ref 31. When comparing these effects, it must be considered that already the subsystem orbitals and transition densities in the case of FDE are influenced by the kinetic-energy functional, and that only approximate kinetic-energy functionals can be applied. In practice, however, such effects are usually very small.40 In ref 30, the overlap contribution to the coupling for the B1u transition of the ethylene dimer was calculated explicitly and found to be roughly 4 orders of magnitude smaller than the Coulomb contribution and 2 orders of magnitude smaller than the exchange-correlation contribution at the rather short intermolecular distance of 4 Å. An additional short-range effect that might be more important for the splitting energies is the interaction with charge-transfer-like configurations (see, e.g., the example in ref 40). Neither the approach originally suggested in ref 31 nor the coupled FDE approach used here includes such effects at the present stage. Finally, it should be noted that the approaches in refs 30 and 31 focused mainly on the inclusion of condensed-matter effects on excitation energy transfer phenomena by means of continuum models (see also the discussion in the Conclusion section). 4. Technical Details All calculations presented in this work are performed with a locally modified version of the ADF program package.49,50 If not indicated otherwise, the following procedure has been applied for structural motifs of LH2: Coordinates for nonhydrogen atoms were taken from the X-ray structure reported in ref 8 (PDB entry 1NKZ). Hydrogen atoms were added using default bond angles and bond lengths. Subsequently, the positions of the hydrogen atoms in the individual pigment molecules were optimized using the Becke-Perdew functional denoted as BP8652,53 and the TZP basis set from the ADF basis set library. The positions of the non-hydrogen atoms were kept

Properties of Natural Light-Harvesting Complexes

J. Phys. Chem. B, Vol. 112, No. 7, 2008 2211

fixed. For the calculation of excitation energies and excitation energy couplings, we proceed as follows: The ground-state densities and molecular orbitals of the pigment molecules and the surrounding fragments under investigation are determined from a frozen-density embedding calculation33 for each subsystem, including three freeze-and-thaw cycles.54 For the pigments, we calculated the local excitations within the frozendensity embedding approach as described in ref 38 which were subsequently coupled according to the approach presented in ref 40. The TZP basis set in combination with the “statistical averaging of (model) orbital potentials’’ (SAOP) potential55-57 is used in all excitation calculations and frozen-density preparation calculations. The exchange-correlation component of the embedding potential is usually approximated using the same potential as that used for the embedded subsystem. Since we want to make use of the density of the environment only, we employ the exchange functional by Becke52 and the correlation functional by Perdew and Wang58 (PW91) for the embedding potential; the SAOP potential is orbital-dependent and can thus not be used for such an embedding potential. The kinetic-energy component of the embedding potential is approximated employing the functional denoted as PW91k in ref 35, which is based on the PW91 exchange functional with the reparameterization by Lembarki and Chermette.59 The choice of this functional follows from the investigation in ref 60. It has proven to yield reliable results in a large number of applications of the FDE scheme for molecular properties.35-39 For the spectra plots, we converted the calculated oscillator strengths f into molar absorption coefficients (ν) according to61

f)

(

)∫

4mec0 ln(10) NAe2

band

(ν) dν

TABLE 1: Calculated Excitation Energies (SAOP/TZP; in Units of eV) of Bchl a for Different Structures and Environmentsa structure RB850 Bchl a (fully opt.) RB850 Bchl a (H opt.) RB850 Bchl a (H opt.) RB850 Bchl a (H opt.) RB850 Bchl a (H opt.*) βB850 Bchl a (fully opt.) βB850 Bchl a (H opt.) βB850 Bchl a (H opt.) βB850 Bchl a (H opt.) βB850 Bchl a (H opt.*) B800 Bchl a (fully opt.) B800 Bchl a (H opt.) B800 Bchl a (H opt.) B800 Bchl a (H opt.) B800 Bchl a (H opt.*) exp. Bchl a 2

environment

B850 R-His31 R-Trp45

B850 β-His31 R′-Tyr44

B800 R-Met1 β-Arg20

Qy

Qx

1.76 1.59 1.60 1.60 1.59 1.77 1.57 1.58 1.58 1.57 1.79 1.61 1.60 1.61 1.55 1.60

2.02 1.92 1.94 1.88 1.91 2.02 1.95 1.97 1.92 1.93 2.07 1.95 1.95 1.90 1.85 2.15

a The starting point for all (partial) structure optimizations carried out here was the X-ray structure with 2.0 Å resolution from ref 8, PDB code 1NKZ. H opt.: hydrogen atoms added and optimized (BP86/TZP; if indicated with an asterisk, also atoms involved in hydrogen bonds have been optimized); fully opt.: hydrogen atoms added, all atoms optimized (BP86/TZP). Environments: B850/B800: all other B850/ B800 Bchl a molecules included. Model compounds of different size were used for the amino acid residues R-His31, R-Trp45, β-His31, R′-Tyr44, R-Met1, and β-Arg20 (see text for details). Experimental values for Bchl a in diethyl ether were taken from ref 2.

(8)

where the factor ln(10) arises if the absorption coefficient is defined in terms of the decadic logarithm, ν denotes the frequency, me is the electron mass, c is the speed of light, 0 is the vacuum permittivity, and A ) ∫band (ν) dν is the absorption coefficient integrated over the absorption band. The absorption coefficient (ν) was modeled by a Gaussian function of halfwidth 0.01 eV. Graphics of the molecular structures were generated with the programs Molden62 and VMD.63 5. Chromophore-Environment Interactions The absorption properties of a natural light-harvesting complex will depend on several factors, involving (i) the absorption properties of the (isolated) pigment molecules, (ii) the environmental effects on the subsystem properties (e.g., axial ligand, protein side chains, neighboring pigments, etc.), and (iii) phenomena related to the aggregation of several chromophores (i.e., excitation energy couplings). In this section, the environmental effects on local fragment properties will be studied by means of the (uncoupled) FDE approach. The following effects on the excitation energies of Bchl a will be analyzed: (a) the structural change from an optimized structure to the crystal structure, (b) the influence of neighboring Bchl a molecules in the aggregate, (c) the effect of axial ligands, and (d) the effect of hydrogen-bonding ligands. In the latter three cases, the neighboring molecule or ligand will be represented in terms of an effective embedding potential obtained from its frozen density; these densities are polarized in all cases by means of freeze-and-thaw iterations (details are described below). The absorption spectrum of a Bchl a molecule consists of three prominent bands: the intense Qy band in the region of

Figure 3. Structure of the Bchl a monomers investigated in this study. All structures were obtained from the X-ray structure from ref 8, PDB code 1NKZ. Positions of added hydrogen atoms have been optimized (BP86/TZP).

750-800 nm with a maximum at 773, the much less intense Qx band between 550 and 600 nm (maximum at 577 nm), and the intense Soret band with a maximum at 358 nm, which extends from 0.4 eV), which might be a hint on the latter effect. On the other hand, only relatively small basis sets were used in that study, so that the energy differences might be subject to changes when larger basis sets are used in connection with the B3LYP functional. In order to analyze whether the neighboring Bchl a molecules have an impact on the (uncoupled) excitation energies of the pigments, we investigated the B850 unit of LH2 by means of uncoupled FDE calculations, as described in detail in ref 35. In these calculations, the effect of the surrounding molecules is incorporated by an effective embedding potential constructed from the environmental density. This density surrounding either an RB850 or a βB850 Bchl a molecule was approximated in the following way: We extracted a tetrameric unit from B850, which is shown in Figure 4. Then, the electron densities were calculated (SAOP/TZP) for each of the monomers within the tetramer without including effects of neighboring molecules. Subsequently, the electron densities for all four monomers were relaxed in a series of ground-state FDE calculations, in which a superposition of the densities of the other three subsystems was employed to model the environment. In each step, the latest approximation available for the density of each system in the environment was used. This cyclic optimization was carried out three times in total. Since the adjacent Bchl a molecules will have a much stronger impact on the electron density of a particular pigment than more distant monomers, the electron densities thus obtained for the inner two Bchl a molecules of the tetramer should already represent a good approximation to the monomer density in the full B850 system. In a last step, we used the densities for these inner two monomers of the tetramer to set up an approximate density of the full B850 system by

Neugebauer

Figure 4. Structure of the tetrameric Bchl a substructure of the B850 unit. The labels A and B for the (R,β)-dimers have been introduced for comparability with the results obtained in ref 14, where the same notation is used.

superposition of these relaxed monomer densities, copied to the positions of all monomers in B850. With this environmental density, we calculated the electron densities of one RB850 and one βB850 Bchl a in the full B850 complex as well as their uncoupled excitation energies in a FDE calculation. The total system in these calculations consists of 2520 atoms, of which 2380 belong to the frozen environment and 140 to the embedded system. It should be noted that this way of approximating the environmental density in an uncoupled FDE calculation is much more sophisticated than, e.g., a simple superposition of nonrelaxed subsystem densities, which nevertheless often yields reasonable results.35,36 The effect of the other Bchl a molecules within the B850 unit on the uncoupled excitation energies (sometimes called “site energies’’) is rather small. The excitation energies of the Qy and Qx transitions of both RB850 and βB850 Bchl a are shifted by only 0.01 and 0.02 eV, respectively, when this part of the environment is included (see Table 1). A similar test was carried out for the B800 system, but because of the larger distance between the monomers, we only performed one relaxation step of the monomer density with respect to the density of all other monomers. Also here, the excitation energies for Qy and Qx transitions are virtually unchanged (the change in the Qy excitation energy is in the order of 0.001 eV but leads to a roundoff effect). In a next step, we investigated the influence of an axial ligand on the (otherwise isolated) Bchl a pigments. For RB850 and βB850, this axial ligand is a histidine residue, which was modeled by the frozen density of an imidazole ring in a FDE calculation. Also, this axial ligand results in only a small change. For RB850 Bchl a, the Qy excitation energy increases by 0.01 eV, and the Qx excitation energy decreases by 0.04 eV. Similar changes are found for βB850 Bchl a. This is consistent with the observations made before that axial ligands only lead to slight modulations of the site energies of Bchl a molecules.17 In the case of B800 Bchl a, the axial ligand reported in refs 7 and 8, which was also used in our study, is a COO--modified R-Met1; i.e., the coordinating group is a carboxy group. This is in contrast to the Met ligand that was assumed in the earlier study reported in ref 12, where the effect of the protein environment on the excited states of Bchl a in B800 was investigated. In this case, we included also the following residues (Asn2-Gln3-Gly4) into the part of the protein modeled by FDE in order to test the influence of a larger protein model. However, also here, the effect is rather small. The Qy excitation energy remains unchanged at 1.61 eV, while the Qx transition decreases by 0.05 eV.

Properties of Natural Light-Harvesting Complexes Finally, the effects of the hydrogen-bonding residues R-Trp45 (RB850), R-Tyr44 (βB850), and β-Arg20 (B800) were investigated. These amino acids were modeled by their side chains only (saturated with hydrogen). Again, the positions of the hydrogen atoms were optimized in the monomers from which the hydrogen-bonded complex was created. Afterward, the positions of the atoms involved in the hydrogen bond were optimized in the hydrogen-bonded complex. For the structures obtained in this way, we performed uncoupled FDE calculations considering the amino acid side chains as the frozen part. For the B800 Bchl a pigment, we find rather large red-shifts of 0.06 and 0.10 eV for the Qy and Qx transitions in comparison to the isolated molecule. This is in good agreement with the corresponding shifts of 0.05 and 0.08 eV obtained from supermolecular B3LYP/6-31G calculations in ref 12. In contrast to this, rather small changes (e0.02eV) were found for RB850 and βB850 in comparison to the isolated pigment molecules. In all of the calculations reported above, the effect of either neighboring molecules, axial ligands, or hydrogen-bonding ligands was included by means of an effective embedding potential in an (uncoupled) FDE calculation. In order to assess the error introduced by the frozen-density approach, the calculation on RB850 Bchl a including the Trp45 ligand was repeated in a supermolecular Kohn-Sham (TD)DFT calculation. The resulting excitation energies for the Qy and Qx transitions agree within 0.01 eV with the FDE values, which shows that FDE does not introduce significant errors for these excitations. 6. Chromophore-Chromophore Couplings The next step in the calculation of the absorption properties of LH2 is to determine the coupled excitation energies. As outlined in section 3, this is often done by calculating CI-like matrix elements between excited monomer states, which are then used to construct an effective Hamiltonian matrix from which coupled excitations can be obtained. In order to validate the coupling constants obtained here by comparison to previous results, we will in this section concentrate on the largest coupling constants between neighboring pigments in the light-harvesting complex, whereas all couplings occurring between the Qy and Qx transitions are implicitly incorporated in the resulting spectra presented in section 7. Couplings were calculated for the combined B850 and B800 units (3780 atoms in total) of LH2, which were separately prepared as described in the last section and then combined (i.e., no additional uncoupled FDE calculations were performed). The intra-subunit couplings obtained by considering the B850 and B800 units independently are virtually the same and are therefore not explicitly mentioned here. The pigments are labeled with a letter (A, B) indicating the protomer unit; i.e., an A-A coupling is an intrapolypeptide coupling in the nomenclature of ref 15, whereas an A-B coupling is an interpolypeptide coupling. The results for the strongest couplings between pigments within the B850 subunit are compared to the (Coulomb-only) transition density cube (TDC) data based on scaled CIS transition densities from ref 14 as well as to the CIS results from ref 15 in Table 2. The first observation that can be made is that the coupling constants calculated here are in reasonable agreement with the ones obtained by Krueger et al.14 and the scaled Coulomb results by Scholes et al.15 This is due to the fact that the couplings are dominated by the Coulomb contribution, which was scaled in refs 14 and 15 to match the experimental transition dipole moment. Scholes et al. assumed an experimental dipole moment of 2.51 au (6.39 D), whereas their CIS/6-31G* calculation resulted in a transition moment

J. Phys. Chem. B, Vol. 112, No. 7, 2008 2213 TABLE 2: Calculated (SAOP/TZP; in Units of cm-1) Coupling Constants for the Qy Transition of (r,β)B850 Bchl a in a Combined B850-B800 Complex (FDEc, Coupled FDE Calculation)a type FDEc, B850 + B800 TDC (scaled Coulomb only)14 CIS15 CIS (scaled Coulomb)15

RB T βA RB T βB RB T RA βB T βA 200 213 550 255

232 238 730 320

-65 -46

-45 -37

a For comparison, also the results from ref 14 obtained from a Coulomb-only coupling with the transition density cube (TDC) method on the basis of scaled CIS transition densities as well as the CIS results (both original data and with scaled Coulomb contribution) from ref 15 are presented.

of 4.03 and 3.98 au for RB850 and βB850 Bchl a, respectively, so that their calculated Coulomb contribution was empirically scaled down to about 40% of the original value. The transition dipole moment calculated in this work (SAOP/TZP) is 2.72 au for RB850 and 2.71 au for βB850 Bchl a in the isolated calculation and decreases slightly to 2.67 for both RB850 and βB850 Bchl a in the uncoupled FDE calculation on the full B850 system, thus demonstrating the reliability of transition dipole moments from TDDFT calculations. The data from ref 15 suggest that short-range contributions to the total coupling are small though not negligible. The shortrange contributions given there are 55 and 60 cm-1, respectively, for couplings of type RB T βB and RB T βA, respectively. However, these couplings were also found to be very sensitive to the inclusion of amino acids with hydrogen-bonding ability, which could not be explained on the basis of their calculation. Additionally, the Coulomb contributions in the results from ref 15 were empirically scaled down as mentioned above, whereas the short-range contribution was not modified. It might therefore be possible that the short-range contribution was overestimated in that work. In the case of coupled frozen-density embedding calculations, short-range contributions are in principle included, although they cannot directly be mapped to the corresponding short-range contributions in a CIS-like treatment or in other DFT-based approaches.30,31 The ground-state calculations in the FDE case include interactions between the monomers, so that the monomer orbitals and densities are already polarized with respect to the surrounding molecules in the aggregate. This can of course lead to changes in the coupling strength, e.g., in a modified Coulomb coupling due to changes in the orbitals. Furthermore, the effective kernel employed in the (uncoupled and coupled) subsystem-TDDFT approach contains both exchange-correlation (XC) and kinetic-energy contributions originating from the effective embedding potential,40 which also give rise to shortrange contributions. However, the current implementation for the coupled FDE approach employs the adiabatic local density approximation (ALDA) for both XC and kinetic-energy contributions, which are usually small and of opposite sign (see ref 40 for an example). In order to analyze the role of shortrange couplings in the results obtained here in more detail, we extracted an (R,β)Bchl a dimer from the B850 subunit, which shows the largest couplings in LH2. The coupling constant obtained for this dimer was 234 cm-1 and thus very similar to the corresponding value in the full B850 + B800 system. The XC and kinetic-energy contributions in the coupled FDE calculation were analyzed by switching the corresponding terms in the kernel off in the dimer calculation. If the kinetic-energy contribution is omitted, the coupling reduces from 234 to 230 cm-1. If additionally the XC contribution is omitted (pure Coulomb coupling), the total coupling increases again to 233

2214 J. Phys. Chem. B, Vol. 112, No. 7, 2008

Neugebauer

TABLE 3: Excitation Energies E (SAOP/TZP; in Units of eV) and Oscillator Strengths f for the Qy and Qx Transitions of a Bchl a Model Dimer (Phytyl Chain Omitted) from Coupled FDE (FDEc) and Supermolecular TDDFT Calculations (super)a calculation

E(Qy)

f(Qy)

E(Qy)

f(Qy)

FDEu (R) FDEu (β) FDEc (R + β)

1.585 1.563 1.541 1.606 1.530 1.596

0.298 0.293 0.552 0.043 0.574 0.043

1.911 1.940 1.907 1.943 1.910 1.925

0.071 0.074 0.032 0.134 0.037 0.074

super (R + β)

a For comparison, also the monomer excitation energies from uncoupled FDE calculations on the dimer (FDEu) are shown.

cm-1. This shows that the combined effect of kinetic-energy and XC contribution is almost negligible in the present approximation for the kernel. An additional effect that may become important at small distances are charge-transfer-like excitations between the subsystems. However, the approach in ref 40 does by construction not account for interactions with charge-transfer excitations. Although this issin the context of the problems for TDDFT to describe such excitations64-70s desirable for calculations on very large systems, which would suffer substantially from the long-range charge-transfer problem,35,70,71 it may lead to inaccuracies in the couplings at very short range. Because of this disagreement concerning the role of the shortrange couplings, we carried out an additional test of the FDE methodology. We performed a supermolecular TDDFT calculation in comparison to a coupled FDE calculation for a model of the (R,β)B850 dimer, in which the phytyl chain was omitted. The data are reported in Table 3. It can be seen that the splitting between the two Qy excitations is nicely reproduced by the coupled FDE calculation. The two states resulting from the coupling appear at 1.541 and 1.606 eV (FDEc) with oscillator strengths of 0.552 and 0.043, respectively, which compares well to the supermolecular result of 1.530 and 1.596 eV (oscillator strengths: 0.574 and 0.043), respectively. The splitting energy is 518 cm-1 in the case of FDEc and 536 cm-1 in the conventional supermolecular TDDFT calculation. The coupling constants calculated according to eq 6 are 242 cm-1 (FDEc) and 252 cm-1 (super), respectively (in the latter case, ω+,- were taken from the supermolecular calculation and ωa,b were taken from the uncoupled FDE calculations), which demonstrates that the additional approximations introduced in the subsystemTDDFT approach have only a minor effect on the coupling constants. It should be mentioned that the situation for the Qx transitions is more complicated. While the excitation energies for the lower of the two resulting states are similar (FDEc, 1.907 eV; super, 1.910 eV), there is a larger discrepancy for the higher-energy state (FDEc, 1.943 eV; super, 1.925 eV). The splitting between the two states in the supermolecular case is even smaller than the splitting between the assumed “site energies”, which are approximated by the FDEu calculations. Therefore, this splitting cannot be analyzed in terms of eq 6, which cannot describe a decrease of the energy gap between the two excitations. Possible reasons for this problem could be (i) inaccuracies in the uncoupled excitation energies (which are not a measurable quantity anyway, since one cannot obtain uncoupled excitation energies from an interacting dimer), or (ii) an interaction with one of the many higher-lying excited states of partial chargetransfer character that can be found in the supermolecular calculation but that are absent in the FDEc calculation. Nevertheless, both supermolecular and FDEc calculations agree

TABLE 4: Calculated (SAOP/TZP; in Units of cm-1) Coupling Constants for the Qy Transition of Bchl a in the Combined B800 and B850 Subunits (FDEc, Coupled FDE Calculation)a type

FDEc, B800 + B850

TDC14

B800B T B800C B800A T B800C B800A T RB850B B800A T βB850B B800A T βB850A

-31 -3 31 3 -3

-27 -3 27 23 5

a For comparison, also the results from ref 14 obtained from a Coulomb-only coupling with the transition density cube (TDC) method on the basis of scaled CIS transition densities are shown.

on the qualitative picture that the interactions between the Qx transitions are much smaller than those for the Qy transitions. In addition, it has been argued that it is not necessary to consider Qx transitions as donor transitions in EET processes because of a rapid Qx-Qy internal conversion.14 After this discussion of the rather strong interactions within the B850 unit, we will now address the couplings among the B800 Bchl a molecules and between the B800 and B850 Bchl a pigments. These couplings are weaker because of the larger distances between the chromophores. Therefore, they should be almost exclusively due to Coulomb coupling. Again, we compare our results to the (Coulomb-only) transition density cube (TDC) data based on the scaled CIS transition densities from ref 14 for the most important couplings. The data are shown in Table 4. It can be seen that there is in most of the cases a reasonable agreement in the computed coupling constants with those from ref 14. However, two differences can be seen: The B800A T βB850B coupling is much weaker in our calculation than in the reference, and the sign differs for the coupling B800A T βB850A. In order to analyze this discrepancy, we calculated the coupling constants in a dipole-dipole approximation according to14

Vd-d )

1 |µD|µA| κ 4π0 R 3

(9)

DA

where µD and µA are the transition dipole moments of the donor and acceptor transition, respectively, and κ is the so-called orientation factor

r A - 3(b r D‚R B)(b r A‚R B) κ)b r D‚ b

(10)

B R is a unit vector from the center of the donor to the center of the acceptor transition, RDA is the distance between these two centers, and b rD,A are unit vectors in the direction of the transition dipole moments. In Table 5, we report the distances RDA and orientation factors κ for the pairs of transitions from Table 4. For comparison, we also give the corresponding values from ref 14. Differences between our data and those from that reference can be due to the different crystal structure used in that study, which leads to slightly different bond distances, and differences in the orientation and magnitude of the transition moments. However, in general, there is a good agreement between the distances, orientation factors, and dipole-dipole couplings obtained here and in ref 14, e.g., for the B800B T B800C, B800A T B800C, and B800A T RB850B couplings. Nevertheless, the comparison suggests that the coupling constant reported for B800A T βB850B may be incorrect in that reference when compared to the B800A T RB850B coupling: Although the orientation factor is smaller by about a factor of

Properties of Natural Light-Harvesting Complexes

J. Phys. Chem. B, Vol. 112, No. 7, 2008 2215

TABLE 5: Calculated (SAOP/TZP; in Units of cm-1) Coupling Constants Vdd for the Qy Transition of Bchl a in the Combined B800 and B850 Subunits within the Dipole-Dipole Approximation, eq 9a this work

ref 14

coupling

RDA

κ

Vdd

B800B T B800C B800A T B800C B800A T RB850B B800A T βB850B B800A T βB850A

21.2 39.9 17.8 21.7 18.5

-1.20 -0.91 0.79 0.22 -0.10

-29 -3 33 5 -4

RDA

κ

Vdd

21.2 39.9 17.6 21.8 18.3

-1.33 -1.04 0.79 0.17 0.13

-26 -3 27 31 4

a

Additionally given are the distances between the pigment centers RDA (Å) as well as the orientation factors κ. For comparison, also the results from ref 14 obtained within the dipole-dipole approximation are shown.

5 in the first case and the distance between donor and acceptor is slightly larger (both here and in ref 14), the coupling constant calculated from the dipole-dipole approximation in ref 14 is even larger than that for B800A T RB850B. Since the magnitudes of all transition dipole moments are almost identical, one would expect that the B800A T βB850B coupling should be much smaller than the B800A T RB850B coupling. This is confirmed by our dipole-dipole results, which yield a coupling of 33 cm-1 for the latter and of 5 cm-1 for the former pair of transitions. As far as the discrepancy for the B800A T βB850A pair of transitions is concerned, we would like to note that our dipole-dipole approximation yields an orientation factor of approximately the same magnitude but opposite sign when compared to the results in ref 14. The sign of the coupling constants is actually arbitrary, since already the sign of the transition densities (or transition moments) is, in principle, arbitrary for each monomer, but it must be used consistently in the calculation of all coupling constants. In our case, we have chosen the same sign of the transition moment for each chromophore within the B850 ring. For this choice, the orientation factors for the B800A T βB850B and B800A T βB850A couplings, and consequently also the coupling constants within the dipole-dipole approximation, differ in sign. 7. Absorption Spectra of the B850 and B800 Units of LH2 The absorption spectra of the B850 unit have been obtained from the coupled FDE calculation described in section 6. In Figure 5, we compare the monomer spectra resulting from the Qy transitions of isolated RB850 and βB850 pigments with those obtained for the (R,β)B850 dimer, and the full B850 unit (i.e., 18 Bchl a molecules). The monomer peaks appear at 779 and 788 nm, respectively, for isolated RB850 and βB850. The coupling in the dimer leads to a considerable splitting with an intense peak at 800 nm and a less intense peak at 769 nm. In the full B850 ring of Bchl a monomers, the intense peak (scaled with a factor of 0.1 in Figure 5) occurs at 812 nm, and an additional peak of low intensity can be found at 759 nm. The main reason for the discrepancy of the Qy excitation energies with the absorption maxima in the natural light-harvesting system (about 850 nm) is small errors in the site energies. When using the empirical site energies from ref 16 of 1.56 eV (RB850), 1.50 eV (βB850), which differ by only 0.04 and 0.08 eV from the energies of the monomers in the uncoupled FDE calculations (see Table 1), the intense peak shifts to 852 nm, while the less intense peak shifts to 785 nm. This is in nice agreement with the spectra modeled in ref 16, where these two peaks were found at 857 and 773 nm, respectively.

Figure 5. Qy contributions to the absorption spectra (SAOP/TZP) of (subunits of) the B850 ring of LH2 (Rhodopseudomonas acidophila). Only the Qy transitions are considered. A Gaussian broadening of 0.01 eV has been applied to the peaks in the calculated absorption spectrum. Note that the spectrum of the full B850 unit was scaled by a factor of 0.1 for better comparability.

Figure 6. Absorption spectra (SAOP/TZP) of the B850 ring and the combined B850 and B800 rings of LH2 (Rhodopseudomonas acidophila). A Gaussian broadening of 0.01 eV has been applied to the peaks in the calculated absorption spectrum.

Figure 6 contains the spectra calculated from both the Qy and Qx contributions of the B850 subunit and the combined B850 and B800 subunits of LH2. It should be noted that these spectra do not include ensemble average effects like static disorder or fluctuations in the environment and just serve the purpose to illustrate the results of the subsystem-TDDFT calculations. The additional signals caused by the B800 ring in the lower panel of Figure 6 can clearly be recognized, and no significant changes occur in the bands already present in the B850 spectrum. 8. Conclusion In this study, it was analyzed how specific interactions with neighboring peptide side chains as well as interactions with other pigment molecules change the absorption properties of chromophores in LH2 of Rhodopseudomonas acidophila. In contrast to earlier studies, this was performed by means of a subsystem approach to DFT, which is capable of describing both types of interactions in a consistent DFT/DFT embedding framework. Both the shifts in the site energies due to hydrogen bonding to the Bchl a molecules and effects of peptide side chains acting as axial ligands to the magnesium atoms as well as splittings and shifts of the absorption bands due to EET couplings could be reproduced. By establishing a relation between coupling constants in the framework of TDDFT to CI-like coupling

2216 J. Phys. Chem. B, Vol. 112, No. 7, 2008 constants, good agreement with earlier calculations of EET coupling constants for LH2 could be observed. Although the results are in accordance with those from ref 14, a disagreement remains for one particular coupling constant both for the results obtained on the basis of transition densities and within the simple dipole-dipole approximation. This application demonstrates the efficiency of the general subsystem-TDDFT approach proposed in ref 40, that is based on FDE and its TDDFT extension.33,34 The limiting factor for the applicability of the method is the size of the largest subsystem, but even individual pigments of the size of a Bchl a molecule (140 atoms), combined in an aggregate with several thousands of atoms, can be studied with reasonable effort. The method shares the advantages of other DFT-based methods as well as the transition density cube (TDC) method from ref 14 in terms of the Coulomb couplings; i.e., it can be applied for closely interacting pigments and weakly allowed transitions. This is due to the fact that the Coulombic interaction between local excitations is essentially treated exactly, apart from the fact that density fitting techniques are applied for the transition densities, which is done in many conventional TDDFT implementations as well. Problems related to a discrete representation of the Coulomb potential induced by the transition density are thus avoided. The present approach may be even more efficient than the TDC method, since the latter uses a double numerical integration of the Coulomb interaction. Because of the more reliable oscillator strengths in TDDFT calculations compared to CIS, a good agreement for the Coulombic couplings between our calculation and the empirically scaled results in ref 14 is obtained. Similar to the DFTbased approaches in refs 30 and 31, also non-Coulombic effects are captured (see section 3 for a discussion of similarities and differences). In addition, this approach is properly integrated into the corresponding ground-state subsystem-DFT approach, so that also ground-state interactions between the monomers are correctly included. It turned out that short-range effects in the present approach, in which they are only included via an ALDA kernel, are very small. Although there is a general agreement on the fact that Coulombic effects are much more important for the strongly allowed transitions of the pigments in LH2, the magnitude of short-range effects found here is much smaller than in a previous study.15 The agreement of the couplings calculated for the model dimer in a supermolecular and subsystem fashion shows that this disagreement is not primarily due to the subsystem approach, and also other DFTbased investigations found that such short-range contributions are typically very small.30 Nevertheless, the results obtained here will at least partly be affected by the missing explicit couplings to charge-transfer excitations and deficiencies in the ALDA kernel (in particular also the kinetic-energy component). As far as environmental effects are concerned, the present study investigated direct effects of neighboring molecules, but no screening effects of a surrounding medium on the coupling constants were considered. Studies concerning this point in the context of light-harvesting phenomena were recently presented by Scholes et al.72 and Curutchet et al.48 They showed that such effects can be very important for excitation energy couplings in condensed matter. As FDE is developed as an explicit model for environmental effects, we will in future work explore the possibilities to include such screening effects in the frozendensity embedding framework. In conclusion, it was demonstrated that general systemenvironment interactions on excited states can be described by

Neugebauer this subsystem DFT approach, so that the range of application could be extended to complex aggregates of pigment molecules. Acknowledgment. The author would like to thank Prof. M. Reiher for helpful discussions and generous support. Funding by a Liebig-Stipendium of the Fonds der Chemischen Industrie (FCI) is gratefully acknowledged. References and Notes (1) Collings, A. F., Critchley, C., Eds. Artificial Photosynthesis, 1st ed.; Wiley-VCH: Weinheim, Germany, 2005. (2) Blankenship, R. E. Molecular Mechanisms of Photosynthesis; Blackwell Science: Oxford, U.K., 2002. (3) van Grondelle, R.; Dekker, J. P.; Gillbro, T.; Sundstrom, V. Biochim. Biophys. Acta 1994, 1187, 1-65. (4) Cogdell, R. J.; Gall, A.; Ko¨hler, J. Q. ReV. Biophys. 2006, 39, 227324. (5) Hu, X.; Ritz, T.; Damjanovic´, A.; Autenrieth, F. Q. ReV. Biophys. 2002, 35, 1-62. (6) Koepke, J.; Hu, X.; Muenke, C.; Schulten, K.; Michel, H. Structure 1996, 4, 581-597. (7) Cherezov, V.; Clogston, J.; Papiz, M. Z.; Caffrey, M. J. Mol. Biol. 2006, 357, 1605-1618. (8) Papiz, M. Z.; Prince, S. M.; Howard, T.; Cogdell, R. J.; Isaacs, N. W. J. Mol. Biol. 2003, 326, 1523-1538. (9) McDermott, G.; Prince, S. M.; Freer, A. A.; HawthornthwaiteLawless, A. M.; Papiz, M. Z.; Cogdell, R. J.; Isaacs, N. W. Nature 1995, 374, 517-521. (10) van Grondelle, R.; Novoderezhkin, V. I. Phys. Chem. Chem. Phys. 2006, 8, 793-807. (11) Fleming, G. R.; Scholes, G. D. Nature 2004, 431, 256-257. (12) He, Z.; Sundstro¨m, V.; Pullerits, T. J. Phys. Chem. B 2002, 106, 11606-11612. (13) Subramanian, V.; Evans, D. G. J. Phys. Chem. B 2004, 108, 10851095. (14) Krueger, B. P.; Scholes, G. D.; Fleming, G. R. J. Phys. Chem. B 1998, 102, 5378-5386. (15) Scholes, G. D.; Gould, I. R.; Cogdell, R. J.; Fleming, G. R. J. Phys. Chem. B 1999, 103, 2543-2553. (16) Scholes, G. D.; Fleming, G. R. J. Phys. Chem. B 2000, 104, 18541868. (17) Cogdell, R. J.; Howard, T. D.; Isaacs, N. W.; McLuskey, K.; Gardiner, A. T. Photosynth. Res. 2002, 74, 135-141. (18) van Oijen, A. M.; Ketelaars, M.; Ko¨hler, J.; Aartsma, T. J.; Schmidt, J. J. Phys. Chem. B 1998, 102, 9363-9366. (19) Freer, A.; Prince, S.; Sauer, K.; Papiz, M.; HawthornthwaiteLawless, A.; McDermott, G.; Cogdell, R.; Isaacs, N. W. Structure 1996, 15, 449-462. (20) Scholes, G. D.; Ghiggino, K. P. J. Chem. Phys. 1994, 101, 12511261. (21) Harcourt, R. D.; Scholes, G. D.; Ghiggino, K. P. J. Chem. Phys. 1994, 101, 10521-10525. (22) Scholes, G. D.; Harcourt, R. D.; Ghiggino, K. P. J. Chem. Phys. 1995, 102, 9574-9581. (23) Scholes, G. D.; Ghiggino, K. P. J. Phys. Chem. 1994, 98, 45804590. (24) Fo¨rster, T. Ann. Phys. 1948, 2, 55. (25) Chang, J. C. J. Chem. Phys. 1977, 67, 3901. (26) Madjet, M. E.; Abdurahman, A.; Renger, T. J. Phys. Chem. B 2006, 110, 17268-17281. (27) Howard, I.; Zutterman, F.; Deroover, G.; Lamoen, D.; Van Alsenoy, C. J. Phys. Chem. B 2004, 108, 19155-19162. (28) Russo, V.; Curutchet, C.; Mennucci, B. J. Phys. Chem. B 2007, 111, 853-863. (29) Curutchet, C.; Mennucci, B. J. Am. Chem. Soc. 2005, 127, 1673316744. (30) Iozzi, M. F.; Mennucci, B.; Tomasi, J.; Cammi, R. J. Chem. Phys. 2004, 120, 7029-7040. (31) Hsu, C.-P.; Hirata, S.; Head-Gordon, M.; Head-Gordon, T. J. Phys. Chem. A 2001, 105, 451-458. (32) Prince, S. M.; Papiz, M. Z.; Freer, A. A.; McDermott, G.; Hawthornthwaite-Lawless, A. M.; Cogdell, R. J.; Isaacs, N. W. J. Mol. Biol. 1997, 268, 412-423. (33) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050. (34) Casida, M. E.; Wesolowski, T. A. Int. J. Quantum Chem. 2004, 96, 577-588. (35) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesolowski, T. A. J. Chem. Phys. 2005, 122, 094115. (36) Neugebauer, J.; Jacob, C. R.; Wesolowski, T. A.; Baerends, E. J. J. Phys. Chem. A 2005, 109, 7805-7814.

Properties of Natural Light-Harvesting Complexes (37) Jacob, C. R.; Neugebauer, J.; Jensen, L.; Visscher, L. Phys. Chem. Chem. Phys. 2006, 8, 2349-2359. (38) Wesolowski, T. A. J. Am. Chem. Soc 2004, 126, 11444-11445. (39) Neugebauer, J.; Baerends, E. J. J. Phys. Chem. A 2006, 110, 87868796. (40) Neugebauer, J. J. Chem. Phys. 2007, 126, 134116. (41) Herek, J. L.; Wohlleben, W.; Cogdell, R. J.; Zeidler, D.; Motzkus, M. Nature 2002, 417, 533-535. (42) Casida, M. E. Time-Dependent Density Functional Response Theory for Molecules. In Recent AdVances in Density Functional Methods Part I; Chong, D. P., Ed. World Scientific: Singapore, 1995. (43) van Gisbergen, S. J. A.; Snijders, J. G.; Baerends, E. J. Comput. Phys. Commun. 1999, 118, 119-138. (44) Rosa, A.; Ricciardi, G.; Gritsenko, O. V.; Baerends, E. J. Excitation energies of metal complexes with time-dependent density functional theory. Structure and Bonding; Springer: Berlin, 2004; Vol. 112. (45) Wesolowski, T. A. One-electron Equations for Embedded Electron Density: Challenge for Theory and Practical Payoffs in Multi-Level Modeling of Complex Polyatomic Systems. In Computational Chemistry: ReViews of Current Trends; Leszczynski, J., Ed.; World Scientific: Singapore, 2006; Vol. 10. (46) Iannuzzi, M.; Kirchner, B.; Hutter, J. Chem. Phys. Lett. 2006, 421, 16-20. (47) Jacob, C. R.; Neugebauer, J.; Visscher, L. J. Comput. Chem., in press. (48) Curutchet, C.; Scholes, G. D.; Mennucci, B.; Cammi, R. J. Phys. Chem. B 2007, 111, 13253-13265. (49) Amsterdam Density Functional program; Theoretical Chemistry, Vrije Universiteit: Amsterdam, The Netherlands (http://www.scm.com). (50) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. J. Comput. Chem. 2001, 22, 931-967. (51) Alden, R. G.; Johnson, E.; Nagarajan, V.; Parson, W. W.; Law, C. J.; Cogdell, R. G. J. Phys. Chem. B 1997, 101, 4667-4680. (52) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100.

J. Phys. Chem. B, Vol. 112, No. 7, 2008 2217 (53) Perdew, J. P. Phys. ReV. B 1986, 33, 8822-8824. (54) Wesolowski, T. A.; Weber, J. Chem. Phys. Lett. 1996, 248, 7176. (55) Schipper, P. R. T.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. J. Chem. Phys. 2000, 112, 1344-1352. (56) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Chem. Phys. Lett. 1999, 302, 199-207. (57) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Int. J. Quantum Chem. 2000, 76, 407-419. (58) Perdew, J. P. In Electronic Structure of Solids; Ziesche, P., Eschrig, H., Eds.; Akademie Verlag: Berlin, 1991. (59) Lembarki, A.; Chermette, H. Phys. ReV. A 1994, 50, 5328-5331. (60) Wesolowski, T. A. J. Chem. Phys. 1997, 106, 8516-8526. (61) Atkins, P. W.; Friedman, R. S. Molecular Quantum Mechanics, 3rd ed.; Oxford University Press: Oxford, U.K., 1997. (62) Schaftenaar, G.; Noordik, J. H. J. Comput.-Aided Mol. Des. 2000, 14, 123-134. (63) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14.1, 33-38. (64) Casida, M. E.; Gutierrez, F.; Guan, J.; Gadea, F.-X.; Salahub, D.; Daudey, J.-P. J. Chem. Phys. 2000, 113, 7062-7071. (65) Tozer, D. J. Chem. Phys. 2003, 119, 12697-12699. (66) Dreuw, A.; Weisman, J. L.; Head-Gordon, M. J. Chem. Phys. 2003, 119, 2943-2946. (67) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. J. Chem. Phys. 2004, 120, 8425-8433. (68) Gritsenko, O.; Baerends, E. J. J. Chem. Phys. 2004, 121, 655660. (69) Maitra, N. T. J. Chem. Phys. 2005, 122, 234104. (70) Neugebauer, J.; Gritsenko, O.; Baerends, E. J. J. Chem. Phys. 2006, 124, 214102. (71) Lange, A.; Herbert, J. M. J. Chem. Theory Comput. 2007, 3, 16801690. (72) Scholes, G. D.; Curutchet, C.; Mennucci, B.; Cammi, R.; Tomasi, J. J. Phys. Chem. B 2007, 111, 6978-6982.