The atomistic modeling of light-harvesting complexes from the physical models to the computational protocol

In photosynthesis, the initial step is accomplished by specialized supercomplexes, photosystems, which comprise two components, light harvesting (LH) pigment–protein complexes, and the reaction center. LH complexes are unique biological systems, which integrate a protein matrix with a generally large aggregate of pigments. The number and type of pigments and the structure of the protein matrix are specific for each different photosynthetic organism.¹ Through the multichromophoric aggregate, the LH complex absorbs sunlight and rapidly funnels the energy toward neighboring complexes and ultimately to the reaction center.^2–4 The energy delivery by LH complexes is explained in terms of multichromophoric excited electronic states called excitons. Each exciton generally involves more than one pigment and is characterized by an energy slightly different from the others even if the aggregate is composed of the same types of molecules. The resulting energy ladder is not only determined by the strength of the interactions between excitations localized on the single pigments but also significantly tuned by the interactions of each pigment with the surrounding residues and the electric fields created by the whole protein and the external environment. Moreover, the excitons are coupled with the nuclear degrees of freedom of the multichromophoric aggregate and the protein.

Due to the complexity outlined above, an accurate atomistic modeling of LH complexes is extremely challenging, as it requires the integration of very different physical models, which span from quantum mechanical (QM) methods to classical models, dynamic descriptions, and statistical analyses. Moreover, for that integration to work, the different models need to be all coupled together, as the functioning of the LH complex is a real multiscale process in space and time.

In the last few decades, there have been important advances in the atomistic modeling of LH complexes,^5–10 which have been largely pushed by the increasing amount of experimental data about their structures, the spectroscopic properties, and the energy transfer processes.^11–17 In parallel, there has been an important activity in the development of computational methods that have allowed for simulating properties and processes of LH complexes embedded in different environments.^18–21 

This Perspective reviews these advances using the typical workflow followed in the computational investigation of a LH complex, which starts from a predicted or experimentally characterized structure, adds the environment, simulates the coupled dynamics of all the parts, calculates the excitonic states, and finally reproduces the optical spectra to be directly compared with experiments. A graphical representation of the workflow with the indication of all the different components is given in Fig. 1. In the presentation and discussion of the different parts of the workflow, we use a minor LH complex present in PSII of plants (CP29) as an example of application. Finally, the optical spectra of the complex will be compared with a recently proposed mutant²² to show the accuracy and completeness that the atomistic modeling can achieve in predicting even small effects as those induced by a single residue mutation.

At the present time, LH complexes of different photosynthetic organisms have been structurally characterized through x-ray diffraction or cryo-EM techniques. These structures represent a fundamental ingredient for the atomistic models; however, they also present important limitations. First of all, they do not include any information on the external environment, which for most of the LH complexes is a membrane. Moreover, fluctuations due to the temperature are not accounted for. Finally, the standard resolution of crystal or cryo-EM structures of LH complexes is not high enough to describe the internal geometry of the pigments with the atomic accuracy needed for safely applying QM calculations on them.

For all these reasons, an atomistic modeling of LH complexes based on only the crystal structure can give a qualitative description of the system at very low temperatures, but, in general, it cannot be representative of the room-temperature situation. To go beyond this too simplistic static picture, we can introduce molecular dynamics (MD) simulations. Molecular dynamics has been shown to be a very effective approach for sampling configurations of the LH complex and its environment at the temperature of interest.²¹ 

Molecular dynamics simulations of light-harvesting complexes usually started from the available x ray or cryo-EM structures deposited in public repositories.²³ Another means of obtaining starting structures for MD simulations when the experimental structures are not available is to employ homology modeling techniques, which have also been employed for LH complexes.^24–26 These techniques can be employed if the target protein sequence is sufficiently similar to the sequences of known structures. In addition, recent algorithms based on neural networks, such as AlphaFold 2,^27,28 have proven able to determine protein structures as accurately as their experimental counterparts, providing another means of obtaining starting structures for molecular simulations. Especially when the starting structure is chosen based on a computational guess, such as in homology modeling or by employing neural networks, MD simulations and their spectroscopic analysis can be employed to validate the structures and gain insights on the architecture of systems that are still to be structurally resolved experimentally. Finally, the starting point being either an experimental structure or a computational guess, it is easy to slightly alter the initial structure so as to introduce, e.g., point mutations. Especially in the context of investigating LH complexes, the simulation of mutants can be significantly valuable as it complements site-directed mutagenesis experiments. We provide such an example in Sec. V.

Due to the complexity and the large dimension of the whole system (comprising the complex and the surrounding environment), a simplified description is commonly used for the potential energy (and the corresponding forces) in terms of Molecular Mechanics (MM) Force Fields (FFs). At the present times, accurate FFs are available to describe protein structures and their thermal fluctuations and for describing water and lipid membranes. However, the application of the available FFs to the pigments is by far less accurate. This is the reason why a reparameterization of the FF for describing the nuclear degrees of freedom of the pigments is generally introduced. In this reparameterization of the FF, QM calculations of internal geometry and vibrational properties are commonly used as reference data.^29–31 

Once the potential energy model of all the components has been defined, the system can be assembled and its trajectory can be determined following the same protocols used in MD simulations of biomacromolecules. In order to achieve a meaningful sampling of the configurational space, ns timescales are often not sufficient. We have observed fluctuations on the 100 ns timescale in LH systems, and we, therefore, recommend μs timescales whenever possible.³² We also recommend running several replicas of the same system obtained using different initial velocities. Different replicas are then subject to standard convergence assessment tests, such as monitoring the kinetic energy of the system, the root-mean-square deviation (RMSD) of different parts of the protein, and the average area available to each lipid in the plane of the membrane. Other quantities, such as selected coordinates or order parameters, may be used to assess equilibration. Averages of these quantities may systematically drift during a transient time, indicating that the system has yet to sample the equilibrium ensemble. When the averages reach a plateau or oscillate around a mean value, the trajectory can be considered locally equilibrated. The transient part should be discarded before further analysis and especially before computing excitation properties.

