Quantum Chemical Modeling of the Photoinduced Activity of Multichromophoric Biosystems.

Multichromophoric biosystems represent a broad family with very diverse members, ranging from light-harvesting pigment-protein complexes to nucleic acids. The former are designed to capture, harvest, efficiently transport, and transform energy from sunlight for photosynthesis, while the latter should dissipate the absorbed radiation as quickly as possible to prevent photodamages and corruption of the carried genetic information. Because of the unique electronic and structural characteristics, the modeling of their photoinduced activity is a real challenge. Numerous approaches have been devised building on the theoretical development achieved for single chromophores and on model Hamiltonians that capture the essential features of the system. Still, a question remains: is a general strategy for the accurate modeling of multichromophoric systems possible? By using a quantum chemical point of view, here we review the advancements developed so far highlighting differences and similarities with the single chromophore treatment. Finally, we outline the important limitations and challenges that still need to be tackled to reach a complete and accurate picture of their photoinduced properties and dynamics.

Introduction

The synergistic relation between the technological progress in spectroscopic techniques and the advancement of the theoretical methods, coupled with the improvement in the available computational capabilities, have produced a tremendous boost in our understanding of multichromophoric systems, especially of biological origin. Ultrafast spectroscopic techniques developed in the last decades have been able to gather a huge amount of information on the response of these systems to light-induced perturbations and, indirectly, on the details of the excited-state dynamics. The possibility to explore events that occur at shorter and shorter time scales has revealed new phenomena and opened new questions: What processes of chemical and biological significance take place in the subpicosecond region?(1)What is the nature of the created electronic states? Are they delocalized over multiple chromophoric units?(2)What pathways are followed to relax/exploit the absorbed energy?(3)What is the role of the environment (both as a static scaffold and as an energy reservoir) in these activated energy transfer processes?(2,3)

Eventually, the advent of two-dimensional electronic spectroscopy (2DES) allowed to explore coherent phenomena. The observation of long living beatings in 2DES maps posed new fascinating questions: Are biological systems using a (robust) quantum-coherence pathway for electronic energy transfer?(4)Did Nature think of it first?(5)

Interpreting the outcomes in such complex spectroscopies is by no means immediate: overlapping and oscillating signals, energetic and structural disorder, and the variable strength of the coupling between the system and its environment make it very difficult to clearly assign the microscopic origin of the observed spectral features. As an example, the interpretation of the quantum beating pattern discovered in 2DES spectra of light-harvesting complexes, initially given as a signature of wavelike exciton transport, was later challenged by theoretical models and newer experiments[6] and remains an open ground for discussion.

New questions call for novel and more accurate theoretical modeling.[7,8] On a parallel track, the continuous evolution of computing power has allowed for more ab initio computational investigations, which can support the interpretation of the experiments. In this quickly evolving landscape of experimental evidence and theoretical interpretation, quantum chemical methods offer a reliable route toward the understanding of complex excited-state dynamics in multichromophoric systems.[9,10] This Focus Review aims at surveying the ample body of theoretical methods that have been developed so far to identify established concepts and to reveal still open challenges. In particular, the focus will be on biological systems, which still represent an extremely challenging problem for modeling as their excited-state relaxation cannot be described without including the multiple effects of the embedding complex biomatrix and its dynamics.

To graphically summarize this specificity, in Figure we represent the multiscale nature of excited-state formation and relaxation in three multichromophoric biosystems, two light harvesting (LH) pigment–protein complexes, namely the Fenna–Matthews–Olson (FMO) complex of green sulfur bacteria and the major LH complex of purple bacteria (LH2), and a DNA duplex. These three examples have been selected as representative of the multichromophoric biosystems family members’ diversity. FMO (as most of the LH complexes) can be seen as a rather rigid system (within the context of biosystems), with chromophores at a relatively large distance, presenting separated excited states with relatively long lifetimes and harmonic potential energy surfaces (PES). In some cases, however, LH complexes display much smaller interchromophores distances (as it happens in LH2), and the resulting excited-states can present characteristics in between standard LH complexes and the other type of multichromophoric biosystems, namely nucleic acids. These are extremely flexible systems, in which different conformations control different relaxation pathways. Localized, delocalized, and charge transfer (CT) states are all present at the same time. States lifetime are ultrashort, and PESs are all but harmonic.

Figure 1

Multiscale description of excited state nature and dynamics in multichromophoric systems. The Fenna–Matthews–Olson (FMO) complex of green sulfur bacteria, the LH2 complex of purple bacteria, and the DNA duplex are shown here as three remarkable, representative examples. “Left-to-right”: path to reduce the system complexity into smaller (tractable) pieces. “Right-to-left”: path that brings from the single site properties to the aggregate properties.

The diversity between LH and DNA/RNA systems is reflected in the completely different fate of the light activation. In nucleic acids, the excess energy is dissipated as quickly as possible to prevent photodamage. On the contrary, the light harvesting apparatus of photosynthetic organisms has evolved to control the transport of excitation energy and its conversion into chemical energy. These behaviors will be reconsidered at the end of the review, where these two groups of systems will be compared in terms of their possible modeling.

The scheme reported in Figure can be read both from left to right and from right to left. On the one hand (left-to-right path), the goal is to investigate the behavior of complex multichromophoric systems to unveil their structural–functional properties. To do so, as usual, it is necessary to reduce the system complexity into smaller (tractable) pieces. The challenge here is to introduce sensible approximations, which preserve the physicochemical properties of the whole system while making it possible to perform calculations on it. This reductionist approach brings one to the computation of very accurate properties of single molecule or small groups of molecules, a task that quantum chemistry is nowadays able to tackle in many cases of interest. On the other hand (right-to-left path), the theoretician should reconnect the various pieces of this multiscale modeling, by introducing, at each step, additional complexity: the single site properties, computed in their environment, and the intersite coupling, are used to access aggregate properties (for example through the Frenkel Exciton model), and transport between exciton states has also to be modeled by including system vibrations. CT states may also be necessarily added to the description, and a sampling over many possible system conformations is usually required to reproduce ensemble characteristics. In this process of connecting the various scales of the problem one may realize that some of the introduced approximations are inadequate: the approximation of considering harmonic electronic potential energy surfaces (PES) can be dramatically wrong in the case of DNA bases. Attention should therefore be paid in the use of techniques in system with different characteristics, which can never be blindly transferred from one to the other.

For the following presentation, we have selected the “right-to-left path” but the different sections can be alternatively read in the reverse order, thus reconstructing the “left-to-right path”.

The Single-Site Problem

The necessary balance between accuracy and scalability of computer simulations generally prevents the possibility of a direct investigation of all the atoms of a large and complex biosystem within a single quantum mechanical region. To overcome this important limitation, the system is commonly partitioned into the site and an environment. In the present context, we define a site to be a fragment, a molecule, or a group of molecules that can be homogeneously described at the QM level of theory. The environment, instead, is what remains beyond the QM site system and is here treated at the classical level. There are many strategies developed so far to integrate QM and classical descriptions for the study of excited-state formation and relaxation. These will be presented and commented in the following subsections for the single-site problem, whereas their generalization to multichromophoric systems is established in the next section, where a discussion of what can be taken and what should instead be left (or further approximated) is also reported.

The Quantum Chemical Methods

The extensive body of quantum chemical methods designed for the computation of excited states[11] spans from single-configurational procedures,[12] such as the CI-singles (CIS), CI-singles and doubles (CISD), and coupled-cluster (CC),[13,14] to multiconfigurational treatments[15] as the complete (or restricted) active space self-consistent field, CASSCF (or RASSCF),[16,17] usually employed in tandem with second-order perturbative corrections (CASPT2 or RASPT2). Another family of methods is included under the denomination propagator approaches,[18] e.g., Green’s function, equation-of-motion, and algebraic–diagrammatic–construction (ADC) methodologies, which do not require to compute the wave functions of individual states to obtain excitation energies and transition probabilities. Among them, the combination of Bethe Salpeter equation with the GW approximation (GW/BSE)[19,20] represents a promising method.[21,22] The quality of these methods relies on the type of reference wave function. Finally, time dependent density functional theory (TDDFT) is nowadays among the most popular approaches, also because of its black-box fashion.[23] DFT/MRCI,[24] which combines DFT with multireference configuration interaction (MRCI), has shown great efficiency for calculating excited states.[25]

