# Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores

Madushanka Manathunga, Xuchun Yang, Hoi Ling Luk, Samer Gozem, Luis

Manuel Frutos, Alessio Valentini, Nicolas Ferré, Massimo Olivucci

► To cite this version:

Madushanka Manathunga, Xuchun Yang, Hoi Ling Luk, Samer Gozem, Luis Manuel Frutos, et al. Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores. Journal of Chemical Theory and Computation, 2016, 12 (2), pp.839 - 850. 10.1021/acs.jctc.5b00945. [hal-01409056](https://hal HAL Id: hal-01409056)

Submitted on 10 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Probing the Photodynamics of Rhodopsins with
Reduced Retinal Chromophores

Madushanka Manathunga, Xuchun Yang, Hoi Ling Luk, Samer Gozem, Luis

Manuel Frutos, Alessio Valentini, Nicolas Ferré, Massimo Olivucci

► To cite this version:

Madushanka Manathunga, Xuchun Yang, Hoi Ling Luk, Samer Gozem, Luis Manuel Frutos, et al. Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores. Journal of Chemical Theory and Computation, 2016, 12 (2), pp.839 - 850. 10.1021/acs.jctc.5b00945. hal-01409056

HAL Id: hal-01409056

https://amu.hal.science/hal-01409056v1

Submitted on 10 Jan 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores

This is the peer reviewed version of the following article:

Original:

Manathunga, M., Yang, X., Luk, H.L., Gozem, S., Frutos, L.M., Valentini, A., et al. (2016). Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores. JOURNAL OF CHEMICAL THEORY AND COMPUTATION, 12(2), 839-850 [10.1021/acs.jctc.5b00945].

Availability:

This version is available http://hdl.handle.net/11365/984822 since 2021-04-13T19:06:49Z

Published:

DOI:10.1021/acs.jctc.5b00945

Terms of use:

Open Access

The terms and conditions for the reuse of this version of the manuscript are specified in the publishing
policy. Works made available under a Creative Commons license can be used according to the terms and
conditions of said license.

For all terms of use and more information see the publisher's website.

(Article begins on next page)
Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores

This is the peer reviewed version of the following article:

Original:

Manathunga, M., Yang, X., Luk, H.L., Gozem, S., Frutos, L.M., Valentini, A., et al. (2016). Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores. JOURNAL OF CHEMICAL THEORY AND COMPUTATION, 12(2), 839-850 [10.1021/acs.jctc.5b00945].

Availability:

This version is available http://hdl.handle.net/11365/984822 since 2021-04-13T19:06:49Z

Published:

DOI:10.1021/acs.jctc.5b00945

Terms of use:

Open Access

The terms and conditions for the reuse of this version of the manuscript are specified in the publishing
policy. Works made available under a Creative Commons license can be used according to the terms and
conditions of said license.

For all terms of use and more information see the publisher's website.

(Article begins on next page)
# Probing the Photodynamics of Rhodopsins with Reduced Retinal Chromophores

Madushanka Manathunga†,%, Xuchun Yang†,%, Hoi Ling Luk†, Samer Gozem†, Luis Manuel Fru-tos#, Alessio Valentini¶,#, Nicolas Ferrè% and Massimo Olivucci†,¶,‡

†Department of Chemistry, Bowling Green State University, Bowling Green, OH 43403, USA, ¶Dipartimento di Bio-tecnologie, Chimica e Farmacia, Università di Siena, via A. Moro 2, I-53100 Siena, Italy, #Departamento de Química Física, Universidad de Alcalá, E-28871 Alcalá de Henares, Madrid, Spain, %Aix-Marseille Université, CNRS, Institut de Chimie Radicalaire, 13397 Marseille Cedex 20, France, ‡Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504 Université de Strasbourg-CNRS, F-67034 Strasbourg, France, %these authors contributed equally.

KEYWORDS: Rhodopsin, Non-adiabatic Dynamics, Semi-classical Trajectories, Conical Intersection, Photoisomerization

**Abstract.** While the light-induced population dynamics of different photoresponsive proteins has been investigated spectroscopically, systematic computational studies have not yet been possible due to the phenomenally high cost of suitable high level quantum chemical methods and the need of propagating hundreds, if not thousands, non-adiabatic trajectories. Here we explore the possibility to study the photodynamics of rhodopsins by constructing and studying quantum mechanics/molecular mechanics (QM/MM) models featuring reduced retinal chromophores. In order to do so we use the sensory rhodopsin found in the cyanobacterium *Anabaena* PCC7120 (ASR) as a benchmark system. We find that the basic mechanical features associated with the excited state dynamics of ASR QM/MM models are reproduced using models incorporating a minimal (i.e. three double-bond) chromophore. Furthermore, we show that ensembles of non-adiabatic ASR trajectories computed using the same abridged models, replicate, at both the CASPT2 and CASSCF levels of theory, the trends in spectroscopy and lifetimes estimated using unabridged models and observed experimentally at room temperature. We conclude that a further expansion of these studies may lead to low-cost QM/MM rhodopsin models that may be used as effective tools in high-throughput *in silico* mutant screening.

## 1. Introduction

Retinylidene proteins (also called retinal proteins or, simply, rhodopsins) are a vast family of membrane proteins of great interest for photochemistry due to their ability to transduce light energy.¹⁻³ Although more than 300 such proteins are reported, Natural Selection has converged on very similar designs sharing common structural features.¹ For instance, all rhodopsins contain a chromophore (the protonated Schiff base of retinal) covalently bound into a pocket formed by seven interconnected transmembrane α-helices. Upon light irradiation the chromophore undergoes either a cis→trans or trans→cis double bond isomerization which triggers various biological functions such as vision in vertebrates and invertebrates, ion transportation, phototaxis and gene expression in archaea and eubacteria and ion-gating in algae. Other, possibly different, functions (e.g. in fungal rhodopsins) have not been clarified so far.¹

The photoisomerization of the rhodopsin chromophore occurs on a sub-picosecond timescale and usually shows a high quantum yield.³,⁴ An understanding of the molecular-level

factors determining such attractive properties is not only of fundamental importance for photobiology, but constitutes a step towards the construction of man-tailored (synthetic) biomimetic molecular devices of similar or augmented functionality and efficiency.⁵,⁶ Furthermore, the engineering of specific rhodopsins is becoming of strategic importance in optogenetics.⁷ Accordingly, the development of computational tools and/or strategies allowing a systematic and effective study of the photoisomerization of rhodopsins from different organisms and/or of entire sets of mutants of a specific rhodopsin is highly desirable.

The photoisomerization of rhodopsins is currently investigated constructing hybrid quantum mechanics/molecular mechanics (QM/MM) models of the protein starting from crystallographic data or suitable homology models. These QM/MM models are then used for subsequent photochemical reaction path or non-adiabatic trajectory computations. Trajectory studies have almost invariably relied on two different assumptions that impose severe limitations on the quality of achievable results. The first limitation is related to the high computational cost associated with the quantum mechanical
Probing the Photodynamics of Rhodopsins with Reduced
Retinal Chromophores

Madushanka Manathunga†,%, Xuchun Yang†,%, Hoi Ling Luk†, Samer Gozem†, Luis Manuel Fru-
tos#, Alessio Valentini¶,#, Nicolas Ferrè% and Massimo Olivucci†,¶,‡

†Department of Chemistry, Bowling Green State University, Bowling Green, OH 43403, USA, ¶Dipartimento di Bio-
tecnologie, Chimica e Farmacia, Università di Siena, via A. Moro 2, I-53100 Siena, Italy, #Departamento de Química
Física, Universidad de Alcalá, E-28871 Alcalá de Henares, Madrid, Spain, %Aix-Marseille Université, CNRS, Institut de
Chimie Radicalaire, 13397 Marseille Cedex 20, France, ‡Institut de Physique et Chimie des Matériaux de Strasbourg,
UMR 7504 Université de Strasbourg-CNRS, F-67034 Strasbourg, France, ^these authors contributed equally.

KEYWORDS: Rhodopsin, Non-adiabatic Dynamics, Semi-classical Trajectories, Conical Intersection, Photoisomeriza-
tion

**Abstract.** While the light-induced population dynamics of different photoresponsive proteins has been investigated spectroscopically, systematic computational studies have not yet been possible due to the phenomenally high cost of suitable high level quantum chemical methods and the need of propagating hundreds, if not thousands, non-adiabatic trajectories. Here we explore the possibility to study the photodynamics of rhodopsins by constructing and studying quantum mechanics/molecular mechanics (QM/MM) models featuring reduced retinal chromophores. In order to do so we use the sensory rhodopsin found in the cyanobacterium *Anabaena* PCC7120 (ASR) as a benchmark system. We find that the basic mechanical features associated with the excited state dynamics of ASR QM/MM models are reproduced using models incorporating a minimal (i.e. three double-bond) chromophore. Furthermore, we show that ensembles of non-adiabatic ASR trajectories computed using the same abridged models, replicate, at both the CASPT2 and CASSCF levels of theory, the trends in spectroscopy and lifetimes estimated using unabridged models and observed experimentally at room temperature. We conclude that a further expansion of these studies may lead to low-cost QM/MM rhodopsin models that may be used as effective tools in high-throughput *in silico* mutant screening.

**1. Introduction**

Retinylidene proteins (also called retinal proteins or, simply,
rhodopsins) are a vast family of membrane proteins of great
interest for photochemistry due to their ability to transduce
light energy. ¹⁻³ Although more than 300 such proteins are
reported, Natural Selection has converged on very similar
designs sharing common structural features. ¹ For instance, all
rhodopsins contain a chromophore (the protonated Schiff
base of retinal) covalently bound into a pocket formed by
seven interconnected transmembrane α-helices. Upon light
irradiation the chromophore undergoes either a cis→trans or
trans→cis double bond isomerization which triggers various
biological functions such as vision in vertebrates and inverte-
brates, ion transportation, phototaxis and gene expression in
archaea and eubacteria and ion-gating in algae. Other, possi-
bly different, functions (e.g. in fungal rhodopsins) have not
been clarified so far. ¹

The photoisomerization of the rhodopsin chromophore oc-
curs on a sub-picosecond timescale and usually shows a high
quantum yield.³,⁴ An understanding of the molecular-level

factors determining such attractive properties is not only of
fundamental importance for photobiology, but constitutes a
step towards the construction of man-tailored (synthetic)
biomimetic molecular devices of similar or augmented func-
tionality and efficiency.⁵,⁶ Furthermore, the engineering of
specific rhodopsins is becoming of strategic importance in
optogenetics.⁷ Accordingly, the development of computa-
tional tools and/or strategies allowing a systematic and effec-
tive study of the photoisomerization of rhodopsins from dif-
ferent organisms and/or of entire sets of mutants of a specific
rhodopsin is highly desirable.

The photoisomerization of rhodopsins is currently investigat-
ed constructing hybrid quantum mechanics/molecular me-
chanics (QM/MM) models of the protein starting from crys-
tallographic data or suitable homology models. These
QM/MM models are then used for subsequent photochemi-
cal reaction path or non-adiabatic trajectory computations.
Trajectory studies have almost invariably relied on two differ-
ent assumptions that impose severe limitations on the quality
of achievable results. The first limitation is related to the high
computational cost associated with the quantum mechanical
(QM) methods required for excited state and, therefore, non-adiabatic trajectory calculations. Such high cost limits the maximum number of trajectories that can be computed and, in principle, the statistical accuracy of the resulting dynamics and mechanism. Indeed, in many trajectory studies it is assumed that the computation of a single trajectory started from the protein ground state equilibrium structure at zero K (called Franck-Condon trajectory) provides information on the excited state dynamics of the entire population at room temperature. A second basic assumption adopted in QM/MM studies is related to the fact that the gradients driving the reaction path and trajectory calculations are obtained using QM method lacking a significant part of dynamic electron correlation such as the CASSCF method⁸. This is usually unavoidable, as the cost of gradients computed at more accurate QM levels such as CASPT2⁹ are unpractical, if not impossible, for a system of the size of the rhodopsin chromophore. Accordingly, rhodopsin studies have been carried out using a CASPT2//CASSCF protocol¹⁰ where the molecular structure evolution is computed at the CASSCF level and successively, the excitation energies and/or energy profiles along the isomerization coordinate are re-evaluated using single-point CASPT2 calculations.

The target of the present study is twofold. The first, substantially technical, target is the assessment of the validity of the assumptions revised above by comparing the results of CASSCF-based and CASPT2-based calculations in abridged but affordable QM/MM rhodopsin models. The second and more applicable target is to investigate the possibility to systematically employ the same affordable models for rapidly getting information on trends of kinetic parameters and photochemical reactivity in large rhodopsin sets at room temperature. In order to achieve these targets we build and employ QM/MM models incorporating structurally modified rhodopsin chromophores with three, rather than six, conjugating double bonds. Such models allow to compute (i) the hundreds of QM/MM trajectories necessary to achieve a description of the molecular population evolution at room temperature (this would be highly unpractical for unabridged models even when employing CASSCF gradients and supercomputer facilities) and (ii) the CASPT2 numerical gradient computations needed to run more accurate trajectories.

As a convenient benchmark system we employ the sensory rhodopsin from the cyanobacterium Anabaena PCC7120 (Anabaena sensory rhodopsin, ASR).⁶ In fact, while in most rhodopsins the function is initiated through a single isomerization (e.g. usually from all-trans to 13-cis in archaea and eu-bacteria rhodopsins and from 11-cis to all-trans in vertebrate and invertebrate rhodopsins), the ASR function is driven by the interconversion of two forms, ASR<sub>AT</sub> and ASR<sub>13C</sub>, featuring the all-trans,15-anti and 13-cis,15-syn chromophores respectively (see Scheme 1, top).¹¹ As already reported in the literature, it is possible to construct QM/MM models of both forms starting from the available crystallographic data. Accordingly, the ASR<sub>AT</sub> model can be employed to study the trans→cis isomerization and the ASR<sub>13C</sub> model to study the inverse reaction. Furthermore, these forms display significant differences in isomerization quantum yields and kinetics thus providing a

basis for evaluating the quality of the constructed computer models. Indeed, the quantum yields were reported to be 0.4 and 0.1 for ASR<sub>AT</sub> and ASR<sub>13C</sub> respectively<sup>12,13</sup>, while kinetics measurements indicate that ASR<sub>13C</sub> isomerizes in 0.15-0.25 ps whereas the opposite ASR<sub>AT</sub> isomerization takes 0.75-4.00 ps.<sup>12-14</sup> A recent study<sup>15</sup> has shown that Franck-Condon trajectories (from now on FC trajectories) computed using unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> QM/MM models produce kinetic data consistent with the corresponding observed quantities.

<figure><img src="image_2.png" /><figcaption>**Scheme 1.** Reduced retinal models of ASR<sub>AT</sub> and ASR<sub>13C</sub>. Constructed by removing a part of the chromophore shown in red.</figcaption></figure>

As anticipated above, the present study is based on the construction of ASR<sub>AT</sub> and ASR<sub>13C</sub> models incorporating abridged retinal chromophores featuring three double bonds (from now on called PSB3 chromophores). These modified models, ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively, are obtained by shortening the native chromophores as shown in Fig. 1. In the gas-phase, the PSB3 chromophore is known to reproduce several features of the full chromophore¹⁶ and for this reason has been extensively used as a benchmark for investigating the performance of several QM methods.<sup>17-20</sup> For instance, these studies have revealed substantial differences between CASSCF and CASPT2 ground state equilibrium structures and photochemical reaction paths. Below, we will assess the magnitude of such differences in ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models where the PSB3 chromophore is mechanically and electrostatically embedded in the protein "molecular" environment and not studied in isolated conditions.

In Section 2 we report a detailed description of the construction of the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models and of the sampling method used to simulate their room temperature populations. Section 3 is divided in three result subsections. In subsection 3.1 we use the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models for comparing the ground state equilibrium structures and the FC trajectories computed at the CASSCF and CASPT2 levels with the reported CASSCF trajectories computed for the full system. In Section 3.2 we compare the result of small sets of room temperature CASSCF and CASPT2 non-adiabatic (Tully surface-hop) trajectories for each abridged model and show that both QM levels yield, in contrast with the FC trajectory results, similar excited state lifetimes. Finally, in subsection 3.3 we study the room temperature photoisomerization kinetics of a 400 non-adiabatic trajectory ensembles for both the ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M model at the CASSCF level. We find that, consistently with the FC trajectory results, that the ASR<sub>13C</sub>-M population decay is significantly faster than the
(QM) methods required for excited state and, therefore, non-adiabatic trajectory calculations. Such high cost limits the maximum number of trajectories that can be computed and, in principle, the statistical accuracy of the resulting dynamics and mechanism. Indeed, in many trajectory studies it is assumed that the computation of a single trajectory started from the protein ground state equilibrium structure at zero K (called Franck-Condon trajectory) provides information on the excited state dynamics of the entire population at room temperature. A second basic assumption adopted in QM/MM studies is related to the fact that the gradients driving the reaction path and trajectory calculations are obtained using QM method lacking a significant part of dynamic electron correlation such as the CASSCF method⁸. This is usually unavoidable, as the cost of gradients computed at more accurate QM levels such as CASPT2⁹ are unpractical, if not impossible, for a system of the size of the rhodopsin chromophore. Accordingly, rhodopsin studies have been carried out using a CASPT2//CASSCF protocol¹⁰ where the molecular structure evolution is computed at the CASSCF level and successively, the excitation energies and/or energy profiles along the isomerization coordinate are re-evaluated using single-point CASPT2 calculations.

The target of the present study is twofold. The first, substantially technical, target is the assessment of the validity of the assumptions revised above by comparing the results of CASSCF-based and CASPT2-based calculations in abridged but affordable QM/MM rhodopsin models. The second and more applicative target is to investigate the possibility to systematically employ the same affordable models for rapidly getting information on trends of kinetic parameters and photochemical reactivity in large rhodopsin sets at room temperature. In order to achieve these targets we build and employ QM/MM models incorporating structurally modified rhodopsin chromophores with three, rather than six, conjugating double bonds. Such models allow to compute (i) the hundreds of QM/MM trajectories necessary to achieve a description of the molecular population evolution at room temperature (this would be highly unpractical for unabridged models even when employing CASSCF gradients and supercomputer facilities) and (ii) the CASPT2 numerical gradient computations needed to run more accurate trajectories.

As a convenient benchmark system we employ the sensory rhodopsin from the cyanobacterium Anabaena PCC7120 (Anabaena sensory rhodopsin, ASR).⁶ In fact, while in most rhodopsins the function is initiated through a single isomerization (e.g. usually from all-trans to 13-cis in archaea and eu-bacteria rhodopsins and from 11-cis to all-trans in vertebrate and invertebrate rhodopsins), the ASR function is driven by the interconversion of two forms, ASR<sub>AT</sub> and ASR<sub>13C</sub>, featuring the all-trans, 15-anti and 13-cis, 15-syn chromophores respectively (see Scheme 1, top).¹¹ As already reported in the literature, it is possible to construct QM/MM models of both forms starting from the available crystallographic data. Accordingly, the ASR<sub>AT</sub> model can be employed to study the trans→cis isomerization and the ASR<sub>13C</sub> model to study the inverse reaction. Furthermore, these forms display significant differences in isomerization quantum yields and kinetics thus providing a

basis for evaluating the quality of the constructed computer models. Indeed, the quantum yields were reported to be 0.4 and 0.1 for ASR<sub>AT</sub> and ASR<sub>13C</sub> respectively<sup>12,13</sup>, while kinetics measurements indicate that ASR<sub>13C</sub> isomerizes in 0.15-0.25 ps whereas the opposite ASR<sub>AT</sub> isomerization takes 0.75-4.00 ps.<sup>12-14</sup> A recent study<sup>15</sup> has shown that Franck-Condon trajectories (from now on FC trajectories) computed using unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> QM/MM models produce kinetic data consistent with the corresponding observed quantities.

<figure><img src="image_2.png" /><figcaption>**Scheme 1.** Reduced retinal models of ASR<sub>AT</sub> and ASR<sub>13C</sub>. Constructed by removing a part of the chromophore shown in red.</figcaption></figure>

As anticipated above, the present study is based on the construction of ASR<sub>AT</sub> and ASR<sub>13C</sub> models incorporating abridged retinal chromophores featuring three double bonds (from now on called PSB3 chromophores). These modified models, ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively, are obtained by shortening the native chromophores as shown in Fig. 1. In the gas-phase, the PSB3 chromophore is known to reproduce several features of the full chromophore¹⁶ and for this reason has been extensively used as a benchmark for investigating the performance of several QM methods.<sup>17-20</sup> For instance, these studies have revealed substantial differences between CASSCF and CASPT2 ground state equilibrium structures and photochemical reaction paths. Below, we will assess the magnitude of such differences in ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models where the PSB3 chromophore is mechanically and electrostatically embedded in the protein "molecular" environment and not studied in isolated conditions.

In Section 2 we report a detailed description of the construction of the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models and of the sampling method used to simulate their room temperature populations. Section 3 is divided in three result subsections. In subsection 3.1 we use the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models for comparing the ground state equilibrium structures and the FC trajectories computed at the CASSCF and CASPT2 levels with the reported CASSCF trajectories computed for the full system. In Section 3.2 we compare the result of small sets of room temperature CASSCF and CASPT2 non-adiabatic (Tully surface-hop) trajectories for each abridged model and show that both QM levels yield, in contrast with the FC trajectory results, similar excited state lifetimes. Finally, in subsection 3.3 we study the room temperature photoisomerization kinetics of a 400 non-adiabatic trajectory ensembles for both the ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M model at the CASSCF level. We find that, consistently with the FC trajectory results, that the ASR<sub>13C</sub>-M population decay is significantly faster than the
decay of the ASR<sub>AT</sub>-M population thus reproducing the trend observed for the corresponding experimental systems.

## 2. Methods

The ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models were constructed by discarding the C1-C10 portion of the retinal chromophore from the corresponding unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> models<sup>6,15</sup> respectively (see Scheme 1 for a schematic illustration of the process). Accordingly, the equilibrium structures of the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models were re-constructed following the protocol reported in Schapiro et al.<sup>15</sup> Their chromophores were then trimmed to obtain the corresponding cavity embedded PSB3 chromophores (see Fig. 1). As usual, at the frontier between the QM subsystem (i.e. the PSB3 chromophore) and MM subsystem (the protein), a hydrogen link atom was employed to saturate the bond between the last QM atom and the first MM atom in the Lys296 residue.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The QM/MM model of the reduced PSB3 retinal chromophore ASR<sub>AT</sub>-M.</figcaption></figure>

The original chromophore cavities, defined by the MM subsystem side-chains with at least one atom not more than 4 Å far from any QM atom, were redefined. The resulting models were then relaxed (only the cavity residues and the QM subsystem are free to move) to obtain the equilibrium geometries of the final ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models. As shown in Fig. 2 for ASR<sub>AT</sub>-M, the cavities hosting the reduced PSB3 chromophores conserve, substantially, the same cavity shape of the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models as the above procedure avoids the collapse of the residues in the vicinity of the removed moieties upon relaxation.

CASSCF/6-31G*/AMBER and CASPT2/6-31G*/AMBER (from now on called CASSCF and CASPT2 respectively) energy and gradient calculations were performed using the MOLCAS software package<sup>21</sup> with an active space comprising 6 π-electrons in 6 π-orbitals corresponding to the full π-systems of the chromophore. The TINKER software package<sup>22</sup> was coupled to MOLCAS to perform the QM/MM calculations and employs the AMBER-FF94 parameters. The description of the electrostatic embedding of the QM subsystem into the MM subsystem<sup>23</sup> was treated using the electrostatic potential fitted methodology.<sup>24</sup>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Structure of the full and modified 4 Å cavities of ASR<sub>AT</sub> and ASR<sub>AT</sub>-M. (A) cavity of the full ASR<sub>AT</sub> model. The QM chromophore is in green, the chromophore bound lysine (Lys 296) is in blue and rest of the MM cavity residues are in red. (B) cavity defined around the reduced retinal chromophore of the ASR<sub>AT</sub>-M model. The MM residues originally belonging to the ASR<sub>AT</sub> cavity and located outside the ASR<sub>AT</sub>-M cavity are shown in black and are kept frozen during the computations.</figcaption></figure>

FC trajectories were computed using the deterministic surface-hop method<sup>25-28</sup> also employed in Schapiro et al.<sup>15</sup>. The room temperature trajectories were instead computed using the stochastic Tully surface-hop method<sup>29,30</sup>. The required initial conditions were computed at the CASSCF level for the CASSCF trajectories and at the MP2 level for the CASPT2 trajectories using a protocol yielding Boltzmann-like distributions. More specifically, 400 initial conditions were generated using the following protocol:

* Room temperature ground state molecular dynamics (40 ns) was performed at the MM level for both the unabridged ASR<sub>13C</sub> and ASR<sub>AT</sub> models using a MM parameterized retinal chromophore and keeping fixed the protein backbone and side-chains outside the 4 Å cavity region (see above).

* Geometries and velocities corresponding to 400 snapshots were extracted from the above simulations every 100 ps and ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M models were created for each snapshot using a specifically written FORTRAN script which substitutes the chromophores with the corresponding abridged PSB3 chromophore. The new 4 Å cavities are then redefined as explained above.

* Each relaxed ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M model was used to propagate a 200 fs ground state trajectory with a step size of 1 fs at the two-root state averaged CASSCF level. The final geometry and velocity vectors were taken as initial conditions for the subsequent CASSCF non-adiabatic trajectory computation. Similarly, a ground state QM/MM trajectory was propagated for 200 fs at the MP2 level for generating the initial conditions for the corresponding CASPT2 trajectories. In all cases the geometry and velocity vectors of the last snapshot of each ground state trajectory were treated as the initial conditions.

## 3. Results and Discussion
decay of the ASR<sub>AT</sub>-M population thus reproducing the trend observed for the corresponding experimental systems.

## 2. Methods

The ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models were constructed by discarding the C1-C10 portion of the retinal chromophore from the corresponding unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> models<sup>6,15</sup> respectively (see Scheme 1 for a schematic illustration of the process). Accordingly, the equilibrium structures of the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models were re-constructed following the protocol reported in Schapiro et al.<sup>15</sup> Their chromophores were then trimmed to obtain the corresponding cavity embedded PSB3 chromophores (see Fig. 1). As usual, at the frontier between the QM subsystem (i.e. the PSB3 chromophore) and MM subsystem (the protein), a hydrogen link atom was employed to saturate the bond between the last QM atom and the first MM atom in the Lys296 residue.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The QM/MM model of the reduced PSB3 retinal chromophore ASR<sub>AT</sub>-M.</figcaption></figure>

The original chromophore cavities, defined by the MM subsystem side-chains with at least one atom not more than 4 Å far from any QM atom, were redefined. The resulting models were then relaxed (only the cavity residues and the QM subsystem are free to move) to obtain the equilibrium geometries of the final ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models. As shown in Fig. 2 for ASR<sub>AT</sub>-M, the cavities hosting the reduced PSB3 chromophores conserve, substantially, the same cavity shape of the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models as the above procedure avoids the collapse of the residues in the vicinity of the removed moieties upon relaxation.

CASSCF/6-31G*/*AMBER and CASPT2/6-31G*/*AMBER (from now on called CASSCF and CASPT2 respectively) energy and gradient calculations were performed using the MOLCAS software package<sup>21</sup> with an active space comprising 6 π-electrons in 6 π-orbitals corresponding to the full π-systems of the chromophore. The TINKER software package<sup>22</sup> was coupled to MOLCAS to perform the QM/MM calculations and employs the AMBER-FF94 parameters. The description of the electrostatic embedding of the QM subsystem into the MM subsystem<sup>23</sup> was treated using the electrostatic potential fitted methodology.<sup>24</sup>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Structure of the full and modified 4 Å cavities of ASR<sub>AT</sub> and ASR<sub>AT</sub>-M. (A) cavity of the full ASR<sub>AT</sub> model. The QM chromophore is in green, the chromophore bound lysine (Lys 296) is in blue and rest of the MM cavity residues are in red. (B) cavity defined around the reduced retinal chromophore of the ASR<sub>AT</sub>-M model. The MM residues originally belonging to the ASR<sub>AT</sub> cavity and located outside the ASR<sub>AT</sub>-M cavity are shown in black and are kept frozen during the computations.</figcaption></figure>

FC trajectories were computed using the deterministic surface-hop method<sup>25-28</sup> also employed in Schapiro et al.<sup>15</sup>. The room temperature trajectories were instead computed using the stochastic Tully surface-hop method<sup>29,30</sup>. The required initial conditions were computed at the CASSCF level for the CASSCF trajectories and at the MP2 level for the CASPT2 trajectories using a protocol yielding Boltzmann-like distributions. More specifically, 400 initial conditions were generated using the following protocol:

* Room temperature ground state molecular dynamics (40 ns) was performed at the MM level for both the unabridged ASR<sub>13C</sub> and ASR<sub>AT</sub> models using a MM parameterized retinal chromophore and keeping fixed the protein backbone and side-chains outside the 4 Å cavity region (see above).

* Geometries and velocities corresponding to 400 snapshots were extracted from the above simulations every 100 ps and ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M models were created for each snapshot using a specifically written FORTRAN script which substitutes the chromophores with the corresponding abridged PSB3 chromophore. The new 4 Å cavities are then redefined as explained above.

* Each relaxed ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M model was used to propagate a 200 fs ground state trajectory with a step size of 1 fs at the two-root state averaged CASSCF level. The final geometry and velocity vectors were taken as initial conditions for the subsequent CASSCF non-adiabatic trajectory computation. Similarly, a ground state QM/MM trajectory was propagated for 200 fs at the MP2 level for generating the initial conditions for the corresponding CASPT2 trajectories. In all cases the geometry and velocity vectors of the last snapshot of each ground state trajectory were treated as the initial conditions.

## 3. Results and Discussion
<figure><img src="image_5.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of ground state optimized geometries of ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) in the protein environment and gas phase computed at both the CASPT2 and CASSCF levels. The bond length is in Angstrom and the -C12-C13=C14-C15- dihedral (in parenthesis) is in degree. See Table of Cartesian coordinates in the Supporting Info.</figcaption></figure>

**3.1 Comparison of CASSCF and CASPT2 geometries and FC trajectories.** Ground state (S<sub>0</sub>) geometry optimizations were performed for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M QM/MM models as well as for the corresponding isolated (i.e. gas-

phase) PSB3 chromophores. The optimizations were carried out using both the CASPT2 and CASSCF levels thus yielding two equilibrium structures for each molecular system investigated. The main geometrical parameters (bond lengths and dihedral angles) are reported in Fig. 3 for both the protein and the isolated chromophores. As previously documented for gas-phase systems<sup>17,18,20</sup>, the results indicate that the two levels of theory yield different bond length alternation (BLA). The BLA value is computed as the difference between the average single and average double bond lengths of the conjugated framework and reflects the higher degree of delocalization (i.e. a BLA value close to zero reflects similar single and double bond lengths) for the CASPT2 level. This effect, which appears insensitive to the protein environment, is dominated by the =C-C=N- moiety as documented by the change in =C-C= and -C=N- bond lengths which display 0.02-0.03 Å variations. On the other hand, a comparison of the isolated and protein-embedded chromophores indicates that the protein environment leads to ca. 10-15 degrees twisted central double bond in both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models and at both the CASSCF and CASPT2 levels.

The ultrafast singlet excited state (S<sub>1</sub>) dynamics of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M has been initially probed by computing the corresponding FC trajectories at both the CASPT2 and CASSCF levels. In both cases we propagated the trajectories for 200 fs. The simulated S<sub>1</sub> lifetime (τ<sub>FC</sub>) is defined as the time required to intercept the S<sub>1</sub>/S<sub>0</sub> conical intersections (CI) located halfway along the isomerization coordinate (typically when the C13=C14 double bond is ca. 90 degrees twisted). τ<sub>FC</sub> is often assumed to represent the time required by the center of the S<sub>1</sub> population for reaching the CI.<sup>10,31</sup> In Table 1 the results are compared to the observed S<sub>1</sub> lifetimes for ASR<sup>14</sup> and with the previously reported FC trajectories<sup>15</sup> for the unabridged ASR models computed at the CASPT2//CASSCF level of theory.

<table><caption>**Table 1.** Comparison of observed<sup>12,14</sup> and computed absorption maxima and excited state lifetimes (using FC trajectories) values for unabridged ASR<sup>15</sup> and modified ASR models.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">ASR<br>Observed</th><th colspan="2">ASR<br>CASPT2//CASSCF</th><th colspan="2">ASR-M<br>CASPT2</th><th colspan="2">ASR-M<br>CASSCF</th></tr><tr><th>τ<br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th></tr></thead><tbody><tr><td>AT</td><td>750/4000</td><td>549</td><td>&gt;400</td><td>533</td><td>&gt;200</td><td>351</td><td>61</td><td>336</td></tr><tr><td>13C</td><td>150</td><td>537</td><td>150</td><td>527</td><td>63</td><td>345</td><td>53</td><td>329</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of ground state optimized geometries of ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) in the protein environment and gas phase computed at both the CASPT2 and CASSCF levels. The bond length is in Angstrom and the -C12-C13=C14-C15- dihedral (in parenthesis) is in degree. See Table of Cartesian coordinates in the Supporting Info.</figcaption></figure>

**3.1 Comparison of CASSCF and CASPT2 geometries and FC trajectories.** Ground state (S<sub>0</sub>) geometry optimizations were performed for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M QM/MM models as well as for the corresponding isolated (i.e. gas-

phase) PSB3 chromophores. The optimizations were carried out using both the CASPT2 and CASSCF levels thus yielding two equilibrium structures for each molecular system investigated. The main geometrical parameters (bond lengths and dihedral angles) are reported in Fig. 3 for both the protein and the isolated chromophores. As previously documented for gas-phase systems<sup>17,18,20</sup>, the results indicate that the two levels of theory yield different bond length alternation (BLA). The BLA value is computed as the difference between the average single and average double bond lengths of the conjugated framework and reflects the higher degree of delocalization (i.e. a BLA value close to zero reflects similar single and double bond lengths) for the CASPT2 level. This effect, which appears insensitive to the protein environment, is dominated by the =C-C=N- moiety as documented by the change in =C-C= and -C=N- bond lengths which display 0.02-0.03 Å variations. On the other hand, a comparison of the isolated and protein-embedded chromophores indicates that the protein environment leads to ca. 10-15 degrees twisted central double bond in both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models and at both the CASSCF and CASPT2 levels.

The ultrafast singlet excited state (S<sub>1</sub>) dynamics of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M has been initially probed by computing the corresponding FC trajectories at both the CASPT2 and CASSCF levels. In both cases we propagated the trajectories for 200 fs. The simulated S<sub>1</sub> lifetime (τ<sub>FC</sub>) is defined as the time required to intercept the S<sub>1</sub>/S<sub>0</sub> conical intersections (CI) located halfway along the isomerization coordinate (typically when the C13=C14 double bond is ca. 90 degrees twisted). τ<sub>FC</sub> is often assumed to represent the time required by the center of the S<sub>1</sub> population for reaching the CI.<sup>10,31</sup> In Table 1 the results are compared to the observed S<sub>1</sub> lifetimes for ASR<sup>14</sup> and with the previously reported FC trajectories<sup>15</sup> for the unabridged ASR models computed at the CASPT2//CASSCF level of theory.

<table><caption>**Table 1.** Comparison of observed<sup>12,14</sup> and computed absorption maxima and excited state lifetimes (using FC trajectories) values for unabridged ASR<sup>15</sup> and modified ASR models.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">ASR<br>Observed</th><th colspan="2">ASR<br>CASPT2//CASSCF</th><th colspan="2">ASR-M<br>CASPT2</th><th colspan="2">ASR-M<br>CASSCF</th></tr><tr><th>τ<br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th><th>τ<sub>FC</sub><br>(fs)</th><th>λ<sub>max</sub><br>(nm)</th></tr></thead><tbody><tr><td>AT</td><td>750/4000</td><td>549</td><td>&gt;400</td><td>533</td><td>&gt;200</td><td>351</td><td>61</td><td>336</td></tr><tr><td>13C</td><td>150</td><td>537</td><td>150</td><td>527</td><td>63</td><td>345</td><td>53</td><td>329</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** QM/MM trajectories of modified models (ASR<sub>13C</sub>-M, ASR<sub>AT</sub>-M) computed at the CASPT2 and CASSCF level as well as full models (ASR<sub>13C</sub>, ASR<sub>AT</sub>) computed at CASSCF and scaled-CASSCF level<sup>10</sup> of theory. The trajectories have been released from the corresponding S<sub>0</sub> equilibrium structure. The circles represent points of surface hop located in the vicinity of a S<sub>1</sub>/S<sub>0</sub> CI point.</figcaption></figure>

The corresponding S₀ and S₁ potential energy profiles are reported in Fig. 4. ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models yield the same trend (but very different magnitudes) in τ<sub>FC</sub> values found for the unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> models as well as for the observed data. In fact, at both levels of theory, τ<sub>FC</sub> for ASR<sub>AT</sub>-M is larger than that of the ASR<sub>13C</sub>-M indicating a longer S₁ lifetime and a slower product formation time.

Notice that at the CASPT2 level, ASR<sub>AT</sub>-M does not decay to S<sub>0</sub> within 200 fs. This establishes a qualitative difference between the CASPT2 and CASSCF trajectories of ASR<sub>AT</sub>-M that is attributed to the enhanced flatness of the S<sub>1</sub> CASPT2 potential energy surface with respect to the CASSCF one.<sup>20</sup> While a flatter isomerization path or barrier would prevent ASR<sub>AT</sub>-M to decay quickly, ASR<sub>13C</sub>-M still decays rapidly probably due to the enhanced 13-cis conformational strain (see Fig. 5) imposing a larger initial S<sub>1</sub> acceleration. More specifically, at the CASPT2 level ASR<sub>AT</sub>-M displays a smaller potential energy decrease (i.e., smaller kinetic energy increase) upon relaxation out of the FC region with respect to ASR<sub>13C</sub>-M. This behavior is related to its initial BLA value closer to the S<sub>1</sub> average value (see Fig. 5A) and preventing a large single-bond/double-bond S<sub>1</sub> relaxation (compare the BLA values and energy relaxation of Fig. 5A and Fig. 5B). Since the only source of kinetic energy available to a FC trajectory originates from the potential energy relaxation, the lower kinetic energy content of the CASPT2 trajectory with respect to the CASSCF trajectory becomes a factor controlling the reactivity of the double bond. The increase in flatness and even the presence of a small barrier along the CASPT2 potential energy surfaces has been previously documented for a gas-phase PSB3 model.<sup>19</sup>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Geometrical changes along the computed FC trajectories. The bond length alteration (BLA) of modified models at CASPT2 and CASSCF levels as well as full models at CASSCF level. The values in red dashed and red solid frames are initial BLA values of modified models at the CASSCF and CASPT2 levels.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** QM/MM trajectories of modified models (ASR<sub>13C</sub>-M, ASR<sub>AT</sub>-M) computed at the CASPT2 and CASSCF level as well as full models (ASR<sub>13C</sub>, ASR<sub>AT</sub>) computed at CASSCF and scaled-CASSCF level<sup>10</sup> of theory. The trajectories have been released from the corresponding S<sub>0</sub> equilibrium structure. The circles represent points of surface hop located in the vicinity of a S<sub>1</sub>/S<sub>0</sub> CI point.</figcaption></figure>

The corresponding S₀ and S₁ potential energy profiles are reported in Fig. 4. ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models yield the same trend (but very different magnitudes) in τ<sub>FC</sub> values found for the unabridged ASR<sub>AT</sub> and ASR<sub>13C</sub> models as well as for the observed data. In fact, at both levels of theory, τ<sub>FC</sub> for ASR<sub>AT</sub>-M is larger than that of the ASR<sub>13C</sub>-M indicating a longer S₁ lifetime and a slower product formation time.

Notice that at the CASPT2 level, ASR<sub>AT</sub>-M does not decay to S<sub>0</sub> within 200 fs. This establishes a qualitative difference between the CASPT2 and CASSCF trajectories of ASR<sub>AT</sub>-M that is attributed to the enhanced flatness of the S<sub>1</sub> CASPT2 potential energy surface with respect to the CASSCF one.<sup>20</sup> While a flatter isomerization path or barrier would prevent ASR<sub>AT</sub>-M to decay quickly, ASR<sub>13C</sub>-M still decays rapidly probably due to the enhanced 13-cis conformational strain (see Fig. 5) imposing a larger initial S<sub>1</sub> acceleration. More specifically, at the CASPT2 level ASR<sub>AT</sub>-M displays a smaller potential energy decrease (i.e., smaller kinetic energy increase) upon relaxation out of the FC region with respect to ASR<sub>13C</sub>-M. This behavior is related to its initial BLA value closer to the S<sub>1</sub> average value (see Fig. 5A) and preventing a large single-bond/double-bond S<sub>1</sub> relaxation (compare the BLA values and energy relaxation of Fig. 5A and Fig. 5B). Since the only source of kinetic energy available to a FC trajectory originates from the potential energy relaxation, the lower kinetic energy content of the CASPT2 trajectory with respect to the CASSCF trajectory becomes a factor controlling the reactivity of the double bond. The increase in flatness and even the presence of a small barrier along the CASPT2 potential energy surfaces has been previously documented for a gas-phase PSB3 model.<sup>19</sup>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Geometrical changes along the computed FC trajectories. The bond length alteration (BLA) of modified models at CASPT2 and CASSCF levels as well as full models at CASSCF level. The values in red dashed and red solid frames are initial BLA values of modified models at the CASSCF and CASPT2 levels.</figcaption></figure>
A

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

B

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 6.** The progression along the twisting angles describing the isomerization coordinate for ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) models at the CASSCF and CASPT2 levels.</figcaption></figure>

