Theory and Simulation of the Ultrafast Double-Bond Isomerization of Biological Chromophores

Ultrafast processes in light-absorbing proteins have been implicated in the primary step in the light-to-energy conversion and the initialization of photoresponsive biological functions. Theory and computations have played an instrumental role in understanding the molecular mechanism of such processes, as they provide a molecular-level insight of structural and electronic changes at ultrafast time scales that often are very difficult or impossible to obtain from experiments alone. Among theoretical strategies, the application of hybrid quantum mechanics and molecular mechanics (QM/MM) models is an important approach that has reached an evident degree of maturity, resulting in several important contributions to the field. This review presents an overview of state-of-the-art computational studies on subnanosecond events in rhodopsins, photoactive yellow proteins, phytochromes, and some other photoresponsive proteins where photoinduced double-bond isomerization occurs. The review also discusses current limitations that need to be solved in future developments.


INTRODUCTION
Nature has developed light-responsive proteins capable of exploiting the energy of a photon to accomplish specific electronic and geometrical changes at the molecular level. Through such changes, the photon energy is stored and then employed to initiate a series of conformational transformations defining the protein photocycle. Such photocycles are implicated in a variety of biological functions such as daylight and nocturnal visual perception in vertebrates and invertebrates, 1 light-oriented growth in plants, 2 phototaxis, 3 light sensing, 4 and ion pumping 5 in eubacteria and Archaea, and ion gating in unicellular algae, 6 to cite a few. In their simplest form, light-responsive proteins consist of a cytoplasmic or membrane protein hosting a prosthetic group, the chromophore (see Figure 1A), responsible for light absorption, photon energy storage, and photocyle initiation. These relatively basic systems are generally classified in families which include rhodopsins, phytochromes, xanthopsins, cryptochromes, phototropins, and BLUF proteins. 9 The latter three families have a flavin chromophore that initiates various intermolecular photophysical and/or photochemical processes (e.g., LOV and CRY in Figure 1A belong to these three families). In contrast, the former three families feature diverse chromophores which, upon illumination, undergo the same "elementary" intramolecular chemical process: a subpicosecond isomerization of one the double bonds of their conjugated frameworks (e.g., RHODOP-SINS and PHY in Figure 1A). In the present review, we mainly focus on these photoisomerizing protein families.
The rhodopsin, phytochrome, and xanthopsins families have served as a virtually endless source of interest, as well as inspiration, for scientists working in different fields such as photochemistry, photobiology, and biophysics. As a consequence, there have been numerous theoretical studies aimed at understanding their spectroscopy, color tuning, and photochemistry, not to mention the light−energy transduction mechanism characterizing their functions. While here, we will be concerned exclusively with double-bond isomerization, we need to warn the reader that the present review is not exhaustive even within this ostensibly narrow topic. Rather, this review reflects the interest of the authors, whose aims are a complete molecular-level understanding of this process and the design of efficient computational tools/strategies providing access to it.
In order to illustrate the impact that the theoretical investigation of double-bond isomerization has on our comprehension of light-responsive proteins, we start with a section (section 2) reviewing the results of complementary spectroscopic and computational studies on the photoreceptor mediating scotopic (i.e., dim light) vision in vertebrates: rod rhodopsin. We start with this example because it is one of the most studied processes in photobiology, both experimentally and theoretically, and therefore serves as a prototypical system. The review will then continue by looking at the current understanding of specific protein families. To do so, it is necessary to start with a section (section 3) revising our present understanding of doublebond isomerization in conjugated organic molecules and the computational methods available to construct and study light responsive protein models. Three more specific sections will then follow focusing on rhodopsins (section 4), photoactive yellow proteins (section 5), and phytochromes (section 6). The involvement of double-bond isomerization in other lightresponsive proteins will be finally discussed in an additional section (section 7). Most of these sections are organized in a way to highlight studies that investigate the structure of CIs first, and then the mapping of excited and ground-state reaction paths, and trajectory studies last. The review is closed by a conclusion section also examining the perspectives of this research field (section 8). We have also included an additional section (section I) that reviews some of the computational methodologies available for studying the thermal and light-induced isomerization of biological chromophores. To get a better appreciation of the work reviewed here, we recommend that the reader reads section I before proceeding to sections 4−7. This will not only serve to get familiar with terms such as "CASPT2", "CASSCF", and "QM/MM", but also to equip the reader with knowledge of the advantages and drawbacks of different approaches, thus allowing them to better understand the results and discussions.

Spectroscopic Studies: Photoisomerization of Rod Rhodopsins Is Vibrationally Coherent
Dim-light vision is initiated by the light-induced isomerization of one of the bonds of the 11-cis-retinal protonated Schiff base (rPSB11). This is the chromophore of rhodopsin: one of the light-responsive proteins found in the retina of vertebrates. The crystallographic structure of bovine rhodopsin is shown in Figure  1B and features seven transmembrane α-helices forming a tight cavity where a lysine residue is covalently linked to the chromophore through a Schiff base group (see Figure 1C). Upon illumination, rPSB11 isomerizes inside the tight protein cavity about the C11C12 double bond, producing the corresponding all-trans isomer (rPSBAT) on a sub-200 fs time scale. 10 This event corresponds to the first step of the rhodopsin photocycle (see Figure 2) and leads to the formation of the photorhodopsin intermediate that, after vibrational relaxation, yields bathorhodopsin. The photon energy stored in the bathorhodopsin intermediate is then used to drive the photocycle and, ultimately, the protein function through the production of the meta II activated state.
A recent spectroscopic study by Miller and collaborators 14 has indicated that the transition from the excited state of rPSB11 to the ground state of rPSBAT is not only ultrafast but also vibrationally coherent. In a contribution 12 on the significance of such a finding, Mathies retraces the progress made through the spectroscopic investigation of double-bond isomerization in bovine rhodopsin (Rh), highlighting Miller's work as the conclusion of a research line started in the early 1990's. In their contributions, both Miller and Mathies provide a mechanistic interpretation of the observed spectroscopic signals in terms of first singlet excited state (S 1 ) decay at a conical intersection (CI) progression along S 1 and ground state (S 0 ) reaction coordinates (and/or paths) intercepting the CI, and trajectories that describe the evolution of excited molecules with time and their transfer from the S 1 to the S 0 potential energy surfaces at the CI region. Such terminology is rooted in theoretical research efforts that have advanced the mechanistic understanding and computer simulation of organic photochemical reactions during the last few decades. Such advancements started to consolidate and became routinely available for computations in the 1990's, therefore roughly at the same time of the early spectroscopic measurements on Rh.
CIs, reaction coordinates, and trajectories constitute key elements for the mechanistic description of photoinduced double-bond isomerizations, as well as for any other photochemical reactions not involving a change of spin symmetry (e.g., singlet photochemical reactions). 15 Indeed, in the Miller 14 and Mathies 12 contributions, the observed coherences in the spectroscopic signals are effectively interpreted as collective nuclear (vibrational) motions of the molecular population toward or away from specific CI structures. Such motion not only ensures that the population moves, as a whole, out of the vertical excitation region of the excited state potential energy surface (see Figure 3A) but also that, in each molecule, the distinct vibrational modes describing such progression maintain a certain phase relationship. More specifically, it is proposed that after light absorption, the population of protein-embedded rPSB11 chromophores initiate a vibrationally phased motion along critical modes such as the bond-length-alternation (BLA) stretching, the hydrogen-out-of-plane (HOOP) wagging, and the torsional deformation of the reactive double bond. These modes would drive the population toward a CI between the S 1 and S 0 potential energy surfaces (in green and red, respectively, in Figure 3A), thus transferring it to S 0 . Such a mechanistic description is consistent with the results of a full body of multiconfigurational quantum chemical computations 16−18 that were initiated at the end of the 1990's. Such computations revealed: (i) the existence of low-lying CIs in rPSB11, (ii) the relationship between electronic origin and structure of the isomerization coordinate of rPSB11, and (iii) the molecular modes dominating the rPSB11 excited-state progression. In the present review article, we will highlight computational work Bathorodopsin may be regarded as the intermediate storing the photon energy. This energy is then used for the production of meta II: a key event within the photocycle constituting the signal-transducing state of rhodopsin (i.e., its biological function). Since the C11C12 isomerization is fully accomplished at the level of the photorhodopsin/ bathorhodopsin intermediates, the present contribution reviews the mechanistic studies focusing on the framed part of the photocycle. focusing on the structure of CIs, the mapping of reaction paths, and trajectory simulations. In fact, each section of this review is organized such that it mirrors this section.

Computational Studies: Rh Isomerizes via a Two-State, Three-Mode Mechanism
It is instructive to discuss at this point how computational studies connect with the experimental results (i−iii) introduced above and, at the same time, provide a motivation for the notions that will be examined in the following six sections. On the other hand, the reader may also skip the present subsection without loss of continuity and revisit it at a later stage after reading sections 3 and 4, if they feel they would prefer a more general introduction to the theoretical terms and concepts used in this subsection.
The impact of the existence of low-lying CIs in rPSB11 is apparent when comparing Figure 3 (panels A and B), illustrating the mechanism proposed in Miller and Mathies' contributions 12,14 with that given in an early report by Schoenlein et al., 10 respectively. In contrast to the early view proposing decay at an avoided crossing (AC), the population transfer from S 1 to S 0 is presently assigned to the efficient internal conversion typical of CIs. Building on the seminal work of Michl and BonacǐcḰ outecky, 19 a CI mediating the isomerization of a reduced rPSB11 model has been computed in 1997 20 shortly after CIs were shown to be implicated in common organic photochemical reactions. In the same report, it was shown that a barrierless coordinate connects the vertical excitation region to the CI located at the bottom of the excited-state potential energy surface. 15 Also, the electronic structure of the excited state driving isomerization was assigned to a B u -like (charge-transfer) rather than A g -like (diradical) wave function symmetry (see also section 3.1): a result that has been used by Miller and collaborators 14 in their analysis of the coherences assigned to the excited-state HOOP mode describing the rPSB11 relaxation.
The second result, the relationship between electronic origin and structure of the isomerization coordinate of rPSB11, is a consequence of the B u -like S 1 electronic structure. Promotion to a B u -like state leads to an instantaneous charge-transfer process effectively unlocking the double bonds and locking the single bonds of rPSB11. This induces a stretching relaxation along C11C12 and other conjugated bonds contributing to a BLA mode (this is quantified by the difference between the average of the double-bond lengths and the average of the single-bond lengths of the conjugated chain). The resulting C9−C10 C11−C12C13−C14 pattern (compare the ground-state and excited-state double-bond pattern in Figure 4A) allows for a facile torsion about the reactive C11C12 bond. In Rh, this torsion is accompanied by torsional motion along the adjacent C9C10 bond. 20,21 Such a mechanism, termed the bicyclepedal coordinate, was first proposed by Warshel 22 and confirmed for rhodopsins from different organisms 23,24 after the development of CASSCF-based quantum mechanics/molecular mechanics models (i.e., full protein models). 25 While such a confined coordinate facilitates isomerization in the tight protein cavity, it has been shown that isolated rPSB11 26 and quantum mechanical/molecular mechanical (QM/MM) models 21 feature essentially the same bicycle-pedal coordinate that, however, is aborted after decay to the ground state (i.e., the C9C10 twisting reverts to the original configuration yielding isomerization at C11C12 exclusively). A more detailed description of these results will be given in section 4.3.
The third finding, the modes dominating the rPSB11 excitedstate progression, is related to the second and refers to an early report by Gonzaĺez-Luque et al., 16 where rPSB11 models were used to investigate the isomerization coordinate. It was shown that the aforementioned relaxation along the BLA mode (incorporating the reactive C11C12 stretch of Figure 4A) precedes the reactive C11C12 torsion. The result pointed to a coherent two-state two-mode motion of the Rh population. Computational evidence for the participation of the HOOP mode in the coherent excited-state relaxation came later when semiclassical trajectory computations became possible. 27   Quantum chemical calculations have, as discussed above, provided a framework for the mechanistic/theoretical understanding of the photochemistry of dim-light vision. However, the reliability of computational studies depends critically on the quality of the computer model used to represent the protein and the light-induced processes in the chromophore. Such models are complex and, currently, not univocally defined. For this reason, they require a validation based on the simulation of observed quantities (e.g., absorption maxima and excited-state lifetimes) or the prediction of quantities that could be experimentally verified. For instance, one of the Rh models 27 introduced above predicted a HOOP mode with a ∼40 fs period (∼830 cm −1 ) and a BLA mode with a ∼20 fs period (∼1670 cm −1 ) and thus qualitatively consistent with the values recently reported by Miller and collaborators for the HOOP (746 cm −1 ) and C11C12 stretch (1679 cm −1 ). Such comparisons support the initial relaxation mechanism provided by the model which points to a C11C12 torsional progression reaching a 40°value at 50 fs.
Validated models are also needed to open new alleys for the comprehension of key properties of double-bond photoisomerization such as excited-state lifetime and quantum yield.
For instance, an effect still under computational investigation is associated with the intersection space structure of rPSB11 35 (IS, see Figure 4B and associated legend). The IS is a (N-2)dimensional nuclear coordinate subspace (N being the vibrational degrees of freedom of the system) formed by CI points. The CI of Figure 3A and 4B is therefore only one point along the IS. Multiconfigurational quantum chemical calculations show that the IS may be approached not only via large, ca. 90°C11 C12 torsions but also via progression along the BLA and HOOP modes. 27 Thus, it could be possible to achieve a CI region at smaller value of the C11C12 torsion when the progression along the HOOP mode approaches the IS (this is schematically represented by the dashed curve in Figure 4B). See section 3.4 for further details.

THEORY OF DOUBLE-BOND ISOMERIZATION
The Rh case study reviewed above shows that a theoretical framework is needed for interpreting or planning experiments designed to reveal the mechanism of double-bond isomerization in light-responsive proteins. Accordingly, in the present section we review the basic theortetical and mechanistic concepts for double-bond isomerization that have emerged from quantum chemical studies. Note that a discussion of the methodologies that may be employed in these studies is available in section I. Here, we continue to focus primarily on protonated Schiff bases as paradigmatic examples of isomerizing chromophores.
To understand how nature has utilized the simple doublebond isomerization motif to drive, very efficiently and effectively, different photobiological processes, two important questions come to mind. First, why has nature selected certain organic chromophores in photobiological systems rather than others? While the answer to this question depends on different factors, including the availability of biosynthetic pathways leading to specific conjugated frameworks, in this section we provide a partial answer by focusing on the factors that determine light sensitivity at the molecular level. Second, how does nature optimize the properties of the chromophore (e.g., to maximize sensitivity to light) by tuning its environment? Indeed, the chromophore environment (i.e., the amino acid composition of the protein cavity surrounding the chromophore) can modulate its maximum absorption wavelength (λ max ), its excited-state lifetime (τ cis→trans ), its photoisomerization "quantum yield" (QY), and the barrier (E a T ) controlling the rate of the competing thermal isomerization.
The events included in the framed region of the photocycle of Figure 2 are associated with specific features of the S 0 and S 1 potential energy surfaces of the chromophore. In the simplified Scheme 1, these S 0 and S 1 N-dimensional potential energy surfaces (with N being the number of independent nuclear degrees of freedom or independent vibrational modes) are reduced to one-dimensional curves along a reaction coordinate. This scheme can be used to assign computable quantities and their relation to experimental observables. Accordingly, light absorption is characterized by the vertical excitation energy ΔE(S 1 −S 0 ) approximating the photon absorption wavelength λ max of Figure 3B. The absorbed photon promotes the chromophore from its S 0 equilibrium structure to the Franck− Condon (FC) point of its S 1 potential energy surface. At this stage, to enhance light sensitivity of the system, it is important to increase the efficiency of the photochemical process. For example, the isomerization QY must be maximized in order to increase the ratio of reactive trajectories relative to nonreactive ones that do not trigger the photocycle. At the same time, any  21 and 27) and representing the center of a coherent population (i.e., the vibrational wavepacket). The trajectory is dominated by progression along three modes: BLA (difference between average single and double-bond lengths of the rPSB11 backbone in Å), HOOP (in degrees), and the bicycle-pedal mode incorporating the C11C12 torsion (in degrees). The initial 5 fs time segment is dominated by BLA, but the HOOP progression is active on a 15 fs time scale. The dashed curve provides a qualitative representation of the projection of the intersection space (indicated with the IS label) of the torsion-HOOP subspace. As we will mention in section 3.2 and discuss in the appendices, the element of the intersection space are geometrically and energetically distinct CI points. competing processes must be minimized. This includes thermal isomerization, which could compete with the photoisomerization process or, in the case of visual rhodopsins, creates the "background noise" in the signal and therefore lowers visual sensitivity. Therefore, an understanding of how to improve light sensitivity requires both calculation of QY and rate of competing thermal isomerizations. To compute QY, one may start by considering the correlation proposed by Weiss and Warshel 36 and Mathies and co-workers. 37 The idea is that, consistently with the Landau−Zener model, 38 a larger velocity of the S 1 population moving toward the CI results in a higher reaction QY. This indicates that variation in QY correlates with the change in S 1 lifetime estimated, for instance, by computing the τ cis→trans time required to reach the CI or, more generally, the point of hop from S 1 to S 0 (i.e., the point when the trajectory algorithm decides to switch or "hop" from one state to another). This simplifies QY calculations considerably because one can get a sense of differences in QY in different systems by simulating the excited state lifetime of a single trajectory. Instead, as illustrated in the inset of Scheme 1, a QY value can be obtained in a more bruteforce way by exploiting the S 0 and S 1 potential energy gradients and nonadiabatic couplings to compute a set (sometimes called ensemble or swarm) of nonadiabatic trajectories. Such calculations start with a simulation of the S 0 conformational population of the protein in the dark (i.e., ideally a representation of the Boltzmann distribution). This provides a finite number (N t ) of S 1 initial structures and velocities whose time evolution is simulated. The QY value can then be computed as the ratio between the number of trajectories successfully producing the isomerized chromophore (N p ) and N t . Such explicit QY computations have been reported 32,39,40 for Rh. Finally, the variation in thermal isomerization rate can be computed from the potential energy barrier E a T corresponding to the energy difference between the transition state (TS) and the protein equilibrium structure. In conclusion, the computational investigation of the photochemical and thermal isomerization of a double bond requires, at the very least, the calculation of the S 0 and S 1 state potential energy surfaces of the chromophore. As discussed below, the S 2 potential energy surface may also be implicated, so this also needs to be computed for some cases.