A single technique appropriate for all the photophysical problems of interest does not exist: there is always a trade-off between price and performance. TDDFT can treat rather large systems, with very accurate results, but significant issues have been raised about its capability of describing CT states[26] and double excitations. Moreover, near degeneracies in the ground state cause difficulties in TDDFT, given the single-reference character of the ground state.[27]

CC-based methods have proven to be the most accurate for small and medium size molecules in situations where the HF reference is a fair approximation for the ground state. The CASSCF/CASPT2 methodology is considered a reference for its accuracy and capability of describing in a balanced way states that are multiconfigurational in nature, and as such it represents the state-of-the-art method for the computation of excited-state properties in generic molecular systems. Nonetheless, it requires an accurate selection of the orbitals to be included in the active space (which is a completely non-black box procedure), and the cost of such computations rapidly increases with the system size.

The Classical Models for the Environment

Many different strategies have been proposed to integrate a classical model of the environment within a quantum chemical description of a solvated or more generally embedded molecule.

The most straightforward and also the most computationally efficient approach is to reduce the environment to an infinite continuum dielectric, which coarse grains all the atomistic degrees of freedom of the environment in an average picture represented by a dielectric constant. Among the many possible realizations of continuum models, the one which projects the polarization of the whole dielectric only on the surface of the cavity chosen to contain the QM system, has obtained the largest success. The polarization is here represented by an “apparent” (e.g., induced) surface charge (ASC) distribution and discretized in terms of point charges by introducing a surface mesh.[28] The ASC distribution, in turn, influences the QM system, effectively describing mutual polarization between the QM system and the environment. In this framework, the definition of the cavity (and its surface) represents the most important adjustable parameter. In the most effective implementation of the model, the cavity is defined in terms of interlocking spheres centered on the atoms (or selected groups of atoms) of the QM molecule.

Because of their simplicity, the implementation of continuum models in electronic structure codes is extremely fast: this, when combined with the easiness of use and the very good ratio between accuracy and computational cost, explains why they are the most used approaches to introduce environment effects in quantum chemistry. Of course, the characteristic of completely neglecting the atomistic nature of the environment can lead to poor descriptions, especially when specific interactions between the QM and selected atoms of the environment are important in determining the property or the process of interest. Continuum models can oversimplify very heterogeneous environments. Finally, its unique advantage of giving an implicitly averaged description of the environment can turn into a limit when its fluctuations occur on the same time scale as the investigated process, preventing a mean field approximation. In all these cases, a much better strategy consists in introducing an atomistic description of the environment.

This detailed (atomistic) description of the environment is achieved through molecular mechanics (MM) force fields, and the resulting QM/MM approach has been formulated in various ways with respect to the type and the number of interactions explicitly included in the Hamiltonian determining the QM molecule.[29] In the “mechanical embedding” QM/MM approach, these interactions are null, as they are all rewritten in terms of additional MM terms. As such, this approach is only applicable if one is not interested in the effect of the environment on the electronic density of the molecule (and all the related properties), and it is therfore completely useless in the present context of photoinduced processes. In the “electrostatic embedding” formulation of QM/MM, the effects of the MM atoms are included in the Hamiltonian of the QM system as a one-electron operator having exactly the same form of the electron–nuclei Coulomb interaction. This formulation, more accurate and complete then the previous one, is nowadays the most commonly used. More refined formulations, which go beyond the electrostatic embedding, are those that include polarization effects in the MM part and/or add explicit operators for the QM-classical dispersive and repulsive interactions. While the latter extension is still not widespread, the polarizable reformulation of QM/MM, known as “polarizable embedding”, is becoming more and more widespread. Different alternative formulations of the polarizable embedding have been successfully applied to describe spectroscopies of embedded molecules and photoinduced processes, where mutual polarization between the QM and the classical system can lead to non-negligible effects.[30−34]

Both QM/continuum and QM/MM approaches, in their different flavors, have been extended to describe excited-state properties and processes. In this context, the dynamical response of the environment is characterized by different time scales: the inertial component of the polarization, which is due to motions and rearrangements of the environment atoms, will present a much slower response time than the electronic component, which is instead characterized by an almost instantaneous reaction of the electronic density. To translate this into classical models, continuum or atomistic, these two fast and slow components have to be explicitly defined and will respond differently to the fast change of electronic state in the QM system during, e.g., an absorption or an emission process: only the fast component is allowed to completely and immediately relax following the quick change in the QM electronic density, while the other component will be kept frozen in its initial configuration. This is exactly what is described by the so-called nonequilibrium regime. In continuum models, this is achieved by separating the optical component of the dielectric constant from the full (static) one: the apparent charges will also split into a dynamical and a inertial component which are allowed to react differently to fast changes in the QM system. The same separation is naturally obtained in polarizable QM/MM: the nuclei are explicitly considered and represent the slow component, while the fast component is accounted for through the polarizability. On the other hand, electrostatic embedding can only account for the slow response, because MM charges are fixed at the positions of the classical nuclei.

The Potential Energy Surfaces

In the Born–Oppenheimer approximation, we can model the process of a chromophore interacting with light with a nuclear wavepacket moving from the bottom of the ground-state (GS) potential energy surface to one of the possible excited-state (ES) PESs. The Franck–Condon principle establishes what are the vibrational levels which are likely to be populated in the process: these must be instantaneously compatible with the nuclear positions of the vibrational level of the molecule in the starting electronic state.

Most commonly, the nuclear degrees of freedom are approximated as harmonic oscillators. In a further simplification, called displaced harmonic oscillator (DHO), the coupling between the electronic and nuclear degrees of freedom is assumed to be linear (Figure ). The excited-state PES has thus the same curvature as the ground-state PES but a shifted equilibrium position. Within this model, all the parameters describing the various PESs can be readily computed at the initial geometry: a frequency computation in the GS minimum gives normal modes and frequencies, while the energy gradient on the ES PES at the Franck–Condon point is directly linked to the displacement between the two PESs minima as well as to the vibronic coupling along each normal mode. This frequency-dependent vibronic coupling can be exactly mapped into a spectral density function,here written on the basis of normal-mode frequencies ω, and the Huang–Rhys factors S, which describe the dimensionless displacement of the excited-state PES along the kth normal mode. The total strength of the vibronic coupling is measured by the reorganization energy .

Figure 2

(a) Drawing of GS and ES PESs, highlighting important points and transitions. The displaced harmonic oscillator (DHO) model is shown in the gray shadowed inset. (b) Relation between the energy gap fluctuation autocorrelation function C(t), the spectral density J(ω), and the DHO model (FT stands for Fourier Transform). The mode responsible for the high frequency peak in the spectral density (highlighted in red) is also depicted on the top of the molecular structure and highlighted in C(t). Similarly, a low frequency mode is highlighted in green and connected to the slow motion of the biomatrix embedding the chromophore. The peaks in J(ω) have position determined by the mode frequency and height determined by the strength of the coupling to the electronic excitation (the two parameters, ω and de, of the DHO model). The effect of the different modes (high and low frequency) on the linear spectrum is shown on the side of the respective parabolae: high frequency modes produce a band structure, while low frequency modes are responsible for the so-called homogeneous broadening.

The spectral density can be evaluated knowing only the GS nuclear trajectory, as the autocorrelation function C(t) of the excitation energy fluctuations (Figure b). If the ground-state trajectory is described classically, one applies a quantum correction to the Fourier transform of the classical autocorrelation function Ccl(t):[35]

This description allows consideration of all the vibrational modes in one function comprising the motion of the environment surrounding the chromophores. Even though the spectral density only describes linear coupling to the vibrational degrees of freedom, it can be used to approximately map any kind of vibronic coupling to an infinite set of harmonic oscillators.[36]

In real systems, PESs can be much more complex than the previously described independent harmonic wells (Figure a). They can have arbitrarily complicated shapes, and the Born–Oppenheimer approximation may be violated in regions where the electronic states become too close: internal conversion (IC) or intersystem crossing (ISC) may take place, respectively, between states of the same or different spin multiplicities. As the evaluation of the entire hypersurface is an impossible task, one usually computes only significant points along the so-called minimum energy path (MEP), such as minima, saddle points, barriers, and crossings between the states. These can be located by suitable optimization algorithms.