The distinct feature of LH complexes is the dense packing of pigments in highly specific arrangements that maximize energy transfer while avoiding concentration quenching.³ At variance with other biomacromolecular systems, then, the analysis of LH complex trajectories should pay particular attention to the relative arrangement of the embedded pigments, i.e., their relative distances and orientations. Pigment–pigment orientation measures that are directly linked to their interaction strength can be based on the Förster orientational factor κ [see Eq. (5) in Sec. III B]. In addition, analyzing the hydrogen bond interactions between the protein and the bound pigments can provide useful insights on interpreting the spectroscopic properties of the complex.²⁵ This kind of analysis provides information on the interaction network within the LH complex and on how it is influenced by the protein dynamics.

Further effective analysis instruments aiming at disentangling the complex conformational landscapes in LH complexes can be borrowed from the ever-growing field of unsupervised learning.³³ In particular, dimensionality reduction methods [such as principal component analysis (PCA)³⁴ and Isomap³⁵] can be employed to consistently reduce the redundancy in the underlying conformational ensemble starting from properly featured protein conformations,³⁶ thereby improving the ease of characterization of protein conformations. In a similar spirit, kinetic-based independent component analysis algorithms, such as time-lagged independent component analysis (tICA),³⁷ can also be used to characterize LH complexes conformations while conserving some degree of information lying in the temporal structure of the trajectory.³⁸ Additionally, clustering techniques can be adopted to group similar conformations together, effectively coarse-graining the explored structures and facilitating the identification of key structural differences between the observed conformers. When running subsequent QM/MM calculations, an a priori division of the explored ensemble into clusters may be useful to reduce the number of calculations while ensuring an accurate sampling of the LH complex conformations. In addition, clustering can be used to connect structural features to spectroscopic observables.^26,36

The electronic properties of the multichromophoric aggregate of LH complexes are commonly described through the Frenkel exciton model. In this model, the electronic excited states of the entire aggregate are expanded on a basis of excitations localized on single pigments. In this basis, the excited states of the aggregate are described by the following exciton Hamiltonian:

where $Ei$ are denoted as site energies and represent the excitation energies of the isolated pigments and V_ij are the excitoniccouplings between different pigments. Here, we assume for simplicity that each pigment contributes with one excitation. The wavefunction of the exciton state K can be expressed in terms of the exciton coefficients $ciK$⁠,

From a strictly computational point of view, the exciton model can be seen as a “divide and conquer” method, where the complexity of the multichromophoric aggregate is reduced by separating the problem into smaller elements. In fact, the matrix elements of the Hamiltonian (1) can be easily computed from quantum chemical calculations on single pigments. However, the exciton model has enormous advantages that go beyond the mere computational efficiency. Strikingly, the exciton coefficients allow for a straightforward interpretation of the aggregate’s excited states and of their energy variations. As we will see in Sec. IV, this model also enables the description and calculation of optical line shapes.

Physically, the exciton model corresponds to neglecting all charge-transfer contributions to the aggregate excited states. This is generally a reasonable assumption as the inter-pigment distances in most of the LH complexes are small enough to have non-negligible excitonic couplings, but at the same time, they are large enough to avoid significant charge transfer. However, when needed, the Hamiltonian (1) can be augmented with charge-transfer states by including their energies and couplings with localized excitations.^39,40

Site energies represent the first main ingredient in the exciton Hamiltonian. Each site energy can be computed as the excitation energy of a pigment embedded in its environment. As such, calculation of these quantities requires an appropriate QM method and a description for the environment. Among the QM methods, a fair compromise between computational cost and accuracy is represented by linear-response time-dependent density functional theory (TD-DFT). Despite its well-known deficiencies, TD-DFT represents a valid and extensively tested method to compute excited states. More accurate QM methods have been employed for light-harvesting complexes, but their use has been mostly limited to benchmark calculations. Multireference methods, such as multireference configuration interaction-DFT (MRCI-DFT) and complete-active-space SCF (CASSCF), possibly in combination with perturbation theory (CASPT2), have been used to complement or assess TD-DFT methods.^41–44 More recently, single-reference accurate methods have been employed. Approximate coupled cluster methods, such as second-order coupled cluster (CC2) and pair natural orbital coupled cluster [domain-based local pair natural orbital similarity transformed equation of motion-coupled cluster singles and doubles (DLPNO-STEOM-CCSD)], have been applied to photosynthetic pigments.^45–47 Other alternatives are post-SCF methods based on equations of motion or response approaches, such as the algebraic-diagrammatic construction (ADC) and the Bethe–Salpeter equation in the GW approximation (GW/BSE). These methods have found applications for photosynthetic pigments.^45,48 Semiempirical methods represent a much less expensive alternative to obtain excitation properties with reasonable accuracy. Among these, time-dependent density functional tight binding (TD-DFTB) and its long-range corrected (LC) version showed promise for describing light-harvesting pigments.⁴⁹ 

