Capturing Nuclear Quantum Effects at Classical Efficiency: a Path-Integral Approach Citation Tao, Xuecheng (2022) Capturing Nuclear Quantum Effects at Classical Efficiency: a Path-Integral Approach. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/s929-5x12. https://resolver.caltech.edu/CaltechTHESIS:03282022-213805038 Abstract Quantum mechanical effects of nuclei are ubiquitous in chemistry. For a typical example, zero-point energies and tunneling effects of the nuclei shift the chemical equilibrium and manipulate the reaction rate. However, theoretical investigation of such nuclear quantum effects in chemical reactions remains a challenge due to the heavy computation cost. To this end, imaginary-time path-integral based approximate methods have been previously introduced, which allows the inclusion of nuclear quantization in real-time chemical dynamics simulations at the efficiency of classical Newtonian dynamics. In the dissertation, we further extend the applicability of those path-integral methods and exploit the methods for practical chemical investigations. Specifically, we introduce novel dynamics approaches based on ring-polymer molecular dynamics methodology to incorporate nuclear quantum effects in the simulations of excited state dynamics and microcanonical scattering processes, and to examine the nuclear quantum effects in Hydrogen/Deuterium sticking to the graphene surface. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: nuclear quantum effects, ring-polymer molecular dynamics, chemical dynamics, rate theory, isotope effects 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: Chan, Garnet K. (chair) Okumura, Mitchio Wei, Lu Miller, Thomas F. Defense Date: 22 March 2022 Record Number: CaltechTHESIS:03282022-213805038 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:03282022-213805038 DOI: 10.7907/s929-5x12 Related URLs: URL URL Type Description https://doi.org/10.1063/1.5005544 DOI Article adapted for Chapter 2 https://doi.org/10.1021/acs.jpca.9b00877 DOI Article adapted for Chapter 3 https://doi.org/10.1063/1.5144307 DOI Article adapted for Chapter 4 https://doi.org/10.1021/acs.jpclett.0c02933 DOI Article adapted for Chapter 5 https://doi.org/10.1063/5.0087156 DOI Article related to the thesis ORCID: Author ORCID Tao, Xuecheng 0000-0003-2907-3839 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 14531 Collection: CaltechTHESIS Deposited By: Xuecheng Tao Deposited On: 20 Apr 2022 19:42 Last Modified: 13 Jul 2022 22:54 Thesis Files PDF - Updated Version See Usage Policy. 4MB Repository Staff Only: item control page