A general protocol (valid for all kinds of valence excited states) employed for that purpose relies on CASSCF determined geometries, with energetics refined by including dynamic electron correlation with CASPT2 single-point computations on the specific nuclear configurations. Attention should be paid here as state swapping may occur when the CASPT2 correction is applied to structures determined at the CASSCF level. Moreover, in order to describe the entire reaction path, the selection of the active orbitals requires them to describe the complete reaction mechanism, as the active space must not change along the reaction coordinate.

If one focuses on states for which TDDFT is known to work properly, the same kind of mapping can be performed much more easily and rapidly at this level of theory. In practice, the use of TDDFT follows the same caveats outlined above on double excitations and charge-transfer states. Importantly, even if the Franck–Condon region is accurately described by TDDFT, at other geometries the excited states might acquire charge-transfer or doubly excited character and would be poorly described.[23,27] Finally, a linear-response approach such as TDDFT cannot correctly describe the regions where one excited state crosses the ground state.[27]

The Excited State Dynamics

The mapping of the PESs provides important structural information, and the computation of spectroscopic signals at selected points along the reaction path gives snapshots of the system evolution. Accounting for the nuclear motion is nonetheless important, as the addition of the nuclear kinetic energy can change the reaction mechanism deduced from the sole mapping: it can for example open new deactivation paths, allowing the molecule to visit regions of the PES unaccessible for a “cold” wave packet. Moreover, the nuclear momentum enters in the evaluation of the terms that describe the nonadiabatic behavior of the system dynamics in proximity of PESs crossings, the so-called nonadiabatic couplings. The proper description of the system relaxation from the initially excited state to the final product state(s), and of the related evolution of the spectroscopic signatures, therefore calls for suitable modeling(s) of the nuclear dynamics.

The nuclear dynamics can be incorporated at various levels of approximation[37,38] (from the most accurate and computationally expensive, to the less accurate and cheap methodologies), employing:

Quantum Dynamics (QD) methods, which explicitly solve the time-dependent Schrödinger equation for the time propagation of the quantum nuclear wave packet on the coupled PESs and can provide numerically exact results. Multiconfiguratinal TD Hartree (MCTDH)[39] is considered as a reference QD method, providing high accuracy PESs and nonadiabatic couplings are given in input. In this respect, a major drawback in these methodologies is that they require precalculated PESs, necessarily restricting their applicability domain to a few electronic states and/or few effective (important) coordinates.

Mixed Quantum Classical Dynamics (MQC) methods are based on a separation between quantum and classical system degrees of freedom. In most of MQC approaches, electrons are treated quantum mechanically, while the nuclei move according to classical equations of motion. To enforce self-consistency, the quantum and classical subsystems are connected through nonadiabatic coupling terms. The MQC split enables the treatment of realistic molecular systems because of the reduced computational costs. Moreover, the MQC treatment allows computing of potential energies, energy gradients, and couplings on-the-fly during the trajectory integration (local approximation), with a significant impact on computational costs because precomputed multidimensional surfaces are not required anymore. The most established methods are Ehrenfest dynamics (ED), trajectory surface hopping (TSH),[40,41] and multiple spawning (MS).[38,42] In the first method, the classical nuclei move on a mean potential given by the weighted average of the PESs according to their coupling; in the second method, a swarm of classical nuclear trajectories evolves on an adiabatic surface with the possibility to jump to other surfaces; in the MS method, the semiclassical dynamics of nuclear wavepackets is captured by classically driven Gaussian functions, and the classical nuclear trajectories are used as an auxiliary grid for a quantum propagation of the nuclei. While the classical localization of nuclear trajectories allows treating realistic systems, it also makes it impossible to provide a description of quantum phenomena depending on global features. These purely nuclear quantum effects, as, for instance, tunneling and coherent spectral oscillations can be partially recovered in the MS treatment.

The Spectroscopic Investigation

The largest part of the knowledge we have about photochemical and photophysical properties of molecules comes from the combination of steady-state and time-resolved techniques. While steady-state spectroscopy delivers information about the energetic positions and probabilities of the observed transitions (typically from ground-state stable configurations/minima), time-resolved spectroscopy is able to provide kinetic information. Among the time-resolved techniques for the study of ultrafast electronic dynamics, pump–probe (PP) spectroscopy is the most widely used.[44] In a PP experiment, a fraction of the molecules is promoted by a first (strong) pulse, called pump pulse, to a resonant electronic excited state. After a time delay τ, a weaker probe pulse is sent to the sample to monitor the pump-induced changes in sample absorption. A difference absorption spectrum (ΔA) is then obtained by subtracting the absorption spectrum of the sample in the ground state to the absorption spectrum of the excited sample. The information contained in ΔA(λ,τ) comes from ground-state bleaching (GSB), stimulated emission (SE), and excited-state absorption (ESA) contributions. Very short laser pulses allow to follow extremely fast phenomena but at the cost of a reduced spectral selectivity.

The 2DES technique can be seen as an extension of PP: in PP, the evolution of the initially excited sample is monitored over a frequency (detection energy) and a time axis (delay between pump and probe pulses); in 2DES one resolves the signal over an additional coordinate, namely the excitation energy. This improvement allows revealing of connections between optical excitations at a given frequency, and the signals they create over a wide range of frequencies, disentangling possible overlapping signals, and eventually giving a direct view of the (potentially) multiple molecular transitions and competitive photo induced processes. Moreover, the introduction of an additional pump pulse, with a controlled time delay with the other pump pulse, makes it possible to circumvent the trade-off between time and frequency accuracy.[43] The relation that holds between PP and 2DES is explained in Figure .

Figure 3

(a) Pulse setup and time delays in 2DES (noncollinear geometry, which allows having a background-free signal). (b) Computer spectroscopic simulations are based on the system response function R. (c) Relation between PP and 2DES and wealth of information gained in 2DES at different waiting times t2. Excitation axis, only resolved in 2D, is labeled with ω1, while detection axis, present in both techniques, is labeled with ω3. In 2DES, two types of peaks are distinguished: diagonal peaks that mirror the linear absorption spectrum and thus highlight the bright transitions from the ground state which lay within the considered spectral window (GSB), and cross-peaks, displaced outside the diagonal. At t2 = 0, cross-peaks reveal possible SE on the red side of the bright transitions, which may red-shift and/or disappear while the system evolves along increasing waiting times (t2 > 0). Peaks of opposite sign (with respect to GSB and SE) on both diagonal or off-diagonal positions can demonstrate the presence of ESA signals and the appearance of new ESA features at increasing waiting-time t2 can indicate the production of new states through, e.g., internal conversion or intersystem crossing. The ratio of the diagonal to antidiagonal widths reflects the degree of inhomogeneous versus homogeneous broadening, which, in contrast with linear absorption, are here resolved independently.[43]

To simulate these spectroscopies and explain their outcomes, a field-matter Hamiltonian is introduced, and the molecular coupling to the radiation field is usually expressed as , where is the dipole operator of the QM system and (t) is the classical incoming field (semiclassical treatment). Time-dependent perturbation theory is usually applied (in terms of H′(t)), and different orders of the perturbative expansion give rise to different kind of signals: linear techniques are studied with first-order signals, while PP and 2DES generate third-order signals.

A density matrix approach for describing the system changes during and between the interaction with the pulses is usually preferred to a wave function based approach, as the density matrix formalism allows including the fluctuations and dissipation processes that occur due to the interaction between the system and the environment. The time dependent polarization induced by the interaction of the ensemble of molecules with the incoming radiation is the quantity of interest for the simulation of spectroscopy, as it gives rise to the measured signal of spectroscopic experiments.

First-order (third-order) polarization connects the incident radiation field to the so-called system response function, whose expressions are given by[45]where is the equilibrium density matrix (i.e., the density matrix prior to any field interaction) and , represents the interaction of the QM system with the laser. The superoperator drives the evolution of the density matrix for a time t, in the absence of an external field: .

The response function R(1) (R(3)) summarizes the evolution of the first-order (third-order) perturbed system in the following way (by reading eqs and 4 from right to left): (i) starts from the equilibrium density matrix , (ii) interacts with the laser via the dipole superoperator μ×, (iii) evolves the perturbed density matrix with according to the field-free molecular Hamiltonian H0 for a time t, (iv) repeats the above two steps as prescribed by the considered techniques, and (v) computes the macroscopic polarization with the perturbed and evolved density matrix.

