Dioxygen Activation by Iron Complexes: The Catalytic Role of

Jan 25, 2018 - Enzymes containing heme, nonheme iron, or copper active sites play an essential role in the dioxygen binding and activation for substra...
2 downloads 4 Views 852KB Size
Subscriber access provided by Thompson Rivers University | Library

Article

Dioxygen Activation by Iron Complexes: The Catalytic Role of Intersystem Crossing Dynamics for a Heme-Related Model Likai Du, Fang Liu, Yanwei Li, Zhongyue Yang, Qingzhu Zhang, Chaoyuan Zhu, and Jun Gao J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b11462 • Publication Date (Web): 25 Jan 2018 Downloaded from http://pubs.acs.org on January 30, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dioxygen Activation by Iron Complexes: The Catalytic Role of Intersystem Crossing Dynamics for a Heme-Related Model

Likai Du*1, Fang Liu1, Yanwei Li2, Zhongyue Yang3, Qingzhu Zhang2, Chaoyuan Zhu4, Jun Gao*1

1

Hubei Key Laboratory of Agricultural Bioinformatics, College of Informatics, Huazhong

Agricultural University, Wuhan, 430070, P.R. China. 2

Environment Research Institute, Shandong University, Jinan, 250100, P. R. China

3

Department of Chemistry and Biochemistry, University of California, Los Angeles, California

90095, United States 4

Institute of Molecular Science, Department of Applied Chemistry and Center for Interdisciplinary

Molecular Science, National Chiao Tung University, Hsinchu 30010, Taiwan *To whom correspondence should be addressed. Likai Du: [email protected]; Jun Gao: gaojun@ mail.hzau.edu.cn

Abstract Enzymes containing heme, non-heme iron or copper active sites play an essential role in the dioxygen binding and activation for substrate oxidation. The conceptual challenges to the quantitative modeling of this primary catalytic step arise from (1) instrinsically electronic non-adiabaticity of the spin flip events of the triplet dioxygen molecule (3O2), mediated by spin-orbit coupling; (2) possible heat dissipation channels, due to the high exothermicity of dioxygen binding processes. Herein, the spin-forbidden dioxygen binding dynamics of a reduced heme model was directly investigated in terms of the non-adiabatic trajectory surface hopping dynamics, involving the coupled singlet, triplet and quintet states. This work reveals the complexity of this elemental reaction, and the binding/dissociation dynamics of iron peroxo species is important to interpret the subsequent H-atom abstraction reaction step. Furthermore, we identify non-adiabatic dynamical effects that could not be observed through traditional calculations of static geometries.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Heme enzymes, such as cytochrome P4501, catalase-peroxidase2, nitric oxide synthase (NOS)3 and heme oxygenase (HO)4 are an important class of proteins that catalyze a variety of essential oxidative reactions, including the biosynthesis of steroids, the degradation of xenobiotics, the metabolism of drugs and so on. These reactions involve O2 binding and activation as the first elemental step at the active sites, which generates iron peroxide complexes that subsequently undergo substrate oxidation, while O2 is reduced to H2O2 or H2O via two-electron or four-electron mechanisms respectively.5-9 The generated iron peroxide complexes are the precursors to other families of iron oxidants, such as the iron hydroperoxo complex (Cpd 0) or iron-oxo complex (Cpd I). The O2 binding processes have been of considerable interest to chemistry community. Both “end-on” (one atom bound to the metal, known as η1) and “side-on” (two atoms bound to the metal, known as η2) are proposed as possible geometries for the binding complex.10-14 End-on O2 complexes always have a bent geometry at the proximal oxygen atom, and myoglobin and hemoglobin are famous examples,15-18 while side-on O2 binding complex is also found in a number of non-heme enzymes19-20, i.e. naphthalene dioxygenase. The influence of O2 binding on the reactivity has been much less explored. It has been known that during oxidations, the molecular system is required to “hop” from the PES corresponding to the initial spin state onto that corresponding to the one of the product state. Two state reactivity (TSR) or multistate reactivity (MSR) models have been developed to better understand how different electronic states contribute to the catalysis.21-28 Since several low-lying spin PESs should be considered on the same foot21,

23, 25, 27-28

, the most possible reaction pathway is usually

determined by connecting the transition states and intermediates with lower energies (Scheme 1).

ACS Paragon Plus Environment

Page 2 of 26

Page 3 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Scheme 1. (a) An overview of dioxygen binding for heme enzymes. The O2 binding potential energy surface (PES) with heme group (top left), and the illustration of two state reactivity (TSR, bottom right) are given. (b) The molecule orbital diagram for the free ground state dioxygen molecule (3O2).

The spin forbidden problem is conceptually resolved by TSR, but the non-adiabatic dynamics detail of how the reactions actually occur, from binding to oxidation, still remains unclear. The reversibly binding/dissociation process of the free oxygen molecule (3O2) at the iron porphyrin site is often said to be non-adiabatic. On the basis of optimized molecular geometries, one could apply the Laudau-Zerner theory to estimate the hopping probabilities between different spin states, and thus the rate coefficients.29 Experimentally, the timescale of spin state change or intersystem crossing (ISC) event of transition-metal complexes ranges from sub-ps to ns-µs time scale.30-36 Non-statistical (dynamic) effects are involved when the spin state change or ISC occurs in a ultrafast time scale, in which intramolecular vibrational energy can not be fully redistributed. The O2 binding at the iron porphyrin site is highly exothermic (typically ~1.0 eV), making the relaxation dynamics a non-trivial problem in understanding the complexation. This, coupled with electronic non-adiabaticity, necessitates the non-adiabatic dynamics simulations, which inform how dioxygen binding and spin state flip happen in a time-resolved manner and at single-molecule level. It also reveals critical features of the spin forbidden processes, and how the proposed TSR or MSR mechanism works at the full atomic resolution. It should be pointed out that the spin state changes as well as other non-adiabatic processes are

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

common in photochemistry. Several approaches, from quantum dynamics in reduced dimensions to semiclassical dynamics have been applied to the excited states non-adiabatic dynamics involving the spin flip events.37-44 The on-the-fly trajectory surface hopping methodology has been successfully extended to the photoinduced intersystem crossing for some organic molecules, such as nucleobases, endoperoxides, acrolein, o-nitrophenol and so on.41, 45-48 Therefore, it is possible to address the non-adiabatic dynamic feature of the spin flip events of the iron-activated O2 species, which are capable of oxidizing a broad range of substrates. Here, we performed direct non-adiabatic dynamics simulations on the oxygen binding to a reduced porphyrin model with on-the-fly trajectory surface hopping methodology. This is also our first attempt to integrate non-adiabatic dynamics to understand catalytic problems. This work mainly focuses on two important issues. 1) For the thermal activated dioxygen binding reaction, can we observe the ultrafast spin flip events in real time by the non-adiabatic dynamics? 2) If the non-adiabatic dynamics works, how about the binding mode of dioxygen on the iron center under the thermal perturbation and spin state changes? How about the possible timescale of the spin crossover events in dioxygen binding process? Thus, we use a very simplified porphyrin model to evaluate the ability of the on-the-fly surface hopping method, which has been carefully validated in a few related works.49-51 The dissociation of dioxygen molecule from the iron center is often observed during the non-adiabatic surface hopping dynamics, beyond the optimized side-on and end-on mode. The timescale of spin flip events can be comparable to the lifetime of transition states (~a few hundred femtoseconds)52-53. The non-adiabatic dynamics simulation with coupled PESs reveals some critical features of the spin forbidden O2 binding process, and how the proposed TSR or MSR mechanism works at the full atomic resolution. We hope this work enriches our understanding of spin crossover dynamics in thermal activated chemical reactions of heme related compounds.