One also needs to explain why the reported FC trajectory of the unabridged ASR<sub>AT</sub> model does not decay within 200 fs even when using the steeper CASSCF gradients. This behavior is attributed to the increase in length of the conjugated framework (from three to six conjugated double bonds) in ASR<sub>AT</sub> making the S<sub>1</sub> potential energy surface flatter as documented in the literature for a series of retinal chromophore models of different lengths.<sup>32</sup> The corresponding CASSCF S<sub>1</sub> potential for ASR<sub>AT</sub>-M (with the three double bond chromophore) is not flat enough to prevent acceleration along the reaction coordinate even if, clearly, the tendency to a slower reactivity remains (i.e. the ASR<sub>AT</sub>-M S<sub>1</sub> lifetime is ca. 10 fs longer than that of ASR<sub>13C</sub>-M).

An additional test reinforcing the similarity between the abridged models and the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models is related to the S<sub>1</sub> reaction coordinate. The reaction coordinate of both ASR<sub>AT</sub> forms corresponds to a bicycle-pedal motion described by the clockwise and counterclockwise twisting of the C13=C14 and C15=N bonds while the C11=C12 remains substantially unreactive. As shown in Fig. 6 this is indeed what is found when the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models are used at the CASSCF and, for ASR<sub>13C</sub>-M only, CASPT2 levels. It is concluded that, mechanistically, the PSB3 in the abridged models follow the same, relatively complex, reaction coordinate of full chromophore models<sup>6</sup> (with different speeds).