From a practical point of view, the necessary ingredients to simulate spectroscopic experiments are the energies of the system’s manifold of states, the transition dipoles between the states, and the line-width broadening. The latter can be provided phenomenologically or by explicitly accounting for the system–bath interaction through the spectral density associated with the various transitions. The cumulant expansion of Gaussian fluctuations (CGF)[45] method is a mean to describe the effect of the environment on the transition line-shape beyond the phenomenological treatment. It is exact for a bath characterized by Gaussian fluctuations (and in the absence of transport of excitation between the various states of the system). Within this model, linear and nonlinear spectra are computed through the line shape function g(t), which is directly obtained from the spectral density:[45,46]

Population relaxation can be included at some approximate level, usually by assuming a separation of time scales between the bath fluctuations and the transport processes[45] or relying on a higher level description of the density matrix evolution which explicitly evolve the coupled electronic and nuclear DOFs.[47]

From Single-Site to Multichromophoric Systems

In the previous section, we have highlighted the achievements of the quantum-chemistry based modeling of systems made of a single QM unit (the site) embedded in a classical environment. To what extent can we use such achievements to describe systems where the number, and dimensions, of the QM units is too large to treat homogeneously the entire system? What should be taken of the site description, and what should instead be abandoned?

On the one hand, a subdivision of the aggregate into a system and environment part can still be performed, and the coupling between the electronic and the nuclear degrees of freedom can still be accounted for in order to describe the system dynamics after photoexcitation. On the other hand, in these systems one has to consider, in principle, a large number of electronic states and of nuclear degrees of freedom. This, coupled to the size of the system, forces the reformulation of some of the previously introduced concepts and the introduction of additional approximations:

Quantum Chemical Description: The multichromophoric nature of the system implies a too large number of atoms to be treated homogeneously, and the electronic Hamiltonian is generally reconstructed in terms of independent but interacting units. However, the nature and the dimension of the most proper units is not unequivocally defined and it might not coincide with the single chromophores.

Environment: The classical models used for single site systems can be generalized to multichromophoric systems, but now their responsive nature becomes much more important, as the QM Hamiltonian explicitly includes interunit couplings which can significantly change if the environment is allowed to polarize. Moreover, the issue of the sampling of the environment is even more delicate, as now each chromophore has its own “local” environment and the interpigment couplings connect all of these environments.

Coupling between Electronic and Nuclear Degrees of Freedom: The excited states of a multichromophoric aggregate are coupled to many nuclear coordinates of several chromophores, which complicates the description of their PESs. The dependence upon the intramolecular nuclear coordinates can still be included by considering the vibronic coupling of each site with the DHO model, restricting the form of the total multichromophoric Hamiltonian (see section ). This assumption greatly simplifies the problem, as many quantities of interest can be computed analytically. However, real potential wells can drastically deviate from the harmonic approximation (as it happens, e.g., in DNA). Moreover, local baths of different units can be correlated, and the interunit couplings can be modulated by large intermolecular motions.

Excited State Dynamics: Contrary to the previously described dynamics, which occurs within a single site, in multichromophoric systems the relevant photophysical processes span different sites. Excitation energy transfer (EET) and/or transfer of electrons (electron transfer, ET) among chromophores have to be modeled. Several theories have been developed, but their range of applicability is generally limited due to the specific regimes that each theory assumes.

Spectroscopy: Spectroscopy data becomes congested of contributions coming from the different chromophore and from their mutual interaction. Theoretical simulations are here even more necessary in order to disentangle different contributions and interpret recorded spectra.

The Quantum Chemical Description

The excited states of a multichromophoric system (such as a pigment–protein LH complex) may be effectively described introducing the Frenkel exciton model: one performs QC computations on single sites, i.e., subregions of the entire system, and computes site–site coupling terms to reconstruct the multichromophoric manifold of electronic states.

In the standard formulation of the model, the electronic excited states of the entire system are expanded on a basis of diabatic states |i⟩, which represent an excitation localized on a single chromophoric unit:It is commonly assumed that the states |i⟩ are Hartree products of the excited state of one chromophore with the ground states of all other chromophores (|0,···,1,0,···⟩). Physically, this corresponds to neglecting the GS exchange interactions between different chromophores. Consistently, the global ground state |g⟩ is represented as the Hartree product of all chromophores’ ground states.

The electronic Hamiltonian of the entire system is then written on the basis of states |i⟩, namely:where are the excitation energies of the isolated chromophores (site energies), and V are the electronic couplings between the different chromophores’ states. Eg is the energy of the ground state, which is taken as the zero of the energy scale whenever one is only interested in the excitation energies. The electronic energies and coefficients of the expansion in eq are obtained by diagonalization of the Hamiltonian matrix. Here, for the sake of simplicity, we have assumed that each chromophore i contributes with only one excited state to the diabatic basis set, but the essential details of the model do not change when considering more than one excited state per chromophore.

The power of the exciton model lies in the possibility of computing its parameters and V only from calculations of single-chromophore properties. In particular, the electronic coupling V between two excited states localized on different chromophores can be very well approximated as the following Coulomb interaction:[48]where ρ0(R) (ρ0(R′)) is the transition density corresponding to the excitation at site i (j). The integral in eq can be evaluated numerically or analytically using a basis set expansion.[49] An alternative derivation through natural transition orbitals has also been proposed.[50,51] Approximated forms of this interaction, including the point-dipole approximation (PDA) and a point-charge expansion, have found widespread use.[52] The PDA coupling is the first nonzero multipole order of the interaction (eq ), and can be expressed aswhere μ is the transition dipole for transition i or j, R is the vector distance between the centers of transitions i and j, with magnitude R. It is well-known that the point-dipole approximation breaks down at close interchromophoric separation, namely, when R is smaller than the dimensions of the chromophores. Nonetheless, because of the simplicity of its implementation, the PDA has been widely employed to compute Coulomb couplings, even outside of its range of applicability.[52]

An alternative simple expression for the Coulomb coupling is the expansion on atomic charges. The densities ρ0(R) can be approximated by a set of point charges {qαTr} centered on the chromophore’s nuclei. Within this expression, the integral in eq is replaced by a double sum:where indices α and β run on the atoms of i and j. The transition charges qαTr can be obtained through a population analysis of the transition density ρ0(R). In analogy to the parametrization of classical force fields, the transition charges can be fitted to match the electrostatic potential generated by the transition density, in the so-called TrEsp method.[53]

It is possible to think at the exciton model as a simplified CIS, in which the molecular orbitals are now occupied and virtual orbitals localized on the sites: CIS configurations will consider transitions both within the sites (locally excited, LE, states) and between the sites (charge transfer, CT, states).

The excitonic ansatz of the standard exciton model (eq ) only considers LE states, namely, states where the excitation is localized on single chromophores. As such, it neglects all the CT states |i+j–⟩, where the hole and electron reside on two different chromophores i and j. Generalizations of the Hamiltonian (eq ) to include all CT configurations are however possible:

Here, denotes the energy of the CT state |i+j–⟩, is the coupling between a CT state and a LE state, and is the coupling between two different CT states.

The Coulomb expression of eq is not valid for couplings involving CT states. The coupling between a LE state and a CT state can still be computed from the properties of the noninteracting chromophores, although it requires computing the Fock matrix elements for each pair of chromophores.[54] Similarly, the energy of a CT state can be computed by knowing the orbital energies of chromophores i and j and Coulomb and exchange integrals.

The number of CT states in a system of N chromophores is in the order of N2 (while the number of exciton states grows linearly with N). This raises the complexity of the CT-extended exciton model, and the number of couplings to compute. The CT-extended exciton model would thus have the same complexity as a CIS performed on the entire system. However, only those couplings between states that have the hole or the electron on the same chromophore are nonzero. Moreover, the CT states where hole and electron reside on distant chromophores can be neglected. This makes the charge-transfer extended exciton model still computationally more tractable than a calculation on the entire system, especially for effectively one-dimensional systems.[55] Moreover, separating locally excited and charge-transfer states can be an effective strategy to correct for the unbalanced treatment of CT states by TDDFT methods.[54]

An alternative approach to the calculation of CT energies and couplings is based on the diabatization of the electronic Hamiltonian of dimeric units. This strategy is particularly effective when selecting only few CT states, or when chromophore pairs are related by symmetry.[55,56]

The Environment Effects