2. Methods and Computation Details As an initial evaluation of the non-adiabatic dynamics on the thermal activated catalytic problems, we decide to consider the active sites of heme enzymes-simplified porphyrin related models. One is the full model FeP(Im), where P stands for porphin and Im for imidazole; the other is its simplified mimic (FeL2NH3), in which the porphin ring is replaced by two bidentate,

ACS Paragon Plus Environment

Page 4 of 26

Page 5 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

N-donor ligands C3N2H5- (abbrev. L), and the axial Im are replaced by NH3. These heme related model complexes have been extensively studied by a series of benchmark works51, 54-57. The structures of the iron porphin model and its simplified mimics are shown in Figure 1.

Figure 1. The reduced porphyrin model (iron porphin) and its simplified mimics are shown.

2.1 Quantum Chemistry Calculations The geometries of iron complexes and their dioxygen binding counterparts were optimized at B3LYP-D3/6-31G* level. Although the B3LYP-D3 is not necessarily the best functional, its performance is generally not the worst in our experience and others50, 58. Further discussions on this issue can be found in the “Implementation Details” section. The used basis set is also compared with larger basis sets (Table S1). And only the lowest singlet, triplet and quintet states of the iron complexes are considered, which are commonly studied in the reaction mechanism studies of heme and related iron complexes.50-51, 59-62 The minimal energy crossing points (MECPs) between different spin states was identified by the gradient following method.63 The quantum chemistry calculations were performed with Gaussian 09 package64. 2.2 The Intersystem Crossing Dynamics The thermal activated spin-forbidden reactions was investigated with the on-the-fly fewest switch surface hopping approach. In the surface hopping method, the nuclear and electronic degrees of freedoms are treated by classical and quantum dynamics, respectively. This approach is possible to treat relatively large and realistic molecular systems with full degree of freedom. The

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 26

original fewest switch surface hopping method65-66 was usually used for describing radiationless transitions between same-spin states for photo-induced processes. In order to treat the intersystem crossing, the electronic Hamiltonian must include the spin-orbit interaction39,

41, 67-72

, which

promotes radiationless transitions between different spin states. Since the SOC value for the iron complexes is relatively small (tens to one hundred of cm-1), we prefer to use a spin-diabatic approach to incorporate the SOC effects45, 71, and each spin multiplet is treated as a single electronic state. Similar treatment has been commonly used to investigate the inter-system crossing photo-induced processes41,

48

. Thus, we only provide a very brief

introduction of the surface hopping method for thermal activated multiple spin states reactivity reactions. The propagation of the quantum amplitudes along the nuclear trajectory R(t) can be given.

dcβ (t ) dt

[

= −ih −1 ∑ cα (t ) H eff , βα − ihv(t ) ⋅ d βα

]

α

(1)

= −ih −1 ∑ cα (t )[H 0 + H so ]βα − ∑ cα (t )v (t ) ⋅ d βα α

α

In Eq. 1, Heff is the effective electronic Hamiltonian, for which the spin-orbit operator is added to the zero-order electronic Hamiltonian (H0) as a perturbation. And v(t) is the nuclear velocities, while dβα is the derivative coupling vector. The spin-orbit matrix element (Hso) is the coupling between electronic state α and β. Only the transitions between different spin states are required for thermal activated reaction involving iron complexes. The derivative coupling contribution would vanish, and Eq. 1 becomes

dcβ (t ) dt

= −ih −1 ∑ cα (t )[H 0 + H so ]βα

(2)

α

The density matrix representation of quantum amplitudes was adopted in our implemented code. The fewest-switches criterion is used to account for the hops caused by spin-orbit coupling, similar as the original one given by Tully.65-66 And the transition probability from state β to state α is given as follows, if only the spin state transitions are considered.

Pαβ dt =

[

(

*

− 2 h −1 Im cα cβ H so ,αβ cα (t )

2

)]dt

(3)

The nuclear degrees of freedom are propagated by the independent classical trajectories R(t) on the currently-occupied electronic spin state, which is computed by numerical integration of

ACS Paragon Plus Environment

Page 7 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Newton’s equations. 2.3 Initial Conditions and Dynamics Simulation Details An ensemble of initial conditions (i.e. coordinates and momenta) needs to be prepared. Using the corresponding harmonic frequencies, the initial coordinates and momenta were generated for the O2 binding complex (S0 minima), from the Wigner distribution function of the first vibrational level for the ground singlet electronic state. Then, the surface hopping dynamics were performed under a few distinct perturbations. (1) Starting from the initial geometries, we performed the surface hopping dynamics at the open shell singlet state. And the initial velocities were included, plus a Boltzmann sampling of thermal energy available at 300 K with a random phase. (2) After the dioxygen molecule binding to the iron site, the high exothermicity in the adsorption process may simply result in excessive kinetic energy. This could cause a highly-distorted molecular geometry or even cause the dissociation of the system, if the heat dissipation channels are missing. To take into account this “hot” molecular motion, the thermal perturbation is considered by adding a Boltzmann sampling of thermal energy of 0.5 eV and 1.0 eV to the initial velocities of the iron and dioxygen molecule with a random phase. (3) Similar as case 1, the surface hopping dynamics is started by electronically flipped to the low-lying quintet spin state. The time step for integration of classical equations was 0.5 fs and of quantum equations, 0.0025 fs. The decoherence correction by Granucci et. al.73 was taken and the parameter is set to α=0.1 Hartree. The dynamics simulations were collected within 500 fs, and extended to 1000 fs for a few specific trajectories. For each type of initial condition, the dynamics simulations were performed for 80 trajectries, and totally 320 trajectories. All dynamics simulations were done with a recently modified version of JADE package74, interfaced with Gaussian 09 packages. The analysis of the dynamics trajectories was also done within JADE package. 2.4 Implementation Details In order to perform the on-the-fly non-adiabatic dynamics for larger molecule systems, the electronic structure calculations should be fast enough to allow the computation of tens of thousands of single points along the dynamics propagation for all trajectories, with reasonable accuracy. The potential energy surfaces (PESs) underlying this dioxygen activation process can be

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

considerably complex. It appears a very accurate electronic structure method75-76, such as CASPT2 or MRCI, is required for the theoretical studies of heme related systems58, 77-79, for which the adiabatic electronic wavefunction in the region of the MECP is a mixture of different spin states. However, the larger size molecules under study restrict the choice of suitable dimension of the active space.50-51, 59-62. Fortunately, in many cases, the wavefunction corresponding to the diabatic electronic states for heme related complex may be well described by single-reference methods, such as DFT and MP21, 22, 59, 80-83. The DFT method is commonly used to locate the reaction energy pathway, and calculate rate coefficients and isotope effects at DFT level.1, 82-83 However, it should be pointed out that the electronic structure character and ordering of spin states are very functional sensitive for iron complexes. Numerous studies have suggested that a given density functional might provide reasonable or even excellent results in some applications, but may fail in another one.50, 58, 84-91 In our experience50 and others58, 84-92, although the B3LYP used in this work is not necessarily the best functional for a specific system, its performance is generally not the worst one. Because the on-the-fly surface hopping dynamic simulations would sample various molecular geometries away from their local minima, the benchmark of their equilibrium structures may not be a good choice. Thus, we suggest that the B3LYP may offer a good compromise between the best and worst conditions. This is also the reason why this functional is still commonly used in the catalytic reaction mechanism studies.59, 80, 86, 93-96 The benchmark of best functional for dynamic simulation beyond their equilibrium structures is still an open question. We also noted the recent improvement of density functional for transition metal systems with Bayesian error estimation97 or neutral network98, which may represent recent efforts to overcome the challenge of current exchange-correlation functionals. The failure of SCF convergence is very common in the quantum chemistry calculations of iron complexes, which results in some technical difficulties for successful dynamics propagation. In this work, an extra quadratically convergent SCF procedure99 is enforced to make sure the success of the dynamics propagation, in case the first-order SCF has not converged. Of course, this would strongly lower the efficiency of the quantum chemistry calculations. Our calculations were mainly performed on Intel (R) Xeon(R) E5-2683 v3. For the full model, one trajectory (1000 fs) would require 3~4 weeks or even more, which is too time consuming, and mainly caused by the SCF