**3.2 Probing the room temperature dynamics with small trajectory sets.** Above we have made the hypothesis that a smaller kinetic energy together with a flatter S<sub>1</sub> potential energy surface could be the cause of the missed/slower ASR<sub>AT</sub>-M reactivity seen in the CASPT2 FC trajectory. As mentioned above, since FC trajectories are 0 K trajectories, it is important to test the above conclusion investigating the CASSCF and CASPT2 dynamics of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M at room temperature (i.e. when one has initial non-zero velocity vectors). In order to do so we focused on a small set of room temperature non-adiabatic trajectories defined by selecting the corresponding initial conditions out of the Boltzmann-like distribution generated using the protocol described in Section 2.

Before computing the trajectory, the quality of the computed Boltzmann-like distributions was evaluated by plotting the potential energy and dihedral angle (C12-C13=C14-C15) values. These are reported in the bar-diagrams of Fig. 7 for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M and show the expected Gaussian-like shape. As reported in Section 2, for both ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M, the distribution employed for the CASPT2 calculations is
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 6.** The progression along the twisting angles describing the isomerization coordinate for ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) models at the CASSCF and CASPT2 levels.</figcaption></figure>

One also needs to explain why the reported FC trajectory of the unabridged ASR<sub>AT</sub> model does not decay within 200 fs even when using the steeper CASSCF gradients. This behavior is attributed to the increase in length of the conjugated framework (from three to six conjugated double bonds) in ASR<sub>AT</sub> making the S<sub>1</sub> potential energy surface flatter as documented in the literature for a series of retinal chromophore models of different lengths.<sup>32</sup> The corresponding CASSCF S<sub>1</sub> potential for ASR<sub>AT</sub>-M (with the three double bond chromophore) is not flat enough to prevent acceleration along the reaction coordinate even if, clearly, the tendency to a slower reactivity remains (i.e. the ASR<sub>AT</sub>-M S<sub>1</sub> lifetime is ca. 10 fs longer than that of ASR<sub>13C</sub>-M).