Capturing Nuclear Quantum Effects at Classical Efficiency: A Path-Integral Approach Thesis by Xuecheng Tao In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2022 Defended March 22, 2022 ii © 2022 Xuecheng Tao ORCID: 0000-0003-2907-3839 All rights reserved except where otherwise noted iii To my family and my friends iv ACKNOWLEDGEMENTS I would like to begin by expressing my gratitude to my advisor, Professor Tom Miller. Tom has a natural flair for leadership. We team around him, and are kept inspired and influenced by him. Tom also makes an effort to create an open and inclusive environment for the group that encourages communication, to recognize people’s talents and their efforts, and to support the growth and career development of his students. As a consequence, my Ph.D. experience with Tom completely reshaped me, and allowed me to explore my role as a member of the scientific community and the even larger society. With Tom starting a new chapter of his life, I sincerely wish you, Tom, all the best for the adventure. Next, I would like to express my thanks to Caltech. Being a closely connected small community, Caltech is considerate, supportive, and makes every effort to take care of its students. I would like to thank the professors at Caltech, especially to my thesis committee—Professors Garnet Chan, Mitchio Okumura, and Lu Wei, for constantly taking the initiative to provide insightful scientific and career advice. A big thank you also goes to the administrative staff of the Division, in particular to Priscilla Boon and Alison Ross, for taking care of the tedious but critical chores for our group. Members from the Miller Group have always been the greatest support to my journey at Caltech, one thing that I cannot appreciate more. We together create a collaborative and warm group atmosphere and take care of each other during hard times. People come and leave, but the atmosphere remains. Thank you to you all for the guidance and help offered, for the discussions, debates, and watercooler chats—Ralph, Brett, Philip, Feizhi, Welborn, Leanne, Ioan, Saleh, Varun, Xinglong, Tamara, Steve, Tomislav, Emiliano, Josh, Connie, Mike, Michiel, Sebastian, Brooke; Zimmer, Jeongmin, Jorge and Dan; Sherry, Roman, Marta, James, Zhuoran, Jiace, Vignesh, Doris, Rui and Bo. I would like to express thank you to everyone for including me in the community and thank you for accepting me as being stubborn and socially awkward. Special thanks to the dynamics team folks, Philip, Roman, and Tomislav—without you, I would have kept struggling in the mud alone. Finally, it is my family that stays with me throughout the whole journey of my graduate study. I would like to thank my mother and grandmother in China for shouldering the family responsibilities for me when I was thousands of miles away v from home to pursue an international education. On the other side of the ocean, it is hard to get familiar with one other and build a lifelong partner relationship, not to mention in the tough isolation period—Linqing, I am blessed that we are gradually learning how to support each other and developing the courage to face unforeseen challenges in the future. There are just so many things to appreciate and I feel very fortunate that I had the best Ph.D. experience with the Miller Group at Caltech for the past years. I sincerely hope that everyone has the chance to work, study, and live a life in the way that one truly enjoys. Miller Group Snack Table, proudly hosted by Xuecheng Tao, 2018–2022. vi ABSTRACT Quantum mechanical effects of nuclei are ubiquitous in chemistry. For a typical example, zero-point energies and tunneling effects of the nuclei shift the chemical equilibrium and manipulate the reaction rate. However, theoretical investigation of such nuclear quantum effects in chemical reactions remains a challenge due to the heavy computation cost. To this end, imaginary-time path-integral based approximate methods have been previously introduced, which allows the inclusion of nuclear quantization in real-time chemical dynamics simulations at the efficiency of classical Newtonian dynamics. In the dissertation, we further extend the applicability of those path-integral methods and exploit the methods for practical chemical investigations. Specifically, we introduce novel dynamics approaches based on ringpolymer molecular dynamics methodology to incorporate nuclear quantum effects in the simulations of excited state dynamics and microcanonical scattering processes, and to examine the nuclear quantum effects in Hydrogen/Deuterium sticking to the graphene surface. vii PUBLISHED CONTENT AND CONTRIBUTIONS 1. Begušić, T., Tao, X. & Miller III, T. F. Equilibrium-Nonequilibrium RingPolymer Molecular Dynamics for Nonlinear Spectroscopy. arXiv preprint arXiv:2202.00887 (2022). Xuecheng contributed to the interpretation of the results in the project and to the writing of the final manuscript. 2. Jiang*, H., Tao*, X., Kammler, M., Ding, F., Wodtke, A. M., Kandratsenka, A., Miller III, T. F. & Bünermann, O. Small Nuclear Quantum Effects in Scattering of H and D from Graphene. Journal of Physical Chemistry Letters 12, 1991–1996. doi:10.1021/acs.jpclett.0c02933 (2021). Xuecheng contributed to the conception of the project, to theory developments and numerical computations in the project, and to the writing of the manuscript. 3. Tao, X., Shushkov, P. & Miller III, T. F. Microcanonical Rates from RingPolymer Molecular Dynamics: Direct-Shooting, Stationary-Phase, and MaximumEntropy Approaches. Journal of Chemical Physics 152, 124117. doi:10 . 1063/1.5144307 (2020). Xuecheng contributed to the conception of the project, to theory developments and numerical computations in the project, and to the writing of the manuscript. 4. Tao, X., Shushkov, P. & Miller III, T. F. Simple Flux-Side Formulation of State-Resolved Thermal Reaction Rates for Ring-Polymer Surface Hopping. Journal of Physical Chemistry A 123, 3013–3020. doi:10.1021/acs.jpca. 9b00877 (2019). Xuecheng contributed to the conception of the project, to theory developments and numerical computations in the project, and to the writing of the manuscript. 5. Tao, X., Shushkov, P. & Miller III, T. F. Path-Integral Isomorphic Hamiltonian for Including Nuclear Quantum Effects in Non-Adiabatic Dynamics. Journal of Chemical Physics 148, 102327. doi:10.1063/1.5005544 (2018). Xuecheng contributed to the conception of the project, to theory developments and numerical computations, and to the writing of the manuscript. viii TABLE OF CONTENTS Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Published Content and Contributions . . . . . . . . . . . . . . . . . . . . . . Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter I: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter II: Path-Integral Isomorphic Hamiltonian for Including Nuclear Quantum Effects in Non-Adiabatic Dynamics . . . . . . . . . . . . . . . . . . 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter III: Simple Flux-Side Formulation of State-Resolved Thermal Reaction Rates for Ring-Polymer Surface Hopping . . . . . . . . . . . . . . . 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Computational Details . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter IV: Microcanonical Rates from Ring-Polymer Molecular Dynamics: Direct-Shooting, Stationary-Phase, and Maximum-Entropy Approaches . 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Computational Details . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter V: Small Nuclear Quantum Effects in Scattering of H and D from Graphene . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Computational Methods . . . . . . . . . . . . . . . . . . . . . . . . iv vi vii vii x xii 1 4 4 4 6 12 17 25 26 32 32 32 33 40 44 48 49 49 50 51 56 57 68 69 71 71 71 73 73 ix 5.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Supplementary Materials . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 77 80 90 x LIST OF ILLUSTRATIONS Number Page 2.1 The isomorphic potentials and the QCLE correlation functions for a model two-level harmonic oscillators system. . . . . . . . . . . . . . 18 2.2 The isomorphic potential for a three-level donor-bridge-acceptor system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Matrix elements of the physical potential for the two-level reactive scattering system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.4 State-to-state thermal reaction rates as a function of temperature. . . . 23 2.5 Time-evolution of the initial phase-space distribution under QCLE dynamics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1 Schematic illustration of a hopping ring-polymer near avoided crossing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 The potential and the quantum transition state theory rates of a gasphase scattering model. . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Arrhenius plots of state-resolved thermal rates versus reciprocal temperature for two reaction channels. . . . . . . . . . . . . . . . . . . . 45 3.4 Dividing surface dependence of thermal rates at 𝛽 = 8. . . . . . . . . 46 4.1 An illustrative example of the “L-curve” that is used to determine the parameter 𝛼 in each MaxEnt calculation. . . . . . . . . . . . . . . . 57 4.2 Microcanonical reaction rate predictions for the Eckart barrier model from the direct shooting methods. . . . . . . . . . . . . . . . . . . . 59 4.3 Stationary phase approximation (SPA) results for the Eckart barrier. . 61 4.4 Maximum entropy (MaxEnt) inversion results for the Eckart barrier. . 63 4.5 Revisiting the microcanonical rate predictions for the Eckart barrier from the direct shooting methods. . . . . . . . . . . . . . . . . . . . 64 4.6 State-resolved microcanonical rates for the two-level model in Eq. 4.30. 66 5.1 Comparing theory with experiment for H and D scattering distributions from graphene. . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.2 Comparison of experiment and theory for H/D sticking probabilities to graphene. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3 Nuclear quantum effects in the H/D-on-graphene sticking. . . . . . . 79 xi 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 The minimum energy pathway (MEP) for the perpedicular approach of an H atom on top of a carbon atom on the graphene flake. . . . . . Isotope differences in the H/D-on-graphene scattering events predicted with theories at different levels. . . . . . . . . . . . . . . . . . Comparing theory with experiment—scattering distribution of H/D collision with graphene. . . . . . . . . . . . . . . . . . . . . . . . . Experimentally derived total survival fluxes. . . . . . . . . . . . . . Phonon density of states spectrum calculated with GFN-xTB and EMFT-REBO potential at 300 K, respectively. . . . . . . . . . . . . Sticking probability of H/D scattering on graphene as a function of normal incidence energy when using the EMFT-REBO potential. . . Interaction time of the scattered H/D atom with the graphene surface from simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . H atom sticking probability as a function of normal incidence energy— testing sensitivity with respect to the initialization of the H/D internal ring-polymer temperature. . . . . . . . . . . . . . . . . . . . . . . . 82 83 84 86 87 88 89 89 xii LIST OF TABLES Number Page 2.1 Parameter values for the physical potential of the two-level reactive scattering system, given in Eq. 2.32. . . . . . . . . . . . . . . . . . . 22 3.1 Parameter values for the potential of two-level reactive scattering model in Eq. 3.32. . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1 Parameters for the two-level model in Eq. 4.30. . . . . . . . . . . . . 67 5.1 Results for the reactive energy term that describes C-H interaction, as defined in Eqs. 5.4 and 5.5. . . . . . . . . . . . . . . . . . . . . . 81 1 Chapter 1 INTRODUCTION Quantum mechanical effects of nuclei are ubiquitous in chemistry. Liquid Helium does not freeze even at the temperature of absolute zero (at the standard atmospheric pressure) [1]; atomically thin membrane material that permits only proton penetration is designed and fabricated [2]; long-range charge transfer in the bacterial proteins cannot occur fast without the help from the fluctuating hydrogen bonds around the reaction center [3]; and even the sense of smell may rely on the quantum mechanical vibrations of odourant molecules [4, 5]; just to name but a few. From the simplest elementary substance to the materials, and to the life processes, it is clear to observe that nuclear quantum effects are not merely a concept at the microscopic scale, but can lead to significant consequences in the macroscopic world. However, theoretical investigation of such nuclear quantum effects in chemical reactions remains a challenge nowadays due to the heavy computation cost. The problem has long been recognized as ‘a curse of dimensionality’ in quantum mechanical simulations, and the most powerful computer soon feel exhausted when increasing the number of particles included in the simulations, if one chooses to stick with a fully exact quantum representation of the atoms and electrons. To this end, imaginary-time path-integral-based approximate methods have been previously introduced. In the approach, the original quantum mechanical degrees of freedom are mapped to an isomorphic classical system where the exact quantum statistical distribution is faithfully preserved. As a result, one is able to describe the effects of nuclear quantization in real-time chemical dynamics by simulating the classical Newtonian dynamics of the path-integral isomorphic system, that is to say, at a classical mechanical computational efficiency, and with support from the whole mature field of classical dynamics simulations. As will be presented, the topics of my graduate research close surround the central idea, to capture the nuclear quantum effects at classical efficiency. In the dissertation, one will see how the applicability of imaginary-time path-integral methods grows as newly developed members join the family; one will also see how the 2 imaginary-time path-integral methods advance the exploration of nuclear quantum effects in various beautiful chemistry contexts. The remaining dissertation is organized as follows. Chapter 2 presents a path-integral approach for including nuclear quantum effects in electronically non-adiabatic chemical dynamics simulations. We derive an isomorphic Hamiltonian such that for any given physical system with multiple electronic states, the Boltzmann sampling of the isomorphic Hamiltonian with classical nuclear degrees of freedom yields the exact quantum Boltzmann distribution for the original physical system. The isomorphic Hamiltonian can be straightforwardly combined with existing mixed quantum-classical dynamics methods, and the numerical results demonstrate the use of the method for the inclusion of nuclear quantum effects, when in combination with either fewest-switches surface hopping or the quantum-classical Liouville equation. The development leads to a variety of promising, novel dynamics methods that retain the simplicity and robustness of both imaginary-time path-integrals for nuclear quantization and the parent mixed quantum-classical method. Chapter 3 presents the further extension of the aforementioned isomorphic Hamiltonian framework to perform state-resolved thermal reaction rate calculations where nuclear quantum effects play an important role. We describe an efficient fluxside formulation of ring-polymer surface hopping in the isomorphic Hamiltonian framework, which is able to calculate multiple state-resolved non-adiabatic thermal reaction rates with only a single free-energy surface calculation, whereas previous non-adiabatic flux-side formulations for surface hopping involve multiple free-energy surface calculations. The method is shown to be robust and straightforwardly implemented, and numerical results reveal that the method leads to better transition-state dividing-surface independence, due to improved preservation of the path-integral statistics. The combined accuracy and simplicity of the method make it amenable to useful application in realistic chemical systems. Chapter 4 presents the approaches to calculate microcanonical reaction rates for processes involving significant nuclear quantum effects using ring-polymer molecular dynamics, both with and without electronically non-adiabatic transitions. After illustrating the shortcoming of the naive direct-shooting method, for which the results are strongly sensitive to the internal ring-polymer temperature that is employed, we investigate alternative strategies based on the expression for the microcanonical rate in terms of the inverse Laplace transform of the thermal reaction rate. It 3 is found that accurate microcanonical rates are obtained from inversions using either the stationary phase approximation or the numerically exact maximum entropy method. We also point out that the use of an internal temperature based on the stationary-phase approximation partially took account of the tunneling effects in the low incidence energy region. The work suggests general strategies for the extraction of microcanonical dynamical quantities from ring-polymer molecular dynamics or other approximate thermal simulations. Chapter 5 presents the study of nuclear quantum effects in H/D sticking to the graphene surface, comparing scattering experiments at near-zero coverage with classical simulation, transition-state calculations, and quantized simulation with ring-polymer molecular dynamics. Experiment shows H/D sticking probabilities that are indistinguishable from one another and markedly smaller than those expected from a consideration of zero-point energy shifts of the chemisorption transition state. The inclusion of dynamical effects and vibrational anharmonicity via ring-polymer molecular dynamics yields results that are in good agreement with the experiment. RPMD also reveals that nuclear quantum effects, while modest, arises primarily from carbon and not from H/D motion, confirming the importance of a C-atom re-hybridization mechanism associated with H/D sticking on graphene. 4 Chapter 2 PATH-INTEGRAL ISOMORPHIC HAMILTONIAN FOR INCLUDING NUCLEAR QUANTUM EFFECTS IN NON-ADIABATIC DYNAMICS This chapter is based on the following publication: 1. Tao, X., Shushkov, P. & Miller III, T. F. Path-Integral Isomorphic Hamiltonian for Including Nuclear Quantum Effects in Non-Adiabatic Dynamics. Journal of Chemical Physics 148, 102327. doi:10.1063/1.5005544 (2018). 2.1 Abstract We describe a path-integral approach for including nuclear quantum effects in nonadiabatic chemical dynamics simulations. For a general physical system with multiple electronic energy levels, a corresponding isomorphic Hamiltonian is introduced, such that Boltzmann sampling of the isomorphic Hamiltonian with classical nuclear degrees of freedom yields the exact quantum Boltzmann distribution for the original physical system. In the limit of a single electronic energy level, the isomorphic Hamiltonian reduces to the familiar cases of either ring polymer molecular dynamics (RPMD) or centroid molecular dynamics Hamiltonians, depending on implementation. An advantage of the isomorphic Hamiltonian is that it can easily be combined with existing mixed quantum-classical dynamics methods, such as surface hopping or Ehrenfest dynamics, to enable the simulation of electronically non-adiabatic processes with nuclear quantum effects. We present numerical applications of the isomorphic Hamiltonian to model two- and three-level systems, with encouraging results that include improvement upon a previously reported combination of RPMD with surface hopping in the deep-tunneling regime. 2.2 Introduction Chemical processes that involve transitions among different electronic states play a central role in photo-induced [6, 7], redox [8, 9], and collisional processes [10, 11]. Widely used mixed quantum-classical (MQC) methods—including Ehrenfest dynamics [12] and surface hopping [13]—have been developed for the simulation of electronically non-adiabatic processes in cases for which the nuclei can be described using classical mechanics. However, nuclear quantum effects are important in many 5 electronically non-adiabatic processes [14–17], creating the need for new methods that robustly and accurately describe the interplay between nuclear and electronic quantum mechanical effects. For chemical dynamics on a single electronic surface, approximate methods based on imaginary-time Feynman path integrals [18, 19] have proven useful for describing nuclear quantization. These methods include ring-polymer molecular dynamics (RPMD) [20, 21] and centroid molecular dynamics (CMD) [22–24], which involve classical molecular dynamics trajectories governed by an isomorphic Hamiltonian that includes the effects of zero-point energy and tunneling. RPMD and CMD exhibit various exact formal properties, including time-reversibility and preservation of the quantum Boltzmann distribution for the physical system, and RPMD additionally recovers semiclassical instanton rate theory in the deep-tunneling regime [25]. The simplicity and robustness of these path-integral-based methods has led to the development of mature technologies [26–29] and enables the study of complex systems [30–33]. These successes motivate the development of path-integral-based methods for describing electronically non-adiabatic dynamics. Previous work includes non-adiabatic extensions of instanton theory [34–36], CMD [37, 38], and RPMD [39–45]. A unifying feature of these previous efforts is that they employ a case-specific development strategy, in which path-integral quantization of the nuclei is specifically tailored for combination with a particular approximation to the electronically non-adiabatic dynamics, such as instanton theory [34–36, 42, 43], surface-hopping [41, 46], linearized semiclassical [47–52], or other approximation. This strategy typically limits each resulting method to the application domain for which the associated non-adiabatic dynamics approximation is valid. The current work employs an alternative strategy to take full advantage of the diversity of previously developed MQC methods for describing non-adiabatic dynamics. We use path integration to obtain a general isomorphic Hamiltonian that incorporates nuclear quantization and that can be easily combined with any MQC method. As will be shown, this leads to a variety of promising, new dynamics methods that retain the simplicity and robustness of both imaginary-time path-integrals for nuclear quantization and the parent MQC method. In the following, we derive the new isomorphic Hamiltonian, and we present applications of it in combination with non-adiabatic dynamics based on either surface hopping [13] or the quantum- 6 classical Liouville equation [53, 54]. These results illustrate the flexibility with which the isomorphic Hamiltonian may be employed, as well as implementations that are readily applicable for the study of complex systems. 2.3 Theory We begin by reviewing the path-integral-based RPMD and CMD methods, which employ an isomorphic Hamiltonian for the description of quantized nuclear dynamics in electronically adiabatic systems. We then extend this approach to obtain an isomorphic Hamiltonian for the description of quantized nuclear dynamics involving multiple electronic states. Isomorphic Hamiltonian for one-level systems: RPMD and CMD For a system obeying the Born-Oppenheimer approximation in the electronic ground state, we consider the Hamiltonian operator 2 𝑝 + 𝑉 (𝑥), 𝐻ˆ = 2𝑚 (2.1) where 𝑥, 𝑝, and 𝑚 are the nuclear position, momentum, and mass, respectively, and 𝑉 (𝑥) is the potential energy surface. Throughout this work, results will be presented for a single nuclear degree of freedom; generalization to multiple dimensions is straightforward. The path-integral discretization of the quantum mechanical canonical partition function for this system is given by [18, 19, 55] 𝑄 ˆ = tr[𝑒 −𝛽 𝐻 ] ∫  𝑛  𝑛∫ iso = lim 𝑑x 𝑑p 𝑒 −𝛽𝐻𝑛 (x,p) , 𝑛→∞ 2𝜋ℏ (2.2) where 𝛽 is the reciprocal temperature, 𝑛 is the number of ring-polymer beads in the path-integral discretization, x = {𝑥 (1) , 𝑥 (2) , . . . , 𝑥 (𝑛) } is the vector of ring-polymer positions such that 𝑥 (1) = 𝑥 (𝑛+1) , and p is the vector of ring-polymer momenta. 𝐻𝑛iso is the ring-polymer Hamiltonian (see Appendix A) 𝐻𝑛iso (x, p) = 𝑛 ∑︁ 𝑝2 𝑛 1 ∑︁ + 𝑈spr (x) + 𝑉 (𝑥 𝛼 ), 2𝑚 𝑛 𝑛 𝛼=1 𝛼=1 𝛼 (2.3) which includes the inter-bead potential 𝑈spr (x) = 𝑛 ∑︁ 1 𝑚 𝑛 𝜔2𝑛 (𝑥 𝛼 − 𝑥 (𝛼+1) ) 2 , 2 𝛼=1 (2.4) 7 where 𝑚 𝑛 = 𝑚/𝑛, 𝜔𝑛 = (𝛽𝑛 ℏ) −1 , and 𝛽𝑛 = 𝛽/𝑛. Approximate real-time quantum dynamics is obtained in the RPMD method [20] by running classical molecular dynamics trajectories associated with the ring-polymer Hamiltonian, which are given by 𝑥¤𝛼 = 𝑝 𝛼 /𝑚 𝑛 (2.5)
1 𝜕 𝑝¤𝛼 = 𝑚 𝑛 𝜔2𝑛 𝑥 (𝛼+1) + 𝑥 (𝛼−1) − 2𝑥 𝛼 − 𝑉 (𝑥 𝛼 ) 𝑛 𝜕𝑥 𝛼 for 𝛼 = 1, . . . , 𝑛. Equation 2.2 can be further reduced with respect to the intra-ring-polymer degrees of freedom, yielding ∫  𝑛  𝑛∫ ¯ iso ¯ 𝑄= 𝑑 𝑥¯ 𝑑 𝑝¯ 𝑒 −𝛽 𝐻 (𝑥,¯ 𝑝) , (2.6) 2𝜋ℏ where 𝐻¯ iso is the centroid Hamiltonian 2 𝑝¯ 𝐻¯ iso ( 𝑥, ¯ 𝑝) ¯ = + 𝑉¯ ( 𝑥) ¯ 2𝑚 which includes the centroid potential of mean force ∫ ∫ iso 1 ∑︁ −𝛽𝑉¯ ( 𝑥) ¯ 𝑒 ∝ lim 𝑑x 𝑑p 𝛿( 𝑥¯ − 𝑥 𝛼 )𝑒 −𝛽𝐻𝑛 (x,p) . 𝑛→∞ 𝑛 𝛼 (2.7) (2.8) Approximate real-time quantum dynamics is obtained in the CMD method [24] by running classical molecular dynamics trajectories associated with the centroid Hamiltonian, which are given by 𝑥¤̄ = 𝑝/𝑚 ¯ 𝜕 ¯ . 𝑝¤̄ = − 𝑉¯ ( 𝑥) 𝜕 𝑥¯ (2.9) Both Eqs. 2.3 and 2.7 provide an isomorphic Hamiltonian for the one-level physical system described by Eq. 2.1, in the sense that classical mechanical trajectories associated with the isomorphic Hamiltonian yield the approximate quantum mechanical time-evolution for the physical system. Moreover, classical Boltzmann sampling of the isomorphic Hamiltonian (i.e., by running the classical trajectories in Eqs. 2.5 or 2.9 in contact with a thermal bath) rigorously preserves the exact quantum Boltzmann statistics associated with the physical system. In the following, we derive both RPMD and CMD versions of the corresponding isomorphic Hamiltonian for physical systems involving multiple electronic surfaces, with the RPMD version presented in the main text and the CMD version in Appendix B. 8 Isomorphic Hamiltonian for multi-level systems Path-integral discretization Consider the Hamiltonian in the diabatic representation for a system with 𝑓 electronic energy levels, 𝑝2 ˆ + 𝑉 (𝑥) 2𝑚  𝑉1 (𝑥) 𝐾12 (𝑥)   2  𝐾12 (𝑥) 𝑉2 (𝑥) 𝑝 +  . = .. 2𝑚  .. .   𝐾1 𝑓 (𝑥) 𝐾2 𝑓 (𝑥)  𝐻ˆ = (2.10) · · · 𝐾1 𝑓 (𝑥)   · · · 𝐾2 𝑓 (𝑥)  ..  . ... .   · · · 𝑉 𝑓 (𝑥)  Discretizing the partition function with respect to both electronic state, 𝑖, and nuclear position, 𝑥, and employing a Trotter factorization such as ˆ ˆ ˆ ˆ 𝑒 −𝛽𝑛 𝐻 = 𝑒 −𝛽𝑛𝑉/2 𝑒 −𝛽𝑛𝑇 𝑒 −𝛽𝑛𝑉/2 + O (𝛽𝑛3 ), we obtain the path-integral representation  𝑛 𝑛 𝑄 = lim 𝑛→∞ 2𝜋ℏ ∫ ∫ Í𝑛 𝑝 2𝛼 × 𝑑x 𝑑p 𝑒 −𝛽( 𝛼=1 2𝑚𝑛 +𝑈spr (x)) 𝜇(x), where " 𝜇(x) = tre 𝑛 Ö (2.11) (2.12) # 𝑒 −𝛽𝑛𝑉ˆ (𝑥 𝛼 ) . (2.13) 𝛼=1 The subscript ‘e’ in Eq. 2.13 indicates the trace taken over only the electronic states. Although path-integral discretization of multi-level systems can also be performed in the adiabatic representation [56], the diabatic representation employed here is particularly convenient. Note that 𝜇, which describes the statistical weight of a given ring-polymer nuclear configuration after thermally averaging over the electronic states, is a familiar and easily evaluated quantity. It is the central object in the Schwieters-Voth non-adiabatic instanton theory [35–37] and mean-field non-adiabatic RPMD [42, 57, 58], both of which provide a thermally averaged (i.e., mean-field) description of the electronically non-adiabatic dynamics. Moreover, as is discussed in Appendix C, 𝜇 is non-negative when evaluated in the limit of large 𝑛, and both 𝜇 and its derivative with respect to the ring-polymer nuclear coordinates can be evaluated using O (𝑛) operations. 9 The Isomorphic Hamiltonian We now address the central goal of this work: Given the physical system associated with the 𝑓 -level Hamiltonian in Eq. 2.10, determine the corresponding 𝑓 -level isomorphic Hamiltonian for which classical Boltzmann sampling of the nuclear degrees of freedom yields the exact quantum Boltzmann distribution for the physical Hamiltonian. It follows from Eq. 2.12 that this requirement is satisfied by an isomorphic Hamiltonian of the form 𝐻ˆ 𝑛iso (x, p) = 𝑛 ∑︁ 𝑝2 𝛼 2𝑚 𝑛 𝛼=1 + 𝑈spr (x) + 𝑉ˆ iso (x), (2.14) where 𝑉ˆ iso is the isomorphic potential energy given by the 𝑓 × 𝑓 matrix that obeys i h ˆ iso (2.15) tre 𝑒 −𝛽𝑉 (x) ≡ 𝜇(x). Special case of a two-level system For a system with two electronic states ( 𝑓 = 2), the isomorphic potential energy has the form " # iso (x) 𝐾 iso (x) 𝑉 1 12 𝑉ˆ iso (x) = . (2.16) iso 𝐾12 (x) 𝑉2iso (x) Given the symmetry of the off-diagonal term, the matrix has only three independent elements at any given ring-polymer configuration. To specify the two diagonal terms, we require that the usual RPMD surfaces be recovered in the regime of zero electronic coupling, such that 𝑛 𝑉𝑖iso (x) = 1 ∑︁ 𝑉𝑖 (𝑥 𝛼 ). 𝑛 𝛼=1 (2.17) The only remaining term is the off-diagonal isomorphic coupling, 𝐾𝑖iso 𝑗 (x), which must satisfy Eq. 2.15, such that  𝛽    2 iso (x)+𝑉 iso (x) 𝑉 2 iso 𝑗 𝐾𝑖 𝑗 (x) = acosh 𝑒 2 𝑖 𝜇𝑖 𝑗 (x)/2 /𝛽2  2 iso iso − 𝑉𝑖 (x) − 𝑉 𝑗 (x) /4, where " 𝑛 Ö 𝜇𝑖 𝑗 (x) = tre 𝛼=1 " exp −𝛽𝑛 𝑉𝑖 (𝑥 𝛼 ) 𝐾𝑖 𝑗 (𝑥 𝛼 ) 𝐾𝑖 𝑗 (𝑥 𝛼 ) 𝑉 𝑗 (𝑥 𝛼 ) (2.18) # !# . (2.19) 10 For the case of a two-level system, 𝜇𝑖 𝑗 (x) = 𝜇(x), where the latter is defined in Eq. 2.13. Eq. 2.18 fully specifies 𝐾𝑖iso 𝑗 (x) to within an absolute sign, which we take to be equal to that of the physical potential coupling evaluated at the ring-polymer centroid position, sgn(𝐾𝑖 𝑗 ( 𝑥)). ¯ For a two-level system, the isomorphic Hamiltonian is given by Eq. 2.14 and Eqs. 2.16-2.18. Inspection of the matrix elements of the isomorphic potential reveals that the diagonal matrix elements (Eq. 2.17) include RPMD-like corrections to the diabatic potential energy surfaces, while the off-diagonal elements (Eq. 2.18) include the effect of nuclear quantization on the pairwise (i.e., two-body) coupling between the electronic states. Before discussing other properties of the isomorphic Hamiltonian, we generalize it to multi-level systems. General case of a multi-level system Following the two-level case, we now present the generalization of the isomorphic Hamiltonian to systems with 𝑓 > 2. We define an 𝑓 × 𝑓 potential energy matrix  𝑉 iso (x) 𝐾 iso (x) · · · 𝐾 iso (x)   1  12 1𝑓  iso (x) 𝑉 iso (x) · · · 𝐾 iso (x)   𝐾 12 2 2𝑓 iso  𝑉ˆ2-body (x) =  .. .. ..  ... .  . .     𝐾 iso (x) 𝐾 iso (x) · · · 𝑉 iso (x)    1𝑓 2𝑓 𝑓 (2.20) for which the diagonal and off-diagonal terms are defined in Eqs. 2.17 and 2.18. And finally, to ensure that Eq. 2.15 is satisfied, we define the isomorphic potential energy to be iso iso 𝑉ˆ iso (x) = 𝑉ˆ2-body (x) + 𝑉many-body (x), (2.21) where     𝜇(x) 1   iso i 𝑉many-body (x) = − ln  h iso ˆ 𝛽  tr 𝑒 −𝛽𝑉2-body (x)   e    and 𝜇(x) is defined in Eq. 2.13. (2.22) Combined with Eq. 2.14, Eqs. 2.20-2.22 present the central result of this work: the isomorphic Hamiltonian for a general multi-level system. We now point out a number of important properties that make the isomorphic Hamiltonian amenable to the description of complex, multi-level systems, much like standard RPMD and CMD are amenable to the description of complex, one-level systems. 11 First, the isomorphic Hamiltonian can immediately be employed with any MQC method for describing nonadiabatic dynamics; by simply running the MQC dynamics on the isomorphic Hamiltonian, nuclear quantum effects are included via the path-integral description. Naturally, the dynamics run on the isomorphic Hamiltonian will inherit the strengths and weaknesses of the MQC method that is employed. As is illustrated in the Results section, the MQC dynamics can either be run directly using the diabatic representation or by diagonalizing it to obtain the corresponding adiabatic states and derivative couplings. Second, by construction, the isomorphic Hamiltonian satisfies the requirement that classical Boltzmann sampling of the nuclear degrees of freedom yields the exact quantum Boltzmann distribution for the physical system. It employs a path-integral discretization that involves no approximation to the quantum statistics of the system. For an (idealized) MQC method for which the equations of motion rigorously preserve the MQC Boltzmann ensemble, then running the corresponding dynamics on the isomorphic Hamiltonian would rigorously preserve the exact quantum Boltzmann distribution; however, we note that most MQC methods do not rigorously preserve the MQC Boltzmann ensemble [59]. Third, as for standard RPMD, evaluation of the matrix elements in the isomorphic Hamiltonian is numerically robust and scales linearly in cost with the number of ring-polymer beads. Quantities that arise the evaluation of the isomorphic h in i iso −𝛽𝑉ˆ2-body (x) Hamiltonian, such as 𝜇(x), 𝜇𝑖 𝑗 (x), or tre 𝑒 (and their derivatives with respect to nuclear position), can be obtained from simple diagonalization of an 𝑓 × 𝑓 matrix or with O (𝑛) operations. Furthermore, the argument of the logartithm in Eq. 2.22 involves a ratio of positive quantities and is thus well behaved. It should be noted that the numerical robustness of the isomorphic Hamiltonian is an important and non-trivial feature; whereas evaluation of the path-integral representation for the underlying density matrix of a many-level system generally gives rise to a numerical sign problem [60], we have expressed the isomorphic Hamiltonian in terms of non-oscillatory quantities. We further note that the isomorphic Hamiltonian obeys various satisfying limits. In the classical mechanical limit for the physical nuclei (i.e., the 1-bead ring polymer limit), the isomorphic Hamiltonian reduces to the original physical Hamiltonian in Eq. 2.10. In the limit of zero coupling among the states in the physical system (i.e., when 𝐾𝑖 𝑗 = 0), the isomorphic Hamiltonian reduces to the standard RPMD Hamiltonian for the diabatic potential energy surfaces. Finally, in the limit for 12 iso (x) = 0, the which the electronic states only couple via separate pairs, 𝑉many-body many-level isomorphic Hamilton simply reduces to the previously discussed twoiso (x) includes the effect of nuclear quantization level result. In this sense, 𝑉ˆ2-body on the pairwise (i.e., two-body) coupling between the electronic states, whereas iso (x) provides a mean-field many-body coupling between the electronic 𝑉many-body states due to nuclear quantization. As will be seen in the results, this many-body coupling is found to be much smaller than the two-body coupling, but inclusion of the many-body term is necessary to rigorously preserve the quantum Boltzmann statistics. Finally, we note that the specification of the matrix elements of the isomorphic potential presented here is not unique. For example, direct inversion of the electronic density matrix within the trace operation of Eq. 2.13 was explored and found to be numerically ill-conditioned. Other alternative choices that satisfy the condition in Eq. 2.15 may be devised, although any revision should both preserve the formal properties listed above and improve upon the numerical results presented in the Results section. We do recognize that a representation-invariant specification of the matrix elements of the isomorphic potential would be a worthy goal for future development. Similarly, we recognize the mathematical possibility that the RHS of Eq. 2.18 may become negative in our specification (although we have found no such case in which this occurs), and we note that the positivity of 𝜇 guarantees the existence of a specification for which the matrix elements of the isomorphic potential are everywhere real. 2.4 Applications The isomorphic Hamiltonian can be used to incorporate nuclear quantum effects in any MQC simulation. To illustrate this, present applications in which the isomorphic Hamiltonian is combined with either quantum-classical Liouville equation (QCLE) or fewest-switches surface hopping non-adiabatic dynamics. Below, we briefly summarize the equations of motion associated with these two MQC methods. QCLE Dynamics The time evolution of a general operator in a multi-level system according to QCLE dynamics is given by [53, 54, 61] 𝜕 𝑂ˆ W (𝑥, 𝑝, 𝑡) = L̂ 𝑂ˆ W (𝑥, 𝑝, 𝑡), 𝜕𝑡 (2.23) 13 where L̂ = i 1 n o n o 𝑖 h Ĥ , • − Ĥ , • − •, Ĥ . ℏ 2 (2.24) In these equations, 𝑂ˆ W (𝑥, 𝑝, 𝑡) is an 𝑓 × 𝑓 matrix that corresponds to the partial Wigner distribution for a given operator with respect to a subset of the degrees of freedom [62, 63], and Ĥ is a generic Hamiltonian in the diabatic representation. Our motivation for using the QCLE approach is to obtain a MQC limit in which the electronic dynamics evolves quantum mechanically and the nuclear dynamics evolves classically. Taking the limit of small ℏ, the partial Wigner distribution ˆ such that the QCLE dynamics reduces to the MQC phase-space distribution 𝑂, retains the same form, except that 𝜕 𝑂ˆ (𝑥, 𝑝, 𝑡) = L̂ 𝑂ˆ (𝑥, 𝑝, 𝑡). 𝜕𝑡 (2.25) Eqs. 2.24 and 2.25 thus cleanly define a MQC limit, where the first term in the RHS of Eq. 2.24 describes the quantum evolution of the electronic states via the commutator, and the second term describes both the classical evolution of the nuclear coordinates and the back-reaction to the quantum subsystem via the symmetrized Poisson bracket. Having taken the classical limit for the nuclei, the Kubo-transformed positionautocorrelation function ∫ 𝛽 h i 1 ˆ ˆ 𝑐˜𝑥𝑥 (𝑡) = 𝑑𝜆 tr 𝑒 −(𝛽−𝜆) 𝐻 𝑥ˆ 𝑒 −𝜆 𝐻 𝑥(𝑡) ˆ (2.26) 𝛽𝑄 0 becomes ∫ h  i 𝑑𝑥 𝑑𝑝 tre 𝑥𝑒 −L̂𝑡 𝑒 −𝛽Ĥ 𝑥 , (2.27) 2𝜋ℏ where we have taken advantage of time-reversal symmetry to ensure that the timeevolved distribution in Eq. 2.27 is conveniently numerically evaluated. 𝑐˜𝑥𝑥 (𝑡) = In this study, we consider the correlation function in Eq. 2.27, with the nuclei clasˆ where 𝐻ˆ is sically evolved either with respect to the physical Hamiltonian (Ĥ = 𝐻, given in Eq. 2.10) or with respect to the CMD version of the isomorphic Hamiltonian (Ĥ = 𝐻ˆ ciso , where 𝐻ˆ ciso is given in Eq. 2.43). The resulting dynamics is used to study two- and three-level systems with a single nuclear degree of freedom. Specifically, we investigate a two-level system comprised of shifted quartic oscillators with constant potential coupling, as well as a three-level system comprised of shifted harmonic oscillators with constant potential coupling. 14 The equations of motion in Eqs. 2.24 and 2.25 are evolved exactly on a numerical grid, using the interaction picture with Heisenberg evolution applied to the quantum subsystem; the resulting time-evolution is both numerically stable and avoids additional approximations to the QCLE dynamics, such as the momentum-jump approximation [64]. The midpoint finite-difference method [65] is used to integrate the partial differential equations. We employ a numerical grid that spans the range of positions for which the classical Boltzmann probability density exceeds 10−12 , 257 grid points in both 𝑥 and 𝑝 directions, and an integration timestep of 2.5 × 10−4 a.u. The matrix elements of the isomorphic potential, 𝑉¯𝑖iso and 𝐾¯ 𝑖iso 𝑗 , are sampled to convergence using path-integral Monte Carlo with 16𝛽 ring-polymer beads. In the Results section, for comparison with the approximate QCLE dynamics described by Eqs. 2.24 and 2.25, we additionally obtain numerically exact quantum mechanical results by propagating the Schrodinger equation in the discrete variable representation (DVR) [66, 67] on a grid. As is necessary, we confirm that the DVR results are identical to the QCLE dynamics in the high-temperature limit. Additionally, for any temperature, we confirm that the DVR results are identical to the QCLE dynamics for the case of a two-level system comprised of linearly coupled harmonic oscillators when Ĥ = 𝐻ˆ and the dynamics is initialized from the multi-level partial Wigner phase-space distribution (Appendix 2.7) [53, 54]. Surface Hopping Dynamics Consider a generic 𝑓 -level system with 𝑑 nuclear degrees of freedom and diabatic Hamiltonian 𝑑 1 ∑︁ 𝑚 𝑗 𝑦¤ 𝑗 + V̂ (y), (2.28) Ĥ = 2 𝑗=1 where V (y) is the diabatic potential energy matrix that depends on the nuclear positions, y = {𝑦 1 , . . . , 𝑦 𝑑 }, and 𝑚 𝑗 is the mass of the 𝑗 th degree of freedom. In fewest-switches surface hopping [13], quantum evolution of the electronic wavefunction 𝜓(y, 𝑡) along a given trajectory obeys 𝑖ℏ 𝜕 𝜓(y, 𝑡) = V̂ (y)𝜓(y, 𝑡), 𝜕𝑡 (2.29) and classical evolution of the nuclear coordinates obeys 𝑚 𝑗 𝑦¥ 𝑗 = − 𝜕 E 𝑘 (y), 𝜕𝑦 𝑗 (2.30) 15 where E 𝑘 is the 𝑘 th adiabatic Born-Oppenheimer surface obtained by diagonalizing the diabatic potential matrix. The nuclear trajectory evolves along a particular BornOppenheimer surface, subject to stochastic hops to other surfaces with probability   2 𝑝 𝑘𝑙 = max − Re((𝑑𝑙 𝑘 · 𝑣)𝑎 𝑘𝑙 )Δ𝑡, 0 (2.31) 𝑎𝑘 𝑘 where 𝑎 𝑘𝑙 is the element of the electronic density matrix in the adiabatic representation, (𝑑𝑙 𝑘 · 𝑣) is the inner product of the first-derivative non-adiabatic coupling with the nuclear velocity vector, and Δ𝑡 is the integration timestep. During hopping events, the total energy associated with the Hamiltonian in Eq. 2.28 is conserved by modifying the component of the velocity along the non-adiabatic coupling vector that connects the two surfaces; hops are forbidden if there is insufficient velocity in this component to ensure energy conservation. We implement forbidden hops without momentum reversal [68, 69], and we neglect decoherence corrections [70, 71], although either could easily be implemented in the current context. In this study, we consider various implementations of fewest-switches surface hopping in a two-level gas-phase scattering system that is a function of a single nuclear coordinate: (i) For the standard case of surface-hopping with classical nuclei (hereafter referred to as SH-classical), we employ Eqs. 2.28-2.31 using the physical Hamiltonian (Ĥ = ˆ given in Eq. 2.10) which includes the physical diabatic potential matrix (V̂ = 𝑉) ˆ 𝐻, as a function of the single nuclear coordinate, such that y = 𝑥. (ii) To quantize the nuclei in the surface hopping dynamics with the CMD version of the isomorphic Hamiltonian (referred to as SH-C-iso), we employ Eqs. 2.28-2.31 using Ĥ = 𝐻ˆ ciso (given in Eq. 2.43), which includes the CMD version of the diabatic potential matrix (V̂ = 𝑉ˆciso , given in Eq. 2.46) as a function of the centroid nuclear coordinate, such that y = 𝑥. ¯ (iii) To quantize the nuclei in the surface hopping dynamics with the RPMD version of the isomorphic Hamiltonian (referred to as SH-RP-iso), we employ Eqs. 2.282.31 using Ĥ = 𝐻ˆ 𝑛iso (given in Eq. 2.14), which includes the RPMD version of the diabatic potential matrix (V̂ = 𝑈spr + 𝑉ˆ iso , given in Eqs. 2.4 and 2.21, respectively) as a function of the ring-polymer coordinates, such that y = x. (iv) Finally, for comparison with an earlier effort to combine RPMD with surface hopping, we also employ the method described in Reference [41] using the “centroidapproximation” defined therein; this method is referred to as SH-RP-nokinks, since 16 it neglects the contribution of the “kinked” ring-polymer configurations that span multiple diabatic surfaces, such that Eq. 2.15 is not obeyed and the quantum Boltzmann statistics are approximated. Note that for all surface-hopping calculations reported here, the dynamics is run in a representation for which the number of electronic states is the same as for the physical system. For results obtained using the various versions of the isomorphic Hamiltonian, the surface-hopping dynamics involves transitions between the adiabatic potential surfaces obtained by diagonalizing the isomorphic diabatic potential energy matrix. Following the implementation in Reference [41], Eq. 2.29 is evolved in the interaction representation using a fourth-order Runge-Kutta integrator [65], and Eq. 2.30 is evolved using the velocity Verlet algorithm [72]. As in previous RPMD simulations, each timestep for the nuclear degrees of freedom involves separate coordinate updates due to forces arising from the adiabatic potential and due to exact evolution of the purely harmonic portion [21, 30]. Matrix elements of the centroid isomorphic potential, 𝑉¯𝑖iso and 𝐾¯ 𝑖iso 𝑗 , are sampled to convergence using path-integral Monte Carlo with either 8𝛽 ring-polymer beads (for 𝛽 ≤ 9) or 24𝛽 ring-polymer beads (for 𝛽 > 9); the larger number of ring-polymer beads was found to be more important for improving statistical sampling of the centroid potential surfaces than for converging the path-integral discretization. The SH-RP-iso results were likewise performed using 8𝛽 ring-polymer beads. For all cases, Eq. 2.30 is integrated with a timestep of 10−4 a.u. Thermal rates in this study are calculated via Boltzmann averaging of the microcanonical reactive probabilities, initializing trajectories outside of the interaction region with a momentum range for which the ratio of the corresponding Boltzmann-weighted microcanonical reactive probability to the total thermal rate is greater than 10−8 a.u. For the SH-classical and SH-C-iso calculations, for which the microcanonical reactive probability changes abruptly at the threshold energy, we discretize this momentum interval at a resolution of 0.01 a.u; for the SH-RP-iso and SH-RP-nokinks calculations, we use a discretization of 0.05 a.u. The microcanonical reactive probabilities are calculated using from 104 to 105 trajectories. In the results section, for comparison with the various surface-hopping implementations, we additionally obtain numerically exact quantum mechanical results via wavepacket propagation, using the split-operator Fourier transform method of Feit and Fleck [73] extended to multiple potential energy surfaces. A wavepacket was 17 initialized in the asymptotic reactant region and evolved forward in time until the scattering event was completed. An absorbing potential was placed in the asymptotic reactant region that eliminated the reflected portion of the scattered wavepacket, while the transmitted component was projected out in the asymptotic product region. The scattering amplitudes were calculated by Fourier transform of the transmitted fraction of the wavepacket, properly normalized, and the squared modulus of the scattering amplitudes is numerically integrated to obtain the quantum rates. To illustrate the full details of our implementation of the SH-RP-iso method, we have provided an example program online.1 2.5 Results We now present numerical results for two possible combinations of the new pathintegral isomorphic Hamiltonian with MQC methods. First, to investigate a welldefined limit for MQC non-adiabatic dynamics in combination with the isomorphic Hamiltonian, we employ the QCLE method, considering both a two-level system of coupled quartic oscillators and a three-level system involving a donor-bridgeacceptor model. Then, to investigate a broadly applicable combination of MQC nonadiabatic dynamics with the isomorphic Hamiltonian, we employ fewest-switches surface hopping to study a model for state-resolved gas-phase reactive scattering. Unless otherwise specified, quantities are reported in atomic units, and we employ a nuclear mass of 𝑚 = 1. QCLE Dynamics Two-level system: Coupled quartic oscillators We begin by considering a two-level system involving a single nuclear coordinate, for which the physical potential energy matrix, 𝑉ˆ (𝑥), is comprised of diagonal elements that are strongly anharmonic quartic oscillators, 𝑉1 (𝑥) = (𝑥 + 𝑥 0 ) 4 /16 and 𝑉2 (𝑥) = (𝑥 − 𝑥 0 ) 4 /16, and the off-diagonal elements, 𝐾12 (𝑥) = Δ, are constant. The lateral shift of the potentials is 𝑥 0 = (32/𝛽) 1/4 , such that the activation energy associated with the crossing of the diabats is consistently 2/𝛽. In studying this system, we will consider (i) numerically exact quantum dynamics, (ii) the classical nuclear limit in which the QCLE dynamics is run using the physical Hamiltonian, 1 See https://github.com/thomasfmiller/SH-RP-iso for the available code. potential energy 0 1 -8 -4 β=8 B -4 β=1 -8 2 0 2 x 0 x 0 4 4 8 8 _ K iso (x) / ∆ _ iso K (x) / ∆ A 1.0 1.5 2.0 1.0 1.5 2.0 -4 -4 β∆ = 1 D -8 -8 β∆ = 0.1 C x 0 x 0 4 β=1 β=8 4 β=1 β=8 8 8 ~ cxx (t) ~ cxx (t) 1.2 1.8 2.4 3.0 3.5 5.0 6.5 8.0 0 0 F E 5 5 t t 10 15 15 β=8 β∆ = 0.1 10 QM QCLE-classical QCLE-C-iso β=1 β∆ = 0.1 0.0 1.2 2.4 3.6 -6.4 -1.2 4.0 9.2 0 0 H G 5 5 t t 10 β=8 β∆ = 1 10 β=1 β∆ = 1 15 15 Figure 2.1: The isomorphic potentials and the QCLE correlation functions for a model two-level harmonic oscillators system. (A,B) Diagonal potential energy matrix elements for the coupled quartic oscillator system at high temperature (𝛽 = 1, panel A) and low temperature (𝛽 = 8, panel B). Matrix elements for the physical potential 𝑉𝑖 (𝑥) and for the CMD version of the isomorphic potential 𝑉¯𝑖iso ( 𝑥) ¯ are shown in dashed and solid lines, respectively. Matrix elements for diabats 1 and 2 are shown in blue and red, respectively. iso , normalized by the off-diagonal coupling (C,D) The off-diagonal matrix element of the CMD version of the isomorphic potential, 𝐾¯ 12 Δ in the physical potential, for weak coupling (𝛽Δ = 0.1, panel C) and for intermediate coupling (𝛽Δ = 1, panel D). High-temperature (𝛽 = 1) and low-temperature (𝛽 = 8) results are shown in black and red, respectively. (E-H) Kubo-transformed position-autocorrelation functions obtained using exact quantum mechanics (QM; black, dots), QCLE dynamics with classical nuclei (QCLE-classical; blue, dashed), and QCLE dynamics with nuclei quantized via the CMD version of the isomorphic Hamiltonian (QCLE-C-iso; red, solid). Results are presented for weak coupling and high temperature (𝛽Δ = 0.1, 𝛽 = 1; panel E), weak coupling and low temperature (𝛽Δ = 0.1, 𝛽 = 8; panel F), intermediate coupling and high temperature (𝛽Δ = 1, 𝛽 = 1; panel G), and intermediate coupling and low temperature (𝛽Δ = 1, 𝛽 = 8; panel H). potential energy ~ cxx (t) ~ cxx (t) 4 18 19 𝐻ˆ (𝑥), and (iii) the case of quantized nuclei in which the QCLE dynamics is run using the CMD version of the isomorphic Hamiltonian, 𝐻ˆ ciso . Methodological and computational details are provided in Section 2.4. Figs. 2.1A-D illustrate the matrix elements of the CMD version of the isomorphic potential, 𝑉ˆciso (Eq. 2.46). In solid lines, panels A and B present the diagonal elements of the isomorphic potential, 𝑉¯1iso ( 𝑥) ¯ and 𝑉¯2iso ( 𝑥), ¯ at high and low temperature, with the physical diabatic potentials 𝑉1 (𝑥) and 𝑉2 (𝑥) shown in dashed lines for comparison. Given that these isomorphic potential matrix elements are identical to the CMD potentials of mean force for the two diabats, they exhibit the familiar features of converging to the physical potential at high temperature (Fig. 2.1A) and exhibiting larger nuclear quantization effects at low temperature (Fig. 2.1B). For weak coupling (𝛽Δ = 0.1) and intermediate coupling (𝛽Δ = 1), Figs. 2.1C and D, respectively, present the off-diagonal matrix elements of the isomorphic iso , at both low (red) and high (black) temperature. Unlike the coupling potential, 𝐾¯ 12 in the physical potential for this model, Δ, the coupling in the isomorphic potential is position dependent, reflecting the changing thermal probability of kinked ringpolymer configurations at different nuclear configurations. In all cases, the inclusion of nuclear quantization via exact path-integral statistics leads to an increase in the effective coupling between the two diabatic surfaces in the vicinity of the diabatic crossing (𝑥 = 0), with more pronounced effects at lower temperature. Figs. 2.1E and F present results for the Kubo-transformed position-autocorrelation function (Eq. 2.26) in the weak-coupling regime (𝛽Δ = 0.1) at high and low temperature, respectively. At the higher temperature (Fig. 2.1E), there is little difference in the QCLE dynamics obtained with classical nuclei (QCLE-classical; blue, dashed) versus with nuclei quantized via the CMD version of the isomorphic Hamiltonian (QCLE-C-iso; red, solid), and both implementations of QCLE are in good agreement with exact quantum mechanics (black, dots) due to the small role of nuclear quantum effects. At low temperatures, however, substantial nuclear quantum effects emerge, as evidenced by the difference between the blue and black curves in Fig. 2.1F. In this low-temperature case, the QCLE-C-iso dynamics exhibit substantial improvement, recovering the exact quantum result at 𝑡 = 0 as a necessary consequence of the path-integral statistics and showing better agreement with the quantum mechanical period of oscillation. 20 Finally, Figs. 2.1G and H present results for the Kubo-transformed position autocorrelation function in the intermediate-coupling regime (𝛽Δ = 1) at high and low temperature, respectively. As before, at high temperature (Fig. 2.1G), the QCLE-classical dynamics differs little from the QCLE-C-iso dynamics; however, both differ substantially from the exact quantum result at longer times. At low temperature (Fig. 2.1H), even larger differences are observed. As is necessary, QCLE-C-iso recovers the exact quantum result at short times, but it deviates from both QCLE-classical and exact quantum results at longer times. As is familiar from standard CMD and RPMD in one-level systems [20, 22], the results in Figs. 2.1E-H highlight that the newly introduced isomorphic Hamiltonian provides a means of exactly incorporating the statistical effects of nuclear quantization while only approximately including the dynamical effects. Moreover, the dynamics obtained from the isomorphic Hamiltonian will reflect the particular shortcomings of the employed MQC method—in this case, QCLE initialized with the MQC phase-space distribution. In Appendix D, we illustrate that a leading source of error for the QCLE-C-iso results in Figs. 2.1E-H is non-preservation of the MQC phase-space distribution in the QCLE dynamics at lower temperatures, where the MQC phase-space distribution differs substantially from the partial Wigner distribution. Three-level system: Donor-Bridge-Acceptor model For systems with more than two levels, a many-body correction appears in the iso isomorphic potential to ensure exact Boltzmann statistics (𝑉many-body in Eq. 2.21 and iso ¯ 𝑉many-body in Eq. 2.46). To investigate the nature of this many-body term, we consider a previously studied model for a three-level donor-bridge-acceptor system [36]. For this system, the physical potential energy, 𝑉ˆ (𝑥), is comprised of diagonal elements that are harmonic oscillators (𝑉1 (𝑥) = (𝑥 + 𝑥0 ) 2 /2, 𝑉2 (𝑥) = 𝑥 2 /2, and 𝑉3 (𝑥) = (𝑥 − 𝑥 0 ) 2 /2), and the off-diagonal elements are constant (𝐾12 (𝑥) = 𝐾23 (𝑥) = Δ, 𝐾13 (𝑥) = 0). The lateral shift of the potentials is 𝑥 0 = 4/𝛽1/2 , such that the activation energy associated with the crossing of the diabats is 2/𝛽. For the case of 𝛽 = 1, the diagonal elements of the physical potential are shown in Fig. 2.2A. Upon computing the matrix elements for the CMD version of the isomorphic potential, 𝑉ˆciso (Eq. 2.46), it is found that the diagonal (not shown) and off-diagonal (Fig. 2.2B) contributions to the two-body isomorphic potential (Eq. 2.47) are qualitatively similar to those illustrated in Figs. 2.1A-D. The many-body contribution to 21 potential energy 4 2 A 0 -8 -4 0 4 8 4 8 4 8 x _ iso K /∆ 1.4 B 1.2 1.0 -8 -4 0 iso log |Vmany-body /∆ | x 0 C -4 -8 -8 -4 0 x Figure 2.2: The isomorphic potential for a three-level donor-bridge-acceptor system. (A) Diagonal matrix elements of the physical potential for the three-level donor-bridge-acceptor system with 𝛽 = 1. (B) Off-diagonal matrix elements of the CMD version of the isomorphic potential, 𝐾¯ iso , normalized by Δ. (C) Manyiso body contribution to the isomorphic potential of the three-level system, 𝑉¯many-body , normalized by Δ. Results are presented for weak coupling and high temperature (𝛽Δ = 0.1, 𝛽 = 1; black, solid), weak coupling and low temperature (𝛽Δ = 0.1, 𝛽 = 8; red, solid), intermediate coupling and high temperature (𝛽Δ = 1, 𝛽 = 1; black, dashed), and intermediate coupling and low temperature (𝛽Δ = 1, 𝛽 = 8; red, dashed). In panel B, the high-temperature results (black lines) are graphically indistinguishable. 22 potential energy 8 4 0 Diabat 1 Diabat 2 Coupling x10 -4 -8 -8 -4 0 x 4 8 Figure 2.3: Matrix elements of the physical potential for the two-level reactive scattering system. The elements include diabat 1 (red), diabat 2 (black), and the off-diagonal coupling (blue, dashed, with 10-fold magnification). iso the isomorphic potential of the three-level system, 𝑉¯many-body , is plotted in Fig. 2.2C, divided by Δ to illustrate the magnitude of this many-body term in comparison to the two-body potential coupling. As is clear from the log-scale in Fig. 2.2C, we find in all studied cases that the many-body contribution is negligible in comparison to the two-body coupling between the electronic states. As a result, the dynamics for this system exhibits very little three-body character, and the computed time correlation functions (not shown) exhibit the qualitative features of those discussed in Figs. 2.1E-H. We thus find that the isomorphic Hamiltonian can be straightforwardly applied in multi-level systems and that, at least for the three-level system studied here, the many-body contribution to the isomorphic potential plays a minor role. Table 2.1: Parameter values for the physical potential of the two-level reactive scattering system, given in Eq. 2.32. Parameter Value 𝐴1 7 𝐴2 −18/𝜋 𝐴3 0.25 𝐵1 −0.75 𝐵2 54/𝜋 Parameter Value 𝑎1 √ 1 𝑎2 3𝜋/4 𝑎3 0.25 𝑥1 −1.6 𝑥3 −2.625 23 -2 4 -3 log k1 A 8 -4 0 3 3.5 4 QM SH-classical SH-RP-nokinks SH-C-iso SH-RP-iso -4 -8 0 1 2 3 4 105 / T -3 4 -4 log k2 B 8 -5 0 3 3.5 4 QM SH-classical SH-RP-nokinks SH-C-iso SH-RP-iso -4 -8 0 1 2 3 4 105 / T Figure 2.4: State-to-state thermal reaction rates as a function of temperature. The results are obtained using surface hopping with classical nuclei (SH-classical; black, dashed) and with nuclei quantized via the SH-RP-nokinks (green), SH-C-iso (blue), and SH-RP-iso (red) methods, as well as with exact quantum mechanics (black, dots). (A) The rate (𝑘 1 ) for the channel that enters on diabat 1 and exits on diabat 2. (B) The rate (𝑘 2 ) for the channel that enters on diabat 2 and exits on diabat 2. The insets expand the axes in the low-temperature region. Unless explicitly shown, the error bars are smaller than the size of the plotted circles. For the inset of panel (B), the SH-RP-iso and SH-RP-nokinks results are within the statistical error at all temperatures. Both temperature and the reaction rate are reported in SI units. 24 Surface-hopping dynamics We finally consider the state-to-state reactive scattering in a two-level model for a gas-phase system with a single nuclear degree of freedom. The physical potential for this system is given by matrix elements 𝑉1 (𝑥) = 𝐴1 + 𝐵1 −𝑎 1 + 𝑒 1 (𝑥−𝑥1 ) 𝑉2 (𝑥) = 𝐵2 𝐴2 +
1 + 𝑒 −𝑎2 𝑥 4 cosh2 𝑎22 𝑥 𝐾12 (𝑥) = 𝐴3 𝑒 −𝑎3 (𝑥−𝑥3 ) (2.32) 2 with parameters given in Table 2.1. Both the diagonal and off-diagonal potential matrix elements are plotted in Fig. 2.3, with reactants at 𝑥 → −∞ and products at 𝑥 → ∞. The basic features of this model resemble the F+H2 co-linear reaction, exhibiting both endothermal and exothermal reactive channels. We consider the thermal reaction rate 𝑘 1 for the channel that enters on diabatic state 1 and exits on diabatic state 2, as well as the thermal reaction rate 𝑘 2 for the channel that enters on diabatic state 2 and exits on diabatic state 2. The state-to-state thermal reaction rates are calculated using methods that include (i) numerically exact quantum dynamics, (ii) surface hopping with classical nuclei (SH-classical), (iii) surface hopping with nuclei quantized via the ring-polymer surface hopping method in Reference [41] that approximates the path-integral statistical distribution (SH-RP-nokinks), (iv) surface hopping with nuclei quantized via the CMD version of the isomorphic Hamiltonian (SH-C-iso), and (v) surface hopping with nuclei quantized via the RPMD version of the isomorphic Hamiltonian (SH-RP-iso). Both the SH-C-iso and SH-RP-iso methods are newly presented in this work. Results were also obtained using classical Ehrenfest dynamics [12], but are excluded due to their poor quality for this model. Computational details are provided in Section 2.4, and an example program that runs the SH-RP-iso trajectories for the system studied here is provided online.2 Figure 2.4A presents results for the thermal reaction rate 𝑘 1 obtained using the various methods as a function of reciprocal temperature, with the inset providing an expanded view of the lowest-temperature results. The large differences between the exact quantum and SH-classical results at low temperature illustrate the strong role of nuclear quantum effects. Although the SH-RP-nokinks method qualitatively recovers the effect of nuclear tunneling in this process, it overestimates the thermal reaction rate at low temperatures by at least an order of magnitude (see inset). Since 2 See https://github.com/thomasfmiller/SH-RP-iso for the available code. 25 SH-RP-nokinks neglects ring-polymer configurations that span the two electronic surfaces, it underestimates the role of the low-lying excited state in suppressing nuclear tunneling; similar errors are observed when standard RPMD on the lower adiabatic surface is used to approximate tunneling through an avoided crossing (see Fig. 2 of Reference [45]). It is clear that both the SH-C-iso and SH-RP-iso results in Fig. 2.4A are in better agreement with the exact quantum results, with the RPMD version of the isomorphic Hamiltonian leading to particularly accurate results. Figure 2.4B presents the corresponding results for the thermal reaction rate 𝑘 2 . Again, large nuclear quantum effects at low temperature are indicated by the difference between the exact quantum and SH-classical results. The inset reveals that for this reactive channel, the SH-C-iso method exhibits the largest errors among the quantized surface hopping methods, overestimating the reaction rate by an order of magnitude in the deep-tunneling regime (𝛽 > 𝛽c ≈ 8 for diabat 2). This result illustrates a well-known shortcoming of CMD for deep-tunneling across asymmetric barriers [74], which is the precise nature of the reaction channel associated with 𝑘 2 . For this process, the SH-RP-iso and SH-RP-nokinks are graphically indistinguishable and are in good agreement with the exact quantum results. We note that this simple model for a gas-phase scattering reaction reveals a significant shortcoming of both the SH-RP-nokinks and CMD-based methods for describing non-adiabatic chemical dynamics. Surface hopping combined with the RPMD version of the isomorphic Hamiltonian (SH-RP-iso) avoids these pitfalls and provides the best accuracy for both reactive channels at all temperatures. 2.6 Summary The current work strives to decouple the methodological challenge of describing electronically non-adiabatic dynamics from that of describing nuclear quantization. For a general physical system with multiple electronic energy levels, we derive a corresponding isomorphic Hamiltonian, such that Boltzmann sampling of the isomorphic Hamiltonian with classical nuclear degrees of freedom yields the exact quantum Boltzmann distribution for the original physical system. The key advantage of this isomorphic Hamiltonian is that it can be combined with existing mixed quantum-classical (MQC) methods for non-adiabatic dynamics, allowing for the straightforward inclusion of nuclear quantum effects. 26 The isomorphic Hamiltonian is presented in two versions, one of which recovers standard ring-polymer molecular dynamics (RPMD) in the limit of a single electronic surface, and the other that recovers standard centroid molecular dynamics (CMD). Numerical results are presented using both the RPMD and CMD versions of the isomorphic Hamiltonian, in combination with either fewest-switches surface hopping or the quantum-classical Liouville equation (QCLE) descriptions of MQC non-adiabatic dynamics. Investigation of a simple model for non-adiabatic gas-phase scattering reveals that a particularly promising approach is to combine surface-hopping dynamics with the RPMD version of the isomorphic Hamiltonian (i.e., the SH-RP-iso method), which exhibits the best accuracy among the studied methods for two different reactive channels at all temperatures. Future work will include applications of the isomorphic Hamiltonian to explore the role of nuclear quantum effects in the non-adiabatic dynamics of complex systems. Methodological extensions of the current work are also of interest, including alternative specification of the matrix elements of the isomorphic Hamiltonian (as discussed in Section 2.3), and combination of the isomorphic Hamiltonian with other MQC methods for describing non-adiabatic dynamics. Also of interest are dimensionality-reduction strategies based on generalization of the isomorphic potential energy in Eq. 2.21 to describe the correlated dynamics of a local subset of electronic states embedded in a mean-field treatment of the environment (akin to quantum embedding strategies for electronic structure [75]). 2.7 Appendix A. Equivalent forms of the ring-polymer Hamiltonian The ring-polymer Hamiltonian is usually introduced [20, 21] by writing the partition function as ∫ ∫ −𝑛 𝑑x 𝑑p 𝑒 −𝛽𝑛 𝐻𝑛 (x,p) , (2.33) 𝑄 = lim (2𝜋ℏ) 𝑛→∞ where 𝐻𝑛 = 𝑛 ∑︁ 𝑝2 𝛼 𝛼=1 2𝑚 + 𝑛𝑈spr (x) + 𝑛 ∑︁ 𝑉 (𝑥 𝛼 ) (2.34) 𝛼=1 and 𝑈spr (x) is defined in Eq. 2.4. The RPMD equations of motion associated with this form of the Hamiltonian are 𝑥¤𝛼 = 𝑝 𝛼 /𝑚 𝑝¤𝛼 (2.35) 𝜕 = 𝑚𝜔2𝑛 𝑥 (𝛼+1) + 𝑥 (𝛼−1) − 2𝑥 𝛼 − 𝑉 (𝑥 𝛼 ) 𝜕𝑥 𝛼
27 or
1 𝜕 𝑥¥𝛼 = 𝜔2𝑛 𝑥 (𝛼+1) + 𝑥 (𝛼−1) − 2𝑥 𝛼 − 𝑉 (𝑥 𝛼 ) 𝑚 𝜕𝑥 𝛼 for 𝛼 = 1, . . . , 𝑛, and the Lagrangian associated with this Hamiltonian is L= 𝑛 ∑︁ 1 2 𝛼=1 𝑚 𝑥¤𝛼2 − 𝑛𝑈spr (x) − 𝑛 ∑︁ 𝑉 (𝑥 𝛼 ). (2.36) (2.37) 𝛼=1 Now, we introduce a new Lagrangian that is obtained by constant scaling of the original, L iso ≡ L/𝑛, (2.38) which yields the corresponding Hamiltonian 𝐻𝑛iso = 𝑛 ∑︁ ( 𝑝 iso ) 2 𝛼 𝛼=1 2𝑚 𝑛 𝑛 1 ∑︁ 𝑉 (𝑥 𝛼 ). + 𝑈spr (x) + 𝑛 𝛼=1 (2.39) The classical equations of motion associated with this Hamiltonian are 𝑥¤𝛼 = 𝑝 iso 𝛼 /𝑚 𝑛 (2.40) = 𝑚 𝑛 𝜔2𝑛 𝑥 (𝛼+1) + 𝑥 (𝛼−1) − 2𝑥 𝛼 − 𝑝¤iso 𝛼
1 𝜕 𝑉 (𝑥 𝛼 ) 𝑛 𝜕𝑥 𝛼 or
1 𝜕 𝑉 (𝑥 𝛼 ) . (2.41) 𝑥¥𝛼 = 𝜔2𝑛 𝑥 (𝛼+1) + 𝑥 (𝛼−1) − 2𝑥 𝛼 − 𝑚 𝜕𝑥 𝛼 Comparison of Eqs. 2.36 and 2.41 confirms that since the two forms of the Hamiltonian (in Eqs. 2.34 and 2.39) are obtained from constant scaling of the same Lagrangian, they yield the same equations of motion. Finally, we can rewrite the exponand in Eq. 2.33 as " 𝑛 # 𝑛 1 ∑︁ 1 ∑︁ 𝑝 2𝛼 −𝛽𝑛 𝐻𝑛 = −𝛽 + 𝑈spr (x) + 𝑉 (𝑥 𝛼 ) 𝑛 𝛼=1 2𝑚 𝑛 𝛼=1 " 𝑛 # 𝑛 ∑︁ 1 ∑︁ 1 𝑉 (𝑥 𝛼 ) = −𝛽 𝑚 𝑛 𝑥¤𝛼2 + 𝑈spr (x) + 2 𝑛 𝛼=1 𝛼=1 " 𝑛 # 𝑛 ∑︁ ( 𝑝 iso ) 2 1 ∑︁ 𝛼 = −𝛽 + 𝑈spr (x) + 𝑉 (𝑥 𝛼 ) 2𝑚 𝑛 𝑛 𝛼=1 𝛼=1 = −𝛽𝐻𝑛iso . We have thus shown that the partition function in Eq. 2.33 can equivalently be rewritten as ∫  𝑛  𝑛∫ iso iso 𝑄 = lim 𝑑x 𝑑piso 𝑒 −𝛽𝐻𝑛 (x,p ) (2.42) 𝑛→∞ 2𝜋ℏ 28 and that the Hamiltonian in Eq. 2.39 yields the usual RPMD equations of motion. In the main text, we employ Eqs. 2.42 and 2.39 for the partition function and the ringpolymer Hamiltonian, respectively, and for succinctness, we drop the superscript “iso” in denoting the bead momenta. B. CMD version of the isomorphic Hamiltonian The CMD version of the isomorphic Hamiltonian is 2 𝑝¯ 𝐻ˆ ciso ( 𝑥, ¯ 𝑝) ¯ = + 𝑉ˆciso ( 𝑥), ¯ 2𝑚 (2.43) where 𝑉ˆciso is the isomorphic potential energy given by the 𝑓 × 𝑓 matrix that obeys h i ˆ iso ¯ tre 𝑒 −𝛽𝑉 (𝑥) ≡ 𝜇( ¯ 𝑥), ¯ (2.44) ∫ 𝜇( ¯ 𝑥) ¯ = lim 𝐶 𝑑x 𝛿( 𝑥¯ − 𝑛→∞ 1 ∑︁ 𝑥 𝛼 ) 𝑒 −𝛽𝑈spr (x) 𝜇(x), 𝑛 𝛼 (2.45) √  𝑚𝑛  (𝑛−1)/2 , and 𝜇(x) is given by Eq. 2.13. Following the logic of the 𝐶 = 𝑛 2𝜋𝛽ℏ 2 main text, we obtain the centroid isomorphic potential energy of the form iso iso ( 𝑥), ¯ ( 𝑥) ¯ + 𝑉¯many-body 𝑉ˆciso ( 𝑥) ¯ = 𝑉ˆ¯2-body (2.46) which includes the two-body contribution  𝑉¯ iso ( 𝑥) ¯ iso ¯  1 ¯ 𝐾12 ( 𝑥)  ¯ iso  𝐾12 ( 𝑥) ¯ 𝑉¯2iso ( 𝑥) ¯ iso ˆ  ¯ 𝑉2-body ( 𝑥) ¯ = .. .. . .    𝐾¯ iso ( 𝑥) 𝐾¯ 2iso𝑓 ( 𝑥) ¯  1𝑓 ¯ · · · 𝐾¯ 1iso𝑓 ( 𝑥) ¯   · · · 𝐾¯ 2iso𝑓 ( 𝑥) ¯   ..  .. . .    · · · 𝑉¯ iso ( 𝑥) ¯  𝑓 (2.47) for which the diagonal terms are the centroid potential of mean force for each diabatic surface, ∫ 1 ∑︁ −𝛽𝑉¯𝑖iso ( 𝑥) ¯ 𝑒 = lim 𝐶 𝑑x 𝛿( 𝑥¯ − 𝑥𝛼 ) (2.48) 𝑛→∞ 𝑛 𝛼 " !# 𝑛 1 ∑︁ ×exp −𝛽 𝑈spr (x) + 𝑉𝑖 (𝑥 𝛼 ) 𝑛 𝛼=1 for 𝑖 = 1, . . . , 𝑓 , and the off-diagonal terms are given by  𝛽    2 𝑉¯𝑖iso ( 𝑥)+ ¯ 𝑉¯ 𝑗iso ( 𝑥) ¯ iso 2 2 𝐾¯ 𝑖 𝑗 ( 𝑥) ¯ = acosh 𝑒 𝜇¯𝑖 𝑗 ( 𝑥)/2 ¯ /𝛽2  2 ¯ − 𝑉¯ 𝑗iso ( 𝑥) ¯ /4, − 𝑉¯𝑖iso ( 𝑥) (2.49) 29 where ∫ 𝜇¯𝑖 𝑗 ( 𝑥) ¯ = lim 𝐶 𝑑x 𝛿( 𝑥¯ − 𝑛→∞ 1 ∑︁ 𝑥 𝛼 ) 𝑒 −𝛽𝑈spr (x) 𝜇𝑖 𝑗 (x), 𝑛 𝛼 (2.50) and 𝜇𝑖 𝑗 (x) is given by Eq. 2.19. Also included in the isomorphic potential is the many-body contribution,       1 𝜇( ¯ 𝑥) ¯   iso  , 𝑉¯many-body ( 𝑥) ¯ = − ln   ˆ iso 𝛽  ¯ ¯   tre 𝑒 −𝛽𝑉2-body (𝑥)     (2.51) which vanishes for the case of a two-level system. C. The positivity and evaluation of 𝜇 In the limit of large bead number, 𝜇 can be expressed as a continuous path integral # " 𝑛 Ö ( 𝛼) ˆ (2.52) lim 𝜇(x) = lim tre 𝑒 −𝛽𝑛𝑉 (𝑥 ) 𝑛→∞ 𝑛→∞ 𝛼=1   ∫ 𝛽  = tre exp (Ô) − 𝑉ˆ (𝑥(𝜏))𝑑𝜏 , 0 where exp (Ô) is the time-ordered exponential, which is needed since 𝑉ˆ (𝑥) may not commute with itself at different imaginary times along the path, 𝑥(𝜏). Application of the generalized cumulant expansion [76] to this time-ordered exponential yields ∞ © ∑︁ ª lim 𝜇(x) = exp ­ (−1) 𝑗 𝐾 𝑗 (x(𝜏)) ® , 𝑛→∞ « 𝑗=1 ¬ where 𝐾 𝑗 is the 𝑗 th -order cumulant ∫ 𝛽 ∫ 𝜏1 ∫ 𝜏𝑛−1 𝐾 𝑗 (x(𝜏)) = 𝑑𝜏1 𝑑𝜏2 · · · 𝑑𝜏 𝑗 0 0 0   ˆ ˆ tr(c) e 𝑉 (x(𝜏1 )) · · · 𝑉 (x(𝜏 𝑗 )) , (2.53) (2.54) and tr(c) e [·] is the cumulant partial trace defined in Eq. 2.9 of Reference [76]. Given that the exponand in Eq. 2.53 is thus a sum of real numbers, it follows that lim𝑛→∞ 𝜇(x) > 0, providing that the cumulant expansion series converge (which is true for most conventional dynamics systems). In practice, for the 𝑛-bead discretization of the path integral, both 𝜇 and its derivatives 𝜕𝜇/𝜕𝑥 (𝛼) are evaluated using Bell’s algorithm [77], which requires only O (𝑛) operations. Details of this algorithm are provided elsewhere [42, 58]. 30 ~c (t) xx 8 A 5 QM QCLE-classical QCLE-C-iso QCLE-Wigner 2 -1 0 5 10 15 10 15 t 〈 x2(t) 〉 6 B 4 2 0 0 5 t Figure 2.5: Time-evolution of the initial phase-space distribution under QCLE dynamics. (A) Kubo-transformed position-autocorrelation functions for two linearly coupled harmonic oscillators, with physical potential energy matrix elements of 𝑉1 (𝑥) = 12 (𝑥 − 𝑥 0 ) 2 , 𝑉2 (𝑥) = 12 (𝑥 + 𝑥0 ) 2 , and 𝐾12 (𝑥) = 1.25, where 𝑥 0 = 2, and 𝛽 = 8. (B) Time-evolution of the second moment of the phase-space distribution with respect to position, ⟨𝑥 2 (𝑡)⟩, for the system in panel A. Results are obtained using exact quantum dynamics (QM; black, dots), QCLE dynamics with nuclei initialized from the classical phase-space distribution on the physical potential (QCLE-classical; blue, dashed), QCLE dynamics with nuclei initialized from the classical phase-space distribution on the isomorphic Hamiltonian (QCLE-C-iso; red, solid), and QCLE dynamics with nuclei initialized from the multi-surface partial Wigner distribution (QCLE-Wigner; black, solid). 31 D. Time-evolution of the initial phase-space distribution under QCLE dynamics Here, we examine a source of error for the QCLE dynamics presented in Fig. 2.1EH of the main text. In particular, we quantify the extent to which the QCLE dynamics preserves the MQC phase-space distribution that arises in the classical limit for the nuclear degrees of freedom (Section 2.4). For a two-level system comprised of linearly coupled one-dimensional harmonic oscillators (see caption), Fig. 2.5A shows results for the Kubo-transformed position-autocorrelation function, and Fig. 2.5B shows the second moment of the time-evolved initial phase-space distribution with respect to position, ⟨𝑥 2 (𝑡)⟩. The results in Fig. 2.5A are similar to those discussed in Fig. 2.1E-H, with substantial errors emerging for both the QCLE-classical and QCLE-C-iso at lower temperature and higher coupling; as is necessary for the system studied in this appendix [64], the QCLE dynamics initialized from the multi-level partial Wigner distribution (QCLE-Wigner in Fig. 2.5A) recovers exact quantum mechanics. As is seen in panel B, the QCLE dynamics exactly preserves the second moment of the initial Wigner phase-space distribution for this system [54, 64], but it does not preserve the initial MQC phase-space distribution associated with either the physical potential (QCLE-classical) or the isomorphic potential (QCLE-C-iso). Indeed, the erroneous features in the time correlation functions in panel A coincide with non-conservation of the MQC phase-space distribution in panel B. Although use of an initial MQC phase-space distribution for the QCLE dynamics emerged (Section 2.4) from our goal of obtaining a classical limit for the nuclear degrees of freedom without double-counting of nuclear quantum effects from the initial distribution, it is clear that the MQC phase-space distribution is not conserved by the QCLE dynamics, leading to erroneous time correlations in both the QCLEclassical and QCLE-C-iso results. 32 Chapter 3 SIMPLE FLUX-SIDE FORMULATION OF STATE-RESOLVED THERMAL REACTION RATES FOR RING-POLYMER SURFACE HOPPING This chapter is based on the following publication: 1. Tao, X., Shushkov, P. & Miller III, T. F. Simple Flux-Side Formulation of State-Resolved Thermal Reaction Rates for Ring-Polymer Surface Hopping. Journal of Physical Chemistry A 123, 3013–3020. doi:10.1021/acs.jpca. 9b00877 (2019). 3.1 Abstract Employing the recently developed isomorphic Hamiltonian framework for including nuclear quantum effects in mixed quantum-classical non-adiabatic dynamics [J. Chem. Phys., 148, 102327 (2018)], we present a flux-side formulation of stateresolved thermal reaction rates for ring-polymer surface hopping (iso-RPSH). An appealing aspect of the new approach is that calculation of multiple state-resolved non-adiabatic thermal reaction rates is enabled with only a single free-energy surface calculation, whereas previous non-adiabatic flux-side formulations for surface hopping involve multiple free-energy surface calculations. The method is shown to be robust and straightforwardly implemented, and numerical results reveal that RPSH in the isomorphic Hamiltonian framework leads to better dividing-surface independence than alternative RPSH methods, due to improved preservation of the path-integral statistics. 3.2 Introduction Ring polymer molecular dynamics (RPMD) [20, 21, 30, 31, 78, 79] provides a simple method for including nuclear quantum effects (NQEs), such as zero-point energy and tunneling, in rate calculations for electronically adiabatic processes. RPMD is based on the imaginary time path-integral formalism [18, 19, 55] and captures NQEs by evolving classical trajectories of a ring-polymer Hamiltonian that exactly preserve quantum Boltzmann statistics. The RPMD thermal reaction rate is rigorously independent of the choice of dividing surface for the reaction [79], numerically exact at high temperature limit and for a parabolic barrier [78], and 33 connected with instanton theory in the deep tunneling region [25]. These properties, in combination with favorable scaling of the computational cost of RPMD with respect to system size, have led it to be widely applied for the study of thermal reaction rates and mechanisms in complex systems [3, 33, 44, 45, 80–90]. Recent effort has focused on extending RPMD for the description of processes involving multiple electronic states [39, 40, 42, 43, 57, 58, 91, 92]. To this end, we recently introduced a multi-state isomorphic ring-polymer Hamiltonian [93] that allows for the inclusion of NQEs in any mixed quantum-classical (MQC) non-adiabatic dynamics method, including surface hopping [13], Ehrenfest dynamics [12], and other widely used methods [53, 54, 61, 94, 95]. The isomorphic Hamiltonian was shown to provide a natural and accurate way to combine RPMD with fewest switches surface hopping (i.e., iso-RPSH) [93], although the previously reported thermal rate calculations were performed using the relatively inefficient method of calculating rates from direct dynamics (i.e., without rare-event sampling). Here, we address this inefficiency by developing the theory for state-resolved iso-RPSH thermal reaction rates in the flux-side formulation. The new approach has the appealing feature of requiring only a single free energy calculation to obtain the quantum transition state theory (QTST) rate, and the subsequent dynamical recrossing correction makes the rate prediction state-resolved and (nearly) dividing-surface independent. The method is illustrated in a model for the non-adiabatic F+H2 reaction and compared to other implementations of ring-polymer surface hopping. 3.3 Method We begin by reviewing the multi-state isomorphic ring-polymer Hamiltonian and its combination with fewest switches surface hopping [93]. The methodology is then extended for the calculation of state-resolved thermal rates in the flux-side formulation. The isomorphic Hamiltonian Consider a general Hamiltonian for a physical system with 𝑁 electronic states in the diabatic representation, 2 𝑝ˆ + V( 𝑞), ˆ 𝐻ˆ = 2𝑚 (3.1) where 𝑝, ˆ 𝑞ˆ and 𝑚 are the nuclear position, momentum, and mass, respectively. The physical potential, V, is expressed as an 𝑁 × 𝑁 Hermitian matrix where the diagonal terms correspond to the diabatic potential energy surfaces and the off-diagonal terms 34 correspond to electronic couplings. Throughout this paper, results are presented for a single nuclear degree of freedom, although extension to higher dimensions is straightforward. By considering the path integral discretization of the canonical partition function, 𝑄, we thus introduce the multi-state isomorphic ring-polymer Hamiltonian, Hiso 𝑛 (p, q) [93], such that h i ˆ 𝑄 = tr 𝑒 −𝛽 𝐻  𝑛  𝑛∫ h i iso = lim 𝑑p 𝑑q tre 𝑒 −𝛽H𝑛 (p,q) , (3.2) 𝑛→∞ 2𝜋ℏ where 𝛽 and 𝑛 are the inverse temperature and the number of imaginary time discretization steps, respectively. The vector q = {𝑞 1 , 𝑞 2 , . . . , 𝑞 𝑛 } represents the positions of the ring-polymer beads, and p represents the corresponding momenta. The symbol ‘tre ’ denotes a partial trace over the electronic subspace. The resulting isomorphic Hamiltonian describes a ring-polymer comprised of 𝑛 replicas, or beads, associated with copies of the original physical system, Hiso 𝑛 (p, q) = 𝑛 ∑︁ 𝑝2 𝛼 2𝑚 𝑛 𝛼=1 + 𝑈spr (q) + Viso 𝑛 (q) (3.3) with 𝑚 𝑛 ≡ 𝑚/𝑛. Neighboring beads are connected by harmonic springs, 𝑛 ∑︁ 1 2 (𝑞 𝛼 − 𝑞 𝛼+1 ) 2 , 𝑈spr (q) = 𝑚 𝑛 𝜔𝑛 2 𝛼=1 (3.4) where 𝜔𝑛 = (𝛽𝑛 ℏ) −1 , 𝛽𝑛 = 𝛽/𝑛, and following Reference [93], the isomorphic potential Viso 𝑛 is an 𝑁 × 𝑁 matrix that satisfies # " 𝑛 h i Ö iso 𝜇(q) ≡ tre 𝑒 −𝛽V𝑛 (q) = tre 𝑒 −𝛽𝑛 V(𝑞 𝛼 ) (3.5) 𝛼=1 and is chosen to have the form diag Viso 𝑛 = V𝑛 coup + V𝑛 + 𝑉𝑛mb I, (3.6) where I is the identity matrix. The diagonal elements of the isomorphic potential correspond to the usual ring-polymer potentials associated with each diabatic state, 𝑛 h diag V𝑛 (q) i 1 ∑︁ [V(𝑞 𝛼 )] 𝑖𝑖 = 𝛿𝑖 𝑗 𝑖𝑗 𝑛 𝛼=1 (3.7) 35 with 𝛿𝑖 𝑗 the Kronecker delta function. The potential coupling between pairs of states in the isomorphic potential is given by " h # i h i  diag diag n 𝛽 V𝑛 + V𝑛 /2  coup  2 𝑖𝑖 𝑗𝑗 V𝑛 = (1 − 𝛿𝑖 𝑗 ) acosh2 𝑒 𝜇𝑖 𝑗 /2 /𝛽2 𝑖𝑗 − h diag V𝑛 i − h diag V𝑛 i 2 𝑖𝑖 o /4 , (3.8) 𝑗𝑗 where 𝜇𝑖 𝑗 (q) = tre " 𝑛 Ö 𝛼=1 " exp −𝛽𝑛 [V(𝑞 𝛼 )] 𝑖𝑖 [V(𝑞 𝛼 )] 𝑖 𝑗 [V(𝑞 𝛼 )] 𝑗𝑖 [V(𝑞 𝛼 )] 𝑗 𝑗 # !# , (3.9) and the sign of the potential coupling is kept to be the same as that of the physical potential evaluated at the ring-polymer centroid position. Lastly, the many-body potential term ensures that path-integral isomorphism (Eq. 3.2) is rigorously satisfied,             1 𝜇(q) mb   diag 𝑉𝑛 (q) = − ln (3.10)  . coup 𝛽  −𝛽 V𝑛 (q)+V𝑛 (q)       tre 𝑒    Like the standard ring-polymer Hamiltonian for a single-level system [20], the Hamiltonian in Eq. 3.3 exploits the quantum-classical ring-polymer isomorphism, in the sense that Boltzmann sampling in the classical limit for the nuclear degrees of freedom on Eq. 3.3 yields the exact quantum Boltzmann distribution for the physical system. iso-RPSH An appealing feature of the isomorphic Hamiltonian (Eq. 3.3) is that it immediately can be combined with any mixed quantum-classical (MQC) dynamics method, thereby including nuclear quantum effects. The case for which the MQC method is trajectory surface hopping [13] leads to the iso-RPSH method [93]. For iso-RPSH, the evolution of the nuclei is described via the extended classical equations of motion 𝑚 𝑛 𝑞¥𝛼 = 𝑚 𝑛 𝜔2𝑛 (𝑞 𝛼−1 + 𝑞 𝛼+1 − 2𝑞 𝛼 ) − 𝜕 iso 𝜀 (q), 𝜕𝑞 𝛼 𝛾 (3.11) th where 𝜀 iso 𝛾 (q) is the 𝛾 adiabatic potential energy surface that is obtained by diagonalizing the isomorphic potential Viso 𝑛 (q), and the evolution of the electronic wavefunction, 𝜓, is described using the time-dependent Schr¥odinger equation along the nuclear trajectories, 𝑖ℏ 𝜕 𝜓(q, 𝑡) = Viso 𝑛 (q) 𝜓(q, 𝑡). 𝜕𝑡 (3.12) 36 Figure 3.1: Schematic illustration of a hopping ring-polymer near avoided crossing. In RPSH dynamics, ring-polymer moves along a single surface until interrupted by instantaneous state transitions. The ring-polymers hop as a whole. As illustrated in Fig. 3.1, the ring-polymer evolves on a particular adiabatic potential energy surface until all ring-polymer beads simultaneously hop to a different adiabatic potential energy surface with probability ) (   ∑︁ 2 (3.13) Re 𝑎 𝛾𝜁 𝑑 𝜁 𝛾 𝑞¤ 𝛼 Δ𝑡, 0 . 𝑝 𝛾→𝜁 = max − 𝑎 𝛾𝛾 𝛼 Here, 𝑎 𝛾𝜁 and 𝑑 𝜁 𝛾 are the elements of the adiabatic electronic density matrix and firstderivative non-adiabatic coupling, respectively, and Δ𝑡 is the propagation timestep. A hop can only occur under the condition of energy conservation, and hops are rejected if the ring polymer has insufficient velocity component along the nonadiabatic coupling vector; exactly as in classical surface hopping [13], the velocity of each nuclear degree of freedom is modified along the direction of the non-adiabatic coupling vector. Although the current paper only presents results obtained using the Tully-surface hopping algorithm in its simplest form [13], modified algorithms [69, 70, 96] including momentum reversal, decoherence corrections, and extension to Liouvillian space can also be straightforwardly and unambiguously implemented in the isomorphic Hamiltonian framework. 37 Flux-side formulation of the iso-RPSH thermal rate By extension of classical surface hopping rate theory [97], the flux-side formulation of the RPSH thermal reaction rate is [41]  𝑛  𝑛∫
RPSH 𝑘 𝑖→ 𝑓 𝑄 𝑟 = lim lim 𝑑p0 𝑑q0 𝜌𝑛 p0 , q0 𝑛→∞ 𝑡→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑞¯𝑡 − 𝑞 ‡ ) 𝛿𝑖𝜃 −𝑡 𝛿 𝑓 𝜃 𝑡 . (3.14) Here, 𝑖 and 𝑓 specify the state upon which the system enter and exit, respectively, 𝑄 𝑟 is the reactant partition function, and 𝑡 → ∞ defines a time period that is sufficiently long to allow for dynamical recrossing. The ring-polymer centroid position at Í time 𝑡 is given by 𝑞¯𝑡 , and 𝑣¯ 0 denotes the initial centroid velocity, 𝑣¯ 0 = 𝛼 𝑞¤ 0,𝛼 /𝑛. 𝑞 ‡ indicates the position of the dividing surface that separates the reactant and product regions in position space, and 𝜌𝑛 (p, q) is the path-integral representation of the equilibrium quantum Boltzmann distribution for the multi-level system. The Heaviside function, ℎ(𝑞 − 𝑞 ‡ ), is non-zero only if the position of the system is in the product region, and the Kronecker 𝛿 function, 𝛿𝑖𝜃 , is non-zero only if the system resides on adiabatic surface 𝑖, with 𝜃 𝑡 specifying the active adiabatic surface upon which the RPSH trajectory resides at time 𝑡. In the earlier work [41], 𝜌𝑛 was approximated using a Boltzmann-weighted sum of adiabatic ring-polymer surfaces, ∑︁ 𝑒 −𝛽𝐻𝑛,𝛾 (p,q) , 𝜌𝑛 (p, q) ≃ (3.15) 𝛾 where 𝐻𝑛,𝛾 (p, q) = 𝑛 ∑︁ 𝑝2 𝛼 2𝑚 𝑛 𝛼=1 𝑛 + 𝑈spr (q) + 1 ∑︁ 𝜀 𝛾 (𝑞 𝛼 ), 𝑛 𝛼=1 (3.16) and 𝜀 𝛾 (𝑞 𝛼 ) is the 𝛾 th eigenvalue of the physical potential V(𝑞 𝛼 ). A quantum transition state theory (QTST) rate [98, 99] was computed for each adiabatic potential energy surface,  𝑛  𝑛∫ nokinks ‡ 𝑘 QTST,𝛾 (𝑞 )𝑄 𝑟 = lim 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛,𝛾 ( p0 ,q0 ) 𝑛→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑝¯0 ), to yield the state-resolved thermal rate, ∑︁ nokinks-RPSH nokinks 𝑘 𝑖→ = 𝑘 QTST,𝛾 (𝑞 ‡ ) 𝜅𝑖→ 𝑓 ,𝛾 (𝑞 ‡ ), 𝑓 𝛾 (3.17) (3.18) 38 where 𝜅𝑖→ 𝑓 ,𝛾 is the recrossing factor [41] associated with trajectories that are initialized at the dividing surface on adiabat 𝛾 and that reside on adiabat 𝑖 as 𝑡 → −∞ and on adiabat 𝑓 as 𝑡 → ∞. Here, the superscript “nokinks” is included to indicate the approximation to the path-integral statistics made in Eq. 3.15. Note that to calculate a single state-resolved rate, this protocol requires the independent calculation of 𝑁 free energy surfaces and 𝑁 recrossing factors. In the current work, the quantum Boltzmann distribution is expressed exactly in terms of the isomorphic Hamiltonian, h i iso 𝜌𝑛 (p, q) =tre 𝑒 −𝛽H𝑛 (p,q) (3.19) T =𝑒 −𝛽p p/2𝑚 𝑛 𝑒 −𝛽𝑈spr (q) 𝜇(q) ∑︁ iso = 𝑒 −𝛽𝐻𝑛,𝛾 (p,q) , (3.20) (3.21) 𝛾 where iso 𝐻𝑛,𝛾 (p, q) = 𝑛 ∑︁ 𝑝2 𝛼 𝛼=1 2𝑚 𝑛 + 𝑈spr (q) + 𝜀 iso 𝛾 (q) , (3.22) iso th and 𝜀 iso 𝛾 (q) is the 𝛾 eigenvalue of the isomorphic potential, V𝑛 (q). The analog of Eqs. 3.17 and 3.18 is still applicable with QTST rate now defined in terms of the adiabats of the isomorphic Hamiltonian, such that  𝑛  𝑛∫ iso ‡ iso 𝑘 QTST,𝛾 (𝑞 )𝑄 𝑟 = lim 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛,𝛾 ( p0 ,q0 ) 𝑛→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑝¯0 ). The resulting expression for the state-resolved thermal rate is ∑︁ iso-RPSH iso iso ‡ 𝑘 𝑖→ = 𝑘 QTST,𝛾 (𝑞 ‡ ) 𝜅𝑖→ 𝑓 ,𝛾 (𝑞 ), 𝑓 (3.23) (3.24) 𝛾 iso ‡ where the 𝜅𝑖→ 𝑓 ,𝛾 (𝑞 ) are the same state-resolved recrossing factors as in Eq. 3.18, except with trajectories run on the adiabats of the isomorphic Hamiltonian. As in Eq. 3.18, for a single state-resolved rate calculation, Eq. 3.24 requires the independent calculation of 𝑁 free energy surfaces and 𝑁 recrossing factors. A central result of the current work is that, by combining Eqs. 3.23 and 3.24 with Eq. 3.20, the state-resolved thermal rate expression can be formulated more simply, iso-RPSH MF iso ‡ 𝑘 𝑖→ = 𝑘 QTST (𝑞 ‡ ) 𝜅𝑖→ 𝑓 (𝑞 ), 𝑓 (3.25) 39 in terms of a state-unresolved (i.e., mean-field) QTST,  𝑛  𝑛∫ MF MF ‡ 𝑘 QTST (𝑞 )𝑄 𝑟 = lim 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛 ( p0 ,q0 ) 𝑛→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑝¯0 ), (3.26) where 𝐻𝑛MF (p, q) = 𝑛 ∑︁ 𝑝2 𝛼 𝛼=1 2𝑚 𝑛 + 𝑈spr (q) − 1 ln 𝜇(q) 𝛽 (3.27) is the familiar mean-field (MF)-RPMD Hamiltonian [42, 57, 58]. In this formulation, state-resolution of the rate appears only via the dynamical recrossing factor, ∫
lim lim 𝑑p0 𝑑q0 𝜌𝑛 p0 , q0 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑞¯𝑡 − 𝑞 ‡ ) 𝛿𝑖𝜃 −𝑡 𝛿 𝑓 ,𝜃 𝑡 𝑛→∞ 𝑡→∞ iso ‡ ∫ 𝜅𝑖→ . 𝑓 (𝑞 ) =
‡ lim 𝑑p0 𝑑q0 𝜌𝑛 p0 , q0 𝛿( 𝑞¯ 0 − 𝑞 ) 𝑣¯ 0 ℎ( 𝑝¯0 ) 𝑛→∞ (3.28) Practical evaluation of the dynamical recrossing factor in Eq. 3.28 involves initialization of the iso-RPSH trajectories from the quantum Boltzmann distribution with the ring-polymer centroid constrained to the dividing surface. First, using MF that 𝜌𝑛 (p, q) = 𝑒 −𝛽𝐻𝑛 (p,q) , nuclear configurations and velocities are sampled from MF the constrained distribution 𝑒 −𝛽𝐻𝑛 ( p0 ,q0 ) 𝛿( 𝑞¯ 0 − 𝑞 ‡ ); then, using Eq. 3.21, the electronic state at the dividing surface, 𝛾, is conditionally sampled from the Boltzmann weight of the isomorphic adiabatic surfaces at that ring-polymer nuclear configuration. The initialized trajectories are propagated both backwards and forwards in time using the iso-RPSH Hamiltonian in the adiabatic representation (Eq. 3.22) and with transitions induced by non-adiabatic couplings that are derived from Eq. 3.3. This leads to characterization of the position-space recrossing, as well as determination of the initial and final electronic state for each trajectory. The primary advantage of the rate expression in Eq. 3.25 is that evaluation of the state-resolved reaction rates requires only the single MF-QTST in Eq. 3.26 and thus only a single free-energy calculation with respect to the well-behaved MF Hamiltonian in Eq. 3.27. This is a distinct feature of the method proposed here. It is worth noting, however, that although the rate expression in Eq. 3.25 provides importance sampling with respect to the ring polymer nuclear configurations, it does not involve importance sampling with respect to the electronic states. That is, 40 the electronic state of the trajectory at the dividing surface is sampled directly from Í iso the conditional probability distribution 𝛾 𝑒 −𝛽𝐻𝑛,𝛾 ( p0 ,q0 ) 𝛿( 𝑞¯ 0 − 𝑞 ‡ ), such that more trajectories will be initialized from electronic states that are thermally accessible at the dividing surface. As a result, the rate expression in Eq. 3.25 will provide a simple and efficient avenue to obtaining the state-resolved reaction rates that are most thermally accessible, whereas the rate expression in Eq. 3.24 may provide a more efficient means of obtaining state-resolved rates that make small contributions to the total thermal reaction rate. Finally, we note that in direct analogy to Eqs. 3.25-3.27, the state-resolved thermal rate for the no-kinks approximation to RPSH can be expressed in terms of a stateunresolved MF-QTST, so that it too can be obtained with only a single free-energy calculation, nokinks-RPSH MF-nokinks ‡ nokinks ‡ (𝑞 ) 𝜅𝑖→ 𝑓 (𝑞 ), 𝑘 𝑖→ = 𝑘 QTST 𝑓 (3.29) in terms of a no-kinks approximation to the MF-QTST,  𝑛  𝑛∫ MF-nokinks p ,q MF-nokinks ‡ ( 0 0) (𝑞 )𝑄 𝑟 = lim 𝑘 QTST 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛 𝑛→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑝¯0 ), (3.30) ! ∑︁ 1 𝐻𝑛MF-nokinks (p, q) = − ln 𝑒 −𝛽𝐻𝑛,𝛾 (p,q) . 𝛽 𝛾 (3.31) where For the evaluation of the recrossing factor in Eq. 3.29, the initialized trajectories are propagated both backwards and forwards in time using the nokinks-RPSH Hamiltonian in the adiabatic representation (Eq. 3.16)) and with transitions induced by non-adiabatic couplings that are averaged over the bead positions (i.e., the “bead approximation” from Reference [41]). 3.4 Computational Details Calculations are presented for the one-dimensional, two-state model of gas-phase F+H2 reactive scattering introduced in Reference [93], which exhibits both substantial nuclear tunneling and electronically non-adiabatic effects. 41 Table 3.1: Parameter values for the potential of two-level reactive scattering model in Eq. 3.32. Parameter Value 𝐴1 7 𝐴2 −18/𝜋 𝐴3 0.25 𝐵1 −0.75 𝐵2 54/𝜋 Parameter Value 𝑎1 √ 1 𝑎2 3𝜋/4 𝑎3 0.25 𝑞1 −1.6 𝑞3 −2.625 Potential energy surfaces and couplings for the system in the diabatic representation are 𝑉1 (𝑞) = 𝐴1 + 𝐵1 −𝑎 1 + 𝑒 1 (𝑞−𝑞1 ) 𝑉2 (𝑞) = 𝐴2 𝐵2 + −𝑎 𝑞 2 2 1+𝑒 4 cosh (𝑎 2 𝑞/2) 𝐾12 (𝑞) = 𝐴3 𝑒 −𝑎3 (𝑞−𝑞3 ) 2 (3.32) with parameters given in Table 3.1. Both diabatic and adiabatic surfaces are plotted in Fig. 3.2A, with the reactant region at 𝑞 → −∞ and the product region at 𝑞 → ∞. For the computed state-resolved thermal reaction rates, we will use the notation that 𝑘 1 indicates the rate for the system that entering on adiabat 1 and exiting to the product region on adiabat 1, and 𝑘 2 indicates the rate for the system that entering on adiabat 2 and exiting to the product region on adiabat 1. Note that the reaction channels are named here differently from those in Reference [93]. We indicate the total (state-unresolved) thermal reaction rate as 𝑘 total , mainly consisting of contributions from above two channels at the investigated temperature range 𝛽 ≥ 1. All the rate results below are reported in the unit of m·molecule−1 s−1 . Three different surface hopping rate calculations are compared: (i) RPSH with the isomorphic Hamiltonian, called iso-RPSH, using Eqs. 3.25-3.27, (ii) RPSH with the no-kinks Hamiltonian, called nokinks-RPSH, using Eqs. 3.29-3.31, and (iii) classical surface hopping using the (shared) 1-bead limit of these expressions. We note that iso-RPSH employs an exact expression for the quantum Boltzmann distribution, nokinks-RPSH employs an approximate for the quantum Boltzmann distribution, and classical surface hopping employs the classical Boltzmann distribution. Exact quantum results are provided using wavepacket propagation the split-operator technique based [73, 100]. 42 Figure 3.2: The potential and the quantum transition state theory rates of a gasphase scattering model. (A) Potential of the model system described as in Eq. 3.32. Diabatic and adiabatic states are shown in solid and dash lines, respectively. The filled area is the region in which following dividing surface independence tests are performed. (B) Optimized positions of the dividing surface as a function of reciprocal temperature is shown in the dotted line. Solid green lines depict the mean-field quantum transition state theory (MF-QTST, Eq. 3.26) rate on a log scale as a function of the dividing surface position for three different temperatures. Position of the reaction barrier on ground adiabat is indicated with a thin solid line. 43 Each of the surface-hopping rate calculations is implemented in two stages—a statistical and a dynamical part. As in previous RPMD simulations [81, 82], the QTST rate for a biomolecular reaction is calculated using the Bennett-Chandler procedure [101, 102]. An auxiliary dividing surface, 𝑞 ∗ , is placed in the reactant asymptotic region, and the QTST rate is then given by 𝑘 QTST (𝑞 ‡ ) 𝑘 QTST (𝑞 ∗ ) 1 𝜌𝑐 (𝑞 ‡ ) = (2𝜋𝛽𝑚) 1/2 𝜌𝑐 (𝑞 ∗ ) 𝑘 QTST (𝑞 ‡ ) = 𝑘 QTST (𝑞 ∗ ) (3.33) with ∫ 𝜌𝑐 (𝑞 ) = 𝑑q 𝑒 −𝛽𝑈spr (q) 𝜇(q) 𝛿( 𝑞¯ − 𝑞 ‡ ) ‡ (3.34) the MF centroid density. For all the results presented in this section, QTST rates are obtained with importance sampling in the path-integral representation. Monte Carlo (MC) samplings with both 8𝛽 and 16𝛽 ring-polymer beads are conducted to make sure imaginary-time discretization convergence is reached. Statistical convergence is achieved for all the simulations with up to 109 MC moves. Optimized dividing surfaces in QTST calculations are obtained with a scan accurate to 0.01 a.u. The dynamical recrossing factor 𝜅 is computed for each method by trajectory simulations. The ring-polymer Hamiltonian is integrated using the velocity Verlet algorithm [72], with the Liouville operator for nuclei motion factorized such that the free ring-polymer movements are solved analytically [30, 31]. The isomorphic Hamiltonian and its derivatives are calculated using Bell’s algorithm [93]. The electronic Schr¥odinger equation is evolved in the interaction representation with an unitary nuclear-position-coupled mid-point propagator [103] for four times during each nuclear timestep. The trajectories are initialized from a thermal distribution at the dividing surface. The thermal ensemble is generated using Metropolis algorithm while velocities are sampled from the Maxwell-Boltzmann distribution. Utilizing a standard reweighting strategy [97], trajectories are first propagated backward in time to determine the correct quantum amplitude at the dividing surface, before they are propogated forward in time. The trajectories are simulated using a timestep of 0.0001 a.u. Dynamics simulation for ring-polymers with both 8𝛽 and 16𝛽 beads have been performed to check the convergence of the path-integral discretization. Rate results below are obtained using 105 trajectories, although we emphasize that to within graphical accuracy, identical results could have been obtained using only 44 103 trajectories. For the calculations reported with the dividing surface placed away from the optimized position, 8×105 trajectories are performed to fully describe the trajectory recrossing at each dividing surface position. 3.5 Results We begin by demonstrating ring-polymer surface hopping (RPSH) in the flux-side formulation for a one-dimensional, two-state model of gas-phase F+H2 reactive scattering (Fig. 3.2A). Given that the flux-side calculations involve a dynamical recrossing correction (Eq. 3.28, for example) to a QTST, Fig. 3.2B illustrates how the QTST estimate in Eq. 3.26 for the total thermal reaction rate varies as a function of the choice of dividing surface over a range of temperatures. Consistent with previous observations of electronically adiabatic RPMD [79], the optimal position of the dividing surface shifts away from the position of the barrier-top for the groundstate adiabat at lower temperatures (𝛽 > 1), due to significant nuclear tunneling effects. However, in a strictly non-adiabatic effect, Fig. 3.2B also shows that the optimal position of the dividing surface shifts away from the barrier-top position at high temperatures (𝛽 < 1), due to the influence of the excited electronic state. Fig. 3.3 presents the state-resolved thermal reaction rate versus inverse temperature for two different reactive channels, obtained using exact quantum mechanics, classical surface hopping, iso-RPSH [93], and nokinks-RPSH [41]. As is clear from comparison of the exact quantum results and classical surface hopping, the nuclear tunneling effects accelerate the reaction rate and becomes more significant at lower temperatures. Flux-side formulation of iso-RPSH dynamics provides accurate stateresolved rates for the full range of temperatures, as seen previously when employed with the side-side formulation of the rate [93]. The key difference using the fluxside formulation in the current work is that substantially smaller statistical errors are achieved with 1000-fold fewer dynamical trajectories than in the previous study. The nokinks-RPSH method overestimates the rate for both channel, but the effect is much more stubstantial for rate 𝑘 1 , which corresponds to the system entering and exiting on the lower electronic adiabat; the deviation of nokinks-RPSH from exact results is due to underestimation of the influence of the low-lying excited state, which acts to suppress the rate (see for example, Fig. 2 in Reference [45]). We now explore the robustness of the RPSH methods to the choice of dividing surface position. The RPSH calculations in Fig. 3.3 were obtained using the fluxside formulation with the optimal choice of dividing surface for the associated 45 Figure 3.3: Arrhenius plots of state-resolved thermal rates versus reciprocal temperature for two reaction channels. More specifically, (A) for adiabat 1 → 1 and (B) for adiabat 2 → 1. Rates are calculated from exact quantum mechanics (QM, thick solid black), classical surface hopping (classical SH, dashed black), ring-polymer surface hopping dynamics with the no-kinks approximation (nokinksRPSH, dotted green) and RPSH in the isomorphic Hamiltonian framework (isoRPSH, solid red). Mean standard deviations for all the data points plotted are smaller than the size of markers. Dividing surfaces in the RPSH calculations are placed at the optimal positions. MF-QTST (i.e., Eq. 3.26 for iso-RPSH and Eq. 3.30 for the nokinks-RPSH). For the inverse temperature 𝛽 = 8, Fig. 3.4A presents total (state-unresolved) thermal 46 Figure 3.4: (Continued on the following page.) 47 Figure 3.4: Dividing surface dependence of thermal rates at 𝛽 = 8. More specifically, (A) for state-unresolved total rate, (B) and (C) for channels adiabat 1 → 1, and adiabat 2 → 1, respectively. Rates are calculated from exact quantum mechanics (QM, thick grey), classical surface hopping (classical SH, dash-dotted black), mean-field quantum transition state theory (MF-QTST, dashed red, Eq. 3.26), mean-field ring polymer molecular dynamics (MF-RPMD, dashed violet), and three implementations of ring-polymer surface hopping as described in the main text (RPSH, solid red, dashed blue and dotted green lines). For data points for which the error bars are not explicitly shown, the standard error is smaller than the size of marker. Optimal positions for the dividing surface are shown in dashed grey. reaction rate for the reaction, using a variety of methods, as a function of the dividing surface position. As is necessary [104], the exact quantum mechanical rate is strictly independent of the choice of dividing surface, and in stark contrast, the MF-QTST of Eq. 3.26 is exponentially dependent on the position of the dividing surface. The classical surface hopping result differ substantially from the quantum result due to the exclusion of nuclear tunneling effects, and as has been seen previously [68], they vary as a function of dividing surface position due to inexact preservation of the mixed quantum-classical Boltzmann distribution by the dynamical trajectories [59, 105], although the degree of variation (∼40%) is small on the scale of this plot. Both RPSH methods likewise vary slightly on this scale as a function of dividing surface position due to inexact preservation of the quantum Boltzmann distribution by the dynamical trajectories. Notably, the rate obtained using MF-RPMD is strictly independent of the dividing surface position due to exact preservation of the quantum Boltzmann distribution by the dynamical trajectories, although this method has the disadvantage of not enabling state-resolved rate calculations. Figs. 3.4B and C present the state-resolved thermal reaction rate for the adiabat 1 → 1 channel (𝑘 1 ) and the adiabat 2 → 1 channel (𝑘 2 ), computed using exact quantum mechanics, iso-RPSH (red), and nokinks-RPSH (green). Again, it is seen that the exact quantum results are strictly independent of the position of the dividing surface, and as in Fig. 3.3, it is seen that nokinks-RPSH (green) substantially overestimates the 𝑘 1 rate while iso-RPSH (red) provides the best estimate of both channels. Given the more narrow range of rates plotted in these panels, the variation of the RPSH methods is more clearly seen with both iso-RPSH and nokinks-RPSH varying by a similar degree (∼25%) as a function of dividing surface position. 48 Also shown in Figs. 4B and C is the result of a rate expression that uses MF-QTST with the exact quantum Boltzmann distribution (Eq. 26), while the recrossing factor is evaluated by trajectories that are initialized from the exact quantum Boltzmann distribution and are propagated using nokinks-RPSH dynamics. Whereas this hybrid method, analogs of which have been suggested in other contexts [57], improves the rate associated with the lower adiabat (𝑘 1 ), it does so at the cost of introducing a much stronger dependence on the position of the dividing surface. This example illustrates the statistical importance of kinked ring-polymer configurations in describing these non-adiabatic reaction rates, as well as the fact that iso-RPSH does a substantially better job of preserving the true quantum Boltzmann distribution than the nokinksRPSH method. 3.6 Summary The development of practical, robust, and accurate quantum methods that incorporate both nuclear quantum effects and allow for electronically non-adiabatic transitions remains an important challenge for the description of complex photochemical, charge-transfer, and energy-transfer processes. In this work, we present an efficient flux-side formulation of ring-polymer surface hopping on the isomorphic pathintegral Hamiltonian (iso-RPSH). The method is straightforwardly implemented, involving no greater complexity than a standard trajectory surface hopping rate calculation, and it is shown to be accurate and relatively insensitive to the position of the transition-state dividing surface. A novel and appealing aspect of the flux-side formulation presented here is that calculation of multiple state-resolved non-adiabatic thermal reaction rates is enabled with only a single free-energy surface calculation. The combined accuracy and simplicity of iso-RPSH make it amenable to useful application in realistic chemical systems. 49 Chapter 4 MICROCANONICAL RATES FROM RING-POLYMER MOLECULAR DYNAMICS: DIRECT-SHOOTING, STATIONARY-PHASE, AND MAXIMUM-ENTROPY APPROACHES This chapter is based on the following publication: 1. Tao, X., Shushkov, P. & Miller III, T. F. Microcanonical Rates from RingPolymer Molecular Dynamics: Direct-Shooting, Stationary-Phase, and MaximumEntropy Approaches. Journal of Chemical Physics 152, 124117. doi:10 . 1063/1.5144307 (2020). 4.1 Abstract We address the calculation of microcanonical reaction rates for processes involving significant nuclear quantum effects using ring-polymer molecular dynamics (RPMD), both with and without electronically non-adiabatic transitions. After illustrating the shortcomings of the naive free-particle direct-shooting method, in which the temperature of the internal ring-polymer modes is set to the translational energy scale, we investigate alternative strategies based on the expression for the microcanonical rate in terms of the inverse Laplace transform of the thermal reaction rate. It is shown that simple application of the stationary-phase approximation (SPA) dramatically improves the performance of the microcanonical rates using RPMD, particularly in the low-energy region where tunneling dominates. Using the SPA as a Bayesian prior, numerically exact RPMD microcanonical rates are then obtained using maximum entropy inversion of the thermal reaction rates, for both electronically adiabatic and non-adiabatic model systems. Finally, the directshooting method is revisited using the SPA-determined temperature for the internal ring-polymer modes, leading to a simple, direct-simulation method with improved accuracy in the tunneling regime. This work suggests a general strategy for the extraction of microcanonical dynamical quantities from RPMD (or other approximate thermal) simulations. 50 4.2 Introduction Ring-polymer molecular dynamics (RPMD) [20, 21] has proven to be a useful tool for the calculation of chemical reaction rates [78, 79], spectra [26, 32], and transport coefficients [30, 31]. The method has been widely applied for the study of electronically adiabatic processes for which nuclear quantum effects play an important role [29, 33, 44, 45, 80, 81, 83, 86, 87, 89, 90, 106–116], and extensions of the method for systems involving electronically non-adiabatic processes are increasingly common [3, 39–43, 57, 58, 91–93, 117]. However, despite the utility of RPMD for calculating quantities in terms of thermal transport coefficients, less work has focused on the extension of the method to non-thermal initial distributions [118] or for the calculation of properties associated with non-thermal ensembles, such as microcanonical reaction rates, which would be of use for both benchmarking and practical applications. Application of RPMD beyond the canonical ensemble immediately encounters the question of how to treat the temperature associated with the intra-bead ring-polymer potential. This temperature is well-defined in thermal applications for which RPMD was initially developed [20, 21], and it has been justified for RPMD with particular non-equilibrium initial conditions [118]. In the context of microcanonical reaction rates, a direct-shooting method based on the free-particle temperature has been proposed [119], in which the internal ring-polymer temperature is fixed based on the microcanonical energy, i.e. 𝑇 = 𝐸/𝑘 B . This protocol has been employed in several model calculations [46, 119], although its reliability has not been systematically examined. The current work addresses the challenge of microcanonical rate calculations using RPMD. In addition to analyzing the previously proposed free-particle directshooting protocol, we introduce alternative stationary-phase and maximum-entropy inversion methods to extract microcanonical rates from thermal reaction rates, the calculation of which is well established using RPMD. Finally, we return to the direct-shooting method for microcanonical RPMD rates, replacing the free-particle temperature with the optimal temperature from the stationary-phase inversion, which is shown to yield greatly improved microcanonical RPMD rates in the low-energy regime. Numerical examples of these microcanonical RPMD methods are presented for both electronically adiabatic and non-adiabatic systems. 51 4.3 Methods Thermal reaction rates from RPMD We begin by briefly reviewing RPMD and its use for the calculation of thermal reaction rates. The theory is presented for a one-dimensional system, and extension to multiple dimensions is straightforward. Consider an electronically adiabatic system with the Hamiltonian 2 𝑝ˆ 𝐻ˆ = + 𝑉 ( 𝑞), ˆ 2𝑚 (4.1) where 𝑉 ( 𝑞) ˆ is the potential energy function. Expressing the quantum canonical partition function, 𝑄, in the path-integral representation yields [18, 19, 55] h i ˆ 𝑄 = tr 𝑒 −𝛽 𝐻  𝑛  𝑛∫ iso 𝑑p 𝑑q 𝑒 −𝛽𝐻𝑛 (p,q) , (4.2) = lim 𝑛→∞ 2𝜋ℏ where 𝛽 and 𝑛 are the reciprocal temperature and the number of imaginary time discretization steps, respectively; q = {𝑞 1 , 𝑞 2 , . . . , 𝑞 𝑛 } denotes the positions of the ring-polymer beads, and p denotes the corresponding momenta. Eq. 4.2 introduces the classical isomorphic ring-polymer Hamiltonian, 𝐻𝑛iso (p, q) = 𝑛 ∑︁ 𝑝2 𝛼 2𝑚 𝑛 𝛼=1 𝑛 + 𝑈spr (q) + 1 ∑︁ 𝑉 (𝑞 𝛼 ), 𝑛 𝛼=1 (4.3) with 𝛽𝑛 = 𝛽/𝑛, 𝑚 𝑛 = 𝑚/𝑛, and neighboring ring-polymer beads are connected via harmonic springs 𝑛 1 𝑚 𝑛 ∑︁ (𝑞 𝛼 − 𝑞 𝛼+1 ) 2 . 𝑈spr (q) = 2 2 𝛽𝑛 𝛼=1 (4.4) Classical sampling of the ring-polymer Hamiltonian faithfully preserves quantum Boltzmann statistics. The classical equations of motion associated with the ring polymer Hamiltonian are given by 𝑞¥𝛼 = 1 1 𝜕𝑉 (𝑞 𝛼 ) (𝑞 𝛼+1 + 𝑞 𝛼−1 − 2𝑞 𝛼 ) − . 2 𝑚 𝜕𝑞 𝛼 𝛽𝑛 (4.5) 52 The calculation of thermal rates from RPMD then simply follows from the application of classical rate theory to the dynamics associated with the ring-polymer Hamiltonian [21, 78, 79]. Specifically, calculation of the thermal RPMD rate in the flux-side formulation yields  𝑛 𝑛 ∫ iso 𝑘𝑄 𝑟 = lim lim 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛 (p0 ,q0 ) 𝑛→∞ 𝑡→“∞” 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) 𝑣¯ 0 ℎ( 𝑞¯𝑡 − 𝑞 ‡ ), (4.6) which correlates the positions and velocities of the ring-polymer beads at time 𝑡 following evolution according to the ring-polymer equations of motion (Eq. 4.5) from an initial distribution in which the ring-polymer centroid is positioned at the dividing surface for the reaction. Here, 𝑞¯ 0 and 𝑞¯𝑡 indicate the ring-polymer centroid position at time zero and 𝑡, respectively, and 𝑣¯ 0 indicates the centroid velocity at time zero. 𝑄 𝑟 denotes the reactant partition function, 𝑞 ‡ indicates the position of the dividing surface that separates the reactant and product, ℎ is the Heaviside function, and 𝛿 is the Dirac delta function. Microcanonical reaction rates from RPMD In the following, we describe three alternative strategies for calculating microcanonical reaction rates from RPMD. The first involves an inverse Laplace transform of the thermal RPMD reaction rates and introduces no approximations beyond that of the thermal RPMD rate theory, although it is numerically the most demanding. The subsequent two methods introduce additional approximations (i.e., the stationary phase approximation and the direct shooting approximation) with the benefit of reduced numerical complexity. Maximum entropy inversion Reaction rates in the microcanonical1 and canonical ensembles are connected via the Laplace transform [120] ∫ +∞ 1 𝑑𝐸 𝑒 −𝛽𝐸 𝑁 (𝐸), (4.7) 𝑘 (𝛽)𝑄 𝑟 (𝛽) = 2𝜋ℏ −∞ which can be inverted to yield [121–123] ∫ 𝛾+𝑖∞ 1 𝑁 (𝐸) = (2𝜋ℏ) 𝑑𝛽 𝑒 Φ(𝛽) (4.8) 2𝜋𝑖 𝛾−𝑖∞ 1 Following previous work [120], we refer to 𝑁 (𝐸) as the microcanonical rate, although it is also known as the cumulative reaction probability. The quantity 𝑟 (𝐸) = 𝑁 (𝐸)/[2𝜋ℏ 𝜌(𝐸)] is also commonly referred to as the microcanonical rate, where 𝜌(𝐸) is the density of states; by estimating 𝜌(𝐸) with statistical mechanical methods, the quantities can be interconverted. 53 with Φ(𝛽) = 𝛽𝐸 + log (𝑘𝑄 𝑟 ) . (4.9) The line integral in Eq. 4.8 is performed along Re[𝛽] = 𝛾, where 𝛾 is greater than the real part of all points for which Φ is singular. However, numerical implementation of this inverse Laplace transform is typically ill-conditioned and sensitive to statistical noise [65], which is unavoidable in simulation-based thermal rate calculations. To ameliorate this problem, we first employ the maximum entropy (MaxEnt) method [124, 125], which utilizes statistic inference and a Bayesian prior to regularize the numerical inversion [89, 125–128]. MaxEnt is implemented by rewriting the integral in Eq. 4.7 in matrix form, 𝜿 = B 𝝂, (4.10) where 𝜿 is the vector of thermal rate input data at discrete temperature points {𝛽𝑖 }, and 𝝂 is the vector of microcanonical rate outputs at discrete energy values {𝐸 𝑗 }. Specifically, the elements of 𝜿 and 𝝂 are 𝜅𝑖 = 2𝜋ℏ 𝑘 (𝛽𝑖 )𝑄 𝑟 (𝛽𝑖 ) and 𝜈 𝑗 = 𝑁 (𝐸 𝑗 ). The matrix B is comprised of the Boltzmann kernel 𝐵𝑖 𝑗 = 𝑒 −𝛽𝑖 𝐸 𝑗 Δ𝐸 𝑗 , where Δ𝐸 𝑗 = 𝐸 𝑗+1 − 𝐸 𝑗 is the integration stepsize. MaxEnt yields the microcanonical rate by maximizing the objective function 𝑄(𝝂; 𝛼) = 𝛼𝑆(𝝂) − 𝜒2 (𝝂)/2 + 𝑉reg (𝝂), (4.11) where the information entropy 𝑆 describes the degree to which solution is faithful to a prior model 𝝀({𝐸 𝑗 }),  ∑︁  𝜈𝑗 , (4.12) 𝑆(𝝂) = 𝜈 𝑗 − 𝜆 𝑗 − 𝜈 𝑗 log 𝜆 𝑗 𝑗 and the likelihood function 𝜒2 describes the accuracy with which the reference thermal rate data is fit, 𝜒2 (𝝂) = (𝜿 − B𝝂) T C−1 (𝜿 − B𝝂). (4.13) Here, C is the covariance matrix for the thermal rate data with elements 𝐶𝑖𝑖 ′ =𝛿𝑖𝑖 ′ 𝜎𝑖2 , (4.14) where 𝛿𝑖𝑖 ′ is the Kronecker delta function, and 𝜎𝑖 is the standard deviation for the 𝑖-th thermal rate datapoint. The parameter 𝛼 balances between accurately fitting the reference data while preserving the prior model. 54 Finally, 𝑉reg (𝝂) penalizes those trial solutions that violate the physical constraints of the microcanonical rate constant, namely that it satisfy 𝑁 (𝐸) ∈ [0, 1]. The lower bound is enforced by conducting a solution search only in the positive subspace while the upper bound is enforced via the functional form ∑︁ 1 𝜁 𝐼 2 (𝜈 𝑗 − 1), (4.15) 𝑉reg (𝝂) = − 2 𝑗 where 𝐼 (𝜈 𝑗 − 1) =    𝜈 𝑗 − 1,  𝜈 𝑗 ≥ 1, (4.16)   0, 𝜈 𝑗 < 1,  and 𝜁 is chosen according to a tolerance criterion. Stationary phase approximation As an alternative to numerically exact inversion, we apply the stationary phase approximation (SPA) [129] to Eq. 4.8. Implementation of the SPA involves finding the stationary point of the phase function Φ in Eq. 4.9 and then approximating the integrand as a Gaussian function along the imaginary axis. Setting the first-order derivative of the phase function to zero yields the energy-temperature correspondence [130, 131] 𝐸 st = − 𝜕 log [ 𝑘 (𝛽)𝑄 𝑟 (𝛽) ] , 𝜕𝛽 𝛽st 𝛽st ∈ R (4.17) where 𝛽st is the stationary temperature that is assumed to dominate the integrand. The resulting SPA microcanonical rate prediction is given by ! −1/2 2𝜋ℏ 𝜕 2 log [ 𝑘 (𝛽)𝑄 𝑟 (𝛽) ] 𝑁SPA (𝐸 st ) = √ 𝜕 𝛽2 2𝜋 𝛽st × 𝑒 𝛽st 𝐸st 𝑘 (𝛽st )𝑄 𝑟 (𝛽st ). (4.18) We note that semiclassical instanton theory can also yield approximate microcanonical rates via the SPA [132–134]. A well-known shortcoming of the SPA is that the calculated microcanonical rate violates the constraint 𝑁 (𝐸) ≤ 1 in the high-energy limit [131]. In this regime, the barrier-crossing dynamics reduces to free particle motion, and the thermal rate becomes [𝑘𝑄 𝑟 ] FP = 1/(2𝜋𝛽ℏ) . (4.19) 55 Substituting Eq. 4.19 into Eqs. 4.17 and 4.18 yields the energy-temperature correspondence relation 𝐸 stFP = 1/𝛽stFP , (4.20) and the corresponding microcanonical rate in the high-energy limit is √ FP 𝑁SPA (𝐸) = 𝑒/ 2𝜋 ≃ 1.084, (4.21) in excess of the correct upper limit. Direct shooting approximation By analogy with classical rate theory, a physically intuitive strategy for approximating microcanonical rates from RPMD is to simply (i) initialize trajectories from the reactant side with specified translational energy, (ii) propagate those trajectories using the microcanonical equations of motion in Eq. 4.5, and (iii) count the proportion of trajectories that reach the product region, such that 𝑁direct (𝐸, 𝛽int ) iso 𝑒 −𝛽int 𝐻𝑛 (p0 ,q0 ) = lim lim 𝑑p 𝑑q 0 0 𝑛→∞ 𝑡→“∞” (2𝜋ℏ) 𝑛−1 𝑒 −𝛽int 𝐸 h i √︁ × 𝛿 𝑝¯0 − 2𝑚(𝐸 − 𝑉𝑎 ) 𝛿( 𝑞¯ 0 − 𝑞 ‡ ) ℎ( 𝑞¯𝑡 − 𝑞 ‡ ), 𝑛𝑛 ∫ (4.22) Í where 𝑝¯ = 𝛼 𝑝 𝛼 is the centroid momentum, and 𝑉𝑎 is the potential energy in the reactant asymptotic region. The centroid kinetic energy for the RPMD trajectories are initialized to match the physical incident energy (as indicated by the 𝛿 function), while the internal modes are thermally sampled from an internal temperature 𝛽int , the appropriate value of which is not obvious. Previously [119], the direct shooting method for calculating microcanonical rates has been applied with the internal temperature set in correspondence to the physical incident energy 𝛽int = 1/𝐸, which we call the free-particle protocol; in the current work, we shall also consider a prescription for the internal temperature that is derived from the SPA. We note that direct shooting is similar in practical implementation to the calculation of non-equilibrium time-correlation functions using RPMD with momentum-kick initial conditions [106, 118], although the theoretical justification is more clearly established for the case of non-equilibrium time-correlation functions than for the calculation of microcanonical rates as considered here. 56 4.4 Computational Details Unless specified, all results are reported in atomic units. Implementation of the direct-shooting approach for microcanonical rates employs Eq. 4.22. Initial configurations for the ring polymer are sampled from the thermal distribution associated with the specified internal temperature (𝛽int ), while dynamical evolution is performed using the standard RPMD integration scheme with a timestep of 0.3. This choice of timestep is confirmed to avoid resonance instabilities, although we note that the best practice for future applications is to employ the Cayley-modification to the RPMD integration [135, 136]. Calculations with up to 144 ring-polymer beads are performed to ensure the convergence of the path-integral discretization. Implementation of the SPA utilizes Eqs. 4.17 and 4.18. Eq. 4.17 is first solved to obtain the stationary temperature from the thermal rate data. Then, the SPA microcanonical rate is obtained using Eq. 4.18. First- and second-order derivatives of 𝑘𝑄 𝑟 are obtained from a standard basis-spline interpolation procedure [137, 138]. Validation of the numerical procedure is performed by comparison with independent SPA results obtained from path-integral Monte Carlo sampling methods (Appendix A). Implementation of the MaxEnt approach closely follows the Bryan algorithm [124, 125]. Calculations are performed with a modified version of an open-source code.2 Quantum mechanical and RPMD thermal rates are the fitting targets in these calculations, while SPA microcanonical rates are employed as the Bayesian prior. For the exact quantum mechanical rate inversion, the diagonal elements of the covariance matrix in Eq. 4.14 are set to 0.4 ‰ of the rate value, as is consistent with the statistical noise in the RPMD thermal rate calculations. While not included here, MaxEnt results were also obtained using the flat Bayesian prior, 𝝀(𝐸) = 1; however, the numerics of these calculations were found to be less stable than those based on the SPA, which requires no additional information beyond the thermal rates that are also used for the fitting target. To ensure that the SPA priors are nonzero and sufficiently smooth, they are filtered with a low threshold value of 10−3 , followed by a simple moving average procedure to suppress the statistical fluctuations. A regulation potential (Eq. 4.15) with 𝜁 = 107 is sufficient to enforce the upper bound of 𝑁 (𝐸) to a tolerance of tight 10−4 in all reported calculations. 2 See https://github.com/jgreitemann/maxent for the available code kindly provided by Jonas Greitemann. 57 Figure 4.1: An illustrative example of the “L-curve” that is used to determine the parameter 𝛼 in each MaxEnt calculation. The optimal value of 𝛼 coincides with the kink in the curve (indicated by a red point). This example corresponds to the Eckart barrier, using RPMD thermal rate data for the fitting target and the SPA-RPMD microcanonical rates for the Bayesian prior. To specify the parameter 𝛼 in the objective function of the MaxEnt calculations (Eq. 4.11), the ‘L-curve’ rule was employed as in many previous studies [89, 126, 127]. As illustrated with a representative example in Fig. 4.1, the balance between fitting accuracy and solution likelihood when plotted as a parametric function of 𝛼 yields a hockeystick-shaped curve. We take the kink of the curve (red point) to correspond to the optimal balance between these attributes. 4.5 Results Microcanonical RPMD rates for electronically adiabatic reactions We begin by analyzing the effectiveness of the direct shooting approach with different choices of ring-polymer internal temperatures, 𝛽int . The symmetric Eckart barrier model for H+H2 reactive scattering [79] is chosen as the test example, with potential energy function . 𝑉 (𝑞) = 𝑉0 cosh2 (𝑞/𝑞 0 ) (4.23) 58 using parameters 𝑚 = 1061, 𝑉0 = 0.425 eV, and 𝑞 0 = 0.734. Analytical solution of the microcanonical rate for this model yields 𝑁 (𝐸) = 𝑓 /( 𝑓 + 𝑔) , where   √ 𝑓 = sinh2 𝜋𝑞 0 2𝑚𝐸/ℏ , and  √︃  2 𝑔 = cosh 𝜋 2𝑚𝑉0 𝑞 20 /ℏ2 − 1/4 , (4.24) and the exact thermal rate is obtained by integrating 𝑁 (𝐸) over Boltzmann kernel, following Eq. 4.7. Free-particle direct shooting Fig. 4.2A plots the microcanonical rate prediction from classical mechanics, exact quantum mechanics, and the direct shooting scheme (Eq. 4.22), as a function of energy. As expected, the step-function shape of the classical result is smoothed due to nuclear quantum effects. It is clear from the figure that the direct-shooting scheme is strongly sensitive to the choice of internal ring-polymer temperature, particularly at low energies for which tunneling plays an important role; irrespective of the employed value of 𝛽int , the direct-shooting scheme reverts to classical behavior in the high-energy regime. Strikingly, almost all nuclear quantum effects are absent using the free-particle protocol (𝛽int = 1/𝐸) for the internal ring-polymer temperature. Fig. 4.2B presents the canonical reaction rates for the Eckart barrier as a function of temperature, including exact quantum and classical results, as well as the standard RPMD calculation of thermal reaction rate (green) [78, 79]. As is well known for such problems, RPMD allows for the direct calculation of thermal reaction rates with good accuracy. it is exact at high temperatures [79] and works well at low energies due to its connection to semiclassical instanton theory [25]. However, the figure also shows the results of the RPMD thermal rate prediction obtained by transforming (via Eq. 4.7) the microcanonical RPMD rates from the free-particle direct-shooting protocol (magenta). Consistent with the results of Fig. 4.2A, this direct-shooting protocol provides an essentially classical description of the thermal reaction rate across the entire range of temperatures. Fig. 4.2B clearly demonstrates that approximation of microcanonical rates via the free-particle direct-shooting method (magenta vs. solid black) is a far greater source of error than the intrinsic approximation of RPMD for calculating thermal rates (green vs. solid black). This 59 Figure 4.2: Microcanonical reaction rate predictions for the Eckart barrier model from the direct shooting methods. (A) Microcanonical reaction rate predictions for the Eckart barrier model. Results are calculated with classical mechanics (classical MD, dashed black), analytical quantum mechanics (QM, solid black) and direct shooting approach with different internal ring-polymer temperatures (direct, solid cyan, orange and magenta). The free-particle direct shooting protocol is labeled 𝛽int = 1/𝐸. (B) Thermal rates obtained by substituting microcanonical rates from various levels of theory (classical MD, QM, and free-particle direct shooting) into Eq. 4.7. For comparison, the standard RPMD thermal rates (green) are also included. figure illustrates the hazards of using direct shooting for RPMD microcanonical rates, and it suggests that better results for microcanonical rate should be achievable on the basis of RPMD dynamics. 60 Stationary phase approximation We now turn our attention to the use of the SPA to calculate RPMD microcanonical reaction rates for the example of the Eckart barrier. Fig. 4.3A presents the calculated stationary temperature 𝛽st as a function of energy, obtained via Eq. 4.17 with input from either exact quantum thermal rates (blue) or standard RPMD thermal rate calculations (red, dashed). For comparison, we also show the temperature associated with the free-particle protocol (𝛽 = 1/𝐸), which differs significantly from the stationary temperature at each energy, as well as an analytical expression for the high-energy stationary temperature [𝛽 = 1/(𝐸 − 𝑉0 )] that is derived in Appendix A. In the low-energy regime, only the RPMD thermal rate data provides a satisfactory description of the stationary temperature obtained from the exact quantum results. Fig. 4.3B presents microcanonical rates from the SPA (Eq. 4.18) using input from RPMD thermal rates. To provide a baseline of accuracy associated with the SPA, we first compare the microcanonical rate from exact quantum mechanics (black, solid) with that obtained via SPA applied to exact quantum thermal rates (blue). The difference in these curves presents a best-case scenario for the accuracy of a method that approximates microcanonical rates via SPA; and it is seen that while the agreement at low energy is excellent, there is substantial deviation associated with energies in the high-energy regime. Encouragingly, essentially identical performance is seen when the microcanonical rates are obtained via application of the SPA to RPMD thermal rates (red, dashed). This indicates that the RPMD thermal rates are a smaller source of error than the SPA in the high-energy regime. Finally, comparison of the results in Fig. 4.3B with Fig. 4.2A makes clear that RPMD offers a much more accurate avenue to the calculation of microcanonical rates than might be concluded from simulations that employ direct shooting. We again note that the SPA errors at high energy in Fig. 4.3B that are well known and due to the neglect of higher-order terms in the phase function [131]. As anticipated form Eq. 4.21, both sets of SPA results in the figure converge to the erroneous high-energy asymptote of 1.084. Maximum entropy inversion As the third alternative for obtaining microcanonical rates from RPMD, Fig. 4.4A presents results for the Eckart barrier obtained using MaxEnt inversion. To establish the baseline error for the MaxEnt procedure, the dashed blue curve presents the 61 Figure 4.3: Stationary phase approximation (SPA) results for the Eckart barrier. (A) The stationary temperature calculated with QM (blue) and RPMD (dashed red) thermal rates, respectively, as a function of incident energy. For comparison, the classical asymptote of the stationary temperature (dashed black), and the freeparticle temperature (magenta) are also included. (B) Microcanonical rates obtained using the SPA with input from QM (blue) and RPMD (dashed red) thermal rates, respectively. For comparison, the microcanonical rates from exact QM and classical MD are also included. The inset expands the low-energy regime. 62 results obtained via inversion of the exact QM thermal rates using the SPA-QM microcanonical rates (Fig. 4.3B) as the Bayesian prior. Also, the dashed red curve presents the MaxEnt results obtained via inversion of the RPMD thermal rates using the SPA-RPMD microcanonical rates as the Bayesian prior. This last result utilizes input from RPMD thermal rates alone. It is clear from Fig. 4.4A that the MaxEnt procedure provides excellent accuracy across the entire range of energies on the scale plotted, avoiding the incorrect highenergy asymptote of the SPA results. Closer examination of the low-energy regime in the inset shows that the agreement persists even in the regime of strong tunneling. Comparison in Fig. 4.4B of the SPA results (SPA-QM and SPA-RPMD) with their MaxEnt refinements reveals that MaxEnt can introduce significant relative errors in the rate predictions at low temperatures for which the rates are very small. Taken together, the results in Fig. 4.4A and B indicate that for this example, the use of MaxEnt inversion helps to improve the quality of the SPA at intermediate and higher energies, but it does little to improve the quality of the SPA in the low-energy regime. Finally, as a self-consistency check, Fig. 4.4C presents the thermal rates obtained by transforming (via Eq. 4.7) the microcanonical rates obtained from the MaxEnt inversion of the RPMD thermal rates (red, dashed). For comparison, the exact quantum, classical, and standard RPMD thermal rates are also included. As expected, the MaxEnt RPMD rates are fully consistent with the standard RPMD thermal rates, and both are in good agreement with the exact QM results. Stationary-temperature direct shooting Given the success of the SPA for extracting microcanonical rates from standard RPMD thermal rates, it is tempting to see whether data obtained from the SPA can be used to improve the direct shooting method. Specifically, we explore the use of the stationary temperature as the ring-polymer internal temperature for initializing and propagating the direct-shooting trajectories, i.e., setting 𝛽int = 𝛽st in Eq. 4.22. This strategy is physically appealing, since the stationary temperature (which is a function of incident energy, see Fig. 4.3A) dictates the delocalizaton of the ringpolymer in its barrier-crossing configuration [25, 130, 132]; also note that the stationary temperature approaches the free-particle temperature at high incident energies. 63 Figure 4.4: Maximum entropy (MaxEnt) inversion results for the Eckart barrier. Methods are labeled with the format “MaxEnt-[thermal rate input type]/[prior type].” (A) MaxEnt solutions for the microcanonical reaction rate as a function of incident energy, with the inset expanding the low-energy regime. Microcanonical rates from classical MD and exact QM are also presented for reference. (B) An expanded view of panel A in the low-energy regime, plotted on a log-log scale. (C) Thermal rates obtained by integrating 𝑁 (𝐸) in panel A over the Boltzmann kernel. For comparison, exact QM, classical MD, and standard RPMD thermal rates are also included. 64 Figure 4.5: Revisiting the microcanonical rate predictions for the Eckart barrier from the direct shooting methods. The plot compares the direct-shooting method with the ring-polymer internal temperature set to either the stationary temperature (red) or the free-particle temperature (magenta). Also included are the exact QM and classical MD results. Fig. 4.5 plots the microcanonical rate for the Eckart barrier, obtained using the stationary-temperature direct-shooting method (red). For comparison, the freeparticle direct-shooting (magenta), classical MD (black, dashed), and exact quantum (black, solid) results are all reproduced from Fig. 4.2A. While stationary-temperature direct shooting remains qualitatively less accurate than the SPA and MaxEnt inversion methods, it nonetheless substantially improves the results of the free-particle direct-shooting approach in the low-energy region where tunneling is important. These results indicate that stationary-temperature direct shooting is a less quantitative tool than SPA or MaxEnt for the calculation of microcanonical rates from RPMD trajectories, but it may nonetheless prove useful in exploratory studies for which a direct trajectory-based simulation approach is needed, or in applications to high-dimensional systems for which obtaining precise thermal reaction rate in the whole temperature region is computationally expensive. Microcanonical RPMD rates for non-adiabatic systems Although we have thus far only discussed the SPA and MaxEnt inversion methods in the context of single-level (i.e., electronically adiabatic) processes, both methods can be naturally extended to multi-level systems. Specifically, given state-resolved 65 thermal reaction rates for a non-adiabatic process, both the SPA and MaxEnt methods can be used to compute state-resolved microcanonical rates for different reaction channels. For the SPA method, state-resolved thermal reaction rates are substituted into Eqs. 4.17 - 4.18, respectively, yielding a single stationary temperature and a single stateresolved microcanonical rate for each reaction channel. For the MaxEnt method, we solve the coupled integral equation ! ! ! 𝜿1→2 B 𝝂1→2 = · (4.25) 𝜿2→2 B 𝝂2→2 for a system with two reaction channels (e.g., diabat 1 → 2, and diabat 2 → 2), with 𝜿, B and 𝝂 defined in Eq. 4.10. The MaxEnt objective function for the two-level system is 𝑄(𝝂1→2 , 𝝂2→2 ; 𝛼) = 𝛼𝑆(𝝂1→2 ) − 𝜒2 (𝝂1→2 )/2 + 𝛼𝑆(𝝂2→2 ) − 𝜒2 (𝝂2→2 )/2 + 𝑉reg (𝝂1→2 , 𝝂2→2 ) (4.26) which sums the information entropy and likelihood function contributions for the state-resolved rates. The regularization potential is likewise generalized, ∑︁ 1 𝜁 𝐼 2 (𝜈1→2, 𝑗 + 𝜈2→2, 𝑗 − 1), (4.27) 𝑉reg (𝝂) = − 2 𝑗 to enforce unitarity (with the detailed balance condition) 𝑁1→2 (𝐸 𝑗 ) + 𝑁2→2 (𝐸 𝑗 ) ≤ 1 ∀ 𝑗 . (4.28) The constraint of non-negativity 𝑁1→2 (𝐸 𝑗 ) ≥ 0, 𝑁2→2 (𝐸 𝑗 ) ≥ 0, ∀ 𝑗, (4.29) is enforced as before by confining the solution search to the positive subspace. As a numerical demonstration for non-adiabatic reaction dynamics, we use a twostate gas-phase reactive scattering model that has been previously introduced [93, 139]. In the diabatic representation, the potential energy functions for this system are 𝐴1 𝑉11 (𝑞) = + 𝐵1 1 + 𝑒 −𝑎1 (𝑞−𝑞1 ) 𝐴2 𝐵2 𝑉22 (𝑞) = + −𝑎 𝑞 2 2 1+𝑒 4 cosh (𝑎 2 𝑞/2) 𝑉12 (𝑞) = 𝑉21 (𝑞) = 𝐴3 𝑒 −𝑎3 (𝑞−𝑞3 ) 2 (4.30) 66 Figure 4.6: State-resolved microcanonical rates for the two-level model in Eq. 4.30. Dashed lines indicate the 1→ 2 diabatic reaction channel and solid lines indicate the 2→ 2 diabatic reaction channel. (A) Microcanonical rates from classical surface hopping (orange) and numerically exact QM wavepacket propagation (black). (B) SPA results for the microcanonical rate, with input from exact QM (blue) and iso-RPSH (red) thermal rates. (C) MaxEnt results for the microcanonical rate, with input from exact QM (blue) and iso-RPSH (red) thermal rates. 67 with parameters specified in Table 4.1 and with reactant and product regions corresponding to 𝑞 → −∞ and 𝑞 → ∞, respectively. Table 4.1: Parameters for the two-level model in Eq. 4.30. Parameter Value 𝐴1 7 𝐴2 −18/𝜋 𝐴3 0.25 𝐵1 −0.75 𝐵2 54/𝜋 Parameter Value 𝑎1 √ 1 3𝜋/4 𝑎2 𝑎3 0.25 𝑞1 −1.6 𝑞3 −2.625 We focus on microcanonical rates in the range of incident energies for which the higher-energy state is unavailable as a product channel (i.e., the only two available reactions channels correspond to the 1 → 2 and 2 → 2 processes on the diabatic states). Calculation of state-resolved thermal rates for the two channels is performed with the flux-side formulation [139] of iso-RPSH [93], with both the methodological details and thermal-rate results both reported elsewhere [139]. The only difference in the current work is that 1000-fold more trajectories are performed to suppress statistical error in the thermal rates for the MaxEnt calculations, although the thermal rate results are graphically indistinguishable from those previously published [139]. For comparison, Fig. 4.6A presents the state-resolved microcanonical rates obtained from numerically exact quantum mechanics [73, 100] and using classical surface hopping [13] as implemented in our previous work [93]. As for the one-level system, the microcanonical rates with classical nuclear dynamics exhibit a sharp increase when the incident energy reaches the barrier height. Although classical surface hopping qualitatively includes the effect of the non-adiabatic transition and performs well in the high-energy regime, it fails to capture the significant nuclear quantum effects in this problem. Fig. 4.6B presents the microcanonical rates obtained via application of the SPA to state-resolved thermal rates from exact QM (blue) and from iso-RPSH (red). Comparison of the QM and SPA-QM results indicate that the SPA approximation is a good approximation in this example. Furthermore, comparison of these curves with the SPA-iso-RPSH results indicates that the iso-RPSH thermal rate data is an even smaller source of error than the SPA. At higher energies, the SPA results exhibit the same pathologies as those discussed in connection with Eq. 4.21, and the better performance of SPA-iso-RPSH versus SPA-QM in this regime is likely due to fortuitous error cancellation. 68 Fig. 4.6C presents the microcanonical rates obtained via application of the MaxEnt method to the state-resolved thermal rates. As for the one-level system described in Fig. 4.4, the MaxEnt method improves the description for the two-level system at high energies but does little to refine the description of the SPA at lower energies. Taken together, these results indicate that the iso-RPSH method can be staightforwardly extended for the accurate calculation of state-resolved microcanonical rates. Moreover, these results show that iso-RPSH provides an accurate description of both the thermal and microcanonical reaction rates of this system in a regime for which both non-adiabatic and nuclear quantum effects play an important role, although the method has been shown to underestimate the asymmetry in the Marcus inverted regime in the golden-rule limit of electron transfer [140]. 4.6 Summary Whereas the ring-polymer molecular dynamics (RPMD) thermal rate theory has proven immensely successful in many chemical application domains, far less attention has been paid to the problem of calculating microcanonical reaction rates using RPMD, which may be of considerable value in the context of gas-phase and surface-molecule scattering processes for which nuclear quantization plays an important role [141–143]. The current work addresses this shortcoming by exploring a variety of strategies to calculate microcanonical reaction rates with RPMD. It is found that the ad hoc strategy of direct shooting of RPMD trajectories is strongly sensitive to the internal ring-polymer temperature that is employed; this is somewhat ameliorated in the tunneling regime via the use of an internal temperature based on the stationary-phase approximation (SPA), but the resulting direct-shooting results remain overly classical in the barrier-crossing energy regime. Far more accurate microcanonical rates are obtained from RPMD thermal rate data via Laplace transform inversion using either the SPA or the numerically exact maximum entropy method. In general, we find that the SPA applied to RPMD thermal rate data provides the best compromise between good accuracy and numerical feasibility, particularly in the low-energy tunneling regime, although we point out that the alternative directshooting and maximum entropy methods described here may also prove useful in particular application cases. While the current paper focuses only on the calculation of microcanonical rates from RPMD (or other approximate) thermal rate data, we note that the same strategy can be applied for the calculation of collisional cross section, and other time correlation functions, spectra, and transport coefficients in the microcanonical ensemble. 69 4.7 Appendix Evaluation of the stationary temperature via path-integral Monte Carlo sampling methods As an alternative method for evaluating the stationary temperature 𝛽st in Eq. 4.17, we apply path-integral Monte Carlo sampling methods in combination with quantum transition state theory (QTST) approximation [98] to the reaction rate [144–146]   𝑘 (𝛽) ≃ 𝑘 QTST (𝛽, 𝑞 ‡o ) = min 𝑘 QTST (𝛽, 𝑞 ‡ ) , (4.31) 𝑞‡ with a dividing surface 𝑞 ‡o that minimizes dynamical recrossing effects. Using the path-integral representation of the QTST rate in the case of a single-surface system  𝑛  𝑛∫ iso ‡ 𝑑p0 𝑑q0 𝑒 −𝛽𝐻𝑛 (p,q) 𝑘 QTST (𝛽, 𝑞 o )𝑄 𝑟 = lim 𝑛→∞ 2𝜋ℏ × 𝛿( 𝑞¯ 0 − 𝑞 ‡o ) 𝑣¯ 0 ℎ( 𝑣¯ 0 ), (4.32) together with Euler’s theorem for homogeneous functions [147], we derive a viriallike expression for the stationary energy-temperature relation, which can be conveniently evaluated using path-integral Monte Carlo sampling methods + * 𝜕𝑉 (𝑞 𝛼 ) 1 ∑︁ 1 1 ∑︁ (𝑞 𝛼 − 𝑞¯ 0 ) 𝑉 (𝑞 𝛼 ) . (4.33) 𝐸 st = + + 𝛽st 2𝑛 𝛼 𝜕𝑞 𝛼 𝑛 𝛼 ‡ 𝑞o Here, ⟨Q⟩𝑞‡ is a constrained ensemble average defined by o ∫ 𝑑q0 𝜌𝑐 (q0 ) Q (q0 ) ∫ ⟨Q⟩𝑞‡ = , o 𝑑q0 𝜌𝑐 (q0 ) (4.34) with 𝜌𝑐 (q0 ) = 𝛿( 𝑞¯ 0 − 𝑞 ‡o ) 𝑒 −𝛽𝑈spr (q0 ) 𝑒 −𝛽𝑛 Í 𝛼 𝑉 (𝑞 𝛼 ) . (4.35) This approximate method for obtaining 𝛽st works well for the numerical examples considered here, although it is more limited than the general numerical differentiation method that is described in the Computational Details section. Now let us proceed by analyzing the high-temperature asymptotic behavior of the stationary temperature. At high temperatures, Eq. 4.33 approaches the classical limit for the stationary energy-temperature relation lim 𝐸 st = 𝛽st →0 1 + 𝑉 (𝑞 ‡o ), 𝛽st (4.36) 70 with 𝑞 ‡o approaching the barrier top. This classical limit can be also obtained upon the substitution of the classical transition state theory rate, [𝑘𝑄 𝑟 ] CTST = 1 −𝛽𝑉 (𝑞‡o ) , 𝑒 2𝜋𝛽ℏ (4.37) into Eq. (4.17), yielding   𝜕 log 𝑘 CTST (𝛽)𝑄 𝑟 (𝛽) 𝐸 st = − 𝜕𝛽 = 𝛽st 1 + 𝑉 (𝑞 ‡o ). 𝛽st (4.38) 71 Chapter 5 SMALL NUCLEAR QUANTUM EFFECTS IN SCATTERING OF H AND D FROM GRAPHENE This chapter is based on the following publication: 1. Jiang*, H., Tao*, X., Kammler, M., Ding, F., Wodtke, A. M., Kandratsenka, A., Miller III, T. F. & Bünermann, O. Small Nuclear Quantum Effects in Scattering of H and D from Graphene. Journal of Physical Chemistry Letters 12, 1991–1996. doi:10.1021/acs.jpclett.0c02933 (2021). 5.1 Abstract We study nuclear quantum effects in H/D sticking to graphene, comparing scattering experiments at near-zero coverage with classical, quantized, and transition-state calculations. Experiment shows H/D sticking probabilities that are indistinguishable from one another and markedly smaller than those expected from a consideration of zero-point energy shifts of the chemisorption transition state. Inclusion of dynamical effects and vibrational anharmonicity via ring-polymer molecular dynamics (RPMD) yields results that are in good agreement with experiment. RPMD also reveals that nuclear quantum effects, while modest, arises primarily from carbon and not from H/D motion, confirming the importance of a C-atom re-hybridization mechanism associated with H/D sticking on graphene. 5.2 Introduction Nuclear quantum effects (NQEs) are ubiquitous in chemistry, typically resulting from quantum resonances, zero-point energy or tunneling. H/D substitution forms the basis of nearly all experiments designed to reveal NQEs, as the resulting energy shifts to the quantum levels influence differential reactive scattering cross-sections [148], alter rate constants or even shift equilibria in chemical reactions [149]. However, NQEs in chemistry are far from being fully understood, particularly in the context of interfacial chemical reactivity for which the classical phonon approximation is by necessity widely made but rarely tested. Within this context, graphene is an ideal test system to examine NQEs in surface chemistry. The combination of light C-atoms and stiff C-C bonds drives graphene’s phonon spectrum to high frequency, where quantum effects might be important. 72 The simplest reaction on graphene is adsorption of an H or a D atom. A physisorption well with a depth of ∼40 meV is found at C-H distances of ∼4 Å [150, 151] and at a C-H distance of ∼1.1 Å, there is an ∼800 meV deep chemisorption well. Since chemisorption involves sp2 to sp3 re-hybridization of a C atom, there is a pronounced barrier between the physisorbed and chemisorbed states [152, 153]. Saturated H/D atom adsorption shows an inverse isotope effect, possibly due to the interplay between the adsorption, reflection and associative desorption [154]. Hydrogen ions also penetrate through a graphene sheet faster than deuterium ions [155]. Additionally, previous computational studies [156, 157] have predicted evident nuclear quantum effects to hydrogen sticking on graphene at low temperature, a prediction that awaits experimental confirmation. Recently, H atom scattering experiments were combined with first-principles theory to show that H collisions at graphene induce concerted in-plane motion of the carbon atom framework and extraordinarily fast energy dissipation leading to chemisorption [106]. The experiments produced nearly mono-energetic H atoms and could be performed at near-zero coverage, removing well-known ambiguities associated with the energy and coverage dependence of C-H bond formation on graphene [158, 159]. In this paper, we present new scattering experiments comparing H and D interactions at graphene, and we analyze the role of NQEs using both transition-state theory, classical molecular dynamics, and a quantized molecular dynamics approach. We employ a tight-binding potential energy surface (PES) reparametrized to match electronic structure data obtained with a hybrid exchange-correlation functional, while maintaining a high level of numerical efficiency. This represents a significant improvement over our previously reported PES [106], correctly reproducing the graphene phonon density of states spectrum (PDOS), a quantity that is critically important to testing the classical phonon approximation. As we will show below, experiment reveals little or no evidence of an H/D isotope effect. This contradicts the prediction of a microcanonical transition-state theory treatment based on the Rice–Rampsberger–Kassel–Marcus (RRKM) rate theory (SI Sec. B), which exhibits a five fold enhancement of stickng for D versus H at the experimentally investigated normal incidence energies. While RRKM is obviously not suited to simulating a scattering experiment, the RRKM predicted H/D isotope effect arises from differences in the zero-point energy of the isotopes at the transition-state, and it might reasonably be expected that these effects are also important in non-equilibrium scattering experiments. To test this idea, ring-polymer 73 molecular dynamics (RPMD) simulations are performed, which, in agreement with experiment, show that the longer interaction time of D with the graphene flake compared to H leads to only an extremely small increased sticking probability for D, a result of classical inertia. Within the RPMD framework, we perform mixed quantum/classical calculations where only selected degrees of freedom are treated quantum mechanically. This approach shows that the largest NQEs in this system are associated with the C-atom motion of the graphene—these NQEs are not present when using a PES that fails to reproduce the high frequency region of graphene’s PDOS spectrum, those that are directly involved in the chemisorption reaction. These results show that the classical approximation for solid phonons leads to errors in the description of the dynamics of transient C-H bond formation in H collisions with graphene, providing guidance concerning the validity of this assumption in other systems. 5.3 Experiment The experimental setup has been described in detail in Reference [160–162]. All D experimental measurements are newly reported in this work. H/D atoms are produced via photodissociation of a supersonic molecular beam of hydrogen/deuterium iodide with a 10 ns pulsed KrF laser producing atoms with incidence energy E𝐼 ∼ 1.0 eV. A portion of the atoms go through a skimmer and two differential pumping stages, enter the ultra-high vacuum chamber and scatter from the graphene sample. The sample is held on a 6-axis manipulator, allowing variation of the incidence angle 𝜗𝑖 . About 0.7 mm above the surface, the scattered H or D atoms are excited to a long-lived Rydberg state (𝑛 = 34) by two spatially and temporally overlapped laser pulses at 121.57 nm and 365.90 nm via a two step excitation. The neutral Rydberg atoms travel 250 mm before they are field-ionized and detected by a multi-channel plate detector. The arrival time is recorded by a multi-channel scalar. The rotatable detector allows data to be recorded at various scattering angles 𝜗𝑠 . The graphene sample is epitaxially grown in situ on a clean Pt(111) substrate by dosing ethylene (partial pressure 3 × 10−8 mbar) at 700 ◦ C for 15 mins [106]. 5.4 Computational Methods The all-atom potential energy surface (PES) employed in this study is obtained from the “Geometry, Frequency, Noncovalent, eXtended Tight Binding” (GFNxTB) method [163]. While an accurate C-H adsorption barrier and binding well can be obtained [164] with a hybrid density functional theory (DFT) [165, 166], the 74 associated computational costs make it impossible to simulate scattering energy and angular distributions. We found however, that a reparameterization of GFN-xTB within the entos software package [167] to best reproduce the minimum energy path of a hybrid functional DFT calculation resulted in an accurate all-atom PES at low computational cost. Details are presented in the SI Sec. A. Notably, we are able to compute the GFN-xTB PES 1000-fold more efficiently than a similar PES using a hybrid DFT functional. We are thus able to run far more trajectories than would otherwise be possible. We use ring-polymer molecular dynamics (RPMD) under the influence of GFNxTB potential energies to model the quantized molecular dynamics of the system. RPMD [20, 21] is a trajectory-based dynamics method, in which the NQEs are taken into account based on Feynman’s imaginary-time path-integral formalism [18]. Although approximate, the method successfully describes both zero-point energy and tunneling effects in simulations at thermal equilibrium [21, 80, 168] and more recently, it was applied to systems with non-equilibrium initial conditions [118]—for instance, the dissociative molecular sticking process [107]—and in the microcanonical ensemble [169]. It is worthwhile mentioning that RPMD reproduces exact quantum dynamics at the short-time and within the harmonic limit [20, 21]. This makes it well suited to the study of ultrafast scattering processes involving surface phonons that are largely harmonic, like H sticking at graphene. RPMD allows for the inclusion of NQEs by propagating classical trajectories of an isomorphic system. The isomorphic system consists of 𝑛 replicas of the physical one, and is constructed such that exact quantum Boltzmann statistics are preserved [18, 19, 55]. If 𝒒 = (𝒓 H1 , 𝒓 C1 , 𝒓 C2 , ...) denotes the column vector of positions of all atoms and 𝑉 (𝒒) the potential energy for a geometry 𝒒, then the RPMD equations of motion (EOM) are 𝒒¥ 𝛼 = 𝜔2𝑛 (𝒒 𝛼−1 + 𝒒 𝛼+1 − 2𝒒 𝛼 ) − 𝒎 −1 · ∇𝒒 𝛼 𝑉 (𝒒 𝛼 ), (5.1) where 𝛼 = 1, 2, ..., 𝑛 is the index for different replicas, 𝜔𝑛 = 𝑛𝑘 B𝑇 is the strength of the harmonic springs that connect neighbouring replicas with 𝑇, being the system temperature, and 𝒎 = diag (𝑚 H1 , 𝑚 C1 , 𝑚 C2 , ...) the mass matrix for all the atoms involved. Note that with Eq. 5.1, all the atoms in the system are described quantum mechanically on the same footing. 75 Not only does RPMD provide an accurate and efficient way to perform quantum simulations, the fact that it captures NQEs via trajectory propagation in classical phase space bridges the gap between classical and quantum mechanics. RPMD recovers exact quantum statistics in the limit 𝑛 → ∞ [18] and reduces to ordinary classical molecular dynamics when 𝑛 = 1. Hence, RPMD can be used for mixed quantum/classical (MQC) calculations where some degrees of freedom are described quantum mechanically and others classically [45, 87, 89, 116, 170]. For a system partitioned into quantum 𝒒 and classical 𝑸 parts, MQC-RPMD evolves the dynamics for the quantum mechanical portion using RPMD EOM (for 𝑛 replicas of the original system), and evolves the dynamics for the classical portion with Newton’s EOM [89], i.e. 𝒒¥ 𝛼 = 𝜔2𝑛 (𝒒 𝛼−1 + 𝒒 𝛼+1 − 2𝒒 𝛼 ) − 𝒎 −1 ∇𝒒 𝛼 𝑉 (𝒒 𝛼 , 𝑸), 𝑛 1 ∑︁ −1 ¥ ¯ 𝑸), 𝑸 = −𝑴 ∇𝑸 𝑉 ( 𝒒, 𝒒¯ = 𝒒𝛼 𝑛 𝛼=1 (5.2) where 𝒎 and 𝑴 are the mass matrices for the quantum and classical degrees of freedom, respectively. As in Reference [106], we perform quantized molecular dynamics simulations by a non-equilibrium RPMD approach [118]. The graphene surface is modelled with a free-standing cluster of 42 carbon atoms. This model is sufficiently large to describe the PDOS spectrum of graphene obtained from calculations using periodic boundary conditions—see Fig. 5.9. The boundary of the carbon cluster is terminated with H atoms, which are held fixed throughout the calculation using a RATTLE scheme [171]. A suspended graphene flake with 80 C-atoms is used to eliminate the edge effects under conditions with larger values of 𝜗𝑖 . RPMD simulations are initialized by separately preparing the initial configurations for the non-interacting graphene sheet and H/D atom. The initial flake geometries are sampled from a thermalized ring-polymer trajectory at 300 K, performed using the Andersen thermostat [172]. The position and velocity of the centroid of the H/D atom ring-polymer is then determined according to the values of scattering energy and incidence angle in the experiment. The internal modes for the H/D-atom ring-polymer are thermalized at the surface temperature 300 K following the nonequilibrium RPMD formulation [106, 118], and it is confirmed that the results are insensitive to this choice of internal temperature—see SI Fig. 5.11. The temporal evolution of the system is then found from either Eq. 5.1 or 5.2 using a time step 76 of 0.5 fs [136]. Propagation continues until the fate of the H or D atom is decided by trajectory analysis. For scattered atoms, the energy loss and outgoing angle is recorded for each trajectory. Convergence of the path-integral discretization was confirmed to be reached with 8 beads. All reported simulations were performed using the entos software package [167] and the trajectory visualization is rendered with the OVITO software [173]. Further details of the trajectory calculations are reported in the SI Sec. C. 5.5 Results Fig. 5.1(a) and (b) show examples of experimental scattering distributions for H and D, respectively, colliding with graphene. Both scattering distributions peak close to the specular angle and energy loss is small but slightly larger for D than for H. Additional scattering distributions for other values of 𝜗𝑖 are shown in Fig. 5.6. The scattered flux comes from atoms reflected at the barrier to chemisorption; it therefore decreases as 𝜗𝑖 is reduced, because ever more projectile atoms pass over the barrier to chemisorption with higher normal energy of incidence [106]. Fig. 5.1(c) and (d) show RPMD simulations of the experimental results of Fig. 5.1 (a) and (b); see also Fig. 5.6. Agreement is excellent. These trajectory calculations also show that the scattered flux is due to H/D atoms that never reach the chemisorption well; those that do remain trapped. Thus, the observed signals are a direct measure of the H/D survival probability. Experimentally derived sticking probabilities are shown in Fig. 5.2(a). See SI Sec. D and Fig. 5.7. The sticking probability increases with the normal component of the incidence energy. Thresholds near ∼0.35 eV are seen for both H and D, a clear sign of the adsorption barrier. No H/D isotope effect can be discerned from the experiment. Fig. 5.2(a) also shows the sticking probabilities from RPMD simulations, which are in good agreement with experiment. As the focus of the article is on the isotope and nuclear quantum effects, a horizontal shift (for both H and D) is introduced before comparing to experiment in Fig. 5.2(a), in order to account for the influence of the Pt substrate, which is neglected in the theory. We note the shift is introduced based on our published analysis and the value of the shift is consistent with that presented previously [106]. A small inverse H/D isotope effect is clearly present—D sticking is ∼10% more likely than H sticking. 77 (a) (b) H (c) D (d) H D Figure 5.1: Comparing theory with experiment for H and D scattering distributions from graphene. (a-b) Experimental distributions. (c-d) Theoretical distributions. Result shown are from RPMD trajectories under the influence of GFN-xTB potential energies. The incidence energy of H/D translation, E𝐼 , was ∼ 1 eV. In all images, the normalized scattering energy, E𝑆 /E 𝐼 , is shown along the radial coordinate and the scattering angle, 𝜗𝑠 , is shown on the polar coordinate. The distributions are normalized to their integrals. The red ticks indicate the specular scattering angle. 5.6 Discussion By comparing the predictions of sticking probabilities from different computational approaches, we gain insight into how classical and quantum effects manifest in the dynamics of H and D sticking to graphene. Fig. 5.2(b) compares the sticking probability results from the classical molecular dynamics (cMD) and RPMD, for H and D atom scattering. cMD simulations predict an increased sticking for D compared to H (Fig. 5.2(b): solid red vs. solid black lines). This results from the ∼1.4x longer interaction time of the D atom compared to the H atom—see Fig. 5.10—a result of classical inertia [174, 175]. The longer interaction time of D allows greater relaxation of the C flake during the trajectory, reducing the height of the effective barrier to sticking. 78 (b) 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 sticking probability (a) 0.0 0.0 / / exp., H / D RPMD, H / D 0.2 0.4 0.6 0.8 normal incidence energy / eV 0.0 0.0 / / RPMD, H / D cMD, H / D 0.2 0.4 0.6 0.8 normal incidence energy / eV Figure 5.2: Comparison of experiment and theory for H/D sticking probabilities to graphene. (a) Experimental sticking probabilities (filled circles with error bars) for H (black) and D (red) are compared to ring polymer molecular dynamics (RPMD, dashed lines and symbols) simulations. (b) Sticking probabilities for H and D from classical (solid lines and symbols) simulations are compared to those from quantized RPMD simulations. The incidence energy was held constant and the incidence angle was varied to control the normal component of incidence energy. The statistical error associated with the trajectory calculations is less than 1.1%. The simulation curves are shifted on the horizontal axis by -0.12 eV to compare with the experiment results. We argue next that Fig. 5.2(b) shows that the NQEs arising from the C-atom motion influence the reaction, enhancing both H and D sticking. First note that NQEs are present for both H and D scattering—the deviation of RPMD from cMD predictions for both isotopes make this clear. To investigate the role of C atoms, we calculated sticking probabilities with the MQC implementation of RPMD [45, 87, 89, 116, 170]—see Fig. 5.3(a). In the MQC-RPMD simulations, either the projectile atom or the graphene flake was described quantum mechanically, while the remainder of the system was described classically. When the projectile is quantized but the graphene flake moves classically (dashed green line), no difference is found in comparison to the conventional cMD (black solid line). However, when the C atom motion is quantized but the projectile atom is treated classically (dashed blue line), no difference is seen compared to full RPMD (dashed black line). We emphasize that the NQEs associated with the C atoms, although modest, are much larger than the statistical uncertainty of the ensemble averages. Fig. 5.3(b) presents the phonon density of states (PDOS) spectrum of a free-standing graphene flake obtained with the GFN-xTB potential energy surface. Calculation details are discussed in SI Sec. E. The PDOS of graphene consists of two major branches, i.e., the out-of-plane vibrations, which are softer in frequencies and are 79 (a) (b) Figure 5.3: Nuclear quantum effects in the H/D-on-graphene sticking. (a). Simulated H-on-graphene sticking probability as a function of normal incidence energy. Results are obtained with classical MD, all-atom quantized RPMD (with Eq. 5.1), and the mixed quantum/classical implementations of RPMD (with Eq. 5.2). (b). Out-of-plane (blue) and total (red) phonon density of states spectrum of the graphene flake from molecular dynamics simulations. Previously published simulation result [176] (filled gray) is also included for reference. relevant to C puckering out (blue branch); and the in-plane vibrations, which are at higher frequencies (the optical branch from 1200–1600 cm−1 ). The NQEs from the C-atom motions that encourage H/D sticking are from the high frequency group of vibrational modes, because these NQEs entirely disappear when we use the same approach with a PES that fails to capture the distribution of those higher vibrational frequencies of the graphene flake (see Figs. S5 and S6). This comparison emphasizes the role that the in-plane C motions play in driving the reaction center C atom to pucker out of the graphene plane and capture the scattering projectile. As a result, when those in-plane optical phonons are modelled with quantized instead of classical simulations, the NQEs from those neighbouring C atoms make it easier for the system to cross the chemisorption barrier and therefore increase the sticking probability. We also note that the importance of those NQEs originating from C atoms reveals a failure of the classical approximation for describing processes that are sensitive to high frequency phonons. While RPMD provides only an approximate description of quantum dynamics, there are reasons to expect that it is well suited to H scattering from graphene. First, for both equilibrium and non-equilibrium processes, RPMD is exact in the limit of short-time, where the relevant timescale is 𝛽ℏ ∼ 25 fs at the temperatures used in these simulations. This is substantially longer than the classical timescale of 80 H/D-surface interaction, which is only 10 fs (Fig. 4A of Reference [106]). Second, RPMD is exact in the classical mechanical limit, giving us confidence that it can accurately describe the small NQEs seen here. Finally, although the exactness of RPMD for both the short-time and classical limits holds regardless of whether the simulated property corresponds to a correlation function of linear or non-linear operators, we note that the RPMD approximation is typically more accurate for time correlation functions (TCFs) of linear operators than for non-linear operators; [21] while this might raise concern because the fluctuation-dissipation description of vibrational energy relaxation involves a TCF of non-linear operators [177], the centroid position of H/D atom that dictates whether sticking occurs is in fact a linear function of position. All in all, these considerations suggest that RPMD provides a reliable description of the physical processes considered here, while allowing for the full-dimensional simulations that are needed to capture key aspects of the sticking mechanism [106]. Regardless of these methodological considerations, we emphasize that the mechanistic picture that emerges is physically intuitive and consistent with our previous interpretation [106]. Although perhaps at first surprising, the results are easily understood when considering the dynamical mechanism for sticking. Ultrafast energy loss of H atom translation is associated with in-plane C-atom excitation adjacent to the reaction center. When the in-plane C-atom motion important to the sticking is more realistically modelled using quantum simulations rather than classical simulations, sticking is influenced. We began this study with the intention to investigate H/D isotope effects in adsorption at graphene; however, it turns out that the H/D isotope effects are negligible. By contrast, our work shows that the quantum motion of C-atoms can have an influence on H/D sticking at graphene. This reflects the fact that the high frequency C-C vibrations—those that are most quantum mechanical—are diretly involved in the formation of a C-H chemical bond during the adsorption process. 5.7 Supplementary Materials A. Modified GFN-xTB Hamiltonian For the purpose of studying the scattering dynamics of H and D from graphene, it was necessary to modify the GFN-xTB Hamiltonian [163], originally developed for near-equilibrium geometries. Following Eq. 16 in Reference [163], the potential 81 energy of the system is 𝐸 GFN-xTB-mod = 𝐸 GFN-xTB + 𝐸 RE (5.3) with 𝐸 RE = atoms ∑︁ CNA 𝑓A . (5.4) A 𝑓A is the element-specific atomization correction energy, and CNA represents the coordination number for atom A [178], defined as follows. atoms  −1 ∑︁  CNA = 1 + 𝑒 −𝜆[𝜇(𝑅A +𝑅B )/𝑟AB −1] (5.5) B≠A 𝑅A/B and 𝑟 AB are the covalent radius of atoms A/B and the distance between atoms A and B, respectively. We fit 𝑓 , 𝜆 and 𝜇 to benchmark density functional theory (DFT) calculations at the B3LYP/cc-pVDZ/JKFIT [165, 166, 179, 180] level, with reference configurations consisting of an H atom interacting with a C42 H16 graphene cluster. Values for the optimized parameters are reported in Table 5.1. Fig. 5.4(a) shows the minimum energy path (MEP) for an H atom’s perpendicular approach to a carbon atom on the graphene flake, comparing the B3LYP and optimized GFN-xTB-mod energies. Fig. 5.4(b) shows the two-dimensional cut through the high dimensional potential energy landscape, with all other degrees of freedom are allowed to relax at each C-H distance. The ability of GFN-xTB-mod to reproduce the B3LYP barrier and well-depth is remarkably good. Furthermore, its use accelerates the dynamics calculations of this work by one-thousand fold, compared to calculations with a hybrid functional. This allows the GFN-xTB-mod approach to be further tested by simulations of experimental scattering distributions, where good agreement is also seen. See the main text. B. Transition state theory estimation of the isotope effect To investigate the isotope effect in the H/D-on-graphene reactive system, we first employ the Rice–Rampsberger–Kassel–Marcus (RRKM) rate theory [181]—a standard microcanonical transition state theory (mTST)—to compare the difference Table 5.1: Results for the reactive energy term that describes C-H interaction, as defined in Eqs. 5.4 and 5.5. Parameters Value Parameters 𝑓H 0.6448 eV 𝜆 𝑓C 0.0261 eV 𝜇 Value 7.141 1.144 82 Figure 5.4: The minimum energy pathway (MEP) for the perpedicular approach of an H atom on top of a carbon atom on the graphene flake. (a) Comparing the one-dimensional surface at the GFN-xTB-mod and B3LYP/cc-pVDZ levels of theory as a function of the C-H distance, with all other degrees of freedom optimized at each C-H distance. (b) Using GFN-xTB-mod, the two-dimensional surface as a function of the C-H distance and the displacement of the interacting C atom from the plane of the graphene flake, again with all other degrees of freedom optimized. All energies are in eV. between the chemisorption rates of H and D atom onto the surface, respectively. The RRKM is applied in a microcanonical manner consistent with that laid out in Reference [182] for surface collisions. We assume that the rate-determining step in the chemisorption is a unimolecular process where the reactant is an activated pre-cursor complex (PC) [183] that corresponds to the H/D atom in a physisorbed complex with the surface at a distance 4Å [150]. Activated PC is also assumed to have sufficiently rapid intramolecular energy distribution such that the statistical mechanical treatment in mTST is applicable. The RRKM microcanonical rate is 𝑘 (𝐸) = 𝑁 ‡ (𝐸 − 𝐸 ‡ ) . ℎ 𝜌PC (𝐸) (5.6) 83 Here 𝑁 ‡ is the transition state sum of states, 𝐸 ‡ is the energy difference between the PC and the transition state geometry (including zero-point energy contributions), and 𝜌 is the PC density of states. Note that the PC consists of a structureless H/D atom and the C atoms that vibrate at the graphene phonon frequencies with no chemical interactions to the H/D atom. As a result, the denominator in the Eq. 5.6 is solely determined by the graphene phonon frequencies and the isotopic difference of the microcanonical rate then becomes † † 𝑘 D (𝐸) 𝑁D (𝐸 − 𝐸 D ) = . 𝑘 H (𝐸) 𝑁 † (𝐸 − 𝐸 † ) H (5.7) H Fig. 5.5 presents the microcanonical rate ratio of the H/D chemisorption, obtained by calculations with Eq. 5.7 on the GFN-xTB potential energy surface. The sum of the states in Eq. 5.7 is calculated with the Beyer-Swinehart algorithm [184] at the level of separable harmonic vibrations. It is predicted that a large isotope effect more than 250% favors D sticking is observed, in contrary to the experimental observation of almost no isotope effect. The inconsistency between the mTST result and the experiment measurement emphasizes the highly non-equilibrium nature of the chemisorption process in the scattering experiment and the ultrafast reaction timescale that prohibits a sufficient energy redistribution of the reactive system, Figure 5.5: Isotope differences in the H/D-on-graphene scattering events predicted with theories at different levels. The ratio of microcanonical transition state theory rates (mTST, Eq. 5.7, with and without the inclusion of the zero-point energy (ZPE) in the transition state complex) and Eckart barrier transmission probabilities are presented to examine the isotope effect. The filled grey area indicates no isotope effect. 84 Figure 5.6: Comparing theory with experiment—Scattering distribution of H/D collision with graphene. (a)-(c) H scattering experiment; (d)-(f) H scattering theory; (g)-(i) D scattering experiment; (j)-(l) D scattering theory. In all figures, the scattering energy, 𝐸 𝑠 , is shown along the radial coordinate as a fraction of 𝐸𝑖 and the scattering angle, 𝜗𝑠 , is shown on the polar coordinate. All theoretical scattering distributions were obtained with RPMD trajectories under the influence of GFN-xTB potential energies. Each heat map is multiplied by the indicated number shown in red. We encourage the reader to compare this figure to Fig. S10 of [106]. Such a comparison shows that the new GFN-xTB PES provides an equally accurate representation of the forces during scattering. 𝐸𝑖 = 0.99 eV for H and 𝐸𝑖 = 0.94 eV for D scattering experiments. 85 which leads to this particular breakdown of the statistical rate theory. Furthermore, the approximation of those non-phonon active modes as harmonic vibrations might also degrade the accuracy of mTST result. To this end, those inadequacies in the mTST modelling urges the need for an appropriate inclusion of the dynamical effects and vibrational anharmonicity in the simulation of H/D-on-graphene scatterings. As a simplified alternative, we also estimate the isotopic difference between H and D atom transmission through the adsorption barrier with the one-dimensional Eckart barrier model [185, 186]. The result is also presented in Fig. 5.5. Significant difference is shown at lower incidence energies between the transmission probabilities that the scattering system tunnels through the adsorption barrier; while at the energies that higher than the barrier, Eckart barrier model predicts exactly the same classical transmission for H/D atoms and no isotope difference is seen. This reduced-dimensional transmission model is extremely intuitive to describe the scattering process, however is also a grossly oversimplified one that lacks the capability to describe the collective graphene motions which are orthogonal to the reaction coordinate. The latter has proven to be essential for quantitatively understanding the reactivity in the H/D scattering on the graphene surface [106]. C. Trajectories and comparison to Experiment Trajectory calculations (cMD or RPMD) were initialized following the procedure in Reference [106] For a given incidence energy 𝐸𝑖 and angle 𝜗𝑖 at which comparison with experiment is to be made, the following steps are performed. (1) The initial configurations of the graphene flake are sampled from a trajectory propagated using either with cMD or RPMD and thermalized at the experiment’s temperature of 300 K. More specifically, the trajectory is preequilibrated for 40 ps and then propagated for another 30 ps, during which initialization geometries are sampled. (2) The collision point of H/D atom on the surface (𝑥 𝑜 , 𝑦 𝑜 , 𝑧 𝑜 = 0 Å) is randomly selected from within the center unit cell of the graphene flake. (3) The azimuthal scattering angle 𝜑𝑖 is drawn from a uniform distribution between [0, 2𝜋). 86 Figure 5.7: Experimentally derived total survival fluxes. (4) The initial height of the projectile atom is set at 4 Å from the surface. The (centroid) position to initialize the projectile atom (𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 ) is then calculated using (𝑥 𝑜 , 𝑦 𝑜 , 𝑧 𝑜 ), 𝜗𝑖 and 𝜑𝑖 , i.e. 𝑥𝑖 = 𝑥 𝑜 + (𝑧𝑖 − 𝑧 𝑜 ) tan 𝜗𝑖 cos 𝜑𝑖 , 𝑦𝑖 = 𝑦 𝑜 + (𝑧𝑖 − 𝑧 𝑜 ) tan 𝜗𝑖 sin 𝜑𝑖 , 𝑧𝑖 = 4 Å. (5.8) (5) The magnitude of the H/D atom’s velocity vector is calculated from 𝐸𝑖 , and its direction is found geometrically by connecting a line between the initial position of the H atom and the collision point. (6) The relative positions and momenta of the internal ring-polymer beads with respect to those of the centroid are thermally sampled at 300 K. Fig. 5.6 presents the scattering distribution at 𝜗𝑖 between 60◦ and 30◦ from experimental measurement and RPMD simulations. They quantitatively agree with each other. D. Experimentally derived survival and sticking probabilities Experimental data like that of Fig. 2 and Fig. 5.6 were used to derive the H and D survival probabilities as function of normal incidence translational energy. The experiments are only sensitive to scattered atoms within 2.8◦ of the plane defined by the incidence atom beam and the surface normal. We integrate over the out of plane scattered flux assuming cylindrical symmetry of the angular distribution with 87 Figure 5.8: Phonon density of states spectrum calculated with GFN-xTB and EMFT-REBO potential at 300 K, respectively. The phonon power spectrum is calculated as the Fourier transform of the velocity autocorrelation function of the carbon atoms, which is obtained from classical molecular dynamics calculations. Reference spectrum are took from Fig. 5 of Reference [176]. respect to a line coincident with the specularly scattered atoms. This assumption was validated in Reference [161]. Fig. 5.7 shows the Survival fluxes at incidence angles from 60◦ to 30◦ —it increases with 𝜗𝑖 , reaching a plateau below 𝐸 n = 0.35 eV, with 𝐸 𝑛 = 𝐸𝑖 cos2 𝜗𝑖 the normal incidence energy. We set the absolute scale of the survival probability to 1 at this plateau. From this, we easily obtain the sticking probabilities shown in Fig. 2. E. Comparing GFN-xTB with EMFT-REBO potential In this work, we have developed a new PES for the H/graphene system using the GFN-xTB method. Using it, we are able to accurately reproduce experimentally measured H and D atom scattering distributions using cMD and RPMD. Previously, using a PES [106] based on reactive empirical bond order (EMFT-REBO) formalism [187] with training data from quantum embedding calculations [188–190]. Fig. 5.8 shows the graphene phonon density of states spectra (PDOS) generated from xTB and REBO potential, respectively, compared with the published periodic DFT calculations at the GGA level [176] (black solid line). The phonon power spectrum is calculated as the Fourier transform of the velocity autocorrelation function of the carbon atoms, which is obtained from classical molecular dynamics calculations. The PDOS spectrum obtained from the GFN-xTB PES is shown in red, while that obtained with EMFT-REBO is shown in green. The three theories model a free 88 Figure 5.9: Sticking probability of H/D scattering on graphene as a function of normal incidence energy. The calculation is performed using the EMFT-REBO PES reported previously [106] . standing graphene sheet. It is immediately clear that GFN-xTB describes the high frequency phonons of graphene much better than does EMFT-REBO. Furthermore, the phonon band gap between 1000 and 1300 cm−1 is not captured by EMFT-REBO. The high frequency phonons so poorly described by EMFT-REBO are attributed to in-plane optical phonons, precisely those so critical to the H sticking process [106]. These deficiencies in the EMFT-REBO PES leads to an absence of significant nuclear quantum effects. Fig. 5.9 shows EMFT-REBO based sticking probabilities for H and D using both cMD and RPMD trajectories. No difference between classical and quantum dynamics is seen. Contrast this with Fig. 3, which show analogous results using GFN-xTB. F. Interaction time of H/D with graphene See Fig. 5.10. G. Robustness with respect to the RPMD simulation protocol In Fig. 5.11, classical molecular dynamics (cMD, in black) and RPMD (in red) are reproduced from Fig. 2, with the label 𝑇int,𝐻 = 300 K indicating the non-equilibrium RPMD method reported in the main text. Also shown are RPMD simulations with the H/D internal ring-polymer temperature is set to 𝑇int,𝐻 = E𝑛 /𝑘 𝐵 = E𝑖 cos2 𝜗𝑖 /𝑘 𝐵 , following the ‘free-particle direct shooting’ protocol [169]. All other simulation 89 Figure 5.10: Interaction time of the scattered H/D atom with the graphene surface from simulations. The symbols present the mean value of the interaction time from trajectory statistics and the shadow area in black and red, respectively, present the distribution width. The interacting period is defined that the distance of H/D projectile and the nearest neighbouring C atom is shorter than 2.5 Å. Figure 5.11: H atom sticking probability as a function of normal incidence energy—testing sensitivity with respect to the initialization of the H/D internal ring-polymer temperature. details in these two sets of RPMD results are kept to be the same. The fact that both RPMD results are in complete agreement indicates that the results reported here are insensitive to the protocol used for initializing the internal temperature of the H/D ring-polymer in the non-equilibrium simulations. This is also consistent with the conclusion from the main text that the NQE associated with the H/D atom is negligible for the sticking dynamics studied here (Fig. 5.3). 90 BIBLIOGRAPHY 1. London, F. On Condensed Helium at Absolute Zero. Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 153, 576–583. doi:10.1098/rspa.1936.0024 (1936). 2. Hu, S. et al. Proton Transport through One-Atom-Thick Crystals. Nature 516, 227–230. doi:10.1038/nature14015 (2014). 3. Kretchmer, J. S., Boekelheide, N., Warren, J. J., Winkler, J. R., Gray, H. B. & Miller III, T. F. Fluctuating Hydrogen-Bond Networks Govern Anomalous Electron Transfer Kinetics in a Blue Copper Protein. Proceedings of the National Academy of Sciences of the United States of America 115, 6129– 6134. doi:10.1073/pnas.1805719115 (2018). 4. Marais, A. et al. The Future of Quantum Biology. Journal of the Royal Society Interface 15, 20180640. doi:10.1098/rsif.2018.0640 (2018). 5. Dyson, G. Some Aspects of the Vibration Theory of Odor. Perfumery and Essential Oil Record 19 (1928). 6. Xu, X., Zheng, J., Yang, K. R. & Truhlar, D. G. Photodissociation Dynamics of Phenol: Multistate Trajectory Simulations Including Tunneling. Journal of the American Chemical Society 136, 16378–16386. doi:10 . 1021 / ja509016a (2014). 7. Levine, B. G. & Martinez, T. J. Isomerization through Conical Intersections. Annual Review of Physical Chemistry 58, 613–634. doi:10.1146/annurev. physchem.57.032905.104612 (2007). 8. Shih, C. et al. Tryptophan-Accelerated Electron Flow through Proteins. Science 320, 1760–1762. doi:10.1126/science.1158241 (2008). 9. Fujita, K., Nakamura, N., Ohno, H., Leigh, B., Niki, K., Gray, H. & Richards, J. Mimicking Protein-Protein Electron Transfer: Voltammetry of Pseudomonas Aeruginosa Azurin and the Thermus Thermophilus Cu(A) Domain at 𝜔Derivatized Self-Assembled-Monolayer Gold Electrodes. Journal of the American Chemical Society 126, 13954–13961. doi:10.1021/ja047875o (2004). 10. Nahler, N. H., White, J. D., LaRue, J., Auerbach, D. J. & Wodtke, A. M. Inverse Velocity Dependence of Vibrationally Promoted Electron Emission from a Metal Surface. Science 321, 1191–1194. doi:10.1126/science. 1160040 (2008). 11. Yao, Y. & Giapis, K. P. Dynamic Molecular Oxygen Production in Cometary Comae. Nature Communications 8. doi:10.1038/ncomms15298 (2017). 12. Ehrenfest, P. Bemerkung über die angenäherte Gültigkeit der klassischen Mechanik innerhalb der Quantenmechanik. Zeitschrift für Physik A Hadrons and Nuclei 45, 455–457. doi:10.1007/BF01329203 (1927). 91 13. Tully, J. Molecular Dynamics with Electronic Transitions. Journal of Chemical Physics 93, 1061–1071. doi:10.1063/1.459170 (1990). 14. Xie, C., Ma, J., Zhu, X., Yarkony, D. R., Xie, D. & Guo, H. Nonadiabatic Tunneling in Photodissociation of Phenol. Journal of the American Chemical Society 138, 7828–7831. doi:10.1021/jacs.6b03288 (2016). 15. Hammes-Schiffer, S. & Stuchebrukhov, A. A. Theory of Coupled Electron and Proton Transfer Reactions. Chemical Reviews 110, 6939–6960. doi:10. 1021/cr1001436 (2010). 16. Huynh, M. H. V. & Meyer, T. J. Proton-Coupled Electron Transfer. Chemical Reviews 107, 5004–5064. doi:10.1021/cr0500030 (2007). 17. Migliore, A., Polizzi, N. F., Therien, M. J. & Beratan, D. N. Biochemistry and Theory of Proton-Coupled Electron Transfer. Chemical Reviews 114, 3381–3465. doi:10.1021/cr4006654 (2014). 18. Feynman, R. P. & Hibbs, A. R. Quantum Mechanics and Path Integrals (McGraw-Hill, 1965). 19. Chandler, D. & Wolynes, P. Exploiting the Isomorphism between QuantumTheory and Classical Statistical-Mechanics of Polyatomic Fluids. Journal of Chemical Physics 74, 4078–4095. doi:10.1063/1.441588 (1981). 20. Craig, I. & Manolopoulos, D. Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 121, 3368–3373. doi:10.1063/1.1777575 (2004). 21. Habershon, S., Manolopoulos, D. E., Markland, T. E. & Miller III, T. F. RingPolymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space. Annual Review of Physical Chemistry 64, 387–413. doi:10 . 1146 / annurev - physchem - 040412 110122 (2013). 22. Cao, J. & Voth, G. The Formulation of Quantum-Statistical Mechanics Based on the Feynman Path Centroid Density. II. Dynamical Properties. Journal of Chemical Physics 100, 5106–5117. doi:10.1063/1.467176 (1994). 23. Voth, G. Path-Integral Centroid Methods in Quantum Statistical Mechanics and Dynamics. Advances in Chemical Physics, Vol XCIII 93, 135–218. doi:10.1002/9780470141526.ch4 (1996). 24. Jang, S. & Voth, G. A Derivation of Centroid Molecular Dynamics and Other Approximate Time Evolution Methods for Path Integral Centroid Variables. Journal of Chemical Physics 111, 2371–2384. doi:10 . 1063 / 1 . 479515 (1999). 92 25. Richardson, J. O. & Althorpe, S. C. Ring-Polymer Molecular Dynamics RateTheory in the Deep-Tunneling Regime: Connection with Semiclassical Instanton Theory. Journal of Chemical Physics 131. doi:10.1063/1.3267318 (2009). 26. Rossi, M., Ceriotti, M. & Manolopoulos, D. E. How to Remove the Spurious Resonances from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 140. doi:10.1063/1.4883861 (2014). 27. Ceriotti, M., More, J. & Manolopoulos, D. E. I-PI: A Python Interface for Ab Initio Path Integral Molecular Dynamics Simulations. Computer Physics Communications 185, 1019–1026. doi:10 . 1016 / j . cpc . 2013 . 10 . 027 (2014). 28. Markland, T. E. & Manolopoulos, D. E. An Efficient Ring Polymer Contraction Scheme for Imaginary Time Path Integral Simulations. Journal of Chemical Physics 129. doi:10.1063/1.2953308 (2008). 29. Marsalek, O. & Markland, T. E. Ab Initio Molecular Dynamics with Nuclear Quantum Effects at Classical Cost: Ring Polymer Contraction for Density Functional Theory. Journal of Chemical Physics 144. doi:10 . 1063 / 1 . 4941093 (2016). 30. Miller III, T. F. & Manolopoulos, D. Quantum Diffusion in Liquid ParaHydrogen from Ring-Polymer Molecular Dynamics. Journal of Chemical Physics 122. doi:10.1063/1.1893956 (2005). 31. Miller III, T. F. & Manolopoulos, D. Quantum Diffusion in Liquid Water from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 123. doi:10.1063/1.2074967 (2005). 32. Habershon, S., Fanourgakis, G. S. & Manolopoulos, D. E. Comparison of Path Integral Molecular Dynamics Methods for the Infrared Absorption Spectrum of Liquid Water. Journal of Chemical Physics 129. doi:10.1063/1.2968555 (2008). 33. Boekelheide, N., Salomon-Ferrer, R. & Miller III, T. F. Dynamics and Dissipation in Enzyme Catalysis. Proceedings of the National Academy of Sciences of the United States of America 108, 16159–16163. doi:10 . 1073 / pnas.1106397108 (2011). 34. Cao, J., Minichino, C. & Voth, G. The Computation of Electron-Transfer Rates: the Nonadiabatic Instanton Solution. Journal of Chemical Physics 103, 1391–1399. doi:10.1063/1.469762 (1995). 35. Schwieters, C. & Voth, G. The Semiclassical Calculation of Nonadiabatic Tunneling Rates. Journal of Chemical Physics 108, 1055–1062. doi:10 . 1063/1.475467 (1998). 93 36. Jang, S. & Cao, J. Nonadiabatic Instanton Calculation of Multistate Electron Transfer Reaction Rate: Interference Effects in Three and Four States Systems. Journal of Chemical Physics 114, 9959–9968. doi:10.1063/1.1371262 (2001). 37. Schwieters, C. & Voth, G. Extension of Path Integral Quantum Transition State Theory to the Case of Nonadiabatic Activated Dynamics. Journal of Chemical Physics 111, 2869–2877. doi:10.1063/1.479569 (1999). 38. Liao, J. & Voth, G. A Centroid Molecular Dynamics Approach for Nonadiabatic Dynamical Processes in Condensed Phases: The Spin-Boson Case. Journal of Physical Chemistry B 106, 8449–8455. doi:10.1021/jp020978d (2002). 39. Ananth, N. Mapping Variable Ring Polymer Molecular Dynamics: A PathIntegral Based Method for Nonadiabatic Processes. Journal of Chemical Physics 139. doi:10.1063/1.4821590 (2013). 40. Richardson, J. O. & Thoss, M. Communication: Nonadiabatic Ring-Polymer Molecular Dynamics. Journal of Chemical Physics 139. doi:10.1063/1. 4816124 (2013). 41. Shushkov, P., Li, R. & Tully, J. C. Ring Polymer Molecular Dynamics with Surface Hopping. Journal of Chemical Physics 137. doi:10 . 1063 / 1.4766449 (2012). 42. Menzeleev, A. R., Bell, F. & Miller III, T. F. Kinetically Constrained RingPolymer Molecular Dynamics for Non-Adiabatic Chemical Reactions. Journal of Chemical Physics 140. doi:10.1063/1.4863919 (2014). 43. Kretchmer, J. S. & Miller III, T. F. Kinetically-Constrained Ring-Polymer Molecular Dynamics for Non-Adiabatic Chemistries Involving Solvent and Donor-Acceptor Dynamical Effects. Faraday Discussions 195, 191–214. doi:10.1039/c6fd00143b (2016). 44. Kretchmer, J. S. & Miller III, T. F. Tipping the Balance between Concerted versus Sequential Proton-Coupled Electron Transfer. Inorganic Chemistry 55, 1022–1031. doi:10.1021/acs.inorgchem.5b01821 (2016). 45. Kretchmer, J. S. & Miller III, T. F. Direct Simulation of Proton-Coupled Electron Transfer across Multiple Regimes. Journal of Chemical Physics 138. doi:10.1063/1.4797462 (2013). 46. Shakib, F. A. & Huo, P. Ring Polymer Surface Hopping: Incorporating Nuclear Quantum Effects into Nonadiabatic Molecular Dynamics Simulations. Journal of Physical Chemistry Letters 8, 3073–3080. doi:10.1021/acs. jpclett.7b01343 (2017). 47. Ananth, N. & Miller III, T. F. Exact Quantum Statistics for Electronically Nonadiabatic Systems Using Continuous Path Variables. Journal of Chemical Physics 133. doi:10.1063/1.3511700 (2010). 94 48. Ananth, N., Venkataraman, C. & Miller, W. H. Semiclassical Description of Electronically Nonadiabatic Dynamics via the Initial Value Representation. Journal of Chemical Physics 127. doi:10.1063/1.2759932 (2007). 49. Sun, X., Wang, H. & Miller, W. Semiclassical Theory of Electronically Nonadiabatic Dynamics: Results of a Linearized Approximation to the Initial Value Representation. Journal of Chemical Physics 109, 7064–7074. doi:10. 1063/1.477389 (1998). 50. Sun, X. & Miller, W. Semiclassical Initial Value Representation for Electronically Nonadiabatic Molecular Dynamics. Journal of Chemical Physics 106, 6346–6353. doi:10.1063/1.473624 (1997). 51. Huo, P. & Coker, D. F. Iterative Linearized Density Matrix Propagation for Modeling Coherent Excitation Energy Transfer in Photosynthetic Light Harvesting. Journal of Chemical Physics 133. doi:10.1063/1.3498901 (2010). 52. Huo, P., Miller III, T. F. & Coker, D. F. Communication: Predictive Partial Linearized Path Integral Simulation of Condensed Phase Electron Transfer Dynamics. Journal of Chemical Physics 139. doi:10 . 1063 / 1 . 4826163 (2013). 53. Kapral, R. & Ciccotti, G. Mixed Quantum-Classical Dynamics. Journal of Chemical Physics 110, 8919–8929. doi:10.1063/1.478811 (1999). 54. Nielsen, S., Kapral, R. & Ciccotti, G. Mixed Quantum-Classical Surface Hopping Dynamics. Journal of Chemical Physics 112, 6543–6553. doi:10. 1063/1.481225 (2000). 55. Parrinello, M. & Rahman, A. Study of an F-Center in Molten KCL. Journal of Chemical Physics 80, 860–867. doi:10.1063/1.446740 (1984). 56. Schmidt, J. R. & Tully, J. C. Path-Integral Simulations beyond the Adiabatic Approximation. Journal of Chemical Physics 127. doi:10.1063/1.2757170 (2007). 57. Duke, J. R. & Ananth, N. Mean Field Ring Polymer Molecular Dynamics for Electronically Nonadiabatic Reaction Rates. Faraday Discussions 195, 253–268. doi:10.1039/c6fd00123h (2016). 58. Hele, T. J. H. An Electronically Non-Adiabatic Generalization of Ring Polymer Molecular Dynamics Master’s Thesis, Exeter College, Oxford University. 2011. 59. Schmidt, J. R., Parandekar, P. V. & Tully, J. C. Mixed Quantum-Classical Equilibrium: Surface Hopping. Journal of Chemical Physics 129. doi:10. 1063/1.2955564 (2008). 60. Chandler, D. Introduction to Modern Statistical Mechanics (Oxford University Press, 1987). 95 61. Nielsen, S., Kapral, R. & Ciccotti, G. Statistical Mechanics of QuantumClassical Systems. Journal of Chemical Physics 115, 5805–5815. doi:10. 1063/1.1400129 (2001). 62. Wigner, E. On the Quantum Correction for Thermodynamic Equilibrium. Physical Review 40, 0749–0759. doi:10.1103/PhysRev.40.749 (1932). 63. Hillery, M., O’Connell, R. F., Scully, M. O. & Wigner, E. P. Distribution Functions in Physics: Fundamentals (Springer, 1997). 64. Mac Kernan, D., Ciccotti, G. & Kapral, R. Surface-Hopping Dynamics of a Spin-Boson System. Journal of Chemical Physics 116, 2346–2353. doi:10. 1063/1.1433502 (2002). 65. Press, W. H. Numerical Recipes 3rd Edition: the Art of Scientific Computing (Cambridge University Press, 2007). 66. Marston, C. & Balintkurti, G. The Fourier Grid Hamiltonian Method for Bound-State Eigenvalues and Eigenfunctions. Journal of Chemical Physics 91, 3571–3576. doi:10.1063/1.456888 (1989). 67. Colbert, D. & Miller, W. A Novel Discrete Variable Representation for Quantum-Mechanical Reactive Scattering via the S-Matrix Kohn Method. Journal of Chemical Physics 96, 1982–1991. doi:10 . 1063 / 1 . 462100 (1992). 68. Jain, A. & Subotnik, J. E. Surface Hopping, Transition State Theory, and Decoherence. II. Thermal Rate Constants and Detailed Balance. Journal of Chemical Physics 143. doi:10.1063/1.4930549 (2015). 69. Hammes-Schiffer, S. & Tully, J. Proton-Transfer in Solution: MolecularDynamics with Quantum Transitions. Journal of Chemical Physics 101, 4657–4667. doi:10.1063/1.467455 (1994). 70. Subotnik, J. E. & Shenvi, N. A New Approach to Decoherence and Momentum Rescaling in the Surface Hopping Algorithm. Journal of Chemical Physics 134. doi:10.1063/1.3506779 (2011). 71. Prezhdo, O. & Rossky, P. Relationship between Quantum Decoherence Times and Solvation Dynamics in Condensed Phase Chemical Systems. Physical Review Letters 81, 5294–5297. doi:10 . 1103 / PhysRevLett . 81 . 5294 (1998). 72. Frenkel, D. & Smit, B. Understanding Molecular Simulation: From Algorithms to Applications (Academic press, 2001). 73. Feit, M., Fleck, J. & Steiger, A. Solution of the Schrodinger Equation by A Spectral Method. Journal of Computational Physics 47, 412–433. doi:10. 1016/0021-9991(82)90091-2 (1982). 96 74. Jang, S., Schwieters, C. & Voth, G. A Modification of Path Integral Quantum Transition State Theory for Asymmetric and Metastable Potentials. Journal of Physical Chemistry A 103, 9527–9538. doi:10.1021/jp992190+ (1999). 75. Manby, F. R., Stella, M., Goodpaster, J. D. & Miller III, T. F. A Simple, Exact Density-Functional-Theory Embedding Scheme. Journal of Chemical Theory and Computation 8, 2564–2568. doi:10.1021/ct300544e (2012). 76. Kubo, R. Generalized Cumulant Expansion Method. Journal of the Physical Society of Japan 17, 1100–&. doi:10.1143/JPSJ.17.1100 (1962). 77. Bell, M. T. Unpublished contribution to an RPMD group meeting in Oxford. 2005. 78. Craig, I. & Manolopoulos, D. Chemical Reaction Rates from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 122. doi:10.1063/1. 1850093 (2005). 79. Craig, I. & Manolopoulos, D. A Refined Ring Polymer Molecular Dynamics Theory of Chemical Reaction Rates. Journal of Chemical Physics 123. doi:10.1063/1.1954769 (2005). 80. Suleimanov, Y. V., Javier Aoiz, F. & Guo, H. Chemical Reaction Rate Coefficients from Ring Polymer Molecular Dynamics: Theory and Practical Applications. Journal of Physical Chemistry A 120, 8488–8502. doi:10 . 1021/acs.jpca.6b07140 (2016). 81. Collepardo-Guevara, R., Suleimanov, Y. V. & Manolopoulos, D. E. Bimolecular Reaction Rates from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 130. doi:10.1063/1.3127145 (2009). 82. Suleimanov, Y. V., Collepardo-Guevara, R. & Manolopoulos, D. E. Bimolecular Reaction Rates from Ring Polymer Molecular Dynamics: Application to H + CH4 → H2 + CH3 . Journal of Chemical Physics 134. doi:10.1063/1. 3533275 (2011). 83. Perez de Tudela, R., Aoiz, F. J., Suleimanov, Y. V. & Manolopoulos, D. E. Chemical Reaction Rates from Ring Polymer Molecular Dynamics: Zero Point Energy Conservation in Mu + H2 → MuH + H. Journal of Physical Chemistry Letters 3, 493–497. doi:10.1021/jz201702q (2012). 84. Suleimanov, Y. V., Perez de Tudela, R., Jambrina, P. G., Castillo, J. F., SaezRabanos, V., Manolopoulos, D. E. & Javier Aoiz, F. A Ring Polymer Molecular Dynamics Study of the Isotopologues of the H + H2 Reaction. Physical Chemistry Chemical Physics 15, 3655–3665. doi:10 . 1039 / c2cp44364c (2013). 85. Li, Y., Suleimanov, Y. V., Yang, M., Green, W. H. & Guo, H. Ring Polymer Molecular Dynamics Calculations of Thermal Rate Constants for the O(3 P) + CH4 → OH + CH3 Reaction: Contributions of Quantum Effects. Journal of Physical Chemistry Letters 4, 48–52. doi:10.1021/jz3019513 (2013). 97 86. Li, Y., Suleimanov, Y. V., Li, J., Green, W. H. & Guo, H. Rate Coefficients and Kinetic Isotope Effects of the X + CH4 → CH3 + HX (X = H, D, Mu) Reactions from Ring Polymer Molecular Dynamics. Journal of Chemical Physics 138. doi:10.1063/1.4793394 (2013). 87. Suleimanov, Y. V. Surface Diffusion of Hydrogen on Ni(100) from Ring Polymer Molecular Dynamics. Journal of Physical Chemistry C 116, 11141– 11153. doi:10.1021/jp302453z (2012). 88. Castillo, J. F. & Suleimanov, Y. V. A Ring Polymer Molecular Dynamics Study of the OH + H2 (D2 ) Reaction. Physical Chemistry Chemical Physics 19, 29170–29176. doi:10.1039/c7cp05266a (2017). 89. Miller III, T. F. Isomorphic Classical Molecular Dynamics Model for an Excess Electron in a Supercritical Fluid. Journal of Chemical Physics 129. doi:10.1063/1.3013357 (2008). 90. Menzeleev, A. R. & Miller III, T. F. Ring Polymer Molecular Dynamics beyond the Linear Response Regime: Excess Electron Injection and Trapping in Liquids. Journal of Chemical Physics 132. doi:10 . 1063 / 1 . 3292576 (2010). 91. Pierre, S., Duke, J. R., Hele, T. J. H. & Ananth, N. A Mapping Variable Ring Polymer Molecular Dynamics Study of Condensed Phase Proton-Coupled Electron Transfer. Journal of Chemical Physics 147. doi:10 . 1063 / 1 . 4986517 (2017). 92. Chowdhury, S. N. & Huo, P. Coherent State Mapping Ring Polymer Molecular Dynamics for Non-Adiabatic Quantum Propagations. Journal of Chemical Physics 147. doi:10.1063/1.4995616 (2017). 93. Tao, X., Shushkov, P. & Miller III, T. F. Path-Integral Isomorphic Hamiltonian for Including Nuclear Quantum Effects in Non-Adiabatic Dynamics. Journal of Chemical Physics 148, 102327. doi:10.1063/1.5005544 (2018). 94. Miller, W. H. & Cotton, S. J. Classical Molecular Dynamics Simulation of Electronically Non-Adiabatic Processes. Faraday Discussions 195, 9–30. doi:10.1039/c6fd00181e (2016). 95. Cotton, S. J. & Miller, W. H. Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations. Journal of Physical Chemistry A 117, 7190–7194. doi:10.1021/jp401078u (2013). 96. Wang, L., Sifain, A. E. & Prezhdo, O. V. Fewest Switches Surface Hopping in Liouville Space. Journal of Physical Chemistry Letters 6, 3827–3833. doi:10.1021/acs.jpclett.5b01502 (2015). 97. Hammes-Schiffer, S. & Tully, J. Nonadiabatic Transition-State Theory and Multiple Potential-Energy Surface Molecular-Dynamics of Infrequent Events. Journal of Chemical Physics 103, 8528–8537. doi:10 . 1063 / 1 . 470162 (1995). 98 98. Voth, G., Chandler, D. & Miller, W. Rigorous Formulation of Quantum Transition-State Theory and Its Dynamical Corrections. Journal of Chemical Physics 91, 7749–7760. doi:10.1063/1.457242 (1989). 99. Hele, T. J. H. & Althorpe, S. C. Derivation of a True (𝑡 → 0+ ) Quantum Transition-State Theory. I. Uniqueness and Equivalence to Ring-Polymer Molecular Dynamics Transition-State-Theory. Journal of Chemical Physics 138. doi:10.1063/1.4792697 (2013). 100. Neuhauser, D., Baer, M., Judson, R. & Kouri, D. The Application of TimeDependent Wavepacket Methods to Reactive Scattering. Computer Physics Communications 63, 460–481. doi:10 . 1016 / 0010 - 4655(91 ) 90270 - U (1991). 101. Bennett, C. H. Algorithms for Chemical Computations (American Chemical Society, 1977). 102. Chandler, D. Statistical-Mechanics of Isomerization Dynamics in Liquids and Transition-State Approximation. Journal of Chemical Physics 68, 2959– 2970. doi:10.1063/1.436049 (1978). 103. Li, X., Tully, J., Schlegel, H. & Frisch, M. Ab Initio Ehrenfest Dynamics. Journal of Chemical Physics 123. doi:10.1063/1.2008258 (2005). 104. Miller, W., Schwartz, S. & Tromp, J. Quantum-Mechanical Rate Constants for Bimolecular Reactions. Journal of Chemical Physics 79, 4889–4898. doi:10.1063/1.445581 (1983). 105. Parandekar, P. & Tully, J. Mixed Quantum-Classical Equilibrium. Journal of Chemical Physics 122. doi:10.1063/1.1856460 (2005). 106. Jiang, H. et al. Imaging Covalent Bond Formation by H Atom Scattering from Graphene. Science 364, 379–382. doi:10.1126/science.aaw6378 (2019). 107. Liu, Q., Zhang, L., Li, Y. & Jiang, B. Ring Polymer Molecular Dynamics in Gas-Surface Reactions: Inclusion of Quantum Effects Made Simple. Journal of Physical Chemistry Letters 10, 7475–7481. doi:10.1021/acs.jpclett. 9b02570 (2019). 108. Kaur, R. & Welsch, R. Probing Photodissociation Dynamics Using Ring Polymer Molecular Dynamics. Journal of Chemical Physics 150, 114105. doi:10.1063/1.5086218 (2019). 109. Yu, Q. & Bowman, J. M. Classical, Thermostated Ring Polymer, and Quantum VSCF/VCI Calculations of IR Spectra of H7 O+3 and H9 O+4 (Eigen) and Comparison with Experiment. Journal of Physical Chemistry A 123, 1399– 1409. doi:10.1021/acs.jpca.8b11603 (2019). 99 110. Novikov, I. S., Suleimanov, Y. V. & Shapeev, A. V. Automated Calculation of Thermal Rate Coefficients Using Ring Polymer Molecular Dynamics and Machine-Learning Interatomic Potentials with Active Learning. Physical Chemistry Chemical Physics 20, 29503–29512. doi:10.1039/c8cp06037a (2018). 111. Wilkins, D. M., Manolopoulos, D. E., Pipolo, S., Laage, D. & Hynes, J. T. Nuclear Quantum Effects in Water Reorientation and Hydrogen-Bond Dynamics. Journal of Physical Chemistry Letters 8, 2602–2607. doi:10.1021/ acs.jpclett.7b00979 (2017). 112. Kreis, K., Kremer, K., Potestio, R. & Tuckerman, M. E. From Classical to Quantum and Back: Hamiltonian Adaptive Resolution Path Integral, Ring Polymer, and Centroid Molecular Dynamics. Journal of Chemical Physics 147. doi:10.1063/1.5000701 (2017). 113. Rossi, M., Ceriotti, M. & Manolopoulos, D. E. Nuclear Quantum Effects in H+ and OH− Diffusion along Confined Water Wires. Journal of Physical Chemistry Letters 7, 3001–3007. doi:10 . 1021 / acs . jpclett . 6b01093 (2016). 114. Kenion, R. L. & Ananth, N. Direct Simulation of Electron Transfer in the Cobalt Hexammine(II/III) Self-Exchange Reaction. Physical Chemistry Chemical Physics 18, 26117–26124. doi:10.1039/c6cp04882j (2016). 115. Ceriotti, M., Cuny, J., Parrinello, M. & Manolopoulos, D. E. Nuclear Quantum Effects and Hydrogen Bond Fluctuations in Water. Proceedings of the National Academy of Sciences of the United States of America 110, 15591– 15596. doi:10.1073/pnas.1308560110 (2013). 116. Menzeleev, A. R., Ananth, N. & Miller III, T. F. Direct Simulation of Electron Transfer Using Ring Polymer Molecular Dynamics: Comparison with Semiclassical Instanton Theory and Exact Quantum Methods. Journal of Chemical Physics 135. doi:10.1063/1.3624766 (2011). 117. Chowdhury, S. N. & Huo, P. State Dependent Ring Polymer Molecular Dynamics for Investigating Excited Nonadiabatic Dynamics. Journal of Chemical Physics 150. doi:10.1063/1.5096276 (2019). 118. Welsch, R., Song, K., Shi, Q., Althorpe, S. C. & Miller III, T. F. NonEquilibrium Dynamics from RPMD and CMD. Journal of Chemical Physics 145. doi:10.1063/1.4967958 (2016). 119. Duke, J. R. & Ananth, N. Simulating Excited State Dynamics in Systems with Multiple Avoided Crossings Using Mapping Variable Ring Polymer Molecular Dynamics. Journal of Physical Chemistry Letters 6, 4219–4223. doi:10.1021/acs.jpclett.5b01957 (2015). 120. Miller, W. Path Integral-Representation of Reaction-Rate Constant in QuantumMechanical Transition-State Theory. Journal of Chemical Physics 63, 1166– 1172. doi:10.1063/1.431444 (1975). 100 121. Slater, N. Theory of Unimolecular Reactions (Cornell University Press, 1959). 122. Davies, J., Green, N. & Pilling, M. The Testing of Models for Unimolecular Decomposition via Inverse Laplace Transformation of Experimental Recombination Rate Data. Chemical Physics Letters 126, 373–379. doi:10.1016/ S0009-2614(86)80101-4 (1986). 123. Firaha, D. S., Doentgen, M., Berkels, B. & Leonhard, K. Pressure-Dependent Rate Constant Predictions Utilizing the Inverse Laplace Transform: A Victim of Deficient Input Data. Acs Omega 3, 8212–8219. doi:10.1021/acsomega. 8b00311 (2018). 124. Bryan, R. Maximum-Entropy Analysis of Oversampled Data Problems. European Biophysics Journal 18, 165–174. doi:10.1007/BF02427376 (1990). 125. Jarrell, M. & Gubernatis, J. Bayesian Inference and the Analytic Continuation of Imaginary-Time Quantum Monte Carlo Data. Physics Reports-Review Section of Physics Letters 269, 133–195. doi:10.1016/0370- 1573(95) 00074-7 (1996). 126. Habershon, S., Braams, B. J. & Manolopoulos, D. E. Quantum Mechanical Correlation Functions, Maximum Entropy Analytic Continuation, and Ring Polymer Molecular Dynamics. Journal of Chemical Physics 127. doi:10. 1063/1.2786451 (2007). 127. Rabani, E., Krilov, G. & Berne, B. Quantum Mechanical Canonical Rate Theory: A New Approach Based on the Reactive Flux and Numerical Analytic Continuation Methods. Journal of Chemical Physics 112, 2605–2614. doi:10.1063/1.480834 (2000). 128. Golosov, A., Reichman, D. & Rabani, E. Analytic Continuation for Quantum Nonadiabatic Rate Constants. Journal of Chemical Physics 118, 457–460. doi:10.1063/1.1535214 (2003). 129. Miller, P. D. Applied Asymptotic Analysis (American Mathematical Society, 2006). 130. Miller, W. Semiclassical Limit of Quantum-Mechanical Transition-State Theory for Nonseparable Systems. Journal of Chemical Physics 62, 1899–1906. doi:10.1063/1.430676 (1975). 131. Miller, W., Zhao, Y., Ceotto, M. & Yang, S. Quantum Instanton Approximation for Thermal Rate Constants of Chemical Reactions. Journal of Chemical Physics 119, 1329–1342. doi:10.1063/1.1580110 (2003). 132. Richardson, J. O. Microcanonical and Thermal Instanton Rate Theory for Chemical Reactions at All Temperatures. Faraday Discussions 195, 49–67. doi:10.1039/c6fd00119j (2016). 101 133. Chapman, S., Garrett, B. & Miller, W. Semiclassical Transition-State Theory for Nonseparable Systems: Application to Collinear H + H2 Reaction. Journal of Chemical Physics 63, 2710–2716. doi:10.1063/1.431620 (1975). 134. Loehle, A. & Kaestner, J. Calculation of Reaction Rate Constants in the Canonical and Microcanonical Ensemble. Journal of Chemical Theory and Computation 14, 5489–5498. doi:10.1021/acs.jctc.8b00565 (2018). 135. Korol, R., Bou-Rabee, N. & Miller III, T. F. Cayley Modification for Strongly Stable Path-Integral and Ring-Polymer Molecular Dynamics. Journal of Chemical Physics 151. doi:10.1063/1.5120282 (2019). 136. Korol, R., Rosa-Raices, J. L., Bou-Rabee, N. & Miller III, T. F. DimensionFree Path-Integral Molecular Dynamics without Preconditioning. Journal of Chemical Physics 152. doi:10.1063/1.5134810 (2020). 137. Dierckx, P. Curve and Surface Fitting with Splines (Oxford University Press, 1995). 138. Jones, E., Oliphant, T., Peterson, P., et al. SciPy: Open Source Scientific Tools For Python. 2001. 139. Tao, X., Shushkov, P. & Miller III, T. F. Simple Flux-Side Formulation of State-Resolved Thermal Reaction Rates for Ring-Polymer Surface Hopping. Journal of Physical Chemistry A 123, 3013–3020. doi:10.1021/acs.jpca. 9b00877 (2019). 140. Lawrence, J. E. & Manolopoulos, D. E. An Analysis of Isomorphic RPMD in the Golden Rule Limit. Journal of Chemical Physics 151. doi:10.1063/ 1.5138913 (2019). 141. Fang, Y., Barber, V. P., Klippenstein, S. J., Mccoy, A. B. & Lester, M. I. Tunneling Effects in the Unimolecular Decay of (CH3)2 COO Criegee Intermediates to OH Radical Products. Journal of Chemical Physics 146. doi:10.1063/1.4979297 (2017). 142. Weston Jr., R. E., Thanh Lam Nguyen, Stanton, J. F. & Barker, J. R. HO + CO Reaction Rates and H/D Kinetic Isotope Effects: Master Equation Models with Ab Initio SCTST Rate Constants. Journal of Physical Chemistry A 117, 821–835. doi:10.1021/jp311928w (2013). 143. Kaufmann, S., Shuai, Q., Auerbach, D. J., Schwarzer, D. & Wodtke, A. M. Associative Desorption of Hydrogen Isotopologues from Copper Surfaces: Characterization of Two Reaction Mechanisms. Journal of Chemical Physics 148. doi:10.1063/1.5025666 (2018). 144. Predescu, C. & Miller, W. Optimal Choice of Dividing Surface for the Computation of Quantum Reaction Rates. Journal of Physical Chemistry B 109, 6491–6499. doi:10.1021/jp040593q (2005). 145. Wigner, E. Calculation of the Rate of Elementary Association Reactions. Journal of Chemical Physics 5, 720–723. doi:10.1063/1.1750107 (1937). 102 146. Wigner, E. The Transition State Method. Transactions of the Faraday Society 34, 0029–0040. doi:10.1039/tf9383400029 (1938). 147. Herman, M., Bruskin, E. & Berne, B. On Path Integral Monte Carlo Simulations. Journal of Chemical Physics 76, 5150–5155. doi:10.1063/1.442815 (1982). 148. Yang, X. & Zhang, D. H. Dynamical Resonances in the Fluorine Atom Reaction with the Hydrogen Molecule. Accounts of Chemical Research 41, 981–989. doi:10.1021/ar700258g (2008). 149. Kwart, H. Temperature-Dependence of the Primary Kinetic Hydrogen Isotope Effect as A Mechanistic Criterion. Accounts of Chemical Research 15, 401–408. doi:10.1021/ar00084a004 (1982). 150. Ghio, E., Mattera, L., Salvo, C., Tommasini, F. & Valbusa, U. Vibrational Spectrum of H and D on the (0001) graphite Surface From Scattering Experiments. Journal of Chemical Physics 73, 556–561. doi:10.1063/1.439855 (1980). 151. Lepetit, B. & Jackson, B. Sticking of Hydrogen on Supported and Suspended Graphene at Low Temperature. Physical Review Letters 107. doi:10.1103/ PhysRevLett.107.236102 (2011). 152. Sha, X. & Jackson, B. First-Principles Study of the Structural and Energetic Properties of H Atoms on a Graphite (0001) Surface. Surface Science 496, 318–330. doi:10.1016/S0039-6028(01)01602-8 (2002). 153. Zecho, T., Guttler, A., Sha, X., Jackson, B. & Kuppers, J. Adsorption of Hydrogen and Deuterium Atoms on the (0001) Graphite Surface. Journal of Chemical Physics 117, 8486–8492. doi:10.1063/1.1511729 (2002). 154. Paris, A. et al. Kinetic Isotope Effect in the Hydrogenation and Deuteration of Graphene. Advanced Functional Materials 23, 1628–1635. doi:10.1002/ adfm.201202355 (2013). 155. Lozada-Hidalgo, M. et al. Sieving Hydrogen Isotopes through Two-Dimensional Crystals. Science 351, 68–70. doi:10.1126/science.aac9726 (2016). 156. Davidson, E. R. M., Klimes, J., Alfe, D. & Michaelides, A. Cooperative Interplay of van Der Waals Forces and Quantum Nuclear Effects on Adsorption: H at Graphene and at Coronene. ACS Nano 8, 9905–9913. doi:10.1021/ nn505578x (2014). 157. Bonfanti, M., Jackson, B., Hughes, K. H., Burghardt, I. & Martinazzo, R. Quantum Dynamics of Hydrogen Atoms on Graphene. II. Sticking. Journal of Chemical Physics 143. doi:10.1063/1.4931117 (2015). 158. Hornekaer, L. et al. Clustering of Chemisorbed H(D) Atoms on the Graphite (0001) Surface Due to Preferential Sticking. Physical Review Letters 97. doi:10.1103/PhysRevLett.97.186102 (2006). 103 159. Casolo, S., Lovvik, O. M., Martinazzo, R. & Tantardini, G. F. Understanding Adsorption of Hydrogen Atoms on Graphene. Journal of Chemical Physics 130. doi:10.1063/1.3072333 (2009). 160. Buenermann, O., Jiang, H., Dorenkamp, Y., Kandratsenka, A., Janke, S. M., Auerbach, D. J. & Wodtke, A. M. Electron-Hole Pair Excitation Determines the Mechanism of Hydrogen Atom Adsorption. Science 350, 1346–1349. doi:10.1126/science.aad4972 (2015). 161. Jiang, H., Dorenkamp, Y., Krueger, K. & Buenermann, O. Inelastic H and D Atom Scattering from Au(111) as Benchmark for Theory. Journal of Chemical Physics 150. doi:10.1063/1.5094693 (2019). 162. Buenermann, O., Jiang, H., Dorenkamp, Y., Auerbach, D. J. & Wodtke, A. M. An Ultrahigh Vacuum Apparatus for H Atom Scattering from Surfaces. Review of Scientific Instruments 89. doi:10.1063/1.5047674 (2018). 163. Grimme, S., Bannwarth, C. & Shushkov, P. A Robust and Accurate TightBinding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All Spd-Block Elements (Z=1-86). Journal of Chemical Theory and Computation 13, 1989–2009. doi:10.1021/acs.jctc.7b00118 (2017). 164. Wang, Y., Qian, H.-J., Morokuma, K. & Irle, S. Coupled Cluster and Density Functional Theory Calculations of Atomic Hydrogen Chemisorption on Pyrene and Coronene as Model Systems for Graphene Hydrogenation. Journal of Physical Chemistry A 116, 7154–7160. doi:10 . 1021 / jp3023666 (2012). 165. Becke, A. Density-Functional Thermochemistry .IV. A New Dynamical Correlation Functional and Implications for Exact-Exchange Mixing. Journal of Chemical Physics 104, 1040–1046. doi:10.1063/1.470829 (1996). 166. Lee, C., Yang, W. & Parr, R. Development of the Colle-Salvetti CorrelationEnergy Formula into A Functional of the Electron-Density. Physical Review B 37, 785–789. doi:10.1103/PhysRevB.37.785 (1988). 167. Manby, F. et al. Entos: A Quantum Molecular Simulation Package. chemrxiv:7762646.v2 (2019). 168. Markland, T. E. & Ceriotti, M. Nuclear Quantum Effects Enter the Mainstream. Nature Reviews Chemistry 2. doi:10 . 1038 / s41570 - 017 - 0109 (2018). 169. Tao, X., Shushkov, P. & Miller III, T. F. Microcanonical Rates from RingPolymer Molecular Dynamics: Direct-Shooting, Stationary-Phase, and MaximumEntropy Approaches. Journal of Chemical Physics 152, 124117. doi:10 . 1063/1.5144307 (2020). 104 170. Collepardo-Guevara, R., Craig, I. R. & Manolopoulos, D. E. Proton Transfer in a Polar Solvent from Ring Polymer Reaction Rate Theory. Journal of Chemical Physics 128. doi:10.1063/1.2883593 (2008). 171. Andersen, H. Rattle: a Velocity Version of the Shake Algorithm for MolecularDynamics Calculations. Journal of Computational Physics 52, 24–34. doi:10. 1016/0021-9991(83)90014-1 (1983). 172. Andersen, H. Molecular-Dynamics Simulations at Constant Pressure and-or Temperature. Journal of Chemical Physics 72, 2384–2393. doi:10.1063/1. 439486 (1980). 173. Stukowski, A. Visualization and Analysis of Atomistic Simulation Data with OVITO-the Open Visualization Tool. Modelling and Simulation in Materials Science and Engineering 18. doi:10 . 1088 / 0965 - 0393 / 18 / 1 / 015012 (2010). 174. Kelly, K. K., Hirschi, J. S. & Singleton, D. A. Newtonian Kinetic Isotope Effects. Observation, Prediction, and Origin of Heavy-Atom Dynamic Isotope Effects. Journal of the American Chemical Society 131, 8382. doi:10.1021/ ja9031083 (2009). 175. Andujar-De Sanctis, I. L. & Singleton, D. A. Racing Carbon Atoms. Atomic Motion Reaction Coordinates and Structural Effects on Newtonian Kinetic Isotope Effects. Organic Letters 14, 5238–5241. doi:10.1021/ol302389f (2012). 176. Tan, X., Shao, H., Hu, T., Liu, G., Jiang, J. & Jiang, H. High Thermoelectric Performance in Two-Dimensional Graphyne Sheets Predicted by FirstPrinciples Calculations. Physical Chemistry Chemical Physics 17, 22872– 22881. doi:10.1039/c5cp03466c (2015). 177. Tanimura, Y. & Ishizaki, A. Modeling, Calculating, and Analyzing Multidimensional Vibrational Spectroscopies. Accounts of Chemical Research 42, 1270–1279. doi:10.1021/ar9000444 (2009). 178. Grimme, S., Antony, J., Ehrlich, S. & Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFTD) for the 94 Elements H-Pu. Journal of Chemical Physics 132, 154104. doi:10.1063/1.3382344 (2010). 179. Wilson, A. K., Woon, D. E., Peterson, K. A. & Dunning Jr, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. IX. the Atoms Gallium through Krypton. Journal of Chemical Physics 110, 7667–7676. doi:10.1063/1.478678 (1999). 180. Polly, R., Werner, H.-J., Manby, F. R. & Knowles, P. J. Fast Hartree-Fock Theory Using Local Density Fitting Approximations. Molecular Physics 102, 2311–2321. doi:10.1080/0026897042000274801 (2004). 105 181. Marcus, R. A. Unimolecular Dissociations and Free Radical Recombination Reactions. The Journal of Chemical Physics 20, 359–364. doi:10.1063/1. 1700424 (1952). 182. Abbott, H. L. & Harrison, I. Microcanonical Transition State Theory for Activated Gas- Surface Reaction Dynamics: Application to H2 /Cu (111) with Rotation as a Spectator. Journal of Physical Chemistry A 111, 9871–9883. doi:10.1021/jp074038a (2007). 183. Bowker, M. The Role of Precursor States in Adsorption, Surface Reactions and Catalysis. Journal of Physics: Condensed Matter 22, 263002. doi:10. 1088/0953-8984/22/26/263002 (2010). 184. Beyer, T. & Swinehart, D. Algorithm 448: Number of Multiply-Restricted Partitions. Communications of the ACM 16, 379. doi:10.1145/362248. 362275 (1973). 185. Eckart, C. The Penetration of a Potential Barrier by Electrons. Physics Review 35, 1303. doi:10.1103/PhysRev.35.1303 (1930). 186. Johnston, H. S. & Heicklen, J. Tunnelling Corrections for Unsymmetrical Eckart Potential Energy Barriers. Journal of Physical Chemistry 66, 532– 533. doi:10.1021/j100809a040 (1962). 187. Brenner, D. W., Shenderova, O. A., Harrison, J. A., Stuart, S. J., Ni, B. & Sinnott, S. B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. Journal of Physics: Condensed Matter 14, 783. doi:10.1088/0953-8984/14/4/312 (2002). 188. Fornace, M. E., Lee, J., Miyamoto, K., Manby, F. R. & Miller III, T. F. Embedded Mean-Field Theory. Journal of Chemical Theory and Computation 11, 568–580. doi:10.1021/ct5011032 (2015). 189. Miyamoto, K., Miller III, T. F. & Manby, F. R. Fock-Matrix Corrections in Density Functional Theory and Use in Embedded Mean-Field Theory. Journal of Chemical Theory and Computation 12, 5811–5822. doi:10.1021/ acs.jctc.6b00685 (2016). 190. Ding, F., Manby, F. R. & Miller III, T. F. Embedded Mean-Field Theory with Block-Orthogonalized Partitioning. Journal of Chemical Theory and Computation 13, 1605–1615. doi:10.1021/acs.jctc.6b01065 (2017).