# Ultrafast carbon monoxide photolysis and heme spin-crossover in myoglobin via nonadiabatic quantum dynamics

Konstantin Falahati, Hiroyuki Tamura, Irene Burghardt, Miquel Huix-Rotllant

► To cite this version:

Konstantin Falahati, Hiroyuki Tamura, Irene Burghardt, Miquel Huix-Rotllant. Ultrafast carbon monoxide photolysis and heme spin-crossover in myoglobin via nonadiabatic quantum dynamics. Nature Communications, 2018, 9 (1), 10.1038/s41467-018-06615-1. hal-02079618

HAL Id: hal-02079618

https://amu.hal.science/hal-02079618

Submitted on 26 Mar 2019

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
ARTICLE

DOI: 10.1038/s41467-018-06615-1

OPEN

Ultrafast carbon monoxide photolysis and heme
spin-crossover in myoglobin via nonadiabatic
quantum dynamics

Konstantin Falahati ID <fnref n="1" />, Hiroyuki Tamura², Irene Burghardt ID <fnref n="1" /> & Miquel Huix-Rotllant ID <fnref n="1" />,³

Light absorption of myoglobin triggers diatomic ligand photolysis and a spin crossover transition of iron(II) that initiate protein conformational change. The photolysis and spin crossover reactions happen concurrently on a femtosecond timescale. The microscopic origin of these reactions remains controversial. Here, we apply quantum wavepacket dynamics to elucidate the ultrafast photochemical mechanism for a heme-carbon monoxide (heme-CO) complex. We observe coherent oscillations of the Fe-CO bond distance with a period of 42 fs and an amplitude of ~1 Å. These nuclear motions induce pronounced geometric reorganization, which makes the CO dissociation irreversible. The reaction is initially dominated by symmetry breaking vibrations inducing an electron transfer from porphyrin to iron. Subsequently, the wavepacket relaxes to the triplet manifold in ~75 fs and to the quintet manifold in ~430 fs. Our results highlight the central role of nuclear vibrations at the origin of the ultrafast photodynamics of organometallic complexes.

<footnote n="1">¹Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt, Frankfurt 60438, Germany. ²Department of Chemical System Engineering, The University of Tokyo, Tokyo 113-8656, Japan. ³Aix Marseille Univ, CNRS, ICR, Marseille, France. Correspondence and requests for materials should be addressed to I.B. (email: burghardt@chemie.uni-frankfurt.de) or to M.H.-R. (email: miquel.huixrotllant@univ-amu.fr)</footnote>
Hemeproteins containing a porphyrin-iron complex (heme) play a major role in the storage and transport of diatomic molecules¹. The diatomic ligand dissociation occurs concurrently to a spin-crossover (SCO) transition from low-spin (LS) to high-spin (HS) of the FeII center². The heme complex with carbon monoxide (CO) in myoglobin is one of the most studied hemeproteins of this kind³⁻¹⁸. Heme-CO is initially in a singlet LS state (S = 0), which transforms to a quintet HS state (S = 2) upon CO dissociation¹⁰. The unbound CO in the myoglobin cavity along with the motions of the remaining heme initiate a “protein quake” that opens a channel for the diatomic molecule to escape³. Over a century ago, Halden and Lorrain discovered that the CO dissociation can be initiated photochemically¹⁹, by a mechanism which is still under debate. The photolysis can be initiated upon absorption to the lowest singlet porphyrin band (¹Q) or the second porphyrin band (¹B)²⁰. Here, we focus on the photochemistry of the heme-CO complex upon excitation to the lowest Q-band of porphyrin.