Pigments in LH complexes are surrounded by a highly anisotropic environment, and they are bound in a protein pocket by establishing specific interactions, such as hydrogen bonds. Only a computational model that describes this anisotropy can reproduce excitation energy variations among the different sites of the complex. The simplest such model is quantum mechanics/molecular mechanics (QM/MM) with electrostatic embedding. In electrostatic embedding QM/MM, the chromophore is treated with the QM level of choice, whereas the environment is described as a set of fixed point charges. Irrespective of the QM description, this model introduces an external one-electron operator analogous to the electron–nucleus Coulomb interaction. The conceptual simplicity and ease of implementation coupled with any QM method have made this model extremely popular for excited-state calculations, also in light-harvesting complexes.^47,50,51 While electrostatic embedding QM/MM can describe the effect of the MM environment on the QM part, it cannot describe how the environment is affected by the presence of the QM electron density or by its changes in response to excitation.

A more refined QM/MM flavor is the polarizable embedding QM/MM.⁵² Polarization has been included in QM/MM using various formulations based, e.g., on Drude oscillators,^53–55 fluctuating charges,^56–58 or induced point dipoles (IPDs).^59–66 

In any of these formulations, the environment and the QM density can mutually polarize each other. From the point of view of the QM method, this formulation introduces a nonlinear term in the electronic Hamiltonian through effective two-electron contributions. These terms render much less straightforward coupling a polarizable embedding to a new QM method. Moreover, the response of the MM part also becomes more complicated. Calculating the MM polarization requires to account for the mutual polarization between MM atoms. Within the IPD formulation, the induced dipoles on the MM part can be computed by solving the following system:^66,67

where E^ext is the electric field resulting from the QM density and the MM charges and A is a matrix containing the atomic polarizabilities and the interactions between all induced dipoles. To avoid the so-called “polarization catastrophe,” e.g., the divergence of the Coulomb interaction between two point dipoles when they get too close, distance dependent damping functions originally proposed by Thole⁶⁸ are commonly introduced.

All these specificities make polarizable QM/MM methods computationally more expensive than standard electrostatic embedding analogs. New linear-scaling implementations are required to reduce the computational overhead of the formulation.⁶⁹ The IPD formulation of polarizable QM/MM methods (from now on QM/MMPol) has been extensively used for calculating site energies of pigments within LH complexes coupled with linear-response TD-DFT methods.^{19,25,70–72}

It is important to disentangle the different effects that the environment has on the site energies. The excitation energy of a chromophore can be modulated by the environment directly, through electrochromic effects, and indirectly by changing the internal geometry of the chromophore. These effects can be separated naturally by computing the excitation energy of the chromophore embedded in the MM(Pol) environment and isolated in vacuo.¹⁹ The difference between the two calculations represents the direct electrochromic effect due to the environment. We stress that the geometry of the pigment must be the same in both calculations in order to obtain a pure electrochromic effect when taking the difference.

The indirect geometrical effect on the site energies can be defined by comparing the in vacuo site energies of each distinct chromophore with a reference geometry. The environment contribution is very sensitive to the fluctuations of the environment. As such, an adequate sampling of environment degrees of freedom via, e.g., MD simulations, is needed to obtain a realistic description of the chromophores embedded in LH proteins. Indeed, as shown in Fig. 2, an analysis based on the experimental structure only, ignoring thermal effects included with the MD, provides a misleading picture of the influence of the environment on the excitation energy of the embedded chromophores. This kind of analysis allows us to separate “geometrical” and “environment” effects on the site energies and to analyze the origin of the specific tuning of the site energy of each pigment.⁷² 

Despite their power, QM/MM calculations of excitation energies suffer from systematic errors stemming mainly from the quantum chemical method and, in part, on the method used for obtaining chromophore geometries. In fact, the expected accuracy of QM methods for excitation energies is at least around 0.2 eV for TD-DFT methods and ∼0.1 eV for more accurate methods, such as ADC(2) and CC2.⁷³ These errors are larger than the typical width of site energies in LH complexes, but they do not prevent a quantitative application of QM/MM methods to computing site energy ladders. In fact, while these errors are large, they are usually systematic, and they do not affect excitation energy differences among chemically identical pigments. For this reason, it has become customary to rigidly shift the calculated spectra, or equivalently the site energies, by a certain amount that takes into account the intrinsic error of the QM level of theory.^25,51,72,74 This correction is justified by assuming that the error is systematic. In this respect, the same shift should be applied to different systems if they contain the same pigments and are treated at the same level of theory.

The situation is more complicated when a LH complex contains different types of pigments, such as Chl a and Chl b, or when multiple excited states need to be considered. There, it is not guaranteed that the QM method will give the same error for all excited states. For example, even the DLPNO-CCSD method gives different errors for the different transitions of Chl a.⁴⁶ A possible solution is to use a different parametric shift for different pigment types. This approach seems quite straightforward, but it introduces a large number of empirical parameters. In some cases, the relative error for two different pigments is small enough not to influence substantially the spectral shape,⁷² and it can be simply ignored. When the relative error cannot be ignored, the safe choice is to determine the shift parameters from a system other than the system of interest, and possibly simpler. For example, the excitation energies of chromophores in solution could be used to determine the relative error that a given QM method makes on different transitions.

A much simpler approach to environment effects on site energies of pigments in LH complexes is the use of continuum models.^75,76 These models lose the fine details of the protein environment, which makes the description necessarily less accurate than an atomistic model. However, continuum models can still represent a useful and inexpensive alternative to QM/MM especially if the QM subsystem is augmented to include the protein residues involved in specific interactions with the pigments.^70,77 A recent alternative is to couple both MMPol and a continuum model to the QM description of a chromophore.^78,79 This strategy has been applied to phycocyanin complexes of red algae.⁸⁰ 

Other approaches are able to go beyond the classical description of the protein environment. The frozen density embedding (FDE) scheme and its extension to excited states (subsystem TD-DFT) divide the entire system into fragments, each described quantum mechanically and experiencing the embedding potential from the other QM fragments.⁸¹ This method has been applied to computing site energies and couplings in LH complexes.⁸² 