ACS Paragon Plus Environment

Page 8 of 26

Page 9 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

convergence problem. For the simplified heme mimics model, each trajectory require ~2 weeks, depending on the efficiency of the SCF procedure. Additionally, a good initial guess of the molecular orbitals is strongly required. All these considerations have been carefully coded in JADE package as a workflow to ensure the successful of dynamics simulations. As an initial illustrative application, we only present non-adiabatic surface hopping molecular dynamics simulation of a dioxygen binding iron complex (FeL2NH3 - O2) to demonstrate the feasibility of this approach. Spin − orbit coupling (SOC) matrix elements between the two interacting multiplets were calculated on the converged unrestricted molecular orbitals at B3LYP/6-31G* level using the spin-orbit Breit-Pauli Hamiltonian.100-101 Further discussions on SOC values for iron complexes can be found in the work of Shaik and others.56, 102-103 For a given pair of spin states, the obtained SOC matrix elements are condensed into an effective SOC value, which is achieved by averaging over the SOC elements. We have calculated the effective SOC value for 1000 randomly sampled structures, and the distribution of the SOC values are very concentrate (Figure S1). The SOC constants for heme related complexes are relatively small, i.e. ~100 cm-1. Thus, the SOC value could be assumed as a constant for a pair of states in the dynamics simulations, although future work is needed to explore better representations of the effective spin-orbit coupling estimations. The effective SOC values for singlet-triplet transitions are assigned to be 51.6 cm-1, and 184.6 cm-1 for between triplets and quintets, which together engender an indirect singlet/quintet mixing.104-105 3. Results and Discussions 3.1 The Structures and Properties of Heme Related Complexes Figure 2 shows the optimized geometries of the Fe(L2)(NH3)O2 complex at their low-lying singlet (S), triplet (T) and quintet (Q) states, while the optimized geometries of FeP(Im)O2 are given in Figure S2. For the singlet state, two possible geometries are obtained, one is referred to the closed shell singlet state, and the other is related to the open shell singlet state. The for the open shell singlet is 0.6289. The open shell singlet state exhibits the lowest energy, followed by the quintet, and triplet states, and the closed shell singlet state has very high energy than its counterpart open shell singlet state.1, 59, 79, 106-107 Therefore, our subsequent studies only include the open shell singlet, triplet and quintet states.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

As shown in Figure 2a, the O2 molecule prefers an end-on superoxo structure after the absorption on the iron site. The equilibrium Fe-O and Fe-N length of the reduced model is very close to the large FePIm model. The critical geometric parameters and energy order for these iron complexes were collected in Table S2. The out-of-plane angle of the porpyrin ring is within 5º, while, for the FeL2, the angle is usually in the range of ~10º. The minimum energy crossing points (MECP) obtained from gradient following method was also given in Figure 2a. The hopping geometry between the singlet and triplet state has a shorter Fe-O bond distance, meanwhile, the hopping geometry between triplet and quintet shows a quite longer Fe-O bond distance. This is consistent with the previous potential energy surface scanning calculations of Md. Ehesan Ali and co-workers108, for which the singlet/triplet spin crossover point is observed at Fe-O distance of ~1.85 Å, and the spin crossover point between higher spin states occurs at ~2.30 Å.

Figure 2. (a) The optimized geometries of the open shell singlet, triplet, and quintet states, and their MECPs. (b) The charge evolution of dioxygen moiety as a function of Fe-O distances. The spin density for the dioxygen binding and dissociation states are also shown for the open shell singlet state.

The charge transfer from the iron center to the dioxygen molecule (O2) is a critical event to facilitate subsequent formation of hydroperoxo complex (Cpd 0) in the dioxygen activation. In Figure 2b, we depict the Mulliken charge population as a function of Fe-O bond distance for each spin state. The calculation is performed by constraining the Fe-O at each distance and fully relaxing other degree of freedoms. And the charge transfer begins to happen at ~3.0 Å the from the iron center, which occurs prior to the spin state crossing points in Figure 2a. At the equilibrium

ACS Paragon Plus Environment

Page 10 of 26

Page 11 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

distance (~1.8 Å), the charge population of the O2 moiety is about 0.5 e, which is still far from 1.0 e. Thus, we prefer to roughly assign this dioxygen binding species as a mixture of FeIII–O2- and FeII–O2 species. This electronic structure feature seems to be consistent with previous theoretical calculations59, 92, 106 and Mössbauer spectra studies109-110. 3.2 Intersystem Crossing Dynamics To explore the dynamic features of the open shell singlet spin state, we have run 80 nonadiabatic surface-hopping dynamics trajectories starting from the open shell singlet state (Case 1), and 65 trajectories have successfully finished. The molecule system was found to stay on the open shell singlet state, and no spin flip event to higher spin states was observed in the limited timescale and number of trajectories. Extending the simulation timescale to 1000 fs for 10 trajectories did not show any spin flip event. Thus, the lifetime of the singlet state species should be longer than the a few picoseconds30-32, 36 Therefore, for thermal equilibrium condition, the static picture of the TSR model is quit reasonable in the assignment of catalytic reaction pathway. The coupled PES dynamics features started near the stable open shell singlet state minima without external perturbations could be viewed as a single PES process within the sub-picosecond regime. Then, the geometric features of dynamics trajectories were depicted by decomposing configuration space into a limited set of discrete states or clusters, which were sampled on the basis of K-means clustering algorithm. Figure 3a shows the overlapped molecular geometries for each cluster. As pointed out in our previous work111, the clustering algorithm is very suitable to map out the dominant long lived, kinetically meta-stable states visited by the system. Some attempts have been performed to vary the number of clusters, and finally, we adopted 4 distinct clusters to roughly represent the time series of dynamics trajectories. The standard deviations of heavy atom coordinates for each cluster after alignments are usually less than ~0.5 Å. Interestingly, the FeL2(NH3)O2 complex shows remarkable planar characters, with only a small out of plane variation (±5º). The cluster 1 mainly represents the vibrational motion near optimized singlet spin state geometry, which covers nearly 65% of dynamics trajectories, as well as the initial steps of most trajectories. The cluster 4 shows a small fraction (~5%) of nearly dissociation states, for which the dioxygen and NH3 moieties only weakly interact with the iron center. The cluster 2 and 3 show transition characters between the dioxygen binding and dissociation states, which covers 20% and 10% of dynamic trajectories.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. The dynamics simulation results, initialized from the stable open shell singlet state geometries. (a) Coarse grained representation of dynamics trajectories by the clustering algorithm. (b) The Fe-O1 v.s. Fe-O2 bond length distribution, and the end-on mode is preferred in dynamics trajectories. (c) The distribution of point charge v.s. spin population of the dioxygen moiety over dynamics trajectories.

Usually, O2 usually binds to a single iron center, with either “end-on” (η1-) or “side-on” (η2-) mode. Figure 3b shows the distribution of r(Fe-O1) v.s. r(Fe-O2) bond lengths. It is clear the end-on mode is much more common over dynamics trajectories. This is very consistent with our optimized geometries and also a few previous reports112-114. Further structural feature of the dynamics trajectories are summarized in Figure S3. Figure 3c provides the charge and spin distribution of the dioxygen moiety, averaged over dynamics trajectories. It reveals that one electron may partially transfer from the iron complex to the dioxygen moiety during the formation of Fe-O2 species. This is in accordance with charge population derived from our potential energy surface scanning of the Fe-O distance in Figure 2b. As a result the oxygen atom of the dioxygen