Vertical Excitation Energies and Conical Intersections
The selection of a suitable conjugated system hosting the reactive double bond may be discussed on the basis of the comparison of linear conjugate hydrocarbons and protonated Schiff bases. Here we focus on systems with five conjugated double bonds. Accordingly, we compare deca-1, 3,5,7,9- The first three electronic singlet states (S 0 , S 1 , and S 2 ) of the reference hydrocarbon decapentaene in its S 0 equilibrium C 2h structure can be classified as having 1A g , 1B u and 2A g symmetry. The 1A g character corresponds to S 0 and has a closed-shell covalent electronic structure. The 1B u state corresponds to S 1 , and is dominated by a singly excited electronic configuration leading to a charge-transfer character, while the 2A g state, corresponding to S 2 , is dominated by double excitations and has a diradical character. The selection rules for electronic transitions dictate that only the 1A g to 1B u and 1B u to 2A g transitions have nonzero oscillator strength.
In the discussion below, the electronic character of the singlet electronic states of our model chromophores will still be classified using the C 2h symmetry labels 1A g , 1B u , and 2A g [i.e., related to a polyene hydrocarbon such as tEt-hexa-1,3,5-triene in its planar S 0 equilibrium geometry (italics stereochemical labels indicate the configuration of single (trans or cis) and double bonds (E or Z), respectively], even for structures with C s symmetry like PSB5 or no group symmetry like any structure located along the PSB5 isomerization reaction coordinate. Alternatively, the same states will be labeled Φ COV/DIR , Φ CT , and Φ XDIR , respectively, to explicitly indicate the character of the dominating electronic configuration (or diabatic state). Notice that along an isomerization path or a trajectory, the character of the potential energy surfaces associated with S 0 , S 1 , and S 2 (i.e., classified according to their energy order which will always be, by definition, S 0 < S 1 < S 2 ) may change. Accordingly, one may have an S 0 potential energy surface with a dominating Φ COV/DIR (i.e., here associated with a 1A g ) character in certain region and a Φ CT (i.e., 1B u ) character in some other region, or a S 1 potential energy surface with a dominating Φ CT character in a region and a Φ COV/DIR or Φ XDIR (i.e., 2A g ) character in a different one.
Note that later on, in section 3.2, we will further distinguish Φ COV/DIR as being either Φ COV or Φ DIR , depending on the geometry, but it is important to recognize that these two labels both describe the same electronic state. To elaborate, Φ COV/DIR indicates that along a double-bond twisting path, the Φ COV closed-shell electronic configuration dominating the S 0 state of the reactant (or product), smoothly changes into Φ DIR openshell configurations representing an orthogonal diradical. This is equivalent to adopting a valence-bond view of a homolytic bondbreaking process where two spin-paired electrons become Scheme 1. Schematic Representation of Photochemical (Full Arrows) and Thermal (Dotted Arrows) Isomerization Paths a a The two paths develop along different reaction coordinates. In fact, the CI is located energetically above the transition state controlling the thermal isomerization (the thermal transition state is indicated with the label TS and corresponds to a first-order saddle point located along the S 0 potential energy surface), features a different geometrical structure, and drives a far-from-equilibrium process. ΔE(S 1 −S 0 ), τ cis→trans , and E a T (in red) are computable quantities. As shown at the top-right, quantum yields (QYs) can also be computed as the ratio of non-adiabatic trajectories reaching the isomerized chromophore (N p ) over the total number of computed trajectories launched (N t ). The competing photoproduct photoisomerization (dashed arrows) yielding the initial reactant may also be computationally investigated as will be discussed in section 4.2.1. localized on distinct fragment π-orbitals, and the bond breaking is described by a change in overlap between these singly occupied orbitals. Therefore, Φ COV and Φ DIR are seen as different components of the S 0 electronic structure. Notice that Φ XDIR is distinct from Φ DIR . In fact, Φ XDIR represents a diradical state obtained via electronic excitation at the untwisted reactant geometry and, as we will see below, may assume a tetraradical character upon double-bond twisting. When using the same valence bond view, Φ CT corresponds to an electron transfer from one fragment π-orbital to the other.
The energies of low-lying CIs are shown in Figure 5A. The diagram indicates that excitation in decapentaene occurs to the S 2 state, after which a 2A g /1B u CI between the S 2 and S 1 potential energy surfaces may be readily accessed from the FC point after photon absorption. After decay to S 1 , the evolution along the S 1 potential energy surface may then lead to a 2A g /1A g CI, where the decay to S 0 would occur. A qualitative analysis of the geometrical and electronic structure of these CIs indicates that (i) the S 0 photoproduct may be reached with UV rather than visible light excitation, (ii) the E to Z isomerization at C5C6 may occur via the 2A g /1B u CI, potentially generating an isomerized photoproduct on the S 1 potential energy surface, and (iii) the decay to S 0 occurs via a 2A g /1A g CI featuring a −(CH) 3 − kink with a tetraradical character. 44 The "kink" geometrical and electronic structures support the possibility to simultaneously generate two different primary photoproducts (a cyclopropyl diradical intermediate and a Z-isomer) in addition to the possibility of regenerating the initial E-isomer. Points (i−iii) explain why linear polyenes would not be favorable chromophores for absorbing photons of visible light selectively generating the isomerization product without going to secondary (e.g., cyclization) photoproducts. In other words, a protein cavity hosting a decapentaene chromophore would have to severely modify, through environmental effects, the three potential energy surfaces to achieve favorable spectroscopic and reactivity properties. Furthermore, in a short polyene hydrocarbon, the charge separation occurring upon vertical excitation (i.e., via the 1A g to 1B u transition) is spread along the chain and does not have a direction, making it more difficult to control the chromophore properties via local electrostatic effects (i.e., via interaction with specific polar residues located in the protein chromophore cavity).
The results discussed above indicate that, in decapentaene, there is no direct access to a low-lying 1B u /1A g CI where decay to S 0 may occur directly. However, an accessible 1B u /1A g CI has been shown to exist in arylethenes such as stilbene (see inset in Figure 5A). With respect to the FC point, such a CI is found 45 lower than the decapentane 1B u /2A g CI and is related to an observed S 1 intermediate called the "P-state". This intermediate features a highly twisted double bond with a pyramidalized sp 3 carbon hosting an anionic center and a positive charge partially delocalized on a benzyl moiety. Thus, E-stilbene decays to S 0 from the initially populated S 1 state, leading to its Z-isomer after having accessed the P-state and the nearby 1B u /1A g CI. However, this decay mechanism does not occur for decapentaene.
We will now compare the energies of low-lying PSB5 CIs to the one described above for decapentaene. Figure 5B shows a 1B u /1A g CI that could be readily accessed from the FC point. In other words, in the absence of an S 1 energy barrier, the events initiated by the photon absorption leads directly to S 0 without involving a 2A g /1B u CI, which is now located at a significantly higher energy. Most importantly, due to the positive charge residing on the chromophore, the 1B u /1A g CI features a single 90°twisted C5C6 double bond with no other major skeletal distortions. This CI is thus geometrically different from the stilbene 1B u /1A g CI characterized by an S 1 state with charge separation rather than by a charge transfer. In fact, the S 1 state near the PSB5 CI may only feature a neutral lone pair localized on the nitrogen atom of the Schiff base functional group with a positive charge far away from it (see Scheme 2). A qualitative analysis of the geometrical and electronic structure of the PSB5 CI indicates that (i) the CI may be reached with visible light Figure 5. Comparison of relative energies of FC and CI energies in the isoelectronic decapentaene and PSB5. The energies computed at the 3root-SA-CASPT2(IPEA = 0)//2-root-SA-CASSCF/6-31G* level of theory (the method acronyms are introduced in section I.1. Note that the double slash notation indicates that energies were computed at the level of theory at the left, here 3-root-SA-CASPT2, and geometries at the second level of theory at the right, here 2-root-SA-CASSCF). (A) Relative stability of the FC states and low-lying CIs of decapentaene. The schematic representation of the geometrical and mixed S 1 /S 0 electronic structure of the CIs is also given. The dashed lines connect FC and CI structures for visibility. The arrows represent relaxation paths from the FC points. The framed structure represents the CI driving decay from an intermediate (P-state) in stilbene and is shown for comparison. The dashed segment connecting the two radical centers indicates that these electrons are singlet coupled. (B) Relative stability of the FC states and low-lying CIs of PSB5. Even in this case, the dashed segments connecting two radical centers indicate that these are coupled (i.e., they appear uncoupled only because they reside in two nearly orthogonal π-systems). In all cases, the Lewis formula used for representing the electronic structure of the CI points represents mixed electronic structures. This is illustrated in Scheme 2 for the case of PSB5. Data from ref 43. excitation and (ii) the isomerization of the central C5C6 double bond occurs via a fully (90°) 1B u /1A g CI potentially generating exclusively an S 0 isomerized photoproduct. Point (i) and (ii) appear to be favorable features for a biological chromophore.
A protein cavity hosting a "PSB5-like" chromophore such as rPSB11, would easily modulate the ΔE(S 1 −S 0 ) and therefore the λ max . In fact, in contrast with decapentaene, the charge transfer occurring upon S 0 to S 1 vertical excitation would result in a positive charge translocation (delocalization) from the NC2 bond to the β-ionone ring. Such a charge translocation is a general feature of protonated and alkylated Schiff bases and is documented in Figure 6A for rPSB11. The molecular environment (i.e., the protein and solvent cavities) interacts with the charge distributions of the S 0 and S 1 states differently and, ultimately, modulates the ΔE(S 1 −S 0 ) value and thus the corresponding λ max . The mechanism through which proteins influence the absorption wavelength λ max of the chromophore is known as color tuning or spectral tuning, a topic that has been widely studied computationally 50−53 and experimentally. 54,55 Other ΔE(S 1 −S 0 ) and ΔE(S 2 −S 1 ) modulating effects are achieved by changing the effective length of the conjugated chain via single bond twisting (e.g., the C6−C7 bond in rPSB11). The magnitude of such an effect has been calculated, 23 and the corresponding values are displayed at the bottom of Figure 6B. An additional effect that also changes ΔE(S 1 −S 0 ) and ΔE(S 2 − S 1 ) in the same direction as increasing the conjugated chain length is by modulating the BLA geometry. This effect can be inferred comparing Figure 6B top with 6B bottom. The top panel shows energies for the S 1 equilibrium structure of a constrained planar (i.e., untwisted) conjugated chain. At a planar geometry, the S 1 electronic structure of PSB5 is dominated by a Φ CT electronic configuration and, therefore, its equilibrium geometry features an almost inverted BLA pattern. Similar modulating, but more complex, effects are expected when looking at the energies of the 1B u /1A g and 2A g /1B u CIs, suggesting that the reaction path leading from the FC point to the CIs may also be modulated by the environment surrounding the chromophore.
Above, we have exclusively focused on the 2A g /1B u CI associated with the twisting about the central C5C6 double bond of PSB5. However, similar CIs also exist that are associated with the twisting about other double bonds in the chromophore. As shown in Figure 7 for the 2A g /1B u CI mediating the C7C8 double-bond isomerization, such CIs have different relative energies. Therefore, in general, a mixture of isomerized Scheme 2. Lewis Resonance Formulas Representing the Electronic Configurations Degenerate at a 1A g /1B u CI in PSB5 Featuring a ca. 90°Twisted C5C6 Double Bond a a At the FC point, the S 0 and S 1 adiabatic electronic states are dominated by the 1A g and 1B u configurations, respectively. Note that 1A g corresponds to Φ COV at the FC point but to Φ DIR at the twisted CI geometry. However, in the CI vicinity, the S 0 and S 1 states display a mixed and interchanging electronic character. The crossing of 1A g and 1B u electronic characters at the CI is represented by the resonance formula at the right.   photoproducts (different stereoisomer of PSB5) may be expected, which will depend not only on their relative energies but also on the details of the reaction coordinates leading to the corresponding CI. It can be hypothesized that the chromophore environment (i.e., again, the protein cavity or the solvent cavity and counterion) controls the isomerization selectivity by stabilizing a specific CI and/or favoring a specific reaction coordinate. Indeed, as we will review in section 3, the cavities of retinal proteins favor single isomerizations about either the C13C14 or the C11C12 double bond of their retinal chromophores. However, secondary isomerization products (e.g., involving the C9C10 and C7C8 double bonds) have been reported.

Structure of the Photochemical and Thermal Reaction Paths
A deeper understanding of the possible effects of the protein cavity on the nonadiabatic S 1 /S 0 transition driving a double-bond photoisomerization calls for the geometrical and electronic characterization of the branching plane (BP) of the CI. 15,58−60 This is the plane (also called the g-h plane) 60 defined by the molecular modes corresponding to the gradient difference and derivative coupling vectors computed at the CI. These vectors are introduced in detail in section I.1. A pictorial representation of the BP and intersecting potential energy surfaces associated with the S 0 and S 1 states of the rPSB11 chromophore is given in Figure 8. One of the vectors is dominated by a combination of twisting deformation about the reactive C11C12 double bond and a HOOP wag. The other vector corresponds to a BLA stretch. 57 An analysis of the S 0 and S 1 charge distributions along the BP reveals a relationship between the electronic structure of the S 0 and S 1 states in the CI region and in the FC region. The main point is that, consistently with a quantum mechanical theorem called the geometric phase, 61 the S 1 state changes its electronic character along a small loop centered at the CI and residing on the BP. More specifically, in rPSB11, its character passes two times (e.g., roughly every 180°) from a Φ CT character typical of the S 1 state at the FC point to the Φ DIR character associated with the S 0 state at the FC point. In Figure 8, these different regions of the S 1 potential energy surface are colored in brown and green, respectively. The geometric phase theorem explains why an opposite electronic state change is observed for S 0 , namely when one moves along the same loop. Therefore, on the S 0 potential energy surface, there are regions dominated by a Φ CT character rather than the Φ DIR (or Φ COV ) character expected considering the S 0 state in the FC region. The qualitative structure of the S 0 and S 1 potential energy surfaces driving the rPSB11 chromophore evolution is given in Figure 9A. In order to drive an ultrafast chemical reaction, it is not enough to have a CI with a suitable geometrical and electronic structure and BP. The S 1 potential energy surface must also feature a barrierless (or nearly barrierless) reaction path connecting the FC point to the CI. The properties of such a path (i.e., reaction coordinate, slope, and electronic character) may be controlled by the interactions with the protein environment. In principle, such interactions may slow down or even stop the isomerization if an energy barrier is created between the FC point and the CI. However, this is not the case for the S 1 potential energy surface illustrated in Figure 9A, which shows a barrierless potential energy surface consistent with that of Rh. 21 In protonated Schiff bases that are free to isomerize, there are specific molecular structure changes describing the path toward the CI (e.g., BLA, HOOP, and C11C12). 21,27 In the following sections, we will see that these changes also describe the excited state relaxation of other biological chromophores whose S 1 state is dominated by a charge-transfer excitation Φ CT of their πsystems. Such a path is described by progression along multiple coordinates (or modes) and cannot be effectively visualized in a 2D representation like the one of Figure 9A. In fact, already a three-dimensional (3D) representation would reveal that a CI is actually not an isolated point, but it is part of a N-2-dimensional "line" of CIs called intersection space or IS (see also Figure 4). While here, we do not further discuss the relationship between the S 1 isomerization coordinate and the IS. In section II, we show that it has been possible to map the S 1 paths connecting the FC structures of different isomers of a minimal rPSB model to the corresponding CIs as well as the position of these CIs within an IS segment. 62 This information is useful for discussing the dynamics of the S 1 deactivation, which does not occur by following a single "static" path but, rather, by following different trajectories possibly reaching different IS points.
In Figure 9B, we display a schematic representation of the S 0 potential energy surface of Rh. As we will discuss in section 3, the features of this surface determine if the thermal isomerization of the Rh chromophore may or may not compete with its photoisomerization by generating "thermal noise." Indeed, consistently with Scheme 1, such competition depends on the activation barrier E a T associated with an S 0 transition state. However, it has been shown that, as a consequence of the potential energy surface topography imposed by the Rh CI, the thermal isomerization may occur via two geometrically and electronically distinct transitions states: TS CT , characterized by the Φ CT configuration, and TS DIR , characterized by the Φ DIR electronic configuration. 56 Note that dominating transition state structure (i.e., the one associated with the lower E a T value) would be also responsible for the thermal back-isomerization of the chromophore to the original reactant. . Schematic representation of the S 0 and S 1 potential energy surfaces plotted along the BP of the CI. The structures and BP vectors correspond to the CI of the PSB11 chromophore of bovine rhodopsin. A loop located on the BP and centered at the CI is shown. A projection of the loop on the S 0 potential energy surface goes from a Φ DIR region (green) to a Φ CT region (brown) and then returns to a Φ DIR region. The Newman projection on the top-right defines the HOOP internal coordinate as corresponding to the dihedral angle representing the displacement from the "ideal" equilibrium position (dashed segments). Adapted with permission from ref 56

Time Evolution of the Excited-State Population
Above, we have focused on the structure of the S 1 potential energy surfaces controlling the chromophore photoisomerization. However, a barrierless reaction (i.e., not controlled by a transition state) can only be satisfactorily described by looking at the evolution of the geometrical and electronic structure of the system as a function of time. In Figure 10, we report the results of a single nonadiabatic trajectory calculation providing information on such process for Rh. 21 Figure 10A shows the time evolution of the S 1 and S 0 energy profiles and the S 1 −S 0 gap. In ref 21, the authors also report the oscillator strength evolution (proportional to the transition charge-transfer character and intensity at the corresponding wavelength). The same calculation provides information on the evolution of the associated geometrical parameters and on the following reaction coordinate, which has already been introduced in Figure 4. In Figure 10B, we display the evolution of the twisting deformation about the reactive C11C12 double bond, unreactive C9C10 double bond, and HOOP parameters. In this way, it has been possible to document the "aborted" bicycle pedal reaction already mentioned in the Introduction. The same calculation also provides information on the evolution of the charge distribution and, therefore, Φ CT or Φ DIR character of the chromophore electronic structure during the S 1 motion and after the hop to S 0 . Notice that such a hop does not necessarily occur at exactly the CI point but, rather, as mentioned in the previous section and detailed in section I.3, the hop occurs along the "volume" surrounding the CIs belonging to an IS "segment".
In Figure 10C, it can be seen how the initial Φ CT character evolves according to the process schematically anticipated in Figure 9. Schematic representation of the structure of the potential energy surfaces driving the ultrafast double-bond photoisomerization and thermal isomerization of rPSB11 in Rh. (A) S 1 potential energy surface driving the photochemical isomerization. The changes in the charge distribution of the chromophore during S 1 evolution and after relaxation to S 0 are also shown with bubbles. The bubbles represent positive (blue) and negative (red) charges on the atoms of the chromophore framework. The curly arrows indicate the intramolecular electron transfer process driving the charge translocation occurring upon photoexcitation and also during the progression and decay of the excited state. The top structure corresponds to the S 1 structure of the chromophore in the vicinity of the CI (the framed region corresponds to the "double cone" seen in Figure 8). The associated electron transfer is the one occurring during the hop from S 1 to S 0 . (B) Mechanism of the S 0 thermal isomerization. This transformation is controlled by two transition states (TS CT and TS DIR ) located along the BLA coordinate and on opposite sides with respect to the CI. TS CT has a charge-transfer character (Φ CT ) and therefore a 1B u -like electronic structure. In contrast, TS DIR has the diradical character (Φ DIR ) of a 1A g -like electronic structure that becomes a closed-shell covalent electronic structure at the reactant (Rh) and product (bathoRh) geometry. The bubbles represent positive and negative charges on the atoms of the chromophore framework at the corresponding transition states. The cross-section along BLA (computed at the 2-root-SA-CASPT2//CASSCF/6-31G* level of theory) provides information on the energy difference between TS CT and TS DIR and shows that in Rh the transition state controlling the isomerization is TS CT . Adapted with permission from ref56.  Notice that at ca. 100 fs, the S 1 state crosses the S 0 state. At this point, the system hops from the S 1 to the S 0 state and continues its motion along the S 0 potential energy surface. (B) Geometrical changes along the S 1 branch (i.e., until the hop point) of the same nonadiabatic trajectory. The plotted geometrical parameters contribute to one of the BP vectors (C10−C11−C12−C13 twisting in black □ and HOOP in blue ○). We also report the coupled C8−C9−C10−C11 twisting (dashed curve), forming the S 1 bicycle-pedal reaction coordinate 22 discussed in the Introduction. Notice that immediately after the hop, the C8−C9−C10− C11 dihedral twisting inverts direction (not shown) and the bicyclepedal motion is ultimately aborted. 21 (C) Changes in the charge distribution along the same trajectory. The parts of the trajectory dominated by a Φ CT electronic structure are marked in brown, while the one dominated by a Φ DIR electronic structure are marked in green. The QM/MM model used is based on a computer time-saving CASSCF-(12,12)/6-31G*/AMBER level of theory, with an active space comprising the full π-system and with methyl groups and part of the β-ionone ring in the MM subsystem to reduce the size of the QM subsystem.  Figure 9A. 21 Notice that a wider region bearing a Φ CT character in S 0 implies a wider Φ DIR dominated region on S 1 . This wider S 1 region would imply a charge evolution significantly different from the one of Figure 10C. In particular, after an initial charge transfer from the Schiff base to the hydrocarbon region, the charge partially gets back to its original position and the fully twisted chromophore would thus have a diradical character.
The trajectory of Figure 10 only describes the motion of the molecule after a vertical transition from its S 0 equilibrium structure to S 1 and with no initial velocities. Accordingly, such a trajectory only provides qualitative information on the dynamics of the photoexcited population (also called vibrational wavepacket), as one misses the statistical description of the motion prepared, for instance, by a laser pulse. Such a population includes a large number of molecules differing in structural displacement and velocities that are consistent with the Boltzmann distribution at the experimental temperature. It is thus the dynamics of the population that one needs to simulate to provide parameters that can be directly compared with the experiment and, most importantly, to determine QY values. As schematically illustrated in Scheme 1, the vibrational wavepacket time evolution can be simulated by computing a finite set of nonadiabatic trajectories each associated with a structure and velocity extracted from a simulated Boltzmann (or Wigner 63,64 ) distribution. Such calculations have actually been reported for Rh by various groups. 32,34 The simulation of the population dynamics with a suitable number of trajectories can be used to simulate transient spectroscopic quantities and compare them with the results of pump−probe transient spectroscopy experiments. For instance, the simulation of transient emission and absorption spectra is schematically illustrated in the inset of Figure 11A, which primarily reports the energy profile computed along five different nonadiabatic trajectories launched from structures and velocities selected from some initial Boltzmann distribution. 32 More specifically, it is possible to trace, along entire sets of trajectories, the evolution of different energy gaps and oscillator strengths between different states, as schematically illustrated in the inset.
In a more recent study by Polli et al., 40 the photoisomerization quantum yield of isorhodopsin (isoRh), a secondary photoproduct of the photoisomerization of bovine rhodopsin featuring a rPSB9 chromophore (see Figure 12A), was simulated using a larger set of semiclassical trajectories. The results provided an overview of different competing processes occurring during the relaxation of the initial isoRh S 1 population. As displayed in Figure 12A, part of the population simply returns to a vibrationally excited form of the original S 0 reactant. Part of the population undergoes a single isomerization event at either the C11C12 double bond, generating a new intermediate with 11-cis and 9-cis double bonds (termed rPSB911) or at the C9 C10 double bond to give bathorhodopsin (containing a rPSBAT chromophore). Finally, a fraction of the population generates Rh (containing a rPSB11 chromophore) via a simultaneous isomerization of the C11C12 and C9C10 double bonds.
In principle, the analysis of the resulting trajectory ensemble not only provides a test for the validity of the mechanistic information extracted by minimum energy path (MEP) computations and single-trajectory probing discussed above but also can be used to directly simulate the evolution of the signals obtained via transient absorption spectroscopy and compare with the experiment. For instance, as shown in Figure  11B, the simulation of time-resolved absorption spectra has been reported for Rh. The matching between experimental and simulated spectra 32 provides a validation of the mechanistic hypothesis originated with previous investigations. In particular, it confirms the involvement of a barrierless relaxation path and decay in the vicinity of a CI (supported by the progressively redshifting wavelengths associated with the observed stimulated emission signal) during the internal rPSB11 photoisomerization process. Such studies employing an ensemble of trajectories therefore can be used to compute more reliable excited-state lifetimes (τ cis→trans ) compared to single-trajectory calculations.
The calculation of hundreds, if not thousands, of nonadiabatic trajectories allows QY computations in the way anticipated in Scheme 1. These are calculated as the N p /N t ratio between the number N p of trajectories reaching the S 0 potential energy surface valley centered on the photoproduct energy minimum (e.g., for Rh, this would be the bathoRh intermediate) and the total number of nonadiabatic trajectories hopping from S 1 to S 0 . However, presently, the accurate computation of QY is limited by the computational cost associated with running many trajectories at an accurate level of theory. Indeed, only indicative values have been reported on the basis of less than 100 trajectories. 32,40 For instance, the QY of rPSB11 isomerization in Rh was estimated using 40 trajectories and found to be 0.61 compared to a 0.64 experimental value. 32

Quantum Yields and Isomerization Stereochemistry
The example above indicates that nonadiabatic trajectories can be employed to investigate the selectivity (see Figure 7 and related text) of a double bond photoisomerization. In fact, selectivity is quantitatively described by QYs of isomerizations controlled by distinct and competing CI channels. For isoRh, a representation of the structures of these CIs (CI 1 and CI 2 ) are given in Figure 12B. CI 1 leads to isoRh reconstitution with a 0.6 QY, indicating that this is the major channel in isoRh photochemistry. However, the same channel also leads to a 0.2 QY of the derivative featuring the PSB911 chromophore and produced by following an aborted bicycle-pedal reaction coordinate ultimately leading to C11C12 isomerization exclusively. In contrast, the CI 2 channel is leading to a very limited (0.04 QY) production of Rh produced via two parallel double-bond isomerizations. This same channel mainly produces bathoRh (featuring the rPSBAT chromophore) through a single C11C12 isomerization with ca. 0.15 QY. Such isomerization is driven by progression along the same aborted bicycle-pedal reaction coordinate but followed in the opposite direction.
As mentioned above, double-bond isomerization selectivity is determined by differences in QYs obtained, as anticipated in Scheme 1, via nonadiabatic trajectory calculations. These calculations are important because they provide direct access to a molecular-level description of the statistical process leading to high or low QYs. However, they do not provide direct access to a QY theory for ultrafast reactions. In fact, while several hypotheses and even experimental verifications 37,65 of these hypotheses have been reported, the molecular factors controlling QYs in rhodopsins are still not established.
According to Weiss and Warshel 36 and Mathies and coworkers, 10 for fast impulsive (rather than diffusive) dynamics, there exists a correlation such that a larger velocity of the S 1 population moving toward the CI correlates with a higher reaction QY. Although the existence and relevance of such relationship, which is based on a Landau−Zener model 38,66−68 of the S 1 decay, is still under investigation, 69,70 it is apparent that a decrease in the time required to reach the CI enforces coherent S 1 dynamics. While a coherent dynamics regime appears to be necessary for the control of the QY value, 71 the molecular-level mechanisms that could exploit such a regime is still substantially unknown.
As mentioned above, QYs can be defined as N p /N t (i.e., the number of trajectories reaching the photoproduct P divided by the total number of trajectories launched). We now examine a set of results aimed at the understanding of the dependence of the photoisomerization "direction" (i.e., either toward the photoproduct or toward the reactant) of single nonadiabatic trajectories on critical molecular coordinates. The central idea is that the time evolution of certain coordinates determines in which way the reactive double bond reconstitutes immediately after the S 1 decay. 27,39 The stereochemistry of such reconstitution depends on: (i) the nature of the electronic wave function and (ii) the overlap between the fragment π-orbitals containing the paired electrons, which, in turn, depends on the chromophore geometry.
In Figure 13, we introduce the molecular modes affecting the π-overlap across the photoisomerizing double bond, and in Figure 14, we introduce the relationships between these modes and the stereochemical directionality of the double-bond photoisomerization. We start by considering a structure where the HOOP and torsion are close to −90°(see the central structure of Figure 14A). At this highly twisted structure, the overlap of the π-systems of the C1C2 and C2C3 fragments (represented by p-orbitals at the C2 and C3 centers respectively) is close to zero. We assume that this structure corresponds to a S 1 point close to a CI and therefore about to hop from S 1 to S 0 . We also assume that the torsion is slow with respect to the HOOP, consistently with the associated frequencies of these modes (see Figure 13 bottom).
We then look at the effect of the vibrational deformation toward negative (see Figure 14A top) and positive (see Figure  14A bottom) HOOP values, which increases the orbital overlap between the two π-systems. One finds that negative and decreasing HOOP values prompt formation of the E reactant (large b-b′ and a-a′ overlap), while positive and increasing values prompt formation of the Z product (large b′-a and b-a′ overlap). Notice that these rules are valid for the M-helical central structure of Figure 14A and associated counterclockwise isomerization motion (i.e., a motion where the negative C1− C2−C3−C4 value is increasing going from −180°to −90°). The rules would invert for the corresponding P-helical structure and clockwise isomerization featuring positive rather than negative HOOP and torsion. We conclude that the isomerization directionality (i.e., toward reactant or product) may be controlled by the HOOP direction of change at the hop point. These results have been computationally verified in several works, 39,40,72,73 where it is clearly shown that the majority of the trajectories follow the above rules not only in terms of velocities (i.e., sign of the HOOP velocity in degrees/fs) but also HOOP acceleration, which is proportional to the HOOP gradient as reported in Figure 15. 39 In Figure 14B, we show how, according to the above rules/ theory, the E (C1−C2−C3−C4 ca. 180°or −180°) to Z (C1− C2−C3−C4 ca. 0°) photoisomerization directionality can be affected by the HOOP frequency and phase. Such variations may be imposed by a change in the retinal chromophore environment, such as those imposed by mutations of the protein, thus influencing the QY magnitude and isomerization selectivity. As stressed above, the hop occurs when the two orbitals are nearly orthogonal (i.e., when the structure is close to a CI as mentioned above). In the reference diagram on the left of Figure 14B, the HOOP torsion oscillates for more than a period before entering the hop region (dashed vertical red line) then it enters it in a situation where the HOOP value is increasing (going to less negative values and thus displaying a positive velocity), while the C1−C2−C3−C4 torsion value is constantly decreasing (i.e., going from the E isomer reactant where the C1−C2−C3−C4 angle is nearly −180°as we are now starting from the E form of the chromophore to the CI that has a C1−C2−C3−C4 angle close to −90°). In this situation, the HOOP has a negative velocity and would hinder the formation of the Z isomer. Thus, the majority of the trajectories are expected to lead to the E isomer reactant even if a non-negligible number of trajectories may lead to the product. Figure 15 shows a computational demonstration of this effect for Rh where the C1−C2−C3−C4 torsional angle discussed above is now the C10−C11−C12− C13 angle of the rPSB11 chromophore and the isomerization proceeds from Z (starting close to 0°) to E.
The directionality of the nonadiabatic trajectory discussed above can be modified by changing either the frequency or phase of the HOOP coordinate relative to the torsion coordinate. These changes are schematically illustrated by the diagrams on the right of Figure 14B. Changing the frequency (see centralright diagram) to a lower value allows the HOOP mode to reach the CI region, while it is becoming more negative, and therefore acquires a positive velocity. In accordance with the above rules, Figure 13. Schematic representation of the geometrical coordinates modulating the overlap between the two π-systems (represented by porbitals at C2 and C3 centers) that interact to form the isomerizing double bond of a retinal chromophore. These coordinates can also be thought of as vibrational modes with a given excited-state frequency. The frequency is low for the torsion mode (i.e., the C1−C2−C3−C4 dihedral angle), larger for the HOOP mode (i.e., the difference between the H−C2−C3−R and the C1−C2−C3−C4 dihedral angles) and even larger for the d CC mode (d CC is part of the BLA mode mentioned above).  this has the consequence of leading to reactive trajectories. The same effect can be achieved by shifting the phase of the HOOP mode (see Figure 14B, bottom-right). In conclusion, one way of controlling QYs would be to "engineer" specific combinations of frequency and phase changes in the HOOP mode.
Recently, research work has been caried out to explore the effect of external forces applied to a chromophore (e.g., such as in a polymer chain incorporating one photochemically reactive unit). 74 The authors report nonadiabatic trajectory simulations, showing a marked and counterintuitive increase in double-bond E to Z photoisomerization quantum yield of a reduced model of the rPSBAT chromophore when pulling forces are applied axially to the terminal C and N atoms of the molecule. 75 The effect is successfully interpreted in terms of the phase relationship seen above.

Summary: Biological Chromophores and Their Modulation
On the basis of the rPSB results discussed in the previous sections, we may conclude that nature has implemented ultrafast double-bond isomerizations through different "selection processes". For instance, the chromophore is selected in such a way to provide a low-lying CI (or IS segment) which could be accessed easily after having absorbed a photon of visible wavelength and could preferentially mediate the isomerization of a specific double bond (e.g., at the center of the conjugated chain). The chromophore must also provide a space-saving reaction coordinate capable to cope with spatial constraints as the isomerization may have to occur in a tight molecular cavity. Other selection processes are related to the "modulation" of the spectroscopic and photochemical properties of the chromophore. Therefore, a chromophore with a spectroscopic/reactive excited state featuring a large charge-transfer character is likely to be very sensitive to the environment electrostatics and, therefore, to the amino acid composition of the hosting protein. These modulable properties ultimately enable the optimization of the biological function including light sensitivity. In fact, above we have reviewed computational evidence showing that the chargetransfer character may be involved in the modulation of both static properties, such as excitation energies (e.g.λ max ), reaction path slopes, conical intersection stability, and thermal isomerization barriers, as well as dynamic properties, such as excited state lifetime and quantum yields. Similarly, the chromophore geometrical "manipulation" (e.g., the pretwisting of single and/ or double bonds) may be involved in excitation energies (e.g., λ max ) modulation as well as in the selection of the reactive double bond.

RETINAL PROTEINS
In sections 2.1, 2.2, and 3.1, we have introduced Rh as a prototypal system for understanding the mechanism of doublebond isomerization in photobiology. As a result, we have extensively discussed the cis-trans isomerization of the C11 C12 double bond in Rh, the geometric and electronic structure of its CI, the reaction path followed during photoisomerization, and we also presented trajectory calculations simulating the isomerization of Rh. However, rhodopsins are a very large family of photoreceptor proteins found in all domains of life, where they serve a diverse array of functions. 76 In humans alone, rhodopsins are known to be responsible for color vision, dim light (scotopic) vision, and nonvisual light perception for regulation of circadian rhythm through a pigment called melanopsin. 77 As also mentioned earlier, not all rhodopsins have the same photochemistry and photophysics. In this section, we will present computational studies on other types of rhodopsins and attempt to draw a comparison between different types of rhodopsins. This section is divided into subsections focusing on conical intersections (section 4.1), reaction paths (section 4.2), and trajectories (section 4.3).

Conical Intersections
In section 3.1, we have discussed the geometrical and electronic structure of the CI of Rh as a reference vertebrate visual pigment. As shown in Figure 16A, the structure and stereochemistry of the CIs driving the rPSB11 isomerization in squid rhodopsin (sqRh) and in human melanopsin (hMeOp), an invertebrate visual pigment and a vertebrate nonvisual pigment, respectively, were found to be similar to the Rh CI. The comparison of these CI structures with the corresponding S 0 equilibrium structures (for , see also below). Green ▼ denote trajectories yielding the all-trans photoproduct (reactive trajectories), while red denote unreactive trajectories that return to the 11-cis configuration. The bar graphs count the number of events in the corresponding left and right halves of the diagrams. Notice that, consistently with the rules discussed in the text and illustrated in Figure 14A Rh see the top-right structure in Figure 16A) provides information on the rPSB11 changes occurring before S 1 decay. It is shown that, in all three cases, rPSB11 undergoes the same clockwise (CW) twisting of the C11C12 reactive double bond accompanied by a nearby counterclockwise (CCW) twisting of the C9C10 double bond. 24 The sense CW and CCW rotation is defined by the motion of the forward or backward carbon atoms of the isomerizing double bond (e.g., C13 and C10 atoms in the Newman projection on the right of Figure 16A).
The CI structures of animal rhodopsins can also be compared with those of microbial rhodopsins. More specifically, while the CI structures of the rPSBAT or rPSB13 chromophores of archaea and eubacteria rhodopsins 78,79 drive a C13C14 rather than C11C12 isomerization, the documented S 1 structural changes are similar to the aborted bicycle-pedal coordinate driving the isomerization in animal rhodopsins. For instance, the sensory rhodopsin from the cyanobacteria Anabaena (ASR), which exists in two forms (ASR 13C and ASR AT ) featuring a 13-cis,15-syn rPSB13, and a 13-trans,15-anti rPSBAT chromophore, respectively, 80,81 displays CIs characterized by a large twisting of the C13C14 and C15NH bonds with a minor twisting of the C11C12 bond (see Figure 16B). 79 The value of the twisting suggests that C15NH will revert, as expected for an aborted process, to the original conformation immediately after the decay to S 0 .
In spite of their differences, the CIs of Rh and ASR appear to drive isomerization processes consistent with alternated CCW and CW twisting deformations along the conjugated chain. In fact, the stereochemistry of the ASR 13C and ASR AT CIs can be described in terms of a C13C14 CCW and a C15NH CW twisting. This is consistent with the hypothesis that, in general, the C9C10 and C13C14 twists in a CCW direction, while the C11C12 and C15NH twists in a CW direction. Such a hypothesis is supported by the CI structure computed for the archaea light-driven proton pump bacteriorhodopsin (bR) that, in its light-adapted form (bR LA ), features a 13-trans,15-trans rPSBAT chromophore. In fact (see Figure 16C), the bR CI is consistent with a bicycle pedal mechanism involving a C11 C12 CW and a C13C14 CCW twist with a minor CW twist along C15NH. In conclusion, while the regiochemistry of the bicycle-pedal isomerization seems to be distinct for animal and microbial rhodopsins, the direction of the twisting motion is found to be consistent at least for the few CIs currently investigated.

Photochemical Reaction Paths in Animal and
Microbial Rhodopsins. The MEP driving the photoisomerization of the rPSB11 or rPSBAT chromophores have been computed for both gas-phase reduced models of the chromophore as well as for QM/MM models of different rhodopsins. We will begin by examining important features of the MEP profile for the photochemical conversion of the visual photocycle intermediate bathoRh to Rh. This back-isomerization process has been studied experimentally 82 (and recently by using time-resolved spectroscopy 83,84 ) as well as computationally (see for instance ref 85), because it is one of the fundamental processes that could limit the sensitivity of a photoreceptor. This has been illustrated in Scheme 1 by a dashed vertical arrow representing a photon absorption from the photoproduct. Its detailed investigation is therefore of interest. What we report below is only a partial account of the work focusing on the reaction paths.
The analysis of the S 1 MEP of this process (see Figure 17A) provides a static description of the photoisomerization mechanism. It also provides a basis for interpreting the results of nonadiabatic trajectory calculations (see section 4.3) in terms of forces (the slopes of the S 1 to S 0 potential energy surfaces) driving the chromophore changes. These forces control the progression of the S 1 population toward the CI and, after decay, its relaxation toward either the primary S 0 photoproduct or the original reactant. 85 Consistently with the data presented in the right part of Figure  17A, the MEP starting from bathoRh initially develops along the central C11C12 double bond turning it into a stretched C11− C12 single bond. As anticipated above, this motion is part of a more complex BLA stretching coordinate. Such relaxation leads rapidly to a mixed coordinate incorporating a HOOP wagging deformation of the two hydrogen atoms linked to the carbons of the reactive double bond. 27 As a reminder, HOOP is defined as the difference between the angle β, corresponding to the H− C11−C12−H dihedral angle and the angle α corresponding to the skeletal C10−C11−C12−C13 dihedral. Due to the S 1 forces initially acting along the C11C12 double bond and the H− C11−C12−H dihedral, the HOOP wagging changes more rapidly than the skeletal C11C12 twisting, reaching its maxima at a point where the skeletal deformation is minimal. This rapid change along the MEP is explained in terms of the S 1 antibonding character between the two π-systems overlapping through the C11−C12 bond. Such a character is caused by the Φ CTdominanted S 1 electronic structure which, effectively, shifts the spin-pairing of the electrons of one position along the chromophore's conjugated chain. Consequently, a force acting on the nuclei develops along the critical BLA stretching, HOOP wagging, and C11C12 twisting modes. 27 As reported above ( Figure 4, Figure 13, and related text), the fastest evolving coordinates are the ones that more rapidly decrease the overlap: the BLA (dominated by the C11C12 expansion) that increases the distance between the centers of the two π-systems and HOOP (dominated by the change in the angle β) that rapidly increases the angle between the orbital centers. In conclusion, while, as evident from Figure 17A the S 1 potential energy constantly decreases along the MEP coordinate of bathoRh, the coordinate itself involves several different modes and is therefore highly curved.
The energy profiles along the S 1 MEP computed for vertebrate, invertebrate, eubacteria, and archaea rhodopsins are all substantially barrierless. 21,24,79,86 However, the MEP computed for ASR AT and bR LA featuring a rPSBAT chromophore show a flat S 1 energy profile immediately following BLA relaxation. 79,86 This flat region has been proposed to be responsible for the observed longer S 1 lifetimes and photoproduct appearance times with respect to rhodopsins hosting either the rPSB11 or rPSB13 chromophore. 79,86 A comparison between the energy profiles computed along an approximated MEP (i.e., a relaxed scan) is given in Figure 17B for the case of ASR AT and ASR 13C .
The analysis of the S 1 MEPs described above suggests that the stereochemistry of the S 1 relaxation of vertebrate, invertebrate, eubacteria, and archaea rhodopsins is the consequence of a small twisting of the initial S 0 equilibrium structures. 21,24,79,86 Such twisting would be responsible for initiating the relaxation in specific CW or CCW directions, ultimately leading to a specific CI. For instance, as shown in Figure 16A, the S 0 equilibrium geometry of Rh includes a rPSB11 with the C11C12 bond twisted by ca. 6°in the CW direction. Such a small twist in the chromophore S 0 equilibrium structure is also observed for sqRh and hMeOp, which both display a 10°twist around the C11 C12 bond. 24 Similarly, the microbial rhodopsins, ASR AT , ASR 13C , and bR LA , show a S 0 equilibrium structure with a C13C14 twist of 8°−11°in the CCW direction (see Figure 17, panels B and C), consistent with the result of the MEP that shows a CCW isomerization coordinate. 23 The S 1 MEPs computed for ASR AT and ASR 13C both describe a twisting of the reactive C13C14 bond occurring in the CCW direction. When considering that ASR is photochromic (see Figure 18A) such that ASR AT and ASR 13C are interconverted via irradiation at two different wavelengths, 80,81 Strambi et al. 79 hypothesized that the ASR chromophore operates like a lightdriven single-molecule unidirectional rotor. 87 The corresponding rotation mechanism has been documented by following the S 1 MEPs reported in Figure 17B and by computing the K and L photocycle intermediates of both the ASR AT and ASR 13C forms. 79 The K and L structures are the result of MEPs that describe the S 0 relaxation from the CIs of the ASR 13C and ASR AT toward the ASR AT and ASR 13C photoproducts, respectively. The mechanistic picture emerging from such calculations is shown in Figure 18B, where the stereochemistry of the relevant double-bond regions is assigned using both the M and P helicity labels and the CW and CCW labels describe the twisting direction associated with either light or thermally induced isomerizations.
The rPSBAT chromophore was studied in the chimeric channelrhodopsin C1C2 more recently by two different research groups. 88,89 Both studies were based on relaxed scan calculations using the CASPT2//CASSCF methodology and reached similar conclusions. Consistently, two pathways were determined with opposite rotating C13C14 bonds. One pathway was found to have a low or no barrier in the excited state, while the other one was showing a larger barrier explaining two different time constants for excited-state lifetimes from transient absorption measurements. 89 4.2.2. Relationship between Thermochemical and Photochemical Isomerization Paths in Rhodopsins. As discussed in section I.1, the S 0 potential energy surface of rPSB11 may host two transition states (TS DIR and TS CT ) controlling the It also provides a definition of the dihedral angles α and β, which define the HOOP and C10−C11−C12−C13 twisting coordinates. The right part describes the coupling of the C10−C11−C12−C13 dihedral angle of the isomerizing double bond with the HOOP angle and C11−C12 length (notice that the x axis here is the torsion and not MEP coordinate; however, since the C10−C11−C12−C13 twisting is monotonically increasing along the MEP, the diagram also loosely represents the HOOP and stretching geometrical progression along the MEP from reactant to the bathoRh CI). The relationship between the progression along the MEP and the overlap between the π-systems (represented by hybrid p-orbitals) describing the bond breaking/ making of the C11C12 bond is also given. same thermal isomerization. This has also been demonstrated for gas-phase reduced models of rPSB11. 43,90 Figure 19 presents the TS DIR and TS CT structures computed in a bovine Rh QM/MM model. 56 TS DIR features the structure expected for a homolytically broken double bond, with two radical centers delocalized along orthogonal π-systems. Such a transition state was also located by Nemukhin and co-workers in a model of Rh employing DFT to describe the chromophore. 91 The TS DIR charge distribution correlates with that of the S 0 rPSB11 reactant with a +0.98e charge localized in the Schiff base region. In contrast, TS CT has most of its charge (+0.90e) located on the βionone region, thus more closely resembling the charge distribution of the vertically excited S 1 rPSB11 (compare the corresponding bubble diagrams in Figure 9B). It was found that TS CT has a computed activation energy of 34 kcal/mol and lies 11 kcal/mol in energy below TS DIR . While these energies depend on the level employed in the calculation, an increase in the level of theory invariably leads to a larger stabilization of TS CT . Therefore, it has been found that TS CT has full control of the thermal isomerization. 56 To understand why both TS CT and TS DIR exist in Rh, and how a single reaction could be mediated by two different transitions states, we need to revisit concepts from section 3.2, where we discuss the effect of the conical intersection on the topography of the ground state. There, we have discussed how the S 0 and S 1 states change their electronic character along a small loop centered at the CI and residing on the BP. In other words, in rPSB11, the S 0 electronic character passes two times (e.g., roughly every 180°) from a Φ DIR character to a Φ CT character. Due to the peaked topography of the CI, and the fact that TS CT and TS DIR lie on opposite sides of the CI along the BLA coordinate, means that TS CT and TS DIR are dominated by different electronic characters, Φ CT and Φ DIR , one associated with the S 1 state at the FC point and one associated with the S 0 state at the FC point, respectively. Therefore, the existence of the low-lying CI in Rh allows the existence of a transition state, TS CT , with the same electronic character as the first excited state of Rh. Such a transition state would not exist if it were not for the existence of the nearby CI.
How does isomerization through TS CT occur? In Figure 20A, we display a graphical representation of the change in π-electron density along the thermal isomerization path for Rh when passing through the lower energy transition state TS CT . It is apparent from inspection of the structure that along the longest path segment leading to TS CT , C11C12 remains in a "locked" Φ DIR dominated state featuring two π-electrons coupled across the breaking bond. However, due to the electronic nature of the region surrounding TS CT , the electronic structure suddenly changes to a Φ CT configuration as it is evident by the gap in πelectron density at C11C12. This leads to unlocking of the C11C12 double bond, similar to how it is unlocked in the S 1 state of Rh prompting photoisomerization. Of course, once the molecule passes the TS CT moving toward bathoRh, the  The sudden change in the electronic character is quantitatively described in Figure 20B. The positive charge on the β-iononering-containing moiety of rPSB11 is shown in terms of an area diagram. It is evident that in the TS CT region, the positive charge on such a moiety suddenly increases and then decreases after the transition state. The same analysis demonstrates that the alternative and higher-energy isomerization path passing through TS DIR displays a change in charge distribution substantially opposite to the one seen for TS CT .
The nature of the transition state controlling thermal isomerization could have important implications for visual sensitivity. The similarity between the Φ CT -dominated electronic structure of FC and TS CT , together with the results indicating that TS CT is controlling the thermal isomerization, provides a direct link between the rate of thermal isomerization (controlled by the energy of TS CT ) and the absorption wavelength (controlled by the FC S 1 energy). This link is illustrated in Figure 21A. The diagram (see the brown color position along the curves indicating a dominating Φ CT character) explains the link between E a T and ΔE(S 1 −S 0 ) and, in turn, between −log k (where k indicates the kinetic constant of the thermal isomerization reaction) and 1/λ max . Accordingly, any opsin red-shifting the absorption [i.e., decreasing the ΔE(S 1 −S 0 ) value] would simultaneously decrease E a T . This behavior, deduced from quantum mechanical calculations, 56 appears to be consistent with the "Barlow correlation," a prediction made by Barlow in the 1950's 92 and observed in a series of experiments 93−95 that demonstrate that indeed a relationship between the λ max value and the thermal activation kinetic constant (k) in different visual pigments exists. The results above provide a basis for the development of a pigment-noise theory. Such a theory would hold that, in a series of related pigments (e.g., the rod visual pigments), the thermal noise (i.e., the thermal isomerization competing with the photochemical isomerization) would decrease when passing from red-light perceiving pigments to blue-light perceiving pigments. In other words, higher visual sensitivity may be achieved by pigments that absorb in the blue. This could be one of the reasons why dim-light vision is tuned to the blue edge of the visible spectrum. 56 The computational results reported above suggest that a Barlow-like relationship can be reproduced by QM/MM models of rhodopsins. This is illustrated schematically in Figure 21B where a negative charge represents the hypothetical effect protein environment in two different visual pigments carrying the Figure 20. Evolution of the electronic structure along TS trajectories in Rh. (A) Schematic representation of the S 0 potential energy surface driving the isomerization of the protein embedded rPSB11 chromophore and of the progression from Rh to bathoRh along a TS trajectory. The three rPSB11 structures display the calculated shape of the π-electron density along the conjugated framework for the trajectory passing through TS CT . (B) Change in S 1 and S 0 potential energy (red and blue lines, respectively) and charge distribution (orange area) along the two TS trajectories passing through TS CT and TS DIR . Each TS trajectory is a composite trajectory defined by releasing two classical trajectories from the optimized TS with minimal velocities toward the reactant and product, respectively, but then reversing the reactant trajectory to generate the plot shown. Figure 21. Relationship between ΔE(S 1 − S 0 ) and E a T in rhodopsins. (A) One-dimensional representation of the relationship between the energy profiles and electronic structure of the S 1 and S 0 isomerization paths when TS CT is controlling the S 0 process. The brown segments along the paths are dominated by a Φ CT configuration. The green segments are instead dominated by a Φ DIR (diradical or covalent) configuration. (B) Schematic representation of the simultaneous effect of the protein environment on the ΔE(S 1 − S 0 ) and E a T values. The negative charge in the "red" and "blue" panels schematically represents the electrostatic effect of the cavity hosting the chromophore. rPSB11 chromophore and absorbing light of shorter and longer wavelength, respectively. The bottom-left panel shows that the negative charge is distant from the β-ionone-ring moiety, and therefore, the S 1 FC and the S 0 TS CT structure are only slightly stabilized with respect to the surrounding regions of the S 1 and S 0 potential energy surface, respectively. On the other hand, the top-left panel shows a situation where the negative charge is closer to the same β-ionone-ring moiety. In this situation, the S 1 FC structure will be stabilized with respect to S 0 , leading to a redshifted absorption maximum. At the same time, TS CT will also be stabilized leading to a lower barrier and faster thermal isomerization rate. By interpreting the change in the position of the negative charge in terms of a change in the chromophore cavity, one achieves a molecular-level understanding of the Barlow relationship. Of course, a behavior not consistent with the Barlow correlation would be observed if the isomerization was controlled by TS DIR .
A number of computational results have been reported that support the mechanism reviewed above. In Figure 22A, the −log k vs 1/λ max correlation has been computed by modeling a set of Rh mutants incorporating two different types of chromophores. These comprise 5 bovine Rh mutants (hosting rPSB11, sometimes referred to as the A1 chromophore) and 6 derivatives featuring the longer wavelength absorbing 3,4-dehydro-retinal (A2) chromophore. The calculated correlation can be compared directly to experimental data collected by Ala-Laurila et al. 94 A correlation of the same type has been obtained by looking at the computed E a T and 1/λ max values for a set of more distant (in the sense of the sequence diversity) set of pigments with respect to Rh mutants. The set comprises four visual pigments from cottoid fish from lake Baikal (see Figure 22B), which are also compared with bovine Rh. In this case, the apparent linear correlation between the two quantities has also been analyzed in terms of two effects. The first is related to a modulation of the chromophore geometry, while the second is a pure modulation of the protein electrostatic field acting on the chromophore. It was demonstrated that the modulation in the electrostatic field is mainly responsible for the correlation between E a T and 1/λ max . The validation of a Barlow-type correlation in the lake Baikal fish flock provides an explanation for the fact that the species with the most blue-shifted absorption reside at the abyssal zone of the lake (400−1500 m depth), where the quantity of available light is at its minimum. In this situation, the visual pigment must have a high thermal isomerization barrier to limit the thermal noise. On the other hand, the species thriving in the maximally illuminated littoral zone of the lake (1−5 m depth) have red-shifted absorption maxima. Figure 22C shows that a Barlow-type correlation may also exist for distant rhodopsins that do not even serve the same biological function. In fact, using the same type of QM/MM model, it has been demonstrated that the computed E a T and 1/λ max values for a small set of vertebrate, invertebrate, and microbial rhodopsins may also be directly proportional, again having the most redabsorbing member of the set featuring the lowest thermal isomerization barrier. In conclusion, the existence of a Barlow correlation (i.e., a link between absorption wavelength and thermal isomerization rate) appears to be present in all types of rhodopsins.  22 This pioneering effort introduced the space-saving reaction coordinate hypothesis. The simulated trajectory pointed to a rPSB11 isomerization mechanism dubbed bicycle-pedal (BCP in Figure 23A) because of the simultaneous twisting of two double bonds resembling the motion of the pedals of a bicycle. The mechanism originally proposed by Warshel involved isomerizations at C11C12 and C15NH and therefore is distinct from the version described in the Introduction based on more recent QM/MM models. In fact, recent studies point to a BCP at C11C12 and C9C10. Furthermore, it is aborted upon decay at the CI, yielding a so- called nonreactive or aborted BCP (NBCP) where only the C11C12 completes the isomerization but not C9C10.
In a seminal study, Warshel employed an empirical force field method (QCFF/PI) 97,98 that provided the S 1 and S 0 gradients necessary to propagate a nonadiabatic trajectory starting on S 1 . At that time, a crystallographic structure of rhodopsin was not available and therefore the steric interactions with the protein matrix were modeled by imposing geometric constraints on the chromophore atoms. It is remarkable to note that Warshel has predicted this reaction to complete in 200 fs before such an ultrafast time scale became experimentally accessible. Other trajectory calculations by Birge and Hubbard were favoring a one bond flip (OBF) isomerization at C11C12. 99 They predicted also a longer excited-state lifetime of ca. 2.0 ps. Using an improved protein-environment model, Warshel and co-workers refined their isomerization model. 100 The original mechanism was modified by Warshel et al. such that only one double bond, C11C12, isomerized successfully and a second one, C9C10, becomes significantly twisted but then reverts back to the original configuration. This is essentially the NBCP mechanism documented by modern QM/MM models based on multiconfigurational quantum chemical methods.
In the course of development of quantum chemical methods and improvement of the computational hardware, more accurate molecular dynamics simulations become accessible. However, even more important was the availability of the crystal structures as a starting point for trajectory calculations. The trajectory calculations are the closest simulations to experiment because it considers the effect of kinetic energy in contrast to MEP calculations. Further, molecular dynamics simulations are free of bias and do not need a priori knowledge of the reaction mechanism as the outcome is determined by the initial conditions. Hence, with the advance of QM/MM trajectories, various proposed mechanisms of the photoisomerization could be tested ( Figure 23). Such mechanisms evolved from chemical intuition or static quantum chemical calculations [e.g., Liu et al. suggested a volume-conserving isomerization of neighboring single and double bonds in the so-called "Hula-twist" (HT in Figure 23A) mechanism]. 101−103 More recently, a "double bicycle-pedal" 86 or "folding-table" 104 (FT in Figure 23A) mechanism was proposed for bR.
The mechanisms above have mainly been invoked to explain the ability of the rPSB chromophore to isomerize in a spacesaving way inside a tight binding pocket that restricts largeamplitude motions of the embedded chromophore. However, in addition to the steric consideration, electronic factors also play an important role in determining the mechanism of isomerization. The physical origin of the preferred excited-state BCP-like relaxation is schematically illustrated in Figure 23B, where it is shown that such a mechanism is driven by the Φ CT electronic structure of the excited-state chromophore. In the same figure, we also show how this mechanism is different for polyenes where the isomerization is driven by an excited state with diradical Φ XDIR electronic structure (see section 3.1).
Rh was the major focus of trajectory calculation for several reasons: (i) the crystal structure became available in 2000 41 and (ii) shortest excited state lifetime leading to short computational times. The exact details of the computational methodology varied in the studies in several aspects. A broad range of the quantum chemical methods was used, starting from restricted open-shell Kohn−Sham (ROKS) method 30 to wave function based methods. The latter methods were the preferred choice and included semiempirical floating occupation molecular orbital configuration interaction (FOMO−CI) 34,105 but used mainly CASSCF. 21,26,27,29,32,106 To compensate for the lack of dynamic electron correlation in CASSCF, Frutos et al. have scaled the forces linearly to qualitatively match those obtained at the CASPT2 level. 21 This approximation was supported by the finding that only the steepness of the excited-state potential was different between the energy profiles described at the CASSCF and CASPT2 levels of theory. More recently, Martínez and coworkers have introduced an empirical correction in a so-called α-CASSCF approach that scales the splitting between state-specific and state-averaged CASSCF. 107 The correction was determined to match the multistate CASPT2 reference energy and gradients.
Despite the differences in the computational setup, the majority of the trajectories were found to reach an S 0 /S 1 degeneracy on the scale of 100 fs. 21,27,29,32 and to form predominantly the bathoRh intermediate. The reaction coordinate was composed of the sequential inversion of BLA and torsion of the C11C12 bond as already determined in MEP or relaxed scan calculations discussed in section 4.2. However, also the 9-cis isomer (isoRh) and a genuine BCP product were found as byproducts in some trajectories, 26,29 which can only be observed in molecular dynamics simulations. Another new insight gained from trajectory simulations was demonstrated by Ishida and co-workers 106 on Rh and isoRh. The results showed vibrational coherence only in the case of Rh but not for isoRh. Further, the authors explained the longer excited-state lifetime of isoRh with respect to Rh. In particular, it was found that the former stays longer on the excited state with respect to Rh because of back-and-forth twisting of the dihedral angle in S 1 , while in Rh this does not occur (see Figure 24). Two amino acids were found to be responsible for this extended S 1 lifetime: Tyr268 and Thr118. Furthermore, the authors proposed a new mechanism: the "wring-a-wet-towel" motion. However, since this mechanism is defined by rotation of C9−C10 and C11−C12 in opposite directions, it appears to be a variant of the BCP mechanism.
Further insight into the isomerization in isoRh was provided by Garavelli and co-workers. 40 Their QM/MM molecular dynamics simulations have located two deactivation pathways starting from rPSB9. One path was leading to a CI that produced a di-cis photoproduct, while the other path reached a different CI that can generate the rPSB11 and rPSBAT photoproducts. The latter pathway was found to have an increased interaction with Tyr268, in line with the work of Ishida and co-workers 106 and, therefore, a longer lifetime on the excited state than the former pathway. Hence, the authors interpreted the presence of two different isomerization directions to be responsible for the biexponential decay of isoRh.

Bacteriorhodopsin and other Microbial
Rhodopsins. The rPSBAT chromophore in bR and microbial rhodopsins, in general, has a longer excited state lifetime than rPSB11 that is thought to originate in the presence of an excited state barrier. The longer simulation time and the accurate assessment of the barrier present a challenge for accurate excited state dynamics simulations. The first trajectories were computed by Birge and co-workers using semiempirical molecular dynamics. 108 Their calculations of the rPSBAT to rPSB13C isomerization in bR were performed at a time before full structural data for bR became available. They predicted that a transoid minimum on S 1 may trap some of the excited-state population before decay. The first high-resolution (1.55 Å) crystal structure of bR became available in 1999, thus providing a more realistic starting point for simulations. In 2001, Warshel and Chu 109 used a QM/MM setup with the empirical valence bond method 110 to describe the S 1 potential of the chromophore in the QM subsystem. In their simulations, they also accounted for the effect of polarization of the protein environment by the QM atoms.
The first simulation of rPSBAT in bR carried out with a QM/ MM model based on ab initio multiconfigurational quantum chemistry was reported two years later by Schulten and coworkers. 111 The QM subystem in their calculations consisted of a reduced model, a three double-bond portion of the retinal from C11 to Cε, and described at the CASSCF level of theory. Consequently, 6 π-electrons in 6 π-orbitals were included in the active space. The remaining atoms of the retinal chromophore, together with the protein, were calculated with a classical empirical force field.
Morokuma and co-workers 112 studied bR and rPSBAT in methanol using the ONIOM scheme for QM/MM molecular dynamics and CASSCF as an electronic structure method. The trajectories were computed in bR as well as in methanol that served a reference system. The NBCP mechanism was found in both environments. A different comparison of the environmental effect on rPSBAT photoisomerization was carried out by Punwong et al. 113 who also studied the gas-phase reaction. The QM subsystem was described at the semiempirical FOMO− CI 105 level of theory. The isomerization in methanol was found to be much longer than in the gas phase. The reason for this longer time is an S 1 barrier which slows down the isomerization. They found that the electrostatic (not the steric) interactions with methanol are the primary cause of this torsional barrier. More specifically, the barrier originates from the solvent being oriented to stabilize the Franck−Condon region of the chromophore, and therefore, isomerization requires the solvent to reorient as the charge distribution changes during torsion on the excited state, slowing down photoisomerization compared to the isolated rPSBAT chromophore. Punwong et al. also discuss the photoisomerization quantum yield in the gas phase and solution. They agreed with the hypothesis of Martínez and coworkers 114 that the CI topography has an important role in determining the branching ratio of product versus reactant formation (i.e., the photoreaction quantum yield). 113 In order to compare to different microbial rhodopsins, Punwong et al. 115 simulated the excited-state dynamics of bR and halorhodopsin (hR) using the same QM/MM model as mentioned above. 113 The authors found that bR has much shorter excited-state lifetime than hR. The different isomerization rates were attributed to different electrostatic stabilization of the positive charge of rPSBAT on S 1 in hR and bR. 115 4.3.3. Comparative Studies of Animal and Microbial Rhodopsins. A comparison of the S 1 dynamics of Rh and isoRh has been reported by Ishida and co-workers, 106 and as seen in sections 3.3 and 4.3.1, by Garavelli and co-workers. 40 In both cases, the slower isomerization of isoRh with respect to Rh has been successfully simulated. Their results indicate that the steric hindrance imposed by the environment on the C9C10 bond of rPSB9, but not on the C11C12 double bond of rPSB11, is the cause for the different isoRh and Rh dynamics. Since isoRh and Rh have the same protein cavity, these findings suggest that the selection of a specific rPSB chromophore isomer is an important factor for the modulation of the rhodopsin function.
While the rPSB9 chromophore is mainly known as a byproduct of the in PSB11 photoisomerization in Rh, distinct rPSB chromophore isomers characterize the dark states of microbial and animal rhodopsins (also called type I and type II rhodopsins). 1 In fact, type I and type II rhodopsins are characterized by the rPSBAT and rPSB11 chromophore, respectively. Both chromophores are responsible for the activation of the protein function that is triggered by the photoisomerization of the C13C14 and C11C12 double bonds. However, the fact that rPSB11 is thermodynamically unstable with respect to rPSBAT indicates that rPSB11 might have been selected by biological evolution to carry out an otherwise inaccessible function. To investigate whether such a function exists, a recent study has compared the photoisomerization mechanisms of rPSBAT and rPSB13 in the microbial rhodopsin ASR (see above and section 4.2.1) and of rPSB11 in Rh using trajectory computations. 23 The simulations unveil electronic effects achievable in rPSB11 but not in rPSBAT that may explain the higher photoisomerization speed reported for animal rhodopsins.
Consistently with previous studies, 21 it was assumed that during the 200 fs simulation time a FC trajectory (see section 3.2) describes the average evolution of the corresponding S 1 population. This has been partially assessed by running small sets of trajectories starting with different initial conditions and comparing them with the single FC trajectory. 21,24,32 As shown in Figure 25A, the ASR AT trajectory remains on S 1 for the entire duration of the trajectory, whereas ASR 13C reaches the CI after ca. 150 fs. The estimated S 1 lifetimes are consistent with the available transient spectroscopy measurements 116 and computations, 117 showing that the ASR AT photoisomerizes at least five times more slowly than ASR 13C . Moreover, ASR 13C (see Figure  25B) seems to react in a time scale that is only slightly slower than that computed for Rh 21,24 (Figure 25C), consistent with the observations. 32,70 The slower photoisomerization of ASR AT is also consistent with its observed larger fluorescence quantum yield. 70 In fact, starting 25 fs after excitation, the authors compute an average oscillator strength in the 0.8−1.6 range (see Figure 24D). Taking the average ΔE S1−S0 values over the 175 fs interval, they predict a fluorescence λ max of 669 nm, which is comparable with the observed values (i.e., a λ max of 700−710 nm with a shoulder around 650 nm). 70 It is apparent from inspection of the energy profiles in Figure  25 (panels A and B) that both S 1 and S 2 are involved in the isomerizations of ASR AT and, partially, of ASR 13C . Indeed, the S 1 and S 2 states become nearly degenerate ca. 25 fs after photoexcitation. However, while ASR AT travels along the degeneracy for the full duration of the trajectory, ASR 13C abandons the degeneracy region ca. 50 fs after entering it and then evolves toward the CI. In contrast, no degeneracy is found for Rh, which evolves toward the CI directly after photoexcitation. A similar behavior is found for other animal rhodopsins that feature the rPSB11 as their chromophore (e.g., sqRh and hMeOp). These data support a relationship between the magnitude of the estimated S 1 lifetime (defined as the time it takes to reach the CI, τ) and the extent of the S 2 /S 1 mixing. Inclusion of more excited states in the calculations does not change this conclusion.
The τ values of ASR, together with those previously reported for Rh, sqRh, and hMeOp 24 provide information on a dynamic factor that may be implicated in the rhodopsin function (e.g., light sensitivity). Whereas an increase of τ may be associated with a decrease in slope of the S 1 potential energy surface accelerating progression toward the CI, in ASR AT and other rPSBAT-hosting rhodopsins the occurrence of the S 2 /S 1 degeneracy changes the dynamics. The S 1 relaxation would thus become slower and more diffusive with the chromophore exploring the rugged potential energy surface shaped by several S 2 /S 1 avoided crossings (see inset in Figure 16B).
The results above indicate that whereas a two-state isomerization model (Scheme 3A) applies to fast (<150 fs) reacting systems incorporating the rPSB11 or even rPSB13 chromophores, a three-state model (Scheme 3B) would apply to slower rhodopsins hosting rPSBAT. A three-state model has been previously reported 120 for the rPSBAT-hosting Archaeal bacteriorhodopsins (bR), which displays slow and complex S 1 dynamics. 120 However, in such early three-state models, the S 2 diradical state ultimately drives the progression toward the CI Figure 25. QM/MM trajectories of ASR AT , ASR 13C , and Rh computed at the two root state average scaled-CASSCF/Amber (black lines) level of theory and corrected at the CASPT2 21 and XMCQDPT2 levels. 118 (A) S 0 (red ◆ ), S 1 (green ▲), and S 2 (blue •) CASPT2//CASSCF/ Amber energy profiles along the ASR AT trajectory. The main out-ofplane (deviation larger than ±5°) dihedral angles of the chromophore S 0 equilibrium structure are given at the top. The vertical dashed arrow represents weak fluorescence.  ( Figure 4C, inset). Although this seems to be the case in polyenes, 121 it is inconsistent with the structure of protonated Schiff base CIs. 122 Instead, the three-state model based on the more recent work of Luk et al. 23 involves the S 2 state only along part of the isomerization path, where it dips down and results in a slowing down of the S 1 isomerization dynamics. Such a threestate model may also apply to other rPSBAT-hosting rhodopsins such as bR and Channelrhodopsin. In fact, this would be consistent with the reported bR 86 and, recently, Channelrhodopsin 88 reaction paths. It would also be consistent with the similar dynamics observed for ASR AT and light-adapted bR (hosting exclusively rPSBAT) 116,123 as well as for ASR 13C and the corresponding rPSB13-containing bR. Finally, the same threestate model may also control the S 1 dynamics of unsubstituted and certain substituted rPSBAT chromophores in solution that display S 1 picosecond lifetimes. 124 The presence of an S 2 /S 1 degeneracy in solvated rPSBAT models has been proposed 125 and held responsible for the observed solution dynamics of rPSBAT. 126 Similar data have been reported for rPSB11, 49,127 suggesting that the Rh cavity has an important role of removing S 2 /S 1 degeneracies along the photoisomerization pathway in order to optimize photoisomerization lifetimes (and, therefore, possibly also optimize photoisomerization quantum yields).

Summary: Conserved CIs and Reaction Coordinate Variability
Above, we have reported on the results of a set of comparative studies carried out to establish which mechanistic aspects of rhodopsin photoisomerization are conserved or undergo variation in different rhodopsins. It is found that the geometrical and electronic structures of the CI's are conserved and consistent with a 1B u /1A g CI (Figure 8) reached via a NBCP type of motion (see Figure 23). The nature of the motion along the reaction coordinate also appears to be conserved and starts, in general, with a BLA relaxation which then immediately couples with HOOP motion. The reactive double-bond skeletal isomerization motion activates after the initial relaxation. We have also pointed out that the stereochemistry of the isomerization is, for the examined rhodopsins, of the CW type for C11C12 and C15 NH and CCW for C9C10 and C13C14. Such stereochemistry appears to be consistent with the documented pretwisting of specific double bonds of the chromophore. We can therefore hypothesize that the degree and direction of pretwisting imposed by the protein cavity determines which double bond is selected for the reaction and in which direction, CW or CCW, it will twist.
The mechanism summarized above needs to be further detailed for the case of the reaction coordinate of an inverse photoisomerization reactions. For bathoRh to Rh phtoisomerization, which starts from a largely CW pretwisted bathoRh intermediate, the reaction coordinate follows a reaction path which is almost exactly inverted with respect to the Rh to bathoRh path. Therefore, notice that for a directly inverse (with "directly" we mean starting from the primary photoproduct of the reaction) process, the C11C12 bond rotates CCW and the C9C10 CW. This situation is markedly different from one of the "direct" and "inverse" reactions occurring in the cyanobacterial rhodopsin ASR. In this case, we have documented reaction coordinates, which are part of a rotary motion driven by the absorption of two-photons and where the isomerizing double bond moves always in the same CCW direction. This is due to the fact that ASR AT and ASR 13C are not directly interconverted but require additional thermal isomerization steps. In fact, they differ in the conformation of their side chain and stereochemistry of the C15NH bond.
In contrast with the above conserved features, the reviewed results indicate that there are differences between the energy profiles and excited state lifetime of rhodopsins featuring the rPSB11 and rPSBAT chromophore. In fact, wild-type rhodopsins with the rPSBAT chromophore (ASR AT and bR LA ) feature much flatter energy profiles. In section 4.3.3, we have provided evidence supporting the involvement of a different number of electronic states for these type of rhodopsins. rPSB11 rhodopsins only involve two electronic states (S 1 and S 0 ), while the rPSBAT mechanism dominating the photoisomerization of microbial rhodopsins also involves the unreactive S 2 state. We have seen that in these last rhodopsins S 2 partially mixes with S 1 leading to a flatter or, better, rugged S 1 potential energy surface consistent with a slower excited state progression toward the CI.
The reviewed results show that each documented CI may come together with two transition states (TS CT and TS DIR ) on the S 0 potential energy surface. The one featuring charge-transfer character, TS CT , is energetically favored and controls the protein thermal isomerization. Interestingly, the corresponding isomerization coordinate is very similar to the NBCP coordinate driving the excited state isomerization discussed above but has a different BLA value which displays a larger degree of single-bond/doublebond inversion. Finally, it has been found that the similar chargetransfer characters (documented by looking at the charge translocation) of the FC and TS CT structures provide a natural and direct explanation for observed and simulated Barlow-like relationships.

PHOTOACTIVE YELLOW PROTEIN
The photoactive yellow protein (PYP) is a 125 amino acid cytoplasmic blue-light photoreceptor controlling the photoavoidance response of the salt-tolerant bacterium Halorhodospira halophile. 128 PYP belongs to the Xanthopsin protein family and is its most studied member. 129 It is also the structural prototype for the PAS (PER-ARNT-SIM) class of signal transduction proteins. 129−132 PYP contains a deprotonated 4-hydroxycinnamic acid (or p-coumaric acid, pCA) chromophore linked covalently to the γ-sulfur of Cys69 via a thioester bond (see Figure 26). Thus, in contrast to the cationic rPSB chromophore found in rhodopsins, PYP has an anionic chromophore where the phenolic oxygen forms short hydrogen bonds with nearby tyrosine (Tyr42) and glutamic acid (Glu46) residues (see Figure  26B). Upon blue-light irradiation, PYP enters a fully reversible photocycle involving several intermediates formed on time scales ranging from a few hundred femtoseconds to seconds (see Figure  27). Although a full understanding of the details of the PYP photocycle is still lacking, it can be divided into three major events: (1) PYP undergoes photoisomerization of the anionic chromophore to transform from pG (i.e., the starting structure) to yield the pR intermediate, which is a structure with a twisted C10C11 bond (see Figure 27B) and a red-shifted absorption, (2) protonation of the anionic chromophore and structural relaxation leads to the pB intermediate, which has a blue-shifted absorption with respect to the starting conformation, and (3) the recovery of the initial structure. The various intermediates are shown in the photocycle in Figure 27.
The detailed photoisomerization mechanism of the PYP chromophore has been debated for years. The light-induced dynamics of several models of pCA have been studied experimentally and computationally both in gas-phase 133 and in solution. 134,135 In the following sections, we focus our

Conical Intersections
The isomerization of the PYP chromophore has been studied using the isolated chromophore (pCK − ) as a model. 137,139,142 These studies suggest the existence of two possible pathways for the decay of pCK − from S 1 to S 0 : α-twist and β-twist (see Figure  28). The corresponding sloped CIs are located in the vicinity of S 1 minima whose structures are shown at the bottom of Figure  28. The α-twist path induces a complete transfer of the negative charge from the phenolic ring to the alkene fragment where the C9−C10 bond length is found to be 1.47 Å. In contrast, the βtwist induces only a partial charge transfer in the same direction, leading to a stretched single bond (which has a bond length of 1.47 Å) rather than to a double-bond twisting. Such preliminary gas-phase studies suggest that protein residues or solvent molecules forming hydrogen bonds with the phenolic oxygen are likely to play an important role in the selection of the isomerization pathway and drive the S 1 relaxation along the twisting coordinate. More recent studies have characterized the S 1 /S 0 CI reached via a twisting about the C10C11 double bond, corresponding to the β-twist of pCK − . The CI has a C9−C10−C11−C12 dihedral angle of 85°, compared to 170°at the equilibrium geometry. 140 The CI is geometrically close to the experimentally observed S 0 intermediate I T (see Figure 27), which also has a twisting angle of 85°in the crystal structure. 143 In the same computational study, 140 it is shown that the optimized intermediate I T has a C9−C10−C11−C12 angle of 80°. However, the experimental assignment of the I T structure is   not without controversy. 144,145 While the computational study by Wei et al., 140 which was performed at the CASSCF level, report on a twisted I T intermediate, DFT calculations 144,146 indicate that the I T structure is structurally closer to the fully isomerized cis product with a C9−C10−C11−C12 angle of 21°. Kaila et al. 144 suggested that such conflict is the consequence of the arbitrary choice of loose restraints in the X-ray data refinement. However, Jung et al. 145 pointed out that the discrepancies might be due to the experimental conditions and crystal preparation, even if the same time-resolved X-ray crystallography were employed. While the exact origin of discrepancies remains to be understood, 145 the most recent femtosecond time-resolved X-ray study shows no evidence for such an ∼90°twisted intermediate but finds an isomerization that leads directly to an ∼30°twisted structure (pR 0 in Figure  27B). 141

Reaction Paths
In order to understand the photochemical properties of the PYP chromophore, the S 1 potential energy surface of the pCK − gasphase model was studied at the CASSCF level. 139 As anticipated above, two excited-state minima were identified (see Figure 28): the α-twisted minimum, which features a 90°twisted bond adjacent to the phenol ring, and the β-twisted minimum, which features a 90°twisted ethylenic bond. In pCK − one finds a negligible barrier for reaching the α-twisted S 1 minimum from the FC region, whereas there is a significant barrier for reaching the β-twisted minimum. Hence, in the gas-phase, the main relaxation channel on S 1 would proceed along the α-twist coordinate. Similar results have been found at the CC2 level of theory. 147 Since isomerization in the protein proceeds exclusively along the β-twist, these results suggest that the PYP apoprotein controls the isomerization path making the β-twist the energetically favorable reaction coordinate. Another study by Gromov et al. at the CC2 level focused on a pCK − chromophore where the phenolic oxygen is stabilized by a number of water molecules. 138 Their study demonstrates that hydrogen bonding with water molecules can modulate the photoreaction pathways through charge stabilization of the phenolic oxygen to also make the β-twist pathway accessible.
As anticipated above, time-resolved X-ray crystallography studies, which provide information on the PYP light-induced dynamics, were employed to resolve the controversy over the PYP isomerization mechanism. 141,143,148 Such studies 143,148 suggested that there are two possible pathways for the photoisomerization. These are the HT mechanism suggested by Liu and Asato 149 and the BCP mechanism 22 originally proposed for the rPSB chromophore in rhodopsin. The bonds involved in the HT mechanism in PYP are the C9−C10 and C10C11 bonds. Instead, the BCP path is characterized by simultaneous C10C11 and C12−S13 twisting (note that C12−S13 has a partial double-bond character because it is a thioamide) and is the only mechanism observed in the timeresolved X-ray study by Pande et al. 141 The HT and BCP mechanisms have been investigated in PYP by computing the corresponding relaxed scans. Unlike animal rhodopsins, which feature a barrierless relaxation with no S 1 intermediate, here a shallow S 1 minimum is located for PYP with a C9−C10−C11−C12 dihedral angle of 155.5°. 140 This is consistent with a less reactive potential energy surface and therefore with the slower (550 fs) photoisomerization compared to the animal rhodopsin (e.g., Rh). However, it is comparable to the microbial rhodopsins with a rPSBAT chromophore (e.g., bR) displaying a flatter S 1 potential energy surface (see section 4.3.3 above). In fact, as reported for bR, the S 1 relaxation continues downhill after the minimum and reaches the CI. The energy profile of the computed path is consistent with the low fluorescence quantum yield observed in PYP, 150,151 which can be assigned to the located S 1 minimum.
Inspection of the computed S 1 relaxation coordinate indicates that, even though the hydrogen bond between the phenolate oxygen of the chromophore and the Glu46 residue is weakened by the decreased negative charge, it still helps to restrain the phenyl ring from twisting along the C9−C10 bond (i.e., the one that would be involved in the HT mechanism). This C9−C10 twisting goes from −3.6°(at the FC point) to 3°(at the S 1 minimum) to 33°(at the CI). In order to continue such deformation on S 0 , the chromophore has to overcome a 7 kcal/ mol barrier to break the hydrogen bond between the chromophore and Glu46. On the other hand, along the S 1 relaxation following a BCP mechanism, the twisting of C10 C11 goes from 170°(at the FC point) to the 155°(at the S 1 minimum) and finally to 85°(at the CI). Simultaneously, the twisting of the C12−S13 bond goes, as expected for such mechanism, in the opposite direction from −177°(at the FC point) to −168°(at the S 1 minimum) to −147°(at the CI). Upon relaxation to S 0 , the twisting about the C12−S13 bond has to overcome a smaller S 0 barrier (4.4 kcal/mol) with respect to the HT and does not require the breaking of the hydrogen bond with Glu46. This suggests that the BCP mechanism is more favorable than the HT mechanism, even though both mechanisms are feasible. The smaller barrier for the BCP mechanism (see Figure 29) shows consistency with the shorter lifetime determined experimentally. The BCP mechanism leads to the formation of I CT (∼1.7 ns) while the HT mechanism leads to the formation of pR 1 (lifetime = ∼3 ns). 143 The conclusions above are based on relaxed scans, which only provide an approximate description of the reaction coordinate. An improved description can be achieved with a minimum energy path or trajectory studies that incorporates the effects of the light-induced dynamics, as discussed below. In order to investigate the photoinduced dynamics of the isomerization of the PYP chromophore, trajectory studies have been carried out in the gas-phase, 137,139 in solution, 138 as well as in the protein. 136 Semiclassical trajectory calculations by Groenhof and co-workers 139 with nonadiabatic transitions are described by Robb's surface hop method (see section I.3) followed the progression of gas-phase pCK − on the S 1 potential energy surface. These molecular dynamics simulations in the gas phase confirm the result of reaction-path calculations. Indeed, it is found that, in the gas phase, the trajectories reach the S 1 minimum along the α-twist coordinate. It is also found that the S 1 /S 0 IS lies far away from the S 1 minimum, making the S 1 decay inefficient. Again, this demonstrates the importance of the environment in controlling the mechanism and time scale of the S 1 decay of the pCA chromophore.
To examine the effect of the hydrogen bonding on the PYP chromophore, a simple water cluster model was employed to perform QM/MM trajectory calculations. 139 The QM subsystem was treated at the CASSCF level. The results demonstrate that the predominant S 1 decay channel involves an α-twist (88%) rather than a β-twist (22%) mechanism. However, this is an important change from the gas-phase where the decay is exclusively via the α-twist. This is explained by the fact that the water molecules hydrogen bonded to the carbonyl oxygen would stabilize the α-twisted minimum, while those hydrogen bonded to the phenolate oxygen would stabilize the β-twisted minimum, thus modulating the excited-state pathways. The transient stabilization of the charge distribution would bring the IS close to these minima, causing a faster decay. A faster decay in solution has also been documented by Martínez and co-workers using ab initio multiple spawning (AIMS) calculations. 34 The results, which are consistent with the CC2 reaction path studies 138,147 mentioned in section 5.2 above, also indicate that the S 1 dynamics and decay time scale are modulated by hydrogen bonding. Naseem et al. 152 further studied the PYP chromophore in solution and found that there is an increase in pK a in the cis form but not in the trans form due to electrostatic repulsion between the coumaryl tail and the phenolate ring moieties. This aids the protonation of the cis form in the next step of the photocycle.
In order to understand the effect of the PYP protein environment, Groenhof et al. 136 carried out QM/MM trajectory computations at the CASSCF level with Robb's surface hopping algorithm. They showed that the S 1 pCA chromophore rapidly decays to S 0 through a β-twist motion instead of the α-twist motion as reported in the gas-phase. The authors have attributed it to hydrogen bonding between the phenolate ion with the nearby Tyr42 and Glu46 residues. Similar to the effects in the water environment discussed above, these hydrogen bonds enhance S 1 decay from the β-twisted minimum. Moreover, these simulations also demonstrated that the same hydrogen bonds remain intact along the isomerization coordinate, which follows a NBCP mechanism consistently with the reaction path studies. At the surface hop, there is a ca. 90°twist at the C10C11 bond while the adjacent C12−S13 bond is twisted by ca. 50°. The simulations show a 200 fs fluorescence lifetime and 30% isomerization quantum yield and are therefore close to the corresponding experimental values (550 fs 141 and 35% quantum yield 151 ).
In contrast to the solution environment, the protein in PYP causes the chromophore to isomerize exclusively along the β-twist coordinate, with none of the excited state population undergoing an α-twist. 136 This indicates that, apart from hydrogen bonding to the phenolate part of the chromophore, there is also a large effect of the protein residues on tuning the isomerization mechanism of the chromophore. In particular, the positively charged guanidinium moiety of Arg52 plays an important role in biasing the relaxation toward double-bond (β) isomerization. Other studies 34 have confirmed that when a small positive point charge (+0.2) is placed in the position where the nitrogen atom (Nε in Figure 26) of the Arg52 residue is, the α-twist motion is suppressed while the β-twist is promoted, while a negative point charge would produce the opposite effect. 34 Such an effect has been experimentally assessed by producing the Arg52Gln mutant. It was shown that Arg52Gln can still enter the photocycle, although with a lower rate and quantum yield. 153 The conclusion is that Arg52 is important but not essential for photoactivation.
In an attempt to simulate the effects described above, QM/ MM trajectory calculations of the PYP Arg52Gln mutant were performed by Groenhof et al. 154 They find that without a positively charged Arg52, the predominant S 1 relaxation follows an α-twist rather than a β-twist coordinate. It was also found that in the Arg52Gln mutant, the hydrogen bond between the carbonyl oxygen of the chromophore and the backbone amino group of Cys69 is broken. However, the chromophore rapidly forms three new hydrogen bonds involving the amino groups Tyr98 and Asp97 as well as a water molecule that enters the chromophore pocket from the solution. These hydrogen bonds stabilize the negative charge on the alkene moiety and lead to a rapid decay of the chromophore to S 0 . In conclusion, the Arg52Gln mutant of PYP has a similar decay mechanism as the chromophore in water. Even though the α-twist path does not produce the cis diastereoisomer of the chromophore, the computations show that it causes variations in the hydrogen bond network and the flipping of the thioester group. However, the same processes were also computationally demonstrated to occur during the β-twist relaxation. It has thus been proposed 154 that the flipping of the thioester group and not the β-twist is the key step in the photocycle activation. This would be consistent with the fact that, as it is experimentally observed, the PYP Arg52Gln mutant is still photoactive. 155

Controversies Regarding Low-Barrier Hydrogen Bonding and Photoisomerization Mechanisms
All the PYP mechanistic studies discussed so far are based on models featuring with protonated Arg52 residue and an anionic chromophore. However, an alternative scenario has been debated where Arg52 is deprotonated and the anionic chromophore is protonated to become a neutral one. More specifically, this would require a low-barrier hydrogen bonding (LBHB) in PYP. LBHB is a nonstandard hydrogen bond network, which was originally proposed to possess covalent bondlike character. 156 As opposed to standard hydrogen bonding, where the hydrogen is covalently bonded to one heteroatom and more weakly interacting with the other, LBHB is a stronger intermolecular interaction where the hydrogen is more equally shared between two heteroatoms. LBHB between the Arg and the chromophore ion pair would cancel out, in part, the negative charge of the chromophore and the positive charge on Arg52, yielding a more neutral charge environment. This controversy was triggered by a neutron diffraction study by Yamaguchi et al., 157 suggesting the existence of a neutral, LBHBbounded chromophore with a neutral Arg52 residue in PYP. In

Chemical Reviews
Review DOI: 10.1021/acs.chemrev.7b00177 Chem. Rev. XXXX, XXX, XXX−XXX fact, the existence of LBHB allows the stabilization of the anionic pCA chromophore without the need of a charged Arg52. However, Saito et al. 158 suggested that there is no LBHB in PYP based on QM/MM calculations and electrostatic pK a calculations. Later, Nadal-Ferret et al. 159 suggested that PYP displays LBHB in solid phase, while it reverts back to a regular hydrogen bonding in solution. Further studies by Graen et al. 160 suggested that a single PYP model in vacuum is not necessarily a good model to describe the experimental results measured in solid phase as they found that a model with multiple copies of PYP, which would better describe the PYP in crystals, show no LBHB. The discrepancy is caused by the fact that the real environment around PYP is sufficiently isotropic due to the high ionic strength in crystals but a single PYP in vacuum is anisotropic, which is not a good model to reflect the solid phase of PYP. Moreover, this study has re-examined the neutron densities found by Yamaguchi et al. 157 and suggested that the large differences in experimental neutron densities at the side chain of Arg52 do not exclude the possibility of a protonated Arg52. Even though most of the studies point toward a fully protonated Arg52 with no LBHB, the effect of LBHB versus deprotonated Arg52 on the S 1 dynamics of PYP has yet to be studied computationally. This is particularly interesting as it has also been proposed, on the basis of ultrafast stimulated Raman studies, that the initial S 1 relaxation involves conversion from LBHB to a normal hydrogen bonding before the photoisomerization occurs. 161,162 The controversy over the photoisomerization mechanism of PYP is also still unresolved. The study by Jung et al. 143 on the picosecond time scale has suggested that the photoisomerization of the double bond can proceed via two structurally distinct cis intermediates through both the HT and BCP mechanisms. However, studies by Pande et al. 141 at the femtosecond resolution suggest that the photoisomerization proceeds exclusively through the BCP mechanism, which is consistent with the simulations from Groenhof et al. 136 Both HT and BCP mechanisms would ultimately lead to pB (i.e., the cis distereoisomer of the chromophore) and generate the signaling state of PYP, even though the BCP mechanism would be the major pathway. However, a recent experimental study 163 of a different PYP from Leptospira biflexa (Lbif) suggested that there is a bifurcation in the photocycle of PYP. Lbif PYP shows an apoprotein different from the most widely studied PYP from Halorhodospira halophile (Hhal), even if it has the same pCA chromophore. The study suggests that the longstanding inconsistency in the existence of a bifurcation in photocycle at room temperature is due to a small population of intermediates only found in the most computationally studied Hhal PYP. These results will likely motivate computational studies of Lbif PYP to see whether any bifurcation exists leading to different isomerization mechanisms.

Summary: Possible Bifurcation Photocycle
The results of different simulations of the PYP photoisomerization are, in general, consistent with the experimental observations. They suggest an S 1 BCP relaxation mechanism leading to a CI. Even though the HT mechanism is not observed in trajectory simulations inside the protein, the existence of a bifurcation cannot presently be ruled out, especially with the recent experimental finding on a different apoprotein that could lead to a detour in the photocycle.
Unlike rhodopsins, PYP has an anionic chromophore with a relatively small apoprotein. This has led to the controversies regarding low-barrier hydrogen bonding. Due to its small protein size, the chormophore is partially exposed to the solvent. This would make a single PYP in vacuum not a good model to reflect the solid phase of PYP. These findings suggest that it is important to have a model that has a good description of the chromophore environment that is consistent with the experimental conditions in order to be reliable.
At any rate, the mechanism that controls the primary steps of the photocycle of PYP is a double-bond photoisomerization, regardless of whether the HT or BCP mechanism is followed. Then, it undergoes a breaking of a hydrogen bond between the chromophore and protein backbone that triggers a proton transfer from the protein to the phenolate oxygen of the chromophore. This leads back to the pB structure and the signaling state of PYP. Since both HT and BCP would lead to the same pB, the BCP mechanism must follow an aborted isomerization, just as in Rh.
PYP is unique among natural photoreceptors because it forms large and stable crystals that are suitable for X-ray diffraction. This system has therefore provided an excellent opportunity not only to simulate the time-resolved isomerization computationally but also to "see" it experimentally with time-resolved X-ray crystallography. PYP therefore has led the way for the application of time-resolved X-ray crystallography in photobiology and has presented a unique and exciting opportunity to compare QM/ MM simulations and X-ray crystallography directly. Despite this progress, there is still more to be done to better understand the details of the isomerization mechanism of the PYP chromophore and the formation of the signaling state. Phytochromes are a family of photoreceptor proteins that bind a linear tetrapyrrole as a chromophore. They are characteristic among photoreceptors for their absorption in the red, which leads to a photoproduct with absorption in the far-red edge of the visible spectrum. Phytochromes are photochromic and can be switched between the red absorbing (P r ) and far-red absorbing (P fr ) states. While originally discovered in higher plants where they play an important role in regulating circadian rhythm, members of the phytochrome family of proteins have also been found in bacteria, cyanobacteria, and fungi where they regulate a variety of biological functions. 164 An example of a prototypal plant phytochrome is shown in Figure 30A. Not all phytochromes share the same architecture, although the photosensory module of most canonical phytochromes consists of three conserved domains: PAS, GAF, and PHY domains. 165 The chromophores of phytochromes are methine-bridged linear tetrapyrroles (known as bilins) that are bound to a cysteine residue of the apoprotein through a thioether group, as shown in Figure 30B. Phytochromobilin (PΦB) acts as the chromophore of plant phytochromes, while in bacterial and fungal phytochromes the chromophore is phycocyanobilin (PCB) or biliverdin (BV). Despite these small differences in the chromophore structures, spectroscopic studies have indicated that phytochromes of plants, bacteria, and fungi all share a common photocycle, shown in Figure 31. In the dark-adapted state, most phytochromes exist in the P r state. Resonance Raman (RR) studies 166,167 indicate that photoexcitation of chromophore at ca. 660 nm initiates a double-bond isomerization of the methine double bond that bridges rings C and D (the C15C16 bond in Figure 30C) series of thermal relaxations of both the chromophore and the protein, the P fr active state is reached, which is called "active" because it triggers the biological function. The protein can then be reverted back to the P r state by exciting in the far red (at ca. 730 nm) to give Lumi-F (although this step may also occur thermally on a much longer time-scale). Back-conversion does not follow the same mechanism as the forward process but rather goes through different intermediates, Lumi-F and Meta-F. The details of the photoinduced Z to E isomerization from P r to Lumi-R, as well as the E to Z isomerization from P fr to Lumi-F, have slowly begun to emerge over the past two decades. The understanding of the photochemistry of phytochromes has been largely driven by experimental spectroscopic studies employing RR, 166,167,[170][171][172][173]174,175 as well as, more recently, from the analysis of the X-ray crystal structures 168,176−178 of the P r and P fr intermediates. Many of these studies have been largely aided by theory as well. For example, theory has been instrumental in the assignment of RR spectra 166,172,173 and in the interpretation of polarization resolved femtosecond pump−probe spectroscopy. 179,180 However, in contrast to the case of rhodopsins and PYP, there are very few computational studies attempting to determine the photoisomerization mechanism of phytochromes (e.g., by optimizing excited state reaction pathways or by calculating semiclassical trajectories). As we have done in previous sections, we will focus exclusively on such studies aimed at modeling the isomerization process in phytochromes. However, due to the limited number of studies, this section will not be divided into subsections related to funnels, paths, and dynamics as we have done in previous sections. There are a number of reasons why such computational studies have been relatively scarce: (1) the large size of the chromophore prevents the use of many of the quantum chemical methods discussed in section I. For example, a CASSCF calculation that includes the full π-electron system in the active space would require an active space of 28 electrons in 25 orbitals, which is not tractable due to the factorial scaling of computational cost with the size of the active space. For this reason, electronic structure methods used to study phytochromes are less computationally demanding and, therefore, often more approximate. For instance, this often means using single reference methods like CIS or TD-DFT that are unable to treat the nonadiabatic regions of potential energy surfaces correctly 181,182 or using approximations that do not account for static and dynamical electron correlations in a balanced way (as explained in section I.1). It is therefore difficult to achieve a solid description of the mechanistically important  9-syn, and 14-anti) conformation. Note that the three chromophores differ only at rings A and D. Only one resonance structure is shown. (C) Abridged chromophore models commonly used in gas-phase studies of PΦB, PCB, or BV. The reactive bond is indicated with an orange curved arrow. Important atoms are numbered to indicate bonds that are free to photoisomerize in the gas phase. Figure 31. Phytochrome photocycle. The first step of the photocycle, the Z to E isomerization of Pr to Lumi-R, occurs on a picoseconds time scale and exhibits multiexponential kinetics with both an ultrafast (2−3 ps) and 1 order of magnitude slower (30−40 ps) component. 169 reaction paths and CIs, for instance, found on the excited-state potential energy surface. (2) Crystal structures of phytochromes have only become available in the past decade. The first crystal structure of a phytochrome was published in 2005 for the bacterium Deinococcus radiodurans, 176 while the first crystal structure of a plant phytochrome was only resolved recently in 2014. 168 This means that early investigations of the photoisomerization of the phytochrome chromophore were confined to the gas-phase or to an implicit solvation model such as the polarizable continuum model (PCM) as an approximation to account for the environmental effect of the protein. 183−185 An added difficulty was that, for a long time, even the conformations of the chromophores in P r and P fr were unknown due to lack of a crystal structure. While there has been, for the most part, a general consensus that the photoisomerization occurs at C15 C16, it was not clear if a single-bond isomerization (e.g., a Hulatwist mechanism similar to the one introduced in section 4.3.1) would have occurred simultaneously or in tandem. In fact, before crystal structures of plant phytochromes became available, it was predicted that PΦB is in the 5-Z,anti-10-Z,syn-15-Z,anti (ZZZasa) configuration on the basis of RR studies interpreted using DFT gas-phase frequency calculations. 172 This conclusion was also supported by SAC−CI calculations, 186 leading several early computational studies to adopt the ZZZasa configuration for computational studies. SAC−CI, short for "symmetry adapted cluster-configuration interaction", is a method that uses an exponential ansatz for the wave function parametrization similar to the EOM-CC (see section I.1), but it uses a different exponential operator than the one used in a canonical coupled cluster. 187−189 It was only confirmed recently that PΦB in plant phytochromes is in the distinct ZZZssa configuration and therefore different from the one originally hypothesized. 168 (3) It has recently been found that, in several phytochromes, both the P r and P fr states do not exist in one homogeneous conformation but in two conformers (or substates) that quickly interconvert (e.g., see ref 190 and references therein). This complicates their modeling, particularly because the two substates are expected to have different photophysics/photochemistry (e.g., see ref 180), and therefore both need to be modeled to achieve a comprehensive understanding.
Unlike rhodopsins and PYP, phytochromes also have a more complex multidomain protein structure, potentially making QM/MM studies difficult. For example, as displayed in Figure  30A, a plant phytochrome monomer typically comprises three domains: GAF, PAS, and PHY. Rhodopsins and PYP are more compact.
Despite the difficulties outlined above, there have been a number of attempts to computationally investigate the photoisomerization mechanism of phytochromes. Most of these studies have been confined to gas-phase chromophore models, which serve as preliminary studies to understand the excited-state behavior of the chromophore in the absence of any external perturbation. This approach is similar, historically, to the case of rhodopsins, where studies were first conducted with gas-phase calculations on reduced models of rPSB, 20,191,192 only to evolve to full QM/MM simulations later on. 21,193 Early studies by Durbeej, Borg, and Eriksson 194−196 of the phytochrome chromophore used a slightly reduced model of PΦB, where they replaced the C3 thioether linkage, the propionic carboxyl groups at C8 and C12, and the methyl groups at C2, C7, C13, and C17 by hydrogen atoms (see Figure  30C) to compute a relaxed scan around the C15C16 isomerizing bond. They used the CIS method to optimize the S 1 state at different constrained C15C16 dihedral angles and TD-DFT to compute the energies of the S 0 and S 1 states along these scans. Their studies indicated that the PΦB has a significantly lower rotation barrier when it isomerizes exclusively along the C15C16 bond, as opposed to a Hula-twist mechanism where an adjacent single bond also isomerizes. 194 However, since such studies investigated the gas-phase chromophore, it is unclear if the same mechanism will apply in the protein where there are steric constraints restricting largeamplitude motions of the chromophore. Durbeej et al. also found that deprotonation of the chromophore at the pyrole nitrogen of ring B or ring C leads to a higher rotation barrier than if the chromophore remains protonated, 195 which indicated that the chromophore does not become deprotonated before isomerizing.
Later, Durbeej 197 used the same reduced model to compare the S 1 isomerization barriers around the C4C5, C10C11, and C15C16 bonds (i.e., around the methine double-bonds bridging rings A−B, B−C, and C−D, respectively). This time, he used TD-DFT instead of CIS for S 1 optimizations. He also tried modeling the chromophore both in the gas-phase and in PCM. He found that the second excited state, S 2 , which has doubly excited character with respect to S 0 , is not accessed along the isomerization coordinate in either the gas phase or PCM models. More interestingly, he found that the chromophore prefers to isomerize around the central C10C11 bond in the gas phase, while the isomerizations of the terminal double-bonds C4C5 and C15C16 have higher barriers in both gas phase and PCM. Since protein-bound PΦB isomerizes selectively along the C15C16 bond, 166,167 this indicates that the protein plays an important role in this selectivity. This is in line with the protein crystal structures, which indicate that the protein cavity forms a tight binding pocket around rings A, B, and C, while ring D is less tightly bound. The gas-phase results were later extended to PCB and BV, 198 indicating that the conclusions above are valid in bacteria, cyanobacteria, and fungal phytochromes as well. A few years later, Zhuang, Wang, and Lan 199 performed nonadiabatic dynamics simulations with on-the-fly surface hopping using the same reduced model as Durbeej to investigate the reactivity of the PΦB chromophore in the gas phase. The dynamics were run with the semiempirical OM2/MRCI method, but single-point calculations at the ADC(2)/TZVP level were used for verification. They found, in agreement with Durbeej,197 that isomerization occurred exclusively around the methine bridge linking rings B and C, with no trajectories isomerizing around C4C5 or C15C16. However, rather than isomerizing around C10C11, most trajectories actually isomerized around the C9−C10 bond with only a few trajectories isomerizing around C10C11. Recently, Falklof and Durbeej 200 investigated the isomerization of protein-embedded BV in a bacteriophytochrome using a full QM/MM model. They again performed a relaxed scan around the C4C5, C10C11, and C15C16 torsions of the BV chromophore inside the protein.
Their calculations confirmed that steric interactions with the protein pocket are the reason for the selectivity of the photoisomerization around C15C16 because they restrict the isomerization at the A−B and B−C bridges. This is the only QM/MM isomerization study of phytochromes, so there is still very limited understanding of the atomic-level protein interactions controlling selectivity and photoisomerization mechanism.
Note that, unlike rPSB, excitation of the PΦB chromophore does not immediately lead to a significant change in dipole Chemical Reviews Review DOI: 10.1021/acs.chemrev.7b00177 Chem. Rev. XXXX, XXX, XXX−XXX moment, 195 meaning that no significant charge transfer occurs upon vertical excitation. The excitation involves a HOMO to LUMO transition, with both orbitals localized on the central B and C rings. 201 The reason for this is most likely that the A and D rings are somewhat twisted (by ca. 45°) 201 with respect to the B and C rings, decreasing the degree of delocalization of the frontier orbitals along the full π-conjugated framework (i.e., forming a locally excited state). However, along the C15C16 isomerization path, a large and sudden change in dipole moment does occur when the chromophore is ca. 60°twisted out of plane, indicating that a charge-transfer character has to develop along the S 1 path before decay. 195 This is similar to protonated Schiff bases, where twisting localizes the HOMO and LUMO π orbitals on either side of the isomerizing bond, inducing a transfer of the positive charge from the Schiff base moiety to the opposite molecular end (see Figure 9A). As mentioned in the caption of Figure 31, pump−probe experiments on BV-binding phytochromes have shown that they exhibit multiexponential decay dynamics, with both fast (2−3 ps) and slow (30−40 ps) components. 169 In 2009, Altoèand coworkers 201 subsequently published a computational study where they attempted to explain the observed dynamics by using a gasphase reduced model of the chromophore. They used the CIS electronic structure method for S 1 optimizations, but reaction paths were computed with CASPT2//CASSCF with a 10electron 11-orbital active space. The reduced (10,11) active space was chosen based on the convergence of the vertical spectroscopic properties upon a systematic increase. They found that the gas-phase chromophore, upon excitation, initially relaxes from the Franck−Condon region to a shallow S 1 minimum. From that minimum, they performed a relaxed scan along four different rotation coordinates: C4C5, C9−C10, C10C11, and C15C16 and found that isomerization of the C9−C10, C10C11, and C15C16 bonds are competitive because they have similar S 1 barriers, while C4C5 is the only bond that is not likely to isomerize. Their calculations also indicate that isomerization at C15C16 does not lead to a CI but rather to a twisted intramolecular charge-transfer (TICT) state that has an S 0 -S 1 energy gap of 17.4 kcal/mol at the CASPT2//CIS level (or, when they attempted to reoptimize it at the CASSCF(10,11) level of theory, an energy gap of 5.4 kcal/mol with CASPT2// CASSCF). This twisted intermediate has its positive charge localized on the D ring, as opposed to the vertical excited state that had its charge localized still on the B and C rings. Such a TICT state would then have to decay to the ground state via a nearby sloped CI (see Figure 32). On the other hand, isomerization around the C9−C10 or C10C11 bonds leads directly to a peaked CI.
On the basis of these observations, the authors 201 hypothesized that the isomerization of the C9−C10 or C10C11 occurs on an ultrafast time scale, and therefore these processes correspond to the fast 3 ps component of the pump−probe experiments (gray panel in Figure 33). Even if the protein restrains the isomerization motion, they propose that deactivation does not lead to an isomerized product but returns back to the original P r reactant. The other path, which is the biologically relevant C15C16 isomerization, leads to a TICT that is longer lived and decays via a sloped CI which is higher in energy with respect to the minimum (see Figure 32 right). Thus, the isomerization through this channel occurs on a slower time scale. The authors assign this path to the experimentally observed slow ∼30 ps component, consistently with pump−probe experiments, indicating that the slow component corresponds to the biologically relevant C15C16 isomerization (white panel in Figure 33). 169