Electronic couplings between localized excitations determine exciton delocalization effects, which reflect on optical spectra of LH complexes and on the exciton dynamics.

The electronic coupling can be decomposed into a sum of a long-range Coulomb contribution and a short-range term.^83,84 The latter term is dependent on the overlap between orbitals of the interacting chromophores and is generally neglected. In the common case of bright excitations, the Coulomb term is dominant even when the two chromophores are in contact. Therefore, the coupling can be calculated as

where $ρi/jtr$ is the transition density of the excitation localized on chromophore i or j.

Equation (4) has a classical interpretation as the electrostatic interaction between two charge densities that represent the excitations of the coupled chromophores. For this reason, the Coulomb coupling is often approximated through a multipole expansion or a distributed monopole expansion. The first multipolar order is the interaction between transition dipoles $μ⃗tr$⁠,

which can be expressed solely in terms of the distance between chromophores R_ij and the orientational factor κ. While this approximation breaks down at the interchromophore distances relevant for LH complexes^84,85 and, thus, it should not be used for computing electronic couplings, it is still useful to uncover and interpret trends or differences in computed couplings or to quickly prescreen inter-pigment interactions.

Contrary to multipole expansions, distributed expansions place a charge on each atom of the chromophore in order to reproduce the transition density and calculate the Coulomb coupling as an interaction between collections of point charges. In the TrEsp method,⁸⁶ these atomic transition charges are determined by fitting the electrostatic potential generated by $ρi/jtr$⁠. This method has proven very effective in approximating the true Coulomb coupling almost exactly.

Importantly, the electronic couplings can also be affected by the surrounding environment. The environment influences the QM density of the chromophore, and therefore, it also affects the transition properties. This results in an indirect effect on the couplings. This effect can be described with any QM/MM embedding (electrostatic or polarizable). If the couplings are calculated in the context of subsystem TDDFT,^82,85 the effect of the surrounding environment is also included.

A second effect arises from the polarizability of the environment. The environment polarizes in response to the electric field generated by an electronic transition on chromophore j and, in turn, distorts the electric potential experienced by the other chromophore. This effect results in a completely new contribution to the Coulomb coupling. In the IPD formulation of polarizable QM/MM models, this contribution can be calculated as the interaction between the transition density $ρitr$ and all the dipoles $μkMMPol$ induced on the environment atoms by the other transition density,⁶² 

where $μkMMPolρjtr$ are obtained by solving Eq. (3) with an electric field generated by $ρjtr$⁠. Due to the symmetric bilinear nature of the interaction, this contribution does not change by exchanging i and j.⁶⁷ This new explicit contribution from the environment is proportional to the magnitude of both transition densities; therefore, it increases with an increase in the Coulomb coupling. In addition, it often has opposite sign to the Coulomb contribution, thus resulting in a reduction of the total coupling with respect to the Coulomb one.

If one approximates the environment as a dielectric, the sum of $VijMMPol$ and $VijCoul$ can be thought as an interaction screened by the dielectric. This model has been extensively used in the literature, for example, coupled to the TrEsp method in the so-called Poisson-TrEsp strategy⁸⁷ and to the Polarizable Continuum Model (PCM).⁸⁸ However, in the highly inhomogeneous environment of LH complexes, the mediated interaction is not always well represented by a dielectric screening. In fact, a different positioning of the environment molecules around the chromophores can lead to substantially different screenings on the Coulomb coupling.^9,40,89

As mentioned above, the Coulomb coupling is often a very good approximation of the total electronic coupling, at least for bright transitions. Things change when we consider optically forbidden transitions and chromophores in van der Waals contact with each other. For optically forbidden transitions, the first multipole term of the Coulomb coupling is zero. This does not guarantee that the Coulomb coupling is negligible, but it usually means that this term is quite small. When two chromophores are in contact, the short-range coupling term is nonzero. Combining the two cases, it is possible that the short-range coupling is of the same magnitude as the Coulomb coupling, if not larger.

Unfortunately, short-range contributions to the coupling are not as easy to compute as the Coulomb term.⁹⁰ Some short-range terms arise from the fact that, at close contact, the electron clouds of the two chromophores experience Pauli repulsion, i.e., distort to satisfy the exclusion principle. As a consequence, perturbation expressions based on the isolated chromophores do not yield accurate results for this contribution.⁸³ 

A more complete method to account for all contributions in electronic couplings is to use a QM calculation of the chromophoric dimer.^40,90,91 In this case, the electronic Hamiltonian has to be transformed to the monomeric “site” basis with a localization procedure. Given that the dimer of chromophores is considered in the QM calculation, this approach gives, in principle, all contributions to the coupling. However, this strategy is clearly more computationally expensive. On the other hand, this strategy allows for deriving both exciton and charge-transfer couplings from the same calculations^40,92 and, thus, describe all types of short-range interactions among chromophores.

The exciton model described in Sec. III is the basis for the simulation of the optical spectra of LH complexes. For example, the same excitonic coefficients defining the exciton states in Eq. (2) can be used to reconstruct the excitonic transition dipoles,

which, together with the excitonic energies, determine the intensity and position of the expected spectroscopic signals. For excitonic systems, however, vertical excitations and transition dipoles are not enough to properly reproduce the spectra. In fact, as for single molecules, vibronic and broadening effects have to be accounted for.

From now on, for the sake of simplicity, we focus on absorption spectra, but similar considerations can be given for other optical spectroscopies and, at least in part, for two-dimensional electronic spectroscopy.²⁰ 