The heme-CO photolysis is an ultrafast process. Recent pump-probe X-ray experiments of myoglobin, with an initial pump pulse exciting the ¹Q state, seem to agree on a two-step kinetic reaction: (i) a first step taking <50–70 fs, attributed to both CO photolysis and partial SCO, and (ii) a completion of the spin transition to the HS state in ~300–400 fs⁵,⁶,²⁰. Despite the numerous studies on heme-CO photolysis, the kinetics and mechanism of dissociation are still under debate, notably regarding the ultrafast nature of the reaction, and the spin and character of the photolytic state. The rate of photolysis has never been experimentally reported, apart from the upper bound of 50–70 fs. As far as the photolytic state is concerned, the most widely accepted hypothesis is that dissociation occurs from a metal-to-ligand charge-transfer (MLCT) state¹⁵⁻¹⁸,²⁰. Still, experiments and theory do not provide a unified picture to date. On the one hand, Franzen, Martin et al., based on time-resolved absorption and Raman experiments, described a “rapid spin state change that must precede photolysis”²⁰. They consider that photolysis is occurring from a triplet metal→porphyrin ring transfer (³MLCT dissociation)²⁰,²¹. In the model by Franzen et al., the ultrafast reaction is due to valence tautomerism. This mechanism implies a rapid interconversion of several quasidegenerate electronic states involving d-transitions, some of which are dissociative for the CO bond. On the other hand, Head-Gordon et al., employing time-dependent density-functional theory, proposed a dissociation in a singlet metal → $\sigma_{CO}^*$ state (¹MLCT dissociation)¹⁵⁻¹⁸. In the model of Head-Gordon et al., CO photolysis occurs in a Marcus-like process, in which the initial ¹Q population is transferred to the photolytic state after crossing a barrier of <0.2 eV. However, this hypothesis seems in contradiction with experimental observations of absence of fluorescence emission, an ultrafast reaction and a unity quantum yield for the heme-CO complex²². The theoretical method used provided an incomplete description of ultrafast heme-CO photolysis, since the coupled electron-nuclear motion was neglected. Recent theoretical studies of similar iron complexes indicate the fundamental role of nuclear motions in explaining the ultrafast nature of intersystem crossing (ISC)²³⁻²⁵. Such electron-nuclear strong coupling has been experimentally observed in time-resolved X-ray absorption in a Fe(II) tris(2,2'-bipyridine) complex²⁶,²⁷ and also in the ferricyanide ion²⁸.

Here, we resolve the mechanism of ultrafast heme-CO photolysis by means of quantum wavepacket dynamics, which accounts for nuclear and electronic coherent motion²⁹. We find that CO photolysis occurs in around 20 fs in the ¹MLCT band, prior to the spin transition. The ultrafast nature of the photolysis is due to strong vibronic couplings and a band of quasi-continuous states. Upon dissociation, the SCO process induces a

sequential transition of the remaining heme first into the triplet ³MLCT and second to the HS quintet manifold ⁵MLCT.

## Results

**Heme-CO model and electronic structure.** The model taken as representative of the active center of the myoglobin protein consists of an iron(II) encapsulated in a porphyrin ring, and axially ligated to a CO and an imidazole. Hereafter this model is referred to as heme-CO for the complex with carbon monoxide and simply heme for the remaining ligands. The LS state of the heme-CO complex has an electronic configuration $(a_2)^2(a_1)^2$ ($e_x^*, e_y^*)^0$ for the porphyrin and $(d_{xy})^2(d_{xz}, d_{yz})^4(d_{z^2})^0$ ($d_{x^2-y^2})^0$ for Fe. In its minimum energy structure, the iron and porphyrin are coplanar and the CO is bound to Fe in an upright position perpendicular to the plane. The Fe-C(O) bond is 1.80 Å, around 0.26 Å shorter than the Fe-N bonds with porphyrin or imidazole, which are 2.03 Å and 2.06 Å, respectively. Upon photolysis, the Fe changes into a quintet state, with electronic configuration $(d_{xy})^1(d_{xz}, d_{yz})^3(d_{z^2})^1$ ($d_{x^2-y^2})^1$. In the HS, the Fe-CO bond is unstable. The iron is displaced from the porphyrin plane, distorted to a square-pyramidal structure. The Fe-N bond length is 2.10 Å for the porphyrin nitrogens and 2.17 Å for imidazole.

The full electronic spectrum of heme-CO in the Franck-Condon region is shown in Table 1. This spectrum has been computed in the gas phase, at the CASSCF(10,9)/CASPT2/ANO-RCC-VDZP level of theory for the minimum energy structure of the S₀ ground state (the active space is shown in Supplementary Fig. 1). More details can be found on the effect on the spectrum of the basis set (Supplementary Table 1), the level of theory (Supplementary Table 2) and the model (Supplementary Tables 3–5). At this geometry, the HS state is found 0.7 eV above the LS state. The bright states are the ¹Q states, which correspond to a $\pi \rightarrow \pi^*$ transition localized on the porphyrin ring. The ¹Qₓ and ¹Qᵧ states have an oscillator strength of 1.37·10⁻² and 9.13·10⁻³ au for Qₓ and Qᵧ, respectively. In our model, these transitions are found at 2.73–2.79 eV, while

<table><caption>Table 1 Electronic spectrum of the heme-CO model in the gas phase</caption><thead><tr><th>State</th><th>Energy</th><th>Character</th><th>Sym.</th></tr></thead><tbody><tr><td>1<sup>5</sup>MC</td><td>0.70</td><td>(d<sub>yz</sub>, d<sub>xy</sub>) → (d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td>E</td></tr><tr><td>2<sup>5</sup>MC</td><td>0.74</td><td>(d<sub>xz</sub>, d<sub>xy</sub>) → (d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td>E</td></tr><tr><td>1<sup>3</sup>MC</td><td>0.89</td><td>d<sub>xy</sub> → d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>B<sub>2</sub></td></tr><tr><td>3<sup>5</sup>MC</td><td>1.40</td><td>(d<sub>xz</sub>, d<sub>yz</sub>) → (d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td></td></tr><tr><td>2<sup>3</sup>MC</td><td>1.43</td><td>d<sub>xz</sub> → d<sub>z<sup>2</sup></sub></td><td>E</td></tr><tr><td>3<sup>3</sup>MC</td><td>1.51</td><td>d<sub>yz</sub> → d<sub>z<sup>2</sup></sub></td><td>E</td></tr><tr><td>4<sup>3</sup>MC</td><td>1.62</td><td>d<sub>xy</sub> → d<sub>z<sup>2</sup></sub></td><td>A<sub>2</sub></td></tr><tr><td>5<sup>3</sup>MC</td><td>1.86</td><td>d<sub>yz</sub> → d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>E</td></tr><tr><td>6<sup>3</sup>MC</td><td>1.90</td><td>d<sub>xz</sub> → d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>E</td></tr><tr><td>1<sup>1</sup>MC</td><td>1.90</td><td>d<sub>xy</sub> → d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>B<sub>2</sub></td></tr><tr><td>1<sup>3</sup>MLCT</td><td>2.15</td><td>(a<sub>1</sub>, d<sub>yz</sub>, d<sub>xy</sub>) → (e, d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td></td></tr><tr><td>2<sup>1</sup>MC</td><td>2.35</td><td>d<sub>yz</sub> → d<sub>z<sup>2</sup></sub></td><td>E</td></tr><tr><td>3<sup>1</sup>MC</td><td>2.36</td><td>d<sub>xy</sub> → d<sub>z<sup>2</sup></sub></td><td>E</td></tr><tr><td>1<sup>3</sup>Q<sub>1</sub></td><td>2.43</td><td>a<sub>2</sub> → e</td><td>E</td></tr><tr><td>5<sup>MLCT</sup></td><td>2.48</td><td>(a<sub>2</sub>, d<sub>yz</sub>, d<sub>xy</sub>) → (e, d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td></td></tr><tr><td>2<sup>3</sup>Q<sub>2</sub></td><td>2.49</td><td>a<sub>2</sub> → e</td><td>E</td></tr><tr><td>1<sup>3</sup>MLCT</td><td>2.56</td><td>(a<sub>1</sub>, d<sub>yz</sub>, d<sub>xy</sub>) → (e, d<sub>x<sup>2</sup>-y<sup>2</sup></sub>, d<sub>z<sup>2</sup></sub>)</td><td></td></tr><tr><td>1<sup>Q</sup><sub>x</sub></td><td>2.73</td><td>a<sub>2</sub> → e/a<sub>1</sub> → e</td><td>E</td></tr><tr><td>1<sup>Q</sup><sub>y</sub></td><td>2.79</td><td>a<sub>2</sub> → e/a<sub>1</sub> → e</td><td>E</td></tr><tr><td>1<sup>MLCT</sup></td><td>2.83</td><td>(a<sub>2</sub>, d<sub>yz</sub>) → (e, d<sub>z<sup>2</sup></sub>)</td><td></td></tr></tbody></table>

Electronic states up to the bright ¹Qₓᵧ state are shown. The excitation energies (in eV) are computed at the CASSCF(10,9)/CASPT2/ANO-RCC-VDZP level of theory at the S₀ minimum energy structure obtained at the B3LYP/LANL2DZ level of theory. In addition, the state nature, the character of the dominant transitions and the symmetry label is shown
experimentally these bands are found at 2.14–2.30 eV³⁰. These
difference are not unexpected, since on the one hand we do
not take into account the electrostatic environment of the
protein, and on the other hand the statistical average of
several protein conformations. Quantum Mechanics/Molecular
Mechanics (QM/MM) calculations decrease the overall spectrum
by around 0.84 eV for the lowest part of the spectrum. More
importantly, the energetic order and the gap between the states
in our model as compared with QM/MM is essentially conserved
(Supplementary Table 6). This is essential to guarantee a realistic
dynamical treatment.

The lowest valence states (up to ca. 2 eV) correspond to dark doubly-degenerate metal-centered (MC) states, involving transitions from the occupied $d_{xy}$, $d_{xz}$ and $d_{yz}$ orbitals to the unoccupied $d_z^2$ and $d_{x^2-y^2}$ orbitals. The crystal field splitting of the iron orbitals corresponds to a square pyramidal complex. The MLCT band is formed by states of mainly MLCT character, although some MC and LMCT character is also observed (Supplementary Table 7 for a detailed analysis of the state character). These bands of states have strong multi-configurational character, and were not correctly represented in previous single-reference studies based on density-functional theory¹⁵,¹⁸.

**Quantum dynamics.** Quantum wavepacket dynamics has been performed using a vibronic model containing the main vibronic coordinates (heme doming, symmetry-breaking, rotational and dissociative modes, see Supplementary Fig. 2 for a plot of the vibrations and Methods section for further details.) The evolution of diabatic state populations during the first 500 fs is shown in Fig. 1a. The results have been obtained by sampling 10 independent quantum dynamics simulations with different initial conditions. The initial conditions are obtained by projecting the heme–CO geometries extracted from molecular dynamics snapshots of a myoglobin protein onto our vibronic model Hamiltonian (Supplementary Table 8.) In the Figure, only the states which are mainly populated are shown (Supplementary Fig. 3 for the dynamics of all states.) Following the dynamics of each spin manifold, a sequential transfer is found to occur, first from singlet to triplet ($S \rightarrow T$) and then from triplet to quintet ($T \rightarrow Q$). The time constant for $S \rightarrow T$ transfer is estimated as 76 ± 15 fs, while the $T \rightarrow Q$ decay time is 429 ± 70 fs, in perfect agreement with the experimental rates of Cammarata et al. and Franzen et al.⁵,⁶,²⁰ (for the timescales for each initial condition, see Supplementary Table 9 and Supplementary Note 1.) In the singlet manifold, we observe an abrupt decay of the ¹Q population, which is mainly transferred to the singlet ¹MLCT manifold. The rate for this transfer is estimated to 26 ± 7 fs. The overall reaction is thus sequential ¹Q → ¹MLCT → ³MLCT → ⁵MLCT (see Fig. 1b). The initial step is a complete transfer from the ¹Q state to the ¹MLCT state, resulting in a negligible ¹Q population after ~100 fs. During the remaining transfer steps the three spin manifolds coexist. This is because the three manifolds are close in energy and strongly mixed through spin-orbit coupling. However, a clearly dominant spin state prevails during different intervals: singlet for $t < 75$ fs, triplet for $75 < t < 425$ fs and quintet at later times.

In Fig. 1c, we show the evolution of the Fe–C(O) and the Fe out-of-plane distances for the first 500 fs (see Supplementary Figs. 4–6 for the corresponding evolution per set of initial conditions). Clearly, we observe damped large amplitude stretching motions for the Fe–C distance, which converges to a quasi-stationary value of 2.2 Å. These motions are coherent, with a period of 42 fs. The Fourier transform of this signal shows a frequency of 800 cm⁻¹, much faster than the Fe–C(O) stretching frequency of 488 cm⁻¹. In 21 fs, the Fe–CO distance

oscillates between the equilibrium distance (1.7–1.8 Å) and 2.5 Å. These are the typical distances at which the ground-state dissociation occurs⁷. After the first oscillation, the wavepacket is essentially in the MLCT bands, which are repulsive for CO such that the CO does not recombine and undergoes continued oscillations. Relaxation of the structure damps these oscillations until an equilibrium value of 2.2 Å is reached. Both the equilibrium value and the coherent frequency might be disrupted by the presence of the protein, an effect not included in the present model. This would introduce a faster decoherence of the Fe and CO interactions than what we observe in our model. As for the evolution of the Fe out-of-plane motion, we observe oscillations of ±0.2 Å around the initial position, in good agreement with recent time-resolved X-ray crystallography data⁴. As expected, the iron center reacts to the elongation of the Fe–CO distance by following the CO. Very rapidly, it oscillates back to the porphyrin plane, where part of the energy is dissipated to histidine. In the second and subsequent oscillations, the oscillations of Fe and CO are out of phase, indicating that the bond is photolyzed. A large dispersion of the Fe–C distances, which is appearing after 0.2 ps, is due to the delocalization of the wavepacket over the ¹,³,⁵MLCT bands, indicating an unbound CO molecule. This dispersion is less marked for the Fe out-of-plane distance due to the more localized nature of the Fe atom, bounded strongly to the heme and the proximal histidine.

Experimental evidence shows that the CO is photolyzed within 70 fs⁵,⁶,²⁰. In order to determine whether dissociation occurs from the singlet or the triplet state, we have considered a reduced model with only the singlet manifold, in which we eliminate the coupling parameters to the triplet and quintet states (Supplementary Figs. 7–9). In this case, an ultrafast transfer from $Q \rightarrow ¹MLCT$ is still observed. The Fe–C(O) distance again exhibits oscillations with an amplitude of 0.7 Å. The fact that the amplitude is smaller than in the full model can be ascribed to the absence spin-orbit couplings which introduce strong mixing of the ¹MLCT surfaces with ³MLCT states, that extend the region of oscillations. From these results, we can infer that the transfer to the ¹MLCT manifold is sufficient to dissociate the CO, although the presence of ³MLCT favors further the photolysis. The dissociation is thus happening between 0.5 and 1.5 periods of the Fe–C(O) oscillation (20–60 fs), when the wavepacket is mainly in the ¹MLCT state.

After photolysis, the SCO mechanism brings the system sequentially to a HS state. After ~75 fs, the ³MLCT band is clearly dominant, although a residual ³MC participation is also observed (Supplementary Fig. 1). From the triplet manifold, a population transfer to the HS states is slowly building up, until it becomes dominant at around 400 fs. At this time, the wavepacket is distributed across the ⁵MLCT and ⁵MC excited states, staying trapped in these states. This is consistent with the slow time constant of ~3 ps observed in myoglobin, which has been attributed to a vibrational cooling due to the protein²⁰,³¹. Such excited state trapping has been observed recently in similar iron complexes²⁷. This trapping can be understood in terms of the potential energy surfaces depicted in Fig. 2. Upon excitation to the ¹Q band, a sequential ultrafast transfer, first to ¹MLCT and then to ³MLCT occurs. These transfers are energetically favored and therefore extremely fast. In the ³MLCT manifold, the wavepacket relaxes in a distribution of ³MLCT and ³MC states. The ³MC state can undergo direct transfer to the ⁵MLCT state, which is, however, energetically less favorable (and thus slower) than the singlet → triplet MLCT. The transition from the lowest ⁵MLCT to the ⁵MC state is energetically unfavoured. On the other hand, the ³MC states are energetically separated from the rest of the states, such that a transfer to the ⁵MC manifold becomes inefficient.
<figure><img src="image_2.png" /></figure>

**Fig. 1** Photodynamics of photolysis and spin-crossover. Quantum photodynamics of heme-CO complex during the first 0.5 ps, with initial conditions averaged over 10 molecular dynamics snapshots. **a** Evolution of diabatic populations for states ¹Q (magenta), ¹MLCT (orange), triplet band (green) and quintet band (blue). The ¹Q population rapidly decays giving rise to ¹MLCT population dominating by 75 fs, at which point the triplet population increases. The quintet population builds up more slowly, and evolves into the dominant state at around 350 fs. **b** Schematic representation of the reaction mechanism and interpretation in terms of time constants. Upon initial excitation to the Q-band, the metal-to-ligand charge transfer (MLCT) state is populated in ~25 fs. In a second step, the system relaxes to the triplet (~75 fs) and to the lowest quintet state (~430 fs). Black arrows indicate the direction of the electron transfer and the main nuclear motions. **c** Evolution of the Fe-C(O) distance (magenta, left axis) and the Fe out-of-plane distance (black, right axis). Large amplitude motions are observed with a period of oscillation of 40 fs. The amplitude of oscillation is initially 0.9 Å and converges towards a value of 2.2 Å. At this distance, the CO is essentially photolyzed. The standard deviation of these geometric values is shown as a shaded area. In the inset, the Fourier transform of the Fe-C(O) oscillations is shown (in cm⁻¹)

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

**Fig. 2** Potential energy surface along a dissociative mode. Plot of the potential energy surfaces along the Fe-C(O) distance. The LS state $S_0$, and the ¹Q states are shown in black and magenta, respectively. The singlet manifold is shown in orange, the triplet manifold in green and the quintet manifold in blue. The metal-centered (MC) states are shown explicitly, whereas the metal-to-ligand charge transfer (MLCT) states are depicted as a quasi-continuous band. Further, a schematic evolution of the wavepacket dynamics is shown. After absorption to the Q-band, the wavepacket undergoes large-amplitude oscillations in the Fe-CO coordinate on the ¹MLCT state. Then, it cascades down acquiring more ³MLCT and ⁵MLCT character as the spin crossover transitions occur. Finally, the wavepacket disperses, and the Fe-C(O) distance oscillations decrease to a value of 2.2 Å

**Vibronic effects.** The initial structure is close to a $C_{4v}$ symmetric geometry, in which the CO is in upright position and the porphyrin ring is close to the square planar conformation. The origin of the heme-CO bond is the typical $\sigma$-donation $\pi$-back-donation mechanism. A partial electron transfer from an occupied $\sigma$ orbital of CO to an empty $d_{z^2}$ orbital of Fe occurs simultaneously to a partial back transfer of electron density from the occupied $d_{yz}$ and $d_{xz}$ of Fe to the empty $\pi^+$ from $\sigma$-donated $\pi$-back-donation, respectively. At this geometry, the $\pi$ orbitals of porphyrin do not overlap with the d-orbitals.

Different types of vibrations can weaken the heme-CO bond, thus liberating the diatomic molecule. In Fig. 3, a schematic representation of the most representative vibrations is shown (for a detailed analysis of all vibrations, see Supplementary Table 10). The most important vibration is the dissociative coordinate (Fig. 3a), which corresponds to a symmetric Fe-CO stretching. This is the main reaction coordinate for CO photolysis. The Fe-CO stretching decreases the overlap between the iron and CO orbitals, thus decreasing both the strength of $\sigma$-donation and the $\pi$-back-donation. Once the CO is released, the remaining complex is then stabilized in the quintet state. This brings the Fe center out-of-plane, stabilized by $\pi$-interactions with the porphyrin antibonding orbitals.

Activation of the dissociative vibration is fundamental for releasing the CO from the heme. This vibration is of $a_1$ symmetry, as it is a vibration along the principal symmetry axis. However, the initial photon is absorbed exclusively in the porphyrin moiety.
<figure><img src="image_4.png" /><figcaption>**b** Dissociative</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**c** Rotational</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**d** Jahn-Teller</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>*π*-back donation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>*σ*-donation</figcaption></figure>

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

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

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

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

**Fig. 3** Vibrational motions and effect on the Fe-CO interaction. Representation of the nuclear motions (left) and the effect on the Fe-CO orbital interaction (right). The $\sigma_{CO}$ and $d_{xz}/d_{yz}$ orbitals (in blue) are occupied, while $e_{x/y}$, $\pi_{x/y/CO}^*$ and $d_{yz}$ orbitals (in red) are unoccupied. Blue arrows indicate the direction of the partial electron transfer, due to the overlap between an occupied and an unoccupied orbital (solid and dashed arrows indicate strong and weak overlap). Magenta arrows indicate the vibrational motion of the specified atoms. **a** In the initial structure, the Fe-CO bond is strongest due to partial electron transfer $d_{xz}/d_{yz} \rightarrow \pi_x^*/\pi_y^*$ ($\pi$-back donation) and $d_{yz}^2 \rightarrow \sigma_{CO}$ ($\sigma$-donation). **b** Dissociative vibration, which corresponds to a stretching of the Fe-CO bond. This structure destabilizes the Fe-CO bond and stabilizes by a partial electron transfer $d_{xz}/d_{yz} \rightarrow e_x^*/e_y^*$ when the Fe atom is out of plane. **c** Rotational vibrations, which correspond to the rotation of the CO moiety. This vibration decreases the overlap between Fe and CO orbitals. **d** Symmetry-breaking vibration, in which the porphyrin is distorted inducing Jahn-Teller and pseudo Jahn-Teller vibronic couplings in the excited state. This vibration has no effect on the Fe-CO bond, but activates the coupling between the porphyrin and the Fe densities

The porphyrin optical Q-state is E-symmetric, thus localized in the perpendicular plane with respect to the principal axis. This implies that the dissociative vibration is not activated in the quasi-C$_4$v initial geometry. The dissociation of CO thus requires an initial energy transfer from the porphyrin plane to the Fe center, which can only be possible by populating symmetry breaking vibrations. In the excited states, strong symmetry breaking of Jahn-Teller (JT) and pseudo Jahn-Teller (pJT) type occurs$^{32,33}$. The JT effect is the sudden increase of inter-state couplings between an E-degenerate state along symmetry breaking vibrational coordinates. These JT couplings induce a double-well minimum in the potential energy surfaces at elongated coordinates, while the energy is maximum at the initially symmetric structure. The pJT effect is analogous to the JT effect

but applies to non-degenerate states, usually requiring more than a single vibration and more than one electronic state. Both mechanisms are characterized by strong vibronic couplings. In the particular case of the heme-CO complex, the E-symmetry of the Q-band can be broken by $b_1$ and $b_2$ vibrations, in the so-called E⊗b JT mechanism$^{32}$. A representative vibration of this type is depicted in Fig. 3d. Such a vibration is not directly coupled to the Fe-CO bond, but is fundamental to induce the interstate coupling between the Q-band and the MLCT band. The relaxation energy induced by such vibrations can be as large as 0.5 eV, which rapidly leads the wavepacket far from the Franck-Condon region. Large amplitude motions have been experimentally observed in similar iron complexes$^{26}$.

The rotation of CO (see Fig. 3c) has been discussed at length in the literature, generating a controversy on whether it has an effect on the photochemistry$^{8,9,16,34}$. The rotation of CO slightly weakens the Fe-CO bond due to a less efficient overlap between the densities of each moiety. Furthermore, this rotation induces a Renner-Teller type of coupling (rovibronic coupling), in which e-symmetry bending modes of CO and $a_1$-type vibrations couple E and $A_1$ states, in what is known as the $(E \oplus A_1) \otimes (e \oplus A_1)$ effect. This mechanism can also couple the Q-band (E symmetry) and the MLCT band ($A_1$ symmetry) to the CO bending (e symmetry) and the dissociative iron-CO stretching ($A_1$ symmetry). The population evolution of the singlet manifold excluding the JT modes—and thus the rotational modes that are responsible of the Q-MLCT coupling—shows a similar population and Fe-C kinetics. This indicates that any symmetry breaking coordinate, related to protein fluctuations, rotations or symmetry breaking Jahn-Teller effects, will accelerate the coupling to the MLCT band.

The porphyrin doming ($A_1$ vibration of the porphyrin ring, see modes 8,12 in Supplementary Fig. 3) has been pointed out as an important motion for the photodynamics$^{5,6,35}$. The doming motion is out-of-plane, and has a slow period of 371 fs. It is of the same symmetry as the dissociative mode, to which it is strongly coupled. The potential energy surface cuts along the doming mode indicate a weak reorganization energy and a small vibronic coupling for all potential energy surfaces. Therefore, even though this motion is necessary to stabilize the high-spin state, it plays a passive role and does not significantly affect the short-time dynamics of photolysis (Supplementary Figs. 10–12). Therefore, we conclude that the doming motion is in fact a consequence of the CO dissociation, that is strongly coupled to the Fe-CO dissociation coordinate, and is activated during the second kinetic step of 280 fs$^{5,6,20}$.

**Spin crossover.** The photochemical spin crossover mechanism is an important reaction in many organometallic complexes, in which ligand chromophores allow a usually fast low-to-high spin transition of the complex. Recent theoretical and experimental studies for similar iron complexes have shown that ligands play a fundamental role in modulating the time constants and the electronic pathway during the photoreaction$^{4–6}$. The spin crossover mechanism is mediated by the relativistic spin-orbit coupling (SOC). Although this is a purely electronic transition, the vibrational effects on the SOC seem fundamental to explain the ultrafast nature of spin crossover in metallic complexes. In the particular case of heme-CO, the main SOC strengths are about 200–300 cm$^{-1}$, which are of similar size as the JT vibronic couplings. Therefore, internal conversion and intersystem crossing are competing mechanisms in the dynamical relaxation mechanism.

In the heme-CO complex, the initially populated Q-band of porphyrin exhibits a negligible SOC with any other state. In
general, SOC is strong when metal density is present in the electronic state. Therefore, the wavepacket is initially dominated by the symmetry breaking vibronic couplings in the singlet manifold, leading to the transition $^1Q \rightarrow ^1MLCT$ and $^1Q \rightarrow ^1MC$. Both the $^1MLCT$ and $^1MC$ have strong spin-orbit coupling strength with $^3MLCT$ and $^3MC$ states, respectively. In general, the SOC is strong when the selection rules $\Delta S = 1$ and $\Delta L = 1$ are obeyed. The latter condition implies that one of the d orbitals changes orientation between the triplet and the singlet state. For example, a strong coupling occurs between the singlet state of character $d_{xy} \rightarrow d_{x^2-y^2}$ with the triplet state of character $d_{xz} \rightarrow d_{x^2-y^2}$. We observe that singlet-triplet transitions occur at a faster rate in the MLCT band than within the MC states (Supplementary Fig. 4). At around 50 fs, the $^3MLCT$ population becomes dominant, at the expense of the singlet manifold populations. This is because there are a larger number of crossing points between MLCT states of different multiplicity than between MC states (see Fig. 2 and Supplementary Fig. 4). This results in a dominant population of $^3MLCT$ and $^5MLCT$ states during the spin-crossover transitions. The $^3Q$ band triplets, which were pointed out as playing a role in the photophysics of heme–CO, contribute negligibly to the dynamics$^{20}$. The $^1MC$ state exhibits a residual population after 100 fs, following prevalent transfer to the $^3MC$ state, with a population which is two times smaller than the $^3MLCT$ population. Finally, the $^3MLCT \rightarrow ^5MLCT$ transfer builds up population in the quintet band. A $^3MC \rightarrow ^5MC$ transfer, however, is not observed, due to the fact that the $^3MC$ and $^5MC$ states do not cross in an accessible region. Our results show rather that the wavepacket remains trapped in the $^5MLCT$ band. The trapping occurs in regions of the energy surface where relaxation to lower states (in this case MC states) is slow due to large gaps and small vibronic couplings, whereas a band of electronic states acts as an energy dissipation force, thus localizing the wavepacket in an excited state. Such observation for heme–CO is in line with the valence tautomerism model of Franzen, Martin et al.$^{20}$ In summary, the existence of quasi-continuous MLCT bands degenerate to the $^1Q$ state and strong vibronic couplings are at the origin of the ultrafast photolysis and spin-crossover in the heme–CO complex.

## Discussion

The mechanism of heme–CO photolysis and spin crossover has been elucidated by means of high-dimensional quantum dynamical calculations. The results indicate that the vibronic mechanism can be summarized in terms of a sequential $^1Q \rightarrow ^1MLCT/MC \rightarrow ^3MLCT \rightarrow ^5MLCT$ transfer. During the first step, which is completed within 20–60 fs, the CO is dissociated, prior to the spin-crossover mechanism. The dissociation occurs in the singlet MLCT excited states. For these states, the $\sigma - \pi$ bond between Fe and CO bond is weakened. Subsequently, the wavepacket relaxes stepwise to the quintet band under the effect of spin-orbit coupling. We assign the three experimental rates observed for the heme–CO complex in myoglobin, recently reported in Ref. 20 and Refs. 5,6 as follows: (i) the 50–70 fs time constant corresponds to the singlet-triplet spin crossover, in which the population is mainly in the $^3MLCT$ excited state; (ii) the 300–400 fs time constant corresponds to the triplet-quintet spin crossover, in which the population is mainly in the $^5MLCT$ excited state; (iii) the 2.4 ps time constant we tentatively assign to the relaxation of the $^5MLCT$ excited state to the final HS state of the heme.

## Methods

**Electronic structure.** Unconstrained structure optimization and normal-mode calculations of $S_0$ in gas phase have been performed with density-functional theory (DFT) using the B3LYP exchange-correlation functional and the LAN2DZ basis

set$^{36–39}$. These computations have been performed with Gaussian09$^{40}$. Excited states have been computed with CASSCF(10,9)/CASPT2$^{41,42}$. The active orbitals are the Fe(II) 3d-orbitals and the Gouterman’s porphyrin $\pi$-orbitals (see Supplementary Fig. 2 and Supplementary Note 2). State-averaged calculations have been performed for 35 singlet states, 50 triplet states, and 30 quintet states. In all calculations, we employed relativistic Douglas-Kroll-Hess (DKH) Hamiltonian with the ANO-RCC-VDZ basis set$^{43–45}$. Larger basis sets do not significantly affect the energetics (Supplementary Table 2). Multi-state CASPT2 has been performed with frozen core and 116 frozen virtual orbitals. Similarly to coupled-cluster$^{46}$, deleting virtual orbitals accelerates computations and keeps the essential physics (Supplementary Table 2). In particular, the spin crossover transition is well represented and the DFT minima are also CASPT2 minima. Spin-orbit couplings have been obtained by expanding the DKH Hamiltonian in the basis of 20 singlets, 20 triplets and 20 quintets$^{45}$. All multi-configuration results have been obtained with Molcas$^{47}$.

**Vibronic Hamiltonian.** The model Hamiltonian for representing the dissociation and the spin crossover is of the form:

$$ \hat{H}(\mathbf{q}) = \sum_{S=0}^{2} \hat{H}_{S}(\mathbf{q}) + \sum_{S,S'=1}^{2} \hat{V}_{SS'}^{SO}, \quad (1) $$

in which $S=0, 1, 2$ is the total spin quantum number, $\mathbf{q}$ is the mass-frequency-weighted vibrational coordinate vector, $\hat{V}_{SS'}^{SO}$ is the spin-orbit coupling. The vibronic Hamiltonian, represented in the basis of ground-state normal modes, is given by

$$ \hat{H}_S(\mathbf{q}) = \hat{T}_S(\mathbf{q}) + \hat{V}_S(\mathbf{q}), \quad (2) $$

containing the harmonic kinetic energy and potential operator. The potential is given by

$$ \hat{V}_S(\mathbf{q}) = \sum_{m_s = -|S|}^{|S|} \sum_{I,J=1}^{20} [\delta_{IJ} E_I^S + V_{IJ}^S(\mathbf{q})] |^S \Psi_I^{m_s} \rangle \langle ^S \Psi_J^{m_s}| + \text{c.c.}, \quad (3) $$

in which $m_s$ is the azimuthal spin and $|^S \Psi_I^{m_s} \rangle$ is the state wavefunction. The $E_I^S$ energies are diabatic at the $S_0$ minimum geometry; at the quasi-C$_{4v}$ geometry these can be considered quasi-diabatic. The $V_{IJ}^{PT2S}(\mathbf{q})$ potential is derived from one-dimensional potential energy surface cuts which were diabatized by fitting to the following functional form,

$$ V_{II}^{S}(\mathbf{q}) = \sum_{I=1}^{5} \frac{1}{\hbar^2} \frac{\partial^2 E_I^S(\mathbf{q})}{\partial q_I^2} |_{\mathbf{q}=0} \mathbf{q}' + \alpha_I^S e^{\beta_I^S} (\mathbf{q} - \mathbf{y}_I^S) + \mathbf{D}_I^S (1 - e^{\beta_I^S} (\mathbf{q} - \mathbf{q}_I^S))^2, \quad (4) $$

$$ V_{II}^{S}(\mathbf{q}) = v_{II}^{S}(\mathbf{q})\mathbf{q}. \quad (5) $$

Here, the inter-state couplings $V_{II}^{S}(\mathbf{q})$ are approximated within a linear vibronic model. The set of parameters are determined by a least square fitting

$$ f(\mathbf{y}^S) = \frac{1}{N} \sum_{I\sigma} (E_I^{\text{mod},S}(\{\mathbf{y}_I^S\}, \mathbf{q}_\sigma) - E_I^{\text{PT2,S}}(\mathbf{q}_\sigma))^2, \quad (6) $$

where $\mathbf{y}^S = \{\frac{\partial^2 E_I^S}{\partial q_I^2}, \alpha_I^S, \beta_I^S, \mathbf{y}_I^S, \epsilon_I^S, \mathbf{q}_I^S, \mathbf{v}_{II}^S\}$, $\mathbf{q}_\sigma$ is the discretized value of a normal mode coordinate and $E_I^{\text{mod},S}(\mathbf{q}_\sigma)$ and $E_I^{\text{PT2,S}}(\mathbf{q}_\sigma)$ are the spin-free adiabatic energies from the model Hamiltonian and the CASPT2 calculation, respectively. The latter are obtained by one-dimensional potential energy surface cuts along each ground-state normal mode direction. The passive normal modes have been kept frozen along the cuts.

The fitting of quasi-continuous bands requires an accurate initial guess. We employ the following protocol: (i) quasi-diabatic surfaces are obtained by maximum overlap of CI coefficients; (ii) pre-fitting is done on the quasi-diabatic surfaces by employing first and second order Taylor expansion or the exponential function, (iii) the adiabatic PT2 surfaces are fitted with the full diagonal Hamiltonian, and (iv) the full Hamiltonian is fitted until the total error is <0.10 eV.

**Vibrational modes.** Vibrational coordinates have been chosen by performing one-dimensional potential energy cuts and selecting the modes with largest reorganization energy. There are a total of 15 explicit vibrational modes in the model: 1 symmetric and 1 asymmetric Fe-CO stretchings, 2 doming modes, 1 CO stretching, 1 in-plane, and 1 out-of-plane porphyrin vibrations, and 4 doubly degenerate modes (Jahn-Teller modes of porphyrin out-of-plane, a porphyrin in-plane mode and a pure CO bending.) These modes have been selected in order to reproduce the symmetry breaking that allows the $Q \rightarrow MLCT$ ultrafast transition, the dissociation and the relaxation of the quintet heme. A schematic plot of these vibrations can be found in Supplementary Fig. 3. Duschinsky rotation effects and
anharmonic couplings between vibrational modes have been neglected in our model; in particular, the latter are negligibly small compared to the vibronic couplings (see Supplementary Table 11 and Supplementary Note 3).

**Quantum dynamics.** The Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) method<sup>48,49</sup> (Heidelberg package, version 8.5.5) was used to perform wavepacket propagation up to 1 picosecond for a vibronic coupling model comprising 179 electronic states of singlet, triplet and quintet multiplicities, and 15 vibrational modes exhibiting strong vibronic coupling. A four-layer representation of the wavefunction was employed, where 8 single-particle functions (spf) were chosen per layer, with a primitive harmonic-oscillator discrete variable representation (DVR) comprising up to 80 DVR points (see Supplementary Note 4 and Supplementary Fig. 13). Initial conditions were constructed in accordance with the relative oscillator strengths of the singlet ligand and metal-centered states at the Franck-Condon geometry.

**Data availability**

The authors declare that the data supplementary the findings of this study are available within the paper and its supplementary information files. Parameters for the model Hamiltonian and potential energy surfaces along each vibrational mode are available from the corresponding authors upon request.

Received: 8 October 2017 Accepted: 14 September 2018

Published online: 29 October 2018

**References**

1. Messerschmidt, A., Huber, R., Poulos, T. & Wieghardt, K. (eds.) *Handbooks of Metalloproteins*, Vol. 1 (John Wiley and Sons, 2001).

2. Chung, L. W., Li, X. & Morokuma, K. Modeling Enzymatic Reactions in Metalloenzymes and Photobiology by Quantum Mechanics (QM) and Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations. Chap. 3. *Quantum Biochemistry* (Wiley-VCH Verlag GmbH, 2010).

3. Ansari, A. et al. Protein states and proteinquakes. *Proc. Natl Acad. Sci. USA* **82**, 5000–5004 (1985).

4. Barends, T. R. et al. Direct observation of ultrafast collective motions in CO myoglobin upon ligand dissociation. *Science* **350**, 445–450 (2015).

5. Levantino, M., Schiro, G., Lemke, H. T. & Cottone, G. Ultrafast myoglobin structural dynamics observed with an X-ray free-electron laser. *Nat. Commun.* **6**, 6772 (2015).

6. Levantino, M. et al. Observing heme doming in myoglobin with femtosecond X-ray absorption spectroscopy. *Struc. Dynam.* **2**, 041713 (2015).

7. Harvey, J. N. DFT computation of the intrinsic barrier to CO geminate recombination with heme compounds. *J. Am. Chem. Soc.* **122**, 12401–12402 (2000).

8. Spiro, T. G. & Kozlowski, P. M. Is the CO adduct of myoglobin bent, and does it matter? *Acc. Chem. Res.* **34**, 137–144 (2001).

9. Rovira, C., Schulze, B., Eichinger, M., Evanseck, J. D. & Parrinello, M. Influence of the heme pocket conformation on the structure and vibrations of the Fe–CO bond in myoglobin: a QM/MM density functional study. *Biophys. J.* **81**, 435–445 (2001).

10. Harvey, J. N. Spin-forbidden CO Ligand recombination in myoglobin. *Farad. Discuss.* **127**, 165–177 (2004).

11. Meier, C. Laser control of vibrational excitation in carboxyhemoglobin: a quantum wavepacket study. *J. Chem. Phys.* **123**, 044504 (2005).

12. Debnath, A., Falvo, C. & Meier, C. State-selective excitation of the CO stretch in carboxyhemoglobin by mid-IR laser pulse shaping: a theoretical investigation. *J. Phys. Chem. A* **117**, 12884–12888 (2013).

13. Kühn, O. Multidimensional vibrational quantum dynamics of CO-heme compounds: ultrafast IVR mediated Fe-CO bond-breaking after CO excitation? *Chem. Phys. Lett.* **402**, 48–53 (2005).

14. Waleh, A. & Loew, G. H. Quantum mechanical studies of the photodissociation of carbonylheme complexes. *J. Am. Chem. Soc.* **104**, 2346–2351 (1982).

15. Dreuw, A., Dunletz, B. D. & Head-Gordon, M. Characterization of the relevant excited states in the photodissociation of CO-ligated hemoglobin and myoglobin. *J. Am. Chem. Soc.* **124**, 12070–12071 (2002).

16. De Angelis, F., Car, R. & Spiro, T. G. Role of ligand bending in the photodissociation of O₂ vs CO-heme: a time-dependent density functional study. *J. Am. Chem. Soc.* **125**, 15710–15711 (2003).

17. Ohta, T., Pal, B. & Kitagawa, T. Excited state property of hardly photodissociable heme-CO adduct studied by time-dependent density functional theory. *J. Phys. Chem. B* **109**, 21110–21117 (2005).

18. Dunietz, B. D., Dreuw, A. & Head-Gordon, M. Initial steps of the photodissociation of the CO ligated heme group. *J. Phys. Chem. B* **107**, 5623–5629 (2003).

19. Halden, J. & Lorrain, S. J. The oxygen tension of arterial blood. *J. Physiol.* **20**, 497–520 (1896).

20. Franzen, S., Kiger, L., Poyart, C. & Martin, J. Heme photolysis occurs by ultrafast excited state metal-to-ring charge transfer. *Biophys. J.* **80**, 2372 (2001).

21. Petrich, J. W., Poyart, C. & Martin, J. L. Photophysics and reactivity of heme proteins: a femtosecond absorption study of hemoglobin, myoglobin, and protoheme. *Biochem. J.* **27**, 4049–4060 (1988).

22. Adar, F., Gouterman, M. & Aronowitz, S. Fluorescence, resonance Raman, and raditionless decay in several hemoproteins. *J. Phys. Chem.* **80**, 2184–2191 (1976).

23. Sousa, C. et al. Ultrafast deactivation mechanism of the excited singlet in the light-induced spin crossover of $[Fe(2,2'-bipyridine)₃]^{2+}$. *Chem. Eur. J.* **19**, 17541–17551 (2013).

24. Eng, J., Gourlaouen, C., Gindensperger, E. & Daniel, C. Spin-vibronic quantum dynamics for ultrafast excited-state processes. *Acc. Chem. Res.* **48**, 809–817 (2015).

25. Pápai, M., Penfold, T. J. & Moller, K. B. Effect of tert-butyl functionalization on the photoexcited decay of a Fe(II)-N-heterocyclic carbene complex. *J. Phys. Chem. C.* **120**, 17234–17241 (2016).

26. Auböck, G. & Chergui, M. Sub-50-fs photoinduced spin crossover in $[Fe(bpy)₃]^{2+}$. *Nat. Chem.* **7**, 629–633 (2015).

27. Lemke, H. T. et al. Coherent structural trapping through wave packet dispersion during photoinduced spin state switching. *Nat. Comm.* **6**, 6772 (2017).

28. Engel, N. et al. Light-induced relaxation dynamics of the ferricyanide ion revisited by ultrafast XUV photoelectron spectroscopy. *Phys. Chem. Chem. Phys.* **19**, 14248–14255 (2017).

29. Manthe, U. Wavepacket dynamics and the multi-configurational time-dependent Hartree approach. *J. Phys. Condens. Matter* **29**, 253001 (2017).

30. Bowen, W. J. The absorption spectra and extinction coefficients of myoglobin. *J. Biol. Chem.* **179**, 235–245 (1949).

31. Mizutani, Y. & Kitagawa, T. Direct observation of cooling of heme upon photodissociation of carbonmonoxy myoglobin. *Science* **278**, 443–446 (1997).

32. Bersuker, I. *The Jahn-Teller Effect and Vibronic Interactions in Modern Chemistry*. Modern Inorganic Chemistry (Springer Verlag, 1984).

33. Bersuker, I. Pseudo-Jahn-Teller effect-a two-state paradigm in formation, deformation, and transformation of molecular systems and solids. *Chem. Rev.* **113**, 1351–1390 (2013).

34. Carme, R. & Parrinello, M. Harmonic and anharmonic dynamics of Fe–CO and Fe–O₂ in heme models. *Biophys. J.* **78**, 93–100 (2000).

35. Klug, D. D. et al. Doming modes and dynamics of model heme compounds. *Proc. Natl Acad. Sci. USA* **99**, 12526–12530 (2002).

36. Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. *Phys. Rev. A.* **38**, 3098 (1988).

37. Lee, C., Yang, W. & Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. *Phys. Rev. B* **37**, 785 (1988).

38. Vosko, S. H., Wilk, L. & Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis. *Can. J. Phys.* **58**, 1200–1211 (1980).

39. Hay, P. J. & Wadt, W. R. Ab initio effective core potentials for molecular calculations. potentials for K to Au including the outermost core orbitals. *J. Chem. Phys.* **82**, 299 (1985).

40. Frisch, M. J. et al. Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford CT (2009).

41. Roos, B. J. The complete active space self-consistent field method and its applications in electronic structure calculations., Vol. 69 In *Adv. Chem. Phys.* 399 (John Wiley and Sons, Ltd., 1987).

42. Andersson, K., Malmqvist, P. & Roos, B. O. Second-order perturbation theory with a complete active space self-consistent field reference function. *J. Chem. Phys.* **96**, 12181–1226 (1992).

43. Roos, B. O., Lindh, R., Malmqvist, P. A., Veryazov, V. & Widmark, P. O. Main group atoms and dimers studied with a new relativistic ANO basis set. *J. Phys. Chem. A* **108**, 2851 (2005).

44. Roos, B. O., Lindh, R., Malmqvist, P. A., Veryazov, V. & Widmark, P. O. New relativistic ANO basis sets for transition metal atoms. *J. Phys. Chem. A* **109**, 6575–6579 (2005).

45. Roos, B. O. & Malmqvist, P. A. Relativistic quantum chemistry: the multiconfigurational approach. *Phys. Chem. Chem. Phys.* **6**, 2919–2927 (2004).

46. Send, R., Kaila, V. R. I. & Sundholm, D. Reduction of the virtual space for coupled-cluster excitation energies of large molecules and embedded systems. *J. Phys. Chem.* **134**, 214114 (2011).
47. Aquilante, F. et al. MOLCAS 8: new capabilities for multiconfigurational quantum chemical calculations across the periodic table. *J. Comput. Chem.* **37**, 506–541 (2016).

48. Meyer, H.-D., Manthe, U. & Cederbaum, L. S. The multi-configurational time-dependent Hartree approach. *Chem. Phys. Lett.* **165**, 73–78 (1990).

49. Wang, H. & Thoss, M. Multilayer formulation of the multiconfiguration time-dependent Hartree theory. *J. Chem. Phys.* **119**, 1289 (2003).

## Acknowledgements

MHR acknowledges financial support from Agence Nationale de la Recherche (grant ANR-16-CE29-0008, BIOMAGNET) and the Alexander von Humboldt-Foundation. K. F. and I.B. gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft via RTG 1986 “Complex Scenarios of Light Control”. I.B. thanks Prof. Gerhard Hummer (MPI Biophysics, Frankfurt) for making MD snapshots available for the present study. All authors thank Prof. Nicolas Ferré, Dr. Padmabati Mondal, Dr. Anthony Kermagoret and the reviewers for useful comments. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program “Investissements d'Avenir” supervised by the Agence Nationale de la Recherche.

## Author contributions

M.H.R. designed the research in collaboration with I.B. M.H.R., I.B., and H.T. developed the vibronic Hamiltonian. M.H.R. and K.F. performed the electronic structure calculations. M.H.R., K.F., and I.B. performed the quantum dynamics. M.H.R. wrote the paper. All authors contributed to the discussions and gave comments on the manuscript.

## Additional information

**Supplementary Information** accompanies this paper at https://doi.org/10.1038/s41467-018-06615-1.

**Competing interests:** The authors declare no competing interests.

**Reprints and permission** information is available online at http://npg.nature.com/reprintsandpermissions/

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018