Summary: Protein Environment Selects the Isomerizing C15C16 Bond
Among photoreceptors that undergo double-bond isomerization, phytochromes and other related bilin-binding photoreceptors have been the least studied so far. There are several reasons for this, but chief among these reasons is the size of the chromophore (making the application of QM methods expensive and, consequently, limiting the understanding of its electronic structure) and the scarse availability of experimental data such as crystal structures until somewhat recently. While QM/MM studies on phytochromes are slowly starting to emerge, QM studies on chromophore models have paved the way to understanding the photochemistry and photophysics of the chromophore.
One of the main results emerging from theoretical simulations is that the chromophore in the gas phase isomerizes preferentially around the central C9−C10 or C10C11 bonds, but the phytochrome protein environment modulates the selectivity, most likely through steric interactions, and causes the isomerization to occur along C15C16 instead. Experimentally, pump−probe experiments on BV-binding phytochromes have shown that they exhibit multiexponential decay dynamics with a slow (30−40 ps) and fast (3 ps) component, so computational studies have also posited that this could be due to competing C15C16 and C9−C10/C10C11 isomerization in the protein.
With new emerging QM methods and more efficient algorithms, computational studies on phytochromes are almost certain to increase in the near future. Because of the size and extent of conjugation of their tetrapyrrole backbone, bilins absorb naturally in the red and can emit in the far-red as well. The BV chromophore is also naturally occurring in mammalian cells due to breakdown of heme. 202 This has important implications for their potential use in optogenetics as well as bioimaging, since mammalian tissue is more transparent at far-red and near-IR wavelengths, the development of optogenetics tools or fluorescent proteins that absorb/emit in the far-red or infrared is highly desirable for deep-tissue processes. Recent studies have demonstrated that this is a practical and promising possibility. 203,204 This is almost certain to further catalyze the use of computer simulations to understand the photoisomerization of BV and related chromophores in different proteins. Figure 32. A TICT is a minimum on the excited-state potential energy surface reached by following a torsional deformation accompanied by a charge transfer. A TICT appears in a vicinity of a sloped conical intersection (right) instead of a peaked conical intersection (left) and can decay to the ground state either by fluorescence or by accessing the nearby sloped CI.