Whereas the excitation energies and dipoles describe the position and intensity of the spectra, the band shape is determined by coupling of the electronic transitions to the vibrational degrees of freedom. The most common approach to include effects of the vibrational motion is the second order cumulant expansion in the exciton basis.⁹³ The cumulant expansion formalism is exact for Gaussian fluctuations of the excitation energy, i.e., for a linear coupling of the excitation energy to the vibrations of the system.⁹⁴ 

The most commonly used flavor of cumulant expansion yields the absorption spectrum as a sum over exciton states,

where each exciton K is characterized by its transition dipole μ_K, excitation energy ω_K, and lifetime τ_K, where the latter depends on the rates of exciton relaxation from state K to the other states. The information on vibronic coupling is contained in the line shape function g_K(t) for each state K. Similarly to the excitonic transition dipoles [Eq. (7)], g_K(t) is reconstructed from the quantities of the individual pigments and their exciton coefficients, namely,

where g_i(t) is the line shape function of the pigment i.

The line shape function can be obtained from the spectral density $Jω$ as

where β = 1/k_bT introduces the temperature dependence of the optical spectra. The intensity of the spectral density J(ω) describes how much each vibrational mode at frequency ω influences the excitation energy. Different approaches have been proposed in the literature for the calculation of spectral densities: the most used are described in Sec. IV A.

The formalism described above is only one of the possible choices when using the cumulant expansion in an exciton system. In fact, it contains two key approximations. First, the excited-state dynamics is assumed to be Markovian and exciton relaxation is described by an exponential decay function with lifetime τ_K, which is usually derived with the Redfield rate equation. Various approaches at overcoming this approximation have been devised, either by correcting the Redfield rates⁹⁵ or by directly considering time-dependent rates in the cumulant expansion.⁹⁶ The second approximation is the secular approximation, which neglects off-diagonal elements in the cumulant expansion and allows for writing the spectra as a sum over exciton states [Eq. (8)]. The most general approach, the full second-order cumulant expansion, relaxes both of these approximations.^97,98

So far, we have introduced only the effects due to the “dynamic” component of excitation energy fluctuations, which defines what is also called “dynamic disorder” and contributes to the homogeneous line shape. However, an additional effect, commonly denoted as “static disorder,” should be considered as well. Here, the motions responsible for the disorder are those that lead to energy fluctuations slower than the absorption process, such as conformational changes of the environment and/or the pigments. Static disorder is usually included by introducing a set of exciton Hamiltonians calculated using site energies and couplings obtained for different configurations of the system (e.g., uncorrelated snapshots extracted from the MD simulation).^72,99–101 The set of eigenstates resulting from this ensemble of excitonic Hamiltonians is then used to generate the corresponding set of spectra through Eq. (8). These spectra are finally convoluted together to obtain the final ensemble-averaged spectrum.

When applying such a strategy, however, one has to pay particular attention to how the different excitonic parameters respond to changes in the pigment geometry. The electronic couplings are quite stable with respect to changes in the internal geometry of the pigments depending more on their relative position and orientation, particularly when bright transitions are involved. Site energies, on the other hand, are known to vary significantly upon modest changes of the molecular geometry and even more when the pigment presents an extended conjugated chain as in the case of (bacterio)chlorophylls and carotenoids.²⁵ If a MM-MD trajectory is used for considering static disorder, site energy changes following alterations of the internal geometry of the pigments can be inflated if an inaccurate force field for the pigments is employed. Furthermore, even if an accurate force field is employed, the fluctuations of the site energies with the fast (e.g., molecular vibrations) degrees of freedom must be removed from the ensemble, to avoid double counting of this contribution when computing optical spectra. A viable route for introducing the static disorder in the simulation of the spectra is presented in Sec. IV B.

The spectral density describes how the nuclear degrees of freedom couple to a specific electronic transition. In the spirit of the second-order cumulant expansion, the spectral density is defined from the time-dependent fluctuations of the site excitation energy,

where $ΔE(t)=E(t)−⟨E⟩$ is the fluctuation of the site energy and $ΔE(t)ΔE(0)g$ denotes its autocorrelation function in the ground-state nuclear ensemble of the pigment.

The real nuclear quantum ensemble is generally inaccessible to practical simulations, but further simplifications can be made. In the case of linear coupling to a harmonic vibrational bath, the spectral density takes the simple form,

where ω_k are the frequencies of the normal modes and S_k are the corresponding Huang–Rhys factors, which quantify the dimensionless displacement of the excited-state potential energy surface with respect to that of the ground state. From the spectral density, one can obtain the reorganization energy as follows:

which measures the energy change of the pigment in its excited state when it moves from the ground-state equilibrium geometry to that of the excited state.

In general, the spectral density of an embedded pigment can be separated into two contributions, one due to the intra-pigment vibrational modes $Jvibω$ and one arising from the motions of the environment $Jenvω$⁠.

The environmental part of the spectral density is commonly represented as a continuous function with low frequency contributions. This function is obtained either from the fitting of the experimental spectra, i.e., the difference between absorption and fluorescence spectra, or from a normal mode analysis of the whole pigment protein complex.¹⁰² The intra-pigment part of the spectral density is composed of discrete peaks at higher frequencies representing normal modes localized on the pigment. This part gives rise to what is normally called vibronic coupling. High-resolution spectroscopy, such as low-temperature fluorescence line narrowing (FLN) spectra,^103–105 can resolve individual vibronic peaks and give information on their frequencies and corresponding Huang–Rhys factors. Therefore, the $Jvibω$ contribution has been often extracted from experiments.