ACS Paragon Plus Environment

Page 12 of 26

Page 13 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

binding iron species would become more electrophilic to yield Fe-OOH complex by the H-atom abstraction. Although the spin flip event is not observed starting from the initial condition of Case 1, the potential energy gaps between different spin states are indeed small, i.e. ~0.2 eV, for a few specific trajectories (Figure S4). This energy gap change can be easily achieved by minor perturbations. For example, the high exothermicity during the dioxygen molecule binding to the iron site may provide such excessive kinetic energy. Especially, we have calculated the binding energy between the FeL2(NH3) complex and O2, which is 0.52 eV and 1.01 eV with and without considering the geometric relaxation effect. Therefore, we continue to explore the non-adiabatic features of this model system under thermal perturbations (Case 2).

Figure 4. (a) The spin state population as a function of time under thermal perturbations of 0.5eV (solid line) and 1.0 eV (dashed line). (b) The distribution of various dioxygen binding modes under the thermal perturbation of 0.5eV. (c) The distribution of various dioxygen binding modes under the thermal perturbation of 1.0 eV.

Figure 4a shows the spin states population of the FeL2(NH3)O2 complex, under two distinct thermal perturbations. The spin flip events are observed in the sub-picoseconds. At the beginning, the molecule stays at the open shell singlet state. The population of higher spin state begins to appear during the first 100 fs. The population of each spin state tends to become stable after ~200 fs. The population of the triplet state is minor, which is very consistent with our common sense that the active iron species are usually in the high spin state or low spin state.115-119 And the quintet

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

state reaches to be 10% v.s. 25% under thermal perturbations of 0.5eV and 1.0eV. The possible dynamics features of the dioxygen binding complexes are shown in Figure 4b. Two primary fates of dynamic trajectories are observed under thermal perturbations, one is the dioxygen binding state (cluster 1) and the other is the dioxygen dissociation state (cluster 4). The preference for binding v.s. dissociation state is critical for the possible H-atom abstraction, and thus the formation of iron hydroperoxo complex (Cpd 0). The rebound of the dioxygen molecule is also observed in a few trajectories (Figure S5), which may be enhanced in protein environments or condense phases, due to the steric effects.14 The end-on binding mode is still dominant in dynamic trajectories (Figure 4b, 4c). And the dissociation states are remarkably increased under thermal perturbations. The excessive kinetic energy should be an essential factor to determine the dioxygen binding/dissociation equilibrium. This confirms our assumption that the excessive kinetic energy could cause a highly-distorted molecular geometry or even cause the dissociation of the system, if the heat dissipation channels are missing. Although this process is quite fast, it is slow enough to imply substantial catalytic reaction steps. Therefore, we suspect the possible heat dissipation channels in enzymes, i.e. the dynamic interaction between the enzyme and its active center, could be very fundamental to modulate the subsequent dioxygen activation. Similar non-adiabatic dynamics effects have been commonly recognized in many photochemistry processes37-39, 42, 120. This is also consistent with the experimental reports on the rebound of dioxygen molecule in enzymes121. The possibility of the temperature dependence and non-adiabatic effects on the dioxygen activation mechanism should be considered experimentally and computationally.

ACS Paragon Plus Environment

Page 14 of 26

Page 15 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. (a) The molecular electrostatic potential (ESP) and spin density for specific dioxygen binding and dissociation trajectories. And the (b) charge and (c) spin population of the dioxygen moiety for a few typical trajectories undergoing the dioxygen binding and dissociation pathways.

The spin and charge density population change can be used as molecular descriptors to identify the activities of the dioxygen moiety. Thus, we explore the spin and charge population evolution of the dioxygen moiety for a few typical trajectories. The spin flip events are possible in these trajectories. As shown in Figure 5, the population significantly varies for the dioxygen binding and dissociation trajectories. Only the dioxygen binding trajectories are possible to undergo the subsequent H-atom abstraction to generate the iron hydroperoxo complex (Cpd 0). Therefore, the binding/dissociation equilibrium ratio between binding and dissociation structures significantly contributes to the reactivity. The spin flip events mainly contribute to modulate the binding/dissociation equilibrium during the elemental dioxygen binding reaction step. In summary, we suppose that the heme enzymes could modulate the dioxygen binding/dissociation ratio to achieve their catalytic functions, which should be considered in the comparison between the experimental and theoretical results.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. (a) The spin state population starting from quintet state, as a function of time. (b) The hopping structures from dynamics trajectories. Each bubble represents one possible position of the dioxygen or NH3 moiety. (c) The overlapped geometries of dynamics trajectories. Each dot represents one possible position of the dioxygen or NH3 moiety.

Finally, the spin flip process between singlet, triplet and quintet states was explored by artificially promoting the molecule system into the higher spin state (quintet). Figure 6a depicts the spin state population, starting from quintet state (Case 3). In this case, the molecule system undergoes a quick drop of the spin state population to lower spin states. After the quick drop within the first 200 fs, the system undergoes significant oscillation. The triplet state plays an important role to meditate the spin state crossover between the singlet and quintet state. The significant oscillation among spin states can be rationalized by the fact that the PESs are very close to each other122. Figure 6b collects the feature of hopping structures between different spin states from dynamics trajectories. Since the PESs are multidimensional, the crosses between spin states may happened in many different points, in contrast to optimized MECP points in Figure 2. Due to the electronic structure changes, the FeL2(NH3)O2 complex excited from the singlet state to the quintet state would undergo repulsion between the iron center and dioxygen moiety. Figure 5c shows the possible distribution of the dioxygen and axial moiety (NH3) near the ligand plane. The dissociation of the dioxygen moiety is a very common pathway in dynamics trajectories. The Fe-O distance is usually in the medium region (2.0~3.5 Å), which would facilitate the spin flip dynamics, and subsequent dioxygen activation reactions59, 107, 123.

ACS Paragon Plus Environment

Page 16 of 26

Page 17 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Conclusions The dioxygen (O2) binding process, as a fundamental step in many catalytic reactions, are very important in cellular respiration, corrosion, and industrial chemistry. The non-adiabatic spin crossover dynamic effects are generally ignored in profiling the PES along the reaction coordinates in catalytic communities. In this work, we try to directly establish the early-stage ultrafast spin crossover dynamics of dioxygen binding species. And such kind of iron peroxo species is the first significant oxygen intermediate before the formation of well known iron oxidants in heme enzymes, such as the iron hydroperoxo complex (Cpd 0) or iron-oxo complex (Cpd I). The on-the-fly surface hopping dynamics simulation is applied to obtain this "movie-style" depiction on the coupled low-lying spin state PESs for a heme related model. By varying the thermal perturbations over hundreds of trajectories, we suggest numerous possible fates of the dioxygen binding iron species. Further results also highlight role of spin flip events on the equilibrium of dioxygen binding and dissociation states, as a critical factor to determine the subsequent reaction rates. We suggest that increasing random thermal fluctuations may easily disrupt established reaction pathways. This provides us complimentary insights on the traditional interpreting of catalytic reaction mechanism. These results imply the emergence of non-adiabatic dynamics as a versatile tool for revealing the catalytic reaction mechanism, in addition to traditional transition state theory for interpreting catalytic findings. Recent work is also going on to study this uniquely fast and reversible O2 binding to heme by incorporating the full heme model, and also the enzymatic environments.