An additional test reinforcing the similarity between the abridged models and the original ASR<sub>AT</sub> and ASR<sub>13C</sub> models is related to the S<sub>1</sub> reaction coordinate. The reaction coordinate of both ASR<sub>AT</sub> forms corresponds to a bicycle-pedal motion described by the clockwise and counterclockwise twisting of the C13=C14 and C15=N bonds while the C11=C12 remains substantially unreactive. As shown in Fig. 6 this is indeed what is found when the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models are used at the CASSCF and, for ASR<sub>13C</sub>-M only, CASPT2 levels. It is concluded that, mechanistically, the PSB3 in the abridged models follow the same, relatively complex, reaction coordinate of full chromophore models<sup>6</sup> (with different speeds).

**3.2 Probing the room temperature dynamics with small trajectory sets.** Above we have made the hypothesis that a smaller kinetic energy together with a flatter S<sub>1</sub> potential energy surface could be the cause of the missed/slower ASR<sub>AT</sub>-M reactivity seen in the CASPT2 FC trajectory. As mentioned above, since FC trajectories are 0 K trajectories, it is important to test the above conclusion investigating the CASSCF and CASPT2 dynamics of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M at room temperature (i.e. when one has initial non-zero velocity vectors). In order to do so we focused on a small set of room temperature non-adiabatic trajectories defined by selecting the corresponding initial conditions out of the Boltzmann-like distribution generated using the protocol described in Section 2.

Before computing the trajectory, the quality of the computed Boltzmann-like distributions was evaluated by plotting the potential energy and dihedral angle (C12-C13=C14-C15) values. These are reported in the bar-diagrams of Fig. 7 for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M and show the expected Gaussian-like shape. As reported in Section 2, for both ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M, the distribution employed for the CASPT2 calculations is
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 7.** The total energy distribution of 400 initial conditions of ASR<sub>AT</sub>-M (A), ASR<sub>13C</sub>-M (B) at the CASSCF level, taking the mean energy as the reference. The distribution of the central dihedral (-C12-C13=C14-C15-) of ASR<sub>AT</sub>-M (C), ASR<sub>13C</sub>-M (D) at the same level of theory. The brown dots represent the initial conditions selected for ten-trajectory study. Note that each dot represents a single initial condition. Two initial conditions out of ten for the ASR<sub>13C</sub>-M case and one for the ASR<sub>AT</sub>-M case had very low energy. These are not shown on histograms.</figcaption></figure>

obtained at MP2 level (the CASPT2 and MP2 S₀ equilibrium geometries are very similar as shown in Section 1 of the supporting information) while the initial conditions for the CASSCF calculations were obtained at CASSCF level. The standard deviation with respect to the mean energy values for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M are 10 and 15 kcal/mol respectively. The Gaussian-like shaped distribution of central dihedral angles show that the ensembles are centered on -170 and 20 degrees for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively. These mean values of the dihedral angles indicate that both models are pre-twisted in the same clockwise direction due to the asymmetric protein environment (i.e. dihedrals of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M are not exactly -180 and 0 degrees) and consistently with the corresponding optimized structures of Fig. 4 and isomerization direction of Fig. 6.

For each model an ensemble of ten trajectories (indicated by the circles in Fig. 7) was computed at both the CASPT2 and CASSCF levels. The evolution of the values of skeletal dihedral angle representing the C13=C14 reacting double bond, are given at the top of Fig. 8. The blue dots show the points where trajectories decay from S₁ to S₀ at CASPT2 level while the green dots show the corresponding CASSCF decay points. In contrast with the results of the FC trajectories (see above), both the CASPT2 and CASSCF trajectories of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M decay to S₀ (i.e. we do not observe qualitative differences between the two QM levels). Furthermore, the diagram in Fig. 8 indicates that the CASSCF and CASPT2 sets decay on similar timescales for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models. On the other hand, consistently with the FC trajectory description in the previous subsection, the room temperature trajectories confirm that the ASR<sub>AT</sub>-M decays slower than ASR<sub>13C</sub>-M. In fact, the corresponding average hop times were 69 fs and 54 fs at the CASSCF level and 102 fs and 52 fs at the CASPT2 level for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively.

It is important to establish if the bicycle pedal reaction coordinate followed by the 0 K FC trajectories is affected by the room temperature condition. Accordingly, in Fig. 8 we also display the twisting deformation of the C13=C14 and C15=N double bonds along the trajectory sets. It is shown that, in both ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M, these dihedrals follow the same clockwise and counterclockwise motion observed in the corresponding FC trajectories irrespective of the level of QM theory employed. In particular, it is shown that the ASR<sub>AT</sub>-M reactive CASPT2 trajectories, which are not described by the corresponding FC trajectory, still describe a bicycle pedal motion. These results indicate that the specific construction of the reduced chromophore cavity (having part of the residues frozen as in the original chromophore cavity) has not perturbed the selectivity and isomerization mechanism of the photoisomerization.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 7.** The total energy distribution of 400 initial conditions of ASR<sub>AT</sub>-M (A), ASR<sub>13C</sub>-M (B) at the CASSCF level, taking the mean energy as the reference. The distribution of the central dihedral (-C12-C13=C14-C15-) of ASR<sub>AT</sub>-M (C), ASR<sub>13C</sub>-M (D) at the same level of theory. The brown dots represent the initial conditions selected for ten-trajectory study. Note that each dot represents a single initial condition. Two initial conditions out of ten for the ASR<sub>13C</sub>-M case and one for the ASR<sub>AT</sub>-M case had very low energy. These are not shown on histograms.</figcaption></figure>

obtained at MP2 level (the CASPT2 and MP2 S₀ equilibrium geometries are very similar as shown in Section 1 of the supporting information) while the initial conditions for the CASSCF calculations were obtained at CASSCF level. The standard deviation with respect to the mean energy values for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M are 10 and 15 kcal/mol respectively. The Gaussian-like shaped distribution of central dihedral angles show that the ensembles are centered on -170 and 20 degrees for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively. These mean values of the dihedral angles indicate that both models are pre-twisted in the same clockwise direction due to the asymmetric protein environment (i.e. dihedrals of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M are not exactly -180 and 0 degrees) and consistently with the corresponding optimized structures of Fig. 4 and isomerization direction of Fig. 6.

For each model an ensemble of ten trajectories (indicated by the circles in Fig. 7) was computed at both the CASPT2 and CASSCF levels. The evolution of the values of skeletal dihedral angle representing the C13=C14 reacting double bond, are given at the top of Fig. 8. The blue dots show the points where trajectories decay from S₁ to S₀ at CASPT2 level while the green dots show the corresponding CASSCF decay points. In contrast with the results of the FC trajectories (see above), both the CASPT2 and CASSCF trajectories of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M decay to S₀ (i.e. we do not observe qualitative differences between the two QM levels). Furthermore, the diagram in Fig. 8 indicates that the CASSCF and CASPT2 sets decay on similar timescales for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models. On the other hand, consistently with the FC trajectory description in the previous subsection, the room temperature trajectories confirm that the ASR<sub>AT</sub>-M decays slower than ASR<sub>13C</sub>-M. In fact, the corresponding average hop times were 69 fs and 54 fs at the CASSCF level and 102 fs and 52 fs at the CASPT2 level for ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M respectively.

It is important to establish if the bicycle pedal reaction coordinate followed by the 0 K FC trajectories is affected by the room temperature condition. Accordingly, in Fig. 8 we also display the twisting deformation of the C13=C14 and C15=N double bonds along the trajectory sets. It is shown that, in both ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M, these dihedrals follow the same clockwise and counterclockwise motion observed in the corresponding FC trajectories irrespective of the level of QM theory employed. In particular, it is shown that the ASR<sub>AT</sub>-M reactive CASPT2 trajectories, which are not described by the corresponding FC trajectory, still describe a bicycle pedal motion. These results indicate that the specific construction of the reduced chromophore cavity (having part of the residues frozen as in the original chromophore cavity) has not perturbed the selectivity and isomerization mechanism of the photoisomerization.
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Dihedral of central double bond (red part of the molecule) for ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) at CASPT2 (red) and CASSCF (black) level. The evolution of dihedral angles along the FC trajectories at CASSCF (FC-CASSCF) and CASPT2 (FC-CASPT2) levels are shown in purple and brown lines respectively. The blue dots show the points where the CASPT2 trajectories decay from S<sub>1</sub> to S<sub>0</sub> while the green dots show corresponding CASSCF trajectory decay points. (C and D) The evolution of dihedrals C13=C14 (red) and C15=N (green) of the two modified retinal models along the CASPT2 trajectories. The same dihedrals of the corresponding CASPT2 FC trajectories are shown in orange and light green lines. The blue and brown dots show the decay points along the C13=C14 (red) and C15=N (green) dihedral evolution respectively.</figcaption></figure>