Normal-mode frequencies can be also computed from standard normal-mode analysis, and Huang–Rhys factors can be computed in the vertical gradient (VG) approximation.^51,106–109 In this approximation, the Huang–Rhys factors are computed from the gradient of the excited PES at the ground-state equilibrium geometry, and the spectral density is constructed as

where q_ξ are pigment normal mode dimensionless coordinates. The damping factors γ_ξ represent the finite correlation time of the normal mode.

Alternatively, the spectral density can be obtained through MD trajectories. Along such trajectories, the excitation energy is then computed at every time step, obtaining the so-called classical autocorrelation function $Cclt$⁠,

where a quantum correction βω is employed to satisfy detailed balance,¹¹⁰ and the classical autocorrelation function is computed as

At variance with the vertical gradient, MD simulations do not rely on the harmonic approximation and, therefore, effectively incorporate some degree of anharmonicity.

In earlier applications, the trajectory of the system was obtained using MM-based MD. The SDs computed in this way, however, were too sensitive to the quality of the underlying force field, which has to match the QM description of the chromophore to give correct results.^31,112–116 More recently, the increase in computer power and algorithms has allowed for the use of Born–Oppenheimer QM/MM(Pol) MDs to generate the trajectory.^101,117,118 Commonly, a DFT description is used for the ground state QM/MM-MD simulations, which is then combined with a TD-DFT modeling of the excitation energies. Recently, semiempirical methods, such as DFTB, have been proposed as a computationally effective alternative to compute spectral densities of pigments in LH complexes.^50,74,119

To give an example of how different models for computing the spectral density influence the corresponding line shape, in Figs. 3(a) and 3(b), we compare the experimentally derived spectral density of CP29 with those obtained from vertical gradient calculations or QM/MM MD simulations. The VG calculations were performed for Chls a604 and a610 of CP29, while the QM/MM MD spectral density refers to Chl a604.⁷⁴ As one can note, comparing the two VG calculations, it is difficult to distinguish Chl a604 from a610. Moreover, both models agree quite well with the experimental estimate, with the QM/MM MD estimates showing a slightly better fit of the peak positions. This suggests that the inclusion of anharmonicity through QM/MM MD improves the description of intermolecular vibrations.

The different SD models have also been used to simulate the absorption line shape. The spectra are computed using the same adiabatic excitation energy $ωeg0=ωeg−λ$ [Fig. 3(c)]. From this comparison, it is clear that the differences between the experimentally derived spectral density (SD exp) and the computed ones have a limited impact on the actual line shape. In fact, despite the apparent differences in Figs. 3(a) and 3(b), the line shapes in Fig. 3(d) are all very similar and differ mainly in the relative intensity of the main 0–0 peak and the vibronic band.

Strikingly, among the vertical gradient calculations, the one on Chl a610 more closely resembles the line shape obtained from the experimental spectral density. This suggests that the low-temperature FLN experiment extracts the vibrational frequencies of Chl a610, which belongs to the terminal emitter state. The line shape calculated for Chl a604 deviates from the experimentally derived one and even more so does the one computed from QM/MM MD, which takes into account temperature effects. This comparison suggests that the experimentally derived spectral density is biased toward the low-energy Chls, and an accurate representation of vibronic coupling requires QM/MM calculations.

On the other hand, the line shapes obtained with the different methods are sufficiently similar to each other. For the purpose of computing exciton spectra, the choice of spectral density is probably not critical.

As reported above, MD is an efficient route to sample the configurational space of the solvated LH complex. The resulting trajectories can be finally used to extract different configurations of the system on which different excitonic Hamiltonians are calculated to account for static disorder in the simulation of the spectra.

We consider here the case of an LH complex containing several identical chlorophylls, and we assume to have performed different MD replicas of the LH complex embedded in its environment. The value of the site energy calculated at a given configuration of the whole system (LH complex + environment) extracted from an MD trajectory can be interpreted as a sum of different contributions, namely,

where i, r, and f label the pigment, the index of the MD replica, and the frame of the replica, respectively. Figure 4(a) shows the high variability of $EMD$⁠, stemming from the dynamics of both the pigment and environment degrees of freedom.

Here, $Egeom,fast$ indicates the sum of the site energy calculated on the isolated pigment in its equilibrium geometry and a correction stemming from the variation of its fast nuclear degrees of freedom, such as the bond lengths, at that specific configuration of the system. In parallel, $Egeom,slow$ is the further change of the site energy due to change of the slower nuclear degrees of freedom of the pigment, such as angle distortions due to protein environment constraints. The total geometric contribution $Egeom=Egeom,fast+Egeom,slow$ is directly accessible from a vacuum QM calculation of the pigment using a structure sampled from the MD simulation of the LH complex. Finally, $Eenv$ is the explicit environment contribution, obtained by subtracting the total geometric component from a QM/MM(pol) calculation as $Ei,r,fenv=Ei,r,fQM/MM(Pol)−Ei,r,fgeom$⁠.

A straightforward separation of the fast and slow components of the geometrical contribution is not available. Still, one can exploit the remarkable linear dependence of the Chl site energy from the bond lengths of the conjugated ring²⁵ to obtain an estimate of the fast component $Egeom,fast$⁠. After featurizing the chlorophyll geometry as a set of bond lengths, the estimate of the fast geometrical component of the site energy $Êi,r,fgeom,fast=xi,r,fw$ can be obtained by minimizing the squared loss within a regularized linear model,

where x_i,r,f denotes the set of chlorophyll bond lengths (padded with 1 to include the intercept in the fitting), w is the coefficient vector of the linear model, and λ is a regularization parameter that prevents overfitting. The determination of $Êgeom,fast$ through the regularized fitting is shown in Fig. 4(b).