Rhodopsin Mimics
To systematically study the correlation between structure, absorption wavelength, and photochemical and thermal reactivity in retinal proteins, the intracellular lipid binding protein family (iLBP), and in particular the human cellular retinol binding protein II (hCRBPII) and human cellular retinoic acid binding protein II (CRABPII), have been employed as targets for protein engineering. 54,205−207 Tolerance to mutations and the large binding cavity of iLBPs allowed flexibility in the redesign of the protein and enabled the binding of a range of different ligands. 208 More specifically, the relative rigidity of the iLBP backbone structure makes it an ideal scaffold for protein engineering because it is more likely that mutations will not result in structural rearrangement, 208 allowing the effects of specific mutations to be compared directly. These systems have been engineered to bind the rPSB chromophore such that it is almost fully embedded within the binding pocket. This is an important property of these systems because otherwise exposure to aqueous media would buffer electrostatic changes and result in the loss of effectiveness of the protein in modulating the photophysics of the embedded chromophore. Furthermore, hCRBPII and CRABPII are water-soluble proteins instead of membrane proteins and thus much easier to crystallize.
Over the past few decades, Borhan, Geiger, and co-workers have prepared and investigated several iLBP-based rhodopsin mimics. These authors exploited the ease of crystallization of these engineered protein and, consequently, have access to detailed structural information via crystallographic methods. Although, most of their reported studies focus on spectroscopy and color tuning, 54 one recent report 209 presented a reengineered CRABPII system that displays both photochemical and thermal isomerization of the embedded retinal chromophore (see Figure 34A). Before this study, it had been difficult to demonstrate that photoisomerization of the protein-bound retinylidene chromophore (usually rPSBAT) occurs in these types of engineered systems. More specifically, Borhan, Geiger and co-workers have reported on an iLBP rhodopsin mimic where the chromophore specifically isomerizes at the protonated Schiff base −CHNH− double bond both thermally and photochemically. This isomerization was characterized at atomic resolution by quantitatively interconverting the isomers in the crystal (in the solid state), both thermally and photochemically. 209 This event is accompanied by a large pK a change of the imine similar to the pK a changes observed in bacteriorhodopsin and visual opsins during isomerization.  Although there is presently no computational study on the type of isomerizing rhodopsin mimics described above, different groups have attempted to reproduce the observed change in absorption maximum upon mutations of an hCRBPII system. 55,210 In one case, the mechanism of the photoisomerization process has also been investigated 211 in a CRABPII mimic, but this is in a different system from the one where the isomerization has been detected experimentally. In this study, the authors investigated the viability of rPSBAT photoisomerization from all-trans to 11-cis and to 13-cis. Since the protein pocket of CRABPII is adapted to bind the rPSBAT chromophore (i.e., the all-trans isomer), the computed isomerized products were predicted to be unstable. Indeed, the optimized geometries of rPSB11 and rPSB13 are significantly distorted, with dihedrals −27°and −9°twisted out of plane, respectively. These deformations extend to other parts of the chromophore backbone. From Figure 34B, it can be seen that the isomerization leads, in both cases, to large skeletal changes that mainly involve the N-containing segment of the chromophore. This happens even though the β-ionone ring sticks out of the protein pocket and would be, in principle, allowed to rotate.
Further computational studies in this area promise to reveal how to construct rhodopsin mimics undergoing photoisomerization and, hopefully, showing photochromism and bistability. This could be done, for instance, by constructing and studying comparative models of rhodopsin mimics prepared in the lab and shown to be either photochemically reluctant or active like the one of ref 209.