The results above suggest that, for the present systems, one
may use the CASSCF level to investigate the S<sub>1</sub> dynamics of a
room temperature ensemble of QM/MM models and obtain
S<sub>1</sub> lifetimes not too different from those obtained by employ-
ing the superior CASPT2 level. It is interesting to note that the
CASSCF FC trajectories (see Fig. 8A and 8B), correctly describe
the S<sub>1</sub> evolution of, roughly, the center of the room tempera-
ture population<sup>33</sup>. This is also true for the ASR<sub>13C</sub>-M reactive
CASPT2 FC trajectory that follows the corresponding room
temperature trajectory set closely. The ASR<sub>AT</sub>-M unreactive FC
trajectory deviates more substantially from the set but still
describes a tendency to decay at longer delay times. Indeed,

notice the presence in the ASR<sub>AT</sub>-M set (see Fig. 8A) of two
CASPT2 trajectories hopping at longer times (ca. 140 and 170
fs) while in ASR<sub>13C</sub> all trajectories hop below 100 fs.

**3.3 Room temperature excited state dynamics at the CASSCF level.** Obviously a ten-trajectory set is insufficient to produce statistically meaningful results and a much larger set of trajectories is required to derived quantities such as quantum yields and converged excited state lifetime parameters. As we will discuss below there are two parts of the population dynamics calculation affecting the statistical significance of the results: the Boltzmann-like distribution computation and the
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Dihedral of central double bond (red part of the molecule) for ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) at CASPT2 (red) and CASSCF (black) level. The evolution of dihedral angles along the FC trajectories at CASSCF (FC-CASSCF) and CASPT2 (FC-CASPT2) levels are shown in purple and brown lines respectively. The blue dots show the points where the CASPT2 trajectories decay from S<sub>1</sub> to S<sub>0</sub> while the green dots show corresponding CASSCF trajectory decay points. (C and D) The evolution of dihedrals C13=C14 (red) and C15=N (green) of the two modified retinal models along the CASPT2 trajectories. The same dihedrals of the corresponding CASPT2 FC trajectories are shown in orange and light green lines. The blue and brown dots show the decay points along the C13=C14 (red) and C15=N (green) dihedral evolution respectively.</figcaption></figure>

The results above suggest that, for the present systems, one
may use the CASSCF level to investigate the S<sub>1</sub> dynamics of a
room temperature ensemble of QM/MM models and obtain
S<sub>1</sub> lifetimes not too different from those obtained by employ-
ing the superior CASPT2 level. It is interesting to note that the
CASSCF FC trajectories (see Fig. 8A and 8B), correctly describe
the S<sub>1</sub> evolution of, roughly, the center of the room tempera-
ture population<sup>33</sup>. This is also true for the ASR<sub>13C</sub>-M reactive
CASPT2 FC trajectory that follows the corresponding room
temperature trajectory set closely. The ASR<sub>AT</sub>-M unreactive FC
trajectory deviates more substantially from the set but still
describes a tendency to decay at longer delay times. Indeed,

notice the presence in the ASR<sub>AT</sub>-M set (see Fig. 8A) of two
CASPT2 trajectories hopping at longer times (ca. 140 and 170
fs) while in ASR<sub>13C</sub> all trajectories hop below 100 fs.

**3.3 Room temperature excited state dynamics at the CASSCF level.** Obviously a ten-trajectory set is insufficient to produce statistically meaningful results and a much larger set of trajectories is required to derived quantities such as quantum yields and converged excited state lifetime parameters. As we will discuss below there are two parts of the population dynamics calculation affecting the statistical significance of the results: the Boltzmann-like distribution computation and the
surface-hop decision. The initial conditions are generated via molecular dynamics simulations (see Section 2) that require an initial random seed and, similarly, during each Tully surface-hop trajectory a random number is required²⁹,³⁰ whenever a hop decision has to be taken.

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 9.** The evolution of the central dihedral (red part of the molecule) of ensembles of 400 trajectories for ASR<sub>AT</sub>-M (A) and ASR<sub>13C</sub>-M (B) at CASSCF level (black). The evolution of dihedral angles along the FC trajectories at CASSCF (FC-CASSCF) level is shown in purple line. The decay of the populations of ASR<sub>AT</sub>-M (C) and ASR<sub>13C</sub>-M (D) containing 400 molecules (red line), Gaussian exponential fitting (black line), exponential fitting reported by Ruckenbauer et al.<sup>16</sup> (blue dashed line).</figcaption></figure>

Here the room temperature, population dynamics for the photoisomerization of ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M is investigated using ensembles of 400 trajectories. The corresponding results, obtained exclusively at the CASSCF level, are displayed in Fig. 9A and 9B. In Table 2 we report the corresponding

quantum yields, computed as fractions of trajectories leading to photoproduct, and kinetic parameters computed by fitting the population decay with the following function:

$$f(t) = a_1 e^{\left(-\frac{(t-t_1)^2}{t_2^2}\right)} + y_0 + \left((1-a_1) - y_0\right) e^{\left(-\frac{(t-t_1)^2}{t_3^2}\right)} \quad (1)$$

