Direct Simulation of Proton-Coupled Electron Transfer Reaction Dynamics and Mechanisms Citation Kretchmer, Joshua S. (2015) Direct Simulation of Proton-Coupled Electron Transfer Reaction Dynamics and Mechanisms. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/Z98913T5. https://resolver.caltech.edu/CaltechTHESIS:01202015-085910230 Abstract Proton-coupled electron transfer (PCET) reactions are ubiquitous throughout chemistry and biology. However, challenges arise in both the the experimental and theoretical investigation of PCET reactions; the rare-event nature of the reactions and the coupling between quantum mechanical electron- and proton-transfer with the slower classical dynamics of the surrounding environment necessitates the development of robust simulation methodology. In the following dissertation, novel path-integral based methods are developed and employed for the direct simulation of the reaction dynamics and mechanisms of condensed-phase PCET. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Charge Transfer, Molecular Dynamics, Non-adiabatic Dynamics, Proton Transfer, Electron Transfer Degree Grantor: California Institute of Technology Division: Chemistry and Chemical Engineering Major Option: Chemistry Thesis Availability: Public (worldwide access) Research Advisor(s): Miller, Thomas F. Thesis Committee: Marcus, Rudolph A. (chair) Okumura, Mitchio Gray, Harry B. Miller, Thomas F. Defense Date: 17 December 2014 Record Number: CaltechTHESIS:01202015-085910230 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:01202015-085910230 DOI: 10.7907/Z98913T5 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 8756 Collection: CaltechTHESIS Deposited By: Joshua Kretchmer Deposited On: 27 Jan 2015 22:41 Last Modified: 04 Oct 2019 00:07 Thesis Files Preview PDF - Final Version See Usage Policy. 4MB Repository Staff Only: item control page