Acknowledges The work is supported by National Natural Science Foundation of China (Nos. 21373124, 21503249), and Huazhong Agricultural University Scientific & Technological Self-innovation Foundation (Program No. 2015RC008), and Project 2662016QD011 and 2662015PY113 Supported by the Fundamental Founds for the Central Universities. The authors also thank the support of Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase) under Grant No.U1501501.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Supporting information The distribution of the effective SOC values, the optimized geometric of reduced heme complexes, the baisis set dependece, and further details of the dynamics trajectories are given in the supporting information.

ACS Paragon Plus Environment

Page 18 of 26

Page 19 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

References 1.

Shaik, S.; Cohen, S.; Wang, Y.; Chen, H.; Kumar, D.; Thiel, W. P450 Enzymes: Their Structure,

Reactivity, and Selectivity—Modeled by QM/MM Calculations. Chem. Rev. 2010, 110 (2), 949-1017. 2.

Klotz, M. G.; Loewen, P. C. The molecular evolution of catalatic hydroperoxidases: evidence for

multiple lateral transfer of genes between prokaryota and from bacteria into eukaryota. Mol. Biol. Evol. 2003, 20 (7), 1098-112. 3.

Groves, J. T.; Wang, C. C. Nitric oxide synthase: models and mechanisms. Curr. Opin. Chem. Biol.

2000, 4 (6), 687-695. 4.

Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H. Heme-Containing Oxygenases. Chem. Rev.

1996, 96 (7), 2841. 5.

Mense, S. M.; Zhang, L. Heme: a versatile signaling molecule controlling the activities of diverse

regulators ranging from transcription factors to MAP kinases. Cell Res. 2006, 16 (8), 681-692. 6.

Solomon, E. I.; Ginsbach, J. W.; Heppner, D. E.; Kieber-Emmons, M. T.; Kjaergaard, C. H.; Smeets,

P. J.; Li, T.; Woertink, J. S. Copper Dioxygen (Bio)Inorganic Chemistry. Faraday Discuss. 2010, 148 (1), 11-39. 7.

Ray, K.; Pfaff, F. F.; Wang, B.; Nam, W. Status of Reactive Non-Heme Metal-Oxygen Intermediates

in Chemical and Enzymatic Reactions. J. Am. Chem. Soc. 2014, 136 (40), 13942-58. 8.

Finkelstein, J. Metalloproteins. Nature 2009, 460 (7257), 813.

9.

Solomon, E. I.; Heppner, D. E.; Johnston, E. M.; Ginsbach, J. W.; Cirera, J.; Qayyum, M.;

Kieber-Emmons, M. T.; Kjaergaard, C. H.; Hadt, R. G.; Tian, L. Copper Active Sites in Biology. Chem. Rev. 2014, 114 (7), 3659-3853. 10. Collman, J. P.; Gagne, R. R.; Reed, C. A.; Robinson, W. T.; Rodley, G. A. Structure of an Iron(II) Dioxygen Complex; A Model for Oxygen Carrying Hemeproteins. Proc. Natl. Acad. Sci. U. S. A. 1974, 71 (4), 1326-1329. 11. Momenteau, M.; Reed, C. A. Synthetic Heme-Dioxygen Complexes. Chem. Rev. 1994, 94 (3), 659-698. 12. Hidai, M.; Mizobe, Y. Recent Advances in the Chemistry of Dinitrogen Complexes. Cheminform 1995, 95 (4), 1115-1133. 13. MacKay, B. A.; Fryzuk, M. D. Dinitrogen Coordination Chemistry:  On the Biomimetic Borderlands. Chem. Rev. 2004, 104 (2), 385-402. 14. Holland, P. L. Metal-dioxygen and metal-dinitrogen complexes: where are the electrons? Dalton Trans. 2010, 39 (23), 5415-5425. 15. Prigge, S. T.; Eipper, B. A.; Mains, R. E.; Amzel, L. M. Dioxygen Binds End-on to Mononuclear Copper in a Precatalytic Enzyme Complex. Science 2004, 304 (5672), 864-7. 16. Alcantara, R. E.; Xu, C.; Spiro, T. G.; Guallar, V. A quantum-chemical picture of hemoglobin affinity. Proc. Natl. Acad. Sci. U. S. A. 2007, 104 (47), 18451-18455. 17. Perissinotti, L. L.; Marti, M. A.; Doctorovich, F.; Luque, F. J.; Estrin, D. A. A microscopic study of the deoxyhemoglobin-catalyzed generation of nitric oxide from nitrite anion. Biochemistry 2008, 47 (37), 9793-802. 18. Perriman, A. W.; Brogan, A. P. S.; Cölfen, H.; Tsoureas, N.; Owen, G. R.; Mann, S. Reversible dioxygen binding in solvent-free liquid myoglobin. Nat. Chem. 2010, 2 (8), 622-626. 19. Karlsson, A.; Parales, J. V.; Parales, R. E.; Gibson, D. T.; Eklund, H.; Ramaswamy, S. Crystal

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Structure of Naphthalene Dioxygenase: Side-on Binding of Dioxygen to Iron. Science 2003, 299 (5609), 1039-42. 20. Ashikawa, Y.; Fujimoto, Z.; Usami, Y.; Inoue, K.; Noguchi, H.; Yamane, H.; Nojiri, H. Structural insight into the substrate- and dioxygen-binding manner in the catalytic cycle of rieske nonheme iron oxygenase system, carbazole 1,9a-dioxygenase. BMC Struct. Biol. 2012, 12 (1), 15. 21. Schröder, D.; Shaik, S.; Schwarz, H. Two-State Reactivity as a New Concept in Organometallic Chemistry. Acc. Chem. Res. 2000, 33 (3), 139-145. 22. Poli, R.; Harvey, J. N. Spin forbidden chemical reactions of transition metal compounds. New ideas and new computational challenges. Chem. Soc. Rev. 2003, 32 (1), 1-8. 23. Hajime Hirao; Devesh Kumar; Lawrence Que, ‡ And; Sason Shaik. Two-State Reactivity in Alkane Hydroxylation by Non-Heme Iron−Oxo Complexes. J. Am. Chem. Soc. 2006, 128 (26), 8590-606. 24. Dr, D. K.; Visser, S. P. D.; Prof, S. S. Multistate Reactivity in Styrene Epoxidation by Compound I of Cytochrome P450: Mechanisms of Products and Side Products Formation. Chem-Eur J. 2005, 11 (9), 2825-35. 25. Du, L.; Gao, J.; Liu, Y.; Liu, C. Water-Dependent Reaction Pathways: An Essential Factor for the Catalysis in HEPD Enzyme. J. Phys. Chem. B 2012, 116 (39), 11837-11844. 26. Ye, S.; Geng, C. Y.; Shaik, S.; Neese, F. Electronic structure analysis of multistate reactivity in transition metal catalyzed reactions: the case of C-H bond activation by non-heme iron(IV)-oxo cores. Phys. Chem. Chem. Phys. 2013, 15 (21), 8017-8030. 27. Sun, Y.; Tang, H.; Chen, K.; Hu, L.; Yao, J.; Shaik, S.; Chen, H. Two-State Reactivity in Low-Valent Iron-Mediated C-H Activation and the Implications on Other First-Row Transition Metals. J. Am. Chem. Soc. 2016, 138 (11), 3715-3730. 28. Shaik, S.; Hirao, H.; Kumar, D. Reactivity of High-Valent Iron–Oxo Species in Enzymes and Synthetic Reagents: A Tale of Many States. Acc. Chem. Res. 2007, 40 (7), 532-542. 29. Harvey, J. N. Understanding the kinetics of spin-forbidden chemical reactions. Phys. Chem. Chem. Phys. 2007, 9 (3), 331-43. 30. Wojciech Gawelda, ‡; Andrea Cannizzo; Vanthai Pham; †, F. V. M.; ChrisPan Bressler, A.; Majed Chergui. Ultrafast Nonadiabatic Dynamics of [FeII(bpy)3]2+ in Solution. J. Am. Chem. Soc. 2007, 129 (26), 8199. 31. G, A.; M, C. Sub-50-fs photoinduced spin crossover in [Fe(bpy)₃]²⁺. Nat. Chem. 2015, 7 (8), 629. 32. Besora, M.; Carreónmacedo, J. L.; Cowan, A. J.; George, M. W.; Harvey, J. N.; Portius, P.; Ronayne, K. L.; Sun, X. Z.; Towrie, M. A combined theoretical and experimental study on the role of spin states in the chemistry of Fe(CO)5 photoproducts. J. Am. Chem. Soc. 2009, 131 (10), 3583-92. 33. Panman, M. R.; Newton, A. C.; Vos, J.; Van, d. B. B.; Bocokić, V.; Reek, J. N.; Woutersen, S. Ultrafast dynamics in iron tetracarbonyl olefin complexes investigated with two-dimensional vibrational spectroscopy. Phys. Chem. Chem. Phys. 2013, 15 (4), 1115-22. 34. Vennekate, H.; Schwarzer, D.; Torres-Alacan, J.; Krahe, O.; Filippou, A. C.; Neese, F.; Vohringer, P. Ultrafast primary processes of an iron-(iii) azido complex in solution induced with 266 nm light. Phys. Chem. Chem. Phys. 2012, 14 (18), 6165-6172. 35. Negrerie, M.; Cianetti, S.; Vos, M. H.; Martin, J.-L.; Kruglik, S. G. Ultrafast Heme Dynamics in Ferrous versus Ferric Cytochrome c Studied by Time-Resolved Resonance Raman and Transient Absorption Spectroscopy. J. Phys. Chem. B 2006, 110 (25), 12766-12781. 36. Gütlich, P.; Goodwin, H. A. Spin Crossover—An Overall Perspective. In Spin Crossover in Transition Metal Compounds I, Gütlich, P.; Goodwin, H. A., Eds. Springer Berlin Heidelberg: Berlin, Heidelberg,