The slow component of the geometrical contribution is then obtained by subtracting the fast component as $Êi,r,fgeom,slow=Ei,r,fgeom−Êi,r,fgeom,fast$⁠. Figure 4(c) shows the distribution of $Êgeom,fast$⁠, $Êgeom,slow$⁠, and $Eenv$⁠. The fast geometrical contribution evidently accounts for the most variance in $EMD$⁠, with the slow geometrical contribution and the environment contribution presenting a comparable variability.

Once all the contributions have been computed for each pigment, replica and configuration, we define the correct site energy to be used in the evaluation of static disorder as

where $⋅x$ denotes the average with respect to x.

The site energies computed with Eq. (19) are shown in Fig. 4(d). The distribution of each site energy is represented with a box plot, and we also report the confidence interval obtained via bootstrapping (red dashed lines). The statistical uncertainty on the mean site energy is around 50 cm⁻¹ for most Chls and ensures that the obtained site-energy ladder is statistically significant.

In Eq. (19), the fast geometrical contribution has been averaged on all the chlorophylls and all the configurations of the different replicas. For the slow contribution instead, the average is performed on all the configurations within each replica separately but not on the pigments, as we have assumed that different chlorophylls will have different distortions due to their specific environment constraints.

With this scheme, we (i) correct the intrinsic error of the MM-MD simulation for the pigment bond lengths, (ii) exclude double-counting of bond length contribution to the energy disorder, and (iii) reduce the “noise” in the excitation energies coming only from the limited sampling of pigment geometries. We also point out that other approaches avoiding double-counting of vibrational motions in the final spectra, relying on geometry optimizations of the pigments, are present in the literature.^120,121

Another approach to separating static and dynamic disorder in MD simulations is based on the calculation of site energies for a large number of closely spaced frames.¹²² If the site energy is mediated over a time window T, the fast fluctuations are reduced by a factor $T$⁠, while slower fluctuations are unaffected. The standard deviation of the means can be extrapolated for large T, yielding the magnitude of static disorder. This approach requires many calculations along a time window of several ns, and it is thus not feasible with standard QM/MM calculations. Using interpolation schemes for the potential energy surfaces¹²² makes this strategy viable. A further strategy involves running separate simulations, both for a single pigment and for an ensemble of pigments. Simulations on single pigments are employed to estimate the fast fluctuations of the site energies. Simulations on the ensemble are then used to extract the static disorder.¹²³ 

In Secs. III A, IV A, and IV B, we have presented some results obtained for the excitonic parameters of CP29. Here, we extend the analysis by comparing CP29 with a mutant and give a rationale for the observed changes in the absorption spectra of the two complexes.

The mutant that we are considering has been investigated experimentally by Guardini et al.²² It presents an asparagine (Asn) residue replacing the histidine (His) as the axial ligand of Chl a603 and, in the following, will be denoted as CP29-H111N for simplicity (we adopt here the same numbering of the referred article). This mutant is particularly interesting as it retains the targeted Chl a603, contrary to other mutations at the same H111 site, and it induces an increase in the red edge of the absorption spectrum, which has been connected to a change in the quenching properties of the complex. Authors suggested that this mutation results in a reduced distance/altered orientation between chlorophylls a603 and a609, induced by the lower steric hindrance of Asn as compared to His, and increased excitonic interactions between the two.²² 

MD simulations for both CP29-WT and CP29-H111N were employed for investigating the rearrangements occurring at the level of the Vio-a603-a609 cluster. The simulation protocol employed for the CP29-WT simulations has been presented elsewhere.³² For the CP29-H111N simulations, we have employed the very same protocol, with the only exception that the His residue H111 was mutated into the N111 residue before starting the minimization.

The reduced steric hindrance of Asn as compared to His results in Chl a603 being closer to its axial ligand in CP29-H111N as compared to CP29-WT. In Fig. 5(a), we show the distance between the Mg atom of Chl a603 and the Cα carbon of its axial ligand. The Mg-Cα distance is reduced by $∼14%$ in CP29-H111N as compared to CP29-WT (the mean distance is 6.64 Å in CP29-WT and 5.74 Å in CP29-H111N), confirming the rationalizations of Guardini et al.²² The a603 position appears to be very robust to local structural changes induced by the protein dynamics in CP29-WT, with a standard deviation of $∼0.13$ Å, and less so in CP29-H111N (standard deviation of $∼0.28$ Å; see also Fig. 5). Thus, the closer Chl a603 in CP29-H111N is more sensitive to local protein fluctuations than in the WT.

The same simulations have also been employed to investigate the relative position and orientation of the a603-a609 pair. When looking at the minimum distance between the two chlorophylls, we see that it is shorter in the mutant than in the WT (the average minimum distance is 3.54 and 3.38 Å in CP29-WT and CP29-H111N, respectively, which translates into a 4% decrease in the mutant compared to the WT). Differences in the a603-a609 cluster between CP29-WT and CP29-H111N, however, are not limited to their position, but also their relative orientation is affected by the altered position of a603 resulting from the mutation. Figure 5(c) shows the squared orientational factor κ_ij (i = a603, j = a609) appearing in Eq. (5), which we analyze here as a qualitative indicator of increased/decreased electronic coupling between the two pigments. The plot in Fig. 5 shows a 22% increase of $κij2$ in CP29-H111N as compared to CP29-WT (the average value of $κij2$ is 1.76 and 2.14 for CP29-WT and CP29-H111N, respectively) and a similar sensitivity to the protein dynamics (standard deviation of 0.25 and 0.22 for CP29-WT and CP29-H111N, respectively). The altered position and orientation of this chlorophyll pair can also be observed when looking at the MD [see Figs. 5(b) and 5(d) for a particular snapshot of the WT and the mutant].