Direct simulation of proton-coupled electron transfer reaction dynamics and mechanisms Thesis by Joshua S. Kretchmer In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2015 (Defended December 17, 2014) ii c 2015 Joshua S. Kretchmer All Rights Reserved iii Acknowledgments First I would like to thank Professor Thomas F. Miller, my advisor. Over the course of my graduate career he has provided me with invaluable guidance as a scientist and as a member of the scientific community. His mentorship has forever shaped how I view the world. I would also like to thank the faculty of the division of chemistry. In particular, I would like to thank my thesis committee, Professor Rudy Marcus, Professor Mitchio Okumura, and Professor Harry Gray, for the their continued support and guidance throughout my graduate career. Their advice regarding both my science and my future has been indispensable. Thank you to the division administrative staff, in particular Priscilla Boon. Without Priscilla, the Miller group would come to a grinding halt. I appreciate the support and friendship of the entirety of the Miller Group. Bin, Artur, Jason, Nick, Taylor, Connie, Mike, Michiel, Joonho, Mark, Romelia, Nandini, Kuba, Frank, Fran, Brett, and Eric, your support and willingness to provide comments for talks and papers, help with scientific questions, and generally be amazing office-mates has made my graduate career all the better. Thank you to Marco, Katie, Tim, and Jeff, my housemates and friends who have been there when I got home on both the good days and the bad. Thank you to my family, my parents and my sister, who have provided love and support throughout my entire life. Without their encouragement I would not be where I am today. And thank you to Danielle. You are always there for me, supporting me and challenging me to be a better person. iv Abstract Proton-coupled electron transfer (PCET) reactions are ubiquitous throughout chemistry and biology. However, challenges arise in both the the experimental and theoretical investigation of PCET reactions; the rare-event nature of the reactions and the coupling between quantum mechanical electron- and proton-transfer with the slower classical dynamics of the surrounding environment necessitates the development of robust simulation methodology. In the following dissertation, novel path-integral based methods are developed and employed for the direct simulation of the reaction dynamics and mechanisms of condensed-phase PCET. Chapter 1 describes the extension of the ring-polymer molecular dynamics (RPMD) method, a path-integral based method, for the direct simulation of PCET reactions across a wide range of physically relevant regimes, including both concerted and sequential PCET and concerted PCET in the fully non-adiabatic, partially adiabatic, and fully adiabatic regimes. Analysis of the RPMD trajectories reveals distinct kinetic pathways associated with the sequential and concerted PCET reaction mechanisms, and it is demonstrated that concerted PCET proceeds by a solvent-gating mechanism in which the reorganization energy is mitigated by charge cancellation among the transferring particles. The RPMD methodology is further verified through comparison to previously derived rate-theories, and it is particularly notable that RPMD accurately predicts the crossover in the thermal reaction rates between different coupling regimes while avoiding any prior assumptions about the reaction mechanism. Analysis of the RPMD trajectories as a function of the electronic coupling reveals an intriguing and distinct “transient-proton-bridge” mechanism for intermediate values of the electronic coupling. This work has been published as Kretchmer, J. S. and Miller, T. F., “Direct simulation of proton-coupled electron transfer across multiple regimes,” Journal of Chemical Physics, 138, 134109 (2013). Chapter 2 describes the use of RPMD to investigate the competition between concerted and sequential PCET in the organometallic complex iron bi-imidazoline, which is a model compound for non-heme iron containing enzymes involved in C-H bond activation. The calculation of the thermal reaction rates for the concerted and sequential PCET mechanism in iron bi-imidazoline predicts that the concerted mechanism is the dominant PCET mechanism. Analysis of the RPMD trajectories v reveal that the favorability of the PCET mechanism in iron bi-imidazoline is associated with (i) the avoidance of the charge-separated intermediates involved in the sequential mechanism, (ii) the mitigation of the barrier associated with solvent reorganization energy due to charge cancellation between the transferring electron and the time-dependent charge on the bi-imidazoline ligands, and (iii) the small barrier associated with PT during PCET due to the short hydrogen-bond in iron bi-imidazoline. Further results from RPMD simulations then predict that the dominant PCET mechanism can be altered by varying (i) the strength of interaction between the electron and the proton through spatial separation of the electron and proton donor sites, (ii) the polarity of the solvent, and (iii) the height of the PT barrier through an increase in the PT distance. This work is in preparation for submission as Kretchmer, J. S. and Miller, T. F. “Understanding and controlling concerted versus sequential proton-coupled electron transfer”. Chapter 3 explores the performance of the recently developed kinetically-constrained (KC) RPMD for the investigation of condensed phase electron transfer (ET) reactions in which the dynamics of the nuclei play a large role in determining the reaction rate and mechanism. The KC-RPMD method extends conventional RPMD to allow for the treatment of general multi-electron, non-adiabatic processes. Through this work, several technical aspects of the KC-RPMD derivation, including the kinetic-constraint and the choice of mass associated with the continuous electronic coordinate, are revisited. In particular, the presence of spurious resonances in the continuous electronic coordinate is corrected through the inclusion of a Langevin thermostat that is attached to the electronic coordinate. The accuracy of the KC-RPMD dynamics is illustrated through the treatment of a wide-range of systems in which nuclear dynamics are important, including systems in which the diabatic coupling depends strongly on the position of the nuclei and systems with varying magnitudes of the friction of the solvent coordinate; KC-RPMD simulations are able to properly capture the transition between the friction-indepedent Marcus regime of ET and both the high-friction Zusman regime and the low-friction weakly-dissipative regime of ET. This work is in preparation for submission as Kretchmer, J. S. and Miller, T. F. “The role of nuclear dynamics during electron transfer in the context of kinetically-constrained ring polymer molecular dynamics”. vi Contents Acknowledgments iii Abstract iv 1 Direct simulation of proton-coupled electron transfer across multiple regimes 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Ring polymer molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 PCET Rate Theories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Concerted PCET in the Fully Adiabatic Regime . . . . . . . . . . . . . . . . 7 1.3.2 Concerted PCET in the Partially Adiabatic Regime . . . . . . . . . . . . . . 7 1.3.3 Concerted PCET in the Fully Non-Adiabatic Regime . . . . . . . . . . . . . . 8 1.3.4 ET rate theories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.4.1 Adiabatic ET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.4.2 Non-adiabatic ET . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4 PCET Model Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 Calculation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.5.1 RPMD Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.5.1.1 RPMD rate calculations for concerted PCET . . . . . . . . . . . . . 13 1.5.1.2 RPMD rate calculations for ET prior to PT . . . . . . . . . . . . . 15 1.5.1.3 RPMD rate calculations for PT prior to ET . . . . . . . . . . . . . 15 1.5.1.4 Two-dimensional FE profiles . . . . . . . . . . . . . . . . . . . . . . 16 1.5.1.5 RPMD transition path ensemble . . . . . . . . . . . . . . . . . . . . 17 1.5.1.6 RPMD rate calculations for pure ET . . . . . . . . . . . . . . . . . 17 1.5.2 PCET Rate Theory Calculations . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.5.3 ET Rate Theory Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.6.1 Sequential versus concerted PCET . . . . . . . . . . . . . . . . . . . . . . . . 22 1.6.2 Reactions across multiple coupling regimes . . . . . . . . . . . . . . . . . . . 29 1.6 vii 1.7 1.6.2.1 ET across electronic-coupling regimes . . . . . . . . . . . . . . . . . 29 1.6.2.2 Concerted PCET across electronic-coupling regimes . . . . . . . . . 32 1.6.2.3 Concerted PCET across vibrational-coupling regimes . . . . . . . . 35 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Appendix A Coefficients for the Electron Potential . . . . . . . . . . . . . . . . . . . . . 39 Appendix B System-Bath Potential Energy Parameters . . . . . . . . . . . . . . . . . . . 39 Appendix C Auxiliary Restraining Potentials . . . . . . . . . . . . . . . . . . . . . . . . 41 1.C.1 Auxiliary Restraining Potential for Concerted PCET . . . . . . . . . . . . . . 41 1.C.2 Auxiliary Restraining Potential for ET prior to PT . . . . . . . . . . . . . . . 41 Appendix D Parameters for Fitting Functions . . . . . . . . . . . . . . . . . . . . . . . . 42 Appendix E Diabatization Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2 Understanding and controlling concerted versus sequential proton-coupled electron transfer 51 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.2 Ring polymer molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.3 Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3.1 Atomistic Representation for PCET . . . . . . . . . . . . . . . . . . . . . . . 57 2.3.2 System-bath Representation for PCET . . . . . . . . . . . . . . . . . . . . . . 62 Calculation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.4.1 Atomistic representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.4.1.1 RPMD rate calculations for concerted PCET . . . . . . . . . . . . . 66 2.4.1.2 RPMD rate calculations for ET prior to PT . . . . . . . . . . . . . 68 2.4.1.3 1D FE profile for PT prior to ET . . . . . . . . . . . . . . . . . . . 69 2.4.1.4 Two-dimensional FE profiles . . . . . . . . . . . . . . . . . . . . . . 69 2.4.1.5 Solvent reorganization energy for concerted PCET . . . . . . . . . . 70 2.4.1.6 Solvent reorganization energy for ET prior to PT . . . . . . . . . . 71 2.4.1.7 Solvent reorganization energy for symmetric ET . . . . . . . . . . . 71 2.4.1.8 RPMD transition path ensemble . . . . . . . . . . . . . . . . . . . . 72 System-bath representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 2.4.2.1 RPMD rate calculations for concerted PCET . . . . . . . . . . . . . 73 2.4.2.2 RPMD rate calculations for ET prior to PT . . . . . . . . . . . . . 74 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.4 2.4.2 2.5 viii 2.6 2.5.1 Understanding concerted versus sequential PCET in iron bi-imidazoline . . . 75 2.5.2 Controlling the dominant PCET Mechanism . . . . . . . . . . . . . . . . . . 81 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Appendix A Coefficients for the electron potential . . . . . . . . . . . . . . . . . . . . . . 89 Appendix B System-bath potential energy parameters . . . . . . . . . . . . . . . . . . . 89 Appendix C Auxiliary restraining potentials . . . . . . . . . . . . . . . . . . . . . . . . . 89 2.C.1 Auxiliary restraining potential for concerted PCET in the atomistic models . 90 2.C.2 Auxiliary restraining potential for ET prior to PT in the atomistic models . . 91 2.C.3 Auxiliary restraining potential for concerted PCET in the system-bath models 91 2.C.4 Auxiliary restraining potential for ET prior to PT in the system-bath models 92 3 The role of nuclear dynamics during electron transfer in the context of kineticallyconstrained ring polymer molecular dynamics 98 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.3 Model Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.4 Calculation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 3.5 3.6 3.4.1 KC-RPMD rate calculation in System A1 . . . . . . . . . . . . . . . . . . . . 108 3.4.2 KC-RPMD rate calculation in System A2 . . . . . . . . . . . . . . . . . . . . 109 3.4.3 KC-RPMD rate calculation in System A3 . . . . . . . . . . . . . . . . . . . . 109 3.4.4 KC-RPMD rate calculation in System B . . . . . . . . . . . . . . . . . . . . . 110 3.4.5 Calculation of reference TST expressions . . . . . . . . . . . . . . . . . . . . . 110 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 3.5.1 Normal to inverted electron transfer . . . . . . . . . . . . . . . . . . . . . . . 111 3.5.2 Friction controlled electron transfer . . . . . . . . . . . . . . . . . . . . . . . . 112 3.5.3 Electron transfer with position dependent coupling . . . . . . . . . . . . . . . 114 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Appendix A Derivation of the penalty function . . . . . . . . . . . . . . . . . . . . . . . 117 3.A.1 1D redox system with constant K and classical nuclei . . . . . . . . . . . . . 117 3.A.2 Multi-dimensional redox system with position-dependent K(R) . . . . . . . . 119 Appendix B Derivation of the Langevin friction coefficient . . . . . . . . . . . . . . . . . 120 ix List of Figures 1.1 Schematic illustration of competing PCET reaction mechanisms. . . . . . . . . . . . . 2 1.2 Representative RPMD trajectories for concerted and sequential PCET. . . . . . . . . 23 1.3 Free-energy profiles along PCET reaction channels. . . . . . . . . . . . . . . . . . . . . 24 1.4 Representative RPMD trajectories illustrating solvent fluctuations for concerted PCET. 27 1.5 ET in the non-adiabatic and adiabatic regimes. . . . . . . . . . . . . . . . . . . . . . . 30 1.6 Concerted PCET in the fully non-adiabatic and partially adiabatic regimes. . . . . . . 33 1.7 Concerted PCET in the partially adiabatic and fully adiabatic regimes. . . . . . . . . 36 2.1 Schematic illustration of PCET in iron bi-imidazoline . . . . . . . . . . . . . . . . . . 52 2.2 Representative RPMD trajectories for concerted and sequential PCET in iron biimidazoline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 2.3 Free-energy profiles along PCET reaction channels in iron bi-imidazoline. . . . . . . . 77 2.4 Representative RPMD trajectories illustrating solvent fluctuations for concerted PCET in iron bi-imidazoline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 79 The time-dependent system dipole of iron bi-imidazoline during the concerted PCET reaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 2.6 Solvent dependence of the PCET mechanism in iron bi-imidazoline. . . . . . . . . . . 82 2.7 Solvent dependence of the PCET mechanism in systems with weaker electron-proton interactions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 2.8 Solvent dependence of the PCET mechanism in system-bath models of PCET. . . . . 85 2.9 Dependence of the PCET mechanism in system-bath models of PCET on the strength of the electron-proton interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 86 Dependence of the PCET mechanism in system-bath models of PCET on the height of the PT barrier. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.1 ET reaction rate coefficients as a function of the driving force . . . . . . . . . . . . . . 112 3.2 Transition from Marcus to Zusman regime of electron transfer . . . . . . . . . . . . . 112 3.3 Transition from the weakly dissipative to the Marcus regime of electron transfer . . . 113 x 3.4 ET rates as a function of the frequency of a donor-acceptor mode . . . . . . . . . . . 115 xi List of Tables 1.1 Values of the electronic coupling, VET , vibrational coupling, VPT , and reorganization energy, λ, for the system-bath model systems for PCET. . . . . . . . . . . . . . . . . . 1.2 Values of the electronic coupling, VET , and reorganization energy, λ, for ET systems that vary between the adiabatic and non-adiabatic regimes. . . . . . . . . . . . . . . . 1.3 25 Reaction rates for the full ET-PT, PT-ET, and concerted PCET mechanisms calculated using RPMD and Eqs. (1.45) and (1.46). . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 22 RPMD rates for the forward and reverse ET and PT reactions in the sequential mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 21 26 Parameters for the donor coulombic well in the intrinsic electron potential energy function of Eq. (1.26). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.6 Parameters for the PCET and ET systems. . . . . . . . . . . . . . . . . . . . . . . . . 40 1.7 Potential energy function parameters for the PCET systems. . . . . . . . . . . . . . . 40 1.8 Potential energy function parameters for the ET systems. . . . . . . . . . . . . . . . . 40 1.9 Parameters for the auxiliary restraining potential in Eqns. (1.51)-(2.78). . . . . . . . . 42 1.10 Parameters for the polynomial fit to the lowest adiabatic electronic state (Eq. (1.42)). 43 2.1 Side-lengths of the simulation cell for the atomistic representation of PCET. . . . . . 65 2.2 Dividing surfaces for the sequential ET reaction prior to PT in the system-bath models. 74 2.3 Coefficients for the electron potential (Eq. (2.34)) in the System bath model. . . . . . 89 2.4 Parameters for the system-bath model of PCET. . . . . . . . . . . . . . . . . . . . . . 90 2.5 Parameters for the auxiliary restraining potential in Eq. 2.72. . . . . . . . . . . . . . . 91 2.6 Parameters for the auxiliary restraining potential in Eq. 2.75. . . . . . . . . . . . . . . 91 2.7 Parameters for the auxiliary restraining potential in Eq. 2.76. . . . . . . . . . . . . . . 92 2.8 Parameters for the auxiliary restraining potential in Eq. 2.79. . . . . . . . . . . . . . . 93 3.1 Parameters for System A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.2 Parameters for System B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 1 Chapter 1 Direct simulation of proton-coupled electron transfer across multiple regimes 1.1 Introduction Proton-coupled electron transfer (PCET) reactions, in which both an electron and an associated proton undergo reactive transfer (Fig. 1.1(a)), play an important role in many chemical and biological processes. 1–4 Key examples include the tyrosine oxidation step of photosystem II 5,6 and the protonpumping mechanism of cytochrome c oxidase. 7,8 Depending on the chronology of the electron- and proton-transfer events and the magnitudes of the electronic and vibrational coupling, a variety of reactive processes can fall under the umbrella of PCET; 9–13 investigation of the dynamics that govern this full range of behavior provides significant experimental and theoretical challenges, and the characterization of transitions between different regimes of PCET remains incomplete. In this study, we extend the ring polymer molecular dynamics (RPMD) method to allow for the direct simulation of PCET reaction dynamics and to characterize condensed-phase PCET reaction mechanisms and thermal rates across a wide range of physically relevant regimes. PCET reactions are typically described (Fig. 1.1(b)) in terms of the following reactant, intermediate, and product species: 1,9,14–16 D—H + A (OU) D− + [H—A] + (OP) [D—H] + A− + (RU) D + H—A (RP) Here, D and A indicate the donor and acceptor molecules, respectively, and the labels O/R and U/P indicate the oxidation state (oxidized or reduced) and the protonation state (unprotonated or protonated) of the acceptor molecule. The reactions can be categorized among two groups, sequential and concerted PCET, depending on whether both the electron and proton transfer in a single chemical step (Fig. 1.1(b)). 9,14–16 Sequential PCET exhibits distinct electron-transfer (ET) 2 e! (a)! H+ De !D p !H!!A p !A e (b)! D ! + [ H!A ] k O p k O p- D!H + A k + k P e D + H!A CPET k k k R p U e U + [ D!H ] + A ! e- Figure 1.1: (a) Schematic illustration of a co-linear PCET reaction, where De /Dp and Ae /Ap are the respective donor and acceptor for the electron/proton. (b) Schematic illustration of sequential and concerted PCET reaction mechanisms, indicating the rate constants for the individual charge transfer processes. The sequential mechanism proceeds along the horizontal and vertical edges of the schematic, whereas the concerted mechanism proceeds along the diagonal. and proton-transfer (PT) reaction events separated by a metastable intermediate species; concerted PCET exhibits the transfer of both particles in a single reactive step, bypassing the formation of the OP and RU species in Fig. 1.1(b). Within these two broad categories for PCET, there exist a range of coupling regimes that depend on the degree of electronic and vibrational non-adiabaticity for the PCET reaction. 9–13 Rate theories have been derived and successfully employed to study concerted PCET reactions in a variety of limiting regimes, including (i) the fully non-adiabatic regime, 1,17–19 in which the reaction is electronically and vibrationally non-adiabatic, (ii) the partially adiabatic regime, 10,12,20,21 in which the reaction is electronically adiabatic and vibrationally non-adiabatic, and (iii) the fully adiabatic regime, 20,21 in which the reaction is both electronically and vibrationally adiabatic. These rate theories, which generally employ Golden Rule and linear response approximations, provide a powerful toolkit for investigating both concerted and sequential PCET reactions in many systems. However, the applicability of any given rate theory is limited to the particular coupling regime for which it was derived, and without prior mechanistic information about a given PCET reaction, it can be difficult to know which formulation to apply in practice. Furthermore, with few exceptions, 12 existing rate theories do not offer scope for the study of PCET reactions with intermediate values for the electronic and vibrational coupling, which exist between the limiting regimes for which the rate 3 theories have been derived. Methods that enable the direct simulation of PCET reactions across all electronic and vibrational coupling regimes, including intermediate regimes, are needed to achieve a unified picture for the dynamics, mechanisms, and driving forces that govern the full range of PCET reactions. Fundamental theoretical challenges in the description of PCET reactions arise due to the coupling of intrinsically quantum mechanical ET and PT dynamics with slower motions of the surrounding environment. New simulation methods are needed to accurately describe this electronproton-environment dynamics and to efficiently and robustly simulate long trajectories that bridge the multiple timescales of these reactions. In this study, we address these challenges by extending the RPMD method to enable the direct simulation of condensed-phase PCET reactions. RPMD 22 is an approximate quantum dynamical method that is based on Feynman’s imaginary-time path integral formulation of statistical mechanics. 23,24 It provides a classical molecular dynamics model for the real-time evolution of a quantum mechanical system that rigorously preserves detailed balance and samples the quantum Boltzmann distribution. 24–26 The RPMD method has been previously employed to investigate a wide range of quantized reactive and dynamical processes, 27–40 ranging from gas-phase triatomic reactions 27 to enzyme-catalyzed hydrogen tunneling. 28 We have demonstrated that RPMD simulations can be extended to accurately and efficiently describe coupled electronic and nuclear dynamics in condensed-phase systems, including excess electron diffusion, 31 injection, 32 and reactive transfer. 33 Prior validation of RPMD for the description of ET reactions throughout the normal and activationless regimes, 33 in combination with prior demonstration of the method for a range of H-transfer processes, 27–30 provides a basis for expecting the method to adequately describe the dynamics of PCET reactions, which will be tested in the current study. Alternative theoretical methods have previously advanced our ability to simulate and understand coupled electronic and nuclear dynamics, 41–52 and promising new methods continue to be introduced. 53 Established methods include Ehrenfest dynamics, 41,42 mixed quantum-classical trajectory surface hopping dynamics, 43–48 the ab initio multiple spawning approach, 49 and semiclassical methods based on the Meyer-Miller-Stock-Thoss mapping. 50–52 However, despite their successes, these methods do not yield a dynamics that rigorously preserves detailed balance, 54,55 a feature that is valuable for the robust calculation of thermal reaction rates 56,57 and for the utilization of rare-event sampling methods. 58,59 Although it is clear that other methods must be part of the toolkit for understanding PCET reactions, we emphasize that the formal properties of the RPMD method are ideally suited to this goal. 4 In this paper, we extend the RPMD method to allow for direct simulation of co-linear, condensedphase PCET reactions across a wide range of physically relevant regimes. In addition to providing validation for the simulation method via extensive comparison with existing PCET rate theories, we analyze the RPMD reactive trajectories to elucidate a variety of mechanisms for the concerted charge-transfer process. The presented analysis offers a unifying picture for PCET across a wide range of physical regimes, and it suggests new PCET regimes that have yet to be characterized. 1.2 Ring polymer molecular dynamics The RPMD equations of motion for N particles that are quantized using n ring-polymer beads are 22,31 (α) v̇i (α) where vi (0) and qi (α) and qi   (α+1) (α−1) (α) = ωn2 qi + qi − 2qi   1 ∂ (α) (α) (α) − U q , q , . . . , q , 1 2 N mi ∂q (α) i (1.1) are the velocity and position of the αth bead for the ith particle, respectively, (n) = qi . The physical mass for particle i is given by mi , ωn = n/(β~) is the intra-bead harmonic frequency, and β = (kB T )−1 is the reciprocal temperature. The potential energy function of the system is given by U (q1 , . . . , qN ). To allow for the straightforward comparison with PCET rate theories, we quantize only the transferring electron and proton in this study and consider the classical (i.e., 1-bead) limit for the N solvent degrees of freedom. Furthermore, we employ a mixed-bead-number path-integral representation that reduces the cost of the potential energy surface calculations by utilizing the more rapid convergence of the path-integral distribution for heavier particles. 60 We thus obtain the modified RPMD equations of motion: v̇e(α)   = ωn2 e qe(α+1) + qe(α−1) − 2qe(α) − v̇p(γ) ! “ ” (α−k) n1ep +1 1 ∂ (α) U qe , qp ,Q , me ∂qe(α) (1.2)   = ωn2 p qp(γ+1) + qp(γ−1) − 2qp(γ) nep −   1 X ∂ U qe((γ−1)nep +l) , qp(γ) , Q , (γ) mp ∂qp l=1 (1.3) 5 and V̇j = − np nep   1 XX ∂ U qe((γ−1)nep +l) , qp(γ) , Q , ne Mj γ=1 ∂Qj (1.4) l=1 where ne is the number of imaginary-time ring-polymer beads for the transferring electron, me is the physical mass for the electron, and qe(α) and ve(α) are the respective position and velocity for the αth ring-polymer bead of the electron; the corresponding quantities for the transferring proton are indicated using subscript “p”. In Eqs. (2.4)-(1.4), it is assumed that nep = ne /np is an integer number, and  k = α − nep  α−1 , nep (1.5) where b. . . c denotes the floor function. As before, the periodic constraint of the ring-polymer is satisfied via qe(0) = qe(ne ) and qp(0) = qp(np ) , and the intra-bead harmonic frequencies are ωne = ne /(β~) and ωnp = np /(β~). The position, velocity, and mass for the jth classical solvent degree of freedom are given by Qj , Vj , and Mj , respectively, and Q = {Q1 , ..., QN }. In the limit of a large number of ring-polymer beads, the RPMD equations of motion yield a time-reversible molecular dynamics that preserves the exact quantum mechanical Boltzmann distribution. 24–26 Eqs. (2.4)-(1.4) introduce no approximation to Eq. (1.1) beyond taking the classical limit of the solvent degrees of freedom. Analogous to the classical thermal rate constant, 61–63 the RPMD thermal rate constant can be expressed as 56,57 kRPMD = lim κ(t)kTST , (1.6) t→∞ where kTST is the transition state theory (TST) estimate for the rate associated with the dividing surface ξ(r) = ξ ‡ , ξ(r) is a collective variable that distinguishes between the reactant and product basins of stability, and κ(t) is the time-dependent transmission coefficient that accounts for recrossing of tran o jectories through the dividing surface. We have introduced r = qe(1) , ..., qe(ne ) , qp(1) , ...qp(np ) , Q1 , ...QN to denote the position vector for the full system in the ring-polymer representation. As is the case for both exact classical and exact quantum dynamics, the RPMD method yields reaction rates and mechanisms that are independent of the choice of dividing surface. 56,57,64 The TST rate in Eq. (1.6) is calculated using 29,33,65,66 ‡ kTST = (2πβ) −1/2 e−β∆F (ξ ) , hgξ ic R ξ‡ dξe−β∆F (ξ) −∞ (1.7) 6 where F (ξ) is the free energy (FE) along ξ, e−β∆F (ξ) = hδ (ξ(r) − ξ)i , hδ (ξ(r) − ξr )i (1.8) ξr is a reference point in the reactant basin, and 29,67–69 gξ (r) = #1/2 " d X 1  ∂ξ(r) 2 i=1 mi ∂ri . (1.9) Here, ri is an element of the position vector r, mi is the corresponding physical mass, and d is the length of vector r. The equilibrium ensemble average is denoted R dr dv e−βH(r,v) (. . .) R R , dr dv e−βH(r,v) R h. . . i = (1.10) and the average over the ensemble constrained to the dividing surface is denoted R dr R h. . . ic = dv e−βH(r,v) (. . . )δ(ξ(r) − ξ ‡ ) R , dr dv e−βH(r,v) δ(ξ(r) − ξ ‡ ) R (1.11) where H(r, v) = N X 1 j=1 np + 2 X1 γ=1 2 Mj Vj2 + ne X 1 α=1 2  2 mb,e ve(α)  2 mb,p vp(γ) + URP (r). (1.12) Here, mb,e and mb,p are the fictitious Parrinello-Rahman masses for the electron and proton, ren o spectively, 25 v = ve(1) , ..., ve(ne ) , vp(1) , ..., vp(np ) , V1 , ..., VN is the velocity vector for the full system in the ring-polymer representation, and n URP (r) = e  2 1 1 X me ωn2 e qe(α) − qe(α−1) ne α=1 2 np  2 1 X1 + mp ωn2 p qp(γ) − qp(γ−1) np γ=1 2 np nep 1 X X  ((γ−1)nep +l) (γ)  + U qe , qp , Q . ne γ=1 l=1 (1.13) 7 The transmission coefficient in Eq. (1.7) is obtained from the flux-side correlation function, 56,57 κ(t) = hξ˙0 h(ξ(rt ) − ξ ‡ )ic , hξ˙0 h(ξ˙0 )ic (1.14) by releasing RPMD trajectories from the equilibrium ensemble constrained to the dividing surface. Here, h(ξ) is the Heaviside function, ξ˙0 is the time-derivative of the collective variable upon initialization of the RPMD trajectory from the dividing surface with the initial velocities sampled from the Maxwell-Boltzmann (MB) distribution, and rt is the time-evolved position of the system along the RPMD trajectory. 1.3 PCET Rate Theories A primary focus of this study is to compare the RPMD method with rate theories that have been derived for the various limiting regimes of PCET. We thus summarize these PCET rate theories below. 1.3.1 Concerted PCET in the Fully Adiabatic Regime For the fully adiabatic regime, both the electronic and vibrational coupling between the concerted PCET reactant and product states are large in comparison to the thermal energy, kB T . The reaction proceeds in the ground vibronic state, and it is appropriately described using the expression of Hynes and coworkers 20,21 ad kCPET = " # −∆G‡ad ωs exp , 2π kB T (1.15) where ωs is the solvent frequency, ∆G‡ad is the free-energy barrier for the reaction calculated from the difference of the ground vibronic energy level at its minimum and at its maximum with respect to the solvent coordinate, and kB is Boltzmann’s constant. 1.3.2 Concerted PCET in the Partially Adiabatic Regime For the partially adiabatic regime, the electronic coupling is large in comparison to kB T , whereas the vibrational coupling is small in comparison to kB T . The reaction proceeds in the ground electronic state, and it is appropriately described using the expression of Cukier 10 and Hynes and 8 co-workers, 20,21 pad kCPET =   −∆G‡ 2π 2 −1/2 Vµν (4πλkB T ) , exp ~ kB T (1.16) where λ is the concerted PCET reorganization energy associated with the transfer of both the electron and proton, ∆G ‡ = λ + ∆G0 4λkB T
2 , (1.17) ∆G0 is the driving force for the concerted PCET reaction, and Vµν is the vibronic coupling. In this regime, the vibronic coupling is equal to the vibrational coupling, VPT , such that Vµν = VPT = E1 − E0 . 2 (1.18) VPT is obtained from the splitting between the vibrational ground state energy, E0 , and first excited state energy, E1 , calculated on the lowest electronic adiabat. Eq. (1.16) assumes that only a single initial and final vibrational state are involved in the concerted PCET reaction. 10,20,21 1.3.3 Concerted PCET in the Fully Non-Adiabatic Regime For the fully non-adiabatic regime, both the electronic and vibrational coupling are small in comparison to kB T . The reaction is appropriately described using the expression of Cukier 17,18 and Hammes-Schiffer and coworkers 1,19 " # ‡ X X −∆G 2π µν −1/2 nad 2 kCPET = Pµ Vµν (4πλkB T ) , exp ~ µ kB T ν (1.19) where µ and ν index the reactant and product vibrational states, respectively, Pµ is the Boltzmann probability of the reactant vibrational state, and ∆G‡µν = λ + ∆G0 + ν − µ 4λkB T
2 , (1.20) 9 where µ and ν are the respective energies of the reactant and product vibrational states relative to their corresponding ground state. In this regime, the vibronic coupling is given by Vµν = hµ|νiVET , (1.21) where hµ|νi is the overlap between reactant and product vibrational wavefunctions, and VET is the electronic coupling. 1.3.4 ET rate theories We also compare RPMD simulations with rate theories that correspond to the electronically adiabatic and non-adiabatic regimes for pure ET. These ET rate theories are summarized below. 1.3.4.1 Adiabatic ET For the electronically adiabatic regime, the electronic coupling between the reactant and product ET states is large in comparison to kB T . The reaction proceeds in the ground electronic state, and it is appropriately described using Eq. (1.15), except with the free-energy barrier, ∆G‡ad , calculated from the difference of the ground electronic energy level at its minimum and at its maximum with respect to the solvent coordinate. 70,71 1.3.4.2 Non-adiabatic ET For the electronically non-adiabatic regime, the electronic coupling is small in comparison to kB T . The reaction is appropriately described using the standard Marcus theory expression, 72–74 2π 2 −1/2 nad kET |VET | (4πλkb T ) exp = ~ where ∆G‡µν = λ + ∆G0 4λkB T   −∆G‡ , kB T (1.22)
2 . (1.23) Here, VET , λ, and ∆G◦ are respectively the electronic coupling, reorganization energy, and driving force associated with the ET reaction. 10 1.4 PCET Model Systems Throughout this paper, condensed-phase PCET is described using a co-linear system-bath model. The model is expressed in the position representation using the potential energy function U (qe , qp , qs , Q) = Usys (qe , qp , qs ) + UB (qs , Q), (1.24) where UB (qs , Q) is the potential energy term associated with the bath coordinates, and Usys (qe , qp , qs ) = Ue (qe ) + Up (qp ) + Us (qs ) +Ues (qe , qs ) + Ups (qp , qs ) +Uep (qe , qp ) (1.25) is the system potential energy. The scalar coordinates qe , qp , and qs describe the positions of the electron, proton and solvent modes, respectively, and Q is the vector of bath oscillator positions. The first term in the system potential energy function models the interaction of the transferring electron with its donor and acceptor sites,  out in   aD qe2 + bD qe + cD , rD ≤ qe ≤ rD           in out aA qe2 + bA qe + cA , rA ≤ qe ≤ rA Ue (qe ) =            1 1    −µe + , otherwise, |qe − rD | |qe − rA | (1.26) where rD and rA are the positions of the electron donor and acceptor sites. This one-dimensional (1D) potential energy function consists of two symmetric coulombic wells, each of which is capped by quadratic functions to remove singularities. The second term in the system potential energy function models the interaction between the transferring proton and its donor and acceptor sites, Up (qp ) = − mp ωp2 2 m2p ωp4 4 qp + q . 2 16V0 p (1.27) Here, ωp is the proton vibrational frequency and V0 is the intrinsic PT barrier height. The next three terms in the system potential energy function model the solvent potential and 11 the electron- and proton-solvent interactions. Specifically, 1 ms ωs2 qs2 , 2 (1.28) Ues (qe , qs ) = −µes qe qs , (1.29) Ups (qp , qs ) = −µps qp qs , (1.30) Us (qs ) = and where ms is the solvent mass and ωs is the effective frequency of the solvent coordinate. The solvent coupling parameters, µes and µps , are of opposite sign due to the opposing charges of the transferring electron and proton. Interactions between the transferring electron and proton are modeled via the capped coulombic potential Ue (qe ) =  µep   − ,   |q  e − qp |  |qe − qp | > Rcut    µ    − ep , Rcut otherwise. (1.31) The potential energy term UB (qs , Q) models the harmonic bath that is coupled to the PCET reaction. The bath exhibits an ohmic spectral density J(ω) with cutoff frequency ωc , 75,76 such that J(ω) = ηωe−ω/ωc , (1.32) where η denotes the friction coefficient. The continuous spectral density is discretized into f oscillators with frequencies 56,77   j − 0.5 ωj = −ωc ln f (1.33) and coupling constants  cj = ωj such that 2ηM ωc fπ  f X  1 M ωj2 UB (qs , Q) = 2 j=1 1/2 , cj qs Qj − M ωj2 (1.34) !2  . (1.35) Here, M is the mass of each bath oscillator, and ωj and Qj are the respective frequency and position for the jth oscillator. We have developed system parameters to model condensed-phase PCET reactions throughout 12 a range of different physical regimes. Specifically, System 1 models the fully non-adiabatic regime, Systems 2a-2f model the transition between the fully non-adiabatic and partially adiabatic regimes, and Systems 3a-3e model the transition between the partially adiabatic and fully adiabatic regimes. Full details of the parameterization are provided in Appendices 1.1 and 1.2. We also employ a system-bath model to investigate pure ET in this study, with a potential energy function U ET (qe , qs , Q) = Ue (qe ) + Us (qs ) + Ues (qe , qs ) +UB (qs Q), (1.36) that is obtained by simply removing the proton-dependent terms in Eqs. (1.24) and (1.25). Systems 4a-4g model the transition between non-adiabatic and adiabatic ET. Full details of the parameterization for the ET reactions are provided in Appendices 1.1 and 1.2. 1.5 Calculation Details Calculations on System 1, Systems 2a-2f, and Systems 4a-4g are performed at T = 300 K; calculations on Systems 3a-3e are performed at the lower temperature of T = 100 K to clearly exhibit the transition between the partially adiabatic and fully adiabatic regimes for PCET. For all systems, the harmonic bath is discretized using f = 12 degrees of freedom. 1.5.1 RPMD Simulations In all simulations, the RPMD equations of motions are evolved using the velocity Verlet algorithm. 78 As in previous RPMD simulations, each timestep for the electron and proton involves separate coordinate updates due to forces arising from the physical potential and due to exact evolution of the purely harmonic portion of the ring-polymer potentials. 79 The electron is quantized with ne = 1024 ring-polymer beads in all systems, while the proton is quantized with np = 32 ringpolymer beads for Systems 1 and 2a-2f and with np = 128 for Systems 3a-3e. The larger number of beads for Systems 3a-3e is necessary due to the lower temperature. Two collective variables are used to monitor the PCET reaction mechanism in the RPMD simulations. The progress of the electron is characterized by a “bead-count” coordinate, fb , that reports 13 on the fraction of ring-polymer beads that are located on the electron donor, ne     1 X fb qe(1) , ..., qe(ne ) = tanh φqe(α) , ne α=1 (1.37) where φ = −3.0/rD . The progress of the proton is characterized using the ring-polymer centroid in the proton position coordinate, np   1 X (γ) q̄p qp(1) , ..., qp(np ) = q . np γ=1 p 1.5.1.1 (1.38) RPMD rate calculations for concerted PCET The RPMD reaction rate is calculated from the product of the TST rate and the transmission coefficient (Eq. (1.6)). The FE profiles that appear in the TST rate expression (Eq. (1.7)) are obtained using umbrella sampling and the weighted histogram analysis method (WHAM), as described below. 80,81 For System 2f, the 1D FE profile used in the rate calculation is obtained in the proton centroid coordinate, F (q̄p ), using the following umbrella sampling protocol. Nine independent sampling trajectories are harmonically restrained to uniformly spaced values of q̄p in the region [−0.20 a0 , 0.20 a0 ] using a force constant of 1.3 a.u. Additionally, 18 independent sampling trajectories are harmonically restrained to uniformly spaced values of q̄p in both the region [−1.10 a0 , −0.25 a0 ] and in [0.25 a0 , 1.10 a0 ] using a lower force constant of 1.0 a.u. to ensure extensive overlap among the sampled distributions. The equilibrium sampling trajectories are performed using path-integral molecular dynamics (PIMD) with mb,e = 2000 a.u. and mb,p = 1836.1 a.u., which allows for a timestep of 0.1 fs. Each sampling trajectory is run for 10 ns, and thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. We note that this choice of the Parrinello-Rahman masses, mb,e and mb,p , allows for a large timestep in the sampling trajectories but has no affect on F (q̄p ) or any other equilibrium ensemble average. 25,26 For all PCET systems other than System 2f, the 1D FE profile used in the rate calculation is obtained in the electron bead-count coordinate, F (fb ), using the following umbrella sampling protocol. 93 independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [−0.92, 0.92] using a force constant of 20 a.u.; seven independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in both the region [−0.991, −0.985] and in [0.985, 0.991] using a higher force constant of 5000 a.u.; nine independent sampling trajectories 14 are harmonically restrained to uniformly spaced values of fb in both the region [−1.0, −0.992] and in [0.992, 1.0] using a higher force constant of 10000 a.u.; 32 independent sampling trajectories are harmonically restrained to the values of fb ∈ { ±0.93, ±0.935, ±0.94, ±0.945, ±0.95, ±0.955, ±0.96, ±0.962, ±0.965, ±0.967, ±0.97, ±0.974, ±0.976, ±0.978, ±0.98, ±0.982} using a force constant of 500 a.u. For Systems 1 and 2a-2e, an auxiliary restraining potential is introduced for the PIMD sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.3. Each sampling trajectory is run for 10 ns using a timestep of 0.1 fs, with mb,e = 2000 a.u. and mb,p = 1836.1 a.u. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. For System 2f, the transmission coefficient (Eq. (1.14)) is calculated using RPMD trajectories that are released from the dividing surface associated with q̄p = 0. A total of 6000 RPMD trajectories are released. Each RPMD trajectory is evolved for 400 fs using a timestep of 1×10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u., and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. 82 For all PCET systems other than System 2f, the transmission coefficient is calculated using RPMD trajectories that are released from the dividing surface associated with fb = 0. A total of 4500 RPMD trajectories are released for Systems 1 and 2a-2e, and at least 10000 trajectories are released for Systems 3a-3e. Each RPMD trajectory is evolved for 300 fs using a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. For Systems 1 and 2a-2e, the same auxiliary restraining potential used in the calculation of F (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.3; throughout this paper, the RPMD trajectories used to calculate the transmission coefficients are not subjected to auxiliary restraining potentials. 15 1.5.1.2 RPMD rate calculations for ET prior to PT For System 1, we calculate the rate for both the sequential and concerted PCET mechanisms. For the ET step in the sequential mechanism, we calculate the forward and reverse ET reaction rates between the OU and RU species (keU and keU− , Fig. 1.1(b)). The symmetry of the system requires that keP = keU− . The 1D FE profile used in the rate calculation for the ET reactions is obtained in the electron bead-count coordinate, FSET (fb ), using the same umbrella sampling protocol described for the calculation of F (fb ); however, in the calculation of FSET (fb ), an auxiliary restraining potential is introduced for the PIMD sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.3. The independent sampling trajectories used to calculate FSET (fb ) are each run for 15 ns. The transmission coefficients (Eq. (1.14)) for the forward and reverse ET reactions are calculated using RPMD trajectories that are released from the dividing surface associated with fb = 0.18. A total of 12000 RPMD trajectories are released. Each RPMD trajectory is evolved for 300 fs using a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. The same auxiliary restraining potential used in the calculation of FSET (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.3. 1.5.1.3 RPMD rate calculations for PT prior to ET For the PT step in the sequential mechanism in System 1, we calculate the forward and reverse PT reactions between the OU and OP species (kpO and kpO− , Fig. 1.1(b)). The symmetry of the system requires that kpR = kpO− . The 1D FE profile used in the rate calculation for the forward and reverse PT reactions is obtained in the proton centroid coordinate, FSPT (q̄p ), using the same umbrella sampling protocol described for the calculation of F (q̄p ). The transmission coefficients (Eq. (1.14)) for the forward and reverse PT reactions are calculated using RPMD trajectories that are released from the dividing surface associated with q̄p = 0.21 a0 . A total of 10500 RPMD trajectories are released. Each RPMD trajectory is evolved for 300 fs with a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. 16 Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. 1.5.1.4 Two-dimensional FE profiles For the purpose of analysis, we calculate the two-dimensional (2D) FE profile for System 1 in the electron bead-count and proton centroid coordinates, F (fb , q̄p ). The 2D FE profile is constructed using PIMD sampling trajectories that are harmonically restrained in both the fb and q̄p coordinates. A total of 4553 sampling trajectories are performed, in which the coordinates fb and q̄p are sampled using a square grid. The coordinate fb is sampled using 93 windows that are harmonically restrained to uniformly spaced values of fb in the region [−0.92, 0.92] using a force constant of 20 a.u.; seven windows are harmonically restrained to uniformly spaced values of fb in both the region [−0.991, −0.985] and in [0.985, 0.991] using a higher force constant of 5000 a.u.; nine windows are harmonically restrained to uniformly spaced values of fb in both the region [−1.0, −0.992] and in [0.992, 1.0] using a higher force constant of 10000 a.u.; 32 windows are harmonically restrained to the values of fb ∈ { ±0.93, ±0.935, ±0.94, ±0.945, ±0.95, ±0.955, ±0.96, ±0.962, ±0.965, ±0.967, ±0.97, ±0.974, ±0.976, ±0.978, ±0.98, ±0.982} using a force constant of 500 a.u. For each value of fb , the coordinate q̄p is sampled using nine windows that are harmonically restrained to uniformly spaced values of q̄p in the region [−0.20 a0 , 0.20 a0 ] using a force constant of 1.3 a.u, and 10 windows that are harmonically restrained to uniformly spaced values of q̄p in both the region [−0.70 a0 , −0.25 a0 ] and in [0.25 a0 , 1.10 a0 ] using a lower force constant of 1.0 a.u. No auxiliary restraining potentials are employed for the calculation of F (fb , q̄p ). Each sampling trajectory is run for 2.5 ns using a timestep of 0.1 fs, with mb,e = 2000 a.u. and mb,p = 1836.1 a.u. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. We additionally calculate the 2D FE profile for System 1 in the electron bead-count and solvent position coordinates, F (fb , qs ), for sampling trajectories corresponding to the concerted PCET reaction. To generate F (fb , qs ), the harmonically restrained sampling trajectories used to calculate F (fb ) for System 1 are utilized. 17 1.5.1.5 RPMD transition path ensemble As we have done previously, 28 we analyze the transition path ensemble 58 for the RPMD trajectories in the current study. Reactive trajectories are generated through forward- and backward-integration of initial configurations drawn from the dividing surface ensemble with initial velocities drawn from the MB distribution. Reactive trajectories correspond to those for which forward- and backwardintegrated half-trajectories terminate on opposite sides of the dividing surface. The reactive trajectories that are initialized from the equilibrium Boltzmann distribution on the dividing surface must be reweighted to obtain the unbiased transition path ensemble. 58,83,84 A weighting term, wα , is applied to each trajectory, correctly accounting for recrossing and for the fact that individual trajectories are performed in the microcanonical ensemble. This term is given by 83 wα = X ˙ i ξ(r) −1 !−1 , (1.39) i ˙ i is where the sum includes all instances in which trajectory α crosses the dividing surface, and ξ(r) the velocity in the dividing surface collective variable at the ith crossing event. The reweighting has a minor effect on the non-equilibrium averages if the reactive trajectories initialized from the dividing surface exhibit relatively little recrossing, as is the case for the systems studied in this paper. Nonequilibrium averages over the RPMD transition path ensemble are calculated by aligning reactive trajectories at time 0, defined as the moment in time when the trajectories are released from the dividing surface. 1.5.1.6 RPMD rate calculations for pure ET The RPMD rates for pure ET are calculated for Systems 4a-4g. For Systems 4a-4e, the 1D FE profile used in the rate calculation is obtained in the electron bead-count coordinate, FET (fb ), using the same umbrella sampling protocol described for the calculation of F (fb ); however, no auxiliary restraining potentials are introduced for the PIMD sampling trajectories. For Systems 4f and 4g, the 1D FE profile used in the rate calculation is obtained in the solvent coordinate, FET (qs ), by reducing the 2D FE profile in the electron bead-count and solvent coordinates, FET (fb , qs ). The 2D FE profile, FET (fb , qs ), is constructed using PIMD sampling trajectories that are harmonically restrained in both the fb and qs coordinates. A total of 5809 sampling trajectories are performed, in which the coordinates fb and qs are sampled using a square grid. The coordinate fb is sampled using 93 windows that are harmonically restrained to uniformly spaced values of fb in 18 the region [−0.92, 0.92] using a force constant of 20 a.u.; seven windows are harmonically restrained to uniformly spaced values of fb in both the region [−0.991, −0.985] and in [0.985, 0.991] using a higher force constant of 5000 a.u.; nine windows are harmonically restrained to uniformly spaced values of fb in both the region [−1.0, −0.992] and in [0.992, 1.0] using a higher force constant of 10000 a.u.; 32 windows are harmonically restrained to the values of fb ∈ { ±0.93, ±0.935, ±0.94, ±0.945, ±0.95, ±0.955, ±0.96, ±0.962, ±0.965, ±0.967, ±0.97, ±0.974, ±0.976, ±0.978, ±0.98, ±0.982} using a force constant of 500 a.u. For each value of the fb coordinate, the qs coordinate is sampled using 37 windows that are harmonically restrained to uniformly spaced values of qs in the region [−9.0 a0 , 9.0 a0 ] using a force constant of 0.03 a.u. Each sampling trajectory is run for 2.5 ns using a timestep of 0.1 fs, with mb,e = 2000 a.u. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. For Systems 4a-4e, the transmission coefficient (Eq. (1.14)) is calculated using RPMD trajectories that are released from the dividing surface associated with fb = 0. A total of 3000 RPMD trajectories are released for Systems 4a-4c, 6000 trajectories for System 4d and 4500 trajectories for System 4e. Each RPMD trajectory is evolved for 300 fs using a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. For Systems 4f and 4g, the transmission coefficient is calculated using RPMD trajectories that are released from the dividing surface associated with qs = 0. A total of 1500 trajectories are released for Systems 4f and 4g. Each trajectory is evolved for 700 fs using a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. 1.5.2 PCET Rate Theory Calculations Expressions for the thermal reaction rates for concerted PCET are provided in Eqs. (1.15)-(1.21). Since the current paper considers only symmetric PCET reactions, the driving force, ∆G◦ , is zero 19 in all cases. The concerted PCET reorganization energy, λ, is calculated using the following result for symmetric systems, 85–87 λ = h∆U ireac , (1.40) where ∆U is the concerted PCET energy gap coordinate, ∆U = URP (−qe(1) , . . . , −qe(ne ) , −qp(1) , . . . , −qp(np ) , qs , Q) − URP (qe(1) , . . . , qe(ne ) , qp(1) , . . . , qp(np ) , qs , Q), (1.41) and h...ireac denotes the equilibrium ensemble average in the reactant basin. The ensemble average is calculated from a 50 ns equilibrium PIMD trajectory, where the electron and proton are initialized and remain in the reactant basin. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u., and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. Values for the reorganization energy in the various systems are presented in Table 1.1. The free-energy barrier for PCET in the fully adiabatic regime, ∆G‡ad in Eq. (1.15), is calculated from the difference of the ground vibronic energy level at its minimum and at its maximum with respect to the solvent coordinate. The adiabatic vibronic states are obtained as a function of the solvent coordinate in the range −4 a0 ≤ qs ≤ 4 a0 . For each value of qs , the system Hamiltonian associated with Usys (qe , qp , qs ) (Eq. (1.25)) is diagonalized using a 2D discrete variable representation (DVR) grid calculation in the electron and proton position coordinates, qe and qp , respectively. 88 The grid spans the range −30 a0 ≤ qe ≤ 30 a0 and −1.5 a0 ≤ qp ≤ 1.5 a0 , with 1024 and 20 evenly spaced grid points for the electron and proton position, respectively. The vibronic coupling in the partially adiabatic regime (Eq. (1.18)) is obtained from the splitting between the ground and first vibrational states calculated for the potential defined by the ground adiabatic electronic state; the ground adiabatic electronic state is calculated for a frozen solvent configuration for which the reactant and product concerted PCET states are degenerate. 20,21 The calculation of the vibronic coupling in the partially adiabatic regime thus requires two tasks that include (i) the calculation of the adiabatic electronic states as a function of the proton coordinate for a frozen solvent configuration and (ii) the calculation of the proton vibrational states for the potential defined by the lowest adiabatic electronic state. To complete task (i), the adiabatic electronic states are obtained as a function of the proton coordinate in the range −1.5 a0 ≤ qp ≤ 1.5 a0 , with qs = 0. 20 For each value of qp , the system Hamiltonian is diagonalized using a 1D DVR grid calculation in the electron position coordinate. The grid spans the range −30 a0 ≤ qe ≤ 30 a0 with 2048 evenly spaced grid points. To complete task (ii), a polynomial of the form Uad (qp ) = 6 X (i) cad |qp |i (1.42) i=0 is fit to the lowest adiabatic electronic state in the range −1.5 a0 ≤ qp ≤ 1.5 a0 . The vibrational energies, E0 and E1 , are calculated for the fitted potential in Eq. (1.42) by diagonalizing the 1D DVR Hamiltonian in the proton position coordinate. The grid spans the range −1.5 a0 ≤ qp ≤ 1.5 a0 with 2048 evenly spaced grid points. The values of the vibrational coupling, and hence the partially adiabatic vibronic coupling, are presented in Table 1.1. The coefficients for the polynomial fit to the lowest adiabatic electronic state (Eq. (1.42)) are presented in Appendix 1.4 (Table 1.10). The vibronic coupling in the fully non-adiabatic regime (Eq. (1.21)) is obtained from the product of the electronic coupling and the overlap of reactant and product vibrational wavefunctions. The vibrational wavefunctions are calculated for the potential defined by the reactant and product diabatic electronic states; the diabatic electronic states are calculated for a frozen solvent configuration for which the reactant and product concerted PCET states are degenerate. 1,17–19 The calculation of the vibronic coupling in the fully non-adiabatic regime thus requires three tasks that include (i) the calculation of the electronic coupling, (ii) the calculation of the diabatic electronic states as a function of the proton coordinate for a frozen solvent configuration, and (iii) the calculation of the vibrational energies and wavefunctions for the potential defined by the reactant and product diabatic electronic states. To complete tasks (i) and (ii) for Systems 1 and 2a-2f, the electronic coupling and diabatic electronic states are obtained as a function of the proton coordinate for qs = 0 using the localization procedure described in Appendix 1.5. The electronic coupling (Eq. (1.61)) is found to be nearly constant over the physical range of qp , so we employ a constant value of VET that corresponds to the qp = 0 value. For systems 2e and 2f, the localization procedure does yield fully localized diabatic states, which contributes to the breakdown of the fully non-adiabatic rate calculation. The values of the electronic coupling are presented in Table 1.1. To complete task (iii), the reactant and product diabatic electronic states (Eqs. (1.59) and (1.60)) are computed for a uniform grid of 2048 points in the range −1.5 a0 ≤ qp ≤ 1.5 a0 , and the reactant and product vibrational energies and wave functions are then obtained by diagonalizing the 1D DVR Hamiltonian in the proton position on this grid. 21 Table 1.1: Values of the electronic coupling, VET , vibrational coupling, VPT , and reorganization energy, λ, for the system-bath model systems for PCET. System 1 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3f a 1.5.3 VET VPT λ 5.0 × 10−6 5.0 × 10−6 5.0 × 10−5 5.0 × 10−4 5.0 × 10−3 2.5 × 10−2 1.0 × 10−1 3.3 × 10−2 2.7 × 10−2 1.8 × 10−2 1.5 × 10−2 1.5 × 10−2 1.8 × 10−6 4.6 × 10−7 3.6 × 10−7 2.1 × 10−7 9.6 × 10−8 4.5 × 10−7 4.8 × 10−6 1.5 × 10−8 2.3 × 10−5 8.5 × 10−4 2.2 × 10−3 2.8 × 10−3 1.84 × 10−2 9.71 × 10−3 9.47 × 10−3 9.78 × 10−3 9.22 × 10−3 9.32 × 10−3 8.47 × 10−3 3.27 × 10−2 3.29 × 10−2 3.34 × 10−2 3.35 × 10−2 3.33 × 10−2 Parameters are given in atomic units. ET Rate Theory Calculations Expressions for the thermal reaction rates for ET are provided in Eqs. (1.15) and (1.22). The free-energy barrier for ET in the electronically adiabatic regime, ∆G‡ad in Eq. (1.15), is calculated from the difference of the ground electronic energy level at its minimum and at its maximum with respect to the solvent coordinate. The adiabatic electronic states are obtained as a function of the solvent coordinate in the range −8.0 a0 ≤ qs ≤ 8.0 a0 . For each value of qs , the system Hamiltonian associated with Eq. (1.36) is diagonalized using a 1D DVR grid calculation in the electron position coordinate, qe . The grid spans the range −30.0 a0 ≤ qe ≤ 30.0 a0 with 2048 evenly spaced grid points. The electronic coupling, VET in Eq. (1.22), is obtained from the splitting between the ground, ε0 (qs ), and first excited, ε0 (qs ), adiabatic electronic state energies, VET = 1 [ε1 (qs = 0) − ε0 (qs = 0)] . 2 (1.43) The ET reorganization energy, λ, is calculated using Eq. (1.40), 85–87 where ∆U is now the ET energy gap coordinate, ∆U = ET URP (−qe(1) , . . . , −qe(ne ) , qs , Q) − ET (1) URP (qe , . . . , qe(ne ) , qs , Q). (1.44) 22 Table 1.2: Values of the electronic coupling, VET , and reorganization energy, λ, for ET systems that vary between the adiabatic and non-adiabatic regimes. System 4a 4b 4c 4d 4e 4f 4g a VET λ 1 × 10−6 1 × 10−5 1 × 10−4 1 × 10−3 4 × 10−3 1 × 10−2 2 × 10−2 7.18 7.45 7.44 7.37 7.26 7.18 7.30 λ is given in units of a.u.×10−2 ; all other parameters are given in atomic units. The ensemble average is calculated from a 50 ns equilibrium PIMD trajectory, where the electron is initialized and remains in the reactant basin. The sampling trajectories employ mb,e = 2000 a.u. and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The values of the reorganization energy are presented in Table 1.2. 1.6 Results The results are presented in two sections. In the first, we analyze the competition between the concerted and sequential reaction mechanisms for PCET. In the second, we study the kinetics and mechanistic features of concerted PCET reactions across multiple coupling regimes, including the fully non-adiabatic (both electronically and vibrationally non-adiabatic), partially adiabatic (electronically adiabatic, but vibrationally non-adiabatic), and fully adiabatic (both electronically and vibrationally adiabatic) limits. 1.6.1 Sequential versus concerted PCET We begin by investigating the competing PCET reaction mechanisms in System 1. Figure 1.2 presents the 2D FE profile for this system along the electron bead-count, fb , and the proton centroid, q̄p coordinates. The FE profile exhibits four basins of stability corresponding to the various PCET reactant (OU), intermediate (OP and RU), and product (RP) species (Fig. 1.1(b)). Distinct channels on the FE surface connect the various basins of stability. Due to the symmetry of the reaction, the two channels associated with the PT step of the sequential pathway (connecting OU to OP and RU to RP) are identical, as are the the two channels associated with the ET step of the sequential pathway (connecting OU to RU and OP to RP). A single channel on the FE surface connects OU to RP, bypassing the intermediate species. 23 25 OP RP -0.4 20 15 0 − qp / a0 0.4 OU RU 10 5 0 -0.8 -0.4 0 0.4 0.8 fb Figure 1.2: Reactive RPMD trajectories reveal distinct concerted (red), sequential PT-ET (purple), and sequential ET-PT (orange) reaction mechanisms for PCET in System 1. The trajectories are projected onto the FE surface in the electron bead-count coordinate, fb , and the proton centroid coordinate, q̄p , with contour lines indicating FE increments of 2 kcal/mol. Also plotted in Fig. 1.2 are representative samples from the ensemble of reactive RPMD trajectories for PCET in System 1. The trajectories cluster within the channels on the FE surface, providing a direct illustration of the concerted (red) and sequential (purple and orange) reaction mechanisms for PCET. Such distinct clustering of the reactive trajectories need not be observed in general systems that undergo PCET; we note that the RPMD method makes no a priori assumptions about the preferred reaction mechanism or the existence of distinct sequential and concerted reaction mechanisms for PCET. We now demonstrate that the concerted PCET mechanism is dominant in System 1 by computing the RPMD reaction rates for both the concerted and sequential processes. Figures 1.3(a) and 1.3(b) illustrate the FE profile and transmission coefficient that together determine the RPMD reaction rate for concerted PCET (Eq. (1.6)). As was previously found for ET reactions, 33 the FE profile exhibits a sharp rise as a function of fb due to the formation of a ring-polymer configuration in which the electron spans the two redox sites (Fig. 1.3(a), inset), and it exhibits more gradual changes in the range of |fb | < 0.97 due to solvent polarization. For the dividing surface fb = 0, the transmission coefficient plateaus at a value of approximately 0.1, indicating that fb is a reasonably good reaction coordinate for the process. These results combine to yield an RPMD rate of kCPET = (2.1 ± 0.7) × 10−20 a.u. for the concerted reaction mechanism in System 1. Figures 1.3(c)-(f) present the components of the RPMD rate calculation for the sequential PCET reaction mechanism in System 1. For the ET step of the sequential mechanism, Figs. 1.3(c) and 24 20 10 0.4 (a) 10 5 0 0 κ(t) F(fb) 30 0.98 1 0.2 0 -1 -0.5 0 0.5 1 fb 20 10 0 15 10 5 0 0.98 1 0 0 0 100 200 300 t / fs 1 (e) κ(t) − FSPT(qp) (d) 0.2 -1 -0.5 0 0.5 1 fb 10 100 200 300 t / fs 0.4 (c) κ(t) FSET(fb) 30 (b) (f) 0.5 0 -1 -0.5 0 0.5 1 − qp / a0 0 100 200 300 t / fs Figure 1.3: (a) The 1D FE profile in the electron bead-count coordinate, F (fb ), utilized in the RPMD rate calculation for the concerted PCET reaction. (b) The corresponding transmission coefficient for the concerted PCET reaction. (c) The 1D FE profile in the electron bead-count coordinate, FSET (fb ), utilized in the RPMD rate calculation for the ET reactions prior to PT in the sequential PCET mechanism. (d) The corresponding forward (red) and reverse (blue) transmission coefficients for the ET reactions prior to PT. (e) The 1D FE profile in the proton centroid coordinate, FSPT (q̄p ), utilized in the RPMD rate calculation for the PT reactions prior to ET in the sequential PCET mechanism. (f ) The corresponding forward (red) and reverse (blue) transmission coefficients for the PT reactions prior to ET. All FE profiles are plotted in kcal/mol. 25 Table 1.3: mechanism. RPMD rates for the forward and reverse ET and PT reactions in the sequential Rate Constant a keU (3.6 ± 2.6) × 10−21 keU− (1.1 ± 0.4) × 10−15 keP (1.1 ± 0.4) × 10−15 kpO (2.9 ± 0.3) × 10−13 kpO− (9.6 ± 1.7) × 10−8 kpR (9.6 ± 1.7) × 10−8 All rates are given in atomic units. The notation for the rate constants is defined in Fig. 1.1(b). 1.3(d) report the FE profile in the electron bead-count coordinate and the forward (red) and reverse (blue) transmission coefficients associated with fb = 0.18. For the PT step of the sequential mechanism, Figs. 1.3(e) and 1.3(f) report the FE profile in the proton centroid coordinate and the forward (red) and reverse (blue) transmission coefficients associated with qp = 0.21 a0 . The oscillations observed in κ(t) for the PT step correspond to the vibrational motion of the transferring proton. These results combine to yield the RPMD rates for the various individual steps in the sequential PCET reaction (Table 1.3). For the reaction mechanism that involves sequential ET followed by PT, the reaction rate is given by 1 kep = keU kpR , R kp + keU− (1.45) which numerically yields kep = (3.6 ± 2.6) × 10−21 a.u. Similiarly, for the reaction mechanism involving sequential PT followed by ET, the reaction rate is given by 1 kpe = kpO keP , keP + kpO− (1.46) which numerically yields kpe = (3.4 ± 1.4) × 10−21 a.u. The computed values for kep and kpe are equal to within statistical error, as is consistent with microscopic reversibility in this symmetric system. Comparison of the reaction rate for the concerted and sequential PCET mechanisms (Table 1.4) reveals that the reaction rate for the concerted mechanism is approximately six times larger than that of the sequential mechanism; the RPMD method thus predicts that the PCET reaction in System 1 proceeds predominantly via the concerted reaction mechanism. We note that although the reaction 26 Table 1.4: Reaction rates for the full ET-PT, PT-ET, and concerted PCET mechanisms calculated using RPMD and Eqs. (1.45) and (1.46). Rate Constant kep kpe kCPET a (3.6 ± 2.6) × 10−21 (3.4 ± 1.4) × 10−21 (2.1 ± 0.7) × 10−20 All rates are given in atomic units. dividing surfaces were selected to minimize trajectory recrossing, the rates reported here for the various sequential and concerted steps are rigorously independent of this choice of dividing surface; the mechanistic analysis provided here thus avoids any TST approximations or prior assumptions about the reaction mechanism. Having established that the concerted mechanism is favored for System 1, we now analyze the RPMD trajectories with respect to the solvent fluctuations and interactions that govern the concerted PCET reaction mechanism. Figure 1.4(a) presents the 2D FE profile in the electron bead-count and solvent coordinates, F (fb , qs ), computed for the concerted pathway as is described in Sec. 1.5.1.4. The FE profile exhibits two basins of stability corresponding to the PCET reactant and product species (OU and RP, respectively), separated by a barrier that corresponds to the dividing surface in the fb coordinate. Also plotted in Fig. 1.4(a) are representative samples from the ensemble of reactive RPMD trajectories (red) and the non-equilibrium average over the ensemble of reactive trajectories (yellow), as described in Sec. 1.5.1.5. As was seen for ET, 33 the reactive RPMD trajectories for concerted PCET follow a Marcus-type solvent-gating mechanism (black arrows), in which solvent reorganization precedes the sudden transfer of both the electron and proton between wells that are nearly degenerate with respect to solvent polarization. Figure 1.4(b) elaborates on this mechanism, schematically illustrating the ring-polymer configurations that accompany the various stages of the concerted PCET reaction. In the reactant OU basin, the system rests at the bottom of the solvent potential well for the reactant vibronic diabat (indicated by a red point in the left panel); for this polarized solvent configuration, the transferring electron and proton experience a potential energy surface that favors occupation of the donor sites (shown for the electron position in the right panel). In the dividing surface (DS) region of the concerted PCET reaction, the solvent fluctuation brings the system to configurations at which the vibronic diabats for the transferring electron and proton are nearly degenerate (shown at left), and the transferring particles undergo tunneling between nearly degenerate wells for the donor and 27 (a) qs / a0 2 OU DS 0 -2 RP -0.8 -0.4 0 fb 0.4 30 25 20 15 10 5 0 0.8 (b) Energy OU DS RP qs qe d(t) / a.u. (c) 0.005 0 -0.005 dp(t) dep(t) de(t) -300 -200 -100 0 100 200 300 t / fs Figure 1.4: (a) Reactive RPMD trajectories (red) and the average over the ensemble of reactive trajectories (yellow) for the concerted PCET reaction in System 1 reveal a Marcus-type solventgating mechanism indicated by the black arrows. The trajectories are projected onto the FE surface in the electron bead-count coordinate, fb , and the solvent position coordinate, qs , with contour lines indicating FE increments of 2 kcal/mol. The regions corresponding to the concerted PCET reactant (OU), product (RP), and dividing surface (DS) are indicated. (b) Illustration of the mechanism for concerted PCET. The left panels present the vibronic diabatic free energy surfaces along the solvent coordinate; the red dot indicates the solvent configuration associated with the OU, RP, and DS regions indicated in part (a). The right panels present the double-well potential that is experienced by the electron in the OU, RP, and DS regions, as well as the ring-polymer configuration in the electron position coordinate at the corresponding points along a typical reactive trajectory. (c) The combined dipole for the transferring particles in the ensemble of reactive RPMD trajectories, dep (t) (black), as well as the individual components from the transferring electron, de (red), and the transferring proton, dp (blue), for the concerted PCET reaction in System 1. 28 acceptor sites (shown at right); also seen in the panel at right is the extended “kink-pair” configuration for the ring-polymer in which the electron spans the two redox sites during the tunneling event. Finally, the figure panels associated with the product RP basin illustrate that as the solvent relaxes to the minimum of the solvent potential well for the product vibronic diabat (left), the transferring electron and proton experience a potential energy surface that favors occupation of the product sites (right). This mechanism observed in the RPMD trajectories is consistent with the mechanisms that are assumed by PCET rate theories in the fully non-adiabatic regime. 1,17–19 Figure 1.4(c) illustrates part of the mechanistic basis for the favorability of the concerted PCET reaction in this system. The figure presents the combined dipole for the transferring particles in the ensemble of reactive RPMD trajectories, dep (t), as well as the individual components from the transferring electron and proton, de (t) and dp (t), respectively. These terms are computed using de (t) = −µes hq̄e (t)itraj , (1.47) dp (t) = −µps hq̄p (t)itraj , (1.48) and dep (t) = de (t) + dp (t), (1.49) where q̄e and q̄p are the ring-polymer centroids for the transferring electron and proton, respectively, and h...itraj denotes the non-equilibrium ensemble average over the time-evolved reactive RPMD trajectories for concerted PCET (Sec. 1.5.1.5). Figure 1.4(c) shows that the orientation of de (t) and dp (t) switch during the reaction on similar timescales, which follows from the fact that the two particles are moving both co-linearly and in concert. However, the figure also shows that the magnitude of dep (t) is at all times smaller than the larger magnitude of the two component dipoles (i.e., |dep (t)| < max(|de (t)|, |dp (t)|)), due to the opposite charge of the two transferring particles. It is thus clear that throughout reactive trajectories for concerted PCET, the degree to which the polar solvent couples to the transferring particles is reduced by the opposing sign of the electron and proton charge. In this sense, the polar solvent creates a driving force for the co-localization of the electron and proton. We emphasize that although we have previously analyzed concerted versus sequential PCET mechanisms in the context of exact quantum simulations, 89 the RPMD simulations presented here constitute the first trajectory-based simulation approach to allow for the detailed, side-by-side comparison of the concerted and sequential PCET mechanisms and thermal reaction rates, with both 29 reaction mechanisms treated on a consistent dynamical footing. However, all of the results presented thus far have been obtained for System 1, which lies in the fully non-adiabatic regime for PCET. In the following section, we extend our analysis to the full range of coupling regimes for PCET. 1.6.2 Reactions across multiple coupling regimes In this section, we employ RPMD simulations to investigate concerted PCET in a range of physical regimes, including the fully non-adiabatic, partially adiabatic, and fully adiabatic regimes. We validate the accuracy of the RPMD method by comparing thermal reaction rates obtained using the simulation method with those obtained using previously developed rate theories, and we investigate the variety of electron and proton tunneling processes that accompany concerted PCET. However, before delving into this analysis of PCET reactions, we first use RPMD to examine the crossover between electronically non-adiabatic (i.e., weak electronic coupling) and electronically adiabatic (i.e., strong electronic coupling) regimes for pure ET; analysis of this more simple process will provide useful context for the subsequent discussion of PCET. 1.6.2.1 ET across electronic-coupling regimes Figure 1.5(a) presents the reaction rates for Systems 4a-4g, computed using RPMD (red), the electronically adiabatic ET rate expression (Eq. (1.15), blue), and the electronically non-adiabatic ET rate expression (Eq. (1.22), black). The results are plotted as a function of the temperature-reduced electronic coupling βVET . For the weak-coupling regime (βVET  1), the non-adiabatic rate expression constitutes the reference result, whereas for the strong-coupling regime (βVET  1), the adiabatic rate expression is the reference. It is clear that the RPMD rate correctly transitions between agreement with the non-adiabatic rate theory results at weak electronic coupling and the adiabatic rate theory results at strong electronic coupling. For systems with weak electronic coupling, we have shown previously that RPMD accurately describes the ET reaction rate throughout the normal and activationless regimes for the thermodynamic driving force, 33 which follows from the method’s exact description of statistical fluctuations 24–26 and its formal connection to semiclassical instanton theory for deep-tunneling processes. 24,90–93 Figure 1.5(a) shows that for symmetric systems, the accuracy of the method extends from the weak-coupling to the strong-coupling limits. It is important to note that the RPMD rates in Fig. 1.5(a) are obtained without prior knowledge or assumption of the electronic coupling regime, and at no point in the RPMD rate calculation is VET required. A natural question, therefore, is whether a posteriori analysis of the trajectories from 30 (b) Adiabatic -14 Non-Adiabatic RPMD (c) Kink Fraction qe 1 Single Kink-Pair 0.5 5 0 -5 Energy -19 (d) 5 0 -5 Energy -9 qe log(kET) / a.u. (a) -4 0 τ / βh 0 -3 -2 0.5 − Multiple Kink-Pairs -1 0 1 -5 0 5 qe 1 log(βVET) Figure 1.5: (a) ET reaction rates as a function of the temperature-reduced electronic coupling, obtained using RPMD (red), the non-adiabatic rate expression for ET (Eq. (1.22), black), and the adiabatic rate expression for ET (Eq. (1.15), blue) for Systems 4a-4g. (b) and (c) At left, the electron position as a function of the ring-polymer bead index for (b) System 4a (log (βVET ) = −2.98) and (c) System 4g (log (βVET ) = 1.32); at right, a schematic illustration of the corresponding doublewell potentials that are experienced by the transferring electron at the dividing surface, as well as the ring-polymer configurations in the electron position coordinate. The orange and purple stripes indicate the positions of the electron donor and acceptor sites, respectively. (d) The fraction of ring-polymer configurations at the dividing surface for ET that contain either a single kink-pair (black) or multiple kink-pairs (red) as a function of the temperature-reduced electronic coupling. 31 the RPMD rate calculation can be used to determine the electronic coupling regime for a given reaction. Figures 1.5(b)-(d) demonstrate that this is indeed the case. Figures 1.5(b) and 1.5(c) present snapshots of the electron position in a ring-polymer configuration at the reaction dividing surface, with the system either in the weak-coupling regime (part (b)) or in the strong coupling regime (part (c)). As described in the Calculation Details, the dividing surface used for the RPMD ET rate calculations is given by fb = 0, which corresponds to configurations for which the electron position evenly spans the two redox sites and for which the solvent is depolarized to accommodate this symmetric charge distribution for the electron. 33 At left in Figs. 1.5(b) and 1.5(c), we plot the electron position as a function of the ring-polymer bead index, α, where τ = β~α/ne . At right, we schematically illustrate the double-well potential that is experienced by the transferring electron at the dividing surface, as well as the ring-polymer configuration in the electron position coordinate. Note that for the weak-coupling regime (Fig. 1.5(b)), the configuration exhibits only a single kink-pair, in which the electron position transits between the redox sites as a function of the ring-polymer bead index; for the strong-coupling regime (Fig. 1.5(c)), the configuration exhibits multiple kink-pairs. It has long been recognized that the thermodynamic weight of ring-polymer kink-pair configurations is related to the eigenstate splitting (i.e., coupling) in symmetric double-well systems. 24,92–98 In particular, the weak-coupling regime corresponds to that for which the thermodynamic weight of ring-polymer configurations with multiple kink-pairs is small in comparison to the thermodynamic weight of ring-polymer configurations with only a single kink-pair; in the strong-coupling regime, configurations with multiple kink-pairs predominate. A straightforward approach to determining the coupling regime from the RPMD reactive trajectories is thus to simply count the fraction of ring-polymer configurations that exhibit multiple kink-pairs during the reactive transition event. For Systems 4a-4g, Fig. 1.5(d) presents the results of this strategy, in which RPMD results are used for the a posteriori determination of the regime of the electronic coupling. For each system, we calculate the fraction of ring-polymer configurations that exhibit either a single kink-pair (black) or multiple kink-pairs (red) in the ensemble from which the RPMD trajectories are initialized in the rate calculation (i.e., the equilibrium ensemble constrained to the dividing surface). Here, a kink is defined as a segment of the ring-polymer for which the electron position spans from the donor region (qe < −0.7σe ) to the acceptor region (qe > 0.7σe ), where σe is the standard deviation of the ringpolymer bead position in the dividing surface ensemble. We note that more sophisticated strategies for identifying the ring-polymer configurations in the transition region may be needed for systems 32 in which the trajectories exhibit extensive recrossing through a given dividing surface, 58,84 although that is not the case for the systems considered here. It is immediately clear from the comparison of Figs. 1.5(a) and 1.5(d) that the onset of multiple kink-pair configurations coincides with the crossover between the adiabatic and non-adiabatic regimes for pure ET reactions at βVET ≈ 1. We have thus shown that RPMD allows for the accurate calculation of the ET reaction rate across multiple regimes, without prior assumption of the electronic coupling regime, and it also enables determination of the coupling regime via simple analysis of the reactive trajectories. 1.6.2.2 Concerted PCET across electronic-coupling regimes We now shift our attention to Systems 2a-2f, which exhibit weak vibrational coupling and which vary in electronic coupling from the weak- to strong-coupling regimes. Figure 1.6(a) presents the thermal reaction rate for concerted PCET in these systems, calculated using the fully non-adiabatic rate theory (Eq. (1.19), black), the partially adiabatic rate theory (Eq. (1.16), blue), and the RPMD method (red). For the weak-coupling regime (βVET  1), the fully non-adiabatic rate expression constitutes the reference result, whereas for the strong-coupling regime (βVET  1), the partially adiabatic rate expression is the reference; the fully non-adiabatic results are discontinued (opencircle) at values of the electronic coupling for which the diabatic-state localization procedure becomes ill-defined (Sec. 1.5.2). As observed for the pure ET reactions, the RPMD method transitions correctly from the weak-coupling reference to the strong-coupling reference, while avoiding any assumptions about the coupling regime and while avoiding explicit calculation of the electronic or vibrational coupling. As for the pure ET reactions, we can analyze the ensemble of reactive RPMD trajectories for concerted PCET to elucidate the associated tunneling processes and to determine the electronic coupling regime for each system. Figures 1.6(b)-(d) present snapshots of a typical electron ringpolymer configuration at the concerted PCET reaction dividing surface, with the system either in the weak-coupling regime (part (b)), the intermediate-coupling regime (part (c)), or in the strong coupling regime (part (d)). In each case, the dividing surface corresponds to configurations for which the electron and proton positions are distributed between the donor and acceptor sites; for such configurations the solvent is depolarized to accommodate this symmetric charge distribution. The left panel in Figs. 1.6(b)-(d) presents the electron position as a function of the ring-polymer bead index; the right panels schematically illustrate the potential that is experienced by the transferring electron at the dividing surface, as well as the ring-polymer configurations in the electron position 33 Fully Non-Adiabatic (c) RPMD 0 -5 1 (d) 0.5 Multiple Kink-Pairs -2 -1 0 log(βVET) 0 -5 0 -3 5 Energy Single Kink-Pair qe Kink Fraction 5 Energy -15 -19 (e) 0 -5 qe log(kCPET) Partially Adiabatic -11 5 Energy (b) qe (a) -7 1 2 3 0 0.5 − τ / βh 1 -5 0 5 qe Figure 1.6: (a) Concerted PCET reaction rates as a function of the temperature-reduced electronic coupling, obtained using RPMD (red), the fully non-adiabatic rate expression (Eq. (1.19), black), and the partially adiabatic rate expression (Eq. (1.16), blue) for Systems 2a-2f. (b)-(d) At left, the electron position as a function of the ring-polymer bead index for (b) System 2a (log (βVET ) = −2.28), (c) System 2d (log (βVET ) = 0.72), and (d) System 2f (log (βVET ) = 2.02); at right, a schematic illustration of the corresponding potentials that are experienced by the transferring electron at the dividing surface, as well as the ring-polymer configurations in the electron position coordinate. The orange, purple, and green stripes indicate the positions of the electron donor site, the electron acceptor site, and transferring proton, respectively. (e) The fraction of ring-polymer configurations at the dividing surface for concerted PCET that contain either a single kink-pair (black) or multiple kink-pairs (red) as a function of the temperature-reduced electronic coupling. 34 coordinate. For the regime of weak electronic coupling (Fig. 1.6(b)), the electronic tunneling event that accompanies the PCET reaction is qualitatively similar to that observed for pure ET (Fig. 1.5(b)); the electron ring-polymer directly transitions between the two redox sites, exhibiting a single kinkpair. The coincident transfer of the proton in this regime simply affects the electron tunneling event by increasing the effective electronic coupling of the donor and acceptor redox sites, such that the concerted PCET mechanism may be described as proton-mediated electron superexchange. However, for the regime of strong electronic coupling (Fig. 1.6(d)), the electron transitions between the two redox sites via a mechanism that is fundamentally different than that observed for the pure ET reactions (Fig. 1.5(c)); in the PCET reaction, the electron collapses to a localized configuration about the position of the transferring proton, such that it adiabatically “rafts” with the proton between the donor and acceptor sites. This concerted PCET mechanism is immediately recognized as hydrogen atom transfer, or HAT. 10,99–101 In both limiting regimes for the electronic coupling (Figs. 1.6(b) and 1.6(d)), the RPMD trajectories reveal concerted PCET reaction mechanisms that are implicit in the associated PCET rate theories (Eqs. (1.19) and (1.16)). However, the RPMD simulations additionally reveal a distinct – and to our knowledge, previously undiscussed – mechanism for concerted PCET in the intermediate coupling regime, in which the tunneling electron partially localizes about three sites: the positions of the electron donor site, the electron acceptor site, and the proton that is simultaneously undergoing transfer (Fig. 1.6(c)). This intermediate mechanism, which might be called “transientproton-bridge” PCET, exhibits hybrid features of the PCET mechanisms from both limiting regimes (Figs. 1.6(b) and 1.6(d)), and it reflects the changing parameters that are employed to modulate the electronic coupling in Systems 2a-2f (Table 1.7); in this sense, it appears to be a physically reasonable mechanism for PCET in systems with intermediate electronic coupling, rather than an artifact of the approximate RPMD dynamics. Nonetheless, the transient-proton-bridge mechanism is certainly one for which no previous PCET rate theory has been been derived, and it remains to be seen whether an unambiguous kinetic signature of this new mechanism can be identified and observed in a physical system. Finally, Fig. 1.6(e) demonstrates that analysis of kink-pair formation in the reactive RPMD trajectories allows for the determination of the electronic coupling regime for the PCET reactions. As in Fig. 1.5(d), we present the calculated fraction of ring-polymer configurations that exhibit either a single kink-pair (black) or multiple kink-pairs (red) for the electron position in the equilibrium 35 path-integral ensemble constrained to the dividing surface. We employ the same definition for a ring-polymer kink as in the previous section, and we note that this definition registers multiple kinkpairs in the ring-polymer configurations in Figs. 1.6(c) and 1.6(d). Comparison of Figs. 1.6(a) and 1.6(e) demonstrates that the onset of multiple kink-pair configurations coincides with the crossover between the fully non-adiabatic and the partially adiabatic regimes for concerted PCET reactions. 1.6.2.3 Concerted PCET across vibrational-coupling regimes We now analyze Systems 3a-3e, which exhibit strong electronic coupling and which vary in vibrational coupling from the weak- to strong-coupling regimes. Figure 1.7(a) presents the thermal reaction rate for concerted PCET in these systems, calculated using the partially adiabatic rate theory (Eq. (1.16), blue), the fully adiabatic rate theory (Eq. (1.15), green), and RPMD (red). The results are plotted as a function of the temperature-reduced vibrational coupling βVPT . For the weak-coupling regime (βVPT  1), the partially adiabatic rate expression constitutes the reference result, whereas for the strong-coupling regime (βVPT  1), the fully adiabatic rate expression is the reference. The figure clearly shows that RPMD correctly transitions from the reference result in the regime of weak vibrational coupling to the reference for the strong-coupling regime. In the previous subsections we have demonstrated that analysis of the RPMD trajectories provides insight into the electron tunneling processes that accompany concerted PCET; here, we show that the same is true for the proton tunneling events. Figure 1.7(b) presents the expectation value for the radius of gyration in the proton coordinate (Eq. (1.50)) calculated either in the ensemble for the PCET reactant basin (blue) or the ensemble constrained to the PCET dividing surface (red). The radius of gyration in the proton coordinate is given by v u np  2 u 1 X (γ) Rg = t qp − q̄p . np γ=1 (1.50) For typical configurations in the dividing surface region, the insets in Fig. 1.7(b) schematically illustrate the potential that is experienced by the transferring proton and the lowest vibrational eigenstate for the proton, as well as the ring-polymer configurations in the proton position coordinate. For the weak-coupling regime, the ring polymer in the proton position “drapes” across the top of the potential barrier in the dividing surface configuration, leading to large values of Rg compared to the values in the reactant basin; for the strong-coupling regime, the proton exhibits more localized configurations. Using the connection between RPMD and semiclassical instanton 36 log(kCPET) / a.u. (a) -8 -12 Fully Adiabatic -16 RPMD -20 Partially Adiabatic -24 (b) 1 〈Rg〉 / a0 0.8 0.6 0.4 0.2 0 -4 -3 -2 -1 0 1 log(βVPT) Figure 1.7: (a) Concerted PCET reaction rates as a function of the temperature-reduced vibrational coupling, obtained using RPMD (red), the partially adiabatic rate expression (Eq. (1.16), blue), and the fully adiabatic rate expression (Eq. (1.15), green) for Systems 3a-3e. (b) The expectation value for the radius of gyration in the proton coordinate (Eq. (1.50)) calculated either in the ensemble for the PCET reactant basin (blue) or the ensemble constrained to the PCET dividing surface (red). The insets schematically llustrate the potential that is experienced by the transferring proton and the lowest vibrational eigenstate for the proton at the concerted PCET dividing surface, as well as the ring-polymer configurations in the proton position coordinate for System 3a (log (βVPT ) = −4.3, bottom-left) and System 3e (log (βVPT ) = 0.95, top-right). theory, 90,91 we recognize that a distended configuration of the ring polymer at the dividing surface indicates that the system is in the deep-tunneling regime. It is thus clear that the transition from strong vibrational coupling to weak vibrational coupling coincides with the onset of deep-tunneling for the proton coordinate. Comparison of Figs. 1.7(a) and 1.7(b) indicates that the onset of the deep-tunneling regime tracks the transition of the RPMD rates between the results calculated using the fully adiabatic and partially adiabatic rate expressions. 37 1.7 Conclusions In this work, we have extended the RPMD method to allow for the direct simulation of condensedphase PCET reactions across a wide range of physically relevant regimes. The results presented here provide validation of the method for the description of PCET reactions, as well as a unifying mechanistic picture for PCET in different regimes of electronic and vibrational coupling. We have demonstrated that the RPMD approach enables direct investigation of the competition between the concerted and sequential mechanisms for thermally activated PCET reactions, with both mechanisms treated on a consistent dynamical footing. For a system in the fully non-adiabatic regime, RPMD trajectories reveal distinct kinetic pathways associated with sequential and concerted PCET reaction mechanisms (Fig. 1.2); it was demonstrated that concerted PCET is favored in this system and that the concerted mechanism proceeds by a solvent-gating mechanism in which the reorganization energy is mitigated by charge cancelation among the transferring particles (Fig. 1.4). This analysis illustrates the potential usefulness of the RPMD method for the determination of PCET reaction mechanisms in systems that would otherwise be ambiguous on the basis of thermodynamic arguments alone, 3 such as the PCET reaction of bis(imidazole) iron tetraphenylporphyrins linked to an ascorbate derivative. 102 We have also validated the RPMD method for the PCET class of reactions by computing the thermal reaction rates for concerted PCET across a wide range of electronic and vibrational coupling regimes, including the fully non-adiabatic (electronically and vibrationally non-adiabatic), partially adiabatic (electronically adiabatic, but vibrationally non-adiabatic), and fully adiabatic (both electronically and vibrationally adiabatic) limits. Comparison of RPMD reaction rates with the results of previously derived PCET rate theories for each limiting regime demonstrates the accuracy of the direct simulation method in each case (Figs. 1.6(a) and 1.7(a)); it is particularly notable that RPMD correctly predicts the crossover in the thermal reaction rates between different coupling regimes, while avoiding the explicit calculation of the electronic or vibrational coupling and avoiding a priori assumptions about the coupling regime for the reaction. Moreover, by utilizing the connections between RPMD rate theory and semiclassical instanton theory, we showed that analysis of ring-polymer configurations in the RPMD transition path ensemble enables the a posteriori determination of the electronic and vibrational coupling regime for the reaction (Figs. 1.6(e) and 1.7(b)). These results indicate the potential utility of the RPMD method both for the calculation of PCET reaction rates in systems that span a broad range of physical regimes and for the determination of 38 the electronic and vibrational coupling regime in systems for which these attributes are unknown. Finally, we have used the RPMD method to identify and characterize a new PCET reaction mechanism that arises in the intermediate regime between fully-nonadiabatic and partially adiabatic PCET (Fig. 1.6(c)). In the limiting regimes of electronic and vibrational coupling, analysis of the reactive RPMD trajectories reveals previously anticipated electron and proton tunneling processes that accompany concerted PCET reactions, including the proton-mediated electron superexchange mechanism for the fully non-adiabatic regime (Fig. 1.6(b)) and the HAT mechanism for the partially adiabatic regime (Fig. 1.6(d)). However, the RPMD simulations also reveal a distinct “transientproton-bridge” mechanism in the crossover between these limiting regimes, in which the transferring electron partially localizes on the positions of the donor site, the acceptor site, and the proton that is simultaneously undergoing transfer (Fig. 1.6(c)). This progression of PCET mechanisms from superexchange to HAT as a function of electronic coupling suggests that concerted PCET may bear analogy to bridge-mediated ET reactions. 71,103–105 Furthermore, the identification of a new mechanism in the RPMD simulations illustrates the capacity of the direct simulation approach to reveal unanticipated, yet physically reasonable, reaction mechanisms for condensed-phase PCET. 39 Appendix A Coefficients for the Electron Potential This appendix presents the coefficients associated with the quadratic functions in Eq. (1.26). The coefficients are provided in Table 1.5 and are chosen such that Ue (qe ) is continuous and its derivative is piecewise continuous over the full range of qe . Table 1.5: Parameters for the donor coulombic well in the intrinsic electron potential energy function of Eq. (1.26). a System aD bD cD in rD out rD 1 2a 2b 2c 2d 2e 2f 3a-3e 4a 4b 4c 4d 4e 4f 4g 0.2266300 0.2266300 0.1824656 0.1523821 0.1379922 0.1298241 0.1067039 0.0743542 0.0932266 0.0929362 0.0923649 0.0910333 0.0891837 0.0867814 0.0834663 2.6933600 2.6933600 2.1684926 1.8109689 1.6399537 1.5428811 1.2681114 0.5740833 1.3904813 1.2049293 1.0121182 0.8059524 0.6684716 0.5657624 0.4753245 4.9465600 4.9465600 3.9826011 3.3259817 3.0118994 2.8336183 2.3289829 0.0207500 3.9598592 2.6659598 1.5130165 0.4947511 -0.0621548 -0.4151582 -0.6820558 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000 -4.0000 -2.0000 -5.4940 -4.5310 -3.5471 -2.5321 -1.8960 -1.4570 -1.1070 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -6.0000 -9.4940 -8.5310 -7.5471 -6.5321 -5.8960 -5.4570 -5.1070 in in The parameters for the acceptor coulombic well are given by aA = aD , bA = −bD , cA = cD , rA = −rD and out out rA = −rD . All parameters are given in atomic units. Appendix B System-Bath Potential Energy Parameters This appendix describes parameters for the system-bath models for condensed-phase PCET and ET. Table 1.6 presents parameters that are common to all PCET and ET model systems; the cutoff and solvent frequencies, ωc and ωs , are obtained from the spectral density of acetonitrile, 106 the masses of the solvent mode and bath oscillators, ms and M , correspond to that of a nitrogen atom, and the friction coefficient, η, is chosen to be consistent with a condensed-phase environment. 33,77 Table 1.7 presents potential energy function parameters for System 1, Systems 2a-2f, and Systems 3a-3e. Table 1.8 presents potential energy function parameters for Systems 4a-4g. 40 Table 1.6: Parameters for the PCET and ET systems. Parameter f ωs ωc M ms me mp η/M ωc a 12 2.3221×10−4 2.3221×10−4 25539 25539 1.0 1836.1 1.0 Parameters are given in atomic units, unless otherwise specified. Table 1.7: Potential energy function parameters for the PCET systems. a System rD rA µe ωp V0 µes µps µep 1 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3e -6.0 -6.0 -6.0 -6.0 -6.0 -6.0 -6.0 -4.0 -4.0 -4.0 -4.0 -4.0 6.0 6.0 6.0 6.0 6.0 6.0 6.0 4.0 4.0 4.0 4.0 4.0 3.6680 3.6680 2.9532 2.4663 2.2334 2.1012 1.7270 1.245 1.245 1.245 1.245 1.245 0.0131 0.011 0.011 0.011 0.011 0.011 0.011 0.009 0.009 0.009 0.009 0.009 23.7 20.0 20.0 20.0 20.0 20.0 20.0 20.0 10.0 5.0 3.0 2.3 -1.20 -1.00 -1.01 -1.06 -1.25 -2.01 -3.30 -1.55 -1.51 -1.43 -1.40 -1.40 6.00 4.95 4.95 4.95 4.95 4.95 1.00 0.0 0.0 0.0 0.0 0.0 1.3 1.3 1.3 1.3 1.3 1.3 1.3 0.15 0.15 0.20 0.30 0.35 V0 , µes and µps are given in units of a.u.×10−3 ; all other parameters are given in atomic units. Table 1.8: Potential energy function parameters for the ET systems. System 4a 4b 4c 4d 4e 4f 4g a rD rA µe µes -7.4940 -6.5310 -5.5471 -4.5321 -3.8960 -3.4570 -3.1070 7.4940 6.5310 5.5471 4.5321 3.8960 3.4570 3.1070 1.5 1.5 1.5 1.5 1.5 1.5 1.5 -0.94 -1.10 -1.30 -1.60 -1.89 -2.24 -2.80 µes is given in units of a.u.×10−3 ; all other parameters are given in atomic units. 41 Appendix C Auxiliary Restraining Potentials This appendix describes auxiliary restraining potentials that are introduced for the PIMD sampling trajectories used in the calculation of 1D FE profiles and in the initial sampling of configurations for the RPMD trajectories. These auxiliary restraining potentials simply prevent the PIMD sampling trajectories from visiting configurations outside of the sequential or concerted PCET channel of interest. 1.C.1 Auxiliary Restraining Potential for Concerted PCET For Systems 1 and 2a-2e, which exhibit distinct channels for sequential and concerted PCET reactions, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the concerted channel. This potential is given by Uaux (q̄p , fb ) =  2  aaux [q̄p − q+ (fb )] ,           2 aaux [q̄p − q− (fb )] ,            0, q̄p > q+ (fb ) q̄p < q− (fb ) (1.51) otherwise where q+ (fb ) = baux fb + caux (1.52) q− (fb ) = baux fb − caux . (1.53) and The coefficients aaux , baux , and caux (Table 1.9) are chosen to restrict the system to the concerted channel. 1.C.2 Auxiliary Restraining Potential for ET prior to PT For System 1, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the ET channel connecting the OU and RU species in the sequential 42 Table 1.9: Parameters for the auxiliary restraining potential in Eqns. (1.51)-(2.78). a System aaux baux caux 1 2a 2b 2c 2d 2e 1.0 1.0 1.0 1.0 1.0 1.0 0.57 0.66 0.67 0.69 0.78 1.11 0.36 0.39 0.40 0.41 0.41 0.38 aaux is given in units of a.u.×10−2 ; all other parameters are given in atomic units. mechanism. This potential is given by USET (q̄p ) =  2   aSET (q̄p − bSET ) ,    q̄p < bSET      0, otherwise. (1.54) The coefficients aSET = 1.0 × 102 a.u. and bSET = −0.21 a0 are chosen to correctly restrict the system to the ET channel. Appendix D Parameters for Fitting Functions This appendix presents the coefficients associated with the polynomials fit to the proton potential energy surfaces in Eq. (1.42). The coefficients are provided in Table 1.10. Appendix E Diabatization Protocol This appendix describes the procedure for transforming the potential energy function from a position basis for the electron (Eq. (1.24)) to a diabatic basis in which the reactant and product states are maximally localized on the position of the electron donor and acceptor, respectively. The diabatization protocol presented below is similar to that used for ET reactions in Ref. 33. We begin by calculating the two lowest adiabatic electronic eigenstates (ψ0 (qe ; qp , qs ) and ψ1 (qe ; qp , qs )) and eigenenergies (ε0 (qp , qs ) and ε1 (qp , qs )) of the system Hamiltonian associated with Usys (qe , qp , qs ) (Eq. (1.25)). The 1D eigenvalue problem is solved at fixed values of the proton coordinate in the range −1.5 a0 ≤ qp ≤ 1.5 a0 and with qs = 0. For each value of qp the system Hamiltonian is diagonalized on a uniform 1D DVR grid of 2048 electron positions in the range −30 a0 ≤ qe ≤ 30 a0 . For each value of qp , reactant and product electronic wavefunctions in the diabatic basis are 43 Table 1.10: Parameters for the polynomial fit to the lowest adiabatic electronic state (Eq. (1.42)). System cad (6) cad (5) cad 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3e System 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3e System 2a 2b 2c 2d 2e 2f 3a 3b 3c 3d 3e 2.340 × 10−5 3.054 × 10−7 1.776 × 10−5 8.203 × 10−6 1.886 × 10−3 4.547 × 10−5 1.122 × 10−5 9.645 × 10−6 1.987 × 10−5 2.553 × 10−4 2.064 × 10−4 (3) cad −2.182 × 10−3 −1.910 × 10−3 7.850 × 10−3 9.529 × 10−2 −1.197 × 10−2 1.201 × 10−3 5.452 × 10−2 4.522 × 10−2 2.438 × 10−2 1.226 × 10−2 1.251 × 10−2 (0) cad −2.948 −2.389 −2.013 −1.842 −1.782 −1.650 −1.037 −1.015 −9.758 × 10−1 −9.576 × 10−1 −9.576 × 10−1 −4.901 × 10−4 −4.783 × 10−4 1.565 × 10−3 1.441 × 10−2 −1.029 × 10−2 5.369 × 10−4 4.872 × 10−3 4.111 × 10−3 2.061 × 10−3 −8.242 × 10−5 1.192 × 10−4 (2) cad −1.191 × 10−1 −1.204 × 10−1 −1.292 × 10−1 −2.042 × 10−1 −1.607 × 10−1 −1.373 × 10−1 −1.493 × 10−1 −1.369 × 10−1 −1.113 × 10−1 −9.804 × 10−2 −9.812 × 10−2 1.545 × 10−1 1.538 × 10−1 1.452 × 10−1 9.150 × 10−2 1.737 × 10−1 1.508 × 10−1 4.441 × 10−2 1.175 × 10−1 2.657 × 10−1 4.571 × 10−1 5.970 × 10−1 (1) cad −4.128 × 10−2 −4.248 × 10−2 −4.301 × 10−2 −2.031 × 10−2 −2.762 × 10−4 4.793 × 10−5 8.250 × 10−4 5.791 × 10−4 4.775 × 10−4 3.088 × 10−4 3.207 × 10−4 a (4) Parameters are given in atomic units. obtained via rotation of the two lowest-energy adiabatic wavefunctions, using φR (qe ; qp , qs ) = cos (θp ) ψ0 (qe ; qp , qs ) − sin (θp ) ψ1 (qe ; qp , qs ) (1.55) φP (qe ; qp , qs ) = sin (θp ) ψ0 (qe ; qp , qs ) + cos (θp ) ψ1 (qe ; qp , qs ), (1.56) and 44 where θp = 1 arctan 2  S10 + S01 S11 − S00  (1.57) and Z 0 Suv = ∗ ψu (qe ; qp , qs ) ψv (qe ; qp , qs ) dqe . (1.58) −∞ Z 0 2 |φR (qe ; qp , qs )| dqe , the probability that the This choice of the rotation angle, θp , maximizes −∞ reactant diabatic state is positioned on the electron donor. Maximization of the probability that the product state is positioned on the electron acceptor yields an identical choice for θp . The corresponding reactant electronic diabat is thus given by VR (qp , qs ) = ε0 (qp , qs ) cos2 θp + ε1 (qp , qs ) sin2 θp , (1.59) and the product electronic diabat is given by VP (qp , qs ) = ε0 (qp , qs ) sin2 θp + ε1 (qp , qs ) cos2 θp . (1.60) The electronic coupling is given by VET (qp , qs ) = (ε0 (qp , qs ) − ε1 (qp , qs )) cos θp sin θp . (1.61) 45 References [1] S. Hammes-Schiffer and A. A. Stuchebrukhov, Chem. Rev. 110, 6939 (2010). [2] S. Y. Reece and D. G. Nocera, Annu. Rev. Biochem. 78, 673 (2009). [3] J. J. Warren, T. A. Tronic, and J. M. Mayer, Chem. Rev. 110, 6961 (2010). [4] J. L. Dempsey, J. R. Winkler, and H. B. Gray, Chem. Rev. 110, 7024 (2010). [5] M. Y. Okamura, M. L. Paddock, M. S. Graige, and G. Feher, Biochim. Biophys. Acta. Bioenerg. 1458, 148 (2000). [6] M. Saraste, Science 283, 1488 (1999). [7] M. Wikstrom, Curr. Opin. Struct. Biol. 8, 480 (1998). [8] H. Michel, Biochemistry 38, 15129 (1999). [9] S. Y. Reece, J. M. Hodgkiss, J. Stubbe, and D. G. Nocera, Phil. Trans. R. Soc. B 361, 1351 (2006). [10] R. I. Cukier, J. Phys. Chem. B 106, 1746 (2002). [11] J. Ulstrup, Charge transfer processes in condensed media (Springer-Verlag, Berlin, 1979). [12] Y. Georgievskii and A. A. Stuchebrukhov, J. Chem. Phys. 113, 10438 (2000). [13] S. Hammes-Schiffer, Energy Environ. Sci. 5, 7696 (2012). [14] R. I. Cukier and D. G. Nocera, Annu. Rev. Phys. Chem. 49, 337 (1998). [15] M. H. V. Huynh and T. J. Meyer, Chem. Rev. 107, 5004 (2007). [16] J. M. Mayer and I. J. Rhile, Biochim. Biophys. Acta. Bioenerg. 1655, 51 (2004). [17] R. I. Cukier, J. Phys. Chem. 98, 2377 (1994). [18] R. I. Cukier, J. Phys. Chem. 100, 15428 (1996). [19] A. Soudackov and S. Hammes-Schiffer, J. Chem. Phys. 113, 2385 (2000). 46 [20] A. Staib, D. Borgis, and J. T. Hynes, J. Chem. Phys. 102, 2487 (1994). [21] D. Borgis and J. T. Hynes, J. Phys. Chem. 100, 1118 (1996). [22] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 (2004). [23] R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals, International series in pure and applied physics (McGraw-Hill, New York,, 1965) p. xiv. [24] D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 (1981). [25] M. Parrinello and A. Rahman, J. Chem. Phys. 80, 860 (1984). [26] B. D. Raedt, M. Sprik, and M. L. Klein, J. Chem. Phys. 80, 5719 (1984). [27] R. Collepardo-Guevara, Y. V. Suleimanov, and D. E. Manolopoulos, J. Chem. Phys. 130, 174713 (2009). [28] N. Boekelheide, R. Salomón-Ferrer, and T. F. Miller III, Proc. Natl. Acad. Sci. 108, 16159 (2011). [29] R. Collepardo-Guevara, I. R. Craig, and D. E. Manolopoulos, J. Chem. Phys. 128, 144502 (2008). [30] Y. V. Suleimaov, R. Collepardo-Guevara, and D. E. Manolopoulos, J. Chem. Phys. 134, 044131 (2011). [31] T. F. Miller III, J. Chem. Phys. 129, 194502 (2008). [32] A. R. Menzeleev and T. F. Miller III, J. Chem. Phys. 132, 034106 (2010). [33] A. R. Menzeleev, N. Ananth, and T. F. Miller III, J. Chem. Phys. 135, 074106 (2011). [34] T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. 123, 154504 (2005). [35] T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. 122, 184503 (2005). [36] S. Habershon, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 131 (2009). [37] I. R. Craig and D. E. Manolopoulos, Chem. Phys. 322, 236 (2006). [38] T. E. Markland, S. Habershon, and D. E. Manolopoulos, J. Chem. Phys. 128, 194506 (2008). [39] Y. V. Suleimanov, J. Phys. Chem. C 116, 11141 (2012). 47 [40] T. E. Markland, J. A. Morrone, K. Miyazaki, B. J. Berne, D. R. Reichman, and E. Rabani, J. Chem. Phys. 136, 074511 (2012). [41] P. Ehrenfest, Z. Phys 45, 455 (1927). [42] A. D. Mclachlan, Mol. Phys. 8, 39 (1964). [43] J. C. Tully, J. Chem. Phys. 93, 1061 (1990). [44] J. C. Tully, Faraday Discuss. , 407 (1998). [45] J.-Y. Fang and S. Hammes-Schiffer, J. Chem. Phys. 106, 8442 (1997). [46] A. Hazra, A. V. Soudackov, and S. Hammes-Schiffer, J. Phys. Chem. B 114, 12319 (2010). [47] A. V. Soudackov, A. Hazra, and S. Hammes-Schiffer, J. Chem. Phys. 135, 144115 (2011). [48] B. Auer, A. V. Soudackov, and S. Hammes-Schiffer, J. Phys. Chem. B 116, 7695 (2012). [49] M. Ben-Nun, J. Quenneville, and T. J. Martinez, J. Phys. Chem. A 104, 5161 (2000). [50] H. D. Meyer and W. H. Miller, J. Chem. Phys. 70, 3214 (1979). [51] H. D. Meyer and W. H. Miller, J. Chem. Phys. 72, 2272 (1980). [52] G. Stock and M. Thoss, Phys. Rev. Lett. 78, 578 (1997). [53] P. Shuskhov, R. Li, and J. C. Tully, J. Chem. Phys. 137, 22A549 (2012). [54] S. Nielsen, R. Kapral, and G. Ciccotti, J. Chem. Phys. 115, 5805 (2001). [55] J. R. Schmidt, P. V. Parandekar, and J. C. Tully, J. Chem. Phys. 129 (2008). [56] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 122, 084106 (2005). [57] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 123, 34102 (2005). [58] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002). [59] L. Maragliano, A. Fischer, E. Vanden-Eijnden, and G. Ciccotti, J. Chem. Phys. 125, 024106 (2008). [60] R. P. Steele, J. Zwickl, P. Shushkov, and J. C. Tully, J. Chem. Phys. 134, 074112 (2011). 48 [61] E. Wigner, Phys. Chem. Abt. B 19, 203 (1932). [62] H. Eyring, J. Chem. Phys. 3, 107 (1935). [63] J. C. Keck, J. Chem. Phys. 32, 1035 (1960). [64] W. H. Miller, J. Chem. Phys. 58, 1664 (1973). [65] D. Chandler, J. Chem. Phys. 68, 2959 (1978). [66] C. H. Bennet, in Algorithms for Chemical Computations, edited by R. E. Christofferson (American Chemical Society, Washington, DC, 1977) p. 63. [67] E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989). [68] G. K. Schenter, B. C. Garrett, and D. G. Truhlar, J. Chem. Phys. 119, 5828 (2003). [69] J. B. Watney, A. V. Soudackov, K. F. Wong, and S. Hammes-Schiffer, Chem. Phys. Lett. 25, 268 (2006). [70] R. A. Marcus, Annu. Rev. Phys. Chem. 115, 155 (1964). [71] A. Nitzan, Chemical dynamics in condensed phases (Oxford university press, 2006). [72] R. A. Marcus, Disc. Faraday Soc. 29, 21 (1960). [73] R. A. Marcus, J. Chem. Phys. 43, 3477 (1965). [74] R. A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 265 (1985). [75] A. O. Caldeira and A. J. Leggett, Annals of Physics 149, 374 (1983). [76] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987). [77] M. Topaler and N. Makri, J. Chem. Phys. 101, 7500 (1994). [78] L. Verlet, Phys. Rev. 159, 98 (1967). [79] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, Annu. Rev. Phys. Chem. 64, 387 (2013). [80] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, and J. M. Rosenberg, J. Comput. Chem. 13, 1011 (1992). 49 [81] S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, J. Comput. Chem. 16, 1339 (1995). [82] H. C. Andersen, J. Comput. Phys. 52, 24 (1983). [83] G. Hummer, J. Chem. Phys. 120, 516 (2004). [84] E. W. Vanden-Eijnden, J. Stat. Phys. 123, 503 (2006). [85] A. Warshel, J. Phys. Chem. 86, 2218 (1982). [86] G. King and A. Warshel, J. Chem. Phys. 93, 8682 (1990). [87] D. W. Small, D. V. Matyushov, and G. A. Voth, J. Am. Chem. Soc. 125, 7470 (2003). [88] D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1991). [89] N. Ananth and T. F. Miller, III, Mol. Phys. 110, 1009 (2012). [90] J. O. Richardson and S. C. Althorpe, J. Chem. Phys. 131, 214106 (2009). [91] S. C. Althorpe, J. Chem. Phys. 134, 114104 (2011). [92] A. Kuki and P. G. Wolynes, Science 236, 1647 (1987). [93] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). [94] M. Marchi and D. Chandler, J. Chem. Phys. 95, 889 (1991). [95] J. O. Richardson and S. C. Althorpe, J. Chem. Phys. 134, 054109 (2011). [96] J. O. Richardson, S. C. Althorpe, and D. J. Wales, J. Chem. Phys. 135, 124109 (2011). [97] A. O. Caldeira and A. J. Leggett, Ann. Phys. 149, 374 (1982). [98] D. Chandler, in Liquids, Freezing and the Glass Transition, Les Houches 51, edited by D. Levesque, J. P. Hansens, and J. Zinn-Justin (Elsevier science publishers, Amsterdam, 1991) pp. 193–285. [99] S. Hammes-Schiffer, Chem. Phys. Chem. 3, 33 (2002). [100] J. H. Skone, A. V. Soudackov, and S. Hammes-Schiffer, J. Am. Chem. Soc. 128, 16655 (2006). [101] R. I. Cukier, Biochim. Biophys. Acta. Bioenerg. 1655, 37 (2004). 50 [102] J. J. Warren, Gray, A. R. Menzeleev, and J. M. Mayer, J. S. Kretchmer, T. F. Miller III, H. B. J. Phys. Chem. Lett. (2013), 10.1021/jz400029w, http://pubs.acs.org/doi/pdf/10.1021/jz400029w . [103] H. B. Gray and J. R. Winkler, Proc. Natl. Acad. Sci. 102, 3534 (2005). [104] M. Bixon and J. Jortner, J. Chem. Phys. 107, 5154 (1997). [105] M. D. Newton, Chem. Rev. 91, 767 (1991). [106] M. Cho, S. J. Rosenthal, N. F. Scherer, L. D. Ziegler, and G. R. Fleming, J. Chem. Phys. 96, 5033 (1991). 51 Chapter 2 Understanding and controlling concerted versus sequential protoncoupled electron transfer 2.1 Introduction Proton-coupled electron transfer (PCET), in which both an electron and a proton undergo reactive transfer, is ubiquitous throughout chemistry and biology. 1–4 PCET reactions are commonly categorized into two groups, concerted and sequential PCET, based on the chronology of the electron- and proton-transfer events. 5–8 The determination of whether a PCET reaction is concerted or sequential provides significant experimental and theoretical challenges, and an understanding of the physical driving forces that lead to one mechanism over the other remains incomplete. In this study, we utilize the ring polymer molecular dynamics (RPMD) method to investigate the competition between concerted and sequential PCET mechanisms in iron bi-imidazoline complexes. Iron(II) and iron(III) complexes of protonated and deprotonated bi-imidazoline ligands, Fig. 2.1a, were first studied as a model for non-heme iron containing enzymes involved in C-H bond activation; FeIII (Hbim) was shown to abstract a hydrogen atom from 9,10-dihydroantrhacene as a mimic of the lipoxygenase enzyme. 9 More recently, Mayer and co-workers have experimentally investigated the single electron transfer (ET), single proton transfer (PT), and PCET reactions 52 (a)! 3+ HN H N N N N N NH N FeII NH N N N H N H FeII(H2bim)! N N N- 2+ HN H N N N FeII N N- N NH FeIII N NH N N N H NH HN 1+ N N FeIII(H2bim)! HN N NH N HN H N N FeIII N N 3+ HN H N N N H HN FeII(Hbim)! HN FeIII(Hbim)! (b)! FeII(H2bim) + ! FeIII(Hbim)! ET! T! E CP ET! FeIII(Hbim) + ! FeII(H2bim)! PT! PT! FeII(Hbim) + ! FeIII(H2bim)! FeIII(H2bim) + ! FeII(Hbim)! Figure 2.1: (a) Lewis structures of the different iron bi-imidazoline complexes. (b) PCET scheme depicting the sequential and concerted PCET mechanisms in iron bi-imidazoline. The sequential mechanism proceeds along the horizontal and vertical edges of the schematic, whereas the concerted mechanism proceeds along the diagonal. 53 present in these iron bi-imidazoline complexes: 10  II 2+  III 3+ Fe (H2 bim)3 + Fe (H2 bim)3 2+  II 3+ ET  III − − * + Fe (H2 bim)3 ) − − Fe (H2 bim)3 (2.1) (2.2)  III 3+  III 2+ Fe (H2 bim)3 + Fe (Hbim)(H2 bim)2 2+  III 3+ PT  III − − * + Fe (H2 bim) ) − − Fe (Hbim)(H2 bim)2 (2.3)  II 2+  III 2+ Fe (H2 bim)3 + Fe (Hbim)(H2 bim)2 2+  II 2+ PCET  III − − * + Fe (H2 bim)3 . ) −− −− − − Fe (Hbim)(H2 bim)2 The rates measured at 298 K in acetonitrile were kET = (1.7 ± 0.2) × 104 , kPT ≈ 2 × 106 , and kPCET = 5.8 ± 0.6 × 103 M−1 s−1 for Eqs. (2.2), (2.3), and (2.4), respectively. Mayer and co-workers explained the interesting result that the rates for ET and PCET were similar as a balance between the larger outer-sphere reorganization for ET and a larger inner-sphere reorganization for PCET. 10 Later theoretical work, using golden-rule rate expressions, explained the similarity between the two rates as a balance between a larger outer-sphere reorganization energy for ET and a smaller vibronic coupling for PCET. 11 The PCET reaction in iron bi-imidazoline, Eq. 2.4, can follow either a concerted or one of two sequential mechanisms. The sequential mechanisms are depicted along the horizontal and vertical edges of the schematic in Fig. 2.1b and involve distinct ET and PT events separated by the  +  3+ metastable charge-separated intermediate FeII (Hbim) + FeIII (H2 bim) . The two mechanisms are distinguished by whether ET proceeds PT or PT proceeds ET. In iron bi-imidazoline, the rates for the two sequential mechanisms are equal based on the symmetry of the system. The concerted mechanism is depicted along the diagonal in Fig. 2.1b and involves the transfer of both the electron and proton in a single reactive step, bypassing the formation of the charge-separated intermediate. Mayer and co-workers concluded that the PCET reaction follows a concerted mechanism by showing that (i) the free-energy cost for forming the charge-separated intermediate, ∆G◦ = 11.5 kcal mol−1 , during either sequential mechanism was almost as large as the observed activation barrier for the PCET reaction, ∆G‡ = 12.3 kcal mol−1 , and (ii) the activation barrier for the uphill ET reaction in the sequential mechanism was calculated to be almost 5 kcal mol−1 higher than the PCET activation barrier. 10 54 The activation barrier for the uphill ET reaction was calculated using the free-energy cost for forming the charge-separated intermediate and assuming that the reorganization energy for the uphill ET reaction was the same as the measured reorganization energy for the single ET reaction, Eq. 2.2. In many cases, the experimental determination of the reorganization energy for the ET reaction in the sequential mechanism is not possible. Consequently, distinguishing between sequential and concerted PCET mechanisms hinges mainly on thermodynamic arguments, leaving the determination of the PCET reaction mechanism ambiguous in many cases. 2,12 Mayer and co-workers state that the large free-energy cost to form the charge-separated intermediates provides the thermodynamic bias for the concerted over sequential mechanism in iron bi-imidazoline. 10 The large free-energy cost is said to arise from the strong interaction between the electron and proton; FeIII H2 bim is ≈0.5 V stronger oxidant than the deprotonated form FeIII Hbim, or similarly FeIII H2 bim is ca. 108 more acidic than FeII H2 bim. This analysis provides only a limited picture of the physical interactions that govern concerted versus sequential PCET, ignoring the role of the solvent in determining the PCET mechanism. Towards this end, we employ the RPMD method to investigate concerted and sequential PCET in iron bi-imidazoline, providing a deeper understanding of the competition between concerted and sequential PCET. RPMD 13 is an approximate real-time quantum dynamical method that is based on Feynman’s imaginary-time path integral formulation of statistical mechanics. 14,15 It provides an isomorphic classical molecular dynamics model for the real-time evolution of a quantum mechanical system. The RPMD method has been previously employed to investigate a variety of quantized reactive and dynamical processes, 16–30 including both ET 22 and PT 16–18 reactions. Recently, we have demonstrated that RPMD simulations can be extended to accurately treat PCET reactions in system-bath models, enabling the direct investigation of concerted and sequential PCET mechanisms. 30 In this paper, we utilize the RPMD method to directly simulate the PCET reaction in both a fully atomistic and system-bath representation of iron bi-imidazoline. We analyze the RPMD trajectories to elucidate the physical interactions that determine the dominant PCET mechanism. We further investigate altering the dominant PCET mechanism through variation of the properties of the chemical system. 55 2.2 Ring polymer molecular dynamics The RPMD equations of motion for a quantized electron, a quantized proton, and N classical particles, including a mixed-bead-number path-integral representation, are 13,20,30,31 v̇e(α)   + qe(α−1) − 2qe(α) = ωn2 e q(α+1) e − v̇p(γ) ! ” “ (α−k) n1ep +1 1 ∇ (α) U qe(α) , qp ,Q , me qe (2.4)   (γ−1) (γ) + q − 2q = ωn2 p q(γ+1) p p p nep −   1 X , ∇q(γ) U qe((γ−1)nep +l) , q(γ) , Q p p mp (2.5) l=1 and V̇j = − np nep   1 XX ∇Qj U qe((γ−1)nep +l) , qp(γ) , Q , ne Mj γ=1 (2.6) l=1 where ne is the number of imaginary-time ring-polymer beads for the transferring electron, me is the physical mass for the electron, and q(α) and ve(α) are the respective position and velocity vectors for e the αth ring-polymer bead of the electron; the corresponding quantities for the transferring proton are indicated using subscript “p”. In Eqs. (2.4)-(2.6), it is assumed that nep = ne /np is an integer number, and  k = α − nep  α−1 , nep (2.7) where b. . . c denotes the floor function. The periodic constraint of the ring-polymer is satisfied via q(0) = qe(ne ) and q(0) = qp(np ) , and the intra-bead harmonic frequencies are ωne = ne /(β~) and e p ωnp = np /(β~), where β = 1/kB T is the inverse temperature. The position, velocity, and mass for the jth classical degree of freedom are given by Qj , Vj , and Mj , respectively, and Q = {Q1 , ..., QN }. Lastly, the potential energy function of the system is given by U (qe , qp , Q). Analogous to the classical thermal rate constant, 32–34 the RPMD thermal rate constant can be expressed as 35,36 kRPMD = lim κ(t)kTST , t→∞ (2.8) where kTST is the transition state theory (TST) estimate for the rate associated with the dividing surface ξ(r) = ξ ‡ , ξ(r) is a collective variable that distinguishes between the reactant and product basins of stability, and κ(t) is the time-dependent transmission coefficient that accounts for recrossing of tra- 56 n o (ne ) (np ) jectories through the dividing surface. We have introduced r = q(1) , q(1) , Q1 , ...QN e , ..., qe p , ...qp to denote the position vector for the full system in the ring-polymer representation. As is the case for both exact classical and exact quantum dynamics, the RPMD method yields reaction rates and mechanisms that are independent of the choice of dividing surface. 35–37 The TST rate in Eq. (2.8) is calculated using 18,22,38,39 ‡ kTST = (2πβ) −1/2 e−β∆F (ξ ) hgξ ic R ξ‡ −∞ dξe−β∆F (ξ) , (2.9) where F (ξ) is the free energy (FE) along ξ, e−β∆F (ξ) = hδ (ξ(r) − ξ)i , hδ (ξ(r) − ξr )i (2.10) ξr is a reference point in the reactant basin, and 18,40–42 gξ (r) = #1/2 " d X 1  ∂ξ(r) 2 i=1 mi ∂ri . (2.11) Here, ri is an element of the position vector r, mi is the corresponding physical mass, and d is the length of vector r. The equilibrium ensemble average is denoted R h. . . i = R dr dv e−βH(r,v) (. . .) R R , dr dv e−βH(r,v) (2.12) and the average over the ensemble constrained to the dividing surface is denoted R h. . . ic = dr R dv e−βH(r,v) (. . . )δ(ξ(r) − ξ ‡ ) R , dr dv e−βH(r,v) δ(ξ(r) − ξ ‡ ) R (2.13) where H(r, v) = N X 1 j=1 np + 2 X1 γ=1 2 Mj Vj2 + ne X 1 α=1 2  2 mb,e ve(α)  2 mb,p vp(γ) + URP (r). (2.14) Here, mb,e and mb,p are the fictitious Parrinello-Rahman masses for the electron and proton, ren o spectively, 43 v = ve(1) , ..., ve(ne ) , vp(1) , ..., vp(np ) , V1 , , VN is the velocity vector for the full system 57 in the ring-polymer representation, and n URP (r) = e  2 1 X 1 me ωn2 e qe(α) − qe(α−1) ne α=1 2 np +  2 1 X1 (γ−1) mp ωn2 p q(γ) − q p p np γ=1 2 (2.15) np nep 1 X X  ((γ−1)nep +l) (γ)  + U qe , qp , Q . ne γ=1 l=1 The transmission coefficient in Eq. (2.9) is obtained from the flux-side correlation function, 35,36 κ(t) = hξ˙0 h(ξ(rt ) − ξ ‡ )ic , hξ˙0 h(ξ˙0 )ic (2.16) by releasing RPMD trajectories from the equilibrium ensemble constrained to the dividing surface. Here, h(ξ) is the Heaviside function, ξ˙0 is the time-derivative of the collective variable upon initialization of the RPMD trajectory from the dividing surface with the initial velocities sampled from the Maxwell-Boltzmann (MB) distribution, and rt is the time-evolved position of the system along the RPMD trajectory. 2.3 Systems The competition between concerted and sequential PCET reaction mechanisms is studied using both all-atom and system-bath representations of the iron bi-imidazoline system. These representations are described below. 2.3.1 Atomistic Representation for PCET The atomistic representation of the PCET reaction in iron bi-imidazoline consists of a single quantized electron, an FeIII (H2 bim) and an FeIII (Hbim) complex solvated in explicit acetonitrile. The single transferring proton on the FeIII (H2 bim) complex is also quantized using RPMD. The rest of the nuclei on both iron complexes and the solvent are treated classically. The atomic coordinates of both iron complexes are obtained from the experimentally obtained crystal structure of FeIII (H2 bim) to maintain symmetry of the PCET reaction 10 . 58 The potential energy function that describes the atomistic representation is given by U (qe , qp , Q) = Ucl (Q) + Up (qp , Q) + Ue (qe , Q) + Uep (qe , qp ), (2.17) where Q is the set of atomic positions for all of the classical nuclei in both the solvent and iron bi-imidazoline complexes. The interactions between all of the classical nuclei, Ucl (Q), are described by a modified version of the Generalized Amber Force Field (GAFF), in which the united-atom approximation was used for hydrogens bonded to carbons; 44 hydrogens bonded to nitrogens are treated explicitly. The modifications to the GAFF are that (i) the parameters for acetonitrile are obtained from the threesite model of Guardia et. al. 45 , (ii) the charges on the iron centers are chosen to be qFe = 1.65 e, such that the rate for the concerted PCET reaction is in agreement with experiment, and (iii) the charges on the bi-imidazoline ligands are obtained through the procedure described below. The steps taken to obtain the atomic charges on the bi-imidazoline ligands are (i) calculate the atomic charges using the CHELPG method on the isolated protonated, H2 bim, and deprotonated, Hbim, ligands in a continuum solvent representation of acetonitrile. 46 The geometries of the ligands were optimized at the RHF/6-31G** level of theory invoking C2v and Cs symmetry, respectively. All electronic structure calculations were performed using the Gaussian09 package. 47 (ii) The CHELPG charges on the hydrogens bonded to carbons are added to the carbon charges in accordance with the united-atom approximation. (iii) The CHELPG charges of the 12 nitrogen atoms directly bonded to the two iron centers are evenly shifted such that the total charges of the FeIII (H2 bim) and FeIII (Hbim) complexes were +3 and +2, respectively, CHELPG qNj = qN + j 3e − qFe , 6 (2.18) CHELPG where qNj is the shifted charge on nitrogen atom j, qN is the charge on nitrogen atom j j obtained from the CHELPG calculation and qFe is the charge on an iron atom obtained by fitting the rate for the concerted PCET reaction to the experimental PCET rate as described above; the charges on both iron atoms are equal. The term 3e − qFe accounts for the difference in the charge on an iron atom used in the simulations and the formal +3 redox state of the iron atoms; this term is divided by six to account for the six nitrogen atoms bonded to each iron atom. Steps (i)-(iii) fully specify the atomic charges of the four bi-imidazoline ligands that do not participate in the hydrogen 59 bond with the transferring proton. The instantaneous atomic charges on the two bi-imidazoline ligands that do participate in the hydrogen bond (specifically the H(2 bim) ligand that includes the transferring proton in the FeIII (H2 bim) complex and the Hbim ligand in the FeIII (Hbim) complex) change between the values corresponding to the H2 bim and Hbim ligands obtained in steps (i)-(iii) depending on the position of the proton as follows 48 (p) (dp) where qj , qj , and qj qj (r) = f (r) = (p) (dp) (1 − f (r))qj + f (r)qj , " # 1 r − r0 1+ p , 2 (r − r0 )2 + l2 (2.19) (2.20) are the instantaneous atomic charge, the atomic charge in the protonated H2 bim ligand, and the atomic charge in the deprotonated Hbim ligand associated with nuclei j. The parameter r0 is given by half the instantaneous distance between the nitrogen atoms that are directly participating in the hydrogen bond with the transferring proton, and l = 0.125 Å is chosen in agreement with previous work. 18 The variable r is given by the distance between the transferring proton and the nitrogen atom participating in the hydrogen bond associated with the same iron complex corresponding to nuclei j. To conserve charge, the charge of the transferring proton, qp , is given by qp = 4e − X qj , (2.21) j where the sum over j runs over all atoms in the system, and the value of 4e corresponds to the total charge of the system. The interaction between the transferring proton and the classical nuclei in the solvent and both iron complexes is given by Up (qp , Q) = Up,coul (qp , Q) + Up,lj (qp , Q) + UHB (qp , QND , QNA ). (2.22) The potentials Up,coul (qp , Q) and Up,lj (qp , Q) correspond to the usual Coulombic and Lennard-Jones interaction between the transferring proton and all of the classical nuclei except for the two nitrogen atoms participating in the hydrogen bond. The potential describing the hydrogen bond between the 60 proton and the two nitrogen atoms is given by an extension of the Azzouz-Borgis model for PT 48 UHB (qp , QND , QNA ) = Urep (QND , QNA ) +    −n(rHD − d)2 D 1 − exp + 2rHD    −n(rHA − d)2 D 1 − exp , 2rHA (2.23) where     σ 12  σ 6   4ε − +ε R R Urep (QND , QNA ) =   0 R < 21/6 σ , (2.24) R ≥ 21/6 σ and where QND and QNA are the positions of the two nitrogen atoms, one associated with the donor and one with the acceptor complex, respectively. The variables rHD and rHA are the distances between the transferring proton and the nitrogen atoms associated with the donor and acceptor complexes, respectively. The parameters D = 93 kcal/mol, n = 11 Å−1 and d = 1 Å are chosen from common experimental values for nitrogen-proton bonds. 11,49 The potential Urep (QND , QNA ) accounts for the core repulsion between the two nitrogen atoms, where R = |QND − QNA | is the distance between the two nitrogen atoms and  = 250 kcal/mol and σ = 2.39 Å are chosen such that the average distance between the two nitrogen atoms corresponds to the distance in the experimental crystal structure of 2.67 Å. 10 The interaction between the transferring electron and the classical nuclei in the solvent and both iron complexes is given by Ue (qe , Q) = Ue,coul (qe , Q) + Ue,rep (qe , Q). (2.25) The potential Ue,coul (qe , Q) describes the scaled pairwise pseudopotential between the electron and the classical nuclei in the solvent and both iron complexes 50 Ue,coul (qe , Q) = X i∈solv i Ue−solv (ri ) + X j Ue−cmplx (rj ), (2.26) j∈cmplx where ri = |qe − Qi | and rj = |qe − Qj |. The atom index i corresponds to nuclei only in the solvent, and j corresponds to nuclei only in the iron complexes. For cases in which the atom index 61 i corresponds to a positively charged atom, i Ue−solv (ri ) = qi e , i 4π0 rcut   −ζsolv qi e , 4π0 ri    −ζsolv i r ≤ rcut , (2.27) i r > rcut and when i corresponds to a negatively charged atom, i Ue−solv (ri ) = −ζsolv qi e . 4π0 ri (2.28) Similarly, when j corresponds to a positively charged atom, j Ue−cmplx (rj )     −ζcmplx qj e j 4π0 rcut = qj e   ,  −ζcmplx 4π0 rj , j r ≤ rcut , (2.29) j r > rcut and when j corresponds to a negatively charged atom, j Ue−cmplx (rj ) = −ζcmplx qj e . 4π0 rj (2.30) The value of the scaling parameter ζsolv = 0.74 is chosen such that the atomistic representation reproduces the previously calculated value for the outer-sphere reorganization energy of the symmetric ET reaction (Eq. 2.2), λo = 17.9 kcal/mol, as described in 2.4.1.7. 10 The value of the scaling parameter ζcmplx = 0.56 is chosen such that the atomistic representation reproduces the experimental driving force at the reactive configuration for the formation of the PCET charge separated 0 intermediate, ∆G◦ = 11 kcal/mol, as described in Sec. 2.4.1.3. 10 The potential Ue,rep (qe , Q) describes an additional pseudopotential used to model the repulsion between the transferring electron and the electron cloud associated with each nuclei 20,51 Ue,rep (qe , Q) =  X A B − 1 , rk4 (C + rk6 ) (2.31) k6∈Fe,H where rk = |qe − Qk | and the atomic index k runs over all nuclei except proton and iron atoms. The values of the parameters in Eq. 2.31 are A = 32.2 kcal/mol Å4 , B = 1956.5 Å6 and C = 276.86 Å6 . The parameters above, which we denote System 1a, fully define an atomistic representation of the PCET reaction in iron bi-imidazoline that represents the original experimental conditions. 10 In addition, we further define sets of parameters, which allows for the investigation of the physical 62 interactions that govern the dominant PCET reaction. First, we define an analogous set of parameters to those above, Systems 1b and 1c, but in which the atomic charges on the acetonitrile solvent molecules are varied by a multiplicative factor to model different solvent conditions of varying polarity. Second, we define a set of parameters, Systems 2a-2c, in which (i) the atomic charges on the acetonitrile molecules are varied by a multiplicative factor and (ii) the value of the scaling parameter ζcmplx = 0.30 is set to model molecular systems in which the implicit strength of interaction between the electron and proton is weaker than in iron bi-imidazoline. 2.3.2 System-bath Representation for PCET In addition to the atomistic representation of the PCET reaction in iron bi-imidazoline presented above, we also employ a co-linear system-bath model for PCET. The system-bath model has been described in detail previously and is thus only summarized below 30 . The model is expressed in the position representation using the potential energy function U (qe , qp , qs , Q) = Usys (qe , qp , qs ) + UB (qs , Q), (2.32) where UB (qs , Q) is the potential energy term associated with the bath coordinates, and Usys (qe , qp , qs ) = Ue (qe ) + Up (qp ) + Us (qs ) +Ues (qe , qs ) + Ups (qp , qs ) +Uep (qe , qp ) (2.33) is the system potential energy. The scalar coordinates qe , qp , and qs describe the one-dimensional (1D) positions of the electron, proton, and solvent modes, respectively, and Q is the vector of bath oscillator positions. The first term in the system potential energy function models the interaction of the transferring electron with its donor and acceptor sites,  out in   aD qe2 + bD qe + cD , rD ≤ qe ≤ rD           in out aA qe2 + bA qe + cA , rA ≤ qe ≤ rA Ue (qe ) =            1 1    −µe + , otherwise, |qe − rD | |qe − rA | (2.34) 63 where rD and rA are the positions of the electron donor and acceptor sites. The second term in the system potential energy function models the interaction between the transferring proton and its donor and acceptor sites, Up (qp ) = − mp ωp2 2 m2p ωp4 4 qp + q . 2 16V0 p (2.35) Here, ωp is the proton vibrational frequency and V0 is the intrinsic PT barrier height. The next three terms in the system potential energy function model the solvent potential and the electron- and proton-solvent interactions. Specifically, 1 ms ωs2 qs2 , 2 (2.36) Ues (qe , qs ) = −µes qe qs , (2.37) Ups (qp , qs ) = −µps qp qs , (2.38) Us (qs ) = and where ms is the solvent mass and ωs is the effective frequency of the solvent coordinate. Interactions between the transferring electron and proton are modeled via the capped coulombic potential Ue (qe ) =  µep   − ,     |qe − qp | |qe − qp | > Rcut    µ    − ep , Rcut otherwise. (2.39) The potential energy term UB (qs , Q) models the harmonic bath that is coupled to the PCET reaction. The bath exhibits an ohmic spectral density J(ω) with cutoff frequency ωc , 52,53 such that J(ω) = ηωe−ω/ωc , (2.40) where η denotes the friction coefficient. The continuous spectral density is discretized into f oscillators with frequencies 35,54   j − 0.5 ωj = −ωc ln f (2.41) and coupling constants  cj = ωj 2ηM ωc fπ 1/2 , (2.42) 64 such that  f X  1 M ωj2 UB (qs , Q) = 2 j=1 cj qs Qj − M ωj2 !2  . (2.43) Here, M is the mass of each bath oscillator, and ωj and Qj are the respective frequency and position for the jth oscillator. We have developed system parameters to model condensed-phase PCET reactions that transition between the concerted and sequential mechanisms. Specifically, Systems 3a-3e vary the strength of the solvent-proton and solvent-electron interactions, Systems 4a-4e vary the strength of the electronproton interaction, and Systems 5a-5e vary the barrier height associated with PT. The parameters associated with the system-bath model of PCET are presented in Appendices 2.1 and 2.2 . 2.4 Calculation details 2.4.1 Atomistic representation The atomistic simulations for all Systems (Systems 1a-1c and Systems 2a-2c) are implemented in the DL POLY molecular dynamics package, and include 410 acetonitrile molecules. 55 In all simulations, the RPMD equations of motion are evolved using the velocity Verlet algorithm. 56 The electron is quantized with ne = 1024 ring-polymer beads, and the proton is quantized with np = 32 ringpolymer beads. As in previous RPMD simulations, each time step for the electron and proton involves separate coordinate updates due to forces arising from the physical potential and due to exact evolution of the purely harmonic portion of the ring-polymer potentials. 57 The temperature is set to the experimental value of 298 K. 10 All pair-wise interactions are truncated at a distance of rpw = 12 Å. Long-range electrostatics, including the Coulombic interactions between classical nuclei, the Coulombic interaction between the proton and the classical nuclei (Up,coul ), and the Coulombic interactions between the electron and the classical nuclei (Ue,coul ), are treated with the force-shifting algorithm, 58 in which the Coulombic potential is multiplied by a damping function S(r), such that both the potential and its derivative smoothly vanish at r = rpw , where here r defines the distance between the two particles participating in the pair-wise Coulombic interaction. Specifically,  2r r2   1− + 2 , r ≤ rpw rpw rpw S(r) =   0, r > rpw . (2.44) 65 Table 2.1: Side-lengths of the simulation cell for the atomistic representation of PCET. a System L1 L2 L3 1a 1b 1c 2a 2b 2c 29.498 28.778 28.179 30.236 29.498 29.103 29.498 28.778 28.170 30.236 29.498 29.103 45.221 44.116 43.198 46.352 45.221 44.116 The parameters are given in Å, where L1 , L2 , and L3 are the three sides of the rectangular simulation cell. All atomistic calculations are performed in a rectangular simulation cell with periodic boundary conditions. The side-lengths of the cell for each System are obtained from 1.5 ns NPT equilibrium simulations run with the Nosé-Hoover barostat and thermostat using a thermostat and barostat relaxation time of 1.0 ps and 2.0 ps, respectively. The side-lengths of the cell for each System are presented in Table 2.1. Several collective variables are used to monitor and characterize the PCET reaction in the atomistic representation. The progress of the electron is characterized by a “bead-count” coordinate, fb , that reports on the fraction of ring-polymer beads that are located on the iron atom associated with the donor complex,   ne   1 X φ (α) (ne ) fb q(1) , . . . , q , Q , Q = tanh f , FeD FeA e e ne α=1 fFeD (2.45) where fFeD (QFeD , QFeA ) = 1 − |QFeA − QFeD | 2 (2.46) and f (α)   (ne ) q(1) , QFeD , QFeA = e , . . . , qe   1 α q − (QFeA + QFeD ) · 2   1/2 (QFeA − QFeD ) |1/2 (QFeA − QFeD )| (2.47) The variables QFeD , and QFeA are the positions of the iron atoms associated with the donor and acceptor complex, respectively, and φ = −3.0. The progress of the proton is characterized by the difference between the distances of the ring- 66 polymer centroid and the position of the two nitrogen atoms participating in the hydrogen bond,   (np ) fp q(1) , QND , QNA = p , . . . , qp |QND − q̄p | − |QNA − q̄p | , (2.48) where np q̄p  (np ) q(1) p , . . . , qp  1 X (γ) q . np γ=1 p = (2.49) The progress of the solvent during the concerted PCET reaction is characterized by the energy gap associated with the transfer of both the electron and proton, ∆U (Q) =    −e X qi qi − ζsolv 4π0 |Qi − QFeA | |Qi − QFeD | i∈solv ! X qi qjD qi qjA , − + |Qi − Qj | |Qi − Qj | (2.50) j∈cmplx where qjD and qjA are the atomic charges associated with nuclei j when the proton is bonded to the (p) donor or acceptor complex, respectively. Thus, qjD = qj (dp) donor complex; qjD = qj (p) and qjA = qj (dp) and qjA = qj if j is associated with the if j is associated with the acceptor complex. Additionally, the progress of the solvent during the ET step in the sequential mechanism and during the single ET reaction (Eq. (2.2)) is characterized by the energy gap associated with just the transfer of the electron, ∆U ET (Q) =   −e X qi qi ζsolv − . 4π0 |Qi − QFeA | |Qi − QFeD | (2.51) i∈solv 2.4.1.1 RPMD rate calculations for concerted PCET The RPMD reaction rate is calculated from the product of the TST rate and the transmission coefficient (Eq. (2.8)). The FE profiles that appear in the TST rate expression (Eq. (2.9)) are obtained using umbrella sampling as described below. For all atomistic systems, the 1D FE profile used in the rate calculation is obtained in the electron bead-count coordinate, F (fb ), using the following umbrella sampling protocol. Forty-five independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.0, 0.88] using a force constant of 10,000 kcal/mol; six independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.90, 0.95] using a 67 higher force constant of 20,000 kcal/mol; two independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.96, 0.97] using a higher force constant of 40,000 kcal/mol; ten independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.98, 0.998] using a higher force constant of 1 × 106 kcal/mol; an additional independent sampling trajectory is harmonically restrained to a value of fb = 0.989 using a force constant of 1 × 106 kcal/mol to ensure extensive overlap among the sampled distributions. The symmetry of the reaction is employed to obtain the full FE profile along fb over the region [-0.998,0.998]. For all Systems, an auxiliary restraining potential is introduced to the sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.C.1. The equilibrium sampling trajectories are performed using path-integral molecular dynamics (PIMD) with mb,e = 5.357 g/mol and mb,p = 3.156 g/mol, which allows for a timestep of 1 fs. Each sampling trajectory is run for at least 1 ns, and thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps. It is important to note, that as always, the choice of the Parinello-Rahman masses, mb,e and mb,p , allows for a large time step in the sampling trajectories, but has no affect of F (fb ) or any other equilibrium ensemble average. 43,57,59 For all atomistic Systems, the transmission coefficient is calculated using RPMD trajectories that are released from the dividing surface associated with fb = 0. At least 3000 trajectories are released for each system. Each RPMD trajectory is evolved for 300 fs using a timestep of 5 × 10−4 fs and with initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 1 ps from long PIMD sampling trajectories that are constrained to the dividing surface using the RATTLE algorithm. 60 The sampling trajectories utilize mb,e = 5.357 g/mol, mb,p = 3.156 g/mol, and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps. For all Systems, the same auxiliary potential used in the calculation of F (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.C.1; throughout this paper, the PRMD trajectories used to calculate the transmission coefficients are not subject to any auxiliary restraining potentials. 68 2.4.1.2 RPMD rate calculations for ET prior to PT For the ET step in the sequential mechanism, we calculate the forward ET reaction reaction rate for the reaction  II 2+  III 2+ Fe (H2 bim)3 + Fe (Hbim)(H2 bim)2 3+  II 1+ ET  III − − * + Fe (Hbim)(H2 bim)2 . ) − − Fe (H2 bim)3 (2.52) The 1D FE profile used in the rate calculation for the ET reactions is obtained in the electron bead-count coordinate, FSET (fb ), using the following umbrella sampling protocol. Eighty-nine independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [−0.88, 0.88] using a force constant of 10000 kcal/mol; six independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the both the region [−0.95, −0.90] and [0.90, 0.95] using a higher force constant of 20000 kcal/mol; two independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the both the region [−0.97, −0.96] and [0.96, 0.97] using a higher force constant of 40000 kcal/mol; ten independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in both the region [−0.998, 0.98] and [0.98, 0.998] using a higher force constant of 1 × 106 kcal/mol; two additional independent sampling trajectories are harmonically restrained to the values of fb = 0.989 and fb = −0.989 using a force constant of 1 × 106 kcal/mol. For all Systems, an auxiliary restraining potential is introduced to the sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.C.2. As before, the PIMD sampling trajectories employ mb,e = 5.357 g/mol, mb,p = 3.156 g/mol, and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for at least 1 ns. The transmission coefficient is calculated using RPMD trajectories that are released from the dividing surface associated with fb = 0.62325. At least 3000 trajectories are released for each system. Each RPMD trajectory is evolved for 300 fs using a timestep of 5 × 10−4 fs and with initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 1 ps from long PIMD sampling trajectories that are constrained to the dividing surface using the RATTLE algorithm. 60 The sampling trajectories utilize mb,e = 5.357 g/mol, mb,p = 3.156 g/mol, and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps. For all Systems, the same auxiliary potential used in the calculation of FSET (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the concerted 69 channel, as described in Appendix 2.C.2 2.4.1.3 1D FE profile for PT prior to ET For the purpose of analysis and for the calculation of the parameter ζcmplx as described in Sec. 2.3.1, we additionally calculate the 1D FE profile, FSPT (fp ), for the PT step in the sequential mechanism corresponding to the reaction  III 3+  II 1+ Fe (H2 bim)3 + Fe (Hbim)(H2 bim)2 2+  II 2+ PT  III − − * + Fe (H2 bim)3 , ) − − Fe (Hbim)(H2 bim)2 (2.53) using the following umbrella sampling protocol. For each atomistic System, fifteen independent sampling trajectories are harmonically restrained to uniformly spaced values of fp in the region [-0.35 Å,0.35 Å] using a force constant of 1000 kcal/mol Å−2 and seven independent sampling trajectories are harmonically restrained to uniformly spaced values of fp in both the region [-0.70 Å,-0.40 Å] and [0.40 Å,0.70 Å] using a force constant of 500 kcal/mol Å−2 . The PIMD sampling trajectories employ mb,e = 5.357 g/mol, mb,p = 3.156 g/mol, and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for 500 ps. In addition, the electron ring-polymer is initialized to the position of the iron atom associated with the donor iron bi-imidazoline complex for each sampling trajectory. 2.4.1.4 Two-dimensional FE profiles We additionally calculate the two-dimensional (2D) FE profile for System 1a in the electron beadcount and proton coordinates, F (fb , fp ). The 2D FE profile is constructed using PIMD sampling trajectories that are harmonically restrained in both the fb and fp coordinates. A total of 1856 sampling trajectories are performed, in which the coordinates fb and fp are sampled using a square grid. The coordinate fb is sampled using forty-five independent sampling trajectories that are harmonically restrained to uniformly spaced values of fb in the region [0.0, 0.88] using a force constant of 10,000 kcal/mol; six independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.90, 0.95] using a higher force constant of 20,000 kcal/mol; two independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.96, 0.97] using a higher force constant of 40,000 kcal/mol; ten independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [0.98, 0.998] using a higher force constant of 1 × 106 kcal/mol; an additional independent sampling trajectory 70 is harmonically restrained to a value of fb = 0.989 using a force constant of 1 × 106 kcal/mol. The symmetry of the reaction is employed to obtain the full FE profile along fb over the region [-0.998,0.998]. The coordinate fp is sampled using fifteen independent sampling trajectories that are harmonically restrained to uniformly spaced values of fp in the region [-0.35 Å,0.35 Å] using a force constant of 1000 kcal/mol Å−2 and seven independent sampling trajectories are harmonically restrained to uniformly spaced values of fp in both the region [-0.70 Å,-0.40 Å] and [0.40 Å,0.70 Å] using a force constant of 500 kcal/mol Å−2 . No auxiliary restraining potentials are employed for the calculation of F (fb , fp ). The PIMD sampling trajectories employ mb,e = 5.357 g/mol, mb,p = 3.156 g/mol, and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for 500 ps. We additionally calculate the 2D FE profile for System 1a in the electron bead-count and concerted PCET energy gap coordinates, F (fb , ∆U ), for sampling trajectories corresponding to the concerted PCET reaction. To generate F (fb , ∆U ), the harmonically restrained sampling trajectories used to calculate F (fb ) for System 1a are utilized. 2.4.1.5 Solvent reorganization energy for concerted PCET For Systems 1a-1c , we calculate the outer-sphere reorganization energy associated with concerted PCET from the equation λCPET = FDCPET (∆UA ) − FDCPET (∆UD ), o (2.54) where ∆UD and ∆UA correspond to the minimum value of ∆U in the FE profiles FDCPET (∆U ) and FACPET (∆U ), corresponding to the electron and proton being associated with the donor or acceptor, respectively. The 1D FE profiles FDCPET (∆U ) and FACPET (∆U ) are calculated along the energy gap coordinate, ∆U , using the following umbrella sampling protocol. The coordinate ∆U is sampled using seventeen independent sampling trajectories that are harmonically restrained to uniformly spaced values of ∆U in the region [-40 kcal/mol,40 kcal/mol] using a force constant of 1.20 × 10−6 kcal/mol (kcal/mol)−2 . The sampling trajectories used to calculate the FE profile FDCPET (∆U ) and FACPET (∆U ) are initialized with the electron and proton associated with the donor or acceptor iron bi-imidazoline complex, respectively. The PIMD sampling trajectories employ mb,e = 5.357 g/mol, mb,p = 3.156 g/mol and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for 1 ns. 71 2.4.1.6 Solvent reorganization energy for ET prior to PT For Systems 1a-1c , we calculate the outer-sphere reorganization energy associated with the ET reaction prior to PT in the sequential mechanism from the equation λETPT = FDETPT (∆UAET ) − FDETPT (∆UDET ), o (2.55) where ∆UDET and ∆UAET correspond to the minimum value of ∆U ET in the FE profiles FDETPT (∆U ET ) and FAETPT (∆U ET ), corresponding to the electron being associated with the donor or acceptor, respectively. The 1D FE profiles FDETPT (∆U ET ) and FAETPT (∆U ET ) are calculated along the energy gap coordinate, ∆U ET , using the following umbrella sampling protocol. The coordinate ∆U ET is sampled using seventeen independent sampling trajectories that are harmonically restrained to uniformly spaced values of ∆U ET in the region [-40 kcal/mol,40 kcal/mol] using a force constant of 1.20 × 10−6 kcal/mol (kcal/mol)−2 . The sampling trajectories used to calculate the FE profile FDETPT (∆U ET ) and FAETPT (∆U ET ) are initialized with the electron ring-polymer at the position of the iron atom associated with the donor or acceptor iron bi-imidazoline complex, respectively; for the calculation of both FE profiles the proton is associated with the donor iron bi-imidazoline complex and an auxiliary restraining potential is introduced to the sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.C.2. The PIMD sampling trajectories employ mb,e = 5.357 g/mol, mb,p = 3.156 g/mol and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for 1 ns. 2.4.1.7 Solvent reorganization energy for symmetric ET For the calculation of the parameter ζsolv as described in Sec. 2.3.1, we additionally calculate the outer-sphere reorganization energy associated with the single ET reaction (Eq. (2.2)) from the equation λETPT = FDETPT (∆UAET ) − FDETPT (∆UDET ), o (2.56) where ∆UDET and ∆UAET correspond to the minimum value of ∆U ET in the FE profiles FDET (∆U ET ) and FAET (∆U ET ), corresponding to the electron being associated with the donor or acceptor, respectively. The 1D FE profiles FDETPT (∆U ET ) and FAETPT (∆U ET ) are calculated along the energy gap coordinate, ∆U ET , using the following umbrella sampling protocol. The coordinate ∆U ET is sampled using seventeen independent sampling trajectories that are harmonically restrained to uniformly spaced values of ∆U ET in the region [-40 kcal/mol,40 kcal/mol] using a force constant 72 of 1.20 × 10−6 kcal/mol (kcal/mol)−2 . The sampling trajectories used to calculate the FE profile FDET (∆U ET ) and FAET (∆U ET ) involve two fully protonated FeIII (H2 bim) complexes, in which all protons are treated classically and the electron ring-polymer is initialized to the position of the iron atom associated with the donor or acceptor iron bi-imidazoline complex, respectively. The PIMD sampling trajectories employ mb,e = 5.357 g/mol and a time-step of 1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 3 ps, and each trajectory is run for 1 ns. 2.4.1.8 RPMD transition path ensemble As we have done previously, 17,30 we analyze the transition path ensemble 61 for the RPMD trajectories in the current study. Reactive trajectories are generated through forward-and backwardintegration of initial configurations drawn from the dividing surface ensemble with initial velocities drawn from the MB distribution. Reactive trajectories correspond to those for which forward- and backward-integrated half trajectories terminated in opposite sides of the dividing surface. The reactive trajectories that are initialized from the equilibrium Boltzmann distribution on the dividing surface must be reweighted to obtain the unbiased transition path ensemble. 61–63 A weighting term, wα , is applied to each trajectory, correctly accounting for recrossing and for the fact that individual trajectories are performed in the microcanonical ensemble. This term is given by 62 wα = X ˙ i ξ(r) −1 !−1 , (2.57) i ˙ i is where the sum includes all instances in which trajectory α crosses the dividing surface, and ξ(r) the velocity in the dividing surface collective variable at the ith crossing event. The reweighting has a minor effect on the non-equilibrium averages if the reactive trajectories initialized from the dividing surface exhibit relatively little recrossing, as is the case for the systems studied in this paper. Nonequilibrium averages over the RPMD transition path ensemble are calculated by aligning reactive trajectories at time 0, defined as the moment in time when the trajectories are released from the dividing surface. 2.4.2 System-bath representation The simulations for the system-bath models of PCET (Systems 3a-3e, Systems 4a-4e, and Systems 5a-5e) are all performed at T = 300 K. The RPMD equations of motion are evolved using the velocity 73 Verlet algorithm, 56 and each time step for the electro band proton involves separate coordinate updates due to forces arising from the physical potential and due to exact evolution of the purely harmonic portion of the ring-polymer potentials. 57 The electron is quantized with ne = 1024 ringpolymer beads, and the proton is quantized with np = 32 ring-polymer beads. The progress of the electron in the system-bath models is monitored using a 1D form of the electron bead-count coordinate, ne     1 X fb qe(1) , ..., qe(ne ) = tanh θqe(α) , ne α=1 (2.58) where θ = −3.0/rD . The progress of the proton in the system-bath models is monitored using the ring-polymer centroid in the proton position coordinate, np q̄p 2.4.2.1  qp(1) , ..., qp(np )  1 X (γ) q . = np γ=1 p (2.59) RPMD rate calculations for concerted PCET As in the atomistic representation for PCET, the RPMD reaction rate in the system-bath models is calculated from the product of the TST rate and the transmission coefficient (Eq. (2.8)). The 1D FE profile used in the rate calculation for the concerted PCET reactions is obtained in the electron beadcount coordinate, F (fb ), using the following umbrella sampling protocol. 93 independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in the region [−0.92, 0.92] using a force constant of 20 a.u.; seven independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in both the region [−0.991, −0.985] and in [0.985, 0.991] using a higher force constant of 5000 a.u.; nine independent sampling trajectories are harmonically restrained to uniformly spaced values of fb in both the region [−1.0, −0.992] and in [0.992, 1.0] using a higher force constant of 10, 000 a.u.; 32 independent sampling trajectories are harmonically restrained to the values of fb ∈ { ±0.93, ±0.935, ±0.94, ±0.945, ±0.95, ±0.955, ±0.96, ±0.962, ±0.965, ±0.967, ±0.97, ±0.974, ±0.976, ±0.978, ±0.98, ±0.982} using a force constant of 500 a.u. For all Systems, an auxiliary restraining potential is introduced for the PIMD sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.C.3. Each sampling trajectory is run for 10 ns using a timestep of 0.1 fs, with mb,e = 2000 a.u. and mb,p = 1836.1 a.u. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The transmission coefficient is calculated using a total of 6000 RPMD trajectories that are 74 Table 2.2: Dividing surfaces for the sequential ET reaction prior to PT in the system-bath models. System fb‡ 1a 1b 1c 1d 1e 2a 2b 2c 2d 2e 3a 3b 3c 3d 3e 8.00 × 10−1 3.76 × 10−1 1.78 × 10−1 −2.03 × 10−2 5.63 × 10−2 8.32 × 10−2 2.25 × 10−3 1.30 × 10−1 1.80 × 10−1 3.00 × 10−1 7.88 × 10−2 5.18 × 10−2 5.18 × 10−2 7.88 × 10−2 6.52 × 10−2 released from the dividing surface associated with fb = 0. Each RPMD trajectory is evolved for 300 fs using a timestep of 1 × 10−4 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u., and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. The same auxiliary restraining potential used in the calculation of F (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the concerted channel, as described in Appendix 2.C.3. 2.4.2.2 RPMD rate calculations for ET prior to PT We calculate the rate for the forward ET reaction in the sequential PCET mechanism. The 1D FE profile used in the rate calculation for the ET reactions is obtained in the electron bead-count coordinate, FSET (fb ), using the same umbrella sampling protocol described for the calculation of F (fb ); however, in the calculation of FSET (fb ), an auxiliary restraining potential is introduced for the PIMD sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.C.4. The independent sampling trajectories used to calculate FSET (fb ) are each run for 10 ns. The transmission coefficients for the forward ET reactions are calculated using RPMD trajectories that are released from the dividing surfaces present in Table 2.2. A total of 6000 RPMD trajectories are released for each System. ach RPMD trajectory is evolved for 300 fs using a timestep of 1×10−4 fs 75 and with the initial velocities sampled from the MB distribution. Initial configurations for the RPMD trajectories are selected every 10 ps from long PIMD sampling trajectories that are constrained to the dividing surface. The sampling trajectories employ mb,e = 2000 a.u., mb,p = 1836.1 a.u., and a timestep of 0.1 fs. Thermostatting is performed by re-sampling the velocities from the MB distribution every 500 fs. The sampling trajectories are constrained to the dividing surface using the RATTLE algorithm. The same auxiliary restraining potential used in the calculation of FSET (fb ) is introduced for the PIMD sampling trajectories to restrict the system to the ET channel, as described in Appendix 2.C.4. 2.5 Results and discussion The results are presented in two sections. In the first, we analyze the competition between the concerted and sequential PCET mechanisms in iron bi-imidazoline using the fully atomistic representation of iron bi-imidazoline. In the second, we vary the strength of different physical interactions to investigate how to alter the dominant PCET mechanism using both the atomistic and system-bath representation of PCET. 2.5.1 Understanding concerted versus sequential PCET in iron bi-imidazoline We begin by examining the competing PCET reaction mechanisms in iron bi-imidazoline under experimental conditions, System 1a. 10 Figure 2.2 presents the 2D FE profile for this system along the electron bead-count, fb , and proton collective variable, fp . The FE profile exhibits four distinct minima corresponding to the reactant, product, and charge-separated intermediate states depicted in Fig. 2.1b. Snapshots of the RPMD simulations of the atomistic representation in each minima are also illustrated in Fig. 2.2. Distinct channels on the FE surface connect the various basins of stability. Overlaid on the FE surface are representative samples from the ensemble of reactive RPMD trajectories for PCET. The trajectories cluster within the channels, illustrating the competing sequential (purple and orange) and concerted (red) reaction mechanisms. The concerted trajectories span directly from the basins corresponding to FeII (H2 bim)+FeIII (Hbim) and FeIII (Hbim)+FeII (H2 bim), bypassing the formation of the intermediate FeII (Hbim)+FeIII (H2 bim) species; the sequential trajectories, in comparison, span along the horizontal and vertical of the FE surface, including the basin corresponding to the intermediate species. Similar results have been seen for the investigation of PCET reactions in system-bath models. 30 76 FeII(Hbim) + ! FeIII(H2bim)! FeIII(Hbim) + ! FeII(H2bim)! FeII(H2bim) + ! FeIII(Hbim)! FeIII(H2bim) + ! FeII(Hbim)! Figure 2.2: Reactive RPMD trajectories reveal distinct concerted (red), and sequential ET-PT (orange) and PT-ET (purple) reaction mechanisms for the PCET reaction in the atomistic representation of iron bi-imidazoline. The trajectories are projected onto the FE surface in the electron bead-count coordinate, fb , and the proton collective variable, fp , with contour lines indicating FE increments of 1 kcal/mol. Snapshots of the RPMD simulations from each of the four PCET basins of stability are also shown. We now demonstrate that the concerted PCET mechanism is dominant for the atomistic representation of iron bi-imidazoline in agreement with experiment by calculating the rates for the concerted and sequential mechanisms. The rate constant for the bimolecular (second-order) concerted PCET reaction may be expressed as 11,64,65 CPET CPET kbi = KA (r)kuni , (2.60) where KA is the equilibrium constant for the formation of the precursor complex at a separation CPET distance between the two iron atoms, r, and kuni is the uni-molecular (first-order) rate constant for the concerted PCET reaction. The equilibrium constant is expressed as KA (r) = Pr exp (−βwr ) , (2.61) where wr is the work to bring the two reacting iron complexes together. The prefactor Pr can be approximated as Pr = 4πNA r2 δr × 10−27 , (2.62) where δr is the range of iron-iron distances over which the rate is appreciable. Here, Pr is given in units of inverse moles per liter, and r and δr are given in angstroms. In this paper, r = 10.3 Å is 77 given by the iron-iron distance in the crystal structure of FeIII Hbim, δr = 0.8 Å, which has shown to provide reasonable results, 11,65 and wr = 1.35 kcal mol−1 has been previously calculated. 10,11 F(fb) / kcal mol -1 (a) 10 10 5 0 0 -1 0.8 -0.5 0.9 0 1 0.5 1 fb F(fb) / kcal mol-1 (b) 20 20 10 15 10 0 -1 -0.5 0.8 0.9 0 1 0.5 1 fb F(fp) / kcal mol-1 (c) 15 10 5 0 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 fp / Å Figure 2.3: (a) The 1D FE profile in the electron bead-count coordinate, F (fb ), associated with the concerted PCET reaction. (b) The 1D FE profile in the electron bead-count coordinate, F (fb ), associated with the ET reaction in the sequential ET followed by PT mechanism. (c) The 1D FE profile in the proton collective variable, F (fp ), associated with the PT reaction in the sequential ET followed by PT mechanism. The uni-molecular rate constant is calculated from the 1D FE profile along fb and the transmission coefficient associated with the concerted PCET reaction according to Eq. (2.8) (Sec. 2.4.1.1). CPET The FE profile used in the calculation of kuni is presented in Fig. 2.3a. As has been shown previously for concerted PCET in system-bath models, 30 the FE profile exhibits a sharp rise as 78 a function of fb associated with the formation of ring-polymer configurations, in which the electron spans the two iron sites (Fig. 2.3a, inset), followed by a more gradual change associated with reorientation of the solvent. The RPMD rate for the uni-molecular concerted PCET reaction is CPET CPET calculated to be kuni = (8 ± 2) × 103 s−1 , which yields a bi-molecular rate of kbi = (5 ± 2) × 102 M−1 s−1 ; the RPMD rate is in decent agreement with the experimentally measured PCET rate, kPCET = 5.8 ± 0.6 × 103 M−1 s−1 (Sec. 2.1). In principle, as stated in Sec. 2.1, there will be a separate rate for the two sequential PCET mechanisms corresponding to either sequential ET followed by PT (orange trajectories, Fig. 2.2) or sequential PT followed by ET (purple trajectories, Fig. 2.2). However, the two sequential PCET mechanisms in iron bi-imidazoline yield the same rate based on the symmetry of the system. Consequently, we focus solely on the calculation of the rate for the mechanism corresponding to sequential ET followed by PT. Figures 2.3(b) and 2.3(c) present the 1D FE profiles along fb and fp associated with the ET (Eq. (2.53)) and PT (Eq. (2.53)) reactions in the sequential ET followed by PT mechanism, respectively. The FE barrier along fp for the PT reaction is on the order of the thermal energy kB T at the experimental temperature, Fig. 2.3(c). Thus, the charge-separated intermediate is only short-lived, and the sequential ET followed by PT mechanism proceeds via rate-limiting ET. The rate constant for the sequential PCET mechanism is thus given by the forward rate for the sequential ET reaction. The rate constant for bi-molecular sequential ET reaction is given by an analogous expression to Eq. (2.60), 11,64,65 ET ET kbi = KA (r)kuni , (2.63) where KA (r) is given by Eq. 2.61, since the reactant species are the same for the sequential ET and concerted PCET reaction. The FE profile along fb in Fig. 2.3(b) along with the transmission ET coefficient for the sequential ET reaction are used to calculate the uni-molecular ET rate, kuni = 1.8 ± 0.8 s−1 , which yields a bi-molecular ET rate, and consequently a bi-molecular rate for the SPET sequential PCET mechanism, of kbi = (1.2 ± 0.5) × 10−1 M−1 s−1 . Comparison of the RPMD reaction rates for the concerted and sequential PCET mechanisms reveals that the concerted mechanism is the dominant mechanism in iron bi-imidazoline, validating the experimental result and analysis. 10 We note that the atomistic representation was not parameterized to obtain a specific relative favorability between the concerted and sequential mechanisms, and no a priori assumptions about the PCET mechanism are necessary for the use of RPMD. Having confirmed that the concerted mechanism is favored in iron bi-imidazoline, we now analyze 79 FeII(H2bim) + ! FeIII(Hbim)! FeIII(Hbim) + ! FeII(H2bim)! Figure 2.4: Reactive RPMD trajectories (red) for the concerted PCET reaction in the atomistic representation of iron bi-imidazoline reveal a Marcus-type solvent-gating mechanism indicated by the black arrows. The trajectories are projected onto the FE surface in the electron bead-count coordinate, fb , and the concerted PCET energy gap coordinate, ∆U , with contour lines indicating FE increments of 2 kcal/mol. Snapshots of the RPMD simulations from the reactant, transition state, and product regions of the concerted PCET reaction are also shown. the RPMD trajectories with respect to the solvent dynamics and the physical interactions that determine the PCET mechanism. Figure 2.4 presents the 2D FE profile along the electron bead-count coordinate, fb , and the concerted PCET energy gap coordinate, ∆U , computed for the concerted pathway as described in Sec. 2.4.1.4. The FE profile exhibits two basins of stability corresponding to the concerted PCET reactant and product species. Overlaid on the FE surface are representative samples from the ensemble of reactive RPMD trajectories (red) for the concerted PCET reaction. The RPMD trajectories exhibit a marcus-type solvent-gating mechanism, which has been seen previously for both ET 22 and PCET in system-bath models 30 , and is assumed in the derivation of PCET rate theories. 4,66–68 This mechanism is illustrated, from left to right, by the black arrows and the snapshots of the atomistic representation along the dynamical trajectories in Fig. 2.4. In the reactant basin, the electron is localized around the donor iron atom and the proton is bonded to the donor nitrogen; the solvent adopts a configuration that stabilizes the charge distribution of the reactant species. A solvent fluctuation then brings the system to configurations at which the vibronic diabats for the transferring electron and proton are nearly degenerate, ∆U ≈ 0, prior to the rapid transfer of both the electron and proton. In the transition state region of the concerted PCET reaction, the ring-polymer electron adopts a “kink-pair” configuration, in which the ring-polymer spans between both iron sites indicative of tunneling from reactants to products; the ring-polymer proton similarly 80 adopts a stretched configuration. Finally, the solvent relaxes, trapping the electron and proton on the acceptor molecule. de(t) 20 d(t) / D 10 dtot(t) dp(t) 0 dL(t) -10 -20 -300 -200 -100 0 100 200 300 t / fs Figure 2.5: The total dipole of the reactive complex projected along the vector connecting the two iron atoms, dtot (t), (black), as well as the individual components arising from the transferring electron, de (t), (blue) transferring proton, dp (t), (red), and the bi-imidazoline ligands that change charge as a function of the position of the transferring proton, dL (t), (green) for the concerted PCET reaction in the atomistic representation of iron bi-imidazoline. Figure 2.5 illustrates part of the system interactions that lead to the favorability of the concerted PCET mechanism in iron bi-imidazoline. The figure presents the total dipole of the reactive complex projected along the vector connecting the two iron atoms, dtot (t), along with the individual components arising from the transferring electron, de (t), transferring proton, dp (t), and the bi-imidazoline ligands, which change charge as a function of the position of the transferring proton, dL (t). These terms are calculated as * de (t) = * dp (t) = + ne −1 X (α) ex (t) ne α=1 np 1 X np γ=1 = (2.64) traj + qp (t)x (γ) * dL (t) , (t) , (2.65) traj + X k∈L qk (t)xk (t) , (2.66) traj and dtot (t) = de (t) + dp (t) + dL (t), (2.67) 81 where the atomic index k runs only over the nuclei in the bi-imidazoline ligands which change charge as a function of the position of the proton and h· · · itraj denotes the non-equilibrium ensemble average over the time-evolved reactive RPMD trajectories for concerted PCET (Sec. 2.4.1.8). The variable xk (t) is calculated as  xk (t) = QFeA − QFeD |QFeA − QFeD |  · (Qk − Qcom ) , (2.68) where Qcom is the position of the center of mass of the reactive species, which includes the nuclei in both iron complexes, the electron ring-polymer, and the proton ring-polymer. The variables x(α) and x(γ) are defined analogously to Eq. (2.68), where the position of the electron ring-polymer bead α, or the proton ring-polymer bead γ, is substituted for Qk , respectively. Figure 2.5 shows that |dp (t)| is small at all times, due to the small distance over which the proton transfers, while |de (t)| and |dL (t)| are large at all times. The orientation of de (t) and dL (t) are opposite in sign and switch on similar timescales. This cancellation of charge leads to the magnitude of the total dipole moment being smaller than the largest magnitude of the component dipoles at all times, i.e. |dtot (t)| < max (|de (t)| , |dp (t)| , |dL (t)|). Consequently, the degree to which the polar solvent couples to the reactive species during the concerted PCET reaction is reduced not by the direct cancellation of charge between the electron and proton, as was seen in system-bath models for PCET 30 , but by a cancellation of charge between the transferring electron and the time-dependent charge distribution of the bi-imidazoline ligands. Thus, there are two factors that contribute to the favorability of the concerted PCET mechanism in iron bi-imidazoline, (i) the previously stated reason that the concerted mechanism avoids the high-energy charge-separated intermediates and (ii) the concerted mechanism mitigates the high-energy barrier associated with solvent reorganization. The interactions in iron bi-imidazoline optimize both of these factors. 2.5.2 Controlling the dominant PCET Mechanism The previous section provided insight into the physical interactions that lead to the favorability of the concerted mechanism in iron bi-imidazoline. In this section we aim to make predictions as to how to vary the physical interactions of the PCET system in order to alter the dominant PCET mechanism. As we saw in the previous section, the solvent plays a large role in determining the dominant PCET mechanism. As such, we begin our analysis by investigating the effect of varying solvent polarity on the PCET mechanism in iron bi-imidazoline. Figure 2.6(a) presents the uni-molecular rates for the concerted and sequential mechanism in iron bi-imidazoline as a function of the polarity 82 (a) (c) 0.02 4 0.015 Concerted ΔGo / kcal mol-1 log(k / s-1) 3 2 Sequential 0.005 1 0 (b) 0 (d) 0.02 0.02 0.015 ΔG± / kcal mol-1 λ / kcal mol-1 0.015 0.01 0.005 0 0.01 0.01 0.005 4 5 6 0 4 dsolv / D 5 6 dsolv / D Figure 2.6: The (a) unimolecular rates, (b) outer-sphere reorganization energy, (c) driving force for the sequential ET prior to PT reaction, and (d) Marcus barrier associated with the concerted (black) and sequential (red) PCET mechanisms in iron bi-imidazoline (Systems 1a-1c) as a function of the polarity of the solvent. of the solvent (Systems 1a-1c); it is only necessary to compare the unimolecular PCET rates, since the pre-factors associated with the calculation of the bi-molecular rate (Eqs. (2.60) and (2.63)) are equivalent for the sequential and concerted PCET mechanism. The polarity of the solvent in the atomistic representation is quantified through the permanent dipole moment associated with the atomic charges on the acetonitrile molecule, dsolv = |qCH3 rCH3−C + qN rC−N | , (2.69) where qCH3 and qN are the atomic point charges associated with the methyl and nitrogen sites on the 83 solvent molecules, respectively, and rCH3−C and rC−N are the methyl-carbon and carbon-nitrogen bond distances, respectively; the atomic point charges vary by a multiplicative constant between Systems 1a-1c as described in Sec. 2.4.1. Though increasing the polarity of the solvent does not lead to a switch of the dominant PCET mechanism in iron bi-imidazoline, the relative favorability between the concerted and sequential PCET mechanisms is clearly affected by the solvent polarity; as the solvent polarity is increased, the sequential mechanism becomes more favorable in comparison to the concerted mechanism. This trend can be understood in terms of the Marcus barrier associated with the solvent reorganization energy for both the sequential and concerted PCET mechanism, ∆G‡ = (λo + ∆G◦ )2 , 4λ (2.70) where λo is the outer-sphere reorganization energy associated with either the concerted PCET reaction or the sequential ET reaction prior to PT; as before, it is sufficient to analyze only the sequential ET reaction, since the sequential PCET mechanism occurs via rate-limiting ET. The reorganization energy associated with the concerted and sequential PCET reactions are presented in Fig. 2.6(b). The reorganization energy associated with both mechanism increases with increasing solvent polarity, but due to the cancellation of charge during the concerted PCET reaction, the outer-sphere reorganization energy for concerted PCET is always lower than that for the sequential PCET mechanism. In addition, the reorganization energy associated with the sequential mechanism increases to a greater extent in comparison to the concerted mechanism, which would in principle lead to a decrease in the favorability of the sequential PCET mechanism with increasing solvent polarity. However, the increasing reorganization energy of the sequential mechanism is compensated by a decrease in the driving force of the sequential ET reaction with increasing solvent polarity, as seen in Fig 2.6(c); the driving force associated with the concerted reaction is always zero since the concerted reaction is symmetric. The decrease in driving force can be explained by an increased stabilization of the charge separated product following the sequential ET reaction. Figure 2.6(d) presents the total Marcus barrier associated with the concerted PCET mechanism (black) and sequential PCET mechanism (red), which accounts for both the solvent reorganization energy and the driving force. The Marcus barrier associated with the sequential mechanism is observed to decrease with increasing solvent polarity due to the decrease in driving force; the Marcus barrier associated with the concerted mechanism is observed to increase with increasing solvent polarity due to the increase in solvent reorganization energy. The trend in the Marcus barrier 84 6 4 log(k / s-1) Concerted 2 0 Sequential -2 -4 3 4 5 dsolv / D Figure 2.7: The uni-molecular rates for the concerted (black) and sequential (red) PCET reactions for systems with weakened electron-proton interactions (Systems 2a-2c) as a function of the polarity of the solvent. explains the increase in favorability of the sequential mechanism in comparison to the concerted mechanism in iron bi-imidazoline with increasing solvent polarity. A major contribution to the robustness of the concerted PCET mechanism in iron bi-imidazoline is the strength of interaction between the electron and proton, which leads to a large driving force for the sequential ET reaction, as observed in Fig. 2.3 and Fig. 2.6(c). An increase in solvent polarity is insufficient to overcome the intrinsically large penalty for separating the electron and proton during the sequential mechanism. As mentioned in the introduction, this electron-proton interaction can be experimentally quantified through the comparison of the redox potential between the protonated and deprotonated forms of the reactive species, which for the remainder of the paper we will define as ∆ε. 2,10 To thus attempt to observe a transition between the dominant PCET mechanism, we turn our attention to Systems 2a-2c, which model systems in which the electron-proton interaction is weakened (Sec. 2.4.1), as in systems in which ∆ε is small. Experimentally this can be achieved through the spatial separation of the electron and proton donor sites, 2,12,69,70 though in our simulations we model this interaction through the parameter ζcmplx . Figure 2.7 presents the uni-molecular rates for the concerted (black) and sequential (red) PCET mechanisms for Systems 2a-2c. The same qualitative trend is observed as in Fig. 2.6(a), in which the relative favorability of the sequential mechanism increases with increasing solvent polarity. However, for Systems 2a-2c, in which the electron-proton interaction is weakened, we can observe that the rates for the concerted and sequential PCET mechanisms become equivalent in highly polar solvents. 85 -17 log(k) / a.u. -18 -19 -20 -21 0.002 0.003 0.004 d / a.u. Figure 2.8: The reaction rates for the concerted (black) and sequential (red) PCET reactions for the system-bath model of PCET (Systems 3a-3e) as a function of solvent polarity. As such, we have observed two possible physical characteristics of a PCET system which can be varied to alter the dominant PCET mechanism, (i) the solvent polarity and (ii) the strength of the electron-proton interaction. We extend our analysis of the competition between concerted and sequential PCET by investigating a range of system-bath models for condensed phase PCET for which we have greater control over the parameters in comparison to the fully atomistic simulations. Figure 2.8 presents the rates for the concerted (black) and sequential (red) PCET mechanisms as a function of the solvent polarity (Systems 3a-3e). The solvent polarity is quantified as d = |µes rD + µps qpmin |, (2.71) where qpmin is the minimum of the proton potential Up in the system-bath model (Eq. (2.35)). Again, we observe the same qualitative trend as in Figs. 2.6 and 2.7, but now there is a complete transition between the concerted mechanism being the dominant mechanism at low solvent polarity and the sequential mechanism being the dominant mechanism at high solvent polarity. Figure 2.9 presents the rates for the concerted (black) and sequential (red) PCET mechanisms as a function of the strength of the electron-proton interaction, µep (Systems 4a-4e). As was discussed in terms of the atomistic simulations, here we can directly observe a transition in the dominant PCET mechanism as a function of the strength of the electron-proton interaction; the sequential mechanism dominates when the interaction is weak, as in when ∆ε would be small, and the concerted mechanism dominates when the interaction is strong, as in when ∆ε would be large. 86 log(k) / a.u. -18 -20 -22 0.8 1 1.2 1.4 1.6 1.8 µep Figure 2.9: The reaction rates for the concerted (black) and sequential (red) PCET reactions for the system-bath model of PCET (Systems 4a-4e) as a function of the strength of the electron-proton interaction, µep . log(k) / a.u. -18 -19 0.01 0.02 0.03 0.04 V0 Figure 2.10: The reaction rates for the concerted (black) and sequential (red) PCET reactions for the system-bath model of PCET (Systems 5a-5e) as a function of the height of the PT barrier, V0 . 87 We conclude by investigating one further physical interaction that can be varied to alter the PCET mechanism. Figure 2.10 presents rates for the concerted (black) and sequential (red) PCET mechanisms as a function of the height of the PT barrier, V0 . The concerted PCET mechanism is observed to be the dominant mechanism when the PT barrier is small, and the sequential PCET mechanism is observed to be the dominant mechanism when the PT barrier is large. This trend is explained by the fact that the sequential mechanism proceeds through a rate-limiting ET mechanism. Therefore, the sequential mechanism should be independent of the height of the PT barrier; deviations from a V0 independent sequential mechanism in Fig. 2.10 can be explained by the fact that varying the magnitude of V0 leads to a change in the equilibrium position of the proton, which in turn affects the electron-proton and proton-solvent interactions. On the other hand, the rate of the concerted mechanism depends strongly on the height of the PT barrier, decreasing with increasing V0 . In fact, the dependence of the PCET mechanism on the PT barrier height provides an additional reason as to why the concerted mechanism in iron bi-imidazoline is so robust; iron bi-imidazoline forms a short hydrogen bond (2.67 Å), which leads to a low PT barrier height. Thus there are three concrete ways of altering the PCET mechanism: (i) by varying solvent polarity, (ii) by varying electron-proton interaction, and (iii) by varying height of the PT barrier. 2.6 Conclusion In this work, we have utilized the RPMD method to investigate the competition between concerted and sequential PCET mechanisms in both fully atomistic and system-bath models of condensed phase PCET. The results presented here provide new insight into the physical interactions that govern the dominant PCET mechanism. We first investigated the competition between the concerted and sequential PCET mechanism in an atomistic representation of the organometallic catalyst iron bi-imidazoline. The thermal reaction rates for the competing PCET mechanisms were calculated using RPMD, and the concerted PCET mechanism was found to be the dominant mechanism in iron bi-imidazoline in agreement with previous experimental work. 10 Analysis of the RPMD trajectories illustrates that the concerted mechanism proceeds through a solvent-gating mechanism in which the magnitude of the solvent reorganization energy is mitigated due to a cancellation of charge during the concerted reaction; the cancellation of charge arises not from the direct interaction of the electron and proton, but instead through the interaction of the electron and the time-dependent charge distribution of the biimidazoline ligands. The favorability of the concerted mechanism in iron bi-imidazoline arises from 88 (i) the mitigation of the solvent reorganization energy, (ii) the avoidance of the charge-separated intermediates formed during the sequential mechanism, and (iii) the short hydrogen-bond distance which yields a low intrinsic barrier for PT. Based on our observations in iron bi-imidazoline we extend our analysis to investigate how to alter the dominant PCET mechanism in both atomistic and system-bath models of condensed-phase PCET. We find that increasing the solvent polarity leads to an increase in the favorability of the sequential mechanism in comparison to the concerted mechanism in iron bi-imidazoline. However, increasing the solvent polarity is not sufficient to observe a full transition between a dominant concerted versus sequential mechanism in iron bi-imidazoline due to the large interaction between the electron and proton. In comparison, a transition in the dominant PCET mechanism is observed when varying the solvent polarity for systems in which the electron-proton interaction is weaker, such as in systems in which the electron and proton donor sites are spatially separated. 2,12,69,70 In addition, directly varying the strength of the electron-proton interaction is also sufficient to alter the dominant PCET mechanism; the concerted mechanism dominates at large electron-proton interactions, while the sequential mechanism dominates at weak electron-proton interactions. Finally, we also observe a transition of the dominant PCET mechanism through increasing the height of the PT barrier during PCET; the concerted mechanism dominates at small PT barrier heights, and the sequential mechanism dominates at large PT barrier heights. As such, we observe three different methods for modifying the PCET system in order to control the dominant PCET mechanism, (i) the solvent polarity, (ii) the strength of the electron-proton interaction, and (iii) the height of the PT barrier. 89 Table 2.3: Coefficients for the electron potential (Eq. (2.34)) in the System bath model. a Coefficient Valuea aD bD cD in rD out rD 0.22663 2.69336 4.94560 -4.0 -8.0 in in The parameters for the acceptor Coulombic well are given by aA = aD , bA = −bD , cA = cD , rA = −rD , and out out rA = −rD . All values are reported in atomic units. Appendix A Coefficients for the electron potential This appendix presents the coefficients associated with the quadratic functions in Eq. (2.34). The coefficients are equivalent for all system-bath models and are provided in Table 2.3. The coefficients are chosen such that Ue (qe ) is continuous and its derivative is piecewise continuous over the full range of qe . Appendix B System-bath potential energy parameters This appendix describes parameters for the system-bath models for condensed-phase PCET. Table 2.4 presents the parameters for Systems 3a-3e, which model the varying polarity of the solvent, Systems 4a-4e, which model the varying strength of interaction between the electron and proton, and Systems 5a-5e, which model the varying PT barrier height, while maintaining the same minimum position of the intrinsic proton potential, Up (Eq. (2.35)). The parameters are chosen to be consistent with previously employed models for condensed phase PCET. 30 Appendix C Auxiliary restraining potentials This appendix describes auxiliary restraining potentials that are introduced for the PIMD sampling trajectories used in the calculation of 1D FE profiles and in the initial sampling of configurations for the RPMD trajectories both in the atomistic and system-bath models of PCET. These auxiliary restraining potentials simply prevent the PIMD sampling trajectories from visiting configurations outside of the sequential or concerted PCET channel of interest. 90 Table 2.4: Parameters for the system-bath model of PCET. Parametera Systems 3a-3e Systems 4a-4e Systems 5a-5e f ωs ωc M ms me mp η/M ωc rD rA µe ωp V0 µes µps µep 12 2.3221 × 10−4 2.3221 × 10−4 25539 25539 1.0 1836.1 1.0 -6.0 6.0 3.668 1.47 × 10−2 3.00 × 10−2 −7.80 × 10−4 − −1.26 × 10−3 3.90 × 10−3 − 6.30 × 10−3 1.1 12 2.3221 × 10−4 2.3221 × 10−4 25539 25539 1.0 1836.1 1.0 -6.0 6.0 3.668 1.47 × 10−2 3.00 × 10−2 −1.20 × 10−3 6.00 × 10−3 1.0 - 1.5 12 2.3221 × 10−4 2.3221 × 10−4 25539 25539 1.0 1836.1 1.0 -6.0 6.0 3.668 6.00 × 10−3 − 1.59 × 10−2 1.50 × 10−2 − 3.50 × 10−2 −1.20 × 10−3 6.00 × 10−3 1.1 a 2.C.1 All values are reported in atomic units. Auxiliary restraining potential for concerted PCET in the atomistic models For Systems 1a-1c and 2a-2c, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the concerted channel.This potential is given by Uaux (fp , fb ) =  2  aaux [fp − q+ (fb )] ,           2 aaux [fp − q− (fb )] ,            0, fp > q+ (fb ) fp < q− (fb ) (2.72) otherwise where q+ (fb ) = baux fb + caux (2.73) q− (fb ) = baux fb − caux . (2.74) and The coefficients aaux , baux , and caux (Table 2.5) are chosen to restrict the system to the concerted channel. 91 Table 2.5: Parameters for the auxiliary restraining potential in Eq. 2.72. System 1a 1b 1c 2a 2b 2c a b aaux a baux b caux b 1000 1000 1000 1000 1000 1000 0.50 0.50 0.50 0.48 0.48 0.47 0.33 0.33 0.33 0.21 0.19 0.15 aaux is given in units of kcal/mol Å−2 . baux and caux are given in units of Å. Table 2.6: Parameters for the auxiliary restraining potential in Eq. 2.75. System 1a 1b 1c 2a 2b 2c a 2.C.2 aSET a bSET b 1000 1000 1000 1000 1000 1000 0.25 0.21 0.20 0.15 0.08 0.08 aSET is given in units of kcal/mol Å−2 . b bSET is given in units of Å. Auxiliary restraining potential for ET prior to PT in the atomistic models For Systems 1a-1c and 2a-2c, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the ET channel in the sequential mechanism corresponding to the reaction given in Eq. 2.53. This potential is given by USET (fp ) =  2   aSET (fp − bSET ) ,    fp < bSET      0, otherwise. (2.75) The coefficients aSET a.u. and bSET are chosen to correctly restrict the system to the ET channel and are provided in Table 2.6. 2.C.3 Auxiliary restraining potential for concerted PCET in the systembath models For Systems 3a-3e, 4a-4e, and 5a-5e, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the concerted channel in the system-bath model. 92 Table 2.7: Parameters for the auxiliary restraining potential in Eq. 2.76. a System aaux baux caux a 3a 3b 3c 3d 3e 4a 4b 4c 4d 4e 5a 5b 5c 5d 5e 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.57 0.56 0.56 0.55 0.55 0.55 0.55 0.56 0.57 0.58 0.57 0.56 0.55 0.55 0.55 0.32 0.32 0.32 0.33 0.33 0.32 0.34 0.35 0.36 0.38 0.33 0.33 0.33 0.33 0.33 aaux is given in units of a.u. ×10−2 ; all other parameters are given in atomic units. This potential is given by Uaux (q̄p , fb ) =  2  aaux [q̄p − q+ (fb )] ,           2 aaux [q̄p − q− (fb )] ,            0, q̄p > q+ (fb ) q̄p < q− (fb ) (2.76) otherwise where q+ (fb ) = baux fb + caux (2.77) q− (fb ) = baux fb − caux . (2.78) and The coefficients aaux , baux , and caux (Table 2.7) are chosen to restrict the system to the concerted channel. 2.C.4 Auxiliary restraining potential for ET prior to PT in the systembath models For Systems 3a-3e, 4a-4e, and 5a-5e, we now discuss the auxiliary restraining potential introduced to restrict equilibrium sampling of the system to the ET channel in the sequential mechanism. This 93 Table 2.8: Parameters for the auxiliary restraining potential in Eq. 2.79. System 1a 1b 1c 1d 1e 2a 2b 2c 2d 2e 3a 3b 3c 3d 3e a aSET a bSET b 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.29 0.29 0.29 0.29 0.29 0.20 0.20 0.20 0.23 0.26 0.30 0.30 0.30 0.30 0.30 aSET is given in units of a.u. ×10−2 . b bSET is given in units of a.u. potential is given by USET (q̄p ) =  2   aSET (q̄p − bSET ) ,    q̄p < bSET      0, otherwise. (2.79) The coefficients aSET a.u. and bSET are chosen to correctly restrict the system to the ET channel and are provided in Table 2.8. 94 References [1] J. L. Dempsey, J. R. Winkler, and H. B. Gray, Chem. Rev. 110, 7024 (2010). [2] J. J. Warren, T. A. Tronic, and J. M. Mayer, Chem. Rev. 110, 6961 (2010). [3] S. Y. Reece and D. G. Nocera, Annu. Rev. Biochem. 78, 673 (2009). [4] S. Hammes-Schiffer and A. A. Stuchebrukhov, Chem. Rev. 110, 6939 (2010). [5] S. Y. Reece, J. M. Hodgkiss, J. Stubbe, and D. G. Nocera, Phil. Trans. R. Soc. B 361, 1351 (2006). [6] R. I. Cukier and D. G. Nocera, Annu. Rev. Phys. Chem. 49, 337 (1998). [7] M. H. V. Huynh and T. J. Meyer, Chem. Rev. 107, 5004 (2007). [8] J. M. Mayer and I. J. Rhile, Biochim. Biophys. Acta. Bioenerg. 1655, 51 (2004). [9] J. P. Roth, S. Lovell, and J. M. Mayer, Inorg. Chem. 38, 2760 (1999). [10] J. P. Roth, S. Lovell, and J. M. Mayer, J. Am. Chem. Soc. 122, 5486 (2000). [11] N. Iordanova, H. Decornez, and S. Hammes-Schiffer, J. Am. Chem. Soc. 123, 3723 (2001). [12] J. J. Gray, Warren, A. R. Menzeleev, and J. M. Mayer, J. S. Kretchmer, T. J. Phys. Chem. Lett. F. Miller III, H. B. (2013), 10.1021/jz400029w, http://pubs.acs.org/doi/pdf/10.1021/jz400029w . [13] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 (2004). [14] R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals, International series in pure and applied physics (McGraw-Hill, New York,, 1965) p. xiv. [15] D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 (1981). [16] R. Collepardo-Guevara, Y. V. Suleimanov, and D. E. Manolopoulos, J. Chem. Phys. 130, 174713 (2009). 95 [17] N. Boekelheide, R. Salomón-Ferrer, and T. F. Miller III, Proc. Natl. Acad. Sci. 108, 16159 (2011). [18] R. Collepardo-Guevara, I. R. Craig, and D. E. Manolopoulos, J. Chem. Phys. 128, 144502 (2008). [19] Y. V. Suleimaov, R. Collepardo-Guevara, and D. E. Manolopoulos, J. Chem. Phys. 134, 044131 (2011). [20] T. F. Miller III, J. Chem. Phys. 129, 194502 (2008). [21] A. R. Menzeleev and T. F. Miller III, J. Chem. Phys. 132, 034106 (2010). [22] A. R. Menzeleev, N. Ananth, and T. F. Miller III, J. Chem. Phys. 135, 074106 (2011). [23] T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. 123, 154504 (2005). [24] T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. 122, 184503 (2005). [25] S. Habershon, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 131 (2009). [26] I. R. Craig and D. E. Manolopoulos, Chem. Phys. 322, 236 (2006). [27] T. E. Markland, S. Habershon, and D. E. Manolopoulos, J. Chem. Phys. 128, 194506 (2008). [28] Y. V. Suleimanov, J. Phys. Chem. C 116, 11141 (2012). [29] T. E. Markland, J. A. Morrone, K. Miyazaki, B. J. Berne, D. R. Reichman, and E. Rabani, J. Chem. Phys. 136, 074511 (2012). [30] J. S. Kretchmer and T. F. Miller, III, J. Chem. Phys. 138, 134109 (2013). [31] R. P. Steele, J. Zwickl, P. Shushkov, and J. C. Tully, J. Chem. Phys. 134, 074112 (2011). [32] E. Wigner, Phys. Chem. Abt. B 19, 203 (1932). [33] H. Eyring, J. Chem. Phys. 3, 107 (1935). [34] J. C. Keck, J. Chem. Phys. 32, 1035 (1960). [35] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 122, 084106 (2005). [36] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 123, 34102 (2005). [37] W. H. Miller, J. Chem. Phys. 58, 1664 (1973). 96 [38] D. Chandler, J. Chem. Phys. 68, 2959 (1978). [39] C. H. Bennet, in Algorithms for Chemical Computations, edited by R. E. Christofferson (American Chemical Society, Washington, DC, 1977) p. 63. [40] E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989). [41] G. K. Schenter, B. C. Garrett, and D. G. Truhlar, J. Chem. Phys. 119, 5828 (2003). [42] J. B. Watney, A. V. Soudackov, K. F. Wong, and S. Hammes-Schiffer, Chem. Phys. Lett. 25, 268 (2006). [43] M. Parrinello and A. Rahman, J. Chem. Phys. 80, 860 (1984). [44] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comp. Chem. 25, 1157 (2004). [45] E. Guardia, R. Pinzon, J. Casulleras, M. Orozco, and F. J. Luque, Mol. Simulat. 26, 287 (2001). [46] C. M. Breneman and K. B. Wiberg, J. Comp. Chem. 11, 361 (1990). [47] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, . Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, “Gaussian?09 Revision D.01,” Gaussian Inc. Wallingford CT 2009. [48] H. Azzouz and D. Borgis, J. Chem. Phys. 98, 7361 (1993). [49] A. Warshel, Computer modeling of chemical reactions in enzymes and solutions (John Wiley, New York, 1991). 97 [50] R. A. Kuharski, J. S. Bader, D. Chandler, M. Sprik, M. L. Klein, and R. W. Impey, J. Chem. Phys. 89, 3248 (1988). [51] D. F. Coker, B. J. Berne, and D. Thirumalai, J. Chem. Phys. 89, 5689 (1987). [52] A. O. Caldeira and A. J. Leggett, Annals of Physics 149, 374 (1983). [53] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987). [54] M. Topaler and N. Makri, J. Chem. Phys. 101, 7500 (1994). [55] W. Smith and T. R. Forester, J. Mol. Graph. 14, 136 (1996). [56] L. Verlet, Phys. Rev. 159, 98 (1967). [57] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, Annu. Rev. Phys. Chem. 64, 387 (2013). [58] C. L. Brooks, B. M. Pettitt, and M. Karplus, J. Chem. Phys. 83, 5897 (1983). [59] B. D. Raedt, M. Sprik, and M. L. Klein, J. Chem. Phys. 80, 5719 (1984). [60] H. C. Andersen, J. Comput. Phys. 52, 24 (1983). [61] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002). [62] G. Hummer, J. Chem. Phys. 120, 516 (2004). [63] E. W. Vanden-Eijnden, J. Stat. Phys. 123, 503 (2006). [64] M. D. Newton and N. Sutin, Ann. Rev. Phys. Chem. 35, 437 (1984). [65] W. Smith and T. R. Forester, Prog. Inorg. Chem. 30, 441 (1983). [66] R. I. Cukier, J. Phys. Chem. 98, 2377 (1994). [67] R. I. Cukier, J. Phys. Chem. 100, 15428 (1996). [68] A. Soudackov and S. Hammes-Schiffer, J. Chem. Phys. 113, 2385 (2000). [69] V. W. Manner, A. G. DiPasquale, and J. M. Mayer, J. Am. Chem. Soc. 130, 7210 (2008). [70] V. W. Manner and J. M. Mayer, J. Am. Chem. Soc. 131, 9874 (2009). 98 Chapter 3 The role of nuclear dynamics during electron transfer in the context of kinetically-constrained ring polymer molecular dynamics 3.1 Introduction Non-adiabatic reactions are ubiquitous throughout chemistry and biology, including such examples as charge transfer, energy transfer, and non-radiative decay. The accurate and efficient simulation of non-adiabatic processes in the condensed phase remains an ongoing challenge for theoretical methods. Recently we have developed an extension to the ring-polymer molecular dynamics method (RPMD), kinetically-constrained (KC) RPMD, 1 which allows for the treatment of general multielectron, non-adiabatic processes. In the current study, we explore the performance of KC-RPMD for the investigation of model condensed phase electron transfer (ET) reactions in which the relative timescale between the electronic and various nuclear degrees of freedom alters both the ET reaction rate and mechanism. Analogous to RPMD in the usual position representation, 2,3 KC-RPMD is an approximate quantum dynamics method based on the path-integral formulation of statistical mechanics. 4 However, in KC-RPMD the path-integral discritization is performed in the electronic state basis instead of in the usual position representation. KC-RPMD further employs a coarse graining procedure to reduce the discrete electronic state variables to a single continuous coordinate, and a “kinetic constraint” modification of the equilibrium distribution to address the known failures of path-integral based estimates for tunneling rates. This kinetically-constrained distribution is rigorously preserved using continuous equations of motion, yielding a real-time model for the non-adiabatic dynamics of a quantum system. The KC-RPMD equations of motion preserve the useful features of conventional position-represention RPMD such as detailed balance, time-reversal symmetry, and invariance of reaction rate calculations to the choice of dividing surface. In addition, KC-RPMD allows for the simulation of multi-electron and photochemical systems, which cannot be treated with conventional 99 RPMD. Previously, KC-RPMD was successfully employed to investigate a range of model non-adiabatic reactions, including a simple avoided-crossing reaction and ET reactions in both the normal and inverted regimes, and in the non-adiabatic and adiabatic regimes. 1 In this study, we explore the performance of the method in modeling condensed phase ET reactions, in which the dynamics of the nuclear degrees of freedom strongly affect the ET rate and mechanism. Specifically, we investigate the ET rate as a function of the solvent friction and as a function of the frequency of a “donor-acceptor” mode that modulates the magnitude of the electronic-coupling. In doing so, we revisit several technical aspects of the KC-RPMD derivation including the kinetic-constraint and the choice of mass associated with the continuous electronic coordinate. In particular, we correct for the presence of spurious resonances in the electronic variable by coupling the electronic variable to a Langevin bath, analogous to the treatment of spurious resonances associated with the internal modes of the ring-polymer in conventional RPMD. 5 The present work demonstrates the ability of KC-RPMD to simulate more complex systems in which dynamical effects play a strong role in determining the reaction rate and mechanism. 3.2 Theory In this section we summarize the derivation of the KC-RPMD method, including the introduction of a Langevin thermostat that is coupled to the auxiliary variable. Consider a general, two-level system in the diabatic representation with a Hamiltonian operator of the form Ĥ = T̂ + V̂ , where T̂ = d X p2j 2mj j=1 (3.1) describes the kinetic energy of a system of d nuclear degrees of freedom and   V0 (R) K(R)  V̂ (R) =   K(R) V1 (R) (3.2) is the potential energy in the diabatic representation as a function of the nuclear coordinates, R. The canonical partition function for the two-level system is Z = Tr[e −β Ĥ Z ]= dR X hR, i|e−β Ĥ |R, ii, i=0,1 (3.3) 100 where β = (kB T )−1 is the inverse temperature and i denotes the diabatic electronic state. By (a) resolving the identity in the product space of the electronic state and nuclear position coordinates and (b) employing the appropriate short-time approximation, we discretize the trace into the ringpolymer representation with n beads and obtain the familiar result 1,6 Z Zn = d{R(α) } X (α) ρRP }, {i(α) }), n ({R (3.4) {i(α) } such that Z = lim Zn . The nuclear position and electronic state of the αth ring-polymer bead  n→∞ are given by R(α) , i(α) and the usual periodic constraint of the ring polymer is satisfied by     R(n+1) , i(n) = R(1) , i(1) . The ring polymer distribution is given by (α) (α) ρRP },{i(α) }) = Ωe−βUint ({R n ({R }) n Y Mi(α) ,i(α+1) (R(α) ), (3.5) α=1 where Ω = n/2 d  Y nmj j=1 2π~2 β . The internal ring-polymer potential is given by n Uint ({R(α) }) = d  2 1 XX (α) (α+1) mj ωn2 Rj − Rj , 2n α=1 j=1 (3.6) where ωn = (βn ~)−1 . The term Mi,i0 (R) denotes the i, i0 element of the matrix   M(R) = e −βn V0 (R) −βn K(R)e −βn K(R)e −βn V1 (R) e −βn V0 (R) −βn V1 (R)   , (3.7) where βn = β/n. Equations (3.4)-(3.7) define the usual path-integral partition function for a twolevel system in the diabatic representation. To obtain the KC-RPMD distribution and equations of motion from the usual path-integral distribution, Eqs. (3.4)-(3.7), we first introduce a discrete collective variable that reports on the existence of kink-pairs in the ring-polymer configuration    −1, i(α) = 0 for all α,    θ({i(α) }) = 1, i(α) = 1 for all α,      0, otherwise. (3.8) We then introduce a continuous dummy variable y that is tethered to θ({i(α) }) via a square-well 101 restraining potential Vr (y, {i(α) }), such that (α) e−βVr (y,{i }) = f (y, θ({i(α) })), (3.9) where 1 f (y, θ) = lim b→∞ 2lθ and lθ =        lθ 1 − tanh b |y − θ| − , 2 L, θ({i(α) }) = 0 (3.10) (3.11)   2 − L, otherwise. Equation (3.10) defines a set of square restraining potentials of width 2-L, 2-L, and L centered at y = −1, y = 1, and y = 0, respectively, and constitutes the first difference between the current and original formulation of KC-RPMD. Now, the relative width between the square restraining potential centered at y = 0 versus centered at y = −1 or y = 1 can be varied. It should be noted, though, that the kinetically-constrained distribution presented below is invariant to the choice of L. The kinetically-constrained distribution is obtained by reducing the ring-polymer probability distribution in Eq. (3.5) with respect to the discrete electronic variables and by introducing a kinetic constraint that penalizes the formation of kink-pairs at ring-polymer configurations for which the diabatic surfaces are non-degenerate, such that KC (α) ρKC }, y) = Ωe−βVeff ({R n ({R (α) },y) , (3.12) where   n X Y (α) 1 KC g({i(α) }, {R(α) })e−βVr (y,{i }) Mi(α) ,i(α+1)(R(α) ), Veff ({R(α) }, y) = Uint ({R(α) }) − ln β (α) α=1 {i and g({i(α) }, {R(α) }) =       } 1, i(α) = 0 for all α, 1, i(α) = 1 for all α,   1      a 2 ηe−a(w(R̄))2 , otherwise π (3.13) is the penalty function. The function w(R) = (V0 (R) − V1 (R)) /K(R) is the scaled difference in n 1 X (α) the diabatic potential surfaces, R̄ = R is the ring-polymer centroid coordinate, and a is a n α=1 unitless convergence parameter chosen sufficiently large to converge the free energy (FE) of kink-pair 102 formation given by ∆F KC = F KC (0) − F KC (−1), where F KC (y) = − Z 1 (α) ln d{R(α) }ρKC }, y). n ({R β (3.14) The parameter η is a multiplicative factor that is chosen to avoid biasing the probability of kink-pair formation at nuclear configurations for which the diabats cross, such that η = 2π h|K(R)|ic |F0 (R) − F1 (R)| |K(R)|2 c , h|F0 (R) − F1 (R)|ic h|K(R)|3 ic (3.15) where F0 (R) and F1 (R) are the force vectors associated with electronic state 0 and 1, respectively. The brackets denote an ensemble average constrained to the intersection of the diabatic surfaces, such that R h(...)ic = dRδ(V0 (R) − V1 (R))(...)e−βV0 (R) R . dRδ(V0 (R) − V1 (R))e−βV0 (R) (3.16) The functional form of Eq. 3.15 differs from Ref. 1, as described in the detailed derivation in App. 3.1. (α) The classical equations of motion associated with the equilibrium distribution ρKC }, y) n ({R are (α) v̇j =− n ∂ V KC ({R(α) }, y) mj ∂R(α) eff j 1 ∂ KC v̇y = − V ({R(α) }, y). my ∂y eff (3.17) These equations exactly preserve the well-defined equilibrium distribution, Eq. 3.12, and as such preserve all of the essential features of conventional RPMD, including detailed balance, timereversibility, invariance of the thermal rate constant on the choice of dividing surface, and the ability to harness the full machinery of classical molecular dynamics techniques. The equations of motion defined in Eq. 3.17 utilize the physical masses for the nuclear degrees of freedom mj . The mass of the auxiliary variable, my , is chosen such that the KC-RPMD transition state theory (TST) rate exactly recovers the Fermi-Golden rule rate 7 kFGR = 2π ~ R dRδ(V0 (R) − V1 (R))|K(R)|2 e−βV0 (R) R . dRδ(V0 (R) − V1 (R))e−βV0 (R) (3.18) Although the derivation of my is analogous, the choice of the Fermi-Golden rule rate as the reference rate theory differs from the original formulation of KC-RPMD, 1 in which the Landau-Zener TST 103 rate expression for non-adiabatic transitions was utilized. 8 We note that the Fermi-Golden rule and Landau-Zener TST rates differ only by a factor of two, but the Fermi-Golden rule expression is more appropriate for the calculation of thermal rate constants in the condensed phase. 7 The resulting expression for the mass of the auxiliary variable is ~2 β 3 my = 2πL2 " h|K(R)|ic |F0 (R) − F1 (R)| |K(R)|2 c h|F0 (R) − F1 (R)|ic h|K(R)|2 ic #2 . (3.19) It is important to take note that the units of my are mass×distance2 /[L2 ]. However, since the auxiliary variable is unitless, L is also unitless, and the units of my are simply mass×distance2 . Up to this point, the formulation of KC-RPMD has been analogous to the original formulation, 1 though now the square restraining potential, f (y, θ) in Eq. 3.10, the multiplicative factor, η in Eq. 3.15, and the mass of the auxiliary variable, my in Eq. 3.19, have different functional forms. Now, however, we discuss the introduction of a Langevin thermostat that is coupled to the auxiliary variable. The auxiliary variable and the nuclear degrees of freedom are only weakly coupled during the dynamics generated by Eq. 3.17. As such, the timescale for thermalization of the auxiliary variable is exceedingly long. This behavior leads to spurious resonances in the auxiliary variable, in which the auxiliary variable transitions between the reactant and product basins without incurring the statistical penalty of kink-formation. To alleviate this issue we follow a protocol analogous to the fix of the spurious resonances associated with the internal modes of the ring polymer observed in the usual form of RPMD. 5 Specifically, we introduce a Langevin thermostat that is coupled to the auxiliary variable such that the KC-RPMD equations of motion are now (α) v̇j =− n ∂ V KC ({R(α) }, y) mj ∂R(α) eff j 1 ∂ KC v̇y = − V ({R(α) }, y) − γy vy + my ∂y eff r (3.20) 2γy my ℵ(t), β where γy is the friction coefficient and ℵ(t) is a Gaussian random force of unit width. The friction coefficient is chosen such that the time-scale for thermalization of the auxiliary variable is the same as the time-scale for the auxiliary variable to cross the reactant basin, such that 1 γy = 2(2 − L) v h u u 1 − 2 ln 2−L t L
a 1/2 ηβ 2 π βmy The detailed derivation of Eq. 3.21 is given in App. 3.2. i − 4hln |K(R)|ic . (3.21) 104 We conclude this section by defining the KC-RPMD reaction rate coefficient. 1 The KC-RPMD rate is calculated using standard methods 9–12 and is separated into a statistical and a dynamical contribution as 13,14 KC−RPMD k KC−RPMD = kTST lim κ(t), (3.22) t→∞ KC−RPMD where kTST is the TST estimate for the rate associated with the dividing surface ξ(r) = ξ ‡ , and κ(t) is the time-dependent transmission coefficient that corrects for dynamical recrossing at the dividing surface. Here, ξ(r) is a collective variable that distinguishes between reactant and product basins of stability, defined as a function of the position vector of the full system in the ring-polymer n o representation, r = {R(α) }, y . The KC-RPMD TST rate is calculated using 3 ‡ 1 KC−RPMD kTST =√ ‡ 2πβ e−β∆F (ξ ) hχξ i R ξ‡ −∞ dξe−β∆F (ξ) , (3.23) where F (ξ) is the FE along ξ relative to a reference value ξ ◦ , such that ‡ e−β∆F (ξ ) = hδ(ξ(r) − ξ ‡ )i , hδ(ξ(r) − ξ ◦ )i and 15–17  nd+1 X χξ (r) =  j 1 mj  ∂ξ(r) ∂rj 2 (3.24) 1/2  . (3.25) The sum in Eq. 3.25 runs over all the nd + 1 degrees of freedom for the ring-polymer representation used here, and mj denotes the mass associated with each degree of freedom. The angle brackets indicate an equilibrium ensemble average R h. . . i = where v = dr R dv e−βH(r,v) (. . . ) R , dr dv e−βH(r,v) R (3.26) nn o o v (α) , vy is the velocity vector for the full system in the ring-polymer representation and H(r, v) is the ring-polymer Hamiltonian associated with the KC-RPMD effective potential. The ensemble average ‡ h. . . i = R dr R dv e−βH(r,v) δ(ξ(r) − ξ ‡ )(. . . ) R dr dv e−βH(r,v) δ(ξ(r) − ξ ‡ ) R (3.27) corresponds to the ensemble average constrained to the dividing surface. The transmission coefficient 105 in Eq. 3.22 is calculated as
hξ˙0 h ξ(rt ) − ξ ‡ i‡ , κ(t) = hξ˙0 h(ξ˙0 )i‡ (3.28) where h(x) is the Heaviside function, and the subscripts 0 and t denote evaluation of the quantity from the trajectory at its initiation and after evolution for time t, respectively. 3.3 Model Systems In this paper we investigate sets of system-bath models that describe condensed phase ET in a wide-range of regimes, in which the potential energy function takes the form V̂ (R) = V̂S (R) + 1VB (R), (3.29) where V̂S (R) defines the potential of the redox system and VB (R) defines the potential of the bath. We use atomic units throughout, unless otherwise noted. System A models a condensed phase ET reaction with constant coupling, for which 18   1 2 2  2 ms ωs (s − s0 ) V̂ (s) =  K K  , 1 ms ωs2 (s − s1 )2 +  2 (3.30) where s corresponds to the local solvent dipole. This solvent coordinate is linearly coupled to a set of f bath modes, x, such that VB (s, x) =  f X  1 M ωj2 2 j=1 cj s xj − M ωj2 !2  , (3.31) and M denotes the mass of the bath modes. The bath exhibits an Ohmic spectral density with cutoff frequency ωc , J(ω) = γωe−ω/ωc , (3.32) where γ denotes the friction coefficient of the bath and controls the strength of coupling between the system and the bath modes. The spectral density in Eq. 3.32 is discretized into f oscillators with frequencies 19  ωj = −ωc ln j − 0.5 f  (3.33) 106 Table 3.1: Parameters for System A. Parameter System A1a System A2a SystemA3b ms ωs s0 s1  K M ωc γ/M ωc f T my η 1836 2.28 × 10−3 -2.40 2.40 0 − 0.236 6.67 × 10−7 1836 2.28 × 10−3 1.0 12 300 K 8.26 × 10−3 6.28 1836 2.28 × 10−3 -2.40 2.40 0 1 × 10−4 1836 2.28 × 10−3 0.3 − 300.0 12 300 K 1.86 × 102 6.28 1836 2.28 × 10−3 -1.198 1.198 0 2.28 × 10−4 1836 2.28 × 10−3 0.013 − 30.0 12 460 K 2.68 × 102 6.28 a The parameters for this System are taken from Ref. 1. The parameters for this System are taken from Ref. 20. Unless otherwise noted, values are reported in atomic units. b and coupling constants  cj = ωj 2γM ωc fπ 1/2 , (3.34) where j = 1 . . . f . Three sets of parameters, System A1, System A2, and System A3, are used to model ET over a wide-range of regimes. System A1 models the transition from the normal to inverted regime of ET. System A2 models the transition from the Marcus regime to the Zusman regime of ET as a function of the friction coefficient, γ. System A3 models the transition from the weakly dissipative regime at low friction to the Marcus regime at intermediate values of the friction coefficient. The parameters for System A are presented in Table 3.1. The quantities η and my are analytically evaluated from Eqs. (3.15) and (3.19) and are also presented in Table 3.1. System B extends System A to include an additional degree of freedom, q, which models the fluctuating distance between the electron donor and acceptor. The coupling then exponentially depends on q, such that  1 2 2  2 ms ωs (s − s0 ) + V (q) V̂ (s) =  K(q)  K(q)  , 1 2 2 ms ωs (s − s1 ) + V (q) +  2 (3.35) 107 Table 3.2: Parameters for System B. Parameter Range of Values mq ωq q∗ q σq bq K∗  my η 51039 4 × 10 − 3 × 10−3 12.0 0.1 8.0 3.0 6.67 × 10−7 0 2.35 × 10−2 − 8.41 × 10−3 0.77 − 6.05 −4 All values are reported in atomic units. where V (q) VHW (q) 1 mq ωq2 (q − q ∗ )2 + VHW (q), 2 "    6 # 12  σq σq   4q − + q , r < 21/6 σq q q =    0, r ≥ 21/6 σq , = (3.36) (3.37) and ∗ K(q) = K ∗ e−bq (q−q ) . (3.38) The parameters for System B are the same as for System A1 except for those noted in Table 3.2. The additional parameters necessary to define System B, and the calculated values of my and η, are also provided in Table 3.2. The parameters are chosen to model ET between two iron atoms solvated in water. 21 3.4 Calculation Details In all simulations, the KC-RPMD equations of motion, Eq. (3.20), of the nuclei are evolved using the velocity Verlet algorithm; 22 the equations of motion of the auxiliary variable, y, are evolved using an extension of the velocity Verlet algorithm to include the Langevin bath. 23 In all simulations, the nuclei are treated classically by restricting the position of the ring-polymer beads for each nuclei to coincide. The coefficient b was found to be converged for all calculations with a value of b = 400. A value of the parameter L = 0.1 was found to yield reasonable numerical simulations when including the Langevin bath connected to the auxiliary variable. In all cases the TST dividing surface is 108 defined as an isosurface of the auxiliary variable, y. 3.4.1 KC-RPMD rate calculation in System A1 In System A1, the driving force parameter, , is varied to model the transition between the normal and inverted regime of ET, such that  ∈ {0, 0.04, 0.07, 0.118, 0.18, 0.236} . For these values of , the calculations were found to converge with values of − log(a) ∈ {−5.3, −5.3, −5.3, −5.3, −5, −5}, respectively. The ring polymer is discretized using n = 32 beads for all values of . The KC-RPMD rates are obtained from the product of the TST rates, Eq. (3.23), and the transmission coefficients, Eq. (3.28). The KC-RPMD TST rates are obtained from F (y), the FE profile in the continuous auxiliary variable. The FE profile is obtained by direct numerical integration. The transmission coefficients are calculated using KC-RPMD trajectories that are released from the dividing surface associated with y ‡ = 0. The values of the mass of the auxiliary variable, my , denoted in Table 3.1 are small in comparison to the mass of the nuclei, which in principle would necessitate an exceedingly small time step to numerically integrate the equations of motion. However, we note that any choice of the value of my that is still small in comparison to the mass of the nuclei, but allows for a larger time step, will yield the same final value for the transmission coefficient. We thus choose a value of my = 10.0 for all values of  when integrating the KC-RPMD equations of motions for the trajectories used to calculate the transmission coefficient; it is important to note that we only use a value of my = 10.0 for the calculation of the dynamical trajectories, but use the true values of my denoted in Table 3.1 in the calculation of the TST rate. The value of the friction coefficient of the Langevin bath connected to y is calculated from Eq. (3.21) using the same value of my = 10.0 as used in the dynamical KC-RPMD trajectories, such that γy = 1.49×10−2 for all values of . We have confirmed that choosing a smaller value of my yields numerically the same results as a value of my = 10.0. For each value of the driving force , a total of 1000 trajectories are released. Each trajectory is evolved for 200 fs using a time step of dt = 0.002 fs with initial velocities sampled from the Maxwell-Boltzmann distribution. The initial configurations for the KC-RPMD trajectories are generated from long KC-RPMD trajectories that are constrained to the dividing surface using the RATTLE algorithm; 24 the constrained trajectories are 200 ps in time and are thermostatted by resampling the velocities from the Maxwell-Boltzmann distribution every 100 fs. 109 3.4.2 KC-RPMD rate calculation in System A2 In System A2, the friction coefficient of the bath, γ, is varied to model the transition from the Marcus to the Zusman regime of ET, such that γ/M ωc ∈ {0.3, 1, 3, 10, 30, 100, 200, 300}. For all values of γ the calculations were found to converge with a value of − log(a) = −1.5 and a value of n = 128. The KC-RPMD TST rates are obtained by numerical integration of F (y). The transmission coefficients are calculated using KC-RPMD trajectories that are released from the dividing surface associated with y ‡ = 0. The true value of the mass of the auxiliary variable, my = 1.86 × 102 , is used for the dynamical KC-RPMD trajectories and for the calculation of γy , such that γy = 1.34 × 10−3 for all values of γ. For each value of the friction coefficient γ, a total of 8000 trajectories are released. Each trajectory is evolved for 200 fs using a time step of dt = 0.005 fs with initial velocities sampled from the Maxwell-Boltzmann distribution. The initial configurations for the KC-RPMD trajectories are generated from long KC-RPMD trajectories that are constrained to the dividing surface using the RATTLE algorithm; 24 the constrained trajectories are 1.6 ns in time and are thermostatted by resampling the velocities from the Maxwell-Boltzmann distribution every 100 fs. 3.4.3 KC-RPMD rate calculation in System A3 In System A3, the friction coefficient of the bath, γ, is varied to model the transition from the weakly dissipative to Marcus regime of ET, such that γ/M ωc ∈ {0.013, 0.026, 0.1, 0.3, 1, 3, 10, 30}. For all values of γ the calculations were found to converge with a value of − log(a) = −1.5 and a value of n = 128. The KC-RPMD TST rates are obtained by numerical integration of F (y). The transmission coefficients are calculated using KC-RPMD trajectories that are released from the dividing surface associated with y ‡ = 0. The true value of the mass of the auxiliary variable, my = 2.68 × 102 , is used for the dynamical KC-RPMD trajectories and for the calculation of γy , such that γy = 6.55 × 10−4 for all values of γ. For each value of the friction coefficient γ, a total of 8000 trajectories are released. Each trajectory is evolved for 200 fs using a time step of dt = 0.005 fs with initial velocities sampled from the Maxwell-Boltzmann distribution. The initial configurations for the KC-RPMD trajectories are generated from long KC-RPMD trajectories that are constrained to the dividing surface using the RATTLE algorithm; 24 the constrained trajectories are 1.6 ns in time and are thermostatted by resampling the velocities from the Maxwell-Boltzmann distribution every 100 fs. 110 3.4.4 KC-RPMD rate calculation in System B In System B, the frequency of the donor-acceptor mode, ωq , that modulates the strength of the  coupling, K(q), is varied such that ωq ∈ 4 × 10−4 , 7 × 10−4 , 1.4 × 10−3 , 2.2 × 10−3 , 3 × 10−3 . For these values of ωq , the calculations were found to converge with values of − log(a) ∈ {−1, −3.3, −4.3, −4.3, −4.3}, respectively. The ring polymer is discretized using n = 128 beads for all values of ωq . The KC-RPMD TST rates are obtained by numerical integration of F (y). The transmission coefficients are calculated using KC-RPMD trajectories that are released from the dividing surface associated with y ‡ = 0. Due to the small values of my presented in Table 3.2, we utilize a value of my = 10.0 for the KC-RPMD trajectories and the calculation of γy for all values of ωq , such  that γy ∈ 1.30 × 10−2 , 1.36 × 10−2 , 1.40 × 10−2 , 1.39 × 10−2 , 1.39 × 10−2 . For each value of the frequency ωq , a total of 1000 trajectories are released. Each trajectory is evolved for 100 fs using a time step of dt = 0.001 fs for ωq = 4 × 10−4 and dt = 0.002 for all other values of ωq ; initial velocities were sampled from the Maxwell-Boltzmann distribution. The initial configurations for the KC-RPMD trajectories are generated from long KC-RPMD trajectories that are constrained to the dividing surface using the RATTLE algorithm; 24 the constrained trajectories are 200 ps in time and are thermostatted by resampling the velocities from the Maxwell-Boltzmann distribution every 100 fs. 3.4.5 Calculation of reference TST expressions Reference values for the thermal reaction rates for System A are evaluated using the Marcus rate expression for non-adiabatic ET and the Zusman extension of Marcus theory to the overdamped regime of ET, the Zusman regime. The Marcus rate expression for non-adiabatic ET with classical solvent is given by 25–28 2π kMT = |K|2 ~ r   β (λ + )2 exp −β , 4πλ 4λ (3.39) where λ is the reorganization energy calculated as λ = V1 (s1 ) − V1 (s0 ), and all other terms have been defined previously. The Zusman expression of ET, which extends Marcus theory to account for frictional effects of the solvent, is given by 29 kZUS = kMT , 4πK 2 τL 1+ ~λ (3.40) 111 where τL = γ/(ωs2 ms ) is the Debye longitudinal relaxation time. 29,30 Reference values for the thermal reaction rates for System B are evaluated using an extension of the Marcus rate expression to account for a position dependent electronic coupling such that, 8,25 Z kDA = dqkMT (q)P (q), (3.41) where kMT (q) = 2π |K(q)|2 ~ r   β (λ + )2 exp −β 4πλ 4λ (3.42) is the Marcus theory expression at a given value of the donor-acceptor distance q, and P (q) = R e−βV0 (q) dqe−βV0 (q) (3.43) is the probability density that at a given q. The rate expression Eq. (3.41) is evaluated via numerical integration for each value of ωq . 3.5 Results The results in this paper are presented in three parts. In the first part we validate the updated form of KC-RPMD by investigating the transition between normal and inverted ET. In the second part we explore the performance of KC-RPMD in calculating the ET rate as a function of the friction coefficient of the bath. In the final part we present results simulating ET as a function of the frequency of the donor-acceptor mode, which modulates the magnitude of the electronic coupling. 3.5.1 Normal to inverted electron transfer We begin by considering numerical results for System A1, which models the transition from the normal to inverted regime of condensed-phase ET. Figure 3.1 depicts the thermal reaction rates calculated using KC-RPMD (red) and classical Marcus theory using Eq. (3.39) (black); Marcus theory constitutes the appropriate reference result since System A1 is in the weak-coupling regime (βK ≈ 7 × 10−4 ) and the solvent coordinate is treated classically. Comparison of the KC-RPMD results and Marcus theory demonstrates that the updated form of KC-RPMD exhibits quantitative agreement with the reference result throughout the normal, activationless, and inverted regime, analogous to what was seen in the original derivation of KC-RPMD. 1 Through introduction of the kinetic-constraint, the KC-RPMD method is able to correct the breakdown of the tunneling rate 112 through asymmetric barriers for instanton based methods that leads to an absence of the turnover in the ET rate for conventional position-based RPMD. 31,32 -10 Log(kET) -15 -20 -25 0 0.05 0.1 0.15 0.2 0.25 -ε Figure 3.1: Thermal reaction rate coefficients for System A1 as a function of the driving force, , obtained using KC-RPMD (red) and classical Marcus Theory (Eq. (3.39), black). 3.5.2 Friction controlled electron transfer We next present results for System A2 and System A3, which model ET as a function of the bath friction. We consider both the transition from the Marcus regime to the Zusman regime of ET and the transition from the weakly dissipative to Marcus regime of ET. The calculation of ET rates as a function of the friction of the bath provides a stringent test of the KC-RPMD methodology and its ability to treat systems in which dynamical effects play a strong role. 5.0e-19 kET 4.5e-19 4.0e-19 3.5e-19 -0.5 0 0.5 1 1.5 2 2.5 Log(γ/msωs) Figure 3.2: Thermal reaction rate coefficients for System A2 as a function of the friction of the bath, γ, obtained using KC-RPMD (red) and the Zusman expression (Eq. (3.40), black). 113 0.14 0.12 kET/kCL 0.1 0.08 0.06 0.04 0.02 0 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 Log(η/msωs) Figure 3.3: Thermal reaction rate coefficients for System A23 as a function of the friction of the bath, γ, obtained using KC-RPMD (red), the Zusman expression (Eq. (3.40), black), and numerically exact quantum dynamics calculations (blue). 20 Figure 3.2 presents the thermal reactions rates for System A2 calculated using KC-RPMD (red) and the Zusman expressing using Eq. (3.40) as a function of the friction of the bath. Both the Zusman expression and the rates calculated using KC-RPMD exhibit a plateau region for values of log(γ/ms ωx ) < 1.5, which corresponds to the Marcus regime of ET; in this regime, the reference rate expression reduces to the usual Marcus expression, Eq. (3.39), and is independent of the strength of the friction coefficient. In the overdamped regime, when γ > 1.5, both the Zusman expression and the rates calculated KC-RPMD exhibit a decrease in the ET rate associated with the diffusive dynamics of the solvent coordinate. 29,33–36 Comparison of the KC-RPMD results and the Zusman expression illustrates that KC-RPMD is able to capture the underlying physics of both regimes of ET. Specifically, the KC-RPMD methodology is able to capture the competition between the timescale of transition between electronic states and the time-scale for motion of the nuclei, illustrating the ability of KC-RPMD to simulate systems in which the coupled electronic and nuclear dynamics determines the reaction rate. To compare with calculations using exact quantum dynamics methods, Fig. 3.3 presents results as a function of the friction of the bath in terms of the dynamical recrossing factor defined as kET /kcl , where kET is the ET rate calculated using either KC-RPMD or exact quantum dynamics, and kcl is a classical TST rate, 20 kcl = ωs −βV0 (s‡ ) e . 2π (3.44) The solvent position s‡ corresponds to the value of the solvent coordinate for which the diabatic states are degenerate. The dynamical recrossing factor provides a measure of the amount of recrossing of the 114 dividing surface in the true system in comparison to the classical TST estimate. Figure 3.3 presents the thermal reaction rates for System A3 calculated using KC-RPMD (red), the Zusman expression for ET (black), and numerically exact quantum dynamics methods (blue) as a function of the friction of the bath. 20 For intermediate values of the friction coefficient, log(γ/ms ωs ) > −1, all three sets of ET rates are observed to be independent of the strength of the friction coefficient corresponding to the Marcus regime of ET. However, the magnitude of the dynamical recrossing factor in the plateau region is larger for the exact quantum results due to an increase in the ET rate associated with nuclear tunneling. 20 For the weakly dissipative regime when log(γ/ms ωs ) < −1, both the rates calculated using KC-RPMD and the exact quantum dynamics method, in comparison to the Zusman expression, show strong friction dependence, though with distinctly different qualitative behavior. In the weakly dissipative regime, the time-scale for relaxation of the solvent-coordinate is long, such that the solvent-coordinate is able to pass through configurations associated with the crossing of the diabats numerous times; these multiple crossing events allow for multiple transitions between electronic states to occur, a phenomenon not captured by the Zusman expression for ET. When the solvent is treated classically, as in the KC-RPMD simulations, the multiple electronic transitions lead to a decorrelation between the initial and final electronic state causing a decrease in the ET rate; as observed in Fig. 3.3 the KC-RPMD simulations are able to properly capture the physics governing this regime of ET. When the solvent is treated quantum mechanically, as in the exact quantum dynamics simulations, 20 the nuclear wavepackets are able to constructively interfere during each pass of the crossing of the diabatic states, leading to an enhancement of the ET rate. 20,34 Like all dynamics methods based on the imaginary-time formulation of statistical mechanics, 3 KC-RPMD is unable to capture such coherence effects due to the loss of phase-information in the approximate dynamics. However, Fig. 3.3 provides a promising result, illustrating that the KC-RPMD method is able to appropriately capture the competition between the time-scale for solvent relaxation and the probability of transitioning between electronic states. 3.5.3 Electron transfer with position dependent coupling We conclude by considering numerical results for System B, which models ET in the presence of a donor-acceptor mode that modulates the magnitude of the diabatic coupling. Figure 3.4 presents the thermal reactions rates calculated using KC-RPMD (red) and the extension of classical Marcus theory to include a donor-acceptor mode (Eq. (3.41), black). The rates calculated using KC-RPMD are in quantitative agreement with those calculated using the reference result for systems with 115 Log(kET) -21 -22 -23 0 0.001 0.002 0.003 ωq Figure 3.4: Thermal reaction rate coefficients for System B as a function of the friction of the friction of the donor-acceptor mode, ωq , obtained using KC-RPMD (red) and the extension of Marcus theory to include a donor-acceptor mode (Eq. (3.41), black). varying strength of the donor-acceptor mode. These results illustrate the ability of the KC-RPMD methodology to treat systems in which the diabatic coupling can depend strongly on the position of the nuclei. As such, KC-RPMD is well suited to investigate complex systems in which the distance between the ET donor and acceptor is a dynamical quantity. 3.6 Conclusions In this work, we have demonstrated the ability of the KC-RPMD methodology to properly simulate systems in which the coupled nuclear and electronic dynamics has a strong impact in determining the reaction rate and mechanism. In doing so, we revisit several key aspects of the KC-RPMD derivation, including the form of the kinetic constraint and the choice of mass associated with the continuous electronic variable. In particular, we introduce a Langevin bath that is coupled to the continuous electronic variable to correct for the presence of spurious resonances in the electronic variable, analogous to the treatment of spurious resonances associated with the internal modes of the ring-polymer in conventional RPMD. 5 The accuracy of the KC-RPMD methodology has been verified in a variety of model systems of condensed-phase electron transfer in which nuclear dynamics plays an important role. We have illustrated that KC-RPMD is able to capture the underlying physics that governs the transition from the weakly-dissipative to the Marcus to the Zusman regime of ET as a function of the friction of the bath. Specifcally, KC-RPMD properly treats the competing time-scales between the relaxation of 116 the solvent and the time-scale for electronic transition. Furthermore, the KC-RPMD methodology is able to correctly simulate systems in which the diabatic coupling can depend strongly on the position of the nuclei. Taken together, the presented work illustrates the ability of KC-RPMD to treat complex systems in which dynamical effects govern the reaction rate and mechanism. 117 Appendix A Derivation of the penalty function In this appendix, we derive the specific form of the multiplicative factor, η, that appears in the penalty function, g, in Eq. (3.13). The multiplicative factor is derived such that the KC-RPMD distribution does not bias the probability of kink-pair formation at nuclear configurations for which the diabatic states are degenerate. Specifically, we consider the relative probability between kinked configurations and configurations in the reactant basin, and equate this to the relative probability between kinked configurations at nuclear configurations within a small distance of the crossing of the diabatic surfaces and configurations in the reactant basin in the unmodified distribution. The unmodified distribution, ρn ({R(α) }, y), is defined analogously to Eq. (3.12), but where the penalty function g({i(α) }, {R(α) }) = 1 for all {i(α) }. For simplicity we first present the derivation for a 1D redox system with constant coupling, K, in the classical limit of the nuclear coordinates. We then outline the derivation for a general multi-dimensional system. 3.A.1 1D redox system with constant K and classical nuclei For a 1D system with classical nuclei, the kinetically constrained ring-polymer distribution (Eq. 3.12) has the form ρKC n (x, y) = Ω X (α) g({i(α) }, x)e−βVr (y,{i }) Γ({i(α) }, x), (3.45) {iα } (α) where Γ({i }, x) = n Y Mi(α) ,i(α+1)(x), and the penalty function in this case takes the form α=1 g({i(α) }, x) =       1, i(α) = 0 for all α, 1, i(α) = 1 for all α,      Ce−a(w(x))2 , otherwise. (3.46) In the kinetically constrained distribution the relative probability between kinked configurations and configurations in the reactant basin is given as P KC (y ‡ ) P KC (y < y ‡ ) ρKC (y ‡ )dy = R y‡n , KC (y) dyρ n −∞ (3.47) 118 where the incremental length dy is necessary to define the unitless probability, y ‡ = 0, and ρKC n (y) = Ω Z dxe−βVeff (x,y) . (3.48) The relative probability in Eq. (3.47) can be simplified by taking the limit of large a and by noting that the denominator is dominated by the statistical weight of unkinked configurations, such that 1 ‡ P KC (y = y ‡ ) Cdy e−βV0 (x ) =R KC ‡ P (y < y ) dx e−βV0 (x) L|w0 (x‡ )| r n/2 2k π X (βK) , a φn (k) (3.49) k=1 where x‡ denotes the point of the intersection of the diabatic surfaces, the prime denotes differenti −1  n 2 . ation with respect to the nuclear coordinate, and φn (k) = 2k 2k n We relate the relative probability in Eq. (3.49) to the relative probability between kinked configurations at nuclear configurations within a small distance of the crossing of the diabatic surfaces and configurations in the reactant basin in the unmodified distribution R x+ dy x− dxρn (x, y ‡ ) P (x− < x < x+ , y ‡ ) = R y‡ , R P (y < y ‡ ) dy dxρn (x, y) (3.50) −∞ where the distance x+ − x− defines a small region around the crossing of the diabatic surfaces. We define this distance as the Landau-Zener length, lLZ , which gives an estimate of the range of nuclear configurations for which the coupling, K, impacts the dynamics, 7 such that x− = x‡ − lLZ /2 and x+ = x‡ + lLZ /2, and lLZ = 2π|K| |F0 (x‡ ) − F1 (x‡ )| . (3.51) Noting again that the denominator is dominated by the statistical weight of unkinked configurations, and making the reasonable approximation that the diabatic surfaces do not change significantly over the Landau-Zener length, such that the integral can be approximated as Z x+ dxρn (x, y ‡ ) = lLZ Z dxδ(x − x‡ )ρn (x, y ‡ ), (3.52) x− Eq. (3.50) simplifies to ‡ n/2 2k P (x− < x < x+ , y ‡ ) e−βV0 (x ) dylLZ X (βK) =R . ‡ P (y < y ) φn (k) dx e−βV0 (x) L (3.53) k=1 To obtain the final expression for multiplicative factor in a 1D redox system with constant K 119 and classical nuclei, we equate the two relative probabilities in Eq. (3.49) and Eq. (3.53), yielding r C= 3.A.2 a 2π. π (3.54) Multi-dimensional redox system with position-dependent K(R) For a general multi-dimensional redox system with position-dependent K(R) and classical nuclei, the multiplicative factor η can be analogously derived. The relative probability in the constrained distribution is now given by P KC (y = y ‡ ) P KC (y < y ‡ ) dyC = dyC = ≈ dyC pπ R a pπ R a pπ R a dRδ(w(R))e−βV0 (R) R L dRe−βV0 (R) Pn/2 (βK)2k k=1 φn (k) dRδ(V0 (R) − V1 (R))|K(R)|e−βV0 (R) R L dRe−βV0 (R) (3.55) Pn/2 (βK)2k k=1 φn (k) dRδ(V0 (R) − V1 (R))β 2 |K(R)|3 e−βV0 (R) R , L dRe−βV0 (R) (3.56) (3.57) where in the last line we assume that terms associated with more than one kink-pair (k = 1) can be neglected. The relative probability in the unconstrained distribution can be defined as dylLZ P (R− < R < R+ , y ‡ ) = ‡ P (y < y ) R dRδ(V0 (R) − V1 (R))|F0 (R) − F1 (R)|ρn (R, y ‡ ) , R y‡ R dy dRρn (R, y) −∞ (3.58) where in multiple dimensions we define the Landau-Zener length as lLZ = 2πh|K(R)|ic , h|F0 (R) − F1 (R)|ic (3.59) and again we have assumed that the diabatic surfaces do not change over the Landau-Zener length. Simplifying Eq. 3.58, again assuming that terms associated with more than one kink-pair (k = 1) can be neglected, and equating the two relative probability yields the final form for the multiplicative constant in multiple dimensions with position dependent coupling r C= a h|K(R)|ic h|F0 (R) − F1 (R)| |K(R)|ic 2π . π h|F0 (R) − F1 (R)|ic h|K(R)|3 ic (3.60) Eq. (3.60) simplifies to (3.54) in the limit of 1D and a position dependent electronic coupling and is appropriately unitless. For the case of a multidimensional system with quantized nuclei, the 120 resulting expression for the multiplicative factor is unchanged if we make the approximation that the ring-polymer position is approximated by its centroid. 1 Appendix B Derivation of the Langevin friction coefficient In this appendix we derive the specific form for the Langevin friction coefficient associated with the auxiliary function that appears in the KC-RPMD equations of motion, Eq. (3.20). The Langevin friction coefficient is chosen such that the time-scale for thermalization of the auxiliary variable matches the time for the auxiliary variable to cross the reactant (or analogously the product) basin in the absence of the Langevin bath. The time-scale for relaxation, τr , is given by the exponential decay constant associated with the decay of the initial kinetic energy of the auxiliary variable when coupled to the Langevin bath, such that τt = 1/(2γy ). 37 The time for the auxiliary variable to cross the reactant basin the absence of the Langevin bath is given by τy = 2−L , |vyr | (3.61) where 2 − L defines the length of the reactant basin associated with the square restraining potential in Eq. (3.10), and vyr defines the average velocity of the auxiliary variable in the reactant basin following initialization of the auxiliary variable at the transition region y = y ‡ with an average initial velocity given by the Maxwell-Boltzmann distribution, 1 1 my (vyr )2 = + hVeff (R, y ‡ )ic − hVeff (R, y = 0)ic 2 2β (3.62) In Eq. 3.62 we have assumed that the motion of the nuclei is negligible on the time-scale of the motion of the auxiliary variable. Simplifying Eq. 3.62 yields the time for the auxiliary variable to cross the reactant basin s τy = (2 − L) 1 − 2 ln βmy   2−L a . 1/2 ηβ 2 − 4hln |K|i ( ) c L π (3.63) Equating the two time-scales τt and τy yields the final form for the Langevin friction coefficient 1 γy = 2(2 − L) v h u u 1 − 2 ln 2−L t L
a 1/2 ηβ 2 π βmy i − 4hln |K(R)|ic . (3.64) 121 References [1] A. R. Menzeleev, F. Bell, and T. F. Miller, III, J. Chem. Phys. 140, 064103 (2014). [2] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 (2004). [3] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, Annu. Rev. Phys. Chem. 64, 387 (2013). [4] R. P. Feynman and A. R. Hibbs, Quantum mechanics and path integrals, International series in pure and applied physics (McGraw-Hill, New York,, 1965) p. xiv. [5] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, J. Chem. Phys. 140, 234116 (2014). [6] D. Chandler, Introduction to modern statistical mechanics (Oxford university press, 1987). [7] V. May and O. Kühn, Charge and energy transfer dynamics in molecular systems (Wiley-VCH, 2011). [8] A. Nitzan, Chemical dynamics in condensed phases (Oxford university press, 2006). [9] D. Frenkel and B. Smit, Understanding molecular simulation: From algorithms to applications (Academic Press, San Diego, 2002). [10] E. Wigner, Phys. Chem. Abt. B 19, 203 (1932). [11] H. Eyring, J. Chem. Phys. 3, 107 (1935). [12] J. C. Keck, J. Chem. Phys. 32, 1035 (1960). [13] D. Chandler, J. Chem. Phys. 68, 2959 (1978). [14] C. H. Bennet, in Algorithms for Chemical Computations, edited by R. E. Christofferson (American Chemical Society, Washington, DC, 1977) p. 63. [15] E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989). [16] G. K. Schenter, B. C. Garrett, and D. G. Truhlar, J. Chem. Phys. 119, 5828 (2003). 122 [17] J. B. Watney, A. V. Soudackov, K. F. Wong, and S. Hammes-Schiffer, Chem. Phys. Lett. 25, 268 (2006). [18] M. Topaler and N. Makri, J. Chem. Phys. 101, 7500 (1994). [19] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 122, 084106 (2005). [20] M. Topaler and N. Makri, J. Phys. Chem. 100, 4430 (1996). [21] H. B. Gray and J. R. Winkler, Proc. Natl. Acad. Sci. 102, 3534 (2005). [22] L. Verlet, Phys. Rev. 159, 98 (1967). [23] E. Vanden-Eijnden and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006). [24] H. C. Andersen, J. Comput. Phys. 52, 24 (1983). [25] R. A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 265 (1985). [26] R. A. Marcus, J. Chem. Phys. 43, 3477 (1965). [27] R. A. Marcus, Disc. Faraday Soc. 29, 21 (1960). [28] R. A. Marcus, J. Chem. Phys. 24, 966 (1956). [29] I. Rips and J. Jortner, J. Chem. Phys. 87, 2090 (1987). [30] A. Garg, J. N. Onuchic, and V. Ambegoakoar, J. Chem. Phys. 83, 4491 (1985). [31] A. R. Menzeleev, N. Ananth, and T. F. Miller III, J. Chem. Phys. 135, 074106 (2011). [32] P. Shushkov, J. Chem. Phys. 138, 224102 (2013). [33] J. T. Hynes, J. Phys. Chem. 90, 3701 (1986). [34] J. N. Onuchic and P. G. Wolynes, J. Phys. Chem. 92, 6495 (1988). [35] D. L. Zusman, Chem. Phys. 49, 295 (1980). [36] M. Morillo and R. I. Cukier, J. Chem. Phys. 91, 281 (1989). [37] M. Tuckerman, Statistical mechanics: Theory and molecular simulation (Oxford University Press, 2010).