where $a_1$ is the weight of the Gaussian decay component, $y_0$ is a fitting parameter, $t_1$ is the latency time (i.e. the time that precedes the decay onset, $t_2$ and $t_3$ are the decay times of the Gaussian and exponential components respectively. As shown in Fig. 9C and 9D the Gaussian plus exponential function $f(t)$ allows a better fit with respect to the more common bi-exponential function¹⁶. In fact, using eq. 1 the decay dynamics is correctly fitted with $\chi^2$ values 0.033 and 0.038 respectively. The conventional bi-exponential decay model would fit the computed decay dynamics less satisfactorily having corresponding $\chi^2$ values of 0.256 and 0.269 which are higher than that of the Gaussian/exponential fit (smaller the $\chi^2$ value, better the fitting).

As also discussed in Schapiro et al.<sup>15</sup> a combination of Gaussian and exponential fitting could be a wise choice when describing the non-exponential decay of the ASR population. The physical meaning of these terms can be viewed as follows. The population remains on $S_1$ for a time $t_1$ before it starts to decay. A part of the population moves toward the CI quickly and decays “impulsively” according to a Gaussian function so that the corresponding “timescale” is the Gaussian decay time ($t_2$). The rest of the population reaches the $S_1/S_0$ CI more diffusively and decays exponentially (corresponds to exponential decay time $t_3$). The contribution of Gaussian ($a_1$) and exponential decay ($1-a_1$) to the total population decay merely depends on the nature of the $S_1$ potential energy surface. It is expected that the steeper is the surface the higher is the Gaussian weight. The “total” excited state lifetime ($\tau$) is calculated as:

$$\tau = t_1 + a_1 t_2 + (1 - a_1 - y_0) t_3 \quad (2)$$

$\tau$ is therefore different from $\tau_{FC}$ used for associating a $S_1$ lifetime to a FC trajectory which, in turn, can be compared to the $t_1 + a_1 t_2$ time. The statistical convergence of the quantum yields and kinetic parameters $a_1$, $\tau$, $t_1$, $t_2$ and $t_3$ obtained by fitting with the function $f(t)$, is evaluated by plotting the corresponding values as a function of the number of trajectories.

For this reason, we have computed, for the ASR<sub>AT</sub>-M model, two independently generated sets of 400 initial conditions (i.e. generating twice the Boltzmann-like distribution and ensuring that the molecular dynamics simulation step uses, for the two cases, different seeds). We have then compared the values obtained by computing parallel sets of 100, 200, 300 and 400 QM/MM non-adiabatic trajectories. In the following such a test is called “full-replica test”. Also, in order to isolate the random effect introduced by the Tully surface-hop method (see Method section), the QM/MM non-adiabatic trajectories were computed twice starting from the same Boltzmann-like distribution. In the following such a convergence test is indicated as “Tully-replica test”.
surface-hop decision. The initial conditions are generated via molecular dynamics simulations (see Section 2) that require an initial random seed and, similarly, during each Tully surface-hop trajectory a random number is required²⁹,³⁰ whenever a hop decision has to be taken.

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 9.** The evolution of the central dihedral (red part of the molecule) of ensembles of 400 trajectories for ASR<sub>AT-M</sub> (A) and ASR<sub>13C-M</sub> (B) at CASSCF level (black). The evolution of dihedral angles along the FC trajectories at CASSCF (FC-CASSCF) level is shown in purple line. The decay of the populations of ASR<sub>AT-M</sub> (C) and ASR<sub>13C-M</sub> (D) containing 400 molecules (red line), Gaussian exponential fitting (black line), exponential fitting reported by Ruckenbauer et al.<sup>16</sup> (blue dashed line).</figcaption></figure>

Here the room temperature, population dynamics for the photoisomerization of ASR<sub>AT-M</sub> and ASR<sub>13C-M</sub> is investigated using ensembles of 400 trajectories. The corresponding results, obtained exclusively at the CASSCF level, are displayed in Fig. 9A and 9B. In Table 2 we report the corresponding

quantum yields, computed as fractions of trajectories leading to photoproduct, and kinetic parameters computed by fitting the population decay with the following function:

$$f(t) = a_1 e^{\left(-\frac{(t-t_1)^2}{t_2}\right)} + y_0 + \left( (1-a_1) - y_0 \right) e^{\left(\frac{-(t-t_1)}{t_3}\right)} \quad (1)$$

where $a_1$ is the weight of the Gaussian decay component, $y_0$ is a fitting parameter, $t_1$ is the latency time (i.e. the time that precedes the decay onset, $t_2$ and $t_3$ are the decay times of the Gaussian and exponential components respectively. As shown in Fig. 9C and 9D the Gaussian plus exponential function $f(t)$ allows a better fit with respect to the more common bi-exponential function¹⁶. In fact, using eq. 1 the decay dynamics is correctly fitted with $\chi^2$ values 0.033 and 0.038 respectively. The conventional bi-exponential decay model would fit the computed decay dynamics less satisfactorily having corresponding $\chi^2$ values of 0.256 and 0.269 which are higher than that of the Gaussian/exponential fit (smaller the $\chi^2$ value, better the fitting).

As also discussed in Schapiro et al.<sup>15</sup> a combination of Gaussian and exponential fitting could be a wise choice when describing the non-exponential decay of the ASR population. The physical meaning of these terms can be viewed as follows. The population remains on $S_1$ for a time $t_1$ before it starts to decay. A part of the population moves toward the CI quickly and decays “impulsively” according to a Gaussian function so that the corresponding “timescale” is the Gaussian decay time ($t_2$). The rest of the population reaches the $S_1/S_0$ CI more diffusively and decays exponentially (corresponds to exponential decay time $t_3$). The contribution of Gaussian ($a_1$) and exponential decay ($1-a_1$) to the total population decay merely depends on the nature of the $S_1$ potential energy surface. It is expected that the steeper is the surface the higher is the Gaussian weight. The “total” excited state lifetime ($\tau$) is calculated as:

$$\tau = t_1 + a_1 t_2 + (1 - a_1 - y_0) t_3 \quad (2)$$

$\tau$ is therefore different from $\tau_{FC}$ used for associating a $S_1$ lifetime to a FC trajectory which, in turn, can be compared to the $t_1 + a_1 t_2$ time. The statistical convergence of the quantum yields and kinetic parameters $a_1$, $\tau$, $t_1$, $t_2$ and $t_3$ obtained by fitting with the function $f(t)$, is evaluated by plotting the corresponding values as a function of the number of trajectories.

For this reason, we have computed, for the ASR<sub>AT-M</sub> model, two independently generated sets of 400 initial conditions (i.e. generating twice the Boltzmann-like distribution and ensuring that the molecular dynamics simulation step uses, for the two cases, different seeds). We have then compared the values obtained by computing parallel sets of 100, 200, 300 and 400 QM/MM non-adiabatic trajectories. In the following such a test is called “full-replica test”. Also, in order to isolate the random effect introduced by the Tully surface-hop method (see Method section), the QM/MM non-adiabatic trajectories were computed twice starting from the same Boltzmann-like distribution. In the following such a convergence test is indicated as “Tully-replica test”.
<table><caption>**Table 2.** The latency time ($t_1$), Gaussian decay time ($t_2$), exponential decay time ($t_3$), excited state lifetime ($\tau$), weight of the Gaussian decay ($a_1$), quantum yield ($\Phi$) of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M models. The hop times of the Franck-Condon trajectories at CASSCF level ($\tau_{FC}$) are also reported.</caption><thead><tr><th rowspan="2">Model</th><th colspan="8">ASR-M CASSCF</th></tr><tr><th>a<sub>1</sub></th><th>t<sub>2</sub> (fs)</th><th>t<sub>3</sub> (fs)</th><th>t<sub>3</sub> (fs)</th><th>τ (fs)</th><th>χ<sup>2</sup></th><th>t<sub>1</sub> + a<sub>1</sub> t<sub>2</sub> (fs)</th><th>Φ (%)</th></tr></thead><tbody><tr><td>AT</td><td>0.465</td><td>44</td><td>27</td><td>42</td><td>80</td><td>0.038</td><td>57 (61)<sup>1</sup></td><td>57</td></tr><tr><td>13C</td><td>0.554</td><td>32</td><td>28</td><td>33</td><td>62</td><td>0.033</td><td>47 (53)<sup>1</sup></td><td>64</td></tr></tbody></table>

1. $\tau_{FC}$ are given in parenthesis.

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 10.** The variation of quantum yields of full replica (A) and Tully replica (B) against the ensemble size for the ASR<sub>AT</sub>-M model. The reference quantum yield of the ensemble size of 400 trajectories is shown in brackets. The latency time ($t_1$), Gaussian decay time ($t_2$), exponential decay time ($t_3$), excited state lifetime ($\tau$), weight of the Gaussian decay ($a_1$) of full replica (C) and Tully replica (D).</figcaption></figure>
<table><caption>**Table 2.** The latency time ($t_1$), Gaussian decay time ($t_2$), exponential decay time ($t_3$), excited state lifetime ($\tau$), weight of the Gaussian decay ($a_1$), quantum yield ($\Phi$) of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M models. The hop times of the Franck-Condon trajectories at CASSCF level ($\tau_{FC}$) are also reported.</caption><thead><tr><th rowspan="2">Model</th><th colspan="8">ASR-M CASSCF</th></tr><tr><th>a<sub>1</sub></th><th>t<sub>2</sub> (fs)</th><th>t<sub>3</sub> (fs)</th><th>t<sub>3</sub> (fs)</th><th>τ (fs)</th><th>χ<sup>2</sup></th><th>t<sub>1</sub> + a<sub>1</sub> t<sub>2</sub> (fs)</th><th>Φ (%)</th></tr></thead><tbody><tr><td>AT</td><td>0.465</td><td>44</td><td>27</td><td>42</td><td>80</td><td>0.038</td><td>57 (61)<sup>1</sup></td><td>57</td></tr><tr><td>13C</td><td>0.554</td><td>32</td><td>28</td><td>33</td><td>62</td><td>0.033</td><td>47 (53)<sup>1</sup></td><td>64</td></tr></tbody></table>

1. $\tau_{FC}$ are given in parenthesis.

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 10.** The variation of quantum yields of full replica (A) and Tully replica (B) against the ensemble size for the ASR<sub>AT</sub>-M model. The reference quantum yield of the ensemble size of 400 trajectories is shown in brackets. The latency time ($t_1$), Gaussian decay time ($t_2$), exponential decay time ($t_3$), excited state lifetime ($\tau$), weight of the Gaussian decay ($a_1$) of full replica (C) and Tully replica (D).</figcaption></figure>
The results of the full-replica and Tully-replica tests are reported in the panels of Fig. 10. In the full-replica test the differences in quantum yields computed for the two independent trajectory sets (see Fig. 10A) oscillates within 5% and appears not to converge (i.e. shift to a zero value). Thus, a sample of 400 trajectories appears to be insufficient to reach a converged quantum yield prediction within, roughly, a 5% precision. In contrast, the Tully-replica test (see Fig. 10C) shows quantum yields with a tendency to converge when moving from 100 to 400 trajectories. Clearly, full convergence requires a larger number of trajectories in spite of the constrained molecular environment (i.e. outside the chromophore cavity) and, therefore, limited number of degrees of freedom with respect to the number of the degrees of freedom of the full protein.

Note that, in the full-replica test the difference in exponential decay time $t_3$ is exceptionally high (~35 fs) when a 100-trajectory ensemble is used pointing to an inability of such set to produce converged time scales. The convergence of all time scales ($t_1, t_2, t_3$ and $\tau$) and weight ($a_1$) improves when moving to the 400-trajectory ensemble leading to a maximum 15 fs uncertainty and 0.4 differences in $a_1$ value. The variations of lifetimes in the Tully-replica test are smaller (below 5 fs) indicating a closer consistency even if there is no regular improvement with increasing of the ensemble size indicating, again, difficulties in the convergence. We conclude that a 400-trajectory ensemble provides only a semi-quantitative description of the population dynamics with expected variations above 5% for quantum yields, 15 fs for the timescales and 0.4 for the Gaussian weight. We also conclude that such variations are mostly due to the randomness inserted by the initial condition computations (i.e. during the computation of the Gaussian-like distribution).

With the limitation discussed above the 400-trajectory ensemble of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M show (see Table 2 and Fig. 9, top) that the first model decays faster ($\tau$ = 62 fs) than the second ($\tau$ = 80 fs). This, again, confirms the faster dynamics of the ASR<sub>13C</sub>-M model with respect to the ASR<sub>AT</sub>-M model and appears to be qualitatively consistent with the trend observed for the ASR<sub>13C</sub> and ASR<sub>AT</sub> FC trajectories. Notice that the difference between the computed timescales (22 fs) is above the 15 fs variation. The faster ASR<sub>13C</sub>-M decay is justified by a shorter latency time ($t_1$). This reflects a faster motion of the initial ASR<sub>13C</sub>-M population out of the FC region and, presumably a more effective coupling with the torsional motion characterizing the evolution towards the CI. Note that the ASR<sub>13C</sub>-M model decay shows a larger Gaussian weight ($a_1$ = 0.554) whereas in ASR<sub>AT</sub>-M, the decay is more exponential like. As we also mentioned above, we attribute the Gaussian time scale to the time scale ($t_2$) associated with the impulsive population decay and the exponential time scale ($t_3$) with a following non-impulsive decay process. The faster population decay of ASR<sub>13C</sub>-M and comparatively slower decay of ASR<sub>AT</sub>-M are apparent when looking at the evolution of the central dihedral angle (see Fig. 9A and 9B). These differences are accompanied by a lesser spreading of hop times in ASR<sub>13C</sub>-M indicating a $S_1$ population decay more coherent than in the ASR<sub>AT</sub>-M case.

In Fig. 9A and 9B we also superimpose, for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models, the corresponding CASSCF FC trajectories. It is clear that the evolution of the $S_1$ population and its impulsive decay (roughly, the evolution of the population center) is correctly described by the FC trajectories. These not only capture the average behavior of the 400 trajectories but also the trend in the decay time. More precisely the $\tau_{FC}$ of the FC trajectories (53 fs and 61 fs for ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M respectively) compare well to the corresponding $t_1+a_1t_2$ timescales of Table 2. This has also the consequence that the higher isomerization speed (or shorter time required to reach the CI) of ASR<sub>13C</sub>-M is determined by its $t_1$ latency time and by its $a_1t_2$ component. In fact, these are the one controlling the (Gaussian) faster and "impulsive" decay timescale.

### 3.4 Have quantum yields a time dependent structure?

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 12.** Time evolution of reactive (red) and non reactive (green) hop points plotted as a function of the variable alpha and beta and of the time. The Newman diagrams provide the definition of the geometrical parameter affecting the overlap between the p-orbitals centered at C13 and C14 and describing the double bond reconstitution. (A) ASR<sub>AT</sub>-M and (B) ASR<sub>13C</sub>-M.</figcaption></figure>

One intriguing result of the population dynamics reported above, is that, in contrast with the $S_1$ lifetime trend, the computed trend for isomerization quantum yields is opposite to the one observed experimentally. While the precise understanding of this result and discrepancy requires an additional research effort, here we report on a heuristic analysis that could provide initial hypotheses. More specifically, while 400 trajectories are not sufficient to deliver quantitative data, the trajectories may be used to get qualitative information on the distribution of reactive and unreactive processes as a function of the values of critical geometrical parameters and, most
The results of the full-replica and Tully-replica tests are reported in the panels of Fig. 10. In the full-replica test the differences in quantum yields computed for the two independent trajectory sets (see Fig. 10A) oscillates within 5% and appears not to converge (i.e. shift to a zero value). Thus, a sample of 400 trajectories appears to be insufficient to reach a converged quantum yield prediction within, roughly, a 5% precision. In contrast, the Tully-replica test (see Fig. 10C) shows quantum yields with a tendency to converge when moving from 100 to 400 trajectories. Clearly, full convergence requires a larger number of trajectories in spite of the constrained molecular environment (i.e. outside the chromophore cavity) and, therefore, limited number of degrees of freedom with respect to the number of the degrees of freedom of the full protein.

Note that, in the full-replica test the difference in exponential decay time $t_3$ is exceptionally high (~35 fs) when a 100-trajectory ensemble is used pointing to an inability of such set to produce converged time scales. The convergence of all time scales ($t_1, t_2, t_3$ and $\tau$) and weight ($a_1$) improves when moving to the 400-trajectory ensemble leading to a maximum 15 fs uncertainty and 0.4 differences in $a_1$ value. The variations of lifetimes in the Tully-replica test are smaller (below 5 fs) indicating a closer consistency even if there is no regular improvement with increasing of the ensemble size indicating, again, difficulties in the convergence. We conclude that a 400-trajectory ensemble provides only a semi-quantitative description of the population dynamics with expected variations above 5% for quantum yields, 15 fs for the timescales and 0.4 for the Gaussian weight. We also conclude that such variations are mostly due to the randomness inserted by the initial condition computations (i.e. during the computation of the Gaussian-like distribution).

With the limitation discussed above the 400-trajectory ensemble of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M show (see Table 2 and Fig. 9, top) that the first model decays faster ($\tau$ = 62 fs) than the second ($\tau$ = 80 fs). This, again, confirms the faster dynamics of the ASR<sub>13C</sub>-M model with respect to the ASR<sub>AT</sub>-M model and appears to be qualitatively consistent with the trend observed for the ASR<sub>13C</sub> and ASR<sub>AT</sub> FC trajectories. Notice that the difference between the computed timescales (22 fs) is above the 15 fs variation. The faster ASR<sub>13C</sub>-M decay is justified by a shorter latency time ($t_1$). This reflects a faster motion of the initial ASR<sub>13C</sub>-M population out of the FC region and, presumably a more effective coupling with the torsional motion characterizing the evolution towards the CI. Note that the ASR<sub>13C</sub>-M model decay shows a larger Gaussian weight ($a_1$ = 0.554) whereas in ASR<sub>AT</sub>-M, the decay is more exponential like. As we also mentioned above, we attribute the Gaussian time scale to the time scale ($t_2$) associated with the impulsive population decay and the exponential time scale ($t_3$) with a following non-impulsive decay process. The faster population decay of ASR<sub>13C</sub>-M and comparatively slower decay of ASR<sub>AT</sub>-M are apparent when looking at the evolution of the central dihedral angle (see Fig. 9A and 9B). These differences are accompanied by a lesser spreading of hop times in ASR<sub>13C</sub>-M indicating a $S_1$ population decay more coherent than in the ASR<sub>AT</sub>-M case.

In Fig. 9A and 9B we also superimpose, for both the ASR<sub>AT</sub>-M and ASR<sub>13C</sub>-M models, the corresponding CASSCF FC trajectories. It is clear that the evolution of the $S_1$ population and its impulsive decay (roughly, the evolution of the population center) is correctly described by the FC trajectories. These not only capture the average behavior of the 400 trajectories but also the trend in the decay time. More precisely the $\tau_{FC}$ of the FC trajectories (53 fs and 61 fs for ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M respectively) compare well to the corresponding $t_1+a_1t_2$ timescales of Table 2. This has also the consequence that the higher isomerization speed (or shorter time required to reach the CI) of ASR<sub>13C</sub>-M is determined by its $t_1$ latency time and by its $a_1t_2$ component. In fact, these are the one controlling the (Gaussian) faster and "impulsive" decay timescale.

### 3.4 Have quantum yields a time dependent structure?

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 12.** Time evolution of reactive (red) and non reactive (green) hop points plotted as a function of the variable alpha and beta and of the time. The Newman diagrams provide the definition of the geometrical parameter affecting the overlap between the p-orbitals centered at C13 and C14 and describing the double bond reconstitution. (A) ASR<sub>AT</sub>-M and (B) ASR<sub>13C</sub>-M.</figcaption></figure>

One intriguing result of the population dynamics reported above, is that, in contrast with the $S_1$ lifetime trend, the computed trend for isomerization quantum yields is opposite to the one observed experimentally. While the precise understanding of this result and discrepancy requires an additional research effort, here we report on a heuristic analysis that could provide initial hypotheses. More specifically, while 400 trajectories are not sufficient to deliver quantitative data, the trajectories may be used to get qualitative information on the distribution of reactive and unreactive processes as a function of the values of critical geometrical parameters and, most
importantly, of time. The results of such analyses are given in Fig. 12 where we display three-dimensional scattered plots of the S<sub>1</sub> to S<sub>0</sub> decay (hop) points as a function of time and of the critical α and β geometrical coordinates for both ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M. Note that these coordinates affect the overlap between the p-orbitals residing on the carbons of the isomerizing double bond and, in turn, the decay time as the hop usually happens when such overlap is close to zero.

The decay points have been colored according to the fate of the corresponding trajectory. Accordingly, hop points corresponding to trajectories evolving towards the photoproduct, and therefore contributing to a higher quantum yield, are marked in red. The hop points leading to the reconstitution of the reactant have hop points marked in green. The diagrams demonstrate that for the ultrafast photo-isomerizations of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M, quantum yields do not build up homogeneously as a function of time. In fact, it can be clearly seen that waves of reactive and non reactive trajectories alternate with periods of ca. 15 fs at least three times before starting to decay, apparently, in a less organized fashion. The first oscillation starts, for both ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M, with a massive production of photoproducts that is rapidly followed by a wave of reactant reconstitution. This first complete reactive-unreactive cycle has a total duration of ca. 30 fs and it is followed by a second period displaying similar features. Notice also that the sets of alternating waves of reactive and non-reactive decays occur at different values of the geometrical parameters and span the time segment dominated by a Gaussian decay. Additional theoretical work is clearly needed to understand the origin of these features and, in turn, of the factors leading to the final quantum yield value.

## 4. Conclusions

Above we have compared the CASSCF and CASPT2 photoisomerization dynamics for QM/MM models of abridged eubacterial rhodopsins hosting the minimal PSB3 chromophores but conserving the original cavity shape. Due to the high cost of CASPT2 gradients the assessment of the validity of the common use of CASSCF gradients for driving trajectory calculations in these systems was based on the comparison of single FC trajectories (i.e. trajectories released from the FC point with no initial kinetic energy) and a set of 10 room temperature non-adiabatic (Tully surface-hop) trajectories. At both the CASSCF and CASPT2 levels, the room temperature trajectories describe the same clockwise aborted bicyclic-pedal mechanism<sup>34</sup>, display not too distant patterns of decay points (although the CASPT2 decay points form a less compact set) and indicate that the ASR<sub>AT</sub>-M has a longer S<sub>1</sub> lifetime with respect to ASR<sub>13C</sub>-M. These results suggest that, at room temperature, the CASSCF gradients provide a reasonable qualitative description of the photoreaction in spite of the undeniable quantitative and qualitative differences reported when comparing CASSCF and CASPT2 potential energy surface structures. These, as also documented in gas-phase studies<sup>35</sup>, include geometrical differences of the FC points, the change in slope of the S<sub>1</sub> potential energy along the isomerization coordinate and the change in topology and topography of the CI region.<sup>20,35</sup>

As far as the validity of FC trajectories (i.e. single trajectories starting from the S<sub>0</sub> equilibrium structure with zero initial velocities) for getting information on the room temperature dynamics is concerned, our results indicate that when the S<sub>1</sub> potential energy surface is flat or feature small barriers, these may not hop within a reasonable timescale (i.e. the average timescale seen in room temperature population dynamics) but may remain in the S<sub>1</sub> state for a long time. In contrast, the FC trajectories seem to provide valid information in cases when the potential energy surface is not flat thus driving the FC trajectories to a prompt hop. We have shown that these "hopping" FC trajectories describe rather closely the motion of the center of the room temperature population evaluated using 400-trajectory ensembles. Of course, after the hop to S<sub>0</sub>, the trajectory necessarily follows one of the possible relaxation branches in the reactive or unreactive direction.

Above we have also tried to assess if QM/MM models of rhodopsins with PSB3 chromophores would be capable to replicate certain qualitative (i. e. mechanistic) features of the S<sub>1</sub> dynamics documented for unabridged rhodopsin models.<sup>6,15</sup> We found that the aborted bicycle-pedal photoisomerization mechanism<sup>34</sup> is replicated. These models also successfully replicated the documented unidirectional clockwise twisting of the C13=C14 double bond. On the other hand, as evident from the comparison of the data in Table 1, the ASR<sub>13C</sub> and ASR<sub>13C</sub>-M FC trajectories in Fig. 4, the trajectory ensemble in Fig. 9 and the comparison of the fitted excited state lifetimes in Table 2, both the absorption maxima and the S<sub>1</sub> lifetime are, with respect to the unabridged and experimental systems, considerably reduced. This is attributed to the decrease in length of the conjugated chain of the chromophore and related increase in the S<sub>1</sub>-S<sub>0</sub> energy gap and slope of the S<sub>1</sub> energy surface. However, in spite of the slope increase, our dynamics studies (FC trajectories and small and large sets of room temperature non-adiabatic trajectories) indicate that ASR<sub>AT</sub>-M has a longer S<sub>1</sub> lifetime with respect to ASR<sub>13C</sub>-M, in line with the trends documented for the full ASR models and with the observed trend. This, in spite of the demonstrated limited convergence of the kinetic parameters with 400 trajectory ensembles.

In contrast, with the conclusions above, the computed trend in isomerization quantum yields indicates that the faster ASR<sub>13C</sub>-M affords a higher quantum yield at variance with the observed experimental trend. On the other hand, the computed trend appears consistent with the reported hypothesis of a direct proportionality<sup>36,37</sup> between the reaction speed and the quantum yield: an idea based on a simple Landau-Zener model<sup>38,39</sup> for the ultrafast decay. In order to further investigate the origin of the computed quantum yield trend we have "decomposed" the ensemble populations in reactive and non-reactive trajectories. The preliminary results of Fig. 12, clearly indicate that the quantum yield value is generated through a relatively complex process displaying regular space and time features. At this stage it is therefore concluded that the familiar Landau-Zener decay picture for ultrafast decay, is therefore too simplistic for even a qualitative discussion of quantum yield trends (at least in double bond isomerizations).
importantly, of time. The results of such analyses are given in Fig. 12 where we display three-dimensional scattered plots of the S<sub>1</sub> to S<sub>0</sub> decay (hop) points as a function of time and of the critical α and β geometrical coordinates for both ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M. Note that these coordinates affect the overlap between the p-orbitals residing on the carbons of the isomerizing double bond and, in turn, the decay time as the hop usually happens when such overlap is close to zero.

The decay points have been colored according to the fate of the corresponding trajectory. Accordingly, hop points corresponding to trajectories evolving towards the photoproduct, and therefore contributing to a higher quantum yield, are marked in red. The hop points leading to the reconstitution of the reactant have hop points marked in green. The diagrams demonstrate that for the ultrafast photo-isomerizations of ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M, quantum yields do not build up homogeneously as a function of time. In fact, it can be clearly seen that waves of reactive and non reactive trajectories alternate with periods of ca. 15 fs at least three times before starting to decay, apparently, in a less organized fashion. The first oscillation starts, for both ASR<sub>13C</sub>-M and ASR<sub>AT</sub>-M, with a massive production of photoproducts that is rapidly followed by a wave of reactant reconstitution. This first complete reactive-unreactive cycle has a total duration of ca. 30 fs and it is followed by a second period displaying similar features. Notice also that the sets of alternating waves of reactive and non-reactive decays occur at different values of the geometrical parameters and span the time segment dominated by a Gaussian decay. Additional theoretical work is clearly needed to understand the origin of these features and, in turn, of the factors leading to the final quantum yield value.

## 4. Conclusions

Above we have compared the CASSCF and CASPT2 photoisomerization dynamics for QM/MM models of abridged eubacterial rhodopsins hosting the minimal PSB3 chromophores but conserving the original cavity shape. Due to the high cost of CASPT2 gradients the assessment of the validity of the common use of CASSCF gradients for driving trajectory calculations in these systems was based on the comparison of single FC trajectories (i.e. trajectories released from the FC point with no initial kinetic energy) and a set of 10 room temperature non-adiabatic (Tully surface-hop) trajectories. At both the CASSCF and CASPT2 levels, the room temperature trajectories describe the same clockwise aborted bicyclic-pedal mechanism<sup>34</sup>, display not too distant patterns of decay points (although the CASPT2 decay points form a less compact set) and indicate that the ASR<sub>AT</sub>-M has a longer S<sub>1</sub> lifetime with respect to ASR<sub>13C</sub>-M. These results suggest that, at room temperature, the CASSCF gradients provide a reasonable qualitative description of the photoreaction in spite of the undeniable quantitative and qualitative differences reported when comparing CASSCF and CASPT2 potential energy surface structures. These, as also documented in gas-phase studies<sup>35</sup>, include geometrical differences of the FC points, the change in slope of the S<sub>1</sub> potential energy along the isomerization coordinate and the change in topology and topography of the CI region.<sup>20,35</sup>

As far as the validity of FC trajectories (i.e. single trajectories starting from the S<sub>0</sub> equilibrium structure with zero initial velocities) for getting information on the room temperature dynamics is concerned, our results indicate that when the S<sub>1</sub> potential energy surface is flat or feature small barriers, these may not hop within a reasonable timescale (i.e. the average timescale seen in room temperature population dynamics) but may remain in the S<sub>1</sub> state for a long time. In contrast, the FC trajectories seem to provide valid information in cases when the potential energy surface is not flat thus driving the FC trajectories to a prompt hop. We have shown that these "hopping" FC trajectories describe rather closely the motion of the center of the room temperature population evaluated using 400-trajectory ensembles. Of course, after the hop to S<sub>0</sub>, the trajectory necessarily follows one of the possible relaxation branches in the reactive or unreactive direction.

Above we have also tried to assess if QM/MM models of rhodopsins with PSB3 chromophores would be capable to replicate certain qualitative (i. e. mechanistic) features of the S<sub>1</sub> dynamics documented for unabridged rhodopsin models.<sup>6,15</sup> We found that the aborted bicycle-pedal photoisomerization mechanism<sup>34</sup> is replicated. These models also successfully replicated the documented unidirectional clockwise twisting of the C13=C14 double bond. On the other hand, as evident from the comparison of the data in Table 1, the ASR<sub>13C</sub> and ASR<sub>13C</sub>-M FC trajectories in Fig. 4, the trajectory ensemble in Fig. 9 and the comparison of the fitted excited state lifetimes in Table 2, both the absorption maxima and the S<sub>1</sub> lifetime are, with respect to the unabridged and experimental systems, considerably reduced. This is attributed to the decrease in length of the conjugated chain of the chromophore and related increase in the S<sub>1</sub>-S<sub>0</sub> energy gap and slope of the S<sub>1</sub> energy surface. However, in spite of the slope increase, our dynamics studies (FC trajectories and small and large sets of room temperature non-adiabatic trajectories) indicate that ASR<sub>AT</sub>-M has a longer S<sub>1</sub> lifetime with respect to ASR<sub>13C</sub>-M, in line with the trends documented for the full ASR models and with the observed trend. This, in spite of the demonstrated limited convergence of the kinetic parameters with 400 trajectory ensembles.

In contrast, with the conclusions above, the computed trend in isomerization quantum yields indicates that the faster ASR<sub>13C</sub>-M affords a higher quantum yield at variance with the observed experimental trend. On the other hand, the computed trend appears consistent with the reported hypothesis of a direct proportionality<sup>36,37</sup> between the reaction speed and the quantum yield: an idea based on a simple Landau-Zener model<sup>38,39</sup> for the ultrafast decay. In order to further investigate the origin of the computed quantum yield trend we have "decomposed" the ensemble populations in reactive and non-reactive trajectories. The preliminary results of Fig. 12, clearly indicate that the quantum yield value is generated through a relatively complex process displaying regular space and time features. At this stage it is therefore concluded that the familiar Landau-Zener decay picture for ultrafast decay, is therefore too simplistic for even a qualitative discussion of quantum yield trends (at least in double bond isomerizations).
The conclusions reported above suggest that, with the described limitations, suitable QM/MM models of rhodopsins (including CASSCF/6-31G*/AMBER models) with reduced retinal chromophores represent potentially useful tools in computational photobiology and should be further investigated. In fact, due to the affordable computational cost, these can be used to rapidly and systematically scan for desirable photoisomerization properties in rhodopsins from different organisms or mutants of a specific rhodopsin. For instance, it could be possible to examine hundreds of ASR mutants to search for novel optogenetic tools displaying suppressed isomerization quantum yields, longer excited state lifetimes and, in turn, enhanced fluorescence properties. This screening would help to limit most demanding computations to the most promising mutants for which QM/MM models with full chromophores would be constructed and studied. These possibilities represent the main motivation for the reported exploratory study.

## ASSOCIATED CONTENT

Supporting text and Cartesian coordinates of the optimized structures. This material is available free of charge via the Internet at http://pubs.acs.org.

## ACKNOWLEDGMENT

We are grateful to Dr. Yoelvis Orozco for assistance with the generation of the initial conditions. M.O. is grateful to the Center for Photochemical Sciences of the Bowling Green State University, the Human Frontier Science Program Organization under Grant RGP0049/385, the National Science Foundation for grants No. CHE-1152070 and CHE-1551416 the Institute for Advanced Studies of the University of Strasbourg for a USIAS fellowship. M.O. is also grateful for partial funding from the EU-FP7 (Marie-Curie PIOF-GA-2012-332233). L.M.F. acknowledges the Spanish MINECO grant CTQ2012-36966 and UAH (Universidad de Alcalà) grants CCG2013/EXP-089 and CCG2014/EXP-083. A.V. acknowledges UAH and Università di Siena for doctoral fellowships. The research was supported in part by the National Science Foundation through XSEDE resources provided by the XSEDE Science Gateways program and by the Ohio Supercomputing Center (OSC).

## References

(1) Spudich, J. L.; Yang, C. S.; Jung, K. H.; Spudich, E. N. *Annu. Rev. Cell. Dev. Biol.* **2000**, *16*, 365-92.

(2) Kandori, H.; Shichida, Y.; Yoshizawa, T. *Biochemistry* (*Moscow*) **2001**, *66*, 1197-209.

(3) Ernst, O. P.; Lodowski, D. T.; Elstner, M.; Hegemann, P.; Brown, L. S.; Kandori, H. *Chem. Rev.* **2014**, *114*, 126-63.

(4) Weingart, O.; Altoè, P.; Stenta, M.; Bottoni, A.; Orlandi, G.; Garavelli, M. *Phys. Chem. Chem. Phys.* **2011**, *13*, 3645-8.

(5) Léonard, J.; Schapiro, I.; Briand, J.; Fusi, S.; Paccani, R. R.; Olivucci, M.; Haacke, S. *Chemistry* **2012**, *18*, 15296-304.

(6) Strambi, A.; Durbeej, B.; Ferré, N.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2010**, *107*, 21322-6.

(7) Deisseroth, K. *Nat. Methods* **2010**, *8*, 26-9.

(8) Roos, B. O.; Bruna, P.; Peyerimhoff, S. D.; Shepard, R.; Cooper, D. L.; Gerratt, J.; Raimondi, M. *Adv. Chem. Phys.* **1987**, *69*, 399-446.

(9) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. *J. Phys. Chem.* **1990**, *94*, 5483-8.

(10) Frutos, L. M.; Andruniów, T.; Santoro, F.; Ferré, N.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2007**, *104*, 7764-9.

(11) Kawanabe, A.; Furutani, Y.; Jung, K. H.; Kandori, H. *J. Am. Chem. Soc.* **2007**, *129*, 8644-9.

(12) Cheminal, A.; Léonard, J.; Kim, S. Y.; Jung, K. -; Kandori, H.; Haacke, S. *Chem. Phys. Lett.* **2013**, *587*, 75-80.

(13) Cheminal, A.; Leonard, J.; Kim, S. -Y.; Jung, K. -H.; Kandori, H.; Haacke, S. *Phys. Chem. Chem. Phys.* **2015**, *17*, 25429-25439.

(14) Wand, A.; Rozin, R.; Eliash, T.; Jung, K. H.; Sheves, M.; Ruhman, S. *J. Am. Chem. Soc.* **2011**, *133*, 20922-32.

(15) Schapiro, I.; Ruhman, S. *Biochim. Biophys. Acta.* **2014**, *1837*, 589-97.

(16) Ruckenbauer, M.; Barbatti, M.; Müller, T.; Lischka, H. *J. Phys. Chem. A.* **2013**, *117*, 2790-9.

(17) Page, C. S.; Olivucci, M. *J. Comput. Chem.* **2003**, *24*, 298-309.

(18) Valsson, O.; Filippi, C. *J. Chem. Theory. Comput.* **2010**, *6*, 1275-1292.

(19) Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. *J. Chem. Theory. Comput.* **2012**, *8*, 4069-80.

(20) Gozem, S.; Melaccio, F.; Lindh, R.; Krylov, A. I.; Granovsky, A. A.; Angeli, C.; Olivucci, M. *J. Chem. Theory. Comput.* **2013**, *9*, 4495-506.

(21) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P. A.; Neogrády, P.; Pedersen, T. B.; Pitonák, M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. *J. Comput. Chem.* **2010**, *31*, 224-47.

(22) Ponder, J. W.; Richards, F. M. *J. Comput. Chem.* **1987**, *8*, 1016-24.

(23) Gozem, S.; Melaccio, F.; Luk, H. L.; Rinaldi, S.; Olivucci, M. *Chem. Soc. Rev.* **2014**, *43*, 4019-36.

(24) Ferré, N.; Ángyán, J. G. *Chem. Phys. Lett.* **2002**, *356*, 331-9.

(25) Preston, R. K.; Tully, J. C. *J. Chem. Phys.* **1971**, *54*, 4297-304.

(26) Tully, J. C.; Preston, R. K. *J. Chem. Phys.* **1971**, *55*, 562-72.

(27) Vreven, T.; Bernardi, F.; Garavelli, M.; Olivucci, M.; Robb, M. A.; Schlegel, H. B. *J. Am. Chem. Soc.* **1997**, *119*, 12687-8.

(28) Weingart, O.; Migani, A.; Olivucci, M.; Robb, M. A.; Buss, V.; Hunt, P. *The Journal of Physical Chemistry. A.* **2004**, *108*, 4685-93.

(29) Granucci, G.; Persico, M. *J. Chem. Phys.* **2007**, *126*, 134114.

(30) Tully, J. C. *J. Chem. Phys.* **1990**, *93*, 1061-71.

(31) Polli, D.; Altoè, P.; Weingart, O.; Spillane, K. M.; Manzoni, C.; Brida, D.; Tomasello, G.; Orlandi, G.; Kukura, P.; Mathies, R. A.; Garavelli, M.; Cerullo, G. *Nature* **2010**, *467*, 440-3.

(32) Sinicropi, A.; Migani, A.; Vico, L. D.; Olivucci, M. *Photochem. Photobiol. Sci.* **2003**, *2*, 1250-5.

(33) Rinaldi, S.; Melaccio, F.; Gozem, S.; Fanelli, F.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2014**, *111*, 1714-9.

(34) Warshel, A. *Nature* **1976**, *260*, 679-83.

(35) Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; Huix-Rotllant, M.; Ferré, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A. *J. Chem. Theory. Comput.* **2014**, *10*, 3074-84.

(36) Weiss, R. M.; Warshel, A. *J. Am. Chem. Soc.* **1979**, *101*, 6131-3.

(37) Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. *Science* **1991**, *254*, 412-415.

(38) Zener, C. *P. R. Soc. Lond. A-Conta.* **1932**, *137*, 696-702.

(39) Kim, J. E.; Tauber, M. J.; Mathies, R. A. *Biophys. J.* **2003**, *84*, 2492-501.
The conclusions reported above suggest that, with the described limitations, suitable QM/MM models of rhodopsins (including CASSCF/6-31G*/AMBER models) with reduced retinal chromophores represent potentially useful tools in computational photobiology and should be further investigated. In fact, due to the affordable computational cost, these can be used to rapidly and systematically scan for desirable photoisomerization properties in rhodopsins from different organisms or mutants of a specific rhodopsin. For instance, it could be possible to examine hundreds of ASR mutants to search for novel optogenetic tools displaying suppressed isomerization quantum yields, longer excited state lifetimes and, in turn, enhanced fluorescence properties. This screening would help to limit most demanding computations to the most promising mutants for which QM/MM models with full chromophores would be constructed and studied. These possibilities represent the main motivation for the reported exploratory study.

## ASSOCIATED CONTENT

Supporting text and Cartesian coordinates of the optimized structures. This material is available free of charge via the Internet at http://pubs.acs.org.

## ACKNOWLEDGMENT

We are grateful to Dr. Yoelvis Orozco for assistance with the generation of the initial conditions. M.O. is grateful to the Center for Photochemical Sciences of the Bowling Green State University, the Human Frontier Science Program Organization under Grant RGP0049/385, the National Science Foundation for grants No. CHE-1152070 and CHE-1551416 the Institute for Advanced Studies of the University of Strasbourg for a USIAS fellowship. M.O. is also grateful for partial funding from the EU-FP7 (Marie-Curie PIOF-GA-2012-332233). L.M.F. acknowledges the Spanish MINECO grant CTQ2012-36966 and UAH (Universidad de Alcalà) grants CCG2013/EXP-089 and CCG2014/EXP-083. A.V. acknowledges UAH and Università di Siena for doctoral fellowships. The research was supported in part by the National Science Foundation through XSEDE resources provided by the XSEDE Science Gateways program and by the Ohio Supercomputing Center (OSC).

## References

(1) Spudich, J. L.; Yang, C. S.; Jung, K. H.; Spudich, E. N. *Annu. Rev. Cell. Dev. Biol.* **2000**, *16*, 365-92.

(2) Kandori, H.; Shichida, Y.; Yoshizawa, T. *Biochemistry* (*Moscow*) **2001**, *66*, 1197-209.

(3) Ernst, O. P.; Lodowski, D. T.; Elstner, M.; Hegemann, P.; Brown, L. S.; Kandori, H. *Chem. Rev.* **2014**, *114*, 126-63.

(4) Weingart, O.; Altoè, P.; Stenta, M.; Bottoni, A.; Orlandi, G.; Garavelli, M. *Phys. Chem. Chem. Phys.* **2011**, *13*, 3645-8.

(5) Léonard, J.; Schapiro, I.; Briand, J.; Fusi, S.; Paccani, R. R.; Olivucci, M.; Haacke, S. *Chemistry* **2012**, *18*, 15296-304.

(6) Strambi, A.; Durbeej, B.; Ferré, N.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2010**, *107*, 21322-6.

(7) Deisseroth, K. *Nat. Methods* **2010**, *8*, 26-9.

(8) Roos, B. O.; Bruna, P.; Peyerimhoff, S. D.; Shepard, R.; Cooper, D. L.; Gerratt, J.; Raimondi, M. *Adv. Chem. Phys.* **1987**, *69*, 399-446.

(9) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. *J. Phys. Chem.* **1990**, *94*, 5483-8.

(10) Frutos, L. M.; Andruniów, T.; Santoro, F.; Ferré, N.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2007**, *104*, 7764-9.

(11) Kawanabe, A.; Furutani, Y.; Jung, K. H.; Kandori, H. *J. Am. Chem. Soc.* **2007**, *129*, 8644-9.

(12) Cheminal, A.; Léonard, J.; Kim, S. Y.; Jung, K. -; Kandori, H.; Haacke, S. *Chem. Phys. Lett.* **2013**, *587*, 75-80.

(13) Cheminal, A.; Leonard, J.; Kim, S. -Y.; Jung, K. -H.; Kandori, H.; Haacke, S. *Phys. Chem. Chem. Phys.* **2015**, *17*, 25429-25439.

(14) Wand, A.; Rozin, R.; Eliash, T.; Jung, K. H.; Sheves, M.; Ruhman, S. *J. Am. Chem. Soc.* **2011**, *133*, 20922-32.

(15) Schapiro, I.; Ruhman, S. *Biochim. Biophys. Acta.* **2014**, *1837*, 589-97.

(16) Ruckenbauer, M.; Barbatti, M.; Müller, T.; Lischka, H. *J. Phys. Chem. A.* **2013**, *117*, 2790-9.

(17) Page, C. S.; Olivucci, M. *J. Comput. Chem.* **2003**, *24*, 298-309.

(18) Valsson, O.; Filippi, C. *J. Chem. Theory. Comput.* **2010**, *6*, 1275-1292.

(19) Gozem, S.; Huntress, M.; Schapiro, I.; Lindh, R.; Granovsky, A. A.; Angeli, C.; Olivucci, M. *J. Chem. Theory. Comput.* **2012**, *8*, 4069-80.

(20) Gozem, S.; Melaccio, F.; Lindh, R.; Krylov, A. I.; Granovsky, A. A.; Angeli, C.; Olivucci, M. *J. Chem. Theory. Comput.* **2013**, *9*, 4495-506.

(21) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P. A.; Neogrády, P.; Pedersen, T. B.; Pitonák, M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. *J. Comput. Chem.* **2010**, *31*, 224-47.

(22) Ponder, J. W.; Richards, F. M. *J. Comput. Chem.* **1987**, *8*, 1016-24.

(23) Gozem, S.; Melaccio, F.; Luk, H. L.; Rinaldi, S.; Olivucci, M. *Chem. Soc. Rev.* **2014**, *43*, 4019-36.

(24) Ferré, N.; Ángyán, J. G. *Chem. Phys. Lett.* **2002**, *356*, 331-9.

(25) Preston, R. K.; Tully, J. C. *J. Chem. Phys.* **1971**, *54*, 4297-304.

(26) Tully, J. C.; Preston, R. K. *J. Chem. Phys.* **1971**, *55*, 562-72.

(27) Vreven, T.; Bernardi, F.; Garavelli, M.; Olivucci, M.; Robb, M. A.; Schlegel, H. B. *J. Am. Chem. Soc.* **1997**, *119*, 12687-8.

(28) Weingart, O.; Migani, A.; Olivucci, M.; Robb, M. A.; Buss, V.; Hunt, P. *The Journal of Physical Chemistry. A.* **2004**, *108*, 4685-93.

(29) Granucci, G.; Persico, M. *J. Chem. Phys.* **2007**, *126*, 134114.

(30) Tully, J. C. *J. Chem. Phys.* **1990**, *93*, 1061-71.

(31) Polli, D.; Altoè, P.; Weingart, O.; Spillane, K. M.; Manzoni, C.; Brida, D.; Tomasello, G.; Orlandi, G.; Kukura, P.; Mathies, R. A.; Garavelli, M.; Cerullo, G. *Nature* **2010**, *467*, 440-3.

(32) Sinicropi, A.; Migani, A.; Vico, L. D.; Olivucci, M. *Photochem. Photobiol. Sci.* **2003**, *2*, 1250-5.

(33) Rinaldi, S.; Melaccio, F.; Gozem, S.; Fanelli, F.; Olivucci, M. *Proc. Natl. Acad. Sci. U. S. A.* **2014**, *111*, 1714-9.

(34) Warshel, A. *Nature* **1976**, *260*, 679-83.

(35) Gozem, S.; Melaccio, F.; Valentini, A.; Filatov, M.; Huix-Rotllant, M.; Ferré, N.; Frutos, L. M.; Angeli, C.; Krylov, A. I.; Granovsky, A. A. *J. Chem. Theory. Comput.* **2014**, *10*, 3074-84.

(36) Weiss, R. M.; Warshel, A. *J. Am. Chem. Soc.* **1979**, *101*, 6131-3.

(37) Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. *Science* **1991**, *254*, 412-415.

(38) Zener, C. *P. R. Soc. Lond. A-Conta.* **1932**, *137*, 696-702.

(39) Kim, J. E.; Tauber, M. J.; Mathies, R. A. *Biophys. J.* **2003**, *84*, 2492-501.
Insert Table of Contents artwork here

<figure><img src="image_27.png" /></figure>
Insert Table of Contents artwork here

<figure><img src="image_27.png" /></figure>