Fluorescent Proteins
Double-bond isomerization in photobiological systems is by no means exclusive to rhodopsins, phytochromes, and the PYP chromophore. Fluorescent proteins (FPs) are another example of proteins where photoinduced double-bond isomerization of the chromophore can take place. Photoinduced processes, including double-bond isomerization, in FPs have been reviewed recently by Acharya et al. 212 FPs have a rich photophysics where several excited-state processes may occur. This includes isomerization, 213a proton transfer, 213b and even electron transfer. 213c,d The interplay of these excited-state processes is a subject of many theoretical and experimental studies in different types of FPs. 212 FPs have revolutionized the field of bioimaging, 214 and since the discovery of the Green Fluorescent Protein (GFP), 215 a palette of FPs with colors spanning almost the entire visible spectrum have been prepared by mutagenesis. 216 However, although the fluorescent properties of FPs have been widely exploited in biotechnologies, the biological function of wild-type FPs is poorly understood. Therefore, whereas double-bond isomerization triggers biological responses in rhodopsins, phytochromes, and PYP, double-bond isomerization in FPs is a secondary and often undesirable process. This can be true, for instance, when it comes to bioimaging applications, because double-bond isomerization competes with fluorescence and may lead to quenching. On the other hand, it has been shown that double-bond photoisomerization in FPs can be exploited for specific purposes, such as in the development of a new generation of photoswitchable FPs such as Dronpa (see Figure 35).
The majority of FP chromophores are made up of an imidazolinone ring conjugated with another aromatic ring via a methine bridge. In GFP and in Dronpa, the aromatic ring is phenolate ( Figure 35B), which is reminiscent of the phenolate already seen in the PYP chromophore.
The methine bridge in the GFP chromophore is the only segment of the π-conjugated framework that is not locked in a ring and can therefore rotate around the single bond (ϕ) or double bond (τ). Rotation along either ϕ or τ usually leads to a twisted excited-state intermediate (TI ϕ and TI τ ). 219,220 In early computational studies of a GFP model chromophore, a minimum energy CI was found along a hula twist coordinate (MECI hula ) that involves the out-of-phase twisting of both ϕ and τ, but in the gas-phase, it lies high in energy relative to TI ϕ and TI τ . More recent studies have found low-lying CIs in the gas phase that are twisted only along ϕ 221,222 and τ 221 (MECI ϕ and MECI τ , respectively). The relative energetics and geometries of such intermediates and minimum energy CIs has been the subject of numerous computational studies. [219][220][221]223,224 However, it is clear that the solvent environment, 220,223,224 protein environment, protonation state of the chromophore, 222 as well as chromophore substituents all have a substantial impact on the relative energetics of these CI points and the pathways leading to them (TI ϕ and TI τ , see Figure 36). Therefore, the photophysics of the chromophore is highly dependent on such factors.
Twisting around either ϕ or τ is a major decay pathway in many FPs, leading to twisted S 1 intermediates (TI ϕ and TI τ ) with nearby low-lying CIs. This is partly also responsible for phenomena such as emission intensity variations, on−off blinking, and photoreversible switching between light and dark states. 225−228 When the chromophore is not held rigidly in the protein binding pocket (e.g., in gas-phase or solution), the barriers for rotation (TS ϕ and TS τ ) are small, which explains why fluorescence quantum yields of GFP chromophore models in vacuo or in solution are several orders of magnitude smaller than in the protein. 229−231 Conversely, restricting the motion of the GFP chromophore by placing it in a rigid environment increases the barriers of rotations (TS ϕ and TS τ ) and leads to brighter fluorescence. 232,233 Steric interaction is not the only effect at play here, since hydrogen-bonding and electrostatic interactions with the proteins likely play a very important role in influencing the ease of rotation around ϕ or τ.
In the vast majority of FPs, the chromophore is in the cis conformation. 234 On the other hand, the trans conformer is usually not fluorescent, with a few notable exceptions. This possibly stems from the fact that the trans form of the chromophore is tilted or twisted along ϕ or τ to relieve a steric clash between the imidazolinone carbonyl oxygen and the phenol ring, making the trans form more conformationally flexible than the cis form. 235 The presence of a long-lived transient dark trans state can be detrimental to the functionality of FPs in singlemolecule studies, 236 because it contributes (or is a major source of) blinking and intensity variations. However, the same phenomenon can be exploited in photoswitcheable FPs if it can be controlled. One notable example, mentioned above, is Dronpa, 237 where the native state is bright (cis form) but can be switched off (to the trans isomer) by illuminating with intense blue light or switched back on again with weak near-UV light. Dronpa is an example of a negative photoswitch because it is bright in its native state and can be turned off with light. There are also several examples of positive photoswitches, such as Padron 238 and kindling FPs, that are dark in their native state (typically, in the trans conformation) but their fluorescence can be turned on momentarily by photoisomerization to the light state. 238 The exact mechanism of photoswitching in such photoswitchable FPs is a subject of debate, 212 in particular with regards to the protonation states of the chromophore. This has complicated the computational investigation of isomerization in such systems. However, the photoswitching occurs on an ultrafast time scale, 239 which indicates that a low-lying CI must be accessed during photoisomerization. For example, a study by Morokuma and co-workers 240 employing a CASSCF/AMBER model of Dronpa found that, regardless of the protonation state of the chromophore, isomerization occurs via a one-bond flip around τ rather than a Hula-twist CI. On the other hand, a more recent study by Morozov and Groenhof 241 found qualitatively different results. They investigated a Dronpa mutant using 12 QM/MM semiclassical trajectories computed at the SA(2)-CASSCF(6,6)/3-21G/AMBER level of theory and found that while some of the trajectories (8 out of 12) did not isomerize at all in 50 ps, four of the trajectories decayed via a Hula-twist mechanism involving twisting of both ϕ and τ. Such a Hula-twist CI was also optimized by the same authors at the XMCQDPT2 level of theory. 241

Twisting about a Double Bond As a Deactivation Channel in Biomimetic or Biological Ring-Locked Double Bonds
As discussed in several of the sections above, twisting around a double bond does not always lead to an isomerized product. Since isomerization quantum yield is never 100%, there are nonreactive (or aborted) isomerization events that lead back to the original reactant. Such aborted isomerization events contribute to internal conversion, since they simply depopulate the excited state at a twisted CI only to regenerate the original reactant.
Sometimes, such an aborted isomerization is the only possibility. For instance, double bonds locked in a ring smaller than an eight-membered ring often cannot isomerize due to steric constraints imposed by the ring. Cyclooctene is the smallest cycloalkene that is stable in the trans conformation at room temperature. 242 trans-Cycloheptene may exist at lower temperatures but at room temperature isomerizes quickly to ciscycloheptene through a low-energy pathway that involves formation of a diradical dimer. 243 The fact that smaller cycloalkenes cannot form isomerized trans products has been exploited in many instances for locking double bonds to prevent their isomerization. One example of this presents itself in cis-locked rhodopsins, which are artificial constructs where the C11C12 bond of the rPSB11 chromophore is locked in the cis conformation by a fivemembered ring (Figure 37). Cis-Locked bovine rhodopsin has been prepared in the lab and has an excited-state lifetime of ca. 85 ps, which is almost 3 orders of magnitude longer than the lifetime of wild-type (unlocked) rhodopsin. 244,245 Since cis-locked rPSB11 in bovine rhodopsin cannot isomerize along the C11C12 bond, computational studies have indicated that it decays by twisting of an adjacent C9C10 double bond instead. 246,247 The use of rings to lock double bonds in a certain configuration has also been useful for the design molecular switches, where all double bonds but one are incorporated in a ring to prevent their isomerization. 409 This makes it possible to construct light-driven molecular switches that isomerize selectively along one specific exocyclic double bond, such as in the case of the N-alkyl-indanylidene-pyrrolinium  (NAIP), 18,249−251 fluorenylidene−pyrroline, 252 and benzylidene−pyrroline rPSB11 mimics. 253 The same concept has also been used to synthesize HBDI mimics called HBDI-like switches. 254 While ring-locked double bonds may not be able to isomerize, they may still undergo internal conversion via a twisted CI similar to the case of an aborted isomerization. One example in proteins is found in the putative fluorophore of bacterial luciferase, which is a derivative of a flavin mononucleotide (FMN). In particular, the fluorophore is an N(5)-ethyl 4a-hydroxy flavin (Et-FlOH). While EtFlOH is strongly fluorescent in the bacterial luciferase protein, 256 the same fluorophore in solution has a very weak fluorescence. 248 A joint computational and experimental study found, using a methyl FlOH model (Me-FlOH, Figure 38A), that the first singlet excited state (which is the fluorescent state) is repulsive and does not have a fluorescent minimum ( Figure  38B). Instead, from the vertical Franck−Condon point, it evolves sequentially along a BLA mode (similar to that observed in other conjugated organic chromophores such as rPSB11) followed by an out-of-plane deformation of the pyrimidine ring. The deactivation BLA and torsional deformation route is barrierless on the excited state and leads to a twisted CI with the pyrimidine twisted almost 90°with respect to the rest of the flavin plane. The documented Me-FlOH deactivation mechanism and CI structure are in some ways similar to the internal conversion of pyrimidine nucleic acid bases (i.e., cytosine, thymine, and uracil), which also follows a BLA−torsional deformation deactivation mechanism. 255,257,258 However, there are some distinct differences between the case of Me-FlOH and pyrimidine bases, as shown in the comparison with uracil in Figure 38B. For example, the BLA pattern is almost opposite in the pyrimidine ring of Me-FlOH and in uracil. Moreover, the uracil CI has a twisted CC double bond (orange dihedral in Figure 38) that involves a high degree of pyramidilization of one of the carbon centers and can be accessed via a facile hydrogen out-of-plane twisting. The Me-FlOH CI, on the other hand, is only partially twisted along the corresponding CC bond but instead is mainly twisted along an adjacent NC bond (red dihedral in Figure 38). This twisting mode involves a large structural deformation due to movement of the entire pyrimidine ring with respect to the original flavin plane, which explains why such a mode is not accessed in a binding cavity and fluorescence is much higher in the protein.