The geometrical analysis of the trajectories strongly suggests that the electronic coupling between the a603-a609 chlorophyll pair will be enhanced in the mutant with respect to WT. Meanwhile, the altered positions of these chlorophylls in the mutant may lead to substantially modified chlorophyll site energies. Assuming for simplicity that the two chlorophylls only make up the lowest or second-lowest exciton in CP29 (which is a reasonable assumption according to previous quantum chemical calculations and spectroscopic results^77,124), both their site energies and their coupling are important for determining the resulting excitonic state. Site energies and couplings have been computed at the QM/MMpol level using TD-DFT M062X/6-31G(d) for both the CP29-WT and CP29-H111N, selecting uniformly spaced structures from the last 1 µs of trajectory for each replica, resulting in a total of 179 and 200 calculations for both chlorophylls for CP29-WT and CP29-H111N, respectively. The phytil tail of the chlorophyll was removed from the QM part (and considered in the polarizable MM system) by cutting the bond after the first aliphatic carbon of the tail. Additionally, in order to recover the environment contribution to the site energy $Eenv$ and extract the geometrical site energy $Egeom$⁠, we performed the same QM calculations without including the effects of the MM subsystem.

The site energy $Esite$ of each chlorophyll was then obtained by means of Eq. (19), which is essential to avoid double counting of the fast internal motions of the pigments when computing the spectrum, while still retaining the pigment identity arising from its particular dynamic environment. The site energies of Chls a603 and a609 are shown in Fig. 5(e). The difference in mean of chlorophylls a603 and a609 between the WT and the mutant is $Ea603WT−Ea603H111N=30±41$ cm⁻¹ and $Ea609WT−Ea609H111N=−18±43$ cm⁻¹ (the mean difference value and the associated standard error pooled among the replicas are reported). These values indicate that we can consider the site energy of both chlorophylls to be practically the same in the mutant and the WT. By contrast, the electronic coupling difference for the a603-a609 pair is estimated as $Va603−a609WT−Va603−a609H111N=−55±13$ cm⁻¹, which is around a 40% change in the absolute coupling of the WT [see also Fig. 5(f)]. Therefore, corroborating the geometrical insights obtained from the orientational factor [Fig. 5(c)], we obtain a significantly different coupling between the WT and the mutant, with the mutant presenting higher excitonic interactions than the WT. This confirms the intuition of Guardini et al.,²² who proposed a strengthening of excitonic interactions for the a603-a609 pair in CP29-H111N.

Considering the local character of the geometrical changes associated with the mutation (i.e., only the relative arrangement of a603 and a609 is substantially altered) and the increased electronic coupling between the two chlorophylls, we can compute the absorption spectrum of CP29-H111N by employing the very same exciton Hamiltonians employed for CP29-WT, with the only difference that the a603-a609 coupling is increased by 55 cm⁻¹ in the mutant. The absorption spectra of CP29-WT and CP29-H111N are then computed with the cumulant expansion theory truncated at second order (see Sec. IV). The spectral density of chlorophylls was modeled as a linear combination of 48 high-frequency modes and an overdamped Brownian oscillator term, taking the parameters from FLN experimental data.¹¹¹ This choice was made to remain consistent with our earlier optical spectra simulations of light-harvesting complex II (LHCII) and light-harvesting complex stress-related protein 1 (LHCSR1).^26,72

The calculated absorption spectra for CP29-WT and CP29-H111N are reported in Fig. 5(g). Notably, the red shift observed by Guardini et al. (from 674.8 to 685.6 nm) is well reproduced by our calculations (from 670.9 to 688.2 nm), demonstrating the effectiveness of our approach. Therefore, our computational model of CP29-H111N validates the interpretation of Guardini et al.,²² namely, that a strengthening of excitonic interactions of the a603-a609 pair is responsible for the red shift observed in the absorption spectrum. We also note that this computational approach does not just reproduce the experimental results, but it complements them by explicitly tracking down the molecular sources, giving rise to the different observables that have been probed experimentally.

In this Perspective, we have presented a general overview of the challenges encountered when light-harvesting complexes are modeled through atomistic approaches and highlighted the peculiarities that characterize these pigment–protein complexes. We have presented an effective computational strategy to simulate their optical spectroscopy, which combines molecular dynamics simulation, quantum chemical methodologies, and line shape theories. The predictive power achieved by such an integrated strategy has been demonstrated on the comparison of absorption spectra of wild-type CP29 and a particular CP29 mutant that has been studied recently by experimental means.

Despite the remarkable power offered by this strategy, there is still much room for improvement in many of the computational aspects discussed in this Perspective.

One possibility is to exploit machine learning methods to substantially reduce the computational cost associated with excitonic properties while retaining an accuracy almost equal to that of expensive quantum chemical methods.^125–132 Machine learning frameworks have also proven to be able to incorporate the description of the environment for complex systems, including pigment–protein complexes.^133,134 This opens a route to accurate calculations of static and dynamics disorders, which require a large number of calculations.

Another important aspect is the extent of MD sampling. In most applications, crystal structures or short MD simulations of the complexes are used, which can lead to significant biases toward the selected initial structure. Instead, to properly explore the whole conformational space of the complex, including all possible configurations of the pigments, enhanced sampling simulations should be used.^135–138 In addition, in this case, machine learning techniques can represent a speedup by providing useful reaction coordinates upon which the sampling is improved.^139–142 

The authors acknowledge funding from the European Research Council, under Grant No. ERC-AdG-786714 (LIFETimeS). The authors thank Ulrich Kleinekathöfer for providing the spectral densities calculated through a DFTB/MM MD simulation.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.