ACS Paragon Plus Environment

Page 20 of 26

Page 21 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

2004; pp 1-47. 37. Harabuchi, Y.; Eng, J.; Gindensperger, E.; Taketsugu, T.; Maeda, S.; Daniel, C. Exploring the Mechanism of Ultrafast Intersystem Crossing in Rhenium(I) Carbonyl Bipyridine Halide Complexes: Key Vibrational Modes and Spin–Vibronic Quantum Dynamics. J. Chem. Theory Comput. 2016, 12 (5), 2335-2345. 38. Penfold, T. J.; Spesyvtsev, R.; Kirkby, O. M.; Minns, R. S.; Parker, D. S. N.; Fielding, H. H.; Worth, G. A. Quantum dynamics study of the competing ultrafast intersystem crossing and internal conversion in the “channel 3” region of benzene. J. Chem. Phys. 2012, 137 (20), 204310. 39. Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. J. Chem. Theory Comput. 2011, 7 (5), 1253-1258. 40. Fedorov, D. A.; Pruitt, S. R.; Keipert, K.; Gordon, M. S.; Varganov, S. A. Ab Initio Multiple Spawning Method for Intersystem Crossing Dynamics: Spin-Forbidden Transitions between 3B1 and 1A1 States of GeH2. J. Phys. Chem. A 2016, 120 (18), 2911-2919. 41. Cui, G.; Thiel, W. Generalized trajectory surface-hopping method for internal conversion and intersystem crossing. J. Chem. Phys. 2014, 141 (12), 124101. 42. Carvalho, F. F. d.; Tavernelli, I. Nonadiabatic dynamics with intersystem crossings: A time-dependent density functional theory implementation. J. Chem. Phys. 2015, 143 (22), 224105. 43. Bai, S.; Barbatti, M. Spatial Factors for Triplet Fusion Reaction of Singlet Oxygen Photosensitization. J. Phys. Chem. Lett. 2017, 8 (21), 5456-5460. 44. Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011–2015. J. Phys. Chem. Lett. 2016, 7 (11), 2100-2112. 45. Mai, S.; Marquetand, P.; González, L. A general method to describe intersystem crossing dynamics in trajectory surface hopping. Int. J. Quantum Chem. 2015, 115 (18), 1215-1231. 46. Richter, M.; Mai, S.; Marquetand, P.; González, L. Ultrafast intersystem crossing dynamics in uracil unravelled by ab initio molecular dynamics. Phys. Chem. Chem. Phys. 2014, 16 (44), 24423-36. 47. Martínez-Fernández, L.; González-Vázquez, J.; González, L.; Corral, I. Time-Resolved Insight into the Photosensitized Generation of Singlet Oxygen in Endoperoxides. J. Chem. Theory Comput. 2015, 11 (2), 406-414. 48. Xu, C.; Yu, L.; Zhu, C.; Yu, J.; Cao, Z. Intersystem crossing-branched excited-state intramolecular proton transfer for o-nitrophenol: An ab initio on-the-fly nonadiabatic molecular dynamic simulation. Scientific Reports 2016, 6, 26768. 49. Berryman, V. E. J.; Boyd, R. J.; Johnson, E. R. Balancing Exchange Mixing in Density-Functional Approximations for Iron Porphyrin. J. Chem. Theory Comput. 2015, 11 (7), 3022-3028. 50. Zhao, H.; Fang, C.; Gao, J.; Liu, C. Spin-state energies of heme-related models from spin-flip TDDFT calculations. Phys. Chem. Chem. Phys. 2016, 18 (42), 29486-29494. 51. Radoń, M. Spin-State Energetics of Heme-Related Models from DFT and Coupled Cluster Calculations. J. Chem. Theory Comput. 2014, 10 (6), 2306-2321. 52. Pratihar, S.; Ma, X.; Homayoon, Z.; Barnes, G. L.; Hase, W. L. Direct Chemical Dynamics Simulations. J. Am. Chem. Soc. 2017, 139 (10), 3570-3590. 53. Yu, P.; Yang, Z.; Liang, Y.; Hong, X.; Li, Y.; Houk, K. N. Distortion-Controlled Reactivity and Molecular Dynamics of Dehydro-Diels–Alder Reactions. J. Am. Chem. Soc. 2016, 138 (26), 8247-8252. 54. Scherlis, D. A.; Estrin, D. A. Structure and spin﹕tate energetics of an iron porphyrin model: An assessment of theoretical methods. Int. J. Quantum Chem. 2002, 87 (3), 158-166.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