CONCLUSIONS AND PERSPECTIVES
Photobiological systems are incredibly complicated, and a comprehensive understanding of their mechanism requires a full understanding of electronic and nuclear processes occurring on a range of time scales. Currently, a single model alone cannot possibly tell a full story. Therefore, photoreceptors have been the subject of numerous experimental and theoretical studies investigating different aspects of their mechanisms. In this review, we have focused on a relatively narrow topic: the QM/ MM modeling of double-bond isomerization in rhodopsins, photoactive yellow proteins, and phytochromes. Even with such a specific topic, it has been impossible to be exhaustive due to the volume of publications tackling these systems. Therefore, as mentioned in the Introduction, the present review is not fully comprehensive but revolves around the background and research interests of the authors. More precisely, it reflects the progress made by a number of computational chemistry schools that started to study photochemical reactions in the 1990's. One of these schools, which focused on the development and application of the CASSCF method, emerged from a collaboration between Fernando Bernardi and Michael A. Robb 259 and has contributed significantly to the computational results reviewed above. However, the theoretical background required for studying such reactions and, in particular, the double-bond photoisomerization of conjugated molecules, was laid out by numerous earlier studies. For instance, the electronic states driving the photoisomerization of a model rPSB had already been described by Josef Michl and Vlasta Bonacǐc-Kouteckýin the late 1980's. 19,260 Michl and co-workers also originally highlighted the pivotal role of CIs in organic photochemical reactions and tested these hypotheses using configuration interaction computations. 17,260 Even earlier than that, 2013 Nobel laureate Arieh Warshel provided the theoretical background for understanding the effect of the environment on the rPSB chromophore by developing, implementing, and applying, for the first time, QM/MM technologies to derive mechanistic information. Thus, effectively, Warshel initiated the transition from computational photochemistry to computational photobiology using semiempirical quantum chemical methods before practical ab initio multiconfigurational quantum chemical tools became available.
In the 1990's, several groups carried out methodological and applicative research in computational photochemistry (e.g., CI optimization methods, reaction path calculations, and trajectory propagation) that made the studies reviewed here possible. Such work is only marginally covered here but has been the subject of a recent review. 17 For instance, Schlegel and co-workers developed an effective method for locating CIs. 261 Yarkony and coworkers 262 developed rigorous methods and theories for locating and characterizing CIs. 60,263 Martínez and co-workers have Chem. Rev. XXXX, XXX, XXX−XXX reported on practical CI optimization methods and nonadiabatic dynamics approaches and used them to investigate several biological chromophores. 114 Morokuma and co-workers have reported an effective branching plane calculation method. 264 Authors such as Schulten 29 and Birge 99 were among the first to perform semiclassical trajectory computations for Rh, bathoRh, 85 and bR. 108,265,266 Needless to say, such progress in photobiology was dependent on progress on the theoretical photochemistry front. For instance, Domcke, 258,267 Koppel, 268 Cederbaum, 269 and Martínez 270 and co-workers, among others, played a vital role in rigorously demonstrating how ultrafast decay occurs at CIs. The demonstration, by quantum dynamics calculations, that few vibrational degrees of freedom are sufficient for irreversible radiationless decay at CIs led to a fundamental change of paradigms in radiationless decay theory. This understanding has been central to photochemistry in general and complemented the above-mentioned computational demonstrations that lowlying CIs are ubiquitous and common features of the excited states of sizable organic and biological molecules.
In conclusion, several authors created the scientific and technical basis for the computational investigation of ultrafast photochemical reactions in light-responsive proteins and, therefore, for the understanding of the biological function of double-bond photoisomerization. In the following subsections, we draw some general conclusions (section 8.1) and highlight directions (section 8.2) that, in the authors' opinion, would be important for evolving from a scientific understanding toward the design of new molecular materials. In fact, until recently, theoretical chemists working on biological photoreceptor proteins have mainly focused on method development with the aim of exploring and uncovering the molecular mechanisms underlying protein function. Today we are witnessing a change of paradigm as some of the same scientists who studied the mechanism of photobiological systems are becoming more and more like engineers who can use computations and lab models to design new functions. However, this transformation requires integration of methods for predicting and analyzing conformational changes, hydrogen bond networking, protein-chromophore interactions, and thermal/photochemical reactivity before revealing new alleys toward the engineering of unprecedented systems.

General Conclusions
The experimental investigations of light-responsive proteins operating through double-bond isomerization have examined different "stages" of the primary photoproduct formation. When referring to Scheme 1, these are light absorption and initial relaxation, reaction coordinate, excited-state intermediate and/ or CI structures, and, finally, reaction time scales and quantum yields. In this review, we have tried to discuss selected groups of computational investigations addressing one or more of these stages in different proteins. Here, we will try to summarize the general mechanistic information that has resulted from these investigations.
8.1.1. Initial Relaxation. The initial event in the "activation" of the examined biological chromophores corresponds to an allowed π−π* excitation which may be described as a chargetransfer state. In cationic chromophores such as rPSB or PΦB, this implies a partial translocation of the positive charge away from the CN Schiff base (in rPSB) or CN pyrrole bond (in PΦB), while in anionic chromophores such as in pCA, this implies a translocation of the negative charge away from the phenolic oxygen. In the framework of the Franck−Condon approximation, such instantaneous change in the electronic structure creates a nonequilibrium nuclear configuration that then relaxes out of the FC point. One of the first results of the computational work presented above has been the description of such a motion that is dominated by the expansion of the double bond and contraction of the single bond (the BLA coordinate). This implies the "unlocking" of the double bonds of the system as the primary effect of photon energy absorption. Notice (see section 3.1) that in neutral and symmetric conjugated systems such as polyenes, these states would not promptly generate a charge translocation and therefore an inversion between double bonds and single bonds but, rather, a charge separation producing structure with almost equivalent single-and doublebond lengths.
The computational investigations reviewed here have also indicated that the BLA relaxation can be effectively coupled to other molecular modes by the chiral protein environment. One example is the protein-induced stabilization of a chromophore conformation with a pretwisted double bond. This would effectively couple the BLA mode with specific torsional molecular modes resulting in a multimodal initial excited-state relaxation. For instance, we have seen that in certain cases, such as that of animal rhodopsins, the initial relaxation effectively couples the BLA and the HOOP hydrogen wagging motion. A few femtoseconds later this bimodal relaxation couples with the twisting deformation of the carbon framework and, effectively, initiates the evolution along the isomerization coordinate. It is thus apparent that the coupling described above is central to determining the stereoselectivity and regioselectivity of the isomerization. This can be sterically or electrostatically imposed by pretwisting a specific double bond of the chromophore and/ or inducing an asymmetric force field (i.e., potential energy surface shape) that results in stereoselectivity/regioselectivity of the isomerization. Such steric and electrostatic effects result in the selection of the C11C12 or C13C14 bonds instead of the C9C10 bond in rPSB, the β instead of α isomerization in pCA, or the C15C16 bond instead of the other double bonds in PΦB. On the other hand, computational studies have shown that the pCA chromophore in gas phase or solution is more likely to isomerize along α, 138,139,147 while PΦB in vacuo is more likely to isomerize along C9−C10 or C10C11. 197−199 Therefore, the protein environment has a large control on the stereoselectivity of photoisomerization for different chromophores. Understanding this control has already been the subject of much interest for the computational community because a better understanding of such control is central to the engineering of artificial proteins.
8.1.2. Reaction Coordinate. Comparison of excited-state reaction paths and trajectories of different isomerizing chromophores (e.g., rPSB, pCA) reveals that the bicycle-pedal or aborted bycicle-pedal mechanism is commonly encountered in photobiology. This is because such a coordinate is space-saving but also favorable if one considers the electronic character of the excited state of these chromophores. In fact, the initial relaxation along BLA naturally unlocks two nearby double bonds. In contrast, a Hula-twist coordinate would require the twisting deformation of a double and a single bond simultaneously, thus increasing the energy cost necessary to achieve such a spacesaving motion. Of course, when a bicycle-pedal motion cannot take place, such as in PΦB or in the HBDI chromophore of fluorescent proteins featuring a single double bond located outside a cyclic framework, the Hula-twist reaction coordinate may become operative if leading to a favorable space saving Chemical Reviews Review DOI: 10.1021/acs.chemrev.7b00177 Chem. Rev. XXXX, XXX, XXX−XXX motion. The Hula-twist reaction coordinate may become competitive with the bicycle-pedal coordinate when there is a penalty in undergoing bicycle-pedal type motion, as in pCA where the bicycle pedal type motion requires breaking the hydrogen bond between the carbonyl oxygen of the chromophore and the backbone amide in a nearby amino acid. Notice that in polyenic hydrocarbon systems, the initial relaxation along the locally excited (charge separated) π−π* state has a high probability of populating a lower energy diradical electronic state with different potential energy surface features.
As we have discussed in section 3.1, such electronic structure leads to a kink-type conical intersection that forces evolution along a Hula-twist-like excited-state reaction coordinate and, ultimately, conical intersection. Even if the bicycle-pedal reaction coordinate seems to be the most favorable, it is often aborted upon decay, leading to the effective isomerization of a single double bond. The selection of the aborted isomerization has not been thoroughly investigated but in rPSB is connected with the amplitude of the twisting deformation. For instance, in Rh a ca. 90°twisting about the reactive C11C12 bond is accompanied by just a 60°twisting about the adjacent C9C10 bond, which, therefore, has the largest probability of reverting back to its original configuration. In the bicycle-pedal reaction coordinate of PYP, the chromophore also has a ca. 90°twisting about the alkene moiety while the adjacent CS bond is only ca. 50°twisted.
8.1.3. Excited-State Dynamics and Decay at the Conical Intersection. Trajectory calculations have largely confirmed the character of the excited-state electronic and geometrical progression indicated by reaction coordinate studies. On the other hand, new information has been provided on the factors that control the excited-state lifetime (i.e., the photoisomerization speed) and isomerization quantum yields. Animal rhodopsins have emerged as the fastest photoisomerizing proteins with estimated excited-state lifetimes as low as 60 fs 24 and experimentally detected lifetime as low as 70 fs in bovine rhodopsin 32 (in fact, in a recent paper Miller and co-workers emphasized the possibility that the initial photoproduct could form in an even shorter time-scale). 14 This speed has been attributed to factors related to the slope of the barrierless potential energy surface connecting the FC point to the CI. However, it has also been documented that a prerequisite for achieving such a speed is the energetic separation between the potential energy surfaces of the first excited state (i.e., the spectroscopic π−π* charge-transfer state) and the higher electronic states such as the nearby diradical excited state. This separation guarantees a continuous acceleration toward the CI without entering a rugged and flat potential energy surface where the nature of the electronic character of the excited state may vary. This last factor has been proposed to explain the different dynamics observed for animal rhodopsins and microbial rhodopsins. 23 The reason for these differences is related, of course, to the less distorted all-trans rPSB isomer in microbial rhodopsins, which features the isomerization of the C13C14 rather than C11C12 bond.
A second reason for observing a slower double-bond photoisomerization dynamics in proteins is the formation of transient excited-state intermediates whose reactivity is restrained by an excited-state barrier or by the existence of a sloped rather than peaked CI. This has been computationally found in the pCA chromophore of PYP and the PΦB chromophore of phytochromes. In PΦB, the intermediate appears to correspond to a TICT state featuring a complete charge translocation and the same 90°twisted geometrical structure of the CI but located energetically above it.
Finally, as documented in section 3.4, the analysis of the result of nonadiabatic trajectory calculations has started to clarify the relationship between the reaction speed and the quantum yield. This appears to be unrelated to a simple proportionality relationship as a ballistic picture of the isomerization motion would suggest. In contrast, the currently emerging picture is focusing on the achievement of a phase relationship between the different vibrational modes populated immediately after photoexcitation and therefore on an extremely short time scale. Such a control would require an ultrafast excited-state dynamics of the excited-state population whose function would be to guarantee the conservation of the molecular level (not necessarily population level) coherence between the modes that contribute to the reaction coordinate (e.g., the HOOP and the skeletal reactive torsional motion). In accordance with this mechanistic model, a productive trajectory will have to display a certain phase relationship at the moment of decay. A faster excited-state motion (or a slower one) may result in a phase relationship, allowing exclusively the return to the initial reactant double-bond configuration and, statistically, in a quantum yield decrease. Presently, such a mechanism has been considered exclusively in animal rhodopsin. Microbial rhodopsins may not utilize it (e.g., when using the rPSBAT chromophore) but instead could utilize other mechanisms for achieving high quantum yields that are applicable to situations where the excited-state motion is slower and less coherent.
Finally, in an ultrafast regime where vibrational coherence is maintained during the entire excited state evolution, photoproduct formation seems to occur according to an oscillatory mechanism, which appears to be connected/determined by the HOOP-like motion occurring along the reaction coordinate. This has been reported for an abridged model of ASR 13C , 73 a microbial rhodopsin isomer for which vibrational coherence has been observed (see Figure 39). 70 This type of behavior and its significance is presently investigated using more realistic models (e.g., a 5 double-bond chromophore model) as highlighted in the graphical TOC of the present review. Most of the computations reviewed in sections 4−7 have been possible thanks to decades of developments leading up to modern electronic structure methods, such as those discussed in section I. Ideally, to study the mechanism of photobiological systems, one requires methods that can (a) quantitatively describe the relative energies of ground and excited-state potential energy surfaces, (b) describe the near-crossing and crossing regions of potential energy surfaces, and (c) be computationally affordable. When one wants to describe transitions through crossings and avoided crossings, accurate nonadiabatic couplings are also needed. CASSCF is one of the methods that most prominently has made possible, and still makes possible, 271 the investigation of biological chromophores. The CASSCF implementations by Michael A. Robb 259 and by the member of the Swedish school led by Bjorn Roos had, as far as the present authors are concerned, a major impact. As a consequence, computer packages such as GAUSSIAN 272 and MOLCAS, 273 allowing the integration of state-average CASSCF gradients and nonadiabatic couplings into QM/MM technologies, have been applied widely. However, CASSCF does not Chemical Reviews Review DOI: 10.1021/acs.chemrev.7b00177 Chem. Rev. XXXX, XXX, XXX−XXX correctly describe reaction barriers and relative energies between electronic states [i.e., (a) above]. For this reason, the CASPT2 method, developed by Roos and co-workers, had, and still has, 274 a paramount importance for correcting CASSCF energies. However, CASPT2 gradients are still too expensive to be systematically employed in the type of studies describe above (c). Furthermore, CASPT2 suffers from some artifacts when it comes to describing near-crossing regions (b), although this has been remedied by more recent extensions to the method. Such drawbacks notwithstanding, this is the reason for the popularity of the CASPT2//CASSCF protocol in photobiology. The issues above indicate that there is a need for electronic structure methods that can yield accurate but affordable QM/ MM models. The further development and efficient implementation of analytical gradients for CASPT2 and multistate multiconfigurational methods such as MRCI, 275−277 MS-CASPT2, 278 XMCQDPT2, 118 QD-NEVPT, 279,280 and XMS-CASPT2 281,282 will certainly impact trajectory calculations and mechanistic studies. In fact, recently, Shiozaki and co-workers have reported both analytical gradients 283 and nonadiabatic couplings 284 computed at the XMS-CASPT2 level. Meanwhile, many research groups continue to push the boundaries of electronic structure theory as methods continue to be developed, improved, and made more efficient. Some of these methods are already used to study photobiological systems, while others could very well be used in the near future. TD-DFT has been widely used, although primarily for geometry optimizations and spectroscopic studies due to its limitations when it comes to describing near-intersection regions. However, several developments addressing these near-intersection issues, such as DFT-configuration interaction methods (e.g., CIC-TDA-TD-DFT, 285 CDFT-CI, 286,287 SI-SA-REKS, 182,288 and spin-flip TD-DFT) 289,290 may result in more mechanistic studies that employ DFT. DFT methods are also being combined with multiconfigurational methods such as CASSCF to account for the missing static correlation energy component. 291 Coupled-cluster approaches such as linear-response CC and EOM-CC are also used due to their high accuracy and systematic improvability, and their application to photobiological systems is likely to increase as more and more efficient implementations become available. 292,293 This is also true for ADC methods. 294−296 Quantum Monte Carlo methods 297,298 are also emerging as players in the photobiology community due to their favorable computational scaling and ability to account for electron correlation. These are just a few examples, out of many, of established or emerging electronic structure methods for studying biological chromophores.
While CASSCF and CASPT2 are used routinely for systems such as rPSB and the pCA chromophore, certain biological chromophores are still too large for such methodologies. This is in spite of recent progress in parallelization and fast algorithms such as the Cholesky decomposition. 299,300 Therefore, we often find that larger chromophores, such as PΦB or carotenoids, are not as well-studied as rPSB or pCA. For example, a CASSCF calculation on PΦB that includes the full π-electron system in the active space would require an active space of 28 electrons in 25 orbitals, which is beyond the current computational limit of CASSCF. Although a reduced active space could be used, 201 this could lead to complications in the course of dynamics because the relative importance of selected configurations can vary in time as the geometry changes. However, development of new approaches built on the CASSCF concept, such as RASSCF, 301 GASSCF, 302 and DMRG, 303,304 could help overcome such problems.
In order to increase the quality of a QM/MM model, not only is the electronic structure method describing the QM subsystem important but also the description of the interaction between QM and MM subsystems. In fact, the search for more accurate and computationally effective QM/MM models is another important technical challenge that computational photochemistry and photobiology faces. This calls for the development of more realistic force fields to describe the protein and its interaction with the QM subsystem. The importance of polarizable embedding (section I.2.1) has been pointed out in several pilot studies applied to the calculation of excitation energies in photoreceptor proteins. 305−307 However, while several polarizable embedding implementations exist and have been used to study spectroscopic properties, there are no mechanistic or dynamics studies using polarizable embedding yet. Further implementations are needed to realize a combination of polarizable embedding with nonadiabatic dynamics. The same can be said for DFT embedding schemes, which are also powerful and accurate but need to be further developed before they can be used in mechanistic/dynamic studies of photorecetors. Finally, another promising approach for combined quantum/classical description is Effective Fragment Potential (EFP). 308 In this approach, a potential for solvent molecules is generated by calculating several parameters from ab initio calculations. The extension to proteins (macromolecule EFP, or mEFP) is very promising for future applications to photoreceptors. 309,310 One brute-force attempt for overcoming the problem of the interaction between the chromophore and its environment is that of extending as much as possible the size of the QM subsystem. This can be achieved in calculations including clusters of amino acids surrounding the chromophore in the QM subsystem, all the way to treating the entire protein quantum chemically. The latter requires exceptional computing power not currently provided through conventional machines. However, hardware developments have made such calculations accessible. Already a few years ago, it has been shown that video gaming machines containing two consumer graphical cards can outpace a state-of-the-art quad-core processor workstation by a factor of more than 180 for Hartree−Fock energy and gradient calculations. 311 Such performance enhancements make it possible to run large-scale Hartree−Fock and DFT calculations, which typically require hundreds of traditional processor cores, on a single graphical processing unit (GPU) workstation. For instance, benchmark molecular dynamics simulations are performed on molecular systems such as an aspartic acid molecule solvated by 147 waters (457 atoms, 2014 basis functions) at the HF level at a rate of 0.7 ps of computed synamics per day on a single desktop GPU. 311,312 Further development and application of this technology, especially for post-Hartree−Fock methods, is expected to have a big impact on the field. 313,314 8.2.2. We Need to Look at Electron, Not Only Nuclear, Dynamics. In Figure 4 and Section 3.4, we have discussed the importance of phase relationships between different molecular modes. Such a relationship is established within a few femtoseconds after photoexcitation but ultimately has an impact on how excited-state dynamics unfold and therefore on the point when molecules hop from the excited state to the ground-state potential energy surface. In order to completely understand the onset of such phase relationships, we must answer questions related to the validity of the Franck−Condon principle and exact separation of nuclear and electronic motions. In fact, the chargetransfer electronic structure achieved upon photoexcitation involves a reorganization of the electron cloud and, possibly, an intimate coupling between such reorganizations and the onset of nuclear motions (e.g., the early BLA and the HOOP motions discussed for visual pigments). This type of coupled electronic and nuclear relaxation have already been computationally studied in ionization processes. 315 The same type of electronic-nuclear coupling may affect the "hop" event itself, influencing the nuclear motion and, possibly, the excited-state lifetime and quantum efficiency upon passage through the CI region. 316 These phenomena are connected with the rapidly advancing field of attosecond spectroscopy that could allow the observation of such coupling phenomena in biological chromophores on the time scale of the electronic response (the attosecond time scale). 317,318 In molecules, attosecond pump−probe spectroscopy enables investigations of the prompt charge redistribution and localization that accompany photoexcitation processes, where a molecule is lifted from the ground Born−Oppenheimer potential energy surface to one or more excited surfaces and where subsequent photochemistry evolves on femtosecond or longer time scales.
8.2.3. We Need to Investigate the Coupling between Isomerization and Other Processes. The engineering of genetically encodable molecular "machines" powered by light requires the understanding and simulation of coupling between double-bond photoisomerization and other chemical processes such as proton transfer. For instance, it is established that, in fluorescent proteins like GFP, 319 the excited-state HBDI chromophore deprotonation is preceding the formation of the fluorescent intermediate while in the Dronpa (section 7.2) proton transfer is coupled to double-bond photoisomerization. 320 Similarly, in rhodopsins, the ground-state deprotonation of the chromophore follows the isomerization process at some stage along the photocycle. 321 While these functionally important processes are usually seen as sequential, little is known about the actual mechanism and time scale of the coupling. For instance, in rhodopsin, the possibility that the isomerization and even the excited-state charge translocation occurring along the reaction coordinate 322 could destabilize the hydrogen bond network and trigger successive events has been proposed and investigated both experimentally and computationally (for a recent work see ref 323) but still with no conclusive picture.
On the other hand, coupling between proton transfer and electron transfer (PCET) has been shown to be implicated in a number of photobiological systems such as, most recently, in blue light using flavin adenine dinucleotide (BLUF) proteins. 324 In this case, free-energy simulations elucidated the role played by the active site conformations before and following photoexcitation.
8.2.4. Quantum Chemistry and Molecular Evolution. In sections 4.2.2 and 4.3.3, we have shown that QM/MM models can be employed for comparative investigations of the thermal and photochemical reactivity of rhodopsins. These studies belong to research lines where quantum chemistry is employed for solving biophysical problems. Further extension of this approach to other light-responsive proteins, such as the ones treated in sections 4−6, would certainly bring new progress. On the other hand, the authors of the present review believe that, in the near future, quantum chemical studies should focus on the possibility to model large sets of light-responsive proteins. This would be of basic importance for the development of computational photobiology, as several biological problems require the comparison of large sets of proteins. For instance, QM/MM models can be constructed for proteins representing the nodes of a phylogenetic tree and, therefore, protein sequences presently extinct (i.e., ancestral proteins). In this way, the unavoidable overlap between, for instance, quantum chemistry and molecular evolution studies should lead to new research lines. This includes the search into the past history of a protein but also information on how to engineer it (a chemical target connected with the preparation of new materials) or, simply, learning about how a certain extant protein has been achieved as a function of time (a biological target connected with evolution). In other words, QM/MM models would allow the user to virtually "resurrect" a protein in the computer, thus providing insight into the evolutionary history of molecules, a window on the past revealing how adaptive specializations may have occurred in nature (see also the Lake Baikal fish rhodopsin example 96 in section 4.2.2). Furthermore, this same information could provide inspiration for the engineering of highly efficient light-converting biomimetic systems such as the one seen in section 7.1. In the more extreme view anticipated in section 7.3, one can learn from such studies how to engineer efficient components for light-driven molecular devices including coupled chromophore/optical cavities or coupled chromophore/plasmons. 325,326 Rhodopsins, PYP's, and phytochromes appear to be excellent targets for understanding how evolution has gradually taken control of such an elementary, but versatile, process as lightinduced double-bond isomerization. For instance, the visual rhodopsins of invertebrates and vertebrates are known to possess Chemical Reviews Review DOI: 10.1021/acs.chemrev.7b00177 Chem. Rev. XXXX, XXX, XXX−XXX extraordinary adaptations for maximizing photosensitivity under low light conditions. To understand how such sensitivity evolved, the laboratory resurrection and functional characterization of a 200 million year old archosaur rod visual pigment has been reported. 327 However, due to the effort and expense associated with laboratory manipulation, this experimental approach is limited to small numbers of inferred sequences, making the detailed description of important evolutionary transitions and the testing of evolutionary hypotheses difficult. Even if a sufficiently large number of ancestors could be resurrected in the lab, the related structural information could be difficult to achieve. In fact, this would require sufficient expression and the availability of crystal structures. These difficulties can be significantly decreased by in silico screening of the many, potentially interesting, ancestor structures. Such screening would be based on the fast production of consistent QM/MM models of rhodopsin proteins as well as of a set of mutants. The study of the corresponding model would allow biologists to limit the lab efforts to key structures carrying the maximum amount of information.
While approaches close to the one highlighted above have been reported by Morokuma, Yokoyama, and co-workers, 328−331 their applications focused on spectroscopy (color perception tuning) rather than photochemical reactivity. Furthermore, the report focused only on one or few QM/MM models. In order to overcome these limitations, hundreds, if not thousands, of QM/ MM rhodopsin models can be generated automatically. Also, these models need to be based on ab initio multiconfigurational quantum chemistry, as they have to be able to also reproduce changes in the mechanism and dynamics of photochemical and thermal activation. A first attempt in such direction has been recently reported. 332 Such technology, discussed in some more detail in section I.2, is called automatic rhodopsin model (ARM) and may be used to compute, for instance, the optical and photochemical properties of inferred ancestral sequences, then compared them to modern visual pigments. The results of benchmark studies carried out using rhodopsins from different organisms (i.e., archaea, eubacteria, invertebrate, and vertebrates) as well as eubacteria and vertebrate rhodopsin mutants, are displayed in Figure 40 (panels A and B, respectively).
In the future, the community may be able to implement research projects that integrate currently distant disciplines such as ab initio quantum chemistry, computational molecular evolution, and ancestral protein resurrection to achieve an unprecedented understanding of the functional divergence of light-responsive proteins. Starting these studies from proteins undergoing ultrafast double-bond isomerization will allow the exploitation of existing knowledge on these systems. Accordingly, more effective computational protocols for the systematic construction of computer models of ancestral pigments need to be developed and benchmarked.