In the single site section we have summarized the most common approaches used in quantum chemistry to include the effects of the environment. As a matter of fact, the same approaches can be (and have been) extended to multichromophoric systems and used to calculate the excitonic parameters (site energies and couplings) needed to build the excitonic Hamiltonian of the embedded system.[10] Both continuum and MM classical models can be applied to describe the changes induced by the environment on the excitation energies of the single chromophoric units, exactly in the same way they are used for single site systems. Some specificities are instead needed to predict the “correct” coupling between embedded chromophores.

The electronic couplings are affected by the environment through two mechanisms. First, the environment can change both the geometrical and the electronic structure of the chromophores and modify their transition properties, i.e., transition dipoles and transition densities. These changes will be “implicitly” reflected in a change of the Coulomb coupling (eq ), which will be generally enhanced due to the electrostatic effect of the environment. The second effect is due to the polarizable environment, which mediates the interaction among chromophores’ excitations.[57] The resulting “explicit” effect generally reduces the magnitude of the direct (Coulomb) coupling. For this reason, it is common to say that the coupling is “screened” by the environment.

All QM/classical methods outlined in section automatically give the first implicit effect within the same formalism used to get the “new” site energies. It is instead more critical to include the explicit effect. The simplest way is to introduce a screening factor s such that V = sVc, where Vc is the direct electronic coupling between chromophores with the implicit effect taken into consideration. The screening factor s can be related to the inverse of the optical dielectric constant ϵ∞, which is generally approximated with the square of the refractive index of the medium. For a refractive index of 1.4 (typical for a hydrophobic region of protein environments), the screening factor is ∼0.5, namely, the electronic coupling is reduced by a factor ∼2. This is indeed what was originally proposed by Förster in his famous model for EET (vide infra) where the PDA coupling of eq was divided by ϵ∞.[58,59] Within this approximation, the screening neither depends on the interacting chromophores nor on their relative orientation and distance. Moreover, at short distances, the dielectric medium can be excluded from the intermolecular region, leading to more complex effects.[60,61] In particular cases, this can also enhance the coupling.[57] To achieve a more accurate description, the apparent surface charges (ASC) formulation of continuum models we have described in section can be introduced. Within this framework, an explicit coupling term has to be added to the direct one. This additional term becomes[49]

Conceptually, the electronic transition in chromophore k drives a response in the polarizable medium, which, in turn, affects the transition in chromophore j. It is important to note that the apparent charges q are calculated using the optical permittivity of the medium, in order to account for the fact that only the electronic component of the polarization can respond (nonequilibrium regime). A very similar strategy is used in the Poisson–TrEsp method,[62] an extension of the TrEsp method described in section .

If one wants to keep the atomistic detail of the environment, MM descriptions can be used as well. However, the MM model has also to be polarizable to generate the analogue of VASC. An expression has been proposed based on the induced dipole (ID) formulation of polarizable embedding. Within this approach, an atomic polarizability is added to the fixed charge to describe each atom of the MM environment and the explicit contribution to the couping becomes[30]where the transition density ρ0 induces a response in the environment which is represented by the induced dipoles μ. As in the case of the QM/continuum, ρ0 here is also calculated self-consistently with the polarization of the environment. A mixed continuum/atomistic strategy has also been proposed:[63] in this case, VASC and VMMPol sum up and both terms are obtained in a fully polarizable scheme.

Before concluding this section on environment effects, it is important to introduce a general problem that will be better analyzed at the end of this Review. Namely, in the case of any atomistic description, the QM–environment interactions, particularly those at short-range, depend critically on the configuration of the environment. Therefore, several configurations of the whole system need to be taken into consideration to get a correct sampling. Obviously, this sampling is not needed when a continuum approach is employed because it implicitly gives a configurationally averaged effect due to the use of macroscopic properties.

The Coupling to Vibrations

In the previous sections, we have considered only the electronic part of the system Hamiltonian. However, vibronic coupling is fundamental in determining the excited-state dynamics of a multichromophoric system.

When the coupling between excitation and vibrational modes is strong, the electronic basis used in eq might not be the optimal choice for writing the total Hamiltonian. When one or few vibrational modes are strongly coupled to the site excitations, it is possible to include all the relevant vibronic states in the exciton Hamiltonian, thus treating them on the same level as the pure electronic states. This vibronic basis bypasses the need of nonperturbative theories to deal with large vibronic couplings because the most coupled modes are already included in the “system” Hamiltonian, while the remaining weakly coupled modes can be safely treated with perturbative theories.[64]

The details of vibronic exciton models are discussed in a recent review, with a focus on molecular aggregates.[65] Here, we present the main characteristics of the vibronic exciton Hamiltonian. For simplicity, let us consider a single high-frequency mode with frequency ω0 coupled to the electronic excitation. The vibronic exciton Hamiltonian can be expanded on the basis of the vibronic states of each chromophore, , where is a vector of vibronic quantum numbers for each chromophore. This basis comprises a single electronic excitation at site i and multiple vibrational excitations, and its dimensions grow exponentially with the number of chromophores.[65]

Normally, no more than two vibrational excitations have to be considered in the basis. In the two-particle basis set, one considers all the vibronic states where the excitation is localized on chromophore i, both chromophores i and j are in their excited vibrational states and v′, while the other chromophores are in their electronic and vibrational ground states. Here, (v′) is the vth (v′th) vibrational state of the excited-state (ground-state) PES of chromophore i (j).

The off-diagonal elements of the vibronic exciton Hamiltonian in the two-particle basis read:where and are the Franck–Condon integrals between vibronic states of chromophores i and k. In the DHO model, these integrals have a simple analytical expression in terms of the Huang–Rhys factors of the harmonic mode. It can be seen that explicit consideration of vibronic states scales the effective electronic interaction between two chromophores by the Franck–Condon integrals of these chromophores. The number of vibronic states is in principle infinite, but the basis is safely truncated when Franck–Condon factors become negligible.

A further approximation to the vibronic exciton Hamiltonian is the one-particle basis, in which only one vibrational excitation is considered for each electronic excitation. For each chromophore i, its vibronic excited states are considered, whereas all the other chromophores are in their ground electronic and vibrational state. This approximation was sufficient to explain the vibronic enhancement of energy transfer dynamics in a cyanobacterial light-harvesting system.[66]

The vibronic exciton Hamiltonian has proven especially useful for polyene molecules, which show large Huang–Rhys factors for a few vibrational modes. Explicit consideration of the carotenoid high-frequency mode in a simplified two-particle model was needed to explain ultrafast carotenoid-to-bacteriochlorophyll energy transfer in the natural antenna complex LH2.[67]

An alternative or complementary approach to treat the coupling with nuclear degrees of freedom, based on the perturbative treatment of vibronic couplings, is based on the formalism of the spectral density outlined in section . Within this formalism, the vibronic coupling of each site j is entirely specified by the fluctuations of the excitation energy . In line with the perturbative approach, one can assume that the coefficients of the expansion in eq are essentially independent of the nuclear degrees of freedom. In this limit, the fluctuations of the exciton energies E are determined by the fluctuations of the site energies and by the excitonic coefficients, and the corresponding spectral density is given bywhere it is also assumed that the site energy fluctuations are uncorrelated. This is usually well justified, because fluctuations in site energy are mostly determined by internal vibrations of the chromophore, which do not influence the site energies of other chromophores. Moreover, fluctuations of the electronic couplings have been assumed to be negligible.[68]

In the limit of N chromophores with identical spectral densities J(ω) = Jsite(ω), the excitonic spectral density J(ω) is reduced by a factor , which is always ≤1 and is smaller the more the excitation is delocalized. The spectral density of an excitation perfectly delocalized over N chromophores is then Jsite(ω)/N. By changing the electronic basis to the adiabatic exciton basis, the vibronic coupling is thus renormalized. Moreover, the energies of the exciton states become correlated because now the vibrations on one chromophore modulate the energies of several exciton states. Finally, off-diagonal vibronic coupling terms appear in the exciton basis. These nonadiabatic vibronic coupling terms are responsible for transitions between different electronic states and therefore mediate the EET process.

The Excited-State Dynamics

The theoretical research on the excited state dynamics in multichromophoric systems has a long history and many different approaches have been proposed so far. Here, it is convenient to introduce the formalism used in the open quantum systems field. The starting point is the definition of the total Hamiltonian as the sum of the electronic Hamiltonian and two additional terms, referring to the so-called “bath” (B) and “system–bath” (SB) interaction, respectively.