55. Radoń, M.; Broclawik, E.; Pierloot, K. DFT and Ab Initio Study of Iron-Oxo Porphyrins: May They Have a Low-Lying Iron(V)-Oxo Electromer? J. Chem. Theory Comput. 2011, 7 (4), 898-908. 56. Strickland, N.; Harvey, J. N. Spin-Forbidden Ligand Binding to the Ferrous−Heme Group:  Ab IniPo and DFT Studies. J. Phys. Chem. B 2007, 111 (4), 841-852. 57. Oláh, J.; Harvey, J. N. NO Bonding to Heme Groups: DFT and Correlated ab Initio Calculations. J. Phys. Chem. A 2009, 113 (26), 7338-7345. 58. Vancoillie, S.; Zhao, H.; Radoń, M.; Pierloot, K. Performance of CASPT2 and DFT for Relative Spin-State Energetics of Heme Models. J. Chem. Theory Comput. 2010, 6 (2), 576-582. 59. Jensen, K. P.; Ryde, U. How O2 binds to heme: reasons for rapid binding and spin inversion. J. Biol. Chem. 2004, 279 (15), 14561-14569. 60. Borowski, T.; Siegbahn, P. E. M. Mechanism for Catechol Ring Cleavage by Non-Heme Iron Intradiol Dioxygenases:  A Hybrid DFT Study. J. Am. Chem. Soc. 2006, 128 (39), 12941-12953. 61. Yoshioka, Y.; Satoh, H.; Mitani, M. Theoretical study on electronic structures of FeOO, FeOOH, FeO(H2O), and FeO in hemes: as intermediate models of dioxygen reduction in cytochrome c oxidase. J. Inorg. Biochem. 2007, 101 (10), 1410. 62. Blomberg, M. R. A.; Borowski, T.; Himo, F.; Liao, R.-Z.; Siegbahn, P. E. M. Quantum Chemical Studies of Mechanisms for Metalloenzymes. Chem. Rev. 2014, 114 (7), 3601-3658. 63. Harvey, J. N.; Aschi, M.; Schwarz, H.; Koch, W. The singlet and triplet states of phenyl cation. A hybrid approach for locating minimum energy crossing points between non-interacting potential energy surfaces. Theor. Chem. Acc. 1998, 99 (2), 95-99. 64. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al, Gaussian 09, Gaussian, Inc.: Wallingford, CT, USA, 2009. 65. Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of H+ with D2. J. Chem. Phys. 1971, 55 (2), 562-572. 66. Hammes‐Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. J. Chem. Phys. 1994, 101 (6), 4657-4667. 67. Richter, M.; Marquetand, P.; González-Vázquez, J.; Sola, I.; González, L. Femtosecond Intersystem Crossing in the DNA Nucleobase Cytosine. J. Phys. Chem. Lett. 2012, 3 (21), 3090-3095. 68. Favero, L.; Granucci, G.; Persico, M. Dynamics of acetone photodissociation: a surface hopping study. Phys. Chem. Chem. Phys. 2013, 15 (47), 20651-61. 69. Martínezfernández, L.; Corral, I.; Granucci, G.; Persico, M. Competing ultrafast intersystem crossing and internal conversion: a time resolved picture for the deactivation of 6-thioguanine. Chem. Sci. 2014, 5 (4), 1336-1347. 70. Mai, S.; Marquetand, P.; González, L. Non-adiabatic and intersystem crossing dynamics in SO2. II. The role of triplet states in the bound state dynamics studied by surface-hopping simulations. J. Chem. Phys. 2014, 140 (20), 204302. 71. Granucci, G.; Persico, M.; Spighi, G. Surface hopping trajectory simulations with spin-orbit and dynamical couplings. J. Chem. Phys. 2012, 137 (22), 22A501. 72. Habenicht, B. F.; Prezhdo, O. V. Ab Initio Time-Domain Study of the Triplet State in a Semiconducting Carbon Nanotube: Intersystem Crossing, Phosphorescence Time, and Line Width. J. Am. Chem. Soc. 2012, 134 (38), 15648-15651. 73. Granucci, G.; Persico, M.; Zoccante, A. Including quantum decoherence in surface hopping. J. Chem. Phys. 2010, 133 (13), 134111.

ACS Paragon Plus Environment

Page 22 of 26

Page 23 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

74. Du, L.; Lan, Z. An On-the-Fly Surface-Hopping Program JADE for Nonadiabatic Molecular Dynamics of Polyatomic Systems: Implementation and Applications. J. Chem. Theory Comput. 2015, 11 (4), 1360. 75. WERNER, H. Third-order multireference perturbation theory The CASPT3 method. Mol. Phys. 1996, 89 (2), 645-661. 76. Werner, H. J.; Knowles, P. J. An efficient internally contracted multiconfiguration–reference configuration interaction method. J. Chem. Phys. 1988, 89 (9), 5803-5814. 77. Meunier, B.; de Visser, S. P.; Shaik, S. Mechanism of oxidation reactions catalyzed by cytochrome p450 enzymes. Chem. Rev. 2004, 104 (9), 3947-80. 78. Saitow, M.; Kurashige, Y.; Yanai, T. Fully Internally Contracted Multireference Configuration Interaction Theory Using Density Matrix Renormalization Group: A Reduced-Scaling Implementation Derived by Computer-Aided Tensor Factorization. J. Chem. Theory Comput. 2015, 11 (11), 5120-5131. 79. Shaik, S.; Chen, H. Lessons on O2 and NO bonding to heme from ab initio multireference/multiconfiguration and DFT calculations. JBIC J. Biol. Inorg. Chem. 2011, 16 (6), 841-855. 80. Dedieu, A.; Rohmer, M. M.; Benard, M.; Veillard, A. Oxygen binding to iron porphyrins. An ab initio calculation. J. Am. Chem. Soc. 1976, 98 (12), 3717-3718. 81. Collman, J. P.; Boulatov, R.; Sunderland, C. J.; Fu, L. Functional analogues of cytochrome c oxidase, myoglobin, and hemoglobin. Chem. Rev. 2004, 104 (2), 561-88. 82. Lawson Daku, L. M.; Aquilante, F.; Robinson, T. W.; Hauser, A. Accurate Spin-State Energetics of Transition Metal Complexes. 1. CCSD(T), CASPT2, and DFT Study of [M(NCH)6]2+ (M = Fe, Co). J. Chem. Theory Comput. 2012, 8 (11), 4216-4231. 83. Ortiz de Montellano, P. R. Hydrocarbon Hydroxylation by Cytochrome P450 Enzymes. Chem. Rev. 2010, 110 (2), 932-948. 84. Gani, T. Z. H.; Kulik, H. J. Unifying Exchange Sensitivity in Transition-Metal Spin-State Ordering and Catalysis through Bond Valence Metrics. J. Chem. Theory Comput. 2017, 13 (11), 5443-5457. 85. Janet, J. P.; Kulik, H. J. Resolving Transition Metal Chemical Space: Feature Selection for Machine Learning and Structure–Property Relationships. J. Phys. Chem. A 2017, 121 (46), 8939-8954. 86. Cramer, C. J.; Truhlar, D. G. Density functional theory for transition metals and transition metal chemistry. Phys. Chem. Chem. Phys. 2009, 11 (46), 10757-10816. 87. Bowman, D. N.; Jakubikova, E. Low-Spin versus High-Spin Ground State in Pseudo-Octahedral Iron Complexes. Inorg. Chem. 2012, 51 (11), 6011-6019. 88. Roy, L. E.; Jakubikova, E.; Guthrie, M. G.; Batista, E. R. Calculation of One-Electron Redox Potentials Revisited. Is It Possible to Calculate Accurate Potentials with Density Functional Methods? J. Phys. Chem. A 2009, 113 (24), 6745-6750. 89. Mortensen, S. R.; Kepp, K. P. Spin Propensities of Octahedral Complexes From Density Functional Theory. J. Phys. Chem. A 2015, 119 (17), 4041-4050. 90. Kepp, K. P. Theoretical Study of Spin Crossover in 30 Iron Complexes. Inorg. Chem. 2016, 55 (6), 2717-2727. 91. Weymuth, T.; Couzijn, E. P. A.; Chen, P.; Reiher, M. New Benchmark Set of Transition-Metal Coordination Reactions for the Assessment of Density Functionals. J. Chem. Theory Comput. 2014, 10 (8), 3092-3103. 92. Rovira, C.; Kunc, K.; Hutter, J.; Ballone, P.; Parrinello, M. Equilibrium Geometries and Electronic Structure of Iron−Porphyrin Complexes:  A Density FuncPonal Study. J. Phys. Chem. A 1997, 101 (47),

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