I. TECHNICAL SECTION: METHODOLOGY OVERVIEW
The acronyms encountered in several figure legends, such as 3root-SA-CASPT2(IPEA = 0)//2-root-SA-CASSCF/6-31G* and 2-root-SA-CASSCF/6-31G*/AMBER, refer to methods employed to calculate ground and excited-state potential energies and gradients. These acronyms are actually an extension of Pople's acronyms 333 used to indicate different "model chemistries". 334 The double slash "//" sign separates the level of theory used for computing the energy of the system (at the left) from the level of theory used for optimizing the geometry or computing the trajectory (at the right). Thus, 3-root CASPT2// 2-root CASSCF means that energies were computed with 3-root CASPT2 as single-point calculations on top of geometries that were obtained with 2-root CASSCF. The 3-root CASPT2 acronym implies that the corresponding zeroth-order wave function is the 3-root CASSCF wave function. The IPEA = 0 label is a required CASPT2 parameter. 273 It is important to have a certain degree of understanding of these methods as the quality of the chemical information extracted from the calculations depends critically on their properties and limitations. More specifically, such acronyms refer to electronic structure methods that cannot be used in a blackbox fashion and necessarily involve a number of approximations that need to be considered carefully to achieve a balanced description of spectroscopy and ground and excited-state reactivity.
A balanced description of the ground and excited state electronic structure appears to be necessary for the chemical process that involves the breaking and making of π-bonds. In fact, as highlighted in several parts of the present review, double-bond isomerization may occur when a π-bond is broken homolytically, a process that involves the unpairing of electrons, or heterolytically, a process that involves a transfer of the bonding electron pair to one side of the broken bond. Therefore, the modeling of isomerization requires a suitable level of electronic structure theory that is capable of describing changes in electron spinpairing and electron densities during isomerization. At the same time, the study of isomerization in photobiological systems requires an understanding of the effect of the protein , invertebrate, and nonvisual vertebrate rhodopsins (blue), eubacterial sensory rhodopsin isomers (green), archaea rhodopsins (yellow), and a chimeric rhodopsin ChR C1C2 , originating from a microbial eukaryote. A proton-pumping eubacterial rhodopsin PR is also reported. The colors connecting the data points indicate related rhodopsins (i.e., related by sequence homology as for sqRh and hMeOp, or because they are isomers as in both the bR and ASR cases, or because they are photocycle intermediates as for Rh 181  environment on such a chemical reaction. These include the steric and electrostatic effects that influence the electronic structure of the reactive chromophore. However, due to the high scaling of the cost of electronic structure methods with the size of the system, the treatment of an entire protein (i.e., a system containing thousands of atoms) is usually impractical if not impossible. Therefore, multiscale methods, where the system is divided into two or more subsystems treated at different levels of theory, have become common tools for studying biological photoreceptors. When the protein and the chromophore are connected via a bond, a suitable approach to treat the frontier of these two subsystems is also required.
The present section comprises four parts reviewing the computational methods used in exploring double-bond isomerization in photoreceptor proteins. The first part (section I.1) will focus on electronic structure methods commonly used for treating the breaking of π-bonds in the ground and excited states of biological chromophores. The second part (section I.2) will describe the details of QM/MM multiscale methods. In the third part (section I.3), the computational tools and the details regarding the calculation of nonadiabatic trajectories will be presented.

I.1. Electronic Structure Methods
The development and testing of electronic structure methods is a very active field of current theoretical chemistry research, with new methods being continuously developed by many groups around the world. Therefore, even if we limit ourselves to methods suitable for studying excited-state isomerization, a comprehensive review would require the presentation of a large arsenal of methods. For this reason, this part will only focus on a specific series of benchmark studies testing the ability of some of the more popular classes of electronic structure methods employed to study the energy changes occurring along the double-bond isomerization coordinate of a specific chromophore. We will not go into any detail regarding the formulation of these methods and their performance, especially when it comes to studying other properties. We refer the readers to refs 335 and 336 for a general background to excited-state electronic structure methods.
The benchmark studies that we will discuss employ a reduced model of rPSB11 that comprises three double-bonds (PSB3, see Figure 41A). Reduced models of rPSB11 have been used in many instances for benchmarking purposes or even as a model to understand the photochemistry of rPSB11 (for instance, see work by Olivucci, Garavelli, Robb and co-workers, 16,20,191,337−341 Elstner and co-workers, 342,343 Buss and co-workers, 33 Birge and co-workers, 344 Thiel and co-workers, 345 Send and Sundholm and co-workers, 346,347 Filippi and co-workers, 348,349 Andrunioẃ and co-wokers, 350,351 and Wong and co-workers 352 ).
PSB3 has many of the same features of the S 0 and S 1 potential energy surfaces as rPSB11. In both PSB3 and rPSB11, excitation from the S 0 to the S 1 electronic state corresponds to the 1A g to 1B u electronic transition described in Section 3.1 and results in a 30−50% positive charge transfer away from the Schiff base moiety. 20,21,353 In the S 1 state, PSB3 and rPSB11 feature a very shallow barrier or barrierless path connecting the FC point to a peaked CI characterized by a ca. 90°twisted isomerizing bond, as shown in Figure 9A. 21,62,192,348 In both systems, the electronic character of the S 0 and S 1 potential energy surfaces varies between diradical and charge transfer along a loop encircling the CI and located along the BP. 57,90 PSB3 also reproduces several important features of the rPSB11 S 0 potential energy surface; specifically, PSB3 features two transition states (TS CT and TS DIR ), the former having charge-transfer character and the latter having diradical character, as shown in Figure 9B. 56,90 Therefore, PSB3 serves as a model system for rPSB11 and has been adopted to investigate the performance of different electronic structure methods. Its reduced size permits the use of methods not affordable for rPSB11.
Recently, we have mapped specific S 0 and S 1 paths of PSB3 and used them to determine the effect of the electron correlation on different regions of the corresponding potential energy surface. As explained below, the electron correlation energy is the electronic energy missed by the Hartree−Fock (HF) method where the interactions between electrons is accounted for only in a mean field way. This work was a collaborative effort involving several research groups and published over several installments. 35  (C) A schematic 2D representation of the S 0 and S 1 potential energy surfaces of PSB3 showing six of the seven pathways used to benchmark electronic structure methods (see text). The coordinate dominating the MEP CIS and MEP TRANS corresponds to a changing mixture of BLA, torsion, and HOOP and therefore cannot be represented completely on the single 2D potential energy plot shown here. The MEPs paths are represented better in a 3D coordinate space plot as discussed in section II. A specific torsion and HOOP mixture is also running along the CI intersection space (the IS curve already mentioned in section 3.2 and reported in Figure 4). The paths are labeled in yellow full (or dashed for hidden segments) lines, while key structures on the potential energy surfaces are labeled with blue (for minima and Franck−Condon points), orange (for transition states), and red (for conical intersections) full circles. The seventh benchmark path, the CI loop, is shown in Figure 42 because it lies on the BP orthogonal to the IS space. seven different paths (see Figure 41). These paths have mechanistic relevance since they involve key structures on the S 0 and S 1 potential energy surfaces of PSB3. The first path, named MEP CT , is a minimum energy path on the ground state that connects the cis-PSB3 equilibrium geometry, TS CT , and trans-PSB3. The second path, MEP DIR , connects cis-PSB3, TS DIR , and trans-PSB3. The third path is the BLA path, defined in Figure  41B, which connects the TS CT and TS DIR transition states but also intercepts a CI point (hereafter referred to as CI BLA to distinguish it from other CI points located along the IS hyperline). Since the BLA mode removes the electronic degeneracy, it is one of the two modes belonging to the BP of the CI. The fourth path, MEP CIS , is an S 1 path starting from the cis-PSB3 FC point and following a minimum energy path leading to a CI in the IS. This CI is geometrically and energetically distinct from CI BLA and is labeled CI CIS . Similarly, MEP TRANS connects the trans-PSB3 FC point to another point on the IS hyperline, CI TRANS . The sixth path, the CI "seam segment", belongs entirely to the IS and connects CI CIS , CI BLA , and CI TRANS . MEP CIS and MEP TRANS span a mixed coordinate space dominated by BLA near the FC region but then dominated by a mixture of torsion and HOOP. CI CIS , CI BLA , and CI TRANS differ mainly along the torsion and HOOP coordinates. The final path, discussed in more detail below, is called the CI loop and is a circular path that lies in the BP and encircles CI BLA .
The purpose of these benchmark studies is to go beyond understanding the performance of methods in reproducing one energy value (e.g., a vertical excitation). Instead, this allows us to look at the differential electron dynamical correlation effects on the potential energy surfaces involved in thermal and photochemical double-bond isomerization. Differential electron correlation energy arises due to changing electron density distribution along the molecular framework (e.g., when going from a delocalized to more localized electronic distribution), requiring different degrees of electron correlation to be accounted for. The seven pathways described above go through regions of the S 0 and S 1 potential energy surfaces with different electronic configurations, and therefore, their correct description requires a treatment of the changing electron correlation on an equal footing. Details regarding the construction of these paths and the results of the benchmarks using different methods are in refs 35, 62, 90, 182, 288, 297, and 354−357. We will focus on the results of two paths: the BLA path and the CI loop of Figure 42. The BLA path provides a stringent test for electronic structure theories because the states involved are dominated by open shell diradical (Φ DIR ) and charge-transfer (Φ CT ) characters, respectively (see Figure 43). Φ CT has a more localized π-electron density with respect to Φ DIR , and therefore, the computation of the corresponding potential energies requires a method accounting for the variation in dynamic electron correlation energy. Along the BLA path, an unbalanced description of these electronic characters will generate an artificial shift in the S 0 and S 1 energies resulting in a change in the position of their crossing, giving a crossing geometry that depends on the accuracy of the method used.
The CI loop is meant to probe whether an electronic structure method is capable of correctly describing the BP associated with a CI. Such loops encircling a CI point are often used to characterize the role of the CI in the photochemical process 15 or to probe the topography of the CI. 57,58,358−360 The BP has been introduced and discussed briefly in section 3.2. At a crossing between two electronic states of the same spin and symmetry, distorting the molecule along one of two specific degrees of freedom is the only way to lift the electronic degeneracy: GDV is the gradient difference vector, DCV is the derivative coupling vector, R represents the nuclear coordinates, E I and E J are the potential energies of the intersecting states, and ψ I and ψ J are the corresponding wave functions. These vectors are related to the g and h vectors; 60 GDV and g are identical, while DCV = h/(E I − E J ) and is therefore related to h by an energy factor. GDV and DCV define the BP (or the g-h) plane. 58 Notice that, in practice, the BP vectors are not uniquely defined but the branching plane is. For instance, near a CI, the vectors may change as a function of geometry even when the change is infinitesimally small. This is equivalent to moving along an infinitesimal loop center on the true CI point.
The To satisfy the condition of degeneracy (E I = E J ), the radicand must be equal to zero, meaning that the off-diagonal coupling terms (H IJ and H JI ) must vanish and that H II and H JJ must be equal. Meeting those two conditions requires two independent nuclear coordinates. 361,362 These are the GDV and DCV. This gives the state crossing a conical shape, as shown in Figures 8 and  42C. The correct description of the BP can be compromised if an electronic structure method does not correctly describe H IJ and H JI . Some methods cannot account for coupling at all (i.e., the coupling terms are zero at any possible geometry). In such cases, the BP gets reduced to just one dimension, and the crossing becomes linear rather than conical. 181 A similar problem occurs in electronic structure methods that apply energy-scaling corrections to H II and H JJ but that do not account for the effect of the scaling on the coupling terms. In such cases, the intersection is displaced to another region of the potential energy surface but the S 0 and S 1 wave functions and the couplings are not correctly described. In the cases described above, the S 0 and S 1 energy profiles along the CI loop yield a qualitatively different shape when compared to those of Figure 42D. In fact, instead of remaining nondegenerate along the entire loop, the two states will cross twice along the loop due to the fact that one of the BP vectors no longer splits the degeneracy. The energy difference between the two states will therefore drop to zero at two points along the loop (see examples in Figure 45). In conclusion, the CI loop distinguishes between methods that yield an unphysical one-dimensional linear crossing from methods yielding a BP with the correct dimensionality.
It is worth mentioning that in non-Hermitian theories, H IJ and H JI are not necessarily equal, introducing an additional condition (and change along an additional geometrical coordinate) for reaching the degeneracy. In such cases, the BP may become a three-dimensional "branching space". 363,364 The CI loop test cannot distinguish between methods that have a BP with a dimensionality larger than 2, such as in non-Hermitian methods. This is also true for cases where crossings between the potential energy surfaces of more than two electronic states occur.
Finally, notice that in the case of crossings between two electronic states of different spin multiplicities (e.g., singlet− triplet crossings), H IJ and H JI are identically zero in the framework of nonrelativistic quantum chemistry and so the corresponding crossings are also always linear. Spin−orbit coupling, however, can couple states of different multiplicity and therefore lifts the degeneracy formally yielding a conical crossing. 365,366 The PSB3 potential energy profiles along the BLA scan and CI loop computed with several electronic structure methods are shown in Figure 44 and Figure 45, respectively. These benchmarks have also been extended to approximate linearresponse coupled-cluster (CC2), 356 and a range of DFT methods, including spin-restricted ensemble-referenced Kohn− Sham (REKS). 90,182,288,357 Furthermore, the quantum Monte Carlo (QMC) method has been benchmarked using the same system. 297 Although QMC was originally a family of ground state methods, its extension to excited states has started somewhat recently. In fact, QMC methods for electronic excited states have started to show promising results in photochemical applications due to its favorable computational scaling and ability to account for electron correlation. 367−369 A large family of methods that is missing from these benchmark studies are semiempirical methods. This includes popular methods used for photochemistry including ZINDO/ S, 370,371 MNDO/MNDOC, 372,373 excited-state variants of PM3, 374 AM1, 375 and orthogonalization-corrected OMx methods (OM1, OM2, and OM3). 376,377 Such semiempirical methods have been used to study the retinal protonated Schiff base chromophore and its models. 343,378 Extensive benchmarking has already been performed by Thiel and co-workers for different semiempirical methods using PSB3 345 as well as for other organic chromophores. 379 The attractiveness of semiempirical methods comes from their low computational cost compared to ab initio quantum chemistry methods, which allows them to be used in larger systems. However, this is often at the cost of reduced accuracy or a requirement that the method be reparameterized for a particular system.
In the benchmark studies reviewed below, the Davidson corrected-MRCISD (MRCISD+Q) method is taken as the reference method due to its superior treatment of both static and dynamic contributions to the total electron correlation energy. To be more specific, we recall that, while the Hartree−Fock (HF) method does account for the electronic correlation due to the exchange interactions between electrons of parallel spins (Fermi correlation), the Coulomb correlation is completely missing. Therefore, the formal definition of electron (Coulomb) correlation energy is 380 where E exact is the exact energy. The exact energy can in theory be obtained at the full configuration interaction (FCI) level of theory, which is characterized by an electronic wave function incorporating all electronic excitations of all orders with respect to the closed-shell electronic configuration. E HF is the Hartree− Fock (HF) energy, which contains only a single closed-shell electronic configuration. Static electron correlation becomes very significant when a state has important contributions from more than one degenerate or nearly degenerate electron configuration. In such cases, the HF wave function is qualitatively wrong because it only accounts for one configuration (since it is  described by a single closed-shell Slater determinant). Static electron correlation can be accounted for by including the few specific electronic configurations that have major contributions to the wave function. For example, electron correlation would play an important role for Φ DIR in PSB3, which involves two electronic configurations that need to be treated on an equal footing ( Figure 43) as they become nearly degenerate during homolytic bond breaking. Dynamical electron correlation, on the other hand, is related to the instantaneous interaction between electron pairs and can only be accounted for by including many configurations, even those with relatively small weights, in the electronic wave function. Dynamic electron correlation would play an important role in both Φ DIR and Φ CT , but it plays a bigger role in Φ CT due to the fact that charge transfer involves a larger rearrangement of electronic density and therefore a larger perturbation with respect to the reference configuration. The different contribution of the dynamical electron correlation to the Φ DIR -dominated and Φ CT -dominated regions of the chromophore potential energy surface provides an important example of the differential electron correlation effect. This effect can modify both the shape and barriers on a single potential energy surface, or the energy gap between different potential energy surfaces (e.g., the excitation energies). This means that an unbalanced description of electron correlation could result in an incorrect description of barriers, potential energy surface shapes, gaps, and crossings. In the following points (i−vii), we briefly describe electronic structure methods benchmarked using PSB3, and the results of their benchmarks along the BLA and CI loop paths. Note that methods (i) and (iii) in particular are popular for the description of isomerization of π-conjugated systems and are mentioned in several instances throughout the review.

Chemical Reviews
(i) The multiconfigurational complete-active-space selfconsistent-field (CASSCF) method is an elegant solution to the missing static electron correlation. It allows the user to select certain electrons and orbitals, termed the "active space," and exactly accounts for their correlation with each other. However, the remaining electrons are left uncorrelated. The active space therefore must be selected carefully to include electrons and orbitals that are considered necessary for the description of the spectroscopic or reactivity properties under investigation. In other words, the selection of the active space is a function of the chemical/physical problem investigated. If a suitable active space is chosen, the CASSCF method should satisfactorily account for static electron correlation. However, the dynamic component is still missed (both for inactive electrons and partly for active space electrons because they are not excited to virtual orbitals).
In the BLA scan ( Figure 44A), this missing dynamic correlation causes CASSCF to have an unbalanced description of the S 0 and S 1 energy surfaces, yielding a destabilized chargetransfer state when compared to the MRCISD+Q result. On the other hand, CASSCF is able to correctly describe the shape of the potential energy surface along the BP and passes the CI loop test ( Figure 45A).
(ii) Multireference configuration interaction (MRCI) is a variational approach that uses the same principle as truncated configuration interaction methods, but instead of using HF as the reference zeroth-order wave function, it uses a multiconfigurational reference such as the CASSCF wave function. Thus, MRCISD (SD stands for single and double) would expand the reference CASSCF wave function by adding single and double excitations. Since the CASSCF wave function already includes the electron correlation within the active space, MRCI effectively introduces correlation of the electrons outside of the active space.
Note that MRCISD is not rigorously size extensive, meaning that its accuracy does not necessarily scale linearly with system size. A simple energy correction, the Davidson correction, 381 can be applied on top of MRCISD to estimate the energy of configuration interaction up to quadruple (indicated with the letter Q) excitations with almost no additional computational cost. This correction also improves the size extensivity of MRCISD. This leads to the MRCISD+Q energy that we use as a reference above, which is considered to yield values close to exact energies obtained with the unaffordable FCI method. For the benchmark studies reviewed here, MRCISD+Q serves as the reference method to measure the performance of other methods discussed here. However, MRCISD+Q is not used to study photobiological chromophores because it is too expensive for systems larger than PSB3.
The inclusion of dynamical electron correlation causes significant changes in the S 0 and S 1 energies along the BLA scan. Compared to CASSCF, both curves are stabilized with respect to the reference cis-PSB3 energy, but the change is disproportionate with the Φ CT -dominated state stabilized more than the Φ DIR -dominated state. This results in a shift in the geometry of the crossing (CI BLA ) along the BLA coordinate and, in turn, to a change in the entire topography of the S 0 potential energy surface, as shown in Figure 46. By topography, we refer to the relative tilt, pinch, and asymmetry of the surfaces forming the CI, which tells us about the shape of the intersection itself. 360 The CASSCF CI BLA is peaked, and therefore the ground-state potential energy surfaces hosts two different transition states mediating cis−trans isomerization. On the other hand, CI BLA in MRCISD+Q has a sloped topography, and therefore its ground state only hosts one transition state that mediates isomerization (TS CT ). Here we distinguish between topography and topology, the latter of which refers to the dimensionality of the corresponding BP and IS (which determines if the crossing is a conical or a linear crossing). MRCISD+Q does not yield a correct topology of the BP as the S 0 and S 1 curves cross twice along the CI loop ( Figure 45D). This is due to the Davidson correction, which is a simple energy correction term not correcting the electronic state coupling. MRCISD, on the other hand, is able to describe the BP correctly like CASSCF ( Figure 45C).
(iii) Similar to MRCI, multireference second-order perturbation theory (MR-PT2) approaches employ a multiconfigurational reference zeroth-order wave function. However, they account for dynamic electron correlation perturbatively instead of using a variational scheme. This makes MR-PT2 methods computationally less expensive than MRCISD, although they may not be as reliable. The CASPT2, 382 MS-CASPT2, 383 XMCQDPT2, 118 and QD-NEVPT2 279,384,385 are all examples of MR-PT2 methods based on a CASSCF zeroth-order wave function. However, due to their different formulations, they behave differently along the BLA coordinate. As shown in Figure  44A, CASPT2 and MS-CASPT2 overstabilize the Φ DIRdominated state, QD-NEVPT2 overstabilizes the Φ CT -dominated state, while XMCQDPT2 slightly overstabilizes both states (the Φ DIR state more than the Φ CT state) with respect to MRCISD+Q. The performance of all four methods could be improved by introducing small modifications (see Figure 44B). For instance, the use of a shift that modifies the energies of the active orbitals in CASPT2 and MS-CASPT2 (the IPEA shift 386 ) significantly improves the agreement with MRCISD+Q. For QD-NEVPT2, a similar improvement is obtained by enlarging the active space of the reference CASSCF wave function by including polarizing π*-orbitals. 90,387 Finally, the agreement of XMCQDPT2 could be improved by using a modified operator that incorporates terms arising due to the nonseparable part Γ ns of the CASSCF state-averaged second-order density matrix Γ. 90 Overall, MR-PT2 methods show a significant improvement over CASSCF. However, it is difficult to know if the modifications mentioned above are generally applicable. For instance, the use of the IPEA shift in CASPT2 is a topic of active debate in the computational chemistry community. 388,389 Also note that while CASSCF is size-extensive (if certain conditions are met), and perturbation theory is formally also size-extensive, MR-PT2 is not rigorously size-extensive. 390 Despite their drawbacks, MR-PT2 methods are widely used in computational photochemistry and photobiology because of their affordable cost coupled with the ability to capture a large portion of electron correlation. 391 Unfortunately, MR-PT2 geometry optimizations are still not practical for large systems (e.g., for biological chromophores discussed in this review). For this reason, a common protocol has been to obtain CASSCF equilibrium structures and then compute the potential energies via MR-PT2 calculations. Such MR-PT2//CASSCF protocols have been widely employed for evaluating the energies of QM/ MM models, as discussed in the next subsection. The CI loop results show that both the regular CASPT2 (SS-CASPT2) and MS-CASPT2 display artifacts near a CI point. CASPT2 suffers from its inability to describe the derivative coupling and therefore yields a linear crossing. MS-CASPT2 corrects for this but behaves erratically when the underlying CASSCF wave function has degeneracies. 90,118,182 MCQDPT2 and NEVPT2 had similar problems, but were remedied using the extended version of MCQDPT2 (XMCQDPT2 118 ) and the quasi-degenerate implementation of NEVPT2 (QD-NEVPT2), respectively. 279,384,385 An extended implementation of MS-CASPT2, termed XMS-CASPT2, has been developed 281 and would be expected to behave correctly in the vicinity of intersections.
(iv) Spectroscopy oriented configuration interaction (SORCI) is effectively a multireference method that combines elements of MRCI and MRPT2 to achieve a balance between accuracy and affordability in a black-box framework. 392 SORCI uses several mathematical approximations such as the resolution of identity and different levels of user-controlled truncations. While SORCI has many of the benefits of multireference methods and can calculate ground and excited states on an equal footing, it is aimed at vertical energy difference calculations (hence the "spectroscopically oriented" title) and not for potential energy surface mapping, particularly because the use of cutoffs results in an unsmooth surface.
Nonetheless, SORCI was tested along the BLA path ( Figure  44C) and found to be capable of producing relatively smooth profiles, with the exception of regions in close proximity to the CI. These artifacts near the CI occur due to the Davidson correction 381 used within SORCI. 392 SORCI captures the relative energies of the Φ DIR and Φ CT states correctly with respect to each other and therefore gives a correct geometry for CI BLA . However, both states are overstabilized relative to the cis-PSB3 geometry. 355 While SORCI was not tested along the CI loop, it would not be expected to yield a conical crossing due to its use of the Davidson correction similar to MRCISD+Q.
(v) Equation-of-motion coupled-cluster (EOM-CC) is a multistate extension of coupled-cluster that allows the description of a variety of target states. Coupled-cluster methods are currently among the most robust approaches for solving electronic structure problems. The exponential ansatz used for coupled-cluster wave functions ensures size extensivity (i.e., ensures that the quality of the calculation does not deteriorate with increasing size of the system, which is a common problem encountered in some electronic structure methods). Furthermore, these methods are relatively black box and can be systematically improved by increasing the level of excitations. The EOM-CC extension to coupled-cluster allows the description of a variety of target states in a single-reference formalism (unlike multireference methods such as MRCI and MR-PT2 that employ a multiconfugurational zeroth-order wave function and then introduce dynamical electron correlation in a second step). The choice of excitation operator used in EOM-CC allows this method to be used to describe different kinds of target states (see Figure 47). 393 For instance, different flavors of EOM-CC allow the calculation of excited states (EOM-EE), 394 ionized states (EOM-IP), 395 electron-attached states (EOM-EA), 396 or spin-flip states (EOM-SF). The spin-flip approach is a powerful way to treat closed-shell and open-shell, as well as doubly excited configurations, on an equal footing if they can all be described as a single excitation from a reference triplet state. 397−400 In ref 354, a number of EOM-CC approaches were tested, including both EOM-EE and EOM-SF. In Figure 44C we show the most successful, which was able to reproduce the MRCISD +Q potential energy profile closely. This required using the spinflip variant of EOM-CC and the inclusion of perturbative triples (dT) corrections on top of single and double excitations. However, the EOM-SF-CCSD(dT) approach was unable to describe the topology of the BP near CI BLA correctly ( Figure  45D). This is most likely due to the triples correction, which employs EOM-SF-CCSD as a reference and corrects the energy to account for triple excitations and therefore behaves similarly to any other energy correction term (e.g., the Davidson correction in MRCISD+Q), yielding an incorrect crossing topology. Unperturbed EOM-SF-CCSD, on the other hand, would be expected to yield a BP spanning three dimensions, due to the nonhermiticity of the coupled-cluster ansatz. 363,364 (vi) The Algebraic-diagrammatic-construction (ADC) method is a polarization propagator approach based on many-body Green's function theory. 401−403 It describes the effect of an external perturbation on electronic structure using Møller− Plesset perturbation theory. The ground state therefore is described at the corresponding level of Møller−Plesset theory [i.e., MP2 for ADC(2) and MP3 for ADC (3)]. Figure 44C shows that ADC (2) gives energies along the BLA that are in very good agreement with MRCISD+Q. However, because of the formulation of ADC methods, no coupling exists between the ground MP and excited ADC states, yielding a linear crossing. 356 (vii) Time-dependent density function theory (TD-DFT) is one of the most widely used approaches in spectroscopy. Its most attractive feature is its ability to account for electron correlation and its relatively low computational cost scaling, which allows it to be used in much larger systems than would be possible with any of the electronic structure methods discussed above. However, there are several drawbacks for the use of TD-DFT in photochemistry. There are a number of well-documented problems with traditional TD-DFT methods, 335,404 and particularly in the description of charge-transfer states, 405 bond-breaking, 406 and in being able to describe in a balanced way states with double-excitation character. 407 Some of these problems are mitigated by using hybrid or long-range separated functionals (to address problems with describing charge-transfer states) or double-hybrid functionals (to better describe states with double excitation). However, due to the large number of different DFT functionals being developed, it is not clear from the outset which functional is best for a given problem. A growing number of benchmarks have helped to identify the functional more suitable for a given class of systems or applications. 407,408 Huix Rotllant et al. have also benchmarked a number of different DFT functionals along the BLA path of PSB3, including hybrid functionals, hybrid meta-GGA functionals, rangeseparated hybrid functionals, and double-hybrid functionals. 288 It was found that the relative energies of the Φ CT and Φ DIR states are sensitive to the type of functional, and the choice of functional is not always obvious. In Figure 44C, we show only one example corresponding to the double-hybrid mPW2PLYP functional. Similarly, Casida, Filippi, and co-workers have benchmarked TD-DFT along torsional coordinates (both double-bond and singlebond isomerization) in reduced rPSB11 models and found that different functionals may yield a qualitatively different potential energy surface along these coordinates. 349 A second problem with TD-DFT is that it is effectively a singlereference method, which may result in the calculation not converging in regions of degeneracy or near-degeneracy. In the case of the BLA scan, linear-response TD-DFT would often break down near the intersection of the Φ CT and Φ DIR states. This is only somewhat mitigated by using a double-hybrid functional.
A third problem with traditional TD-DFT approaches is their inability to account for coupling between the S 0 and S 1 state, resulting in a linear crossing (as revealed by the CI loop in Figure  45A) instead of a conical crossing. 182 Such artifacts in the case of TD-DFT has been explored in detail by Levine, Martı́nez, and co-workers. 181 Despite these drawbacks, TD-DFT remains a powerful approach in photochemistry and photobiology, particularly for calculating spectroscopic properties and even excited-state processes. Attempts are being made to extend DFT methods to be able to treat nonadiabatic regions of the potential energy surface correctly. Examples of this include configurationinteraction corrected Tamm-Dancoff approximation (CIC-TDA) TD-DFT, 285 constrained density functional theory− configuration interaction (CDFT-CI), and the state-interaction state-averaged spin-restricted ensemble-referenced Kohn−Sham (SI-SA-REKS) approach based on ensemble DFT. 288,409 Another very promising approach is the application of the spin-flip method to TD-DFT. 410 SF-TD-DFT is capable of yielding the correct 2D topology of the branching space around the conical intersection 182,289 However, SF-TD-DFT often suffers from severe problems with spin-contamination, especially far from the Franck−Condon region. A spin-adapted variant of SF-TD-DFT (SA-SF-TD-DFT) has recently been developed that resolves the issue with spin-contamination. Development of analytical gradients may allow this method to be used affordably to run ab initio excited-state or multistate molecular dynamics. 290 As interest in photochemical and photobiological problems continues to grow, there is increasing demand for accurate yet affordable electronic structure methods that can treat extended regions of the ground and excited-state potential energy surfaces. Electronic structure methods have made considerable progress in recent years. However, remarkably, classical methods that have existed for over 20 years, particularly those building on multiconfigurational methods and coupled-cluster approaches, remain among the most widely used ab initio methods for photobiology today. While coupled cluster methods are mainly used for highly accurate studies of small molecules, it has been possible to employ multiconfugurational and multireference approaches to larger systems such as biological chromophores, particularly through the use of the MR-PT2//CASSCF protocol where geometries are optimized at the CASSCF level of theory and energies obtained at the MR-PT2 level of theory. 391 Many of the studies we discuss in the present review employ the MR-PT2//CASSCF method. The reason is that CASSCF can capture multiconfigurational character, which is often a necessity for describing photochemical processes such as bond-breaking and isomerization and can also capture the correct topology of conical intersections, allowing it to be used for nonadiabatic trajectory calculations. At the same time, to capture energetics correctly, it is necessary to account for dynamical electron correlation, which can be achieved relatively affordably with MR-PT2.
On the other hand, while the CASPT2//CASSCF protocol is used routinely for systems such as rPSB11 and chromophores of similar size, it requires the user to be knowledgeable about the drawbacks of CASSCF and CASPT2 and therefore is certainly not a black-box approach. This protocol may also become too expensive when looking at larger chromophores, in particular those with more electrons that need be included in the active space. In particular, this might apply to chromophores with a larger π-conjugated framework or with heteroatoms with lone pairs that may be involved in the photochemistry. Since the computational cost of CASSCF scales factorially with the size of active space, large active spaces become rapidly intractable. It is worth mentioning here that two families of methods are being developed to overcome this shortcoming. The first approach is the restricted active space (RAS) method 301 that divides the total active space into three subspaces. One subspace is similar to CAS and allows all possible excitations, while the other two are typically restricted to single and double excitations from one subspace to another. A more general implementation was achieved by developing the general active space (GAS) method where the number of such spaces is arbitrary. 302 The second strategy is based on the density matrix renormalization group (DMRG) in combination with CASSCF. 411 This method allows an active space of up to 50 electrons and orbitals. Several programs are already available to do such calculations (e.g.. BLOCK by the Chan group, 412,413 QCMaquis by the Reiher group, 414,415 and CheMPS2 by Wouters et al. 416 ). Furthermore, DMRG can be combined with MR-PT2 methods such as NEVPT2 or CASPT2 to yield more quantitative results. 303,304,417 However, while a large size active space is becoming feasible in a foreseeable future, this also means that the selection of a relevant active space will be more difficult. For this purpose, the Reiher group has developed an algorithm that takes over the task for the user. 418

I.2. Hybrid Quantum Mechanics/Molecular Mechanics Approach
The computational investigation of light-responsive proteins necessitates an approach that can treat the chromophore in the excited state and account for the protein and/or solvent environment as well. This is most commonly achieved by using a QM/MM model. 419 The QM/MM approach belongs to a family of multiscale methods, which are suitable for the investigation of the properties of a large molecular system. Multiscale methods employ two or more levels of theory with distinct accuracy and computational cost. Specifically, in QM/ MM, the use of an accurate, but relatively expensive, quantum mechanical (QM) level of theory for treating the group involved in electronic excitation and less accurate, but affordable, molecular mechanical (MM) approach for the remaining larger part makes it feasible to calculate protein-chromophore complexes. The key steps in constructing a basic two level QM/MM model are (i) partitioning the protein into QM and MM subsystems and (ii) choosing the scheme that describes the interaction (i.e., coupling) between the two subsystems. Moreover, a special situation arises in cases where the boundary between the QM and MM subsystems is located at one or more covalent bonds, as in the example of rhodopsin mentioned in the Introduction. 420 In such cases, a suitable scheme describing the covalent link between the two subsystems needs to be adopted as well.
QM/MM models of photoreceptor proteins have reached an evident level of maturity and are capable of treating the absorption and emission spectra of different types of proteins (e.g., with backbones dominated by α-helics and β-sheet structures), incorporating different chromophores (e.g., cationic, anionic, and neutral chromophores). 409 The CASPT2// CASSCF/MM approach in particular has been used to study a number of different photoactive systems as shown in Figures 48  and 49.  I.2.1. QM/MM Embedding Schemes. Besides the accuracy of the electronic structure method used to describe the QM subsystem (see section I.1), the QM and MM subsystems coupling is of central importance. For a detailed discussion of different QM/MM embedding schemes, we refer the reader to ref 424). There are three main coupling schemes. (1) Mechanical embedding: The QM subsystem does not "feel" the electrostatic field generated by the point charges of the MM subsystem. The only nonbonding interaction coupling the two subsystems is the short-range van der Waals interaction. (2) Electrostatic embedding: in addition to the mechanical embedding, the QM subsystem is also perturbed by the point charges of the MM subsystem. This interaction can lead to differential effect on different electronic states that might result in stabilization or destabilization of the ground and excited states as shown in Figure 6. However, only the QM subsystem experiences an electrostatic effect. The MM subsystem is described by point charges whose values are not affected by the QM subsystem electron density. (3) Polarizable embedding: for a mutual polarization between the QM and MM subsystems, the MM environment should also respond to the changes in the electron density of the QM region. This requires a MM subsystem that features a polarizable charge model.
Although the polarizable embedding scheme is, in principle, the most accurate, it is rarely used. The reason is that the few available polarizable force fields are still being benchmarked and developed. 425,426 In practice, presently, the electrostatic embedding is the most commonly used QM/MM coupling scheme. Although recent attempts to use polarizable embedding have been reported, their application in the simulation of ultrafast dynamics is limited. 427 The same holds also for QM/DFT embedding, which is only available for single point calculations at the moment. 428 The inherent intricacy of excited-state dynamics studies is also the reason why complex QM/MM models (e.g., large membrane embedded photoreceptor models) are mostly used for calculations of vertical excitation energies. 423 In some cases, the order of the electronic states is perturbed so much by the electrostatic embedding that a subsequent molecular dynamics or reaction path mapping is not possible due to the presence of several intruder electronic states between the spectroscopic state and the ground state. If these intruder states are physical (i.e., they would be experimentally observed), one can, in principle, correctly account for them in a nonadiabatic trajectory computaions (see section I.3) or, when using an ad hoc surface hopping strategy, 429 in a MEP computation. However, certain intruder electronic states are unphysical as they would be removed at higher levels of theory (e.g., CASSCF generated intruder states may be removed at the CASPT2 level). In these cases, certain authors have switched to a simpler mechanical embedding that avoids these technical problems and allows one to obtain a qualitative description of the system reactivity. 29  The three energy terms arise from one QM calculation of the inner system and two MM calculations of the total (inner + outer) system and the inner system, respectively. The last term is subtracted from the total MM energy, thus eliminating double counting. This QM/MM energy expression was coined subtractive. The main advantage of the subtractive scheme is the avoidance of explicitly calculating the QM-MM coupling term, and therefore, regular QM and MM programs can be used without any modifications. The first realization of the subtractive scheme was named integrated molecular orbital/molecular mechanics (IMOMM). 436 Further development extended this method to combine a high level QM method with a lower-level one (IMOMO), 434 eventually leading to a fully generalized multilayer scheme called ONIOM. 433,437 I.2.3. Link Atom Scheme. Partitioning a molecular system into two subsystems is trivial for a solute surrounded by solvent molecules. In such a scenario, there is no covalent bond that connects the QM and MM subsystems. However, in many photoreceptor proteins the chromophore is covalently bound to the apoprotein. Hence a QM/MM model has to deal with bonds at the boundary between the subsystems. In the following, we discuss the case of Rh, where the rPSB11 chromophore is covalently bound to the Lys296 residue (see Figure 50). Since we are concerned with the photochemistry of rhodopsin, the QM subsystem needs to include the whole conjugated π-electron system. Hence, the boundary must be located on a bond of the Lys296 side-chain beyond the Schiff base functional group. The C δ -C ε was selected to minimize the computational cost without altering significantly the electronic structure of the QM subsystem. However, the bond has to be cleaved homolytically or heterolytically, leading to undesirable changes in electronic structure. Three different approaches have been proposed to overcome this problem: (1) link atom schemes, 431 (2) boundary atoms, 438 and (3) localized orbital schemes. 439 The studies presented in this review make use of the link-atom scheme exclusively, so we will limit our discussion to this particular Figure 50. A link atom scheme for Rh. The rPSB11 chromophore is covalently bound to the lysine residue Lys296 (in blue) embedded in the protein cavity (see also Figure 49B as a reference). The common partitioning puts the boundary on the C δ -C ε single bond. The hydrogen link atom (shown in orange) is placed along this bond. scheme and refer the interested reader to the review by Senn and Thiel 424 for a description of the other two schemes.
In the link-atom scheme, a new atom is introduced to cap the cleaved bond. In most cases it is a hydrogen atom, and the bond that is truncated is selected to be a single bond of small or negligible polarity (like the Lys296 C δ -C ε bond in bovine rhodopsin). The introduction of a new atom that is not present in the real system introduces a few problems, however. One problem is the addition of three degrees of freedom that an atom has and that needs to be optimized in a geometry optimization or that can introduce some artifacts in trajectory computations. This disadvantage was resolved by constraining the bond between the QM frontier atom and the link atom to be aligned along the MM bond (see Figure 47) connecting the QM and MM subsystems. 440 Such a constraint removes all additional degrees of freedom except the link atom bond length. Another caveat that results from the link atom scheme is the overpolarization of the link atom by the adjacent MM frontier atom. If not addressed, this could lead to large artifacts that even lead to an incorrect ordering of excited states. The most common solution to this problem is the charge-shifting scheme that was put forward by Sherwood and co-workers. 441 In this approach, the charge of the closest MM frontier atom is distributed onto neighboring MM atoms bound to it. For instance, the C δ charge (e.g., the AMBER point change) of the Lys296 in Rh is set to zero and its value is evenly added to the remaining MM atom of the Lys residue (leaving the total residue charge unchanged).
I.2.4. Protocols for QM/MM Model Building. There are many ways to construct QM/MM models, so here we will discuss one example to illustrate how such a model construction could be done. A photoreceptor QM/MM model is often derived from a crystallographic or nuclear magnetic resonance (NMR) structure deposited in the Protein Data Bank, 442 or from a homology model. Usually, these structures do not include protons, may be missing a number of residues, and do not contain information on the state of ionizable residues or on the distribution of nearby counterions. Therefore, such missing parts must be incorporated into the model using some QM/MM model-building workflow (see Figure 51). This may require structural manipulations, including the relaxation of the system after addition of the hydrogen to the crystal structure. This relaxation is usually carried out at the MM level of theory via geometry optimization and/or molecular dynamics (MD) simulations. There is no universal way to do this, and the result depends on the protocol, the software, and the parametrization of the chromophore used. This often results in different laboratories having slightly different models. Furthermore, to have a more complete model, the entire protein may be placed in a solvent (or a membrane patch if it is a membrane protein). Finally, it may be important to add ions in the simulation to neutralize the charges of the protein. Once the classical model is completed, a minimization and a subsequent MD trajectory are often necessary. Upon equilibration, such a trajectory could be used to generate one single QM/MM model, or different snapshots taken at different times can be used as starting points to generate a series of QM/MM models. The QM/MM model(s) can then be optimized at the QM/MM level of theory and used for subsequent excitation energy, reaction path, or trajectory calculations.
Another option for generating a set of QM/MM models that can be used to run an ensemble of trajectory calculations is to use a sampling protocol to generate initial geometries and velocities according to a specific distribution such as the Wigner or the Boltzmann distribution. In this way, one has the opportunity to investigate, by running a group of trajectories in parallel, the dynamics of an entire population and even compute photoisomerization QYs. As discussed below, the accuracy of these computations will depend on the number of trajectories run. Due to the computational cost limitations, QM/MM trajectories computed at the CASSCF/MM level are usually limited to lifetimes that do not exceed a few picoseconds. However, this is enough for investigating the light-driven primary event in the ultrafast activation of biological photoreceptors.
Recently, there have been attempts to automate and, partially, standardize the building of QM/MM models for the purpose of comparison of mutants or of photoreceptors of the same family. 332 Presently, these are monomer QM/MM models of retinal proteins whose construction only requires an experimental structure or homology model and information about the protonation states of the ionizable residues and ions that counterbalance the charges. The so-called automatic rhodopsin protocol (ARM) generates a QM/MM model consisting of a frozen protein environment with a flexible protein binding pocket formed by the side-chains of the chromophore cavity (i.e., consistent with the model in Figure 49B and 50). The quality of the models has been assessed by comparing calculated and observed excitation energies for a series of rhodopsins from different organisms. As detailed in section 8.2.4, it is shown that ARM models are capable of providing information on excitation energy trends.

I.3. Nonadiabatic Molecular Dynamics Simulations and Trajectory Surface Hopping
In the course of photoisomerization, a molecule progresses on the excited state potential energy surface(s) and ultimately decays to the S 0 state before completing the isomerization event or regenerating the reactant. The simulation of photoisomerization dynamics therefore requires the description of transitions between the potential energy surfaces associated with different electronic states, also referred to as hopping events. In the context of photoisomerization, such a transition occurs in the proximity of CIs located along the IS, where the adiabatic (Born−Oppenheimer) approximation breaks down.
Such nonadiabatic transitions can be described in several different ways, ranging from full quantum dynamics (where the nuclear motion is treated quantum mechanically) to semiclassical methods (where nuclei are treated classically but the surface hop is treated in a special way to account for nonadiabatic events). 444 Subotnik et al. have recently examined popular surface hopping algorithms and presented modern advances in such algorithms. 445 Here, in the context of photobiology, we will focus our discussion on three trajectory surface hopping (TSH) algorithms commonly used for systems the size of biological chromophores.
One approximate algorithm implemented by Robb and coworkers 136,337,340,446 describes the coupling between two states, k and j, using the numerical approximation: where Ψ k and Ψ j are the time-dependent electronic wave function of the states, t 0 the initial time, and Δt the time-step. The coupling is used to estimate the interaction of the two states and to decide if a hop between the states occurs ( Figure 52). A transition between Ψ k and Ψ j occurs if a specified threshold is reached.
After a hopping event, the energy is conserved by correcting the momentum in the direction parallel to the derivative coupling vector. 340,447,448 The second is the Tully's fewest switches algorithm. 449 This approach requires the coupled propagation of both the nuclei and the electrons motion. As for the previous method, the former is performed by integrating Newton's equations of motion for the nuclei moving on a given potential energy surface on the femtosecond time scale, while the latter is carried out by integrating the time-dependent Schrodinger equation: on the attosecond time scale. R stands for the nuclear position (representing the molecule geometry), r stands for the electronic coordinates, d kj = ⟨Ψ k (r;R)|∇ R Ψ j (r;R)⟩ are the derivative coupling terms, and V kj = ⟨Ψ k (r;R)|H(r;R)Ψ j (r;R)⟩ are the electronic Hamiltonian elements for a given geometry R. The wave function Ψ(r, R, t), is expanded as a linear combination of stationary states ϕ j (r;R), yielding the time-dependent expression Ψ(r, R, t) = ∑ j c j (t)ϕ j (r;R). The d kj terms are key parameters for inducing the surface hopping as they promote the mixing of the electronic states.
The equation above is integrated by computing the Rḋ kj terms. By applying the chain rule, the gradient of the wave function is transformed into a wave function time-derivative, yieldinġ . These terms can be calculated analytically, but there is also a possibility to determine them numerically along the trajectory. It is required that at each step of the nuclear motion the integration step defined by Δt coincides with the exact calculation of the wave function. However, the time-step of the electronic integration is typically at least 1 or 2 orders of magnitude smaller than Δt. Therefore, to get accurate results in the integration of the time-dependent Schrodinger equation, an interpolation between consecutive values of Ṙ· d kj (R) is required. An interpolation/extrapolation scheme for the determination of time-derivative coupling terms has been put forward by Hammes-Schiffer and Tully. 450 Additionally, in order to preserve the consistency of the population density matrix, a decoherence correction proposed by Granucci and Persico 451 is required. By integrating the equation above one obtains, for each time step, the populations of each electronic state |c i | 2 . Depending on the time variation of these populations, the fewest switches algorithm 449 stochastically decides whether the trajectory should continue to propagate on the same potential energy surface (e.g., the S 1 potential energy surface) or should jump onto another one (e.g., to the S 0 potential energy surface). After the hop, the next energy gradient for the nuclear motion is calculated using the potential energy surface associated with the new state while the nuclear velocities are scaled to preserve the total energy.
The third algorithm is based on the Zhu-Nakamura theory. 452 The surface hopping approach based on the Zhu-Nakamura theory provides a way to compute the probability for a nonadiabatic transition for a wide spectrum of energy gaps. 453 This approach is computational feasible and straightforward to implement because it can be also applied if the derivative coupling terms are not available. 454 However, it has several important advantages compared to some of the other surface hopping approaches, such as the treatment of surface hops that are energetically forbidden. The Zhu-Nakamura surface hopping scheme has been applied to several one-dimensional problems before it was recently extended for multidimensional problems. 455 The key step in the Zhu-Nakamura approach is based on the assumption that hops can take place in the close vicinity to a CI and that the derivative coupling direction can be approximated by the eigenvector of the Hessian of the potential energy surface driving the dynamics ∂ 2 ΔV/∂R i ∂R j corresponding to its maximum eigenvalue.
As mentioned earlier, the surface-hopping methods presented above belong to the group of semiclassical methods where the electrons are treated quantum mechanically and the nuclei are described classically. However, due to the breakdown of the Born−Oppenheimer approximation in the vicinity of conical intersections, the quantum nature of the nuclear motion can become important. For example, quantum effects such as quantum decoherence, zero-point energy effects, tunneling, geometric phase effects, and superexchange cannot be described classically. As a result, there have been numerous attempts at extending classical (e.g., surface-hopping) approaches to account for some quantum effects (e.g., see refs 456−460).
While semiclassical approaches are affordable and work well for many types of problems, their extension to treating quantum effects is not without difficulty. Therefore, full-dimensional quantum wavepacket propagation methods are still widely Figure 52. Change in the adiabatic electronic wave functions for two consecutive steps of a semiclassical surface-hop trajectory calculation. The first time-point is marked by t 0 and the time-step by Δt. Ψ k and Ψ j are time dependent electronic wave function of the k and j states, while ϕ k and ϕ j denote stationary diabatic states (i.e., one can think of ϕ k and ϕ j as reference, time-independent electronic configurations whose mixing provide a description of the adiabatic states. The change in mixing coeficients at each time step makes possible the description of the time evolution of the adiabatic states). Adopted from ref 273. Copyright 2016 American Chemical Society. applied to accurately study dynamics of small molecules with relatively few degrees of freedom. Examples of such approaches includes multiconfiguration time-dependent Hartree (MCTDH) and variational multiconfiguration Gaussian wavepackets (vMCG) 461,462 (see ref 463 for a discussion of these methods). However, such methods are not usually applied for systems the size of biological chromophores reviewed in this work due to the exponential growth of their computational cost with increasing number of degrees of nuclear freedom. However, progress in this direction, such as the development of multilayer techniques, 464 is allowing MCTDH to be applicable to larger systems. For example, multilayer MCTDH was recently applied to several hundreds of degrees of freedom to study the FMO complex. 465 Similar advancements have been reported for the vMCG method. 466 One approach that combines elements of semiclassical and quantum mechanics is the ab initio multiple-spawning (AIMS) approach developed by Martı́nez and co-workers. 467,468 In this method, a set of Gaussians, simulating a nuclear vibrational wavepacket, follow classical trajectories while their expansion coefficients are propagated according to the time-dependent Schrodinger equation. Nonadiabatic transitions are then described by expanding the basis set of the Gaussian functions after the transition, leading to "spawning" of new populations that may either continue to evolve on the excited state or decay to the ground state. The AIMS approach improves on semiclassical methods by describing nonadiabatic transitions quantum mechanically, including coherence and decoherence effects. However, the AIMS computational cost may increase in cases of many nonadiabatic transitions (i.e., depending on the topology of the potential energy surface) and can become intractable for extended molecular systems. One workaround for the rapidly increasing number of spawned function is stochastic selection, 469 where the number of basis function can be reduced for efficiency. In general, the AIMS method has been constantly extended over the past decade by introducing full ab initio spawning that is free from any empirical parameters 468 and recently including the treatment of intersystem crossings. 470

II. TECHNICAL SECTION: RELATIONSHIP BETWEEN
REACTION COORDINATE AND INTERSECTION SPACE The MEP coordinate described in section I.1 is expected to operate in all rPSB chromophores. Since the MEP is dominated by three modes (i.e., BLA, HOOP, and torsional coordinate), it cannot be realistically visualized in 2D. This has forced us to illustrate the relationship between S 1 MEPs, S 0 MEPs, and the CI only schematically as seen in Figure 41C. In such a 2D representation, the relationship between an IS segment (containing the CI CIS , CI BLA , and CI TRANS ) cannot be understood. In fact, the IS appears parallel to the mixed torsion and HOOP mode usually associated with the isomerization coordinate. In oder to clarify the relationship between S 1 MEPs and the IS, in this section we present such a detailed mapping for the case of the minimal three double-bond model of rPSB, PSB3 introduced in section I. 1.
It has been possible to map the evolution of PSB3 along its S 1 MEPs in a 3D coordinate space. 62 The 3D schematic representation of the relationship between the BLA, HOOP, and torsional coordinate change in PSB3 is given in Figure 53 for both cis to trans (starting from PSB3 CIS ) and trans to cis (starting from PSB3 TRANS ) isomerization of the central C2C3 double bond. In the same figure, we show that the segments of the S 1 MEPs hitting the CIs is a combination of skeletal double-bond torsion and HOOP while the BLA value remains substantially constant in this region of the MEP. On the other hand, a change in the BLA value would lift the S 1 /S 0 energy degeneracy at the CI. Therefore, the MEP and BLA coordinates have components belonging to the two-dimensional BP seen in Figure 42. In contrast, a combination of skeletal double-bond torsion and HOOP opposite to the MEP hitting the CI would span the IS curve. As anticipated in Figure 4 and related text of the Introduction, the IS is a N-2-dimensional space (i.e., a hyperline with respect to the N-dimensional space of the full set of independent vibrational degrees of freedom of the system). The IS is orthogonal to the BP and consists entirely of CI points. S 1 MEPs launched from different stereoisomers reach different CIs along the IS, each mediating different relaxations. This is illustrated in Figure 53, where excited-state MEP calculations launched from PSB3 CIS and PSB3 TRANS reach the distinct CI TRANS and CI CIS located along the IS. CI TRANS and CI CIS are not minimum energy CI points (i.e., energy minima along the IS). In fact, the search for a minimum energy CI 62 leads to CI BLA , which is located between CI TRANS and CI CIS and mainly characterized by a different value of the HOOP coordinate. More specifically, CI BLA features the carbon atoms of the isomerizing double bond in a substantially sp 2 -hybridized state and, therefore, has a nearly zero HOOP value. Note that Figure  53 is a representation of the MEP CIS , MEP TRANS , and CI seam pathways displayed in Figure 43, although in coordinate space.
The geometrical progression described by the MEP is consistent with that of the S 1 trajectory of Figure 4. More specifically (follow, for instance, the PSB3 CIS progression toward CI CIS ), the geometrical changes are driven by the need to rapidly decrease the π-overlap across the isomerizing double bond. This happens because, in the FC region, S 1 has a Φ CT character setting the orbital interaction across the bond into an antibonding state. Therefore, the fastest and most effective way to decrease the πoverlap is the expansion of the C2−C3 bond length. This is immediately followed by the HOOP deformation and finally by the skeletal twisting. The initial change in HOOP has the effect of changing the hybridization status of the two carbon atoms of the isomerizing double bond. This is the main reason why at both CI TRANS and CI CIS such carbons are largely sp 3 hybridized. The hybridization then changes as one moves to CI BLA , where the same atoms have sp 2 hybridization again.

ASSOCIATED CONTENT Special Issue Paper
This paper is an additional review for Chem. Rev. 2017, 117, issue 16, "Ultrafast Processes in Chemistry".

Author Contributions
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Notes
The authors declare no competing financial interest. where he is in charge of the Laboratory for Computational Photochemistry and Photobiology. In 2015, he was appointed USIAS Fellow at the University of Strasbourg, France. His interests span both methodological and applicative aspects of computational photochemistry and photobiology; two research fields that he contributed to shape.