is commonly represented by the exciton formulation of eq or 11, whereas the “bath” indicates the nuclear degrees of freedom within the single chromophoric units and with the surrounding environment (see Figure ). Now, if the bath is treated quantum-mechanically, one enters in the formalism of the quantum master equation (QME), where it is commonly modeled as an infinite set of quantum harmonic oscillators. The bath and the system are assumed to be distinguishable and all the needed information can be extracted from the reduced density operator ρ(t) obtained after tracing out the bath degrees of freedom. As a matter of fact, with the term QME we refer to many different theoretical approaches that can be classified according to the completeness and the accuracy of the resulting picture. Starting from the first step of this ladder, one encounters the well-known Förster model for EET.[58,59]

Figure 4

System, bath, system–bath terms in the total Hamiltonian . The GS energy Eg is safely assumed to be zero. and represent, respectively, nuclear momenta and coordinates operators; g are the coupling constants between nuclear and electronic degrees of freedom (linear coupling). The perturbative terms (PT) that enter in the various transport theories are highlighted in the colored boxes.

This framework assumes that the electronic couplings between site excitations are small compared to the coupling of the electronic degrees of freedom with the bath (this assumption is also called “weak coupling” limit, see Figure ). In this limit, the dynamics can be described with a system of kinetic equations, whose rate constants describe the transfer of the excitation energy from one site (the donor j) to the other (the acceptor k), and are given by the Fermi Golden Rule expression:where V is the electronic coupling and F(E) and A(E) are the Franck–Condon weighted density of states (FCWD) for the emission of the jth chromophore and absorption of the kth chromophore. The integral in eq , referred to as the spectral overlap, can be recast using the properties of the Fourier transform:where λ is the reorganization energy of the donor state, g(t) and g(t) are the line shape functions of the donor and acceptor, and ω is the difference between their vertical excitation energies.

The Förster formulation describes an energy transfer mechanism occurring by incoherent “hops” between chromophores.[69−72] However, when the electronic coupling increases, the perturbative Fermi Golden Rule expression is not valid anymore,[68,72] as the electronic coupling tends to delocalize the excitation over many sites, such as it is observed in DNA or light-harvesting pigment–protein complexes, giving rise to exciton states.

In some conditions, however, it is still possible to adopt a generalization of the Förster formulation: this is when the multichromophoric system can be divided into weakly connected aggregates, each containing multiple strongly coupled chromophores (their coupling being strong enough to delocalize the excitation). Within this framework, also known as Generalized Förster (GF), the EET rate can be recovered from a generalization of the Förster eq by substituting the local excitation j and k with new “exciton” states, each delocalized in different aggregates.[73−75] As a result, the coupling V is replaced by the effective coupling between excitons L and M, i.e., as , whereas the effective spectral overlap is expressed through the excitation energy difference ℏω = (E–E), and through the line shape functions (g(t)) and reorganization energies (λ) of the exciton states.

A reversed perturbative formulation is instead needed when it is not possible to partition the system into weakly interacting components. In the “strong coupling” regime, the electronic coupling is included explicitly, while the system–bath interactions are treated as the perturbation (see eq ). As a result, the excited-state dynamics of the system consists in relaxation between exciton states instead of hopping between localized excitations. The corresponding relaxation rate can be expressed according to the Redfield theory as[68,71]where M and L are two exciton states and J(ω) is the spectral density of site j evaluated at the energy gap between exciton states M and L.

The standard Redfield (sR) rate eq describes a relaxation process where one quantum of electronic energy ℏω is released to the vibrational degrees of freedom through the vibronic coupling represented by the spectral density; as such it cannot treat transfer processes where multiple energy quanta are involved.[71] A possible correction is represented by the modified Redfield (mR) theory, where the diagonal part of is included in the zeroth-order Hamiltonian while the perturbation part includes only the off-diagonal part of (see 4). By treating the diagonal part of nonperturbatively, bath reorganization and pure dephasing are no more approximated in this formalism. Moreover, multiphonon effects are included. Within this framework, the relaxation rate becomes[68,71,76,77]where V describes the interaction between the donor and acceptor exciton in terms of the reorganization energies and the line shape tensors and their first and second derivatives. This definition is valid for arbitrary delocalization of the donor and acceptor states. For example, if the two involved states are localized at the jth and kth sites, then V(t) is constant and reduces to the square of the interaction energy |V|2, corresponding to a weak coupling between the localized sites j and k. From this analysis, and the comparison of eq with the expression of the Förster rate (eq ), a parallelism appears evident, even though here the “coupling” term V(t) is time-dependent and cannot be factorized out of the integral. Modified Redfield theory has been successfully employed to model exciton dynamics in multichromophoric systems.[78,79] Given its validity for arbitrary mixing of states, eq has also been applied to model charge separation in the reaction center of Photosystem II.[80]

In contrast to the sR, here the system–bath is not supposed to be weak, but the corresponding displacements of the equilibrium positions of the nuclear modes are taken to be independent of the exciton wave function. This means that eq is valid if the exciton delocalization is controlled by static disorder rather than by a relaxation of the nuclear degrees of freedom (a phenomenon known as dynamic localization).[68,76,81,82] In the other cases, excited-state dynamics is likely better modeled by Förster theory, assuming complete localization of the excitation, or by GF theory, assuming partial localization of the excitation within excitonic aggregates. Modified Redfield and Generalized Förster theories can be employed together to model strongly mixed states as well as weakly connected clusters.[83] This is the case, for example, of the charge-separation process in Photosystem II.[84]

As said, the applications of these theoretical models to multichromophoric systems has a long history that is based on the integration of these theories with empirical parameters generally taken from experiments and/or simplified molecular models. For example, a large part of the application of these theories to light-harvesting pigment–protein complexes has been based on the use of empirical values for the site energies in combination with a dipole–dipole approximation of the electronic couplings still based on empirical estimates of the transition dipoles.[85] Moreover, empirical formulations have also been used to define the line shape functions (g(t)) and the reorganization energies (λ).[78,82,84]

Only in the past few years there has been a shift toward a fully computational estimate of the required parameters. This evolution has been made possible by the combination of quantum mechanical descriptions of the excitonic parameters (site energies and couplings) with classical molecular dynamics simulations for the derivation of line shape functions and reorganization energies (see section ). Moreover, the simulation has been further improved so to explicitly account for the effects of the environment on all these parameters; in particular, this has been achieved by the important progresses seen in the development of accurate hybrid QM/classical approaches that we have described in section .

All the methods described so far, however, are based on some perturbative formulation which can be safely used only in specific ranges of application. As a matter of fact, a number of methods have been devised to treat the quantum dynamics beyond the perturbative formulations. Among various higher-order methods developed for decades, the hierarchical equations of motion (HEOM) approach[86−88] has gained popularity recently and its results are considered as benchmark data by many researchers.[82,83] The HEOM approach is based on the idea that higher-order terms of the system–bath interactions can be accounted for by introducing a hierarchy of coupled auxiliary operators, an assumption valid as long as the spectral densities can be expressed as a sum of Drude–Lorentz contributions. This hierarchy of equations is in principle infinite and must be truncated to a finite depth for numerical calculations. Recent alternative formulations and implementations have increased the numerical robustness and scalability of the method.[89,90]

A simple alternative to the HEOM method is the explicit treatment of selected nuclear degrees of freedom within the “System” Hamiltonian , in the vibronic exciton model described in section , whereas the other degrees of freedom are treated within the Redfield model.[91] This “Vibronic Redfield” approach was shown to reproduce the dynamics obtained by HEOM,[64] and was employed to model the primery charge separation event in the reaction center of Photosystem II.[92,93]

The vibronic models, the HEOM methods, and the other QME ones, are still mostly limited to the harmonic oscillator model for the vibronic interaction or even to a specific form for the spectral density. In the case of HEOM, the computational complexity increases with the number of peaks in the spectral density function,[90] which limits the use of complex spectral densities.

To go beyond these limitations, it could be necessary to consider the fully anharmonic dynamics of the nuclear degrees of freedom. Such improvement can be achieved by dropping the quantum mechanical description of nuclear motions, within the mixed quantum-classical approaches such as surface hopping and mean field approximations. Within this framework, the bath coordinates are classical variables (Rcl(t)), and their time dependence arises from Newtonian dynamics dictated by the electronic PESs, whereas the dynamics of the excitonic degrees of freedom is given by the time-dependent Schrödinger equation. The extension of exciton Hamiltonian to mean field (Ehrenfest) models has been presented some years ago and applied to light-harvesting complexes.[94] On the contrary, only recently, the same extension has been presented within a surface hopping dynamics. A first surface-hopping description of dynamics and 2D spectra of an excitonic aggregate was presented a few years ago based on a simplified resolution of the time-dependent Schrödinger equation.[95] In such a formulation, in fact, the quantum subsystem is propagated in time assuming a constant Hamiltonian during the time step. This is conveniently solved in the (local) site basis, which involves a Hamiltonian diagonalization but does not require an explicit calculation of the nonadiabatic coupling vectors. The same approach was successively used to simulate dynamics and 2D spectra of LH2.[96]

The first complete implementation of excitonic Surface Hopping was presented by Martinez and co-workers.[97,98] A multitiered parallel exciton framework was used to carry out an “on-the-fly” nonadiabatic dynamics trajectory using Tully’s fewest switches surface hopping (FSSH) algorithm. In such a formulation, the exciton Hamilonian was using a dipole–dipole approximation of the electronic coupling and the effects of the environment were neglected. Going further along the same line, more recently Menger et al.[99] have presented a combination of surface hopping nonadiabatic dynamics with an exciton model that includes the interaction between the chromophores and an external environment through a hybrid QM/MM formulation using an electrostatic embedding scheme. In this implementation, the excitonic couplings are modeled using atom centered transition charges instead of point dipoles. In both implementations a TDDFT level of QM theory was used, where in the first implementation GPU-accelerated algorithms[100] were used.

Notably, MQC approaches require an expression for the gradients of V with respect to the nuclear coordinates. Although analytical expressions have been derived for the couplings,[51] approximate expressions eqs and 10 are easy to differentiate, especially if one neglects the geometric dependence of the magnitude of the transition dipoles or of the transition charges.

The results of the MQC approaches can offer important information for multichromophoric systems that are not accessible through QME and HEOM approaches. On the other hand, these schemes treat nuclear degrees of freedom completely classically and as such yield correct equilibrium population only in the high-temperature limit. Explicit quantum mechanical treatment should be important for coupled vibrations whose frequency is greater than the thermal energy. Nonetheless, numerical comparisons with the HEOM method have shown that FSSH performs surprisingly well in a range of parameters relevant for multichromophoric processes.[101−103]

Spectroscopy in Multichromophoric Systems

Optical absorption, time-resolved fluorescence, pump–probe, and 2DES, to name a few, have been extensively applied to explore multichromophoric systems, evidencing clear signatures of a cooperative optical response (band splitting, redistribution of the optical intensity, EET, and ET dynamics between the different components of the aggregate, etc.).

Circular dichroism (CD) spectra are extremely informative for multichromophoric systems as they can be used to obtain information about the three-dimensional disposition (distances and relative orientations) of the interacting chromophoric units. Regarding nonlinear spectroscopy, 2DES has been widely employed, especially in the IR–VIS spectral range. This has allowed studying light harvesting complexes and reaction centers,[104−107] leading to a number of breakthroughs in the field. 2DES cross-peaks give valuable information about excitonic coupling between transitions, and they reveal energy transfer pathways when growing along increasing waiting times (t2 > 0).[108] The EET dynamics can also be monitored by the appearance and disappearance of ESA signals along t2. Oscillations of the diagonal and cross-peaks along t2 are here of different origin: they can reflect either electronic, vibrational, or vibronic coherences.

One of the first and more studied examples was 2DES on Fenna–Matthews–Olson (FMO) protein complex from green sulfur bacteria,[104] allowing revealing of the EET cascade at 77 K. 2D spectroscopy was also recently employed to track the EET flow through the entire photosynthetic system of green sulfur bacteria from the chlorosome to the reaction center.[109] Besides FMO, LH complexes of purple bacteria were also extensively studied with 2D techniques in various spectral regions.[70] Here the presence of carotenoids, with their ultrashort state lifetimes and strong coupling to intramolecular vibrational modes, enriches the photo physics of the system, producing a number of overlapping spectroscopic features which make their interpretation highly nontrivial.[110]

When they were first discovered, long living oscillations in 2D spectroscopy of the FMO complex were interpreted as a signature of quantum coherence between donor and acceptor electronic states,[105] leading researchers to speculate about robust wave-like transport through quantum coherence, but doubts concerning this interpretation were later raised. This has produced a wealth of studies aimed at identifying the origin of these oscillations,[111−114] boosting the development of accurate theoretical models while remaining the subject of intense debate.[6]

From a computational perspective, the main ingredients of spectroscopy simulations within the exciton model are readily obtained from the coefficients of eq . For example, the absorption intensity of an excitonic state K depends on the excitonic transition dipole , where μ are the transition dipoles of the chromophores and the coefficients c provide the redistribution of the chromophore dipole strength. When moving to the intensity of CD signals, instead, the exciton rotatory strength has to be introduced, namely , where denotes the imaginary part, and the and are the operators of electric and magnetic moment vectors, respectively. As for the exciton electric transition dipoles, the magnetic ones also can be obtained using the site properties (the intrinsic magnetic moments, mint) and the exciton coefficients. Now, however, an additional term depending on the position R of the chromophore has to be considered to account for the gauge invariance, namely

In most cases, however, the intrinsic magnetic moments of the chromophores are negligible, and the following approximate Rosenfeld equation, involving only electric transition dipoles, can be introduced:where R is the vector distance between the centers of chromophores i and j.

In the spirit of the perturbative treatment of vibronic coupling, the line shape functions g(t) of the exciton states are obtained from the excitonic spectral densities of eq . The multichromophoric spectrum is thus computed as a sum over all exciton states (SOS), each with its energy, position, and line shape.[68] The same treatment can be extended to nonlinear spectroscopies.[46] To describe ESA signals in nonlinear spectroscopies, one has to consider higher-lying states of the multichromophoric system. Usually, these are easily constructed on the basis of doubly excited states |ij⟩ by diagonalization of a double-exciton Hamiltonian analogous to that of eq .[46,115] Higher-lying excited states of the chromophores themselves, which may fall in the same energy window contributing to the overall ESA signal,[110] can also be added, providing their energy location and dipole moments are accurately estimated by QC computations.

The SOS approach[45,46] treats the exciton states as independent: the response functions (eqs and 4) are expressed in terms of exciton states, and the population transfer between the states described by a rate equation such as the Redfield equation. When the excited-state dynamics is more complex, the SOS approach is not correct anymore. Indeed, the linear and nonlinear responses of the system depend on the details of the excited-state dynamics. Nonlinear response functions can be calculated with the HEOM method by explicitly considering the successive application of dipoles operators on the reduced density matrix,[115] as well as with other nonperturbative methods, such as the semiclassical path integral description of the density matrix dynamics.[116] In the case of mixed quantum-classical dynamics, recipes have been developed for calculating nonlinear response in the Ehrenfest[94] and Surface Hopping[95,117] approaches.

Can the Exciton Model Always Be the Answer?

In the previous section, we have presented and discussed the possible theoretical approaches and the related computational strategies that are most commonly applied to describe (large) multichromophoric systems. The main assumption at the basis of such methods is that the entire system can be mapped into a simplified Hamiltonian, built from the independent but interacting chromophores. This is indeed a sensible approximation for most of the multichromophoric biosystems; in particular, it has shown to be a very effective approach for describing many LH complexes.[10] These, in fact, are generally rather “rigid” systems, with spatially separated chromophores, whose first excitation is usually well separated from higher-lying states. The approximation of harmonic PESs and linear vibronic coupling, fundamental for the application of most QME techniques, has been verified for chromophores in LH complexes.[35,119,120] Finally, the chromophores in LH complexes generally have long intrinsic lifetimes,[121] thus, the internal dynamics of the excited chromophores can be safely neglected. Of course, there are some exceptions to this general picture of LH systems, especially those containing carotenoids. The latter in fact present a complex ultrafast dynamics with a rapid decay of their bright S2 state.[122,123] Nonetheless, in most LH complexes, the electronic excitations on carotenoids are generally well localized due to the large energy gap with the neighboring chromophores, which has allowed a simplified treatment of photoinduced dynamics for these systems.[67,110] For example, the peridinin–chlorophyll protein, which contains eight polar carotenoids (peridinins) and one chlorophyll per monomer, has been investigated with modified Redfield theory, showing that the exciton model is a powerful technique even for multichromophoric systems containing carotenoids.[124]