8914-8925. 93. Ghosh, A. First-Principles Quantum Chemical Studies of Porphyrins. Acc. Chem. Res. 1998, 31 (4), 189-198. 94. Johansson, M. P.; Sundholm, D.; Gerfen, G.; Wikström, M. The Spin Distribution in Low-Spin Iron Porphyrins. J. Am. Chem. Soc. 2002, 124 (39), 11771-11780. 95. Jensen, K. P.; Ryde, U. Comparison of the chemical properties of iron and cobalt porphyrins and corrins. Chembiochem 2003, 4 (5), 413–424. 96. Yang, T.; Quesne, M. G.; Neu, H. M.; Cantú Reinhard, F. G.; Goldberg, D. P.; de Visser, S. P. Singlet versus Triplet Reactivity in an Mn(V)–Oxo Species: Testing Theoretical Predictions Against Experimental Evidence. J. Am. Chem. Soc. 2016, 138 (38), 12375-12386. 97. Wellendorff, J.; Lundgaard, K. T.; Jacobsen, K. W.; Bligaard, T. mBEEF: An accurate semi-local Bayesian error estimation density functional. J. Chem. Phys. 2014, 140 (14), 144107. 98. Janet, J. P.; Kulik, H. J. Predicting electronic structure properties of transition metal complexes with neural networks. Chem. Sci. 2017, 8 (7), 5137-5152. 99. Bacskay, G. B. A quadratically convergent Hartree—Fock (QC-SCF) method. Application to closed shell systems. Chem. Phys. 1981, 61 (3), 385-404. 100. Chiodo, S. G.; Leopoldini, M. MolSOC : A spin–orbit coupling code ☆. Comput. Phys. Commun. 2014, 185 (2), 676-683. 101. Heß, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. A mean-field spin-orbit method applicable to correlated wavefunctions. Chem. Phys. Lett. 1996, 251 (5–6), 365-371. 102. Danovich, D.; Shaik, S. Spin−Orbit Coupling in the OxidaPve AcPvaPon of H−H by FeO+. Selection Rules and Reactivity Effects. J. Am. Chem. Soc. 1997, 119 (7), 1773-1786. 103. Shiota, Y.; Yoshizawa, K. A spin–orbit coupling study on the spin inversion processes in the direct methane-to-methanol conversion by FeO+. J. Chem. Phys. 2003, 118 (13), 5872-5879. 104. Li, Z.; Suo, B.; Zhang, Y.; Xiao, Y.; Liu, W. Combining spin-adapted open-shell TD-DFT with spin–orbit coupling. Mol. Phys. 2013, 111 (24), 3741-3755. 105. Cheng, L.; Gauss, J. Perturbative treatment of spin-orbit coupling within spin-free exact two-component theory. J. Chem. Phys. 2014, 141 (16), 164107. 106. Chen, H.; Ikeda-Saito, M.; Shaik, S. Nature of the Fe−O2 Bonding in Oxy-Myoglobin: Effect of the Protein. J. Am. Chem. Soc. 2008, 130 (44), 14778-14790. 107. Kepp, K. P. O2 Binding to Heme is Strongly Facilitated by Near-Degeneracy of Electronic States. ChemPhysChem 2013, 14 (15), 3551-3558. 108. Ali, M. E.; Sanyal, B.; Oppeneer, P. M. Electronic Structure, Spin-States, and Spin-Crossover Reaction of Heme-Related Fe-Porphyrins: A Theoretical Perspective. J. Phys. Chem. B 2012, 116 (20), 5849-5859. 109. Tsai, T. E.; Groves, J. L.; Wu, C. S. Electronic structure of iron–dioxygen bond in oxy‐Hb‐A and its isolated oxy‐α and oxy‐β chains. J. Chem. Phys. 1981, 74 (8), 4306-4314. 110. Oshtrakh, M. I. Study of the relationship of small variations of the molecular structure and the iron state in iron containing proteins by Mössbauer spectroscopy: biomedical approach. Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 2004, 60 (1), 217-234. 111. Liu, F.; Du, L.; Zhang, D.; Gao, J. Direct Learning Hidden Excited State Interaction Patterns from ab initio Dynamics and Its Implication as Alternative Molecular Mechanism Models. Scientific Reports 2017, 7 (1), 8737. 112. Roelfes, G.; Vrajmasu, V.; Chen, K.; Ho, R. Y. N.; Rohde, J.-U.; Zondervan, C.; la Crois, R. M.;

ACS Paragon Plus Environment

Page 24 of 26

Page 25 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Schudde, E. P.; Lutz, M.; Spek, A. L.; Hage, R.; Feringa, B. L.; Münck, E.; Que, L. End-On and Side-On Peroxo Derivatives of Non-Heme Iron Complexes with Pentadentate Ligands:  Models for Putative Intermediates in Biological Iron/Dioxygen Chemistry. Inorg. Chem. 2003, 42 (8), 2639-2653. 113. Ohta, T.; Liu, J.-G.; Naruta, Y. Resonance Raman characterization of mononuclear heme-peroxo intermediate models. Coordin. Chem. Rev. 2013, 257 (2), 407-413. 114. Li, L.; Chang, Z.; Pan, Z.; Fu, Z. Q.; Wang, X. Modes of heme binding and substrate access for cytochrome P450 CYP74A revealed by crystal structures of allene oxide synthase. Proc. Natl. Acad. Sci. U. S. A. 2008, 105 (37), 13883-8. 115. Cheesman, M. R.; Ferguson, S. J.; Moir, J. W. B.; Richardson, D. J.; Zumft, W. G.; Thomson, A. J. Two Enzymes with a Common Function but Different Heme Ligands in the Forms as Isolated. Optical and Magnetic Properties of the Heme Groups in the Oxidized Forms of Nitrite Reductase, Cytochrome cd1, from Pseudomonas stutzeri and Thiosphaera pantotropha. Biochemistry 1997, 36 (51), 16267-16276. 116. Solomon, E. I.; Decker, A.; Lehnert, N. Non-heme iron enzymes: contrasts to heme catalysis. Proc. Natl. Acad. Sci. U. S. A. 2003, 100 (7), 3589. 117. Vanwetswinkel, S.; van Nuland, N. A.; Volkov, A. N. Paramagnetic properties of the low- and high-spin states of yeast cytochrome c peroxidase. J. Biomol. NMR 2013, 57 (1), 21-26. 118. Kepp, K. P. Heme: From quantum spin crossover to oxygen manager of life. Coordin. Chem. Rev. 2017, 344 (Supplement C), 363-374. 119. Schuth, N.; Mebs, S.; Huwald, D.; Wrzolek, P.; Schwalbe, M.; Hemschemeier, A.; Haumann, M. Effective intermediate-spin iron in O2-transporting heme proteins. Proc. Natl. Acad. Sci. U. S. A. 2017, 114 (32), 8556-8561. 120. Liu, F.; Du, L.; Lan, Z.; Gao, J. Hydrogen bond dynamics governs the effective photoprotection mechanism of plant phenolic sunscreens. Photochem. Photobiol. Sci. 2017, 16 (2), 211-219. 121. Lundberg, M.; Morokuma, K. Protein environment facilitates O2 binding in non-heme iron enzyme. An insight from ONIOM calculations on isopenicillin N synthase (IPNS). J. Phys. Chem. B 2007, 111 (31), 9380-9. 122. de Groot, B. L.; Daura, X.; Mark, A. E.; Grubmüller, H. Essential dynamics of reversible peptide folding: memory-free conformational dynamics governed by internal hydrogen bonds. J. Mol. Biol. 2001, 309 (1), 299. 123. Nakashima, H.; Hasegawa, J. Y.; Nakatsuji, H. On the reversible O2 binding of the Fe-porphyrin complex. J. Comput. Chem. 2006, 27 (4), 426-433.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic

ACS Paragon Plus Environment

Page 26 of 26