However, there is another class of multichromophoric biosystems that represents a real challenge for these methods. This important exception is represented by nucleic acids: here, orbital overlap between neighboring chromophore units is not negligible; moreover, ground and excited states display qualitatively different, possibly crossing, topographies, breaking most of the previously introduced assumptions and activating both photophysical and photochemical processes. A face to face comparison between LH-like and DNA-like multichromophoric systems is shown in Figure .

Figure 5

Different “faces” of the modeling when applied to LH-like and DNA-like multichromophoric systems. On the right, a perfect face (“Mona Lisa” by L. da Vinci, adapted with permission from ref (118). Copyright 2015 American Chemical Society), which represents the modeling of LH complexes, with approximations which can be controlled and tuned. On the left, a peculiar face (cubist-style portrait of Pablo Picasso, image provided by shutterstock.com; image ID 533658418): one identifies various components (eyes, nose, mouth, etc.), but the face as a whole is not easy to read. This represents the present challenges in the modeling of DNA-like multichromophoric systems.

In DNA/RNA systems, the “natural” units, i.e., the nucleobases, already present complex excited-state PESs with multiple crossings between states,[125] which contribute to rich photoinduced dynamics for the isolated chromophores.[126] The energy of UV light, absorbed by DNA due to strongly allowed ππ* transitions, is deposited into the ground state via ultrafast decay pathways, which involve multiple crossings with other singlet and triplet states.[125,127] Exactly because of these unique characteristics, the potential damage (carcinogenic mutations, cell lethality, etc.) induced by UV light absorption is significantly lowered in nucleic acids by highly efficient nonradiative decay to the ground state. On the other hand, these same characteristics prevent a simple application of the techniques described above to model photoinduced dynamics in polynucleotides.

Isolated nucleobases have been extensively investigated by quantum chemical studies, both in gas phase or solution and in DNA-like environments.[125] Both mapping of the PES landscapes of the photochemically relevant states (by critical points, real crossings, and MEPs optimizations)[127−129] and running nonadiabatic MQC dynamics[130−133] or QD[134,135] have been largely used in this context. A full arsenal of QM methods has been also employed, spanning from semiempirical methods[136] to single reference (TDDFT, ADC(2)) and multiconfigurational/multireference–perturbative approaches (CASSCF/CASPT2).

Such a detailed investigation, however, has been shown to be insufficient, as the correct picture of DNA cannot be fully recovered from the properties of the isolated nucleobases. When one moves from single nucleobases to oligo- and polynucleotides, in fact, new and slower (by several orders of magnitude) decay processes appear.[137−140] On the one hand, the already mentioned deactivation pathways from the bright ππ* states are still operative although slightly slowered by the structural constraints of the rigid architecture embedding the nucleobases.[141,142] On the other hand, a wealth of new dynamics arises, originating from “collective” excited states, including Frenkel excitons and CT states, whose existence have been proven both experimentally[138,143,144] and theoretically.[145−147]

To describe this collective behavior, the QM description has been extended to several nucleobases but, as such, the level of the QM approaches has been necessarily lowered. Semiempirical methods such as ZINDO have been used[148] but TDDFT has been mainly applied, including QM regions as large as eight stacked adenines,[146] providing information about the extension of the excitation (mainly over two nucleobases) and showing that the spectral intensity is affected by Frenkel excitons and states with partial CT character in equal amounts. More accurate, but expensive, multireference perturbative approaches have been also employed although mostly limited to QM dimers.[142,147,149,150]

DNA has been also modeled through exciton approaches, showing the sensitivity of the exciton delocalization upon DNA conformations and nucleobases composition.[151−157] Different partitions of DNA structures into sites have also been assessed.[158] Moreover, the extension of the exciton Hamiltonian to CT states, in combination with explicit vibronic dynamics and Redfield theory, has shown that CT states are populated from bright states within a few tens of femtoseconds.[159]

In nucleic acids, however, the application of these exciton models has to be limited to static descriptions or to very short time windows. In fact, at longer time scales, the nuclei move away from the Franck–Condon region and explore both anharmonic regions and crossings with other states. The flexibility of the DNA/RNA scaffold translates into a large variability of the long-range structure of the nucleobases assembly, opening and/or quenching different relaxation paths: indeed, understanding the (static and dynamic) structural factors that control the dissipation of the excess electronic energy in DNA, without producing deleterious reactions, is one of the most fascinating challenges in the field. Photochemical reactions may also take place, implying an additional challenge to the theoretical description of these systems, which must go beyond the one employed in LH complexes; for instance, a proton transfer pathway between base pairs was shown to be key in the relaxation of the excited duplex.[160] Dimerization reactions also occur for specific configurations, showing that the evolution of delocalized Frenkel exciton states may branch into different intra- and interbase decay paths, including CT states.[150,161]

All these specificities clearly make nucleic acids an extremely challenging multichromophoric system where exciton approaches show all their intrinsic limitations.

Challenges and Future Directions

The trivial, “brute-force” approach, which enlarges the size of the QM region so as to consistently describe all the photophysically/chemically relevant processes, is still an unfeasible route for the modeling of large multichromophoric systems: therefore, alternative and more approximate methods have been proposed. Here we have presented and discussed the state of the art of these methods, but we have also shown their intrinsic limitations which still prevent a complete and accurate description of the various photoinduced processes characterizing the different biologically relevant systems. As such, new and alternative strategies are required.

Among them, a promising one is to merge the brute-force QM strategy with the excitonic approach, using larger, still feasible, QM building blocks for the exciton Hamiltonian. As a result, the enormous difficulty of the single full QM treatment is split within several easier QM calculations which can be performed independently and therefore are perfectly suitable for modern multicore architectures. Moreover, by cleverly enlarging the QM units, possible effects of CT states would be automatically included without the need of introducing a posteriori corrections. Finally, by combining this approach with some of the already proposed extensions which remove harmonic PESs constraints and/or integrate the exciton Hamiltonian to TSH nonadiabatic dynamics, a very powerful strategy will be achieved.

To turn this strategy into a feasible computational approach, however, important problems have to be faced. On the one hand, the dimensions of the extended Hamiltonian will grow up significantly, and algorithms will likely need a reformulation to allow exploiting GPU parallelization; on the other hand, the calculation of the forces for the dynamics will call for new smart and fast methodologies such as machine learning methods.[162] Finally, environment models will also need to be improved, both in accuracy and efficiency when applied to problems involving multiple excited states and their nonadiabatic dynamics.

On top of all that, a fundamental problem remains to be faced, namely the correct inclusion of fluctuations. Indeed, in multichromophoric systems, the coupled fluctuations in intra- and intermolecular motions play a much more important role than for a single molecule: in addition to vibronic coupling, fluctuations in the van der Waals and the hydrogen bonding interactions between each chromophore and the components of the bio matrix, in the electrostatic fields and local polarization response, and large scale structural deformations are here extremely important, as they can significantly affect both the nature of the excitons (tuning local excitations, interchromophoric interactions, etc.) and their dynamics.

Taking into account all these effects means to cover many time scales and necessarily requires the use of classical dynamics simulations, possibly in conjunction with enhanced sampling techniques.[163] In this context, a major role will be played by the selected MM force field: if this should not constitute a critical issue for the description of the biomatrix/solvent motions, as FF are optimized exactly for this purpose, the same accuracy is not expected to be retained for the embedded chromophores, where structures coming from classical simulations are likely to introduce artifacts and/or inaccuracies in the evaluation of electronic properties. Even if the highly conjugated nature of the most common chromophores present in multichromophoric biosystems represents a real challenge for MM formulations,[164] recently proposed examples of new FFs optimized for the specificities of each single chromophore have been shown to give a much better description with respect to more traditional FFs.[165−167] On the same line, an effective integration of classical MD based on such optimized FFs and Born–Oppenheimer QM/MM dynamics can represent a valid strategy to obtain accuracy in both the sampling of the complex conformational space of the biosystem and the description of intramolecular motions of the embedded chromophores.[168]

In conclusion, to achieve a complete and accurate description of the complex photoinduced activity in multichromophoric biosystems, a long path has still to be walked. Nonetheless, the smart integration of classical and quantum chemical (static and dynamic) approaches combined with the development in high performing multicore architecture, GPU accelerated algorithms, and machine learning techniques draws an optimistic horizon.