MOLECULAR PHOTOCHEMISTRY AT THE NANOSCALE AND DEVELOPMENTS TOWARDS MODELING NONADIABATIC DYNAMICS ON MANY ELECTRONIC STATES By Michael Paul Esch A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Chemistry—Doctor of Philosophy 2020 ABSTRACT MOLECULAR PHOTOCHEMISTRY AT THE NANOSCALE AND DEVELOPMENTS TOWARDS MODELING NONADIABATIC DYNAMICS ON MANY ELECTRONIC STATES By Michael Paul Esch Nonadiabatic molecular dynamics (NAMD) simulation methods are useful for gauging the presence of and characterizing pathways for efficient nonradiative recombination (NRR), which can greatly reduce the efficacy of optoelectronic materials. This utility is herein applied to lead- halide perovskites, which have received great attention for use in solar energy conversion, since they readily absorb visible light and have a low propensity for NRR, to further the understanding of why these materials do not readily undergo NRR. More specifically, we investigate the energetic accessibility of conical intersections (CIs) in an archetypal lead-halide perovskite, CsPbBr3, through simulation of the dynamics of an excited molecule sized cluster model (Cs4PbBr6) and time-independent optimizations of a nanoparticle model (Cs12Pb4Br20) of the CsPbBr3 surface. Using knowledge of the perovskite potential energy surfaces (PESs) gained from NAMD simulations of the electronically excited cluster model, geometries of both cluster and nanoparticle model minimal energy CIs (MECIs) were optimized. Energies of the observed MECIs were corrected for both dynamic electron correlation and spin–orbit coupling and were observed to be greater than the optical band gap of bulk CsPbBr3. This suggests that this particular set of intersections is energetically inaccessible after optical excitation and does not facilitate efficient NRR in this material. Characterization of the PESs surrounding the MECIs suggests that the ionic nature of the bonds in CsPbBr3 could be a contributing factor to the high energy of the observed CIs. Thus, the investigation of other, ionic materials for optoelectronic applications would be a promising direction for future research. Though NAMD methods have been successfully applied to study photophysical processes in a range of fields, such as molecular photochemistry, spectroscopy and materials science, the application of these methods has been limited to processes in which a small number of electronic states are populated when they use properties of the PESs to refine trajectory dynamics. This limitation arises from the immense computational cost of computing a large number of PESs at each simulation time step. Here, we introduce the collapse to a block (TAB) correction to Ehrenfest dynamics and its adaptation for dense manifolds of states (TAB-DMS), which model electronic coherence loss in a state-pairwise fashion. The employed state-pairwise formulation is demonstrated to be necessary for maintaining the most physical description of coherence loss when a model exists in a coherent superposition of three or more electronic states. Additionally, the TAB-DMS adaptation employs the history of the time-dependent Ehrenfest wavefunction to efficiently generate a small number of approximate electronic states, eliminating the need to compute a large number of electronic PESs. We demonstrate through the computation of TAB and TAB-DMS ensemble probabilities of transmission for a series of one-dimensional model problems, that TAB-DMS retains the accuracy of TAB simulations when only a modest number of approximate eigenstates are computed. The accuracy of the TAB-DMS methodology is also demonstrated to be systematically improvable through the computation of additional approximate electronic states. These results indicate pairing TAB-DMS with ab-initio electronic structure theory and efficient GPU propagation schemes would be fruitful next steps towards modeling real photophysical processes occurring over many electronic states. Harold Esch (1933-2015) and Mary Schneider (1931-2018) In loving memory of iv ACKNOWLEDGMENTS I would first like to acknowledge the tremendous support from my family, especially my wife Laura. The research presented here would not have been possible without their love and encouragement. I would like to acknowledge the outstanding instructional skills of my advisor, Benjamin G. Levine, with emphasis on computational method development and scientific writing and communication. His efforts have immensely polished my skills as a computational researcher. I would also like to thank him for his compassion as a mentor. During my time in his lab, I have lost two close members of my family. Their deaths had put an immense strain on me, and I am grateful for Benjamin’s understanding and patience during my recovery. For serving as my committee, I would like to thank Robert I. Cukier, Katharine C. Hunt and Thomas W. Hamann for their guidance and scientific insights. I would also like to thank both our current and former group members, Yinan, Garrett, Scott, Wei-tao, Dmitry, Brandon, Dylan, Andy, Fangchun, Adam and Katy for all of their support, advice, and the opportunities for collaboration. I would lastly like to thank my undergraduate advisor Christopher Lawrence. I appreciate the time he took to introduce me to computational chemistry and teach me to code, and all the advice he gave me. v TABLE OF CONTENTS LIST OF TABLES ....................................................................................................................... viii LIST OF FIGURES ..................................................................................................................... xiii CHAPTER 1: INTRODUCTION ....................................................................................................1 1.1 A Molecular Photochemistry Perspective on NRR in Semiconductor Nanomaterials ................................................................................................................2 1.2 Decoherence Corrections to Ehrenfest Dynamics for Simulations Incorporating Many Electronic States ..................................................................................................4 CHAPTER 2: A CONICAL INTERSECTION PERSPECTIVE ON THE LOW NONRADIATIVE RECOMBINATION RATE IN LEAD HALIDE PEROVSKITES ........................................................................................................9 2.1 Abstract ..........................................................................................................................9 2.2 Introduction ..................................................................................................................10 2.3 Computational Methods ...............................................................................................12 2.3.1 Cluster Model Ab Initio Multiple Spawning Simulations ............................12 2.3.2 Quantum Chemical Study of Cluster Model CIs ..........................................14 2.3.3 Optimization of Nanoparticle MECIs ...........................................................16 2.4 Results and Discussion ................................................................................................16 2.4.1 Photodynamics of the Cluster Model ............................................................16 2.4.2 Analysis of Cluster Model PESs ...................................................................19 2.4.3 FOMO-CASCI PES Validation and Active Space Selection for NP Model ......................................................................................................31 2.4.4 Analysis of NP Model PESs .........................................................................34 2.5 Conclusions ..................................................................................................................40 2.6 Acknowledgments........................................................................................................41 CHAPTER 3: STATE-PAIRWISE DECOHERENCE TIMES FOR NONADIABATIC DYNAMICS ON MORE THAN TWO ELECTRONIC STATES ........................42 3.1 Abstract ........................................................................................................................42 3.2 Introduction ..................................................................................................................43 3.3 Computational Methods ...............................................................................................46 3.3.1 DCE Algorithm .............................................................................................46 3.3.2 Decoherence Corrections ..............................................................................48 3.3.2.1 Collapse To A State (TAS) ............................................................48 3.3.2.2 Independent Pairwise Stochastic Decoherence (IPSD) .................50 3.3.2.3 Collapse To A Block (TAB) ..........................................................52 3.3.3 Model Problems and Initial Conditions ........................................................56 3.4 Results and Discussion ................................................................................................60 3.4.1 Two-State Model ..........................................................................................60 3.4.2 Three-State Model ........................................................................................61 vi 3.4.3 Five-State Model ...........................................................................................64 3.5 Conclusions ..................................................................................................................66 3.6 Acknowledgments........................................................................................................67 CHAPTER 4: DECOHERENCE-CORRECTED EHRENFEST DYNAMICS ON MANY ELECTRONIC STATES ........................................................................................68 4.1 Abstract ........................................................................................................................68 4.2 Introduction ..................................................................................................................69 4.3 Computational Methods ...............................................................................................71 4.3.1 Ehrenfest Propagation Scheme .....................................................................71 4.3.2 Decoherence Corrections ..............................................................................72 4.3.2.1 To-a-Block (TAB) Decoherence Correction .................................73 4.3.2.2 Approximate Electronic Eigenstates for Dense Manifolds of States (DMS) .............................................................................74 4.3.2.3 Wave Function Collapse Algorithm ..............................................76 4.3.3 Model Problems and Initial Conditions ........................................................79 4.4 Results and Discussion ................................................................................................83 4.4.1 Convergence of TAB-DMS with Respect to Number of Approximate Eigenstates ....................................................................................................83 4.4.2 Comparison to Exact Quantum Dynamics and Width Parameter Sensitivity .....................................................................................................86 4.5 Conclusions ..................................................................................................................91 4.6 Acknowledgments........................................................................................................92 CHAPTER 5: CONCLUSIONS AND OUTLOOK ......................................................................93 APPENDIX ....................................................................................................................................99 BIBLIOGRAPHY ........................................................................................................................151 vii LIST OF TABLES Table 2.1: Optimized cluster model MECI and local S1 minimum energies. ................................26 Table 2.2: Largest singlet–triplet SOCs calculated at the FC point and optimized MECI geometries, and the largest singlet–triplet SOCs including the intersecting electronic states (S0 and S1). .......31 Table 2.3: Energies of the NP model at optimized geometries. Corrections are computed as described in the text. ......................................................................................................................36 Table 3.1: The slopes, y-intercepts, and initial populations of the two-state model......................58 Table 3.2: The slopes, y-intercepts, and initial populations of the three-state model....................58 Table 3.3: The slopes, y-intercepts, and initial populations of the five-state model. ....................58 Table 3.4: Fractions of trajectories that have collapsed to individual electronic states at the end of the simulations of the three-state model. We define a trajectory to have collapsed to a state if it its final population on that state is greater than 0.98. A method that properly treats the coherence between parallel states 1 and 2 would collapse with the probabilities reported in the target column. ................................................................................................................................64 Table 4.1: Number of diabatic states in each one-dimensional model, the energy spacing between positively sloped states and the values of the band splitting parameters. ......................................81 Table A.1: Absolute energies of singlet electronic states calculated at the B3YLP optimized FC geometry. The singlet state energies were calculated at the EOM-CCSD/SBKJC Polarized(p, 2d), SA5-CASSCF(2/5)/SBKJC VDZ ECP, MS5-CASPT2-c19/SBKJC Polarized(p, 2d), MS5- CASPT2-c0/SBKJC Polarized(p, 2d), MRCI/SBKJC Polarized(p, 2d) with the same active space and number of electronic states and the Davidson correction, TD-PBE0/SBKJC VDZ ECP and FOMO-CASCI(2/5)/SBKJC VDZ ECP levels of theory. ............................................................105 Table A.2: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the B3LYP optimized FC geometry. .................................................................................................107 Table A.3: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the B3YLP optimized FC geometry in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. .............................................................................................................................107 Table A.4: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized S1 minimum local to the FC region calculated at the same level of theory as the optimization. ..........................................................................................................................108 viii Table A.5: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized S1 minimum local to the FC region. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)- c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory........................................................................................109 Table A.6: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. ..........109 Table A.7: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI. .........................................................111 Table A.8: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ......................................................................................111 Table A.9: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized BrD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ..........................................................................................................................................112 Table A.10: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. ..........113 Table A.11: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized BrD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ............................................................................................................................114 Table A.12: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. ..........115 Table A.13: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI. .........................................................116 Table A.14: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ......................................................................................117 ix Table A.15: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ..........................................................................................................................................117 Table A.16: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. ..........119 Table A.17: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD PES region S1 minimum calculated at the same level of theory as the optimization. ................................................................................................................................120 Table A.18: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ............................................................................................................................121 Table A.19: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. ........121 Table A.20: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI. .......................................................123 Table A.21: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ......................................................................................123 Table A.22: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized DBrS CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ..........................................................................................................................................124 Table A.23: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization. ......125 Table A.24: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI. .....................................................126 x Table A.25: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ....................................................................127 Table A.26: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsBrD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ..........................................................................................................................................127 Table A.27: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization. ......129 Table A.28: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD PES region S1 minimum calculated at the same level of theory as the optimization. ................................................................................................................................130 Table A.29: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsBrD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ............................................................................................................................130 Table A.30: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the same level of theory as the optimization. ...........131 Table A.31: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI. ..........................................................133 Table A.32: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ......................................................................................133 Table A.33: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized Bal CI. Energies are calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ..........................................................................................................................................134 Table A.34: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal PES region S1 minimum calculated at the same level of theory as the optimization. ................................................................................................................................135 xi Table A.35: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized Bal PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. ............................................................................................................................136 Table A.36: Absolute energies of the singlet electronic states at the PBE0 optimized ground state minimum energy geometry calculated at the TD-PBE0/SBKJC VDZ ECP and FOMO- CASCI(6/11)/SBKJC VDZ ECP levels of theory. ......................................................................137 Table A.37: Absolute energy of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP FC region S1 minimum, calculated at the same level of theory as the optimization. ................................................................................................................................139 Table A.38: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. ..........141 Table A.39: Absolute energy of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP BrD region S1 minimum, calculated at the same level of theory as the optimization. ................................................................................................................................143 Table A.40: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. ..........145 Table A.41: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. ........147 xii LIST OF FIGURES Figure 2.1: Cs4PbBr6 cluster model at the Franck-Condon (FC) point optimized at the B3LYP/SBKJC Polarized(p, 2d) level (left) and Cs12Pb4Br20 NP FC point optimized at the PBE0/SBKJC VDZ level (right). Cs, Br, and Pb atoms are represented in green, purple, and black, respectively. All molecular images in this work were created with MacMolPlt11. ......................12 Figure 2.2: The S1–S0 energy gap of the twenty initial TBFs as a function of time drawn from AIMS simulations of the Cs4PbBr6 cluster model. ........................................................................17 Figure 2.3: Thick lines show the populations of S0 (red) and S1 (dark blue) as well as the aggregate S2–S4 populations (orange) averaged over all AIMS calculations. The populations of individual simulations are represented by thin lines, with S0, S1 and aggregate S2–S4 populations in pink, light blue, and light orange, respectively. ..............................................................................................18 Figure 2.4: The positions of all TBFs (initial and spawned) projected onto two coordinates: the largest Br–Pb distance and the largest Cs–Pb distance. Geometries at which these TBFs spawned new basis functions to the S0 electronic state are marked by pink circles. ....................................19 Figure 2.5: The fractionally occupied SANOs at the cluster model FC point, computed at the SA- CASSCF level of theory, and generated with an isovalue of 0.01. All subsequent SANOs displayed have also been generated with an isovalue of 0.01........................................................22 Figure 2.6: Optimized Cs4PbBr6 cluster model S1–S0 MECI geometries. MECI geometries are organized by PES region, with each geometry labeled with the level of theory at which the optimization was carried out. .........................................................................................................24 Figure 2.7: Optimized cluster model MECI and local S1 minimum energies. Energy calculations were performed at the SA-CASSCF (solid black lines), MS-CASPT2-c19 (solid red lines), MS- CASPT2-c0 (dashed orange lines) MRCI, (dashed purple lines), and the FOMO-CASCI (dotted blue lines) levels of theory, as described above. The energy calculations performed at the MS- CASPT2-c0 and MRCI levels of theory were performed at geometries optimized at the MS- CASPT2-c19 level of theory. All other calculations were performed at geometries that were optimized at their respective levels of theory. ...............................................................................25 Figure 2.8: a) Scan of the S0 (black) and S1 (red) PESs of the cluster model as the dissociated (most distant) Br atom at the BrD CI is displaced along the vector connecting it with the Pb center. All other degrees of freedom remain frozen. The PES is computed at the SA-CASSCF level. b) The S0 (black) and S1 (red) Mulliken charges of the dissociated Br atom at the same geometries as in the PES in a)...................................................................................................................................28 the SA-CASSCF-optimized BrD MECI Figure 2.9: Fractionally occupied SANOs at geometry. .......................................................................................................................................29 xiii Figure 2.10: Fractionally occupied SANOs of the cluster model computed at the FOMO-CASCI level of theory at the FC geometry. ...............................................................................................32 Figure 2.11: Fractionally occupied SANOs computed at the FOMO-CASCI(6/11)/SBKJC VDZ ECP level of theory at the NP FC geometry. .................................................................................34 Figure 2.12: Cs12Pb4Br20 nanoparticle optimized geometries. MECI optimizations were performed at the FOMO-CASCI level of theory, while the FC geometry was optimized at the PBE0 level. ...............................................................................................................................................35 Figure 2.13: Energies of the NP model at optimized geometries. Corrections are computed as described in the text. The experimental optical band gap of CsPbBr3 (purple dashed line) was taken from reference 125. ................................................................................................................36 Figure 2.14: Fractionally occupied SANOs averaged over the S0 and S1 electronic states computed at the FOMO-CASCI level of theory at the NP optimized DBrS CI. ............................................38 Figure 2.15: a) Scan of the PES surrounding the NP model BrD CI. The dissociated Br atom is displaced along the vector between it and the closest Pb center while all other degrees of freedom are frozen. Energies are calculated at the FOMO-CASCI level of theory. b) The S0 and S1 Mulliken charges on the dissociated Br atom as a function of the same coordinate. ....................................39 Figure 2.16: Fractionally occupied SANOs calculated at the NP BrD MECI geometry. Orbitals are computed at the FOMO-CASCI level as described in the text. ...............................................40 Figure 3.1: The adiabatic PESs of the a) two-state, b) three-state, and c) five-state models. ........59 Figure 3.2: The ensemble averaged absolute coherences as a function of time for the two-state model as computed by TAS (dashed red), IPSD (blue dash/dot) and TAB (light blue dots). The target exponential decay is shown with a solid black line for comparison. ...................................61 Figure 3.3: The ensemble averaged absolute coherence between states 1 and 2 of the three-state model as a function of time computed via TAS (red line), IPSD (dark blue line), and TAB (light blue line). .......................................................................................................................................62 Figure 3.4: Ensemble averaged absolute coherences for the five-state model computed with TAB[(a)–(d)], TAS[(e)–(h)], and IPSD [bottom, (i)–(l)]. Each column of panels corresponds to the coupling between state pairs with equal force differences; in (a), (e), and (i), in (b), (f), and (j), etc. Target exponential decays in the absolute coherences [as computed with Equation (3.29)] are shown as solid black lines. ...........................................................................................65 Figure 4.1: Adiabatic PES crossing regions for model 9_10000 (panel a) and model 9_10000_split (panel b). ........................................................................................................................................81 Figure 4.2: Adiabatic PES crossing regions for model models 17_5000 (panel a), model 17_2500 (panel b), model 17_1250 (panel c) and model 17_625 (panel d). ................................................82 xiv ,1ii−,2ii− Figure 4.3: TAB-DMS transmission probability (solid black dots) for models 9_10000 (panel a), 9_10000_split (panel b), 17_2500 (panel c), and 17_1250 (panel d) shown as a function of the number of approximate eigenstates. Solid red lines correspond to full TAB calculations. The width parameter α is set to 216.0 bohr-2 for all simulations.....................................................................85 Figure 4.4: TAB ensemble (blue dots), Ehrenfest (dashed red line), and numerically exact (solid black line) transmission probabilities for models 9_10000 (panel a), 9_10000_split (panel b), 17_2500 (panel c), and 17_1250 (panel d). TAB transmission probabilities with intuitive parameterization are highlighted as light blue triangles. ...............................................................87 Figure 4.5: Illustration comparing exponential (red dashed) and Gaussian (black) approximations to the absolute quantum coherence as a function of time, with equivalent time constants. ..........89 Figure 4.6: Error in the probabilities of transmission predicted by traditional Ehrenfest trajectories (black dots), intuitively parameterized TAB ensembles (light blue triangles), TAB ensembles with an α parameter of 24.0 bohr-2 (dark blue dots), and TAB ensembles with an α parameter of 216.0 bohr-2 (red dots), shown as a function of the δ parameter of each 17 state model. ........................90 Figure A.1: B3LYP/SBKJC Polarized(p, 2d) optimized ground state minimum energy geometry. The absolute energy of this structure at the optimization level of theory is -84.9289162 Hartree. ...................................................................................................................105 Figure A.2: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory. .......................................................................106 Figure A.3: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5- CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..............................................................106 Figure A.4: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO- CASCI(2/5)/SBKJC VDZ ECP level of theory. ..........................................................................106 Figure A.5: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the B3LYP optimized FC geometry calculated at the SA5-CASSCF(2/5)/cc- PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. ..................................107 Figure A.6: Local S1 minimum to the FC geometry optimized at the SA5-CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. ....................................................108 Figure A.7: Local S1 minimum to the FC geometry optimized at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. .......................108 the B3YLP optimized FC geometry calculated at the B3LYP optimized FC geometry calculated at the B3LYP optimized FC geometry calculated at xv Figure A.8: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................109 Figure A.9: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization......................................................................................110 Figure A.10: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. ..........................................................................................................................................110 Figure A.11: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ........................................................................................................111 Figure A.12: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized BrD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..........................112 Figure A.13: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................113 Figure A.14: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization......................................................................................113 Figure A.15: Local S1 minimum in the cluster model BrD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) the CIOpt software package. .......................................................................................................................................114 Figure A.16: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................115 Figure A.17: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization......................................................................................116 Figure A.18: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. ............................................................................................................................................116 theory using level of xvi Figure A.19: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ........................................................................................................117 Figure A.20: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..........................118 Figure A.21: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................118 Figure A.22: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization......................................................................................119 Figure A.23: Local S1 minimum in the cluster model CsD PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. ...............120 Figure A.24: Local S1 minimum in the cluster model CsD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) the CIOpt software package. .......................................................................................................................................120 Figure A.25: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the DBrS region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................121 Figure A.26: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization......................................................................................122 Figure A.27: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. ............................................................................................................................................122 Figure A.28: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the DBrS region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ........................................................................................................123 Figure A.29: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized DBrS CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..........................124 theory using level of xvii Figure A.30: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ..............................................................................................................125 Figure A.31: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization......................................................................................126 Figure A.32: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. ............................................................................................................................................126 Figure A.33: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ........................................................................................................127 Figure A.34: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsBrD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..........................128 Figure A.35: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. ..............................................................................................................128 Figure A.36: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP CsBrD CI calculated at the same level of theory as the optimization. ..............................................................................................129 Figure A.37: Local S1 minimum in the cluster model CsBrD PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. ...............129 Figure A.38: Local S1 minimum in the cluster model CsBrD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) the CIOpt software package. .......................................................................................................................................130 Figure A.39: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the Bal region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................131 Figure A.40: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the same level of theory as the optimization......................................................................................132 theory using level of xviii level of theory using Figure A.41: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. ............................................................................................................................................132 Figure A.42: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the Bal region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................133 Figure A.43: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized Bal CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. ..........................134 Figure A.44: Local S1 minimum in the cluster model Bal PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. ...............135 Figure A.45: Local S1 minimum in the cluster model Bal PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) the CIOpt software package. .......................................................................................................................................136 Figure A.46: Ground state minimum energy geometry of the NP model optimized at the PBE0/SBKJC VDZ ECP level of theory. ....................................................................................137 Figure A.47: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the PBE0 ground state minimum energy geometry calculated at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory. ........................................................................139 Figure A.48: Local S1 minimum in the NP model FC PES region optimized at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory using the CIOpt software package..................139 Figure A.49: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the BrD region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................141 Figure A.50: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. ...............................................................................143 Figure A.51: Local S1 minimum in the NP model BrD PES region optimized at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory using the CIOpt software package..................143 Figure A.52: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the CsD region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................145 xix Figure A.53: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. ...............................................................................147 Figure A.54: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the DBrS region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. .........................................................................................................................147 Figure A.55: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. ...............................................................................149 Figure A.56: a) Cluster model SA-CASSCF PES scan as the dissociated Cs atom at the CsD CI is displaced along the vector between it and the Pb center while all other degrees of freedom are frozen. Energies are shifted such that the S0 energy at the FC point is zero. b) Cluster model Mulliken charges of the dissociated Cs atom as it is displaced along the dissociated Cs–Pb vector around the SA-CASSCF optimized CsD CI with all other degrees of freedom frozen. ..............149 Figure A.57: a) NP model FOMO-CASCI PES scan as the dissociated Cs atom at the CsD CI is displaced along the vector between itself and the average position of the four Pb atoms while all other degrees of freedom are frozen. Energies are shifted such that the S0 at the FC point is zero. b) NP model Mulliken charges of the dissociated Cs atom as it is displaced along the vector between it and the average position of the four Pb atoms around the FOMO-CASCI optimized CsD CI with all other degrees of freedom frozen. .......................................................................150 xx CHAPTER 1: INTRODUCTION The presence of efficient pathways for nonradiative recombination (NRR), the fundamental process of converting electronic energy to heat, often determines the efficacy of molecules and materials for use in photophysical applications. For example, a material with long-lived electronic excitations would be better suited to optoelectronic applications, such as photovoltaic energy conversion and light emission, than a material which readily undergoes NRR after light absorption. In contrast, devices designed for rapid heat generation through photothermal energy conversion would employ materials which readily undergo efficient NRR. Determination of excitation lifetimes is possible with sophisticated experimental means developed over the past several decades, however, elucidating the reasons why a particular material does, or does not, readily undergo NRR remains difficult. Computation can provide additional insight when the reasons for, or the barriers to, efficient NRR are unknown. These insights can be provided through either time-independent, or time-dependent simulations. In time-independent computations, electronic potential energy surfaces (PESs) are calculated for fixed sets of nuclear coordinates. This formulation is particularly useful for numerical optimizations to determine quantities such as the energetic accessibility of pathways for efficient NRR, which only depend on the characteristics of the electronic PESs. At an additional computational cost, time-dependent modeling of coupled electronic and nuclear motion through nonadiabatic molecular dynamics simulations (NAMD) can be used to obtain an atomistic picture of nuclear motion and timescales of electronic transitions in photophysical processes. Being able to obtain these observables can also greatly benefit fields where obtaining experimental data is quite difficult, such as atmospheric chemistry. However, NAMD methods 1 which compute electronic PESs for trajectory propagation are currently limited to model processes involving only a small number of electronic states, and quickly become intractable as the number of included states increases. In this work, we will describe our efforts to investigate the energetic accessibility of pathways of NRR in lead-halide perovskite nanomaterials and our development of a novel NAMD method, designed for time-dependent simulations involving many electronic states. 1.1 A Molecular Photochemistry Perspective on NRR in Semiconductor Nanomaterials Over the past decade, lead-halide perovskites have received a tremendous amount of attention, and are particularly promising to incorporate into optoelectronic devices for photovoltaic energy conversion, due to the substantial, rapid increase of reported device efficiencies.2-4 A large portion of this attention was gained due to the low propensity of lead-halide perovskites to undergo NRR.5, 6 However, the question of why lead halide perovskites have a low propensity for NRR is still an active area of research.5-12 A structural understanding of why the rate of NRR is limited in lead-halide perovskites could inform the design of yet more efficient materials. In Chapter 2, we present our work to answer this question for an archetypal lead-halide perovskite, CsPbBr3, where the investigation performed relied heavily on principles drawn from molecular photochemistry. It is widely accepted in molecular photochemistry that conical intersections (CIs), which are seams of degeneracy between two electronic PESs, can provide pathways for efficient nonradiative decay when energetically accessible after electronic excitation.13-17 Additionally, characterization of CIs, associated with specific surface structures and defects, such as dangling bonds in luminescent silicon nanostructures, has provided great insight towards semiconductor photodynamics.18-25 Our research group has also argued that the study of CIs is more generally 2 applicable to other semiconductor nanomaterials, based on the proposal that CIs are local features of electronic PESs.26 In this previous work, the two conditions used to define CI locality were that the associated electronic transition occurs through the transfer of an electron from one spatially local orbital to another, and the degeneracy of the involved electronic states is brought about by local geometric distortions.26 These conditions have been met in models of silicon nanocrystals18- 20, silicon surfaces with defects on next-nearest neighbor silicon atoms27, and linear polyenes26. Knowing the viability of lead-halide perovskites for optoelectronic applications, we hypothesized that the lowest energy points on CI seams, or minimal energy conical intersections (MECIs), were energetically inaccessible after optical excitation in CsPbBr3. The first steps of our investigation, fully described in Chapter 2, entailed NAMD simulations of a cluster model of the CsPbBr3 surface, and subsequent optimization of model MECI geometries. The successful use of cluster models to gauge energies of bulk material and large nanoparticle MECIs is rooted in their locality.19, 26, 27 More specifically, when the geometric and electronic distortions associated with a particular MECI are completely captured by a model, the relative energy of the MECI to the ground state geometry of the system will not change with increasing model size, effectively representing an extrapolation to the bulk limit. However, it is impossible to know whether a given cluster model adequately captures all contributing effects without performing additional computations with larger models. At the nanoscale, the types of tractable electronic structure calculations have previously been greatly limited due to the poor scaling of computational cost with respect to system size, but creative hardware adaptation and clever methodological developments have helped to increase the repertoire. Regarding hardware, graphical processing units (GPUs) have been leveraged by electronic structure software packages for more efficient computations.28-33 The development of 3 efficient GPUs has largely taken place in the video game industry over the past several decades, from the increased demand for improved graphics and more complex physics in modern games. More specific to efficiently modeling electronic excited states of nanoscale systems, is the development of size-intensive two-step complete active space configuration interaction (CASCI) methods. Two-step CASCI methods capture the effects of static electron correlation, similar to state-averaged multiconfigurational self-consistent field (SA-MC-SCF) methods, which are robust for, and routinely used, to model photophysical processes in molecule sized systems. However, SA-MC-SCF methodologies do not provide a size intensive description of electronic excitations and can become numerically unstable at the nanoscale.34, 35 Size-intensive two-step CASCI methods do not suffer from numerical instability though, and are often less expensive than comparable SA-MC-SCF methods. Thus, two-step CASCI methods provide a promising alternative as efficient means of nanoscale electronic structure computations. 1.2 Decoherence Corrections to Ehrenfest Dynamics for Simulations Incorporating Many Electronic States Despite the wealth of information time-independent calculations can provide, data exclusively available from NAMD simulations, such as directly computed electronic excited state lifetimes, are also important for advancing the fields of materials science, molecular photochemistry, and spectroscopy. One such example, pertaining to methylammonium lead-iodide perovskites, was a study that demonstrated, through the use of NAMD simulations, that doping the perovskite with chlorine impurities extends the excited state lifetime.12 However, the computational tractability of NAMD simulations that use knowledge of the electronic PESs in trajectory propagation are currently limited to a few states, even when employment of GPUs and 4 efficient electronic structure methodologies are considered. This is due to the vast computational cost arising from computing properties of a large number of electronic PESs at every trajectory time step, for hundreds, or potentially thousands, of trajectories. Though many interesting photophysical processes that involve only a handful of electronic states can be modeled with this class of NAMD method, a systematic means of applying the refinements in them to dynamics that occur over many electronic states would have widespread benefits in areas such as strong field physics, the study of plasmonic materials, and the modeling of radiation damage to both molecules and materials. To avoid the computational cost of characterizing electronic PESs altogether, one could turn to Ehrenfest (mean-field) dynamics.36-39 Ehrenfest dynamics is an independent trajectory (IT) method, which replaces the nuclear wavefunction by a swarm of non-interacting classical trajectories that propagate on a quantum mechanically defined PES. In Ehrenfest dynamics, this PES is the expectation value of the electronic Hamiltonian bracketed by the time-dependent electronic wavefunction. This definition of the quantum PES circumvents the need to explicitly determine the electronic states and avoids the computational cost associated with their determination. Additionally, this definition allows trajectories to smoothly transition between electronic states in regions of strong coupling. The IT formalism also allows for straightforward interpretation of simulation results, due to the classical representation of the nuclei. However, there are well-known qualitative errors associated with propagating on the mean-field PES of Ehrenfest dynamics.40 These errors arise because trajectory propagation does not account for the loss of coherence between populations on nonparallel electronic PESs.41-44 This is particularly problematic for dynamics involving the traversal of multiple PES crossing regions, which are likely to occur after high energy electronic excitations.45 5 To account for the loss of coherence between electronic populations on nonparallel PESs while preserving the IT dynamics formulation, decoherence corrections have been applied to both Ehrenfest dynamics and the popular fewest switches surface hopping algorithm, which similarly suffers from this issue. Decoherence corrections destroy electronic coherences by forcing individual trajectories towards propagating on a single electronic PES over time.45-62 As a result, ensemble electronic coherences decay at a rate determined by the applied correction, while the ensemble averaged electronic populations remain unchanged. Decoherence corrected methodologies have proven themselves robust for modeling photophysical processes over a small number of electronic states, but they require characterization of electronic PESs in order to be applied to trajectory dynamics. Thus, limiting the photophysical processes they can be employed to simulate. Two areas that would benefit from computationally inexpensive decoherence corrections to Ehrenfest dynamics, in addition to those mentioned above, are ultrafast spectroscopy and atmospheric chemistry. In ultrafast spectroscopy, the use of shaped pulses, can induce high energy excitations in the target that would fall outside the range of states considered in routine decoherence-corrected Ehrenfest (DCE) dynamics. Even if the excitation lifetime is short, the dynamics of the target immediately after the excitation would be influenced by the corresponding PESs, and not considering them in DCE simulations would lessen the correspondence between the simulations and the experiment. In atmospheric chemistry, there are also many interesting photophysical processes involving high energy electronic excitations, such as electron ejection and recombination with its parent ion, that similarly fall outside the range of states routinely considered in DCE dynamics. However, tractable NAMD simulations are even more vital to advancing the 6 field of atmospheric chemistry due to the potentially great difficulty of obtaining experimental data. One conceptually simple solution is to provide a smaller number of efficiently calculated approximate electronic PESs to a decoherence correction to refine the Ehrenfest propagation. However, computation of an approximate electronic basis that can be employed by a decoherence correction without loss of accuracy is not trivial. Additionally, the decoherence correction used in trajectory propagation must not require knowledge of all electronic PESs to be applied to the Ehrenfest dynamics, must have all necessary values be readily attainable from data structures available during Ehrenfest propagation, and describe electronic coherence loss in a state-pairwise fashion. A state-pairwise description of electronic coherence loss is one where the rate two electronic populations lose coherence is exclusively dependent on properties of the involved PESs. This rate will be referred to as the state-pairwise decoherence time. In IT dynamics, when state- pairwise decoherence times are computed to refine the coherent propagation in Ehrenfest dynamics or trajectory surface hopping, they are often approximated to depend on the difference in the instantaneous gradients of the involved PESs, as is done in this work.43, 44, 56 Conveniently, the number of state-pairwise decoherence times, , where N is the total number of electronic states, is always equal to the number of independent coherences in the electronic density matrix. This suggests that a state-pairwise description of coherence loss will always have a sufficient number of independent variables for accurate treatment. In Chapter 3, we describe our developed collapse to a block (TAB) correction to Ehrenfest dynamics, which employs a state- pairwise description of coherence loss and meets the other criteria, specified above, for using approximate electronic states to refine Ehrenfest propagation. Chapter 4 details the development of the algorithm used for generating an approximate electronic basis that can be used by the TAB 7 ()1/2NN− decoherence correction, or other compatible refinements to Ehrenfest dynamics, in simulations involving dense manifolds of electronic states, without loss of accuracy. 8 CHAPTER 2: A CONICAL INTERSECTION PERSPECTIVE ON THE LOW NONRADIATIVE RECOMBINATION RATE IN LEAD HALIDE PEROVSKITES Reproduced with permission from “A Conical Intersection Perspective on the Low Nonradiative Recombination Rate in Lead Halide Perovskites”, Michael P. Esch, Yinan Shu and Benjamin G. Levine, J. Phys. Chem. A 2019, 123, 2661-2673. Copyright 2019 American Chemical Society. 2.1 Abstract The utility of optoelectronic materials can be greatly reduced by the presence of efficient pathways for nonradiative recombination (NRR). Lead halide perovskites have garnered much attention in recent years as materials for solar energy conversion because they readily absorb visible light, are easy to synthesize, and have a low propensity for NRR. Here we report a theoretical study of the pathways for NRR in an archetypal lead halide perovskite: CsPbBr3. Specifically, we have identified a set of conical intersection (CIs) in both a molecule-sized cluster model (Cs4PbBr6) and nanoparticle model (Cs12Pb4Br20) of the CsPbBr3 surface. The energies of the minimal energy CIs, corrected for both dynamical electron correlation and spin–orbit coupling, are well above the bulk band gap of CsPbBr3, suggesting that these intersections do not provide efficient pathways for NRR in this material. Analysis of the electronic structure at these intersections suggests that the ionic nature of the bonds in CsPbBr3 may play a role in the high energy of these CIs. The lowest-energy intersections all involve charge transfer over long distances, whether it be across a dissociated bond or between neighboring unit cells. 9 2.2 Introduction The function of an optoelectronic material is often dependent on its ability to maintain long-lived electronic excitations. Nonradiative recombination (NRR), which is the fundamental process by which electronic energy is converted to heat, can greatly reduce the efficiencies of optoelectronic materials for photovoltaic energy conversion, light emission, and photocatalysis. Having garnered much attention for their low rate of NRR,63, 64 lead halide perovskites have emerged as a promising new class of optoelectronic materials for these applications. They are particularly promising for photovoltaic conversion, where reported efficiencies have increased rapidly from 3.8%2 in 2009 to over 20%3, 65 in more recent years. In addition, colloidal CsPbX3 (X = Cl, Br, I) perovskite nanocrystals have bright emission, and modification of the halide composition is a convenient means of tuning their emission energies.66 Recently much attention has been focused on increasing the stability of these materials67-71 and circumventing the need for toxic heavy metals.72 There are still fundamental questions, even about pure lead halide perovskites, that remain unanswered, however. Among the most important is, “why do lead halide perovskites have such a low propensity for NRR?”63, 64, 73 A deep understanding of how the structure and composition of these materials limits their nonradiative rates could inspire the design of yet more efficient materials for optoelectronic applications. Nonradiative recombination processes are inherently difficult to study, however, especially in materials as structurally complex as lead halide perovskites. Simulation is well-suited to address these challenges, because it offers an approximate look at coupled electron–nuclear dynamics. Computational studies have now provided several explanations for the low rate of NRR. It has been suggested that the coupling of nuclear motion to charge carrier motion in these materials has a strong influence on recombination 10 kinetics.74, 75 Others have suggested that the most common defects do not introduce deep levels in the band gap.10, 11 Doping strategies have been proposed to increase the excited state lifetime.12, 76 Sophisticated time-dependent calculations have enabled the assignment of transient absorption bands to specific electronic motions that precede recombination.77 In this work, we approach this question from a different theoretical angle, drawing inspiration from the field of molecular photochemistry. Our group has recently argued that conical intersections (CIs) associated with specific surface structures and defects can introduce efficient pathways for NRR in semiconductors.78 Conical intersections are points of degeneracy between potential energy surfaces (PESs) that allow for rapid transitions between the intersecting states. Though very challenging to observe experimentally,79 CIs are now widely accepted to have predictive value in the study of molecular photochemistry.13, 15-17, 80 Recent theoretical and computational advances have enabled the identification and characterization of conical intersections in full semiconductor nanoparticles. Already, the study of conical intersections in materials has yielded numerous insights into the photophysics of luminescent silicon nanostructures.18, 19, 21, 78, 81-84 In this work we characterize CIs between the ground and first excited states of lead halide perovskites, with the goal of providing insight into their low rates of NRR. First, ab initio nonadiabatic molecular dynamics (NAMD) simulations were performed to model the photodynamics of a molecule-sized cluster model (Cs4PbBr6). High-level quantum chemical calculations were then used to characterize the PESs explored during the dynamics simulations. Lastly, we employ graphics processing unit (GPU) accelerated electronic structure calculations to identify and characterize CIs in a nm Cs12Pb4Br20 nanoparticle (NP) model. It is 11 1.2 1.2 0.6 our expectation that the insights yielded by this study will guide the development of future materials with low rates of NRR. 2.3 Computational Methods Our computational study will follow the following outline. First we identify potential nonradiative recombination pathways by applying ab initio NAMD simulations to Cs4PbBr6, a cluster model of CsPbBr3, shown in Figure 2.1 (left panel). Second, we apply highly accurate static quantum chemical methods to characterize these recombination pathways and to quantify errors in the low-level PESs used in the NAMD study. Finally, we identify minimal energy conical intersections (MECIs) in a larger NP model, Cs12Pb4Br20, also pictured in Figure 2.1 (right panel). Figure 2.1: Cs4PbBr6 cluster model at the Franck-Condon (FC) point optimized at the B3LYP/SBKJC Polarized(p, 2d) level (left) and Cs12Pb4Br20 NP FC point optimized at the PBE0/SBKJC VDZ level (right). Cs, Br, and Pb atoms are represented in green, purple, and black, respectively. All molecular images in this work were created with MacMolPlt.1 2.3.1 Cluster Model Ab Initio Multiple Spawning Simulations NAMD simulations were performed using the ab initio multiple spawning (AIMS) method.85,86 Described in detail in references 85-87, the approximations in AIMS enable the solution of the full time-dependent Schrӧdinger equation including all nuclear and electronic degrees of 12 freedom explicitly. In short, the molecular wave function is expanded in a basis of frozen Gaussian trajectory basis functions (TBFs). The average positions and momenta of these TBFs evolve according to classical equations of motion and, thus, are ripe for interpretation in classical terms. The expansion coefficients are solved by integration of the time-dependent Schrӧdinger equation. The basis of TBFs is dynamically expanded via a “spawning” procedure to ensure that the basis is sufficient to describe nonadiabatic transitions between electronic states at conical intersections and avoided crossings. In the present work, the electronic wave functions and PESs were computed on the fly at the state-averaged complete active space self-consistent field (SA-CASSCF) level of theory88-90 using the SBKJC VDZ basis set and corresponding scalar relativistic effective core potentials (ECPs).91 An active space of two electrons in five orbitals (2/5) was employed with a state average over the lowest five electronic states (SA5). This active space was selected because, as described in detail in Section 3, the lowest four vertical excitations at the optimized ground-state minimum (Franck-Condon; FC) geometry were in good agreement with those computed at the equation-of- motion coupled cluster with single and double excitations (EOM-CCSD)92-94 level of theory. The FC geometry was optimized at the B3LYP95-97/SBKJC Polarized(p, 2d)98, 99 level, and the EOM- CC calculations used the SBKJC Polarized(p, 2d) basis set. The geometry optimization and vertical excitation energy calculations were all performed using the MOLPRO electronic structure software package.100-102 Twenty AIMS simulations were initialized on the first excited singlet (S1) state with initial nuclear positions and momenta sampled randomly from the ground-state vibrational Wigner distribution calculated in the harmonic approximation at the B3LYP/SBKJC Polarized(p, 2d) level of theory. The maximum time step of the AIMS adaptive multiple time step integrator was set to 13 100 atomic units (~2.42 fs), and the nonadiabatic coupling matrix elements were calculated using the norm-preserving interpolation method.103 The TBF Gaussian widths associated with the Br, Cs, and Pb atoms were determined using an optimization procedure similar to that developed by Thompson and Martínez.104 The scheme is described in detail in the Appendix. In short, the overlap of the initial nuclear wave function with the ground-state vibrational wave function was maximized, subject to the constraint that identical Gaussian widths are associated with the degrees of freedom associated with identical atoms. (That is, identical widths are assigned to all three Cartesian directions for all Br atoms, etc.) The widths of the TBFs were calculated to be 36.69, 20.96, and 24.63 bohr-2 for Br, Cs, and Pb, respectively. All AIMS simulations were performed using the FMSMolPro package.86 2.3.2 Quantum Chemical Study of Cluster Model CIs Geometries with small S1–S0 energy gaps were drawn from the AIMS simulations to serve as initial guesses for MECI geometry optimizations. These optimizations were performed using the CIOpt105 software package. Initial optimizations were carried out at the SA5- CASSCF(2/5)/SBKJC VDZ level of theory. Subsequently, more accurate optimizations including dynamic electron correlation were carried out at the multistate complete active space second-order perturbation theory (MS-CASPT2) level.106-109 The MS-CASPT2 optimizations were performed with five mixed states (MS5), an active space of two electrons in five orbitals, a level shift of 0.2 Hartree,110 the SBKJC Polarized(p, 2d) basis set, and 19 of the molecular orbitals (MOs) left uncorrelated (“cored”; abbreviated c19). Ten orbitals remain correlated (including all active orbitals) in these calculations. Local S1 energy minima were also optimized at these levels of theory. We conclude that a particular geometry is a conical intersection based on the small energy 14 gap alone. Truhlar and Mead have made a compelling argument that points of small energy gap in polyatomic systems are more likely to be conical intersections than avoided crossings.111 Two additional sets of single point calculations were performed at the MS-CASPT2-c19 optimized geometries to provide a more accurate description of dynamic electron correlation. The first set was similar MS-CASPT2 calculations without any uncorrelated (“cored”) MOs (c0). The second set were performed at the multireference configuration interaction with single and double excitations (MRCI) level with the Davidson correction112-114 using the same state average and active space as above. All electrons are correlated in the MRCI calculations. State averaged natural orbitals (SANOs) were computed at all optimized geometries. These orbitals were computed by diagonalization of the S1–S0 state-averaged first-order reduced density matrix. To investigate the possibility that spin–orbit coupling (SOC) might influence the results, calculations of the singlet–triplet SOC were performed at the cluster model FC geometry and optimized MECI geometries. Wave functions for the SOC calculations were similarly generated at the SA-CASSCF level of theory. However the basis set employed was changed to cc-PVDZ- PP115, 116 with the core electrons being replaced by fully relativistic ECPs containing SO potentials.115, 117, 118 The SOCs between the first five singlet states and first four triplet states were calculated using the MOLPRO electronic structure software package.119 Because of the large number of geometry optimizations, we only present geometric parameters, SANOs, and other computed quantities in the main manuscript, as they are required to support the primary conclusions. However, the Appendix includes the geometries and SANO figures for all optimized structures. 15 2.3.3 Optimization of Nanoparticle MECIs Minimal energy conical intersection optimizations of the Cs12Pb4Br20 NP were carried out using the GPU-accelerated implementation of the floating occupation molecular orbital complete active space configuration interaction (FOMO-CASCI)120, 121 method in the TeraChem software package.28, 29, 31, 33, 122, 123 An active space of 6 electrons in 11 orbitals was chosen. The active space selection process is described in subsection 2.4.3. The SBKJC VDZ basis and ECPs were used. The FOMO thermal smearing temperature was set to 0.25 Hartree. The FOMO- CASCI(6/11)/SBKJC VDZ MECI optimizations and S1 energy minimizations of the NP model were performed using the CIOpt software package. For a consistent comparison between the cluster and NP models, cluster model MECIs were also reoptimized using FOMO-CASCI with an active space of 2 electrons in 5 orbitals and the SBKJC VDZ basis set. This active space was chosen based on the agreement with the previously performed EOM-CCSD calculations at the FC geometry. The FOMO thermal smearing temperature was set to 0.25 Hartree, as above. 2.4 Results and Discussion 2.4.1 Photodynamics of the Cluster Model As described above, we ran AIMS simulations of the electronically excited Cs4PbBr6 cluster model to identify possible nonradiative recombination pathways. The S1–S0 energy gap of each initial AIMS TBF as a function of time after initial excitation is presented in Figure 2.2. These energy gaps are computed at the centroids of the TBFs, effectively tracking the PESs of the classical trajectories underlying the TBFs. During the first 1000 fs, all 20 initial TBFs exhibit 16 similar behavior; the energy gaps are oscillating and decreasing towards zero. After 1000 fs the dynamics bifurcate, with one group of TBFs having S1–S0 energy gaps near zero, while the second set has S1–S0 energy gaps oscillating around ~2 eV. The set with near-zero gap suggests the existence of S1–S0 CIs. Figure 2.2: The S1–S0 energy gap of the twenty initial TBFs as a function of time drawn from AIMS simulations of the Cs4PbBr6 cluster model. The populations of S0 and S1 as a function of time (Figure 2.3) further suggest the existence of CIs. An increase in the S0 population begins at 1400 fs after excitation, when it can be seen that many TBFs are exploring regions of near-zero gap in Figure 2.2. By the end of the 5 ps simulation window, the population of S0 is 0.61, indicating fast nonradiative decay. It is important to note that this is the decay time of our model cluster and does not correspond to the expected decay time of an exciton in a bulk perovskite. 17 Figure 2.3: Thick lines show the populations of S0 (red) and S1 (dark blue) as well as the aggregate S2–S4 populations (orange) averaged over all AIMS calculations. The populations of individual simulations are represented by thin lines, with S0, S1 and aggregate S2–S4 populations in pink, light blue, and light orange, respectively. We now turn our attention to the nuclear dynamics that bring about degeneracy in the S0 and S1 electronic state energies. In Figure 2.4 we present projections of the positions of all TBFs (initial and spawned) onto two coordinates of interest: the largest Cs–Pb distance and the largest Br–Pb distance. Geometries where new TBFs were spawned to S0 are marked with pink circles. Specifically, these are geometries, where the nonadiabatic coupling from S1 to S0 reaches a maximum, triggering the creation of a new trajectory on S0. During the course of the AIMS simulations it is very common for a Br and/or Cs atom to dissociate from the cluster. As indicated by the locations of the pink circles, these dissociated geometries are also where spawning to the S0 electronic state most often takes place; no spawning is observed when the longest Br–Pb and Cs–Pb distances are both less than 6 Å (compared to equilibrium values of 3.086 and 4.048 Å, respectively). Clearly, a large displacement from the FC region is required for efficient nonradiative decay to the ground state, hence the relatively long (>1 ps) delay before the onset of population transfer to S0. 18 Figure 2.4: The positions of all TBFs (initial and spawned) projected onto two coordinates: the largest Br–Pb distance and the largest Cs–Pb distance. Geometries at which these TBFs spawned new basis functions to the S0 electronic state are marked by pink circles. 2.4.2. Analysis of Cluster Model PESs To characterize the decay pathways observed in the AIMS simulations described above, as well as to quantify any errors in the shapes of the SA-CASSCF PESs used in these simulations, we have performed a static quantum chemical study of our cluster model. Because of the large variation in spawning geometries observed in our simulations, it is useful to subdivide the PES into several regions. To this end, we define six regions that were explored during the dynamics simulations: the FC region, the Cs dissociation region (CsD), the Br dissociation region (BrD), the double Br stretch region (DBrS), the CsBr dissociation region (CsBrD), and the ballooned region (Bal). The BrD region is where the most distant Br atom has dissociated from the cluster, but the most distant Cs atom has not. Similarly, the CsD region is where the most distant Cs atom has dissociated from the cluster, but the most distant Br atom has not. The DBrS region is where two Br atoms trans to one another on the cluster model have Pb–Br distances that are significantly stretched compared to the ground state minimum distance (3.08 Å), but not so stretched as to consider fully dissociated (<7 Å). The CsBrD region is where both the most distant Br and Cs 19 atoms have dissociated from the cluster. The Bal region is where several Br–Pb and Cs–Pb distances are significantly stretched but not yet fully dissociated (<7 Å). Lastly, the FC region surrounds the ground state minimum—where all Cs–Pb and Br–Pb bonds remain near their equilibrium lengths. To qualify the accuracy of SA-CASSCF in the FC region, we will consider vertical excitation energies at the FC point as computed at the SA-CASSCF level of theory compared to those computed with EOM-CCSD. We choose EOM-CCSD as a reference because it is known to accurately predict the vertical excitation energies of singly excited electronic states and does not depend on the choice of an orbital active space. The EOM-CCSD S1–S0 vertical excitation energy was calculated to be 4.57 eV with a set of triply degenerate S2–S4 electronic states just above S1 at 4.66, 4.67, and 4.67 eV relative to the S0 energy. (These states are not precisely degenerate, because the optimized cluster geometry deviates very slightly from Td symmetry.) SA-CASSCF predicts the S1–S0 vertical excitation energy to be 5.20 eV and the nearly degenerate S2–S4 states at 5.31, 5.32 and 5.32 eV. Thus, the gap between S1 and the triply degenerate S2–S4 set is 0.1 eV in both cases. The overestimation of the vertical excitation energies by SA-CASSCF is likely due to the lack of dynamic correlation, a common error in SA-CASSCF calculations. When dynamic correlation is added via MS-CASPT2-c0, the energies of states S1–S4 are in better agreement with the EOM-CCSD results—4.99, 5.08, 5.08, and 5.09 eV, respectively, while the S1–S2 gap remains 0.1 eV. Additionally, to gauge whether vertical excitation character of the cluster model is similar to the excitation character of bulk lead halide perovskite models, we compare the orbitals involved in the SA-CASSCF excitation to the orbitals corresponding to the valence band maximum and conduction band minimum known in bulk lead halide perovskite models. In bulk CsPbBr3, the 20 valence band maximum is primarily comprised of Br p orbitals and the conduction band minimum is comprised of Pb p orbitals.124 Similarly it is known that, in similar methylammonium (MA) lead triodide and CsPbI3 perovskites, the valence band maximum is comprised of the s orbital of the Pb atoms in an antibonding superposition with the p orbitals of the I atoms, and the conduction band minimum is comprised of Pb p orbitals and I p orbitals, again in an antibonding superposition. Comparing these to the fractionally occupied SANOs at the cluster model FC point displayed in Figure 2.5, we see that the n = 1.50 orbital, corresponding to the valence band maximum, has Pb s orbital and Br p orbital antibonding character, consistent with expectations. However, the n = 0.50 orbital, which is a diffuse Rydberg-like molecular orbital surrounding the Cs atoms, is inconsistent with the expected conduction band character. This is not unexpected in our model though, because the model calculations were performed in the gas phase. With the n = 0.50 orbital being so diffuse, the density of the excited electron minimally overlaps the density on the cluster, reducing repulsion and making this the preferred orbital for electronic excitation. Recognizing the source of this discrepancy, these results do not suggest that our model is incapable of accurately describing the energies of CIs present in the model NRR pathways, because the electronic structure at the CIs (analyzed in detail below) much more closely resembles that of the bulk material. That is, the diffuse Rydberg-like excitation is observed only in the FC region. 21 Figure 2.5: The fractionally occupied SANOs at the cluster model FC point, computed at the SA- CASSCF level of theory, and generated with an isovalue of 0.01. All subsequent SANOs displayed have also been generated with an isovalue of 0.01. To characterize the various decay paths observed in the dynamics, MECI and S1 minimization calculations were conducted in each PES region described above, excluding the FC region. Snapshots from the AIMS simulations were used as initial guesses. As described above, these optimizations were carried out at both the SA-CASSCF and MS-CASPT2-c19 levels of theory. The optimized MECI geometries are presented in Figure 2.6. Energies of the optimized MECIs and corresponding local S1 minima can be viewed graphically in Figure 2.7 and are listed in Table 2.1. Throughout this work, all energies are referenced to the ground-state energy at the ground-state minimum structure. Single-point MS-CASPT2-c0 and MRCI energy calculations were performed at the MS-CASPT2-c19-optimized MECI and S1 minimum geometries to provide more accurate energetics, again reported in Figure 2.7 and Table 2.1. Before discussing the optimized MECI and S1 minimum energies in detail, we will introduce the primary metric of comparison to our results, the optical band gap of bulk CsPbBr3, which is reported to be 2.36 eV.125 Direct comparison of the cluster conical intersection energies to the optical band gap of the material enables us to assess whether the intersection is energetically accessible. Intersections that fall below the optical gap can, in principle, be accessed after vertical excitation, possibly resulting in efficient NRR. With the assumption of fast cooling of charge 22 carriers to the band gap, MECIs whose energies are significantly above the gap would not be accessible and, thus, are irrelevant to the photophysics of the bulk material. The application of this argument to the cluster model MECIs is justified by past experience that CIs in larger materials typically arise from local electronic and geometric distortions, the energies of which can be accurately represented with small models.19, 21 The validity of this approach will be discussed further in subsection 2.4.3 with the results of the Cs12Pb4Br20 NP MECI optimizations. 23 Figure 2.6: Optimized Cs4PbBr6 cluster model S1–S0 MECI geometries. MECI geometries are organized by PES region, with each geometry labeled with the level of theory at which the optimization was carried out. 24 Figure 2.7: Optimized cluster model MECI and local S1 minimum energies. Energy calculations were performed at the SA-CASSCF (solid black lines), MS-CASPT2-c19 (solid red lines), MS- CASPT2-c0 (dashed orange lines) MRCI, (dashed purple lines), and the FOMO-CASCI (dotted blue lines) levels of theory, as described above. The energy calculations performed at the MS- CASPT2-c0 and MRCI levels of theory were performed at geometries optimized at the MS- CASPT2-c19 level of theory. All other calculations were performed at geometries that were optimized at their respective levels of theory. 25 Table 2.1: Optimized cluster model MECI and local S1 minimum energies. Optimization Method B3LYP SA-CASSCF MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 FOMO-CASCI MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 FOMO-CASCI SA-CASSCF MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 FOMO-CASCI SA-CASSCF MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 SA-CASSCF MS-CASPT2-c19 Geometry FC S0 Min FC S1 Min BrD CI BrD S1 Min CsD CI CsD S1 Min DBrS CI CsBrD CI CsBrD S1 Min Bal CI Bal S1 Min Level of Theory S0 Energy (eV) S1 Energy (eV) 0.00 0.00 0.00 0.00 0.00 0.77 0.44 0.96 1.03 3.31 3.02 4.16 4.02 3.16 2.85 3.90 3.83 4.20 3.80 4.45 4.61 3.95 1.02 2.70 3.41 3.46 2.54 2.24 3.49 3.23 4.18 3.92 5.29 5.10 3.67 3.81 3.18 4.49 4.23 4.09 3.68 5.26 5.05 3.59 2.49 4.00 3.60 5.20 4.99 4.50 4.80 4.87 3.86 3.50 3.89 4.21 3.34 3.04 4.27 4.42 3.17 2.94 4.18 4.27 4.21 3.82 4.52 4.73 3.98 3.81 3.75 4.42 4.70 2.54 2.26 3.79 4.00 4.20 3.94 5.79 5.80 3.67 4.11 3.80 5.59 5.66 4.11 3.70 5.78 5.75 3.68 3.24 4.83 4.94 SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI FOMO-CASCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI FOMO-CASCI MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI FOMO-CASCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI FOMO-CASCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI 26 MECI structures were found in all five regions of the PES (BrD, CsD, DbrS, CsBrD, and Bal). S1 minima distinct from the MECIs were optimized in each of these regions except for the DBrS region, where no distinct minimum was found at either the SA-CASSCF or MS-CASPT2- c19 level of theory. An additional S1 minimum is found in the FC region. At both the SA-CASSCF and MS-CASPT2-c19 levels of theory, the local MECI observed to have the lowest energy is the DBrS CI with a predicted energy of 2.54 eV at the SA-CASSCF level of theory and 2.25 eV at the MS-CASPT2-c19 level of theory. Comparison of the best estimates of the intersections energies (3.49 eV and 3.23 eV at the MS-CASPT2-c0 and MRCI levels, respectively) to the experimental band gap of bulk CsPbBr3 (2.36 eV) suggests that none of these intersections are accessible upon excitation of the bulk material. In addition, these intersections fall in the violet-to-ultraviolet energy range, suggesting that they would not introduce recombination pathways in CsPbBr3 NPs for the majority of visible-light applications. Looking to the CI geometries in Figure 2.6, significant increases in Br–Pb and Cs–Pb bond distances are observed in all MECI structures, in agreement with the motions observed in our AIMS simulations (Figure 2.4). To better understand the nature of the electronic structure surrounding these intersections, scans of the PESs surrounding the BrD and CsD intersections were performed. A one-dimensional slice of the PES surrounding the SA-CASSCF-optimized BrD CI can be seen in Figure 2.8a. Additionally, the charges of the most distant (dissociated) Br atom in both the S0 and S1 electronic states are reported in Figure 2.8b at every step of the PES scan. The point of intersection in Figure 2.8a is the BrD MECI optimized at the SA-CASSCF level. The remaining points in the scan are constructed by displacing the most distant Br atom along the vector connecting it to the Pb center. All other degrees of freedom remain frozen. 27 Figure 2.8: a) Scan of the S0 (black) and S1 (red) PESs of the cluster model as the dissociated (most distant) Br atom at the BrD CI is displaced along the vector connecting it with the Pb center. All other degrees of freedom remain frozen. The PES is computed at the SA-CASSCF level. b) The S0 (black) and S1 (red) Mulliken charges of the dissociated Br atom at the same geometries as in the PES in a). a) b) Clearly changing the Br–Pb distance lifts the S0 and S1 degeneracy. The shapes of the S0 and S1 PESs strongly suggest the crossing of an ionic bond dissociation curve and a covalent bond dissociation curve, and the charges of the dissociated Br atom are consistent with this characterization. At shorter Br–Pb distances the charge of the Br atom is zero on S1, and the S1 PES is flat, consistent with a covalent state. The S0 state similarly resembles an ionic state, with 28 a strongly attractive potential (consistent with Coulomb’s law) and a charge of ca. –1 on the most distant Br atom. Once the BrD CI geometry is passed, the dissociated Br S0 and S1 charges swap, which is consistent with the well-known picture of an ionic bond dissociation curve crossing a covalent bond dissociation curve. The CsD CI was also found to be crossing between a covalent and an ionic state, and the PES scan around this intersection can be seen in Figure A.56. These intersections are classic biradicaloid intersections,80 in which breaking of a chemical bond leads to two nearly degenerate orbitals. Taking the BrD CI as an example, the two orbitals are a p orbital on the dissociated Br atom and a p orbital on the central Pb atom, displayed in Figure 2.9. The biradicaloid nature of these intersections is similar to that observed in various silicon defects.78 Figure 2.9: Fractionally occupied SANOs at the SA-CASSCF-optimized BrD MECI geometry. Because the Br atom is very electronegative, the energy of this p orbital is quite low, favoring the ionic state—in which it is doubly occupied—when this atom remains near the remainder of the cluster. Only when it reaches a long distance, where the Coulombic energy between the Br– and the remaining cationic fragment becomes small, does the covalent state, in which this Br p orbital is only singly occupied, become energetically favorable. The large degree 29 of dissociation required to achieve degeneracy results in the relatively high energy of the intersection. The calculated properties of lead halide perovskites, such as band gaps, are known to depend strongly on the inclusion of scalar relativistic and spin–orbit effects.70, 126, 127 Above we have included scalar relativistic effects via ECPs, but neglected spin–orbit effects. It is absolutely vital to gauge their importance. To estimate how our calculated MECI energies would change if spin–orbit effects were included in our computations, we calculated the singlet–triplet SOC at the cluster model FC point and the five SA-CASSCF-optimized MECIs. Table 2.2 reports three SOCs for each geometry: the largest singlet–triplet SOCs including each of the intersecting electronic states (S0 and S1), as well as the largest calculated singlet–triplet SOC including any states. (The large number of spin–orbit couplings prevents us from reporting all of them here.) These largest values represent rough upper bounds to the energetic error of our scalar-relativistic calculations of each MECI attributable to neglect of SOC. Note that the largest SOCs vary significantly between the different MECI geometries. This is due to the nature of the electronic structure. Those states in which the Pb p orbitals are occupied have the largest SOC. The largest SOC involving the intersecting states is 0.28 eV between S1 and T3 of the DBrS MECI, the lowest-energy MECI predicted by our calculations. We roughly estimate the lower bound on the energy of this MECI after SOC correction according to where is the largest SOC involving the intersecting states. Taking our best estimates of the scalar relativistic energy of the lowest-energy MECI (DBrS; 3.49 eV and 3.23 eV at the MS- CASPT2-c0 and MRCI levels, respectively), we estimate the spin–orbit-corrected energy of the lowest (DBrS) MECI to be between 2.95 and 3.21 eV, well above the bulk band gap, 2.36 eV. 30 MECI,full-relativisticMECI,scalar-relativisticSOEEV−SOV Predicted bounds on the MECIs in the other four regions of the PES are even higher in energy. Thus, we judge that our conclusion that all identified MECIs are energetically inaccessible upon excitation across the band gap of the bulk material is insensitive to inclusion of the SOC. Table 2.2: Largest singlet–triplet SOCs calculated at the FC point and optimized MECI geometries, and the largest singlet–triplet SOCs including the intersecting electronic states (S0 and S1). Geometry SOC (eV) Singlet State Triplet State S0 S1 S4 S0 S1 S4 S0 S1 S3 S0 S1 S2 S0 S1 S2 S0 S1 S2 T1 T1 T2 T4 T4 T1 T3 T3 T4 T1 T3 T3 T3 T4 T4 T1 T3 T3 0.0004 0.0002 0.0079 0.0039 0.0680 0.0681 0.0010 0.0002 0.0107 0.0209 0.2815 0.3852 0.0012 0.0039 0.0124 0.0009 0.0046 0.0080 FC S0 Min BrD CI CsD CI DBrS CI CsBrD CI Bal CI 2.4.3 FOMO-CASCI PES Validation and Active Space Selection for NP Model Having identified and characterized the recombination pathways of the smaller cluster model, we will turn our attention to the larger NP model. To this end, we will use the FOMO- CASCI method. FOMO-CASCI is preferable to SA-CASSCF for treatment of larger systems for two reasons. First, the large number of nearly degenerate orbitals can lead to instabilities in the SA-CASSCF equations, even for relatively small cluster models of semiconductors.128 Second, SA-CASSCF does not yield a size-intensive description of electronic excitations.129 FOMO- 31 CASCI provides solutions to both of these problems. FOMO-CASCI is also less computationally expensive than SA-CASSCF in most cases. To validate the FOMO-CASCI method for this system, we will first compare the MECIs of the smaller cluster model optimized at the FOMO-CASCI level to the SA-CASSCF-optimized MECIs presented in Subsection 2.4.2. For the reoptimization of cluster model MECIs using FOMO-CASCI, an active space of two electrons in five orbitals was chosen because the predicted vertical excitations energies (4.87 eV for the S1–S0 transition) and the excitation character (see SANOs in Figure 2.10) are in excellent agreement with those predicted by EOM-CCSD and MS- CASPT2 (4.57 eV and 4.99 eV, respectively, Figure 2.5). Figure 2.10: Fractionally occupied SANOs of the cluster model computed at the FOMO-CASCI level of theory at the FC geometry. The optimized FOMO-CASCI cluster model MECI geometries are displayed in Figure 2.6, and energies of the optimized MECIs are listed in Table 2.1 and can be viewed graphically in Figure 2.7. As seen in Figure 2.6, the FOMO-CASCI optimized MECI geometries in the BrD, CsD, and CsBrD PES regions are similar in energy and structure to those optimized at the SA- CASSCF and MS-CASPT2-c19 levels of theory. The similarity suggests that the FOMO-CASCI method is a good choice for the NP model MECI optimizations. Note that, like SA-CASSCF, 32 FOMO-CASCI significantly underestimates the MECI energies relative to the dynamically correlated MS-CASPT2-c0 and MRCI results (Figure 2.7 and Table 2.1). In subsection 2.4.4, we will discuss how we correct for this error in the NP calculations. No DBrS intersection was identified at the FOMO-CASCI level of theory. This may raise concern, as the DBrS intersection was the lowest-energy MECI at the SA-CASSCF and MS- CASPT2 levels of theory. However, as will be described in subsection 2.4.4, the DBrS intersection is predicted in the large NP cluster and remains the lowest-energy intersection. We hypothesize that the DBrS region of the seam connects to the BrD region, making its optimization difficult in the cluster model case. Having validated the FOMO-CASCI method for the small cluster model, we turn our attention to the selection of the active space for the Cs12Pb4Br20 NP model. Typically we choose active spaces by comparison to a more trustworthy dynamically correlated wave function-based level of theory (e.g. CASPT2 and/or EOM-CCSD). This is unfeasible for the larger NP model, so instead we choose to compare vertical excitation energies to the linear response time-dependent density functional level of theory with the PBE0 functional (TD-PBE0) with the SBKJC VDZ basis and ECPs. TD-PBE0 was chosen because it predicts excitation energies of the smaller cluster model that are in excellent agreement with EOM-CCSD; the TD-PBE0 vertical excitation energies of S1–S4, computed at the PBE0 S0 minimum energy geometry, were found to be 4.48, 4.58, 4.58 and 4.58 eV at the TD-PBE0 level. These compare well with the EOM-CCSD values (4.57, 4.66, 4.67, and 4.67 eV). For the NP model, a 6-electron/11-orbital active space was chosen, because, among the various active spaces we considered, the vertical excitation energies and excited electronic state characters of the NP model were found to be in the best agreement with the TD-PBE0 level. 33 Specifically, there is good agreement between the 4.72 eV FOMO-CASCI S1–S0 excitation energy and the 4.00 eV TD-PBE0 S1–S0 excitation energy, considering the lack of dynamic electron correlation at the FOMO-CASCI level. All of the other active spaces considered overestimate the TD-PBE0 result more dramatically. The SANOs for the FOMO-CASCI calculations are shown in Figure 2.11. There is excellent agreement between the n = 1.57 SANO and the expected valence band maximum character (Br p orbitals and Pb s orbitals). Similar to the smaller cluster model, the n = 0.43 SANO is Rydberg-like and not in good agreement with the expected conduction band character. As argued above, this result is not unexpected, given that these calculations are performed in the gas phase. This discrepancy does not by itself raise doubt regarding the accuracy of computed NP model MECI energies, as the Rydberg-like excitation is only observed in the FC region of the PES. Figure 2.11: Fractionally occupied SANOs computed at the FOMO-CASCI(6/11)/SBKJC VDZ ECP level of theory at the NP FC geometry. 2.4.4 Analysis of NP Model PESs Using initial guesses generated using knowledge of the cluster model MECIs, we have performed both MECI optimizations and excited-state optimizations in the BrD, CsD, and DBrS regions of the PES of the NP model. We focus on these three regions because the other two regions of the PES yielded very high-energy MECIs in our study of the cluster model. In the NP model, 34 an MECI was found in each of the three regions studied. The optimized MECI geometries are shown in Figure 2.12. The energies of the MECIs and excited-state minima are listed in Table 2.3 and can be viewed graphically in Figure 2.13. The FOMO-CASCI energies are approximately corrected for dynamic electron correlation using the computed energies of the cluster model at the FOMO-CASCI and higher (MS-CASPT2-c0, MRCI) levels of theory according to Here “higher-level” indicates either MS-CASPT2-c0 or MRCI. This correction is justified . by the local nature of the electronic and nuclear distortions predicted at the conical intersections. Because no DBrS could be found in the cluster model at the FOMO level, is replaced by in this case. Figure 2.12: Cs12Pb4Br20 nanoparticle optimized geometries. MECI optimizations were performed at the FOMO-CASCI level of theory, while the FC geometry was optimized at the PBE0 level. 35 correctedFOMOFOMOhigher-level(NP)(NP)(cluster)(cluster)EEEE=−+FOMO(cluster)ESA-CASSCF(cluster)E Figure 2.13: Energies of the NP model at optimized geometries. Corrections are computed as described in the text. The experimental optical band gap of CsPbBr3 (purple dashed line) was taken from reference 125. Table 2.3: Energies of the NP model at optimized geometries. Corrections are computed as described in the text. Geometry FC S0 Min FC S1 Min BrD CI Optimization Method PBE0 FOMO-CASCI FOMO-CASCI FOMO-CASCI FOMO-CASCI BrD S1 Min DBrS CI FOMO-CASCI FOMO-CASCI CsD CI FOMO-CASCI Level of Theory S0 Energy (eV) S1 Energy (eV) FOMO-CASCI FOMO-CASCI FOMO-CASCI FOMO-CASCI with MS-CASPT2-c0 correction FOMO-CASCI with MRCI correction FOMO-CASCI FOMO-CASCI FOMO-CASCI with MS-CASPT2-c0 correction FOMO-CASCI with MRCI correction FOMO-CASCI FOMO-CASCI with MS-CASPT2-c0 correction FOMO-CASCI with MRCI correction 0.00 0.92 2.78 3.78 3.64 2.27 2.50 3.45 3.19 4.59 5.09 5.25 4.72 3.36 2.78 3.88 4.03 2.48 2.51 3.76 3.97 4.60 5.14 5.35 36 The geometries of the BrD, CsD, and DBrS conical intersections are quite similar to those reported in the smaller cluster model. After approximate correction for dynamic correlation, the energetic ordering is similar. Corrected at the MRCI level, the DBrS intersection is the lowest in energy at 3.19 eV above the ground-state minimum, followed by the BrD intersection at 3.64 eV, and finally the CsD MECI at 5.25 eV. These numbers compare respectively to 3.23, 3.83, and 3.46 eV in the smaller cluster. The DBrS intersection is the lowest in both cases. Clearly the CsD MECI is disfavored in the larger NP model, suggesting it is less likely to be easily accessed in realistic materials. Most importantly, none of these intersections are accessible at the bandgap energy, and thus efficient nonradiative recombination pathways involving conical intersections do not appear to exist on the pristine surface of CsPbBr3. As with the smaller cluster, approximate spin–orbit correction does not change the story; all three intersections remain well above the bulk band gap (Figure 2.13, blue dashed lines). Though not accessible at the bulk band gap energy, nanostructured materials with larger band gaps may still be able to access the intersection regions. As such, it is worth discussing the topography of the PES in the region surround these intersections. Though MECIs were found in all three regions of the PES in which we searched, only the BrD region has a local S1 minimum distinct from the intersection seam. The S1 minimum in the BrD region is found to be only 0.03 eV below the MECI at the FOMO-CASCI level. Thus, efficient transitions to the ground state would be expect in all three regions of the PES, as the MECI is the lowest or nearly lowest point on S1 in these regions. 37 Figure 2.14: Fractionally occupied SANOs averaged over the S0 and S1 electronic states computed at the FOMO-CASCI level of theory at the NP optimized DBrS CI. To better understand the electronic structure at the lowest-energy MECI, we consider the fractionally occupied SANOs of the NP DBrS CI, which are presented in Figure 2.14. The excitation character seen here is different from what was observed in the cluster model DBrS MECI. The electronic excitation is surprisingly nonlocal. The electron is excited from a Pb atom with a relatively undistorted coordination environment to a more distorted Pb atom (around which two Pb–Br bonds are significantly stretched, with Pb–Br distances of 4.32 and 4.38 Å compared to the equilibrium Pb–Br distance, 3.09 Å). Specifically, the electron is excited from an s orbital on the relatively undistorted Pb atom (n = 1.49) to the p orbital of the more distorted Pb atom (n = 0.50). Both orbitals have significant density on the surrounding Br p orbitals, which form an antibonding superposition with the orbital on the Pb atom. Like the majority of conical intersections identified in semiconductors,78 this is a biradicaloid intersection. However, this marks a never-before-identified class of biradicaloid intersection in which the two radical centers exist on identical atoms of neighboring unit cells of the material. The excitation that brings about the intersections thus actually requires charge transfer between unit cells. Such charge transfer has an energetic cost, which likely contributes to the fact that this intersection is not accessible upon band gap excitation of the bulk material. 38 Figure 2.15: a) Scan of the PES surrounding the NP model BrD CI. The dissociated Br atom is displaced along the vector between it and the closest Pb center while all other degrees of freedom are frozen. Energies are calculated at the FOMO-CASCI level of theory. b) The S0 and S1 Mulliken charges on the dissociated Br atom as a function of the same coordinate. a) b) In order to determine if the BrD intersection of the NP has the same biradicaloid character as observed in the cluster model, we computed PESs and atomic charges along a dissociation coordinate just as we did for the cluster model in Figure 2.8, displacing the dissociated Br atom along the vector connecting it to the nearest Pb atom. The resulting PES and atomic charges are presented in Figure 2.15a and Figure 2.15b respectively. Just as in the cluster model, this intersection appears to be between a covalent and ionic state. In the PES, the intersection occurs 39 between a flat (covalent-like) and sloped (ionic-like) surface. The charge on the dissociated Br atom changes from –1.0 (0.0) to 0.0 (–1.0) as S0 (S1) passes through the intersection, from left to right. The SANOs computed at the NP BrD MECI (Figure 2.16) also support the conclusion that this is a biradicaloid intersection. Figure 2.16: Fractionally occupied SANOs calculated at the NP BrD MECI geometry. Orbitals are computed at the FOMO-CASCI level as described in the text. 2.5 Conclusions In this work we investigate possible pathways for nonradiative recombination in an archetypal lead halide perovskite, CsPbBr3. Using both a molecular cluster model and a model nanoparticle, we identified and characterized several distinct conical intersections that could potentially introduce pathways for nonradiative recombination. Most importantly, none of these intersections were predicted to be energetically accessible upon excitation across the band gap of the bulk material. Thus, they will not provide efficient pathways for nonradiative recombination of the bulk material, consistent with the fact that these materials are known to exhibit low recombination rates. Analysis of the electronic structures of the intersecting states at the intersection points provides insights into the reason why these intersections are not energetically accessible. Such 40 information can guide the design of future materials. All three intersections identified in the nanoparticle molecules required charge transfer over a significant distance in one of the intersecting states. Such transfer is energetically costly and certainly results in intersections that are higher in energy than the more local excitations that we have previously observed at intersections in silicon. Specifically, in the lowest-energy intersections (the double Br stretch or DBrS intersection), the excitation at the intersection geometry involves the transfer of an electron between Pb atoms in neighboring unit cells. The next lowest intersection (the Br dissociation or BrD intersection) required charge transfer across a dissociated Pb–Br bond. The degree of dissociation required to reach the intersections is large due to the ionic character of the Pb–Br bond. These results suggest that semiconductors with relatively ionic bonds are advantageous due to their lack of low-energy conical intersections. As discussed in detail in a recent review by the authors,78 halides are particularly advantageous, because it is difficult to form defects with more covalent bonds in these materials. 2.6 Acknowledgments We are grateful for the financial support of the NSF through grant CHE-1565634 and Michigan State University through a Strategic Partnership Grant. We are also grateful for a generous allocation of supercomputer time from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the NSF under Grant ACI-1548562 (allocation CHE-140101). 41 CHAPTER 3: STATE-PAIRWISE DECOHERENCE TIMES FOR NONADIABATIC DYNAMICS ON MORE THAN TWO ELECTRONIC STATES Reproduced from “State-pairwise decoherence times for nonadiabatic dynamics on more than two electronic states”, Michael P. Esch and Benjamin G. Levine, J. Chem. Phys. 152, 234105 (2020), with permission of AIP Publishing. 3.1 Abstract Independent trajectory (IT) nonadiabatic molecular dynamics simulation methods are established as powerful tools for modeling processes involving transitions between electronic states. Incorporation and refinement of decoherence corrections into popular IT methods, e.g. Ehrenfest dynamics and trajectory surface hopping, is an important means of improving their accuracies. In this work, we identify a new challenge in the development of such decoherence corrections: that when a system exists in a coherent superposition of three or more electronic states, coherences may decay unphysically when the decoherence correction is based on decoherence times assigned on a state-wise basis. As a solution, we introduce decoherence corrected Ehrenfest schemes based on decoherence times assigned on a state-pairwise basis. By application of these methods to a set of very simple one-dimensional model problems, we show that one of these state- pairwise methods (“collapse to a block”; TAB) correctly describes the loss of coherence between all pairs of states in our multistate model problems, whereas a method based on a state-wise description of coherence loss does not. The new one-dimensional models introduced here can serve as useful tests for other decoherence correction schemes. 42 3.2 Introduction Accurate and efficient simulation methods for modeling nonadiabatic processes are essential to progress in a range of fields involving electronic excited states, e.g. molecular photochemistry,16, 130 spectroscopy,131, 132 and materials science.133-136 The development of methods for simulating nonadiabatic molecular dynamics is not only important but also challenging. Over decades, successful methods have ranged from simple and elegant to thoughtfully complex. Out of early work in this area arose two dominant independent trajectory (IT) approaches: Ehrenfest (mean-field) molecular dynamics36, 38, 137, 138 and trajectory surface hopping.40, 139 IT methods replace the nuclear wave function with a swarm of classical trajectories that need not communicate with one another. The simplicity of this approach has led to the broad adoption of IT methods and provides an elegant starting point for further refinements. However, both Ehrenfest and trajectory surface hopping suffer from well known problems describing electronic coherence. Individual Ehrenfest trajectories propagate on mean-field potential energy surfaces (PESs) associated with coherent superpositions of electronic states, for example.40 This mean-field nature results in unphysical reaction outcomes in some cases. Elimination of unphysical trajectories by windowing has been proposed as a solution to this problem.140 Perhaps the most broadly adopted solution, however, is trajectory surface hopping.40 In its most widely used forms, trajectories propagate on only one adiabatic PES at a time, eliminating the unphysical mean-field behavior associated with Ehrenfest dynamics. However, the popular fewest switches algorithm also suffers from well known issues treating electronic coherences, arising because the propagation of the electronic wave function does not take into account the loss of coherence between populations propagating on nonparallel PESs.44, 141-143 Solving this problem has been the focus of much attention. One alternative 43 approach is to replace IT schemes by more sophisticated Gaussian wave packet methods, in which the nuclear time-dependent Schrodinger equation is solved in a small, classically-inspired, time- dependent basis.85, 87, 144-147 Such methods naturally incorporate decoherence effects in a systematically improvable first-principles fashion, but lose some of the appealing simplicity of IT methods. Rigorous path-integral methods have also been introduced to incorporate both electronic decoherence and decoherence associated with the environment.148-151 Many modifications of both surface hopping and Ehrenfest dynamics have been designed to correct for decoherence effects, as well,45, 47, 50, 52, 53, 55-57, 59, 61, 152-160 and in many cases, such methods can provide a more accurate physical description of molecular dynamics on multiple electronic states, while maintaining much of the desirable simplicity of their parent methods. In this work, we will focus on such decoherence corrections to IT methods. It is known that the loss of coherence between populations on a pair of electronic states is driven by the differences in the shapes their respective PESs. With this in mind, in the context of an IT method, the rate of decoherence can be approximated in terms of the instantaneous difference in the forces on those two states.44, 143 Though this state-pairwise nature of the decoherence time is known, most decoherence corrections to IT methods define the decoherence time in a state-wise fashion. That is, for a given trajectory in the ensemble, and at a particular time, a single decoherence time is assigned to each state rather than to each pair of states. This state-wise decoherence time may be defined either relative to the mean-field PES or relative to a particular pointer state. In the latter case, which is common in surface hopping schemes, the definition of the decoherence time may be identical to the state pair-wise definition, but decoherence between a pair of states, ij, is not explicitly treated when neither i nor j is the current pointer state for a given trajectory. 44 In systems where the dynamics of interest are confined to two electronic states (including many interesting photochemical and photophysical systems as well as the most common model problems used to test nonadiabatic molecular dynamics methods), a state-wise treatment is likely to be indistinguishable from a state-pairwise one. However, in systems with more than two electronic states, intuition suggests that a state-wise approach may fail in many cases. In a system with many electronic states, some pairs of states may have PESs that are very nearly parallel (e.g. a manifold of Rydberg states arising from a common cation core). However, there will be other states that are not parallel to this set. Imagine, for example, that states I and J are nearly parallel, but state K has a very dissimilar PES. In an exact treatment, the population on state I would maintain coherence with the population on state J for a long time, but coherence with the population on state K would be lost very quickly. In this scenario, it appears that assigning a single decoherence time to state I (describing its coherence loss with both states J and K) would be impossible. In more general terms, simply counting the number of independent coherences in the electronic density matrix ( , where N is the total number of electronic states) and comparing to the total number of decoherence times in a state-wise treatment (N), indicates that a state-wise approach lacks the flexibility to completely describe the independent dynamics of the coherences when many states are present. In such cases, it appears that defining the decoherence time as a state-pairwise property may enable a more accurate treatment. In this work, we will investigate whether this is the case by applying three IT methods to a set of simple model problems. All three methods employed here are variants of Ehrenfest dynamics, though the ideas introduced in this work would also be applicable to other IT methods (e.g. surface hopping). One of the IT methods employed here (collapse to a state, TAS) is based on a state-wise description of decoherence, while the other two (independent pairwise stochastic 45 ()12NN− decoherence, IPSD, and collapse to a block, TAB) utilize a state-pairwise description of decoherence. All three methods are applied to three simple one-dimensional models, two of which involve more than two electronic states. Rather than focusing on reproduction of transition probabilities, which are indirect measures of the coherences at best, we choose to directly analyze the off-diagonal elements of the density matrix to determine the extent to which these methods provide a physically correct description of coherence loss. In section 3.3, we describe our decoherence corrected Ehrenfest algorithm and the three decoherence corrections that are employed. We also describe the three model systems that we use for comparison. In section 3.4, we compare and analyze the behavior of our three methods in the three simple test cases. In section 3.5, we draw conclusions and discuss future prospects. 3.3 Computational Methods The details of our computational study are presented in three subsections. In subsection 3.3.1, we describe the general DCE framework employed in this work. In subsection 3.3.2, three stochastic decoherence corrections that can be plugged into this framework are introduced. These three schemes will be compared by application to three one-dimension models described in subsection 3.3.3. 3.3.1 DCE Algorithm This subsection presents the steps that comprise a single classical time step of the DCE algorithm used in this work. Each of the three decoherence corrections compared in this work is plugged into this algorithm (at step 8). The following procedure defines a single classical time step from time t to : 46 tt+ 1) The electronic Hamiltonian, is constructed, where is the vector of classical positions at time t. 2) The Ehrenfest force vector is calculated according to , (3.1) where is the time-dependent electronic wave function at time t. The classical coordinates are integrated to using the Velocity Verlet algorithm.161 3) The electronic wavefunction, , is numerically propagated to the classical half time step, , according to the electronic time-dependent Schrödinger equation, This propagation is carried out with a second-order symplectic split operator integrator.162 . (3.2) An electronic time step, Δte, shorter than the classical time step (Δt) is used. 4) The Hamiltonian at the end of the classical time step, , is constructed and diagonalized to obtain the adiabatic electronic states. 5) Employing the same algorithm as described in Step 3, the electronic wavefunction is propagated from time to according to the time-dependent Schrödinger equation, 6) The Ehrenfest force, , is calculated . (3.3) 7) The forces and are used in the velocity Verlet algorithm to integrate . (3.4) the velocities of each classical degree of freedom to t + Δt. 47 ()()tHR()tR()()()()()()MFtttdtdtt=−ΗFR()t()tt+R()t/2tt+()ditdt=H()()tt+HR/2tt+tt+()dittdt=+H()MFtt+F()()()()()()MFttttttdttdtttt++++=−++ΗFR()MFtF()MFtt+F 8) One of the three decoherence corrections described in subsection 2.2 is applied, potentially altering , , and the classical velocities at t + Δt. Steps 1 through 8 are repeated until the specified simulation length is reached. 3.3.2 Decoherence Corrections Here we describe the three decoherence corrections that are compared in this work. 3.3.2.1 Collapse To A State (TAS) The key feature of the TAS correction scheme is that decoherence times are assigned on a state-wise basis. TAS is very similar to the mean-field with stochastic decoherence method of Schwartz,50 replacing the collapse probabilities with those from Jaeger and Prezhdo’s decoherence induced surface hopping method.55 We find this replacement provides collapse probabilities that are consistent with the electronic amplitudes prior to collapse. Below we describe the algorithm in detail. In TAS, the probability that the population on an individual electronic state, k, loses coherence with all other electronic populations during a classical time step is defined to be , where τk is the mean-field decoherence time of state k. This decoherence time is defined , (3.5) (3.6) where for classical degree of freedom η, αη is a decoherence parameter (which can be thought of as the width of a Gaussian wave packet in units of inverse length squared), , (3.7) 48 ()tt+()MFtt+F1ktckPe−=−()2,,22avgavgkMFkFFћ−−=()(),,,2kkavgkFtFttF++= and . (3.8) Equations (3.7) and (3.8), respectively, define elements of the force on adiabatic electronic state k and the Ehrenfest (mean field) force, averaged over a single classical time step. To determine whether the population on state k has lost coherence with all other populations during a classical time step, a random number, γ, is drawn from the uniform distribution on interval [0,1). If γ is less than , the coherence is considered to be lost in that classical time step. This procedure is repeated for all electronic populations. When the electronic population on state k is determined to have lost coherence with the remainder of the states, the probability of collapse into state k is where is the population of adiabatic state k computed from the electronic wave , (3.9) function at the end of the classical time step. A second random number, ζ, is drawn from the uniform distribution on interval [0,1), and if ζ is less than the electronic wavefunction collapses into state k. In this event the electronic amplitude of state k, is set to , (3.10) and the electronic amplitudes of all other electronic states are set to zero. However, if ζ is greater than the amplitude of state k is set to zero and all other electronic amplitudes (l ≠ k) are scaled according to . (3.11) 49 ()(),,,2MFMFavgMFFtFttF++=ckP()kMFkkPtt=+()kktt+kMFP()()()'kkkcttcttctt++=+kMFP()()()()'1/21llkkcttctttt++=−+ This procedure preserves the relative populations and the relative phase angles of the amplitudes on the remaining electronic states. If multiple electronic populations lose coherence with all others within a single classical time step, all collapses are carried out in a randomly determined order. Each collapse changes the electronic amplitudes. Thus, as each collapse is realized, all electronic populations and amplitudes are updated, and is recalculated for the remaining states with the new electronic populations. Once all collapses have taken place for that classical time step, the velocities of classical degrees of freedom are rescaled so that the total energy of the trajectory is conserved. The Ehrenfest forces are also recalculated using the new wavefunction for subsequent propagation. Note that in the present work, only one-dimensional models are studied. In multidimensional cases rescaling the velocity in the nonadiabatic coupling direction, if available, may be preferable to simple rescaling of the velocity vector, but in the present case they are equivalent. In addition, we note that frustrated collapses are not addressed in this initial work because the models we study below never require their consideration. Similar solutions to those applied in other IT methods would likely suffice,163-165 but we leave investigation of frustrated collapses to a future work. 3.3.2.2 Independent Pairwise Stochastic Decoherence (IPSD) The key feature of the IPSD correction scheme is that decoherence times are assigned on a state-pairwise basis, rather than state-wise (as in TAS). The added flexibility of this approach enables populations on pairs of states with nearly parallel PESs to remain coherent for longer, while populations on pairs of states with very different PESs may lose coherence faster. In IPSD, the loss of coherence between pairs of states is treated by introducing stochastic collapse events 50 kMFP that move all population from one adiabatic state into another. Below we describe the algorithm in detail. In IPSD, the probability that a pair of electronic populations on states i and j lose coherence during a classical time step is defined to be (3.12) where the state-pairwise decoherence time, τij, is defined to be . (3.13) Similar to Equation (3.6), αη is the decoherence parameter for classical degree of freedom η, and the time step averaged adiabatic forces are defined in Equation (3.7). This probability is calculated for all pairs of electronic populations at the end of the classical time step. To determine if a pair of populations loses coherence during that classical time step, a random number, γ, is drawn from the uniform distribution on interval [0,1). If γ is less than , that pair of electronic populations is considered to have lost coherence. For a pair of populations that has lost coherence, a second probability, , (3.14) is calculated to determine whether population will be collapsed into state i or j. A second random number, ζ, is drawn from the uniform distribution on interval [0,1). If ζ is less than the electronic amplitude of state i is set to (3.15) and the amplitude of state j is set to zero. If ζ is greater than , the electronic amplitude of state i is set to zero and the amplitude of state j is set to 51 1ijtcijPe−=−()2,,22avgavgijijFFћ−−=cijP()()()iiijiijjttPtttt+=+++ijP()()()()()()1'2iiiijjicttcttttttctt++=++++ijP . (3.16) All amplitudes not corresponding to states i and j remain unchanged. If multiple pairs of populations are determined to lose coherence in the same classical time step, the order of collapses is randomized. Each time a collapse is performed, the electronic amplitudes and populations are adjusted before additional collapses are performed. After all collapses are performed, the velocity is scaled such that the total energy is conserved, and the Ehrenfest forces are recalculated using the new wavefunction. 3.3.2.3 Collapse To A Block (TAB) Compared to IPSD, the TAB scheme includes a more sophisticated treatment of collapse based, again, on state-pairwise decoherence times. In TAB, a set of subsets of electronic states are generated at the end of each classical time. Then, pairwise decoherence times are used to determine the probability of collapsing into a coherent superposition of each subset. Thus, the collapse procedure in TAB modifies all electronic amplitudes simultaneously, compared to IPSD which modifies only two electronic amplitudes at a time. Below we describe the algorithm in detail. First, we build an approximate electronic density matrix at the end of the time step, , assuming exponential decay of each off-diagonal element over the time step. Specifically, this density matrix has diagonal elements and off-diagonal elements 52 (3.17) (3.18) ()()()()()()1'2jjiijjjcttcttttttctt++=++++dρ()dciiiitt=+()ijtdcijijtte−=+ where is the fully-coherent electronic density matrix associated with the Ehrenfest trajectory, and τij is the state-pairwise decoherence time defined exactly as it was in IPSD (Equation (3.13)). Now we define a procedure by which the electronic wave function may stochastically collapse into a subset of the occupied adiabatic states. Specifically, steps 1 and 2 of the procedure below will generate a set of coherent density matrices, , and associated weights, , such that . (3.19) Each represents a pure electronic state that spans a subset of the adiabatic electronic states occupied in . This subset may include anywhere from all of the states occupied in (in which case will be identical to ) down to just a single state (in which case will contain only a single non-zero element). Next, steps 3 and 4 define a procedure for stochastically collapsing the electronic wavefunction into the state described by one of the generated . In this way, each trajectory will remain in a pure electronic state at all times, but the ensemble average will have coherences that decay exponentially, as supposed in the definition of . The specific procedure for generating and collapsing into is as follows: 1) The density matrix, , is calculated from according to Equations (3.17) and (3.18). A copy of , which we will call ρtmp, is created to be iteratively updated in step 2. 2) We loop over the following steps (with loop index a) to build each and compute the associated : a. A matrix b is computed from the current ρtmp, with all elements defined 53 cρblockaρblockaPdblockblockaaaP=ρρblockaρcρcρblockaρcρblockaρblockaρdρblockaρdρcρdρblockaρblockaP . (3.20) b. The smallest, non-zero element of b, bkl, is found. In this work, the absolute value of an element must be greater than a user-chosen threshold of 10-10 to be considered non-zero. c. Now we construct, β, the set of electronic states to be populated in . Set β is initiated with states k and l (identified in the previous step). For all states , if then state p is added to β. Having populated β, we consider all pairs of states, r , (3.21) and s, such that ; ; ; and . If then one of states r and s is randomly chosen and excluded from β. Once a state, , (3.22) s, is excluded in this way, additional pairs including s are no long considered, so as to avoid excluding states unnecessarily. d. An electronic density matrix that represents a coherent superposition of the states in β, , is constructed. The elements (diagonal and off-diagonal) of are defined to be , (3.23) where and . Any element where and/or is set to zero. The sum over the populations of all electronic states included in β in the denominator of Equation (3.23) is to ensure conservation of population during 54 tmpijijcijb=blockaρ,pkl0 and 0plpkbb,rkl,sklrs0rsb=blockaρblockaρ()(),cmnblockamnciiitttt+=+mn,blockamnmn collapse. This procedure ensures that the ratios of populations and phases of coherences of pairs of states in β are the same in as in . e. The weight of each block in Equation (3.19) is computed according to f. The temporary density matrix, ρtmp, is updated according to . . (3.24) (3.25) g. Repeat, incrementing loop index a, until every element of b is zero (including diagonal elements). The end result is a set of coherent density matrices that sum to the approximate density matrix as stated in Equation (3.19). 3) A random number, γ, is drawn from the uniform distribution on interval [0,1) to determine which coherent subset of electronic states the electronic wavefunction collapses into. Taking advantage of the fact that We will collapse into when 4) Now, is assigned , . (3.26) . (3.27) (3.28) and the electronic wavefunction is reconstructed from the new . The velocities of the classical degrees of freedom are scaled to conserve total energy and the Ehrenfest forces are recalculated from the new wavefunction. It is noteworthy that the decomposition in Eq. (3.19) is not unique. Our approach is a greedy algorithm. By zeroing the smallest element of b at each iteration, the algorithm identifies 55 blockaρcρ,tmpblockklablockaklP=tmptmpblockblockaaP−ρρρ1blockaaP=blockaρ and blockblockbbbabaPPcρcblockaρρcρ the largest subset of states that can be coherently coupled, and assigns it the largest possible probability. In some rare cases (e.g. in some pathological multidimensional cases or due to accumulated numerical error), we note that Equation (3.24) may yield a negative probability of collapse. We will address this in two cases: Numerically zero negative probabilities and numerically non-zero negative probabilities. We define a negative to be numerically zero if its absolute value is less than 10-9. If a numerically zero is calculated, after ρtmp is updated in Step 7, is simply set to zero. Numerically non-zero probabilities are only observed in systems with more than one nuclear degree of freedom, and therefore did not pose a problem in the current study. Nonetheless we note that a possible addition to TAB to eliminate the negative probability while preserving the current basis, is applying a linear least squares procedure to solve Equation (3.19), constraining each to be non-negative. 3.3.3 Model Problems and Initial Conditions In this work we apply the above methods to very simple models comprising one- dimensional, uncoupled, linear adiabatic PESs. The fact that there is no coupling between electronic states may seem strange, but it allows us to focus on the primary goal of this study: accurately describing the loss of coherence between pairs of electronic populations in systems that exist in a coherent superposition of three or more electronic states. With this goal in mind, we will investigate the coherences themselves rather than population transfer probabilities, which are only an indirect measure of the effects of coherence. Three models were used, with adiabatic PESs pictured in Figure 3.1 and PES surface parameters and initial electronic populations listed in Tables 3.1–3.3. First, a two-state model 56 blockaPblockaPblockaPblockaP (Figure 1a; Table 1) was used to demonstrate the behavior of all three methods in an “easy” case. A three-state model (Figure 1b; Table 2) and five-state model (Figure 1c; Table 3) are used to test the ability of the three methods to accurately describe coherence loss in systems with three or more states. For each of the nine combinations of method and model, an ensemble of 10000 trajectories was computed, differing only in their initial random seed. Classical and electronic time steps (Δt and Δte) of 2.5 and 0.00625 atomic time units (au) were used, respectively. The two- and three- state model trajectories were propagated for 1000 au and the five-state model trajectories were propagated for 1200 au. The decoherence parameter, α, was assigned to be 4.7 bohr-2, which corresponds roughly to the width of the ground state vibrational wave function of a hydrogen atom in an organic molecule.85, 104 A mass of 22500 au (comparable to a carbon atom) was chosen for numerical convenience. Lastly, the particle’s initial position was set to zero in all models, and the momentum was set to 30.000, 36.742, and 47.434 atomic units for the two-, three-, and five-state models, respectively. Initial electronic state populations are also reported in Tables 3.1–3.3. (All amplitudes are initially positive and real, though the initial complex phase of the amplitudes does not influence the results presented below.) 57 Table 3.1: The slopes, y-intercepts, and initial populations of the two-state model. State 0 1 Slope (Hartree/Bohr) -0.015 0.000 Intercept (Hartree) 0.020 0.030 Initial Population 0.25 0.75 Table 3.2: The slopes, y-intercepts, and initial populations of the three-state model. State 0 1 2 Slope (Hartree/Bohr) -0.030 0.000 0.000 Intercept (Hartree) Initial Population 0.050 0.060 0.070 0.500 0.250 0.250 Table 3.3: The slopes, y-intercepts, and initial populations of the five-state model. State 0 1 2 3 4 Slope (Hartree/Bohr) -0.040 -0.030 -0.020 -0.010 0.000 Intercept (Hartree) 0.030 0.040 0.050 0.060 0.070 Initial Population 0.20 0.20 0.20 0.20 0.20 58 Figure 3.1: The adiabatic PESs of the a) two-state, b) three-state, and c) five-state models. 59 3.4 Results and Discussion Like many decoherence corrections, all three methods applied in this work assume electronic coherences decay at an exponential rate determined by the difference in the gradients of the involved electronic states. Using this assumption, the absolute electronic coherence, , of linear, uncoupled PESs is analytically represented , (3.29) where is the initial electronic coherence and τij is the state-pairwise decoherence time defined in Equation (3.13). (The decoherence times are constant due to the linearity of the PESs.) Below, we apply TAS, IPSD, and TAB to the three models described above to test whether the coherences decay according to this expected exponential decay. Hereafter, we will call this expected function the target. Although the complex phases of the coherences are also important to electronic dynamics, here we focus only on the decay of the absolute electronic coherences. The complex phases are determined by the Ehrenfest propagation and unaffected by the collapse procedures, and thus one should not expect the choice of decoherence correction to affect them. Throughout, the reported ensemble averaged absolute coherences are computed according to , (3.30) where N is the total number of trajectories in the ensemble and I indexes the trajectories. 3.4.1 Two-State Model For the two-state model, the absolute coherences are shown with the target decay calculated using Equation (3.29) in Figure 2.2. It can be seen that TAS, IPSD, and TAB all predict coherences in excellent agreement with the target exponential decay (black line, Equation (3.29)). This 60 ()ijt()()0ijtijijte−=()0ij()()1,ijijIItNt−= suggests that all three methods describe coherence loss correctly in systems with two electronic states. Figure 3.2: The ensemble averaged absolute coherences as a function of time for the two-state model as computed by TAS (dashed red), IPSD (blue dash/dot) and TAB (light blue dots). The target exponential decay is shown with a solid black line for comparison. 3.4.3 Three-State Model In systems with more than two populated electronic states, the rate of loss of coherence between populations on two electronic states should not depend on the presence of a third electronic state. By extending to models with more than two populated electronic states, we can test whether the three proposed decoherence corrections are sensitive to the existence of additional 61 coherent populations. To gauge the existence and extent of this sensitivity, two of the electronic states in our three-state model were chosen to be parallel. For a pair of states with perfectly parallel PESs, the absolute coherence between the two electronic populations should be constant in time.59 Figure 3.3: The ensemble averaged absolute coherence between states 1 and 2 of the three-state model as a function of time computed via TAS (red line), IPSD (dark blue line), and TAB (light blue line). The coherences between parallel states 1 and 2 as computed by TAS, IPSD, and TAB are shown as a function of time in Figure 3.3. TAB correctly predicts that the coherence between states 1 and 2 is constant with time. However, the coherence between these parallel states decreases with time when computed by TAS or IPSD. The reason for this unphysical behavior is that both TAS 62 and IPSD allow collapses that change the population of one parallel state without proportionally changing the population of the other. Though the total ensemble averaged populations of the parallel states are not changed by such collapses, the ensemble averaged coherences are. TAB, on the other hand, is designed such that all collapses preserve these quantities, on average. (Very small deviations from 0.25 are visible at times before 200 au. These deviations are never larger than twice the standard error in our average over trajectories, and thus we believe these are statistically insignificant.) Another related behavior can be seen in Table 3.4, which presents the total fraction of trajectories that collapse to states 0, 1, and 2 before the end of a simulation for the three-state model. (We define a trajectory to have collapsed to a state if it its final population on that state is greater than 0.98.) In principle, no trajectories should collapse to states 1 or 2, because they are parallel and therefore populations on these states should never lose coherence with one another. This fact is indicated by the values in the target column. TAB meets this expectation; though roughly half of trajectories collapse to state 0 (indicating a physically correct loss of coherence between the population on nonparallel state 0 and populations on states 1 and 2), no trajectories completely collapse to state 1 and 2. On the other hand, collapse to states 1 and 2 is observed in TAS and IPSD, indicating, again, an unphysical loss of coherence between these parallel states. Taken together, Table 3.4 and Figure 3.3 support TAB as the preferred method for incorporating state-pairwise decoherence into nonadiabatic simulations compared to IPSD. 63 Table 3.4: Fractions of trajectories that have collapsed to individual electronic states at the end of the simulations of the three-state model. We define a trajectory to have collapsed to a state if it its final population on that state is greater than 0.98. A method that properly treats the coherence between parallel states 1 and 2 would collapse with the probabilities reported in the target column. Decoherence Correction Fraction on state 0 Fraction on state 1 Fraction on state 2 Total fraction on states 3.4.3 Five-State Model TAS 0.5005 0.1671 0.1694 0.8370 IPSD 0.5000 0.0853 0.0823 0.6676 TAB 0.4983 0.0000 0.0000 0.4983 Target 0.5000 0.0000 0.0000 0.5000 Given the definition of decoherence time in Equation (3.13), one would expect that the rate of coherence loss between pairs of states would be determined entirely by the difference between the slopes of their PESs. To test whether our three methods achieve this behavior, we turn our attention to the five-state model. As can be seen in Table 3.3, the five states in this model each have a slope that is an integer multiple of 0.010 au. Thus, four unique decay rates are expected, corresponding to pairs of states differing in slope by 0.010, 0.020, 0.030, and 0.040 au, respectively. For example, one would expect coherences , , , and to all decay with the same lifetime, because the difference in the slopes of the two associated PESs is 0.010 in each case. Figure 3.4 shows the TAB, TAS and IPSD ensemble averaged absolute coherences as a function of time. For each method, four panels are shown corresponding to the four possible force differences between paris of states. Solid black curves show the target exponential behavior. For TAB, the agreement between the coherences and the target decays is excellent (Figure 3.4a–3.4d). However, agreement is poor for TAS and IPSD [panels (e)–(h) and (i)–(l), respectively]. For TAS [panel (b)], the most notable discrepancy is that and decay much faster than and (Figure 3.4e), despite the fact their force gap suggests that these terms should decay at the same 64 1021324310432132 rate. In fact, at short time, and decay even faster than , , and , despite the fact that the latter coherences correspond to state pairs with larger force gaps. These discrepancies can be attributed to the use of the mean-field force to determine decoherence times in TAS [Equation (3.6)], as opposed to the state-pairwise decoherence times used in TAB and IPSD [Equation (3.13) ]. We see no obvious way to reconcile these errors without employing a state-pairwise description of decoherence. Figure 3.4: Ensemble averaged absolute coherences for the five-state model computed with TAB[(a)–(d)], TAS[(e)–(h)], and IPSD [bottom, (i)–(l)]. Each column of panels corresponds to the coupling between state pairs with equal force differences; in (a), (e), and (i), in (b), (f), and (j), etc. Target exponential decays in the absolute coherences [as computed with Equation (3.29)] are shown as solid black lines. IPSD coherences (Figure 3.4i–3.4l) are computed to be faster than the target exponential behavior in all cases. This is for reasons very similar to those described in the context of the failure 65 1043203142,1ii−,2ii− of IPSD to accurately describe coherence loss in the three-state model. A pairwise collapse involving a distant pair of states (e.g., 0 and 4) will be fast, but such collapses affect the coherences between the collapsing states and their nearer neighbors (e.g. 0 and 1, or 3 and 4). This results in an unintended loss of coherence between these nearer pairs at a similarly fast rate. The more sophisticated collapse procedure in TAB avoids this pitfall. 3.5 Conclusions Here we have assessed three DCE methods—TAS, IPSD and TAB—which utilize distinct means of enforcing electronic coherence loss. The primary development reported here is the use of a state-pairwise description of decoherence (in contrast to the state-wise description of decoherence employed is most decoherence correction schemes for surface hopping and Ehrenfest dynamics). TAS employs the typical state-wise description of coherence loss, while IPSD and TAB use a state-pairwise description. All three methods describe coherence loss through a stochastic collapse of the electronic wave function. IPSD and TAB differ in that the collapse in TAB modifies all electronic amplitudes simultaneously, while that in IPSD only modifies two amplitudes at a time. It is found that all three methods provide a suitable description of coherence loss in a simple two-state model, but only TAB, the more sophisticated state-pairwise scheme, provides a correct description of coherence loss in two models containing more than two states. Many interesting photochemical problems and most of the standard analytical model problems used to test nonadiabatic molecular dynamics simulations involve dynamics on only two electronic states. As such, the need for a state-pairwise description of coherence loss has not previously attracted much attention. As more theorists turn their attention to systems with large numbers of electronic states,160, 166-168 we propose that the simple one-dimensional model problems 66 presented here be applied as a standard test for correct treatment of coherence loss in systems with more than two electronic states. In addition, investigating whether the adaptation of TAB into other IT methods (e.g. surface hopping and coherent switches with decay of mixing) would improve performance in systems with more than two electronic states would be a productive direction for future research. 3.6 Acknowledgments The authors would like to thank the Air Force Office of Scientific Research (AFOSR) for support through grant FA9550-17-1-0411. We would like to thank Andy Durden, Dmitry Fedorov, Dylan Hardwick, Joe Subotnik, and Don Truhlar for fruitful conversations related to this work. 67 CHAPTER 4: DECOHERENCE-CORRECTED EHRENFEST DYNAMICS ON MANY ELECTRONIC STATES 4.1 Abstract Decoherence-corrections increase the accuracy of mixed quantum-classical nonadiabatic molecular dynamics methods, but they typically require explicit knowledge of the potential energy surfaces of all occupied electronic states. This requirement renders them impractical for applications in which large numbers of electronic states are occupied. The authors recently introduced the collapse to a block (TAB) decoherence correction [J. Chem. Phys. 152, 234105 (2020)], which incorporates a state-pairwise definition of the decoherence time to accurately describe dynamics on more than two electronic states. In this work, TAB is extended by introduction of a scheme for efficiently computing a small number of approximate eigenstates of the electronic Hamiltonian, eliminating the need for explicit knowledge of a large number of potential energy surfaces. This adaptation of TAB for dense manifolds of states (TAB-DMS) is systematically improvable by increasing the number of computed approximate eigenstates. Application to a series of one-dimensional model problems demonstrates that TAB-DMS can be accurate when even a very modest number of approximate eigenstates are computed (four in all models tested here). Comparison of TAB simulations to exact quantum dynamical simulations indicates that TAB is quite accurate so long as the decoherence correction is carefully parameterized. 68 4.2 Introduction Nonadiabatic molecular dynamics simulations have shed light on the dynamics of electronically excited systems important to a range of fields, such as molecular photochemistry,16, 130 spectroscopy,131, 132 and materials science.62, 133, 134, 136 Popular methods for this purpose include those based on trajectory surface hopping,40, 139 Ehrenfest molecular dynamics,36, 38, 137, 138 coupled Gaussian trajectory basis functions,85, 87, 146 and path integrals.150, 169 For many applications in which only a small number of states are populated, these methods have been successfully applied. However, there are many interesting chemical problems that involve very large numbers of electronic states. Developers of nonadiabatic molecular dynamics methods have begun to direct their attentions toward such systems,166-168, 170-172 with an eye towards applications such as strong field physics, atmospheric chemistry, plasmonic materials, optoelectronic materials, and radiation damage to molecules and materials. The human and computational effort required to compute a large number of PESs places severe practical constraints on the application of nonadiabatic molecular dynamics methods in these cases. Approaches like surface hopping can be adapted to large numbers of electronic states in some specific cases where the independent electron approximation is acceptable,170 but in the most general case, only Ehrenfest dynamics can be applied without explicit knowledge of the PESs of the system. Ehrenfest molecular dynamics can be applied to systems with dense manifolds of electronic states because it is formulated in terms of a single, time-dependent electronic wave function. This wave function defines a single mean-field PES along which the trajectory propagates. Unfortunately, qualitative errors are associated with propagating on this mean-field PES.40 Whereas a real quantum mechanical system that exists in a superposition of electronic states can bifurcate, yielding multiple chemical outcomes, an Ehrenfest trajectory cannot. Instead, the mean- 69 field nature of Ehrenfest molecular dynamics may result in observation of a subset of the real outcomes in some cases and wholly unphysical outcomes in others. In addition, an improper treatment of the loss of coherence between populations whose motion is governed by nonparallel PESs can lead to an incorrect prediction of transition rates and probabilities.44, 141-143 Many researchers have developed variants of the Ehrenfest approach that fix these qualitative problems,50, 153-155, 172-177 and related decoherence corrections have been carefully studied in the context of trajectory surface hopping methods, as well.45, 47, 52, 53, 55-57, 59, 152, 159, 178, 179 In the context of Ehrenfest molecular dynamics, there is an important practical consequence to these corrections. Standard decoherence-corrected Ehrenfest simulation methods require knowledge of the PESs associated with all populated states, rendering them impractical for calculations on systems with dense manifolds of electronic states. In this work we present a variant of Ehrenfest dynamics that is corrected for decoherence and does not require explicit knowledge of all electronic states. It is derived from the collapse to a block (TAB) correction recently introduced by the authors.177 TAB is based on a state-pairwise description of electronic coherence loss. The state-wise nature of the decoherence times employed in most corrections has been shown to introduce errors in both corrected Ehrenfest177 and surface hopping180 simulations of systems with more than two electronic states. TAB, however, has been shown to accurately describe coherence loss in such systems. As such, it is well suited for further development as a method for modeling dynamics in dense manifolds of electronic states. In this work, the TAB decoherence correction is modified to utilize a small number of efficiently computed approximate eigenstates in place of the full electronic eigenspectrum of the system. This is achieved by leveraging the history of the Ehrenfest time dependent wave function, in similar spirit to the multiple cloning in dense manifolds of states (MCDMS) method.172 70 In section 4.3, we present the TAB and TAB for dense manifolds of states (TAB-DMS) algorithms and the analytic models employed in our study. In section 4.4, we analyze the convergence of computed transition probabilities with respect to the number of approximate eigenstates employed in TAB-DMS. In addition, we analyze the accuracy of TAB compared to exact quantum dynamical simulations, paying particular attention to the parameterization of the decoherence time. In section 4.5, conclusions are drawn, and future prospects are discussed. 4.3 Computational Methods Three mixed quantum-classical methods will be compared in this work: standard Ehrenfest molecular dynamics, TAB, and TAB-DMS. The details of these methods and the models used to test them will be presented in three subsections. Subsection 4.3.1 describes the Ehrenfest propagation scheme common to all three mixed quantum-classical methods. In section 4.3.2 we describe the decoherence corrections used in TAB and TAB-DMS. In section 4.3.3, we present the one-dimensional models and initial conditions used for comparison of the methods. 4.3.1 Ehrenfest Propagation Scheme The Ehrenfest propagation scheme is identical to that described in detail in reference 177. Here we provide essential details only, referring the reader to the more explicit presentation of the algorithm in this previous work. Throughout this work, a single classical time step is defined as spanning the range from time to . Classical positions at time t are represented . The time-dependent electronic amplitudes, , are propagated using a second-order symplectic split operator integration scheme.162 A smaller electronic timestep, , is used for this propagation. The electronic amplitudes are propagated for the first half of the classical time step 71 0t0tt+()tR()tcet (from to ) using the electronic Hamiltonian computed at the beginning of the time step, , while they are propagated for the second half of the classical time step ( to ) using the electronic Hamiltonian computed at the end of the time step, . The classical degrees of freedom are propagated via the velocity Verlet algorithm.161 Ehrenfest forces, , are computed on the classical time steps. For standard Ehrenfest simulations, the above completely describes the algorithm. TAB and TAB-DMS include one additional step: a decoherence correction that is applied at the end of each classical time step ( ). The TAB and TAB-DMS decoherence corrections may modify the quantum amplitudes, , and the classical momentum vector, . They are described in detail in subsection 2.2. 4.3.2 Decoherence Corrections Here we describe the TAB and TAB-DMS decoherence corrections. Subsection 4.3.2.1 reviews the decoherence correction used in TAB (previously presented in reference 177). The primary difference between TAB and TAB-DMS is the eigenbasis used. Standard TAB uses the exact eigenstates of the electronic Hamiltonian, but in dense manifolds of states this becomes intractable. Subsection 4.3.2.2 presents a procedure for efficiently determining a set of approximate eigenstates that can be employed in dense manifolds of states. When employed in TAB, the resulting method is TAB-DMS. Details of the wave function collapse algorithm of both TAB and TAB-DMS are described in subsection 4.3.2.3. 72 0t02tt+()()0tHR02tt+0tt+()()0tt+HR()()0MFtFR0tt+()0tt+c()0tt+P 4.3.2.1 To-a-Block (TAB) Decoherence Correction At the heart of the TAB decoherence correction is the assumption that electronic coherences will decay exponentially at a rate governed by the difference in the gradients of the associated PESs. To enforce this behavior, the electronic wave functions of individual trajectories undergo a stochastic collapse procedure at the end of the classical time step, such that the sum of the electronic density matrices of an ensemble of trajectories yields a density matrix consistent with this assumed exponential decay. Specifically, we will define the target density matrix, , to have diagonal elements, and off diagonal elements, , . (4.1) (4.2) Here, is the fully coherent Ehrenfest density matrix, and is the decoherence time between a pair of adiabatic electronic states, i and j, defined . (4.3) Here, η indexes classical degrees of freedom, αη is a user selected decoherence parameter (akin to the width of a Gaussian wave packet in units of inverse length squared), and is an element of the force vector of adiabatic state i averaged over the classical time step, The algorithm detailed in subsection 4.3.2.3 constructs a set of density matrices, , and . (4.4) associated weights, , that sum to according to 73 dρ()0dciiiitt=+()0ijtdcijijtte−=+()ctt+ρij()2,,22avgavgijijFFћ−−=,avgiF()(),0,0,2iiavgiFtFttF++=blockaρblockaPdρ . (4.5) Here a indexes each . Each represents a pure superposition of a subset of the adiabatic electronic states that are populated in . This subset may be the entire set of adiabatic electronic states that are populated in (in which case ), it may be a single electronic state (in which case contains only a single nonzero element), or it may represent any other subset. Each is defined such that the total electronic population, the ratios between populations of populated states, and the phase angles between pairs of electronic amplitudes are conserved upon collapse from to any . Finally, the sum of the weights is unity . (4.6) Given these properties, stochastically collapsing the electronic wave functions of individual trajectories into the pure states described by with probabilities will yield ensemble density matrix , as desired. The algorithms for generating each and their associated weights for subsequent collapse are described in detail in subsection 4.3.2.3. 4.3.2.2 Approximate Electronic Eigenstates for Dense Manifolds of States (DMS) Both TAB and TAB-DMS are based on the collapse procedure described in subsection 4.3.2.1 and further detailed in subsection 4.3.2.3. The sole difference between the two is in the electronic basis used in this procedure. TAB makes use of the exact adiabatic electronic basis. This may prove intractable when one wishes to model dynamics in dense manifolds of states. With this case in mind, we introduce the following procedure for computing a smaller set of approximate 74 dblockblockaaaP=ρρblockaρblockaρ()0ctt+ρdρblockca=ρρblockaρblockaρdρblockaρ1blockaaP=blockaρblockaPdρblockaρ eigenstates. The TAB-DMS algorithm arises from simply replacing the exact eigenspectrum used in TAB by this approximate eigenspectrum. The approximate eigenstates are computed from the history of the time-dependent electronic wave function, . In this sense the approach employed here is similar in spirit to filter diagonalization.181, 182 During coherent propagation of TAB-DMS trajectories, the electronic wave function of the previous N electronic time steps is saved. More specifically, it is split into two real vectors representing the real, , and imaginary, , components of according to . (4.7) These columns are concatenated into a real, rectangular matrix, . (4.8) We note that these vectors span a subspace akin to a Krylov subspace of H acting on c, thus one would expect diagonalization in this space to provide a reasonable approximation to the eigenstates of H populated in c. The columns of X are not, in general, orthogonal. Thus, for numerical convenience we orthogonalize X by QR decomposition. The resulting rectangular matrix, Q, may then be used to transform both the expansion vector defining the electronic wave function and the electronic Hamiltonian matrix into the orthogonalized Krylov-like subspace (denoted by superscript k): and (4.9) (4.10) . 75 ()tc()Rtc()Itc()tc()()()RIttit=+ccc()()()()()()000011,,ReIeRIttNtttNttttt=+−−+−−++Xc,cc,c()()Tktt=cQcTk=HQHQ The approximate adiabatic electronic states are obtained by diagonalization of Hk and are used as the electronic basis during the wave function collapse step of the TAB-DMS algorithm, replacing the exact adiabatic eigenstates employed in TAB. 4.3.2.3 Wave Function Collapse Algorithm This subsection presents the detailed wave function collapse algorithm of TAB and TAB- DMS. The main focus of this section will be the generation of the set of . We note that there is not a unique set of weighted density matrices that satisfies Equation (4.5). Instead, we devise the below greedy decomposition algorithm, which identifies, at each iteration, the largest set of coherently coupled electronic states that can be constructed and assigns it the largest possible weight. To initialize our algorithm, is generated according to Eqs. (4.1)-(4.4) in the appropriate basis: the exact eigenstates of H in the case of TAB and the approximate eigenstates described in subsection 4.3.2.2 in the case of TAB-DMS. (Similarly, state indices i, j, k, l, m, n, p, r, and s below correspond to states in the exact and approximated adiabatic bases in the cases of TAB and TAB-DMS, respectively.) A copy of , , is created to be iteratively updated. At each iteration this matrix will be updated to represent the residual, the portion of not represented by the set of constructed so far. We now loop over the follow steps, with loop index a, to generate each and associated : 1) Elements of the matrix b are computed according to , (4.11) 76 blockaρdρdρtmpρdρblockaρblockaρblockaPtmpijijcijb= where is an element of the current , and is an element of the fully coherent density matrix associated with the trajectory, . 2) The smallest non-zero element of b, , is identified. This may be either an off-diagonal ( ) or diagonal ( ) element. In the present study, a user-defined threshold of 10-9 is used to determine whether elements of b are numerically zero. 3) Next, β, the largest possible set of coherently coupled electronic populations that includes states k and l, is generated. The procedure for generating β is as follows: a. We initiate β to include states k and l. b. For each electronic state, p, where , if , (4.12) then state p is added to β. c. For all pairs of states r and s, where and , if , (4.13) either state r or state s is randomly chosen to be removed from β. When state s has been excluded, additional pairs including s are no longer considered. 4) Nonzero elements of the coherent density matrix, , are computed , (4.14) where . All other elements of , for which and/or , are assigned to be zero. 5) The associated weight, , is computed according to 77 tmpijtmpρcij()0ctt+ρklbklkl=,pkl0 and 0plpkbbrs0rsb=blockaρ()()0,0cmnblockamnciiitttt+=+,mnblockaρmnblockaP . (4.15) Due to the accumulation of numerical error, it is possible for a to be numerically zero but with a negative sign. In this work, if , then is assigned a value of exactly zero. 6) The residual density matrix, , to be used in the subsequent iterations is updated according to . (4.16) Steps 1 through 6 are repeated, incrementing a at each iteration, until every element of b (diagonal or off-diagonal) is numerically zero. Each iteration zeros at least one element of b, thus this algorithm is guaranteed to converge in at most iterations. The result of this algorithm is a set of weights and density matrices that satisfies Equation (4.5). Now one is chosen for the electronic structure to collapse into. To this end, a random number, γ, is drawn from the uniform distribution [0,1) to determine which will be assigned as the new for subsequent propagation. If then , (4.17) . (4.18) After this assignment, is reconstructed in the basis of propagation using the new . 78 ,tmpblockklablockaklP=blockaP7100blockaP−−blockaPtmpρtmptmpblockblockaaP−ρρρ()12NN+blockaρblockaρ()0ctt+ρ and blockblockbbbabaPP()0cblockatt+ρρ()0tt+c()0ctt+ρ In order to continue propagation, the Ehrenfest energy and force vector is recomputed from the new wave function at . Then, the classical momentum vector at is scaled such that the total energy of the trajectory is conserve through the collapse. In the one-dimensional models studied in the present work, simple rescaling of the momentum vector is equivalent to scaling the momentum vector in the direction of the nonadiabatic coupling vector. Future work will address the question of momentum rescaling in multidimensional systems. Additionally, we note that there are never energy-forbidden collapses of the wave function in the present models. Consideration of such forbidden collapses is left to future work, but existing treatments of energy forbidden transitions in nonadiabatic molecular dynamics simulations would likely suffice.163, 164, 183 4.3.3 Model Problems and Initial Conditions In this work we apply TAB, TAB-DMS, traditional Ehrenfest, and numerically exact quantum mechanical simulations to six one-dimensional model problems, pictured in Figures 4.1 and 4.2. We define the electronic Hamiltonian in the diabatic basis. In each model there is a single negatively sloped diabatic state with energy , (4.19) where is a positive number and x is the nuclear coordinate. We will refer to this state as diabat 1. The remaining diabatic states form a band of parallel states that cross diabat 1. They are defined to have energies (4.20) where M is the total number of diabatic states in the model, is a positive number defining the slope of the remaining diabatic states, and δ is a positive number defining the energy spacing 79 0tt+0tt+111Hwx=−1w1M−()21for 10.05) are observe when only two approximate eigenstates are used. 84 Figure 4.3: TAB-DMS transmission probability (solid black dots) for models 9_10000 (panel a), 9_10000_split (panel b), 17_2500 (panel c), and 17_1250 (panel d) shown as a function of the number of approximate eigenstates. Solid red lines correspond to full TAB calculations. The width parameter α is set to 216.0 bohr-2 for all simulations. It is noteworthy that TAB-DMS becomes full TAB in the limit that the number of approximate eigenstates reaches the full dimension of the electronic Hilbert space, therefore convergence of TAB-DMS to the TAB limit is guaranteed. However, in all models above, the full dimension of the electronic space is either 9 or 17, considerably larger than four. This encourages optimism that TAB-DMS is a very efficient approximation to full TAB. It is also noteworthy that use of four approximate eigenstates remains sufficient even when the electronic band is split into two (model 9_10000_split, Figure 4.3b). The ability of a modest 85 number of approximate eigenstates to describe passage through multiple bands was also illustrated in our previous work on the fully quantum mechanical MCDMS method.172 4.4.2 Comparison to Exact Quantum Dynamics and Width Parameter Sensitivity Our previous work on TAB focused on analysis of the coherences rather than observable properties such as branching ratios. In fact, the development of TAB was motivated by the hypothesis that a true state-pairwise definition of the decoherence time is essential to an accurate description of coherence loss in systems with more than two electronic states. In this section, we address a more practical question; under what conditions does TAB accurately predict chemical outcomes in agreement with exact quantum dynamic simulations? To this end, we compare transmission probabilities computed via TAB with those computed via exact quantum dynamical simulations. In this context, we investigate the dependence of the TAB transmission probability on the Gaussian width parameter, α, on which the decoherence times depend. The sensitivity of other mixed quantum-classical dynamical simulations to such parameters has been extensively studied, and is generally understood to be significant.50, 55, 184-186 Below our reference numerically exact quantum dynamical simulations were performed with the initial vibrational wave function chosen to be a Gaussian wavepacket with α = 6.0 bohr-2 (with average position and momentum equal to those of the TAB simulations; described in subsection 4.3.3). One might expect the best agreement between TAB and exact quantum results to occur when the TAB simulations are also performed with α = 6.0 bohr-2. Thus, in the below discussion TAB simulations based on an α parameter of 6.0 bohr-2 will be referred to as having “intuitive parameterization.” 86 Figure 4.4: TAB ensemble (blue dots), Ehrenfest (dashed red line), and numerically exact (solid black line) transmission probabilities for models 9_10000 (panel a), 9_10000_split (panel b), 17_2500 (panel c), and 17_1250 (panel d). TAB transmission probabilities with intuitive parameterization are highlighted as light blue triangles. Shown in Figure 4.4 are the TAB transmission probabilities (dark blue dots) for models 9_10000 (panel a), 9_10000_split (panel b), 17_2500 (panel c), and 17_1250 (panel d) as a function of the value of α. The points corresponding to the intuitive parameterization are highlighted as light blue triangles. Also shown in Figure 4.4 are transmission probabilities computed from standard Ehrenfest trajectories (dashed red lines) and from numerically exact quantum dynamical simulations (solid black lines). The transmission probability for standard Ehrenfest simulations is computed by projection of the mean-field electronic wave function onto 87 diabat 1 at the end of the simulation. (Procedures for computing the transmission probabilities of the TAB, TAB-DMS, and exact quantum dynamical simulations are described in subsection 4.4.1.) In all four models, excellent agreement between TAB transmission probabilities and numerically exact transmission probabilities can be seen for larger values of α. However, simulations based on smaller values of α, including the intuitive parameterization, underestimate the exact probabilities of transmission by 0.1-0.3. We can understand this apparent discrepancy by considering the formal justification of the decoherence time in Equation (4.3). Like many who have used this definition before us, a key assumption of our method is that the coherences decay exponentially. The original derivation of the lifetime in Equation (4.3), however, was based on the assumption of Gaussian wave packets moving on non-parallel PESs. The coherence between such wave packets does not decay exponentially, but instead as a Gaussian function.44 Figure 4.5 presents an illustration of an exponential function (dashed red line) and a Gaussian function (solid black line) with the same decoherence time. As is well known, at short times the exponential function decays faster than the Gaussian function, while at longer time it decays more slowly. Thus, using the intuitive parameterization may result in the coherences being destroyed too quickly in the case of narrowly avoided crossings (like those studied here). This slows population transfer and biases the simulations towards adiabatic propagation (low transmission probability). In our TAB simulations this problem is alleviated by increasing α. 88 Figure 4.5. Illustration comparing exponential (red dashed) and Gaussian (black) approximations to the absolute quantum coherence as a function of time, with equivalent time constants. Turning our attention back to Figure 4.4, model 17_1250 (panel d) is particularly interesting because there is a significant different between the exact quantum mechanical results and the Ehrenfest results. Note that given a large enough value of α, TAB provides a better approximation to the exact results than does Ehrenfest. Of the models discussed so far, model 17_1250 has the densest band of states, suggesting that overcoherence in traditional Ehrenfest trajectories causes more significant errors in the predicted transmission probabilities in system with denser bands of electronic states. 89 Figure 4.6: Error in the probabilities of transmission predicted by traditional Ehrenfest trajectories (black dots), intuitively parameterized TAB ensembles (light blue triangles), TAB ensembles with an α parameter of 24.0 bohr-2 (dark blue dots), and TAB ensembles with an α parameter of 216.0 bohr-2 (red dots), shown as a function of the δ parameter of each 17 state model. To determine the extent of this correlation we consider a series of models differing only in the density of the band of states (17_5000, 17_2500, 17_1250, and 17_625). Figure 4.6 presents the error in the computed transmission probabilities relative to exact quantum dynamical simulations as a function of δ, the energy spacing between states. Remember that lower δ corresponds to denser bands of states, and the total number of states in the bands is held constant. 90 The absolute error in the Ehrenfest probability of transmission (black dots) is largest for models with denser bands of states (0.07 for the densest model, 17_625). By comparison, TAB simulations with an α parameter of 216.0 bohr-2 (red dots) consistently underestimate the exact quantum dynamics results by only ~0.02, independent of δ. As above, strong dependence on the choice of α is observed. Development of a systematic scheme for its selection will be the goal of a future work. 4.5 Conclusions Here we have introduced an extension of the TAB mixed quantum-classical scheme for dense manifolds of electronic states, TAB-DMS. TAB-DMS is a variant of Ehrenfest molecular dynamics that incorporates a decoherence correction whereby population stochastically collapses into a linear combination of approximate eigenstates. These approximate eigenstates are determined from the history of the time-dependent mean-field electronic wave function. This efficient procedure eliminates the need to explicitly diagonalize the electronic Hamiltonian, and therefore is particularly well suited for systems with a large number of electronic states. Increasing the number of approximate eigenstates allows convergence to the full TAB (exact eigenstate) limit, and by application to several one-dimension models we have shown that excellent results can be achieved with a modest number of eigenstates (only four in all cases studied here). TAB-DMS maintains the primary advantages of TAB; the incorporation of a state-pairwise description of decoherence enables the decay of the coherences between pairs of electronic states to be accurately described in systems with more than two electronic states. Additionally, the accuracies of the probabilities of nonadiabatic transmission through a band of parallel states was evaluated by comparison of TAB results to exact quantum dynamical 91 simulations. TAB is found to be accurate across a range of densities of states, so long as the α (Gaussian width) parameter is carefully chosen. It is found that the most intuitive choice of α does not provide the most accurate results, owing to the assumption of exponential decay of the coherences. These results suggest directions for development of a systematic scheme for parameterization. Ehrenfest methods based on on-the-fly ab initio treatments of the electronic structure are powerful tools for modeling coupled electron-nuclear dynamics.187, 188 The primary goal of this work is to develop a scheme that can be linked to ab initio electronic structure methods, enabling even more accurate and efficient simulation of nonadiabatic molecular dynamics in dense manifolds of electronic states from first principles. Toward this end, coupling to high-performance graphics processing unit time-dependent configuration interaction tools189 is underway. 4.6 Acknowledgments The authors would like to thank the Air Force Office of Scientific Research (AFOSR) for support through grant FA9550-17-1-0411. The authors would like to thank Andy Durden, Dmitry Fedorov, Dylan Hardwick, and Srini Iyengar for helpful discussions. 92 CHAPTER 5: CONCLUSIONS AND OUTLOOK In this work, we have further demonstrated that insights towards semiconductor nanomaterial photophysics can be obtained using ideas from molecular photochemistry and NAMD simulations, and we have designed and implemented a decoherence-corrected Ehrenfest method capable of simulating photophysical processes over many electronic states. Our computational investigation of semiconductor nanomaterial photophysics was focused on determining the energetic accessibility of pathways for efficient NRR in CsPbBr3 perovskites after optical excitation. Though lead-halide perovskites are known to have a low propensity for undergoing NRR, research to understand what structural properties of lead-halide perovskites are the underlying cause of this observation is ongoing, because this insight could inform the design of yet more efficient optoelectronic materials.5-12 In the field of molecular photochemistry, CIs are known to provide pathways for efficient nonradiative decay.13-17 Our goal was to gauge the energetic accessibility of CIs in large CsPbBr3 nanoparticles and bulk CsPbBr3 after optical excitation to determine whether they could provide pathways for efficient NRR. Described in chapter 2, the first step towards achieving this goal was identifying pathways for NRR in a molecule sized cluster model of CsPbBr3 using NAMD simulations. Key geometries from the dynamics were then used to inform time-independent optimizations of structures corresponding to the lowest energy CIs in both the cluster model and a larger nanoparticle model of the CsPbBr3 surface. Energies of the optimized MECIs were corrected for both dynamic electron correlation and spin–orbit coupling and were found to be well above the optical band gap of bulk CsPbBr3, suggesting these observed CIs do not provide pathways for efficient NRR in bulk CsPbBr3 or large nanoparticles. Further investigation of the PESs 93 surrounding the intersections suggested that the ionic nature of the bonds in CsPbBr3 contributed to the high energy of the observed intersections. Thus, our results suggest investigation of ionic semiconductor nanomaterials would be a promising direction towards the design of yet more efficient optoelectronic devices. It is important to note the characterization of CIs in the nanoparticle model was made practical through the employment of GPU accelerated two-step CASCI methods.28-33 Though largely consumer driven, the increase in power of GPUs over the past several decades is a boon to computational sciences as a whole due to a general increase in efficiency from these hardware developments. For computing electronic excited state energies for nanoscale systems, two-step CASCI methods have also been very beneficial to the field material science due to greater numerical stability, size-intensive electronic excitation energies, and generally lower cost than similar MC-SCF methods.34, 35 However, the combined benefits of computational tractability are yet insufficient to model photophysical processes occurring over many electronic states using PES refined IT-NAMD dynamics. Thus, motivating the subsequently described developments of chapters 3 and 4. In chapters 3 and 4, we have introduced the novel decoherence corrected Ehrenfest methods, TAB and TAB-DMS. The TAB decoherence correction enforces electronic coherence loss through a stochastic collapse of the electronic wavefunction into a coherent superposition of a subset of the populated electronic states at the end of classical simulation time steps. The probability the wavefunction collapses into a particular superposition is computed using state- pairwise decoherence times, such that ensemble averaged coherences of populations propagating on near-parallel PESs will be preserved for a long time, while coherences decay exponentially at rates exclusively dependent on the difference in the involved PES gradients. The importance of 94 the state-pairwise formulation of coherence loss, and the simultaneous treatment of all electronic amplitudes, was demonstrated in chapter 3. Of the three methods tested, one which assigned decoherence times in a state-wise fashion and another which sequentially treated coherence loss between pairs of electronic states using state-pairwise decoherence times, only TAB ensembles correctly predicted absolute coherences over time in agreement with the target functions for the suite of simple, multistate models employed. The TAB-DMS methodology is an extension of the TAB scheme for simulations of photophysical processes over many electronic states and was introduced and tested in chapter 4. In the TAB-DMS methodology, the computational cost of computing an eigenbasis with a rank matching that of the full electronic Hilbert space is avoided by using the history of the time- dependent Ehrenfest wavefunction to generate a significantly smaller set of approximate eigenstates. This approximate eigenbasis is then used to enforce electronic coherence loss between populations propagating on nonparallel approximate PESs. Specifically, electronic coherence loss is enforced by applying the wavefunction collapse algorithm of TAB to the density matrix associated with the Ehrenfest wavefunction computed in the approximate basis. The direct application of the TAB wavefunction collapse algorithm to this density matrix is possible because the algorithm can be applied to any electronic density matrix associated with a pure electronic state, so long as all values necessary for computing the state-pairwise decoherence times can be obtained. Conveniently, approximate electronic state populations and approximate PES gradients can be readily obtained from data structures present in coherent Ehrenfest propagation in the TAB- DMS formalism. Comparing the ensemble averaged probabilities of transmission computed using simple one-dimensional models, we demonstrated that TAB-DMS recovers the accuracy of TAB 95 dynamics when the number of approximate eigenstates used is only a modest fraction of the full dimensionality of the electronic Hilbert space. The accuracy of TAB dynamics was also evaluated through comparison to numerically exact quantum dynamics. Our results have shown that well- parameterized TAB simulations predicted probabilities of transmission in excellent agreement to those computed using exact dynamics. Additionally, TAB was shown to outperform traditional Ehrenfest dynamics for models in which the trajectory passed through the densest bands of parallel electronic states. Suggesting pairing TAB-DMS with efficient GPU accelerated ab-initio electronic structure theory and electronic propagation schemes is a promising for future IT-NAMD simulations of photophysical processes over many electronic states. Looking to other future developments, we hypothesize the already quite good accuracy of well parameterized TAB simulations can be further improved by adding some flexibility to the functional form of electronic coherence loss. To summarize from chapter 4, TAB ensembles consistently underestimated exactly computed probabilities of transmission for each of our simple one-dimensional models, though only slightly for the best parameterized ensembles. This was discussed to be the result of employing the commonplace assumption that electronic coherences decay in an exponential fashion, but using a lifetime corresponding to a Gaussian decay, which was derived under the assumption Gaussian wave packets were moving along nonparallel PESs.44 More specific to the trajectory dynamics, is the possibility that TAB coherences were being destroyed too quickly as trajectories approached weakly avoided PES crossings, slowing the population transfer as a result of these assumptions. To amend the trajectory dynamics, we suggest altering the functional form of coherence loss. The present focus is modifying the form of the state- pairwise decoherence time in order to preserve the numerical advantages of using exponential functions to model coherence loss. 96 As previously discussed, it is known that electronic coherences only decay when the corresponding populations evolve on nonparallel PESs.44, 59 Thus, every term in an altered TAB state-pairwise decoherence time must be inversely proportional to the difference in the corresponding PES gradients. However, this does not mean additional functional dependencies cannot be included in this lifetime, as seen in the decoherence times of Subotnik et. al in Augmented Fewest Switches Surface Hopping (AFSSH).56 Though, for present considerations we will forgo looking further into employing the AFSSH decoherence times in TAB dynamics, even though they may yield significant improvements, since the AFSSH lifetimes are not presently implementable in the TAB-DMS framework. Instead, we look to the energy-based schemes of Granucci et. al. and the coherent switches with decay of mixing method by Truhlar et. al. for inspiration.49, 52, 53 Although these methods do not allow for the preservation of coherences between populations propagating on parallel nondegenerate PESs, they have seen regular and successful use modeling nonadiabatic dynamics. One particular property they have, rates of coherence loss between populations on degenerate PESs being zero, we hypothesize could be useful to incorporate in TAB dynamics to help alleviate the adiabatic bias in population transfer. More specifically, as an energy-based decoherence corrected trajectory approaches a weakly avoided PES crossing, the rate of coherence loss is slowed, allowing for the coherence to build up more quickly. Subsequently, causing greater population transfer between adiabatic states. Therefore, we propose investigation of including a dependence on the energy difference of the involved PESs in the TAB state-pairwise decoherence time, in addition to the difference in the involved PES gradients. A potential form of this type of state-pairwise decoherence time may be , (5.1) 97 2,,avgavgijijijEEFF−−− where Ei and Ej are the energies of adiabatic states i and j. This form of the state-pairwise decoherence time not only preserves coherences between parallel PESs, but will also slow rates of coherence loss as trajectories approach near-degenerate regions around PES crossings. This may allow for greater adiabatic population transfer in TAB simulations, improving the accuracy. To all readers, Thank you. Reading this is no small investment of time, and I hope the information presented here will be helpful in your future endeavors. Michael Paul Esch 98 APPENDIX 99 Atomic Width Optimization Scheme for AIMS TBFs The Full Multiple Spawning algorithm, which is referred to as AIMS when the electronic PESs are calculated on the fly using ab-initio electronic structure theory, has nuclear trajectory basis functions (TBFs) that are a product of 3n one dimensional Gaussians, where n is the number atoms in the system, in Cartesian coordinate space with the constraints that widths of Gaussians for a particular atom are the same in the x, y and z directions and that atoms of the same element have identical widths. To the best of our knowledge there are no known atomic widths for Cs, Br and Pb published for perovskite systems, so we optimized these widths such that the overlap of the cluster model FMS nuclear wavefunction with zero momentum and ground state vibrational wavefunction at the minimum energy geometry calculated at the harmonic approximation was maximized. This section describes in detail how the overlap of the FMS nuclear wavefunction and the harmonic oscillator wavefunction is calculated without losing generality so this width optimization method can be applied to other systems. In calculating this overlap, we will first address the ground state harmonic oscillator wavefunction (A.1) where is the ith normal coordinate, is the normalization constant and the Gaussian widths are defined as (A.2) where is the ith vibrational frequency and is the corresponding reduced mass. As seen in Equation (A.1), the harmonic oscillator wavefunction can be written as a product of Gaussians like the FMS nuclear wavefunction. However, to calculate the overlap of these two wavefunctions they need to be cast into the same coordinate space. Thus, we will transform the FMS nuclear 100 ()213621236,,...,iinqHOnHOiqqqNe−−−=iqHONiiiiћ=ii wavefunction into normal coordinate space to calculate its overlap with the harmonic oscillator wavefunction. Before addressing the transformation of the FMS nuclear wavefunction, we will substitute Equation (A.2) into Equation (A.1) and equivalently rewrite the harmonic oscillator wavefunction in matrix notation as where the vector q (A.3) (A.4) contains all normal coordinates, the matrix W (A.5) contains all normal mode vibrational frequencies, and the matrix (normal mode mass matrix) (A.6) contains the reduced masses corresponding to each normal mode vibration, in order to simplify the calculation of the overlap of the wavefunctions. Each in Equation (A.6) is calculated as 101 ()11222TTNMNMHOHONe−=WqMMqq1236nqqq−=q123600000000n−=WNMM123600000000NMn−=Mi (A.7) where is the eigenvector of the mass weighted Hessian corresponding to the ith normal vibrational mode and is the 3n rank mass matrix in Cartesian Coordinates. Since both W and are diagonal, Equation (A.3) can be simplified to with a normalization constant of (A.8) (A.9) where is the determinant of the product of matrices W and . The FMS nuclear wavefunction with zero momentum in Cartesian Coordinates, atomic units and written in matrix notation is defined as (A.10) where is a vector containing the displacement of each atom from its equilibrium position (A.11) and is the diagonal matrix 102 TiiciUU=MiUcMNMM()2TNMHONe−=WMqqq()436detNMHOnN−=WM()detNMWMNMM()()()000TFNFNNe−−−−=xxAxxxx()0−xx()1101101100330330330nnnnnnxxyyzzxxyyzz−−−−=−−−xxA (A.12) where is the Gaussian width for the ith atom. The normalization constant of the FMS nuclear wavefunction in Cartesian coordinates is calculated to be (A.13) The FMS nuclear wavefunction is first transformed into normal coordinate space using the 3n-6 eigenvectors of the mass weighted Hessian corresponding to the normal vibrational modes. After the transformation, is (A.14) where U is the 3n by 3n-6 matrix containing the 3n-6 eigenvectors of the mass weighted Hessian corresponding to the normal vibrational modes. In normal coordinate space, the normalization constant of the FMS nuclear wavefunction with zero momentum is (A.15) With both wavefunctions expressed in normal coordinates, the overlap can be written as (A.16) Using standard Gaussian integral definitions 103 111333000000000000000000000000000000000000nnn=Ai()4det2FNnN=AFN()()TTFNFNNe−=qUAUqq()436det2TFNnN−=UAU()2TTNMFNHOSNNed−+=WMqUAUqAq (6.17) To calculate the widths of the Gaussians that maximize , a Nelder-Mead (simplex) minimization of the negative natural logarithm of Equation (6.17) was performed. To calculate cluster model widths for Cs, Br and Pb we used the information from the B3LYP/SBKJC Polarized(p,2d) harmonic frequency calculation performed at the optimized ground state minimum energy geometry. Although this method for optimizing atomic widths for FMS TBFs is completely general there is one limitation to this method, which is that the number of unique elements in a system cannot exceed the number of normal mode vibrational frequencies. When the number of unique atoms is greater than the number of normal mode frequencies, like in a heteronuclear diatomic, one of the atomic widths becomes a free variable and the optimization is not able to converge to a meaningful solution. 104 ()1236det2nFNHOTNMSNN−=+AWMUAU()SA B3LYP Optimized Cluster Model FC geometry Figure A1: B3LYP/SBKJC Polarized(p, 2d) optimized ground state minimum energy geometry. The absolute energy of this structure at the optimization level of theory is -84.9289162 Hartree. Table A1: Absolute energies of singlet electronic states calculated at the B3YLP optimized FC geometry. The singlet state energies were calculated at the EOM-CCSD/SBKJC Polarized(p, 2d), SA5-CASSCF(2/5)/SBKJC VDZ ECP, MS5-CASPT2-c19/SBKJC Polarized(p, 2d), MS5- CASPT2-c0/SBKJC Polarized(p, 2d), MRCI/SBKJC Polarized(p, 2d) with the same active space and number of electronic states and the Davidson correction, TD-PBE0/SBKJC VDZ ECP and FOMO-CASCI(2/5)/SBKJC VDZ ECP levels of theory. Level of Theory EOM-CCSD SA-CASSCF MS-CASPT2-c19 MS-CASPT2-c0 MRCI TD-PBE0 FOMO-CASCI S0 Energy -83.9545752 -83.4237955 -83.5166954 -83.9039095 -83.9427159 -84.9999461 -83.4011538 S1 Energy -83.7864772 -83.2325636 -83.3333620 -83.7385824 -83.7662422 -84.8354565 -83.2222163 S2 Energy -83.7832002 -83.2284427 -83.3298928 -83.7344416 -83.7619967 -84.8317651 - S3 Energy -83.7831240 -83.2283287 -83.3297927 -83.7343283 -83.7618810 -84.8316987 - S4 Energy -83.7830434 -83.2282150 -83.3296932 -83.7342209 -83.7617704 -84.8316010 - Geometry of the B3LYP optimized FC point listed in Cartesian Coordinates: Cs 0.3391847229 0.4348353933 0.3346285579 Cs 5.0104872469 5.0763602944 0.3445116108 Cs 0.4042902225 5.0354636667 5.0746523923 Cs 5.0677370465 0.3372899343 4.9458226215 Br 2.7006904235 -0.3546407121 2.6323253968 Br 2.6574089744 2.7743380598 -0.4025851082 Br -0.3673765795 2.7261433353 2.7261060571 Br 5.7769314565 2.7054518213 2.6402300694 Br 2.7234656718 5.7897155150 2.7360796194 Br 2.7401197276 2.6713958829 5.7424915120 Pb 2.7048180433 2.7178863174 2.6820480609 105 Figure A2: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory. the B3YLP optimized FC geometry calculated at Figure A3: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5- CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. the B3LYP optimized FC geometry calculated at Figure A4: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO- CASCI(2/5)/SBKJC VDZ ECP level of theory. the B3LYP optimized FC geometry calculated at 106 Figure A5: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the B3LYP optimized FC geometry calculated at the SA5-CASSCF(2/5)/cc- PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. Table A2: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the B3LYP optimized FC geometry. n 0 1 2 3 4 Sn Energy -2764.85630122 -2764.66721318 -2764.66317230 -2764.66306862 -2764.66296763 Table A3: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the B3YLP optimized FC geometry in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.85630123 -2764.66788653 -2764.66788652 -2764.66788652 -2764.66721323 -2764.66439003 -2764.66427062 -2764.66424769 -2764.66421327 -2764.66356607 -2764.66355949 -2764.66350685 -2764.66346568 -2764.66345725 -2764.66301098 -2764.66292316 -2764.66283747 107 SA-CASSCF Optimized Cluster Model FC S1 Minimum Figure A6: Local S1 minimum to the FC geometry optimized at the SA5-CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. Table A4: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized S1 minimum local to the FC region calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.3953516 S1 Energy -83.2813017 S2 Energy -83.2790737 S3 Energy -83.2778992 S4 Energy -83.2773458 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized S1 minimum local to the FC region listed in Cartesian Coordinates: Cs 0.0971386674 0.2017733212 0.1061290133 Cs 5.2470834595 5.3177387009 0.1160259420 Cs 0.1613203652 5.2781483057 5.3376046927 Cs 5.3127454609 0.0898944399 5.1915067792 Br 2.6966922252 -0.1642252934 2.6230301281 Br 2.6622995954 2.7684252129 -0.2180034133 Br -0.1744982994 2.7231563414 2.7093103443 Br 5.5871097722 2.7086034480 2.6306489604 Br 2.7231199382 5.5973093641 2.7179756760 Br 2.7416947310 2.6711829753 5.5430827091 Pb 2.7030510756 2.7222326870 2.6990012483 MS-CASPT2-c19 Optimized Cluster Model FC S1 Minimum Figure A7: Local S1 minimum to the FC geometry optimized at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. 108 Table A5: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized S1 minimum local to the FC region. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)- c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.5004631 -83.8686083 -83.9047122 S2 Energy -83.3735181 -83.7462178 -83.7735672 S1 Energy -83.3882014 -83.7609598 -83.7881183 S3 Energy -83.3704435 -83.7427095 -83.7700476 S4 Energy -83.3675522 -83.7393938 -83.7667075 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d) optimized S1 minimum local to the FC region listed in Cartesian Coordinates: Cs 0.2661583783 0.3853662995 0.2532018246 Cs 5.0761321293 5.1128435414 0.2892357464 Cs 0.0587730439 5.3798236667 5.4430839043 Cs 5.1767312279 0.2746928814 5.0350809893 Br 2.7295218631 -0.1033006836 2.6184240904 Br 2.6875855237 2.7333667225 -0.6267251060 Br -0.1122242681 2.6920988044 2.7475051555 Br 5.7728602636 2.6455356750 2.6183018465 Br 2.7305559435 5.4499994346 2.7410924762 Br 2.7461774878 2.6365453499 5.5220160668 Pb 2.6254854931 2.7072676857 2.8150935396 SA-CASSCF Optimized Cluster Model BrD CI Figure A8: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A6: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.30194549 S2 Energy -83.26519353 S4 Energy -83.24584052 S3 Energy -83.2505347 S1 Energy -83.3010363 109 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI listed in Cartesian Coordinates: Cs 0.0321902683 0.4952567584 0.8379192116 Cs 3.1089782203 6.4040538986 -1.0564520499 Cs 0.3598703877 4.6366990171 5.0836136677 Cs 2.4424059965 -0.3831075599 5.6183057001 Br 3.3128301287 0.0916803364 2.2035473672 Br 2.2003486952 3.0054978398 -0.5369632299 Br -0.6357845834 1.2067448471 4.3272513902 Br 12.2224298000 3.0240014049 3.1997258510 Br 2.4303044738 5.8239991018 2.3774013714 Br 3.5840020070 2.9703406454 5.1320718980 Pb 3.7773241198 2.9810239219 2.1974006851 Figure A9: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. Figure A10: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. 110 Table A7: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI. n 0 1 2 3 4 Sn Energy -2764.75302087 -2764.74503710 -2764.71855776 -2764.70359145 -2764.69916923 Table A8: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.75302209 -2764.74541376 -2764.74541376 -2764.74541374 -2764.74538659 -2764.71846691 -2764.71846690 -2764.71846690 -2764.71846536 -2764.70355702 -2764.70355701 -2764.70355701 -2764.70355181 -2764.69895170 -2764.69895169 -2764.69895168 -2764.69895061 MS-CASPT2-c19 Optimized Cluster Model BrD CI Figure A11: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. 111 Table A9: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized BrD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.4056410 -83.7511516 -83.7949362 S1 Energy -83.4049966 -83.7469968 -83.7803357 S2 Energy -83.3611537 -83.6890124 -83.7244948 S3 Energy -83.3519300 -83.6804193 -83.7166380 S4 Energy -83.3471433 -83.6761647 -83.7119625 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized BrD CI listed in Cartesian Coordinates: Cs 0.0054586516 0.4665022097 0.8887315056 Cs 2.3416136410 6.1904779780 -0.9193687206 Cs 0.6210920643 4.5831457429 5.0629361979 Cs 2.2492351760 -0.0802150154 5.3412985330 Br 3.1981935397 0.0429389800 2.0823563656 Br 1.7603864588 2.9378228328 -0.5228649962 Br -0.7731083498 1.3683318031 4.2095654578 Br 13.9987247045 2.9269262078 3.7995357303 Br 2.3682449187 5.9844204407 2.4005683310 Br 3.5970844593 2.9224780877 4.9782401009 Pb 3.4679805885 2.9133557581 2.0628968290 Figure A12: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized BrD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. 112 FOMO-CASCI Optimized Cluster Model BrD CI Figure A13: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the BrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A10: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -83.28499080 S1 Energy -83.28473388 Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized BrD CI listed in Cartesian Coordinates: Cs 0.3829312634 0.0464543236 0.3657236566 Cs 2.1911747705 6.3801302067 -0.9437791916 Cs 0.7528598343 5.0315635633 5.4940751673 Cs 2.1662559754 -0.3163030109 5.6340151774 Br 3.3542257656 -0.1821849452 2.3090994453 Br 2.1956272021 2.8776776674 -0.6843694062 Br -0.0775016061 1.8902505343 3.6402220608 Br 12.5461410642 2.9177401480 3.3178411899 Br 2.4592023954 5.8996769765 2.5303307001 Br 3.6057134804 2.8955678796 5.3486420318 Pb 3.2582666309 2.8156088152 2.3720250734 Figure A14: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. 113 MS-CASPT2-c19 Cluster Model BrD S1 Minimum Figure A15: Local S1 minimum in the cluster model BrD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. Table A11: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized BrD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.4118839 -83.7606907 -83.8017373 S1 Energy -83.4087214 -83.7502250 -83.7858919 S2 Energy -83.3593941 -83.6913689 -83.7265083 S3 Energy -83.3466760 -83.6820707 -83.7171755 S4 Energy -83.3437247 -83.6794888 -83.7136609 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized BrD PES region S1 minimum listed in Cartesian Coordinates: Cs 0.2280743487 0.6532918709 0.9717299195 Cs 3.0396076554 6.1751714886 -0.8468714175 Cs 0.5608805781 4.4555035207 4.9048528521 Cs 2.3800373141 -0.1152078890 5.3802598269 Br 3.3339633706 0.0317389681 2.1905359759 Br 2.2103423973 2.9823054800 -0.5004639605 Br -0.7219741122 1.2239287629 4.3236787840 Br 12.1217863589 3.0319835960 3.1814803305 Br 2.3913705934 5.9087400307 2.4076760471 Br 3.5832607321 2.9636408811 5.1503323610 Pb 3.6875502194 2.9450935383 2.2206104992 114 SA-CASSCF Optimized Cluster Model CsD CI Figure A16: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A12: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.2692858 S1 Energy -83.2689911 S2 Energy -83.2267461 S3 Energy -83.2267415 S4 Energy -83.2227539 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI listed in Cartesian Coordinates: Cs -0.2655524567 0.7953755307 -0.3518873715 Cs 9.1827768944 -3.8163187602 8.8291500029 Cs 4.6437835922 5.6910280158 -0.3335022994 Cs -0.2852081424 5.7243879758 4.7938409881 Br 2.2202325237 0.3241448163 2.2142925175 Br 2.2037280229 3.2317454235 -0.6863937996 Br -0.6558706143 3.1625922626 2.3269879973 Br 5.0941406250 3.1852242273 2.2199149547 Br 2.2788877264 6.0810095144 2.3367323179 Br 2.2361219696 3.1723185374 5.1695089384 Pb 2.3543409035 3.0640816091 2.4206670051 115 Figure A17: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. Figure A18: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized BrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. Table A13: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI. n 0 1 2 3 4 Sn Energy -2764.70835297 -2764.69915348 -2764.66556017 -2764.66555472 -2764.66315390 116 Table A14: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.70835306 -2764.70835127 -2764.70835127 -2764.70835127 -2764.69915348 -2764.66617592 -2764.66617578 -2764.66617577 -2764.66617514 -2764.66508682 -2764.66508682 -2764.66508667 -2764.66508666 -2764.66302219 -2764.66302218 -2764.66302218 -2764.66300734 MS-CASPT2-c19 Optimized Cluster Model CsD CI Figure A19: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A15: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.3769131 -83.7401925 -83.7733696 S1 Energy -83.3761500 -83.7377623 -83.7687867 S2 Energy -83.3244520 -83.6861071 -83.7142632 S3 Energy -83.3244188 -83.6860736 -83.7142311 S4 Energy -83.3233355 -83.6850965 -83.7132660 117 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsD CI listed in Cartesian Coordinates: Cs -0.1162410615 0.9318393931 -0.2044820649 Cs 4.4492803257 5.5619492547 -0.1366136901 Cs -0.1229164940 5.6287639365 4.6396315296 Cs 9.1924886546 -3.8278186979 8.8400738536 Br 2.1956621182 0.3244229539 2.2792962287 Br 2.1561592290 3.2648349007 -1.1086475061 Br -0.4759506880 3.1827424013 2.3142861975 Br 4.9467198693 3.1754755135 2.2559054649 Br 2.2320549249 6.2684893277 2.3308492318 Br 2.2304489795 3.1455408487 5.1752821634 Pb 2.3196752149 2.9593493773 2.5537290045 Figure A20: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. FOMO-CASCI Optimized Cluster Model CsD CI Figure A21: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the CsD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. 118 Table A16: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -83.25582985 S1 Energy -83.25500648 Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsD CI listed in Cartesian Coordinates: Cs -0.2920571376 0.8247060368 -0.2933536325 Cs 4.6213824105 5.7136327271 -0.2806982565 Cs -0.2589590024 5.6946843299 4.7299222163 Cs 9.1649632627 -3.7982755276 8.8122637174 Br 2.2390092696 0.2702110122 2.2329314732 Br 2.2281096115 3.2033503677 -0.8028496993 Br -0.7965197465 3.1470442879 2.3410984635 Br 5.1508485014 3.1753052390 2.2502792488 Br 2.2948739396 6.2027078272 2.3404172172 Br 2.2704072209 3.1476555090 5.1840867899 Pb 2.3853227145 3.0345673458 2.4252137116 Figure A22: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. 119 SA-CASSCF Optimized Cluster Model CsD S1 Minimum Figure A23: Local S1 minimum in the cluster model CsD PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. Table A17: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsD PES region S1 minimum calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.3864385 S1 Energy -83.2838949 S2 Energy -83.2713023 S3 Energy -83.2699703 S4 Energy -83.2697453 Geometry of the optimized SA5-CASSCF(2/5)/SBKJC VDZ ECP CsD PES region S1 minimum listed in Cartesian Coordinates: Cs 0.0733355479 0.2472495692 0.0975225304 Cs 5.1146817885 5.3717471293 0.0609506972 Cs 0.0878063852 5.3105072519 5.1923561241 Cs 5.4936186286 -0.0185316258 5.3787971430 Br 2.6333068635 -0.0433504922 2.5946312788 Br 2.5700462919 2.8308194979 -0.3050688764 Br -0.3006425922 2.7975478031 2.6268350484 Br 5.4690916194 2.8331301672 2.5750046430 Br 2.5784515410 5.7238764064 2.6079339918 Br 2.6414874146 2.7985835794 5.4590297807 Pb 2.6761975561 2.7640098664 2.6513189567 MS-CASPT2-c19 Optimized Cluster Model CsD S1 Minimum Figure A24: Local S1 minimum in the cluster model CsD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. 120 Table A18: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.4173879 -83.7786336 -83.8154363 S1 Energy -83.3789549 -83.7413036 -83.7698338 S3 Energy -83.3303841 -83.6926654 -83.7210739 S2 Energy -83.3304365 -83.6927185 -83.7211220 S4 Energy -83.3290902 -83.6917143 -83.7201122 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsD PES region S1 minimum listed in Cartesian Coordinates: Cs -0.0813495902 0.8116028735 0.0759564691 Cs 4.6703711319 5.5078542485 0.0246099181 Cs 0.0477736882 5.3276420897 4.6540670923 Cs 7.7275507264 -2.3513360807 7.4425881392 Br 2.4319079034 0.1384430363 2.3254001720 Br 2.3706602923 3.0758128736 -0.3838899525 Br -0.9941952328 3.0610850660 2.3848982562 Br 5.3214431661 3.0455127148 2.3533116031 Br 2.4071832653 6.1758304880 2.3969341086 Br 2.4430690905 3.0014171893 5.1881976841 Pb 2.6629795846 2.8217420190 2.4773079777 SA-CASSCF Optimized Cluster Model DBrS CI Figure A25: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the DBrS region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A19: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.33030808 S2 Energy -83.25394366 S1 Energy -83.33028948 S3 Energy -83.25045512 S4 Energy -83.21819391 121 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI listed in Cartesian Coordinates: Cs 0.5998160779 0.6381183756 0.7853763804 Cs 4.7488583550 5.8808054060 -0.3556763517 Cs 0.7242967082 4.5894856820 4.9999006104 Cs 2.0222949488 -0.4346440667 5.6961203877 Br 3.6046577005 -0.3196154790 2.5427477536 Br 3.0495460311 2.8453132721 -0.5399521423 Br -0.7445217183 1.3494904550 4.1300368557 Br 8.1936332461 4.3757125654 1.0632725173 Br 3.0532405538 5.9258630622 2.6820222302 Br 3.6896684393 2.6843036647 5.7237888226 Pb 3.8734100124 2.7213608795 2.6561812490 Figure A26: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. Figure A27: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. 122 Table A20: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI. n 0 1 2 3 4 Sn Energy -2764.79272277 -2764.78099350 -2764.70379673 -2764.70023991 -2764.67764051 Table A21: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized DBrS CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.79296161 -2764.78709406 -2764.78709325 -2764.78683860 -2764.78683179 -2764.71102411 -2764.71102008 -2764.71092957 -2764.71090126 -2764.68836879 -2764.68836868 -2764.68819948 -2764.68816497 -2764.67802139 -2764.67802134 -2764.67801403 -2764.67706670 MS-CASPT2-c19 Optimized Cluster Model DBrS CI Figure A28: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the DBrS region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. 123 Table A22: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized DBrS CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.4344013 -83.7756854 -83.8240352 S1 Energy -83.4337358 -83.7646051 -83.7957141 S2 Energy -83.3615560 -83.6910275 -83.7252978 S3 Energy -83.3554430 -83.6851918 -83.7196985 S4 Energy -83.3158349 -83.6443528 -83.6804743 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized DBrS CI listed in Cartesian Coordinates: Cs 0.7078229310 0.7111119627 0.9548000825 Cs 4.5826338891 5.8354158265 -0.2908160733 Cs 0.8376502450 4.4376351270 4.8709672432 Cs 2.0133341445 -0.2477587367 5.5273243498 Br 3.5387930347 -0.3456195752 2.5945788349 Br 3.0042468405 2.8476569760 -0.7916154301 Br -0.7430128757 1.3422332692 4.1471241499 Br 8.3112949766 4.4029593243 1.0931333170 Br 3.0841272537 6.0011774778 2.6982793320 Br 3.6534180589 2.6422626017 5.7738091351 Pb 3.8245919107 2.6291195289 2.8062343442 Figure A29: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized DBrS CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. 124 SA-CASSCF Optimized Cluster Model CsBrD CI Figure A30: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A23: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.2699688 S1 Energy -83.2695015 S2 Energy -83.2275289 S3 Energy -83.2271096 S4 Energy -83.2236320 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI listed in Cartesian Coordinates: Cs -0.7424645932 5.3093834728 -0.6707208480 Cs 4.5882097537 5.5371669692 -0.6459403147 Cs -0.2792133684 6.5728030706 5.0034081245 Cs 7.0375844100 -3.9807732165 6.3513510558 Br 9.7424187310 -12.5314666179 7.3262189940 Br 1.8978184815 5.8336378228 -2.8707926158 Br -1.3415929453 4.0403890137 2.6923602461 Br 4.4510927769 2.8904970656 1.8697018280 Br 1.9476432500 6.0492491620 2.0168281026 Br 2.2291211979 4.0484268160 5.5593423163 Pb 1.8318253949 3.2938625503 2.9178183461 125 Figure A31: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization. Figure A32: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. Table A24: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI. n 0 1 2 3 4 Sn Energy -2764.72083600 -2764.71594157 -2764.67797186 -2764.67767413 -2764.67507184 126 Table A25: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.72083627 -2764.72073689 -2764.72073689 -2764.72073687 -2764.71594253 -2764.67846935 -2764.67846862 -2764.67846757 -2764.67845352 -2764.67733494 -2764.67733456 -2764.67733444 -2764.67732967 -2764.67523997 -2764.67523994 -2764.67523987 -2764.67493692 MS-CASPT2-c19 Optimized Cluster Model CsBrD CI Figure A33: Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A26: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsBrD CI. Energies are calculated at the MS5- CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.3725607 -83.7094236 -83.7553696 S1 Energy -83.3718665 -83.6909103 -83.7294617 S2 Energy -83.3196832 -83.6386296 -83.6745586 S3 Energy -83.3193799 -83.6383230 -83.6742271 S4 Energy -83.3190619 -83.6380353 -83.6741183 127 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsBrD CI listed in Cartesian Coordinates: Cs 0.9001884270 -1.6511087726 -0.5461170586 Cs 6.1812021176 -0.3204819309 -0.9538132757 Cs -1.4499659060 9.8145137165 7.8735786356 Cs 5.9916873370 0.3139165581 4.0120356893 Br 3.6669609522 -0.9685408400 1.6158243532 Br 2.9874386637 0.6027944339 -2.0655994583 Br 0.0883523678 0.9479550258 1.4951981196 Br 8.3160739271 0.0170912836 1.5966611624 Br -3.2583287077 19.5364456357 13.7453534524 Br 3.3815224324 2.5329048319 3.6809406010 Pb 2.6982901538 1.6747894476 1.1095053230 Figure A34: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsBrD CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. FOMO-CASCI Optimized Cluster Model CsBrD CI Figure A35: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the CsBrD region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. 128 Table A27: Absolute energies of the singlet electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsBrD CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -83.26619082 S1 Energy -83.26604766 Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized CsBrD CI listed in Cartesian Coordinates: Cs -0.7336907790 5.3147694955 -0.6386474082 Cs 4.5913947315 5.5268721480 -0.6319402505 Cs -0.2676922182 6.5655372536 5.0071232723 Cs 7.0970864185 -4.1843420489 6.3351420547 Br 9.6483028267 -12.2319721373 7.2963702662 Br 1.8993111274 5.7977051805 -2.8841638943 Br -1.3528831017 4.0340708132 2.7086523364 Br 4.4547653552 2.8657573165 1.8666717839 Br 1.9816917515 6.0702533681 2.0329266647 Br 2.2043045529 4.0043927686 5.5569278094 Pb 1.8398524271 3.3001319276 2.9005126002 Figure A36: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(2/5)/SBKJC VDZ ECP CsBrD CI calculated at the same level of theory as the optimization. SA-CASSCF Optimized Cluster Model CsBrD S1 Minimum Figure A37: Local S1 minimum in the cluster model CsBrD PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. 129 Table A28: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD PES region S1 minimum calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.2836087 S3 Energy -83.2374719 S1 Energy -83.2725812 S2 Energy -83.2399448 S4 Energy -83.2349863 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized CsBrD region S1 minimum listed in Cartesian Coordinates: Cs -0.5506338756 5.9697092444 -0.4139602624 Cs 5.0308854811 4.4173925455 -0.4562802283 Cs -0.0287573535 6.0293607601 5.4517697443 Cs 5.9165927583 -1.1827028424 4.2830708524 Br 9.2785917869 -11.0008161344 7.2316417392 Br 2.1983684403 4.9219565024 -2.5053353581 Br -1.2601573549 4.1492137609 2.5965251145 Br 4.5446309101 2.0149787549 2.3746256161 Br 2.5164248207 5.7257463582 2.4207438653 Br 1.9155282095 3.1060454672 5.6733831147 Pb 1.8309497015 2.9122579118 2.8936191220 MS-CASPT2-c19 Optimized Cluster Model CsBrD S1 Minimum Figure A38: Local S1 minimum in the cluster model CsBrD PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. Table A29: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized CsBrD PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.3998858 -83.7387985 -83.7870701 S2 Energy -83.3429888 -83.6648132 -83.7005838 S1 Energy -83.3769085 -83.6984527 -83.7346164 S3 Energy -83.3396961 -83.6611672 -83.6967291 S4 Energy -83.3396961 -83.6598055 -83.6956673 130 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized CsBrD region S1 minimum listed in Cartesian Coordinates: Cs -0.2351765501 5.7019411545 -0.3105143228 Cs 4.7278368661 4.3419015868 -0.2821946851 Cs 0.0870935111 5.8906187438 5.2476346793 Cs 5.8763039933 -1.0527939259 4.2155294872 Br 9.2786190339 -11.0010654438 7.2317362191 Br 2.2190014242 5.0030357220 -2.3516688450 Br -1.2694244308 4.2288705275 2.6296326973 Br 4.4523650399 2.0889064277 2.3292201427 Br 2.5214836976 5.7065944639 2.4249121895 Br 1.9470959647 3.1422118463 5.5739855229 Pb 1.7872250298 3.0129211946 2.8415299706 SA-CASSCF Optimized Cluster Model Bal CI Figure A39: Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized MECI in the Bal region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. Table A30: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.2733404 S1 Energy -83.2725247 S2 Energy -83.2387540 S3 Energy -83.2325133 S4 Energy -83.2211943 131 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI listed in Cartesian Coordinates: Cs -0.4617386917 2.9199453262 -2.9935844646 Cs 4.7318823408 3.2987308911 -0.2571350572 Cs 0.6561939748 5.2035622430 5.7824280965 Cs 6.1687129274 3.7196251222 5.2582188645 Br 4.5244562904 -2.6587060206 4.3639827117 Br 1.2525052968 1.5891976533 0.0762411700 Br -0.9514632505 2.5403112151 3.6008990213 Br 7.3028504976 2.7969371043 2.0681128949 Br 2.8539359763 3.6129411968 3.1182381564 Br 3.6374542989 4.7034212902 7.5151660287 Pb 1.3457447334 1.2951766124 2.8087888815 Figure A40: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the same level of theory as the optimization. Figure A41: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the SA-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs. 132 Table A31: Absolute energies of the singlet electronic states calculated at the SA5- CASSCF(2/5)/cc-PVDZ-PP level of theory with fully relativistic Stuttgart/Cologne ECPs at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI. n 0 1 2 3 4 Sn Energy -2764.72490975 -2764.72406804 -2764.69245106 -2764.68668082 -2764.67824157 Table A32: Eigenvalues of the spin-orbit matrix generated using the first five singlet states and first four triplet states calculated at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal CI in the molecular orbitals calculated at the SA5-CASSCF(2/5)/cc-PVDZ-PP with fully relativistic Stuttgart/Cologne ECPs level of theory. N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Eigenvalue (Hartree) -2764.72491017 -2764.72411782 -2764.72411782 -2764.72411781 -2764.72406942 -2764.69247675 -2764.69247675 -2764.69247675 -2764.69247302 -2764.68667105 -2764.68667105 -2764.68667100 -2764.68667025 -2764.67823245 -2764.67823245 -2764.67823245 -2764.67822844 MS-CASPT2-c19 Optimized Bal CI Figure A42: Geometry of the FOMO-CASCI(2/5)/SBKJC VDZ ECP optimized MECI in the Bal region of the cluster model PES. The optimization of the MECI was performed using the CIOpt software package. 133 Table A33: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized Bal CI. Energies are calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.3813495 -83.7106403 -83.7569656 S1 Energy -83.3805418 -83.6915325 -83.7312786 S2 Energy -83.3451308 -83.6548811 -83.6915805 S3 Energy -83.3388573 -83.6495516 -83.6862347 S4 Energy -83.3307575 -83.6416177 -83.6782958 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized Bal CI listed in Cartesian Coordinates: Cs -0.0855034504 3.3547533713 -2.9367396561 Cs 4.4852779451 3.0322398796 -0.0994750871 Cs 0.7618721619 5.0806143820 5.7663072385 Cs 5.8071177396 3.7649778008 5.0958119745 Br 5.3247137332 -2.7440612545 4.3835191476 Br 1.1710861074 1.4719847466 0.1671205576 Br -0.9562554937 2.6686293557 3.6406560074 Br 6.8587759887 2.1948392331 2.1424296814 Br 2.7372395382 3.9333791662 2.9168173718 Br 3.5747222721 4.7454222418 7.3583546645 Pb 1.3814874080 1.5183637091 2.9065651367 Figure A43: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized Bal CI calculated at the SA5-CASSCF(2/5)/SBKJC Polarized(p, 2d) level of theory. 134 SA-CASSCF Optimized Cluster Model Bal S1 Minimum Figure A44: Local S1 minimum in the cluster model Bal PES region optimized at the SA5- CASSCF(2/5)/SBKJC VDZ ECP level of theory using the CIOpt software package. Table A34: Absolute energies of the singlet electronic states at the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal PES region S1 minimum calculated at the same level of theory as the optimization. Level of Theory SA-CASSCF S0 Energy -83.2918032 S1 Energy -83.2886197 S2 Energy -83.2550008 S3 Energy -83.2473327 S4 Energy -83.2387471 Geometry of the SA5-CASSCF(2/5)/SBKJC VDZ ECP optimized Bal region S1 minimum listed in Cartesian Coordinates: Cs 0.8894638639 2.5066729217 -3.7524813422 Cs 3.7539231233 0.9766943132 0.9699087215 Cs 2.3340570101 4.0034352189 5.4247922753 Cs 7.5514166279 3.9187672029 5.7653438634 Br 7.0381808883 1.1386357820 3.0664678812 Br 0.2056961060 1.7981598959 -0.0983494273 Br -0.7454936103 4.1932730820 3.1950028000 Br 5.2490207029 3.6028999783 2.9627071770 Br 1.0835332773 0.6943568599 3.6564972691 Br 4.8551563465 4.5454877940 7.8209852580 Pb -1.1844199415 1.6427595867 2.3304818278 135 MS-CASPT2-c19 Optimized Bal S1 Minimum Figure A45: Local S1 minimum in the cluster model Bal PES region optimized at the MS5- CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) level of theory using the CIOpt software package. Table A35: Absolute energies of the singlet electronic states calculated at the MS5-CASPT2(2/5)- c19/SBKJC Polarized(p, 2d) optimized Bal PES region S1 minimum. Energies are calculated at the MS5-CASPT2(2/5)-c19/SBKJC Polarized (p, 2d), MS5-CASPT2(2/5)-c0/SBKJC Polarized (p, 2d), and Davidson corrected MRCI using the same number of electronic states and active space levels of theory. Level of Theory MS-CASPT2-c19 MS-CASPT2-c0 MRCI S0 Energy -83.4252856 -83.7568589 -83.8102052 S1 Energy -83.3976913 -83.7263629 -83.7612020 S2 Energy -83.3642980 -83.6922617 -83.7265445 S3 Energy -83.354953 -83.6824482 -83.7166224 S4 Energy -83.3510712 -83.6795907 -83.7138062 Geometry of the MS5-CASPT2(2/5)-c19/SBKJC Polarized(p, 2d) optimized Bal region S1 minimum listed in Cartesian Coordinates: Cs 0.8823939733 2.4876302615 -3.7072544003 Cs 3.5661712439 1.1422333113 1.1006846316 Cs 2.3665466649 3.9176863079 5.3169210023 Cs 7.3589515309 3.9124720292 5.7206598074 Br 7.1344461905 1.0364034366 3.0787940527 Br 0.2882749193 1.8281752674 -0.0928206712 Br -0.6517503987 4.2109135256 3.2148360369 Br 5.2121101478 3.6533003179 2.9689220483 Br 1.1743446412 0.6708766419 3.7054457963 Br 4.7775520498 4.4938242335 7.6770550038 Pb -1.0485066015 1.6676272614 2.3581125303 136 PBE0 Optimized Nanoparticle FC Point Figure A46: Ground state minimum energy geometry of the NP model optimized at the PBE0/SBKJC VDZ ECP level of theory. Table A36: Absolute energies of the singlet electronic states at the PBE0 optimized ground state minimum energy geometry calculated at the TD-PBE0/SBKJC VDZ ECP and FOMO- CASCI(6/11)/SBKJC VDZ ECP levels of theory. Level of Theory S1 Energy -285.63822620 TD-PBE0 FOMO-CASCI -280.08586514 S0 Energy -285.78527761 -280.25924118 137 PBE0/SBKJC VDZ ECP optimized ground state minimum energy geometry listed in Caresian Coordinates: Pb 2.0094997031 2.3801990255 2.4130958099 Br 5.0345338491 2.4805497638 2.5925691555 Br 1.6021438509 5.2272904531 2.6985187061 Br 2.2708430645 2.3806202869 5.5485445531 Pb 2.6445984311 2.7984296677 8.4346449324 Br 2.0208423787 5.7537435306 7.8787122369 Br 2.4983525255 -0.6971825003 2.8843825143 Pb 2.7502685703 -3.6279026880 2.8815990728 Br 5.8757293613 -3.3008719895 2.4932315403 Br 1.9051956637 1.6931005119 -0.3936524223 Br -0.8510382042 1.5812372236 3.0697518311 Br 2.7903325792 2.6167178839 11.4922604402 Br 5.5129830052 2.9637497301 7.9591725642 Br 2.8963438135 -0.1322969172 8.4316446745 Pb 3.3852909046 -3.2095539402 8.9031743774 Br 0.3602890591 -3.3101167189 8.7239822753 Br -0.4808084237 2.4717276723 8.8231986110 Br 3.3745379466 -6.5829861396 3.4381138278 Br 3.1236973878 -3.2100271581 5.7676485270 Br -0.1181006546 -3.7933850868 3.3564530475 Br 2.6049826656 -3.4460859323 -0.1761560070 Br 3.7928121179 -6.0566580863 8.6178068535 Br 3.4900578250 -2.5219540338 11.7097736875 Br 6.2457437793 -2.4103477132 8.2462319164 Cs 4.9021034432 -0.3473169709 0.2787726197 Cs 5.7200639523 0.3066430032 10.7276294839 Cs 0.9019653970 -6.0481903926 6.2962699241 Cs -0.3247932962 -1.1355961612 0.5879261004 Cs 0.4558438701 5.2665657353 11.0309379468 Cs 4.9394898262 -6.0958973882 0.2858165855 Cs 6.0824199247 -5.1852973909 5.8378686320 Cs 4.4928023019 5.2187396206 5.0198475684 Cs 0.4928889885 -0.4819714524 11.0374471926 Cs 0.2789999705 -0.7351507083 5.8501344525 Cs 5.1155432271 -0.0942496844 5.4659379606 Cs -0.6873837809 4.3561929395 5.4785778462 138 Figure A47: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the PBE0 ground state minimum energy geometry calculated at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory. FOMO-CASCI Optimized Nanoparticle FC S1 Minimum Figure A48: Local S1 minimum in the NP model FC PES region optimized at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory using the CIOpt software package. Table A37: Absolute energy of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP FC region S1 minimum, calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -280.22539119 S1 Energy -280.13558402 139 Geometry of the optimized FOMO-CASCI(6/11)/SBKJC VDZ ECP FC region S1 minimum listed in Cartesian Coordinates: Pb 1.9757352796 2.4899123001 2.3130500896 Br 5.1836062671 2.5009800059 2.5802376240 Br 1.5753134379 5.4175028679 2.7936713735 Br 2.2767828810 2.4019860681 5.5956525875 Pb 2.6088406296 2.7796247160 8.4275597571 Br 2.1543459308 5.5861407088 7.9674719930 Br 2.4993085375 -0.7368314084 2.8452412698 Pb 2.7194086336 -3.6602333484 2.8545455213 Br 5.9796097672 -3.3199745013 2.4664977996 Br 1.8737701678 1.6484149846 -0.5240461840 Br -0.9399443834 1.6129396451 3.0695795000 Br 2.8619402640 2.7003418781 11.2854961395 Br 5.4358747602 2.9923082809 8.0034760457 Br 2.8821529428 -0.0801218302 8.4763645778 Pb 3.3993500610 -3.2946003633 9.0119800618 Br 0.2496534911 -3.3015247046 8.7619769490 Br -0.2547239141 2.3635809096 8.6709804927 Br 3.4536979385 -6.7647992621 3.4588488472 Br 3.1173731385 -3.2302037116 5.7449530163 Br -0.1768014120 -3.7930170901 3.3617274598 Br 2.6307984385 -3.4604110875 -0.3732293551 Br 3.8115067120 -6.1860105096 8.5974246431 Br 3.4949315238 -2.3980152269 11.8953235488 Br 6.3815515785 -2.3917572344 8.2321274486 Cs 4.9382100096 -0.3450339721 0.1421857222 Cs 5.8449950730 0.2283019941 10.8478047831 Cs 0.8369819145 -6.1586661225 6.3123148715 Cs -0.3978895993 -1.1403771689 0.4917811568 Cs 0.1389037079 5.5983196140 11.3939637263 Cs 4.9509190871 -6.1392031034 0.2485542414 Cs 6.1456485617 -5.2826251984 5.8514764961 Cs 4.5843333934 5.3450810877 4.9159089544 Cs 0.3649911124 -0.5321617609 11.2887415735 Cs 0.0766715962 -0.8163109298 5.8647371302 Cs 5.2948223196 -0.0942017259 5.3852633630 Cs -0.8635948334 4.5331131908 5.4322258150 140 FOMO-CASCI Optimized Nanoparticle BrD CI Figure A49: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the BrD region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. Table A38: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -280.15706720 S1 (Energy) -280.15698939 141 Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD CI listed in Cartesian Coordinates: Pb 1.9918750195 2.4617263737 2.3041438883 Br -0.9583607093 1.6075568775 3.0314447316 Br 5.1167729823 2.4941191048 2.5294250795 Br 1.5663825354 5.4009717942 2.7250638586 Br 2.2838282779 2.4528223068 5.5709331949 Pb 2.7110457587 2.8756217455 8.4605185323 Br 2.9105986140 -0.0976170668 8.5084958991 Br 1.9627313121 1.6919263515 -0.6111788983 Br 1.9667600089 5.9334944763 7.9222010284 Br 2.7781608383 2.6917127306 11.6727465814 Br 5.5783448296 3.0192130945 7.9946040838 Pb 3.4278744029 -3.3080853950 9.0268732637 Br 6.3682571655 -2.4376136743 8.2273905685 Br 0.2525865581 -3.3660217501 8.7514995344 Br 2.9686666968 -3.3178796792 5.7046849786 Pb 2.6376912778 -3.7157521942 2.8485517771 Br 2.5168318210 -3.5094896295 -0.3466646630 Br 2.3753195609 -0.7658017566 2.8163965046 Br 3.3326723457 -6.7813419551 3.4480481800 Br -6.3862698248 -4.1803153723 4.4335952770 Br 3.9174935881 -6.2916957419 8.6250414843 Br 3.4885872627 -2.5154084652 11.8550855628 Br -0.5786317699 2.4877053393 8.8500574706 Br 6.0254670864 -3.3288232859 2.4402679788 Cs 0.3505558062 5.3209937477 11.1131043504 Cs 4.9509615057 -6.1346169037 0.2648224271 Cs -0.7850319610 4.4883642429 5.4224515193 Cs 6.1043131927 -5.2753075534 5.8198436706 Cs 0.4497320544 -0.5165182092 11.2047196546 Cs 4.9222199694 -0.3764272125 0.1470774296 Cs 0.9481992299 -6.3907369397 6.5636901187 Cs 4.5330934331 5.3524263270 4.9683290442 Cs -0.4061601083 -0.9042081679 0.2406295075 Cs 5.7733244040 0.2756269714 10.8379152125 Cs 5.2418035064 -0.0930332647 5.3967399936 Cs 0.0899313382 -0.5779212673 6.0790361790 142 Figure A50: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD CI calculated at the same level of theory as the optimization. FOMO-CASCI Optimized Nanoparticle BrD S1 Minimum Figure A51: Local S1 minimum in the NP model BrD PES region optimized at the FOMO- CASCI(6/11)/SBKJC VDZ ECP level of theory using the CIOpt software package. Table A39: Absolute energy of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP BrD region S1 minimum, calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -280.17596539 S1 Energy -280.16795508 143 Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized BrD region S1 minimum listed in Cartesian Coordinates: Pb 2.0110022996 2.5180987798 2.2842989204 Br -1.0679201927 1.8873463134 2.8549560016 Br 5.1228883900 2.3259291357 2.6322165251 Br 1.7701556925 5.4376051617 2.7166784494 Br 2.1327128675 2.4000886345 5.6064788338 Pb 2.5546951812 2.8284251975 8.4832699184 Br 2.7843408843 -0.1162597629 8.5181989980 Br 2.1387788393 1.7768231364 -0.6479013156 Br 1.9079670155 5.9109682344 7.8985185246 Br 2.6968543429 2.6457980593 11.6919645310 Br 5.4446255839 2.9474817231 7.9640091948 Pb 3.4483275092 -3.2910264143 9.0380367401 Br 6.2656254544 -2.3722462015 8.1001261512 Br 0.2180394029 -3.4521454679 9.0455840226 Br 2.6926229750 -3.3218836816 5.7627351313 Pb 2.2879288539 -3.7600534153 2.8778397058 Br 2.1604465034 -3.5110755629 -0.3031193893 Br 2.0430666685 -0.7959594097 2.8568234211 Br 2.9524068603 -6.7917316467 3.5321938631 Br -3.8860767118 -4.0897052902 3.9705321218 Br 4.1201566836 -6.2329994691 8.6508324301 Br 3.7024872025 -2.4611708462 11.8600969008 Br -0.6480798724 2.5041082256 8.9017966951 Br 6.1962009369 -3.4883212043 2.2706327135 Cs 0.3528966431 5.3366612133 11.1107351624 Cs 4.6292078924 -6.1335676135 0.3325381837 Cs -0.7247586818 4.6402934870 5.3380775498 Cs 5.9458986619 -5.2032085832 5.6357109386 Cs 0.5527901724 -0.5463456569 11.3863938572 Cs 4.8044539300 -0.5796130333 0.2848884773 Cs 1.0195002661 -6.4817627403 6.8573915432 Cs 4.6093274433 5.2279657234 5.0438657292 Cs -0.4261460912 -0.6554268226 0.0572444282 Cs 5.7971437688 0.3628514861 10.7379814727 Cs 4.9709860983 -0.1830049624 5.3710910395 Cs -0.1129254314 -0.6132705471 6.1248674135 144 FOMO-CASCI Optimized Nanoparticle CsD CI Figure A52: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the CsD region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. Table A40: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -280.09048753 S1 Energy -280.09019333 145 Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized CsD CI listed in Cartesian Coordinates: Pb 2.1152151537 2.4358266257 2.3665103214 Br 1.9764772945 1.6257598607 -0.5721481139 Br -1.1608056245 1.5064453288 3.1541986738 Br 4.9932076886 2.7629296622 2.2568939908 Br 1.6413363531 5.4268675819 2.7569892971 Br 2.3194850275 2.4604495732 5.5538320381 Pb 2.7077522168 2.8232269466 8.3460125833 Br 5.5472813668 3.2348791524 8.1219726189 Br 2.8947067408 -0.2559558432 8.5225259207 Br 2.0257991042 5.7949712112 7.8402148084 Br 2.7549411528 2.4983891557 11.8823720974 Pb 3.3976037534 -3.1750816786 8.8777963444 Br 3.4738906012 -2.6554424828 11.6514170052 Br 6.1645595268 -2.7111636888 8.5691752602 Br 0.4739995829 -3.3936504070 8.8457934627 Br 3.0681794494 -3.2576824913 5.8823783645 Pb 2.7908129748 -3.5618949025 2.8746722524 Br -0.1224081694 -3.7527840717 3.3718540882 Br 3.7642229521 -5.9561427384 8.7219165124 Br 2.5745474573 -3.4592585812 -0.2382009355 Br 2.4688652942 -0.7080227982 2.7917081499 Br 3.5492126985 -6.9267647891 3.5788715552 Br -0.6974266087 2.4873206360 8.8751516245 Br 6.0314146639 -3.5203633957 2.2561183004 Cs 0.5139872247 5.1844141748 11.0383053393 Cs 4.7921326846 -6.2249304328 0.3003896692 Cs -0.7291790864 4.3921974649 5.3748509847 Cs 6.1295641702 -5.1050730348 5.6345011768 Cs 0.4355552127 -0.3251003181 11.1574950860 Cs 5.0009076866 -0.3402229174 0.2990380156 Cs 0.8858964419 -6.1705464225 6.1451086832 Cs 4.5839826798 5.2395985042 4.9710117428 Cs -0.3541327589 -0.9987226375 0.5414681940 Cs 5.7346112625 0.4770282326 10.6902996747 Cs 0.2624825371 -0.6830887098 5.7898331136 Cs 16.2612372913 1.4247932256 4.8644761137 146 Figure A53: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized CsD CI calculated at the same level of theory as the optimization. FOMO-CASCI Optimized NP DBrS CI Figure A54: Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized MECI in the DBrS region of the NP model PES. The optimization of the MECI was performed using the CIOpt software package. Table A41: Absolute energies of the singlet electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. Level of Theory FOMO-CASCI S0 Energy -280.16734851 S1 Energy -280.16685486 147 Geometry of the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized DBrS CI listed in Cartesian Coordinates: Pb 1.9458174450 2.4582899903 2.2181926129 Br -0.9465973363 1.7244424394 2.9765909107 Br 5.2056953680 2.5066378557 2.5161792068 Br 1.5648633602 5.4438844829 2.7022495687 Br 2.2477624249 2.4067571994 5.6523027083 Pb 2.6530779182 2.8500385504 8.4319725636 Br 2.9101571489 -0.2253989202 8.4961445646 Br 1.8653405259 1.5698072233 -0.5442641962 Br 1.9641152618 5.9016669957 7.8613079185 Br 2.7441649078 2.5735182465 11.8601518479 Br 5.5719807825 2.9318756774 7.9593422612 Pb 3.3843534582 -3.3122624061 8.9036618572 Br 6.1880631848 -2.5491657544 8.3628640152 Br 0.5126045417 -3.4378989663 8.9593773168 Br 3.0281607132 -3.2864782872 5.9414357295 Pb 2.7669626567 -3.7545772249 2.7642821831 Br 2.6603081702 -3.5105397839 -0.3611972149 Br 2.4976135218 -0.7678826327 2.8754263710 Br 3.5854292595 -6.8360457486 3.4906507884 Br -1.4771113899 -3.5746587879 3.8224878110 Br 3.8214367934 -6.0261743048 8.7368875553 Br 3.5443961734 -2.6121286793 11.6954121992 Br -0.6418720513 2.5419582984 8.9084493052 Br 7.0431217765 -3.3209989095 2.2872693984 Cs 0.4899034598 5.2944931484 11.1194102975 Cs 5.2932577765 -5.9383187635 0.3638544993 Cs -0.7500585015 4.5215450888 5.4241643558 Cs 6.2240096327 -5.1675384850 5.5397535735 Cs 0.4781404960 -0.3446820166 11.1949429743 Cs 5.0146749176 -0.4761099680 0.2275024193 Cs 0.5195028202 -5.9854326857 5.7865098433 Cs 4.5229846390 5.3130668914 4.9189881065 Cs -0.5144484658 -1.3953970809 0.7777459039 Cs 5.7086243031 0.4274419726 10.8528861076 Cs 5.2419338558 -0.0154386296 5.3344908664 Cs -0.0059465312 -0.7850900179 5.7499927761 148 Figure A55: Fractionally occupied state averaged natural orbitals averaged over the S0 and S1 electronic states at the FOMO-CASCI(6/11)/SBKJC VDZ ECP optimized DBrS CI calculated at the same level of theory as the optimization. CsD Potential Energy Surface Scans and Dissociated Cs Charges Figure A56: a) Cluster model SA-CASSCF PES scan as the dissociated Cs atom at the CsD CI is displaced along the vector between it and the Pb center while all other degrees of freedom are frozen. Energies are shifted such that the S0 energy at the FC point is zero. b) Cluster model Mulliken charges of the dissociated Cs atom as it is displaced along the dissociated Cs–Pb vector around the SA-CASSCF optimized CsD CI with all other degrees of freedom frozen. a) b) 149 Figure A57: a) NP model FOMO-CASCI PES scan as the dissociated Cs atom at the CsD CI is displaced along the vector between itself and the average position of the four Pb atoms while all other degrees of freedom are frozen. Energies are shifted such that the S0 at the FC point is zero. b) NP model Mulliken charges of the dissociated Cs atom as it is displaced along the vector between it and the average position of the four Pb atoms around the FOMO-CASCI optimized CsD CI with all other degrees of freedom frozen. a) b) 150 BIBLIOGRAPHY 151 BIBLIOGRAPHY 1 B. M. Bode, and M. S. Gordon, Journal of Molecular Graphics and Modelling 16 (1998) 133. 2 A. Kojima et al., Journal of the American Chemical Society 131 (2009) 6050. 3 W. S. Yang et al., Science 348 (2015) 1234. 4 Q. Jiang et al., Advanced Materials 29 (2017) 1703852. 5 S. D. Stranks et al., Science 342 (2013) 341. 6 A. L. Abdelhady et al., The Journal of Physical Chemistry Letters 7 (2016) 295. 7 D. A. Egger et al., Advanced Materials 30 (2018) 1800691. 8 J. P. H. Rivett et al., Nature Communications 9 (2018) 3531. 9 L. Wang et al., Advanced Optical Materials 6 (2018) 1700975. 10 J. Kim et al., The Journal of Physical Chemistry Letters 5 (2014) 1312. 11 W.-J. Yin, T. Shi, and Y. Yan, Applied Physics Letters 104 (2014) 063903. 12 R. Long, and O. V. Prezhdo, ACS Nano 9 (2015) 11143. 13 D. R. Yarkony, Reviews of Modern Physics 68 (1996) 985. 14 J. Michl, 1939-, and V. Bonačić-Koutecký, (1990) 15 B. G. Levine, and T. J. Martínez, Annual Review of Physical Chemistry 58 (2007) 613. 16 F. Bernardi, M. Olivucci, and M. A. Robb, Chem. Soc. Rev. 25 (1996) 321. 17 W. Domcke, and D. R. Yarkony, Annual Review of Physical Chemistry 63 (2012) 325. 18 Y. Shu, and B. G. Levine, The Journal of Physical Chemistry C 120 (2016) 23246. 19 W.-T. Peng et al., Chemical Science 9 (2018) 681. 20 Y. Shu et al., The Journal of Physical Chemistry Letters 8 (2017) 4091. 21 Y. Shu, B. S. Fales, and B. G. Levine, Nano Letters 15 (2015) 6247. 22 Y. Shu et al., The Journal of Physical Chemistry C 119 (2015) 26683. 23 Y. Shu, and B. G. Levine, The Journal of Physical Chemistry C 118 (2014) 7669. 152 24 Y. Shu, and B. G. Levine, The Journal of Chemical Physics 139 (2013) 081102. 25 B. G. Levine et al., Annual Review of Physical Chemistry 70 (2019) 21. 26 B. G. Levine, W.-T. Peng, and M. P. Esch, Physical Chemistry Chemical Physics 21 (2019) 10870. 27 W.-T. Peng, and B. G. Levine, The Journal of Physical Chemistry C 123 (2019) 16588. 28 I. S. Ufimtsev, and T. J. Martinez, Journal of Chemical Theory and Computation 5 (2009) 2619. 29 A. V. Titov et al., Journal of Chemical Theory and Computation 9 (2013) 213. 30 E. G. Hohenstein et al., The Journal of Chemical Physics 143 (2015) 014111. 31 B. S. Fales, and B. G. Levine, Journal of Chemical Theory and Computation 11 (2015) 4708. 32 E. G. Hohenstein, The Journal of Chemical Physics 145 (2016) 174110. 33 B. S. Fales, E. G. Hohenstein, and B. G. Levine, Journal of Chemical Theory and Computation 13 (2017) 4162. 34 Y. Shu, and B. G. Levine, The Journal of Chemical Physics 139 (2013) 074102. 35 Y. Shu, E. G. Hohenstein, and B. G. Levine, The Journal of Chemical Physics 142 (2015) 024102. 36 W. H. Miller, and C. W. McCurdy, The Journal of Chemical Physics 69 (1978) 5163. 37 H. D. Meyera, and W. H. Miller, The Journal of Chemical Physics 70 (1979) 3214. 38 D. A. Micha, The Journal of Chemical Physics 78 (1983) 7138. 39 J. M. Cohen, and D. A. Micha, The Journal of Chemical Physics 97 (1992) 1038. 40 J. C. Tully, The Journal of Chemical Physics 93 (1990) 1061. 41 O. V. Prezhdo, and P. J. Rossky, The Journal of Chemical Physics 107 (1997) 5863. 42 B. J. Schwartz, and P. J. Rossky, The Journal of Chemical Physics 101 (1994) 6902. 43 E. R. Bittner, and P. J. Rossky, The Journal of Chemical Physics 103 (1995) 8130. 44 B. J. Schwartz et al., The Journal of Chemical Physics 104 (1996) 5942. 45 M. Thachuk, M. Y. Ivanov, and D. M. Wardlaw, The Journal of Chemical Physics 109 (1998) 5747. 153 46 J.-Y. Fang, and S. Hammes-Schiffer, The Journal of Chemical Physics 110 (1999) 11166. 47 J.-Y. Fang, and S. Hammes-Schiffer, The Journal of Physical Chemistry A 103 (1999) 9399. 48 M. D. Hack, and D. G. Truhlar, The Journal of Chemical Physics 114 (2001) 9305. 49 C. Zhu, A. W. Jasper, and D. G. Truhlar, The Journal of Chemical Physics 120 (2004) 5543. 50 M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, The Journal of Chemical Physics 123 (2005) 234106. 51 A. W. Jasper et al., Accounts of Chemical Research 39 (2006) 101. 52 G. Granucci, and M. Persico, The Journal of Chemical Physics 126 (2007) 134114. 53 G. Granucci, M. Persico, and A. Zoccante, The Journal of Chemical Physics 133 (2010) 134111. 54 J. E. Subotnik, and N. Shenvi, The Journal of Chemical Physics 134 (2011) 024105. 55 H. M. Jaeger, S. Fischer, and O. V. Prezhdo, The Journal of Chemical Physics 137 (2012) 22A545. 56 B. R. Landry, and J. E. Subotnik, The Journal of Chemical Physics 137 (2012) 22A513. 57 T. Nelson et al., The Journal of Chemical Physics 138 (2013) 224111. 58 J. E. Subotnik, W. Ouyang, and B. R. Landry, The Journal of Chemical Physics 139 (2013) 214107. 59 J. E. Subotnik et al., Annual Review of Physical Chemistry 67 (2016) 387. 60 L. Wang, A. Akimov, and O. V. Prezhdo, The Journal of Physical Chemistry Letters 7 (2016) 2100. 61 D. V. Makhov et al., The Journal of Chemical Physics 141 (2014) 054110. 62 D. A. Fedorov, and B. G. Levine, The Journal of Physical Chemistry Letters 10 (2019) 4542. 63 S. D. Stranks et al., Science 342 (2013) 341. 64 A. L. Abdelhady et al., J. Phys. Chem. Lett. 7 (2016) 295. 65 Q. Jiang et al., Advanced Materials 29, 1703852 (2017) 7. 66 L. Protesescu et al., Nano Letters 15 (2015) 3692. 67 G. Niu, X. Guo, and L. Wang, Journal of Materials Chemistry A 3 (2015) 8970. 154 68 S. Dastidar et al., Nano Letters 16 (2016) 3563. 69 L. N. Quan et al., Journal of the American Chemical Society 138 (2016) 2649. 70 J. Gebhardt, Y. Kim, and A. M. Rappe, The Journal of Physical Chemistry C 121 (2017) 6569. 71 L. Liu et al., The Journal of Physical Chemistry Letters 9 (2018) 1164. 72 L. Debbichi et al., Advanced Materials 30 (2018) 1707001. 73 E. D. A. et al., Advanced Materials 30 (2018) 1800691. 74 J. P. H. Rivett et al., Nat. Commun. 9, 3531 (2018) 8. 75 L. L. Wang et al., Adv. Opt. Mater. 6, 1700975 (2018) 8. 76 J. Liu, and O. V. Prezhdo, The Journal of Physical Chemistry Letters 6 (2015) 4463. 77 T. S. Nguyen, and J. Parkhill, J. Phys. Chem. A 120 (2016) 6880. 78 B. G. Levine et al., Annual Review of Physical Chemistry 70 (2019) review in advance. 79 D. F. Yuan et al., Science 362 (2018) 1289. 80 J. B.-K. Michl, V., Electronic Aspects of Organic Photochemistry (Wiley, New York, 1990), 81 Y. Shu et al., J. Phys. Chem. Lett. 8 (2017) 4091. 82 Y. Shu et al., J. Phys. Chem. C 119 (2015) 26683. 83 Y. Shu, and B. G. Levine, J. Phys. Chem. C 118 (2014) 7669. 84 Y. Shu, and B. G. Levine, J. Chem. Phys. 139, 081102 (2013) 081102. 85 M. Ben-Nun, J. Quenneville, and T. J. Martínez, The Journal of Physical Chemistry A 104 (2000) 5161. 86 B. G. Levine et al., Chem. Phys. 347 (2008) 3. 87 B. F. E. Curchod, and T. J. Martinez, Chem. Rev. 118 (2018) 3305. 88 B. O. Roos, P. R. Taylor, and P. E. M. Si≐gbahn, Chem. Phys. 48 (1980) 157. 89 H. J. Werner, and P. J. Knowles, The Journal of Chemical Physics 82 (1985) 5053. 90 P. J. Knowles, and H.-J. Werner, Chemical Physics Letters 115 (1985) 259. 91 W. J. Stevens et al., Canadian Journal of Chemistry 70 (1992) 612. 92 H. Sekino, and J. Bartlett Rodney, International Journal of Quantum Chemistry 26 (1984) 255. 155 93 C. Hampel, K. A. Peterson, and H.-J. Werner, Chemical Physics Letters 190 (1992) 1. 94 T. Korona, and H.-J. Werner, The Journal of Chemical Physics 118 (2003) 3006. 95 A. D. Becke, The Journal of Chemical Physics 98 (1993) 5648. 96 C. Lee, W. Yang, and R. G. Parr, Physical Review B 37 (1988) 785. 97 S. H. Vosko, L. Wilk, and M. Nusair, Canadian Journal of Physics 58 (1980) 1200. 98 P. Labello Nicholas, M. Ferreira Antonio, and A. Kurtz Henry, Journal of Computational Chemistry 26 (2005) 1464. 99 P. Labello Nicholas, M. Ferreira Antonio, and A. Kurtz Henry, International Journal of Quantum Chemistry 106 (2006) 3140. 100 F. Eckert, P. Pulay, and H. J. Werner, Journal of Computational Chemistry 18 (1998) 1473. 101 R. Lindh, Theoretica chimica acta 85 (1993) 423. 102 H. J. Werner et al., 103 G. A. Meek, and B. G. Levine, The Journal of Physical Chemistry Letters 5 (2014) 2351. 104 A. L. Thompson, C. Punwong, and T. J. Martínez, Chemical Physics 370 (2010) 70. 105 B. G. Levine, J. D. Coe, and T. J. Martínez, The Journal of Physical Chemistry B 112 (2008) 405. 106 J. Finley et al., Chemical Physics Letters 288 (1998) 299. 107 K. Andersson, P. Å. Malmqvist, and B. O. Roos, The Journal of Chemical Physics 96 (1992) 1218. 108 P. Celani, and H.-J. Werner, The Journal of Chemical Physics 112 (2000) 5546. 109 P. Celani, and H.-J. Werner, The Journal of Chemical Physics 119 (2003) 5044. 110 B. O. Roos, and K. Andersson, Chemical Physics Letters 245 (1995) 215. 111 D. G. Truhlar, and C. A. Mead, Physical Review A 68, 032501 (2003) 112 R. Langhoff Stephen, and R. Davidson Ernest, International Journal of Quantum Chemistry 8 (1974) 61. 113 P. J. Knowles, and H.-J. Werner, Chemical Physics Letters 145 (1988) 514. 114 H. J. Werner, and P. J. Knowles, The Journal of Chemical Physics 89 (1988) 5803. 156 115 K. A. Peterson et al., The Journal of Chemical Physics 119 (2003) 11113. 116 J. G. Hill, and K. A. Peterson, The Journal of Chemical Physics 147 (2017) 244106. 117 B. Metz, H. Stoll, and M. Dolg, The Journal of Chemical Physics 113 (2000) 2563. 118 I. S. Lim et al., The Journal of Chemical Physics 122 (2005) 104103. 119 A. Berning et al., Molecular Physics 98 (2000) 1823. 120 G. Granucci, M. Persico, and A. Toniolo, The Journal of Chemical Physics 114 (2001) 10608. 121 P. Slavíček, and T. J. Martínez, The Journal of Chemical Physics 132 (2010) 234102. 122 E. G. Hohenstein et al., J. Chem. Phys. 143, 014111 (2015) 014111. 123 E. G. Hohenstein, J. Chem. Phys. 145, 174110 (2016) 174110. 124 J. Endres et al., The Journal of Physical Chemistry Letters 7 (2016) 2722. 125 M. Kulbak et al., The Journal of Physical Chemistry Letters 7 (2016) 167. 126 P. Umari, E. Mosconi, and F. De Angelis, Scientific Reports 4 (2014) 4467. 127 J. Even et al., The Journal of Physical Chemistry Letters 4 (2013) 2999. 128 Y. Shu, and B. G. Levine, J. Chem. Phys. 139, 074102 (2013) 074102. 129 Y. Shu, E. G. Hohenstein, and B. G. Levine, J. Chem. Phys. 142, 024102 (2015) 024102. 130 B. G. Levine, and T. J. Martinez, Ann. Rev. Phys. Chem. 58 (2007) 613. 131 W. Domcke, and D. R. Yarkony, Ann. Rev. Phys. Chem. 63 (2012) 325. 132 M. S. Schuurman, and A. Stolow, Annu. Rev. Phys. Chem. 69 (2018) 427. 133 S. Kilina, D. Kilin, and S. Tretiak, Chem. Rev. 115 (2015) 5929. 134 L. J. Wang, R. Long, and O. V. Prezhdo, in Annual Review of Physical Chemistry, Vol 66, edited by M. A. Johnson, and T. J. Martinez (Annual Reviews, Palo Alto, 2015), pp. 549. 135 B. G. Levine et al., in Annual Review of Physical Chemistry, Vol 70, edited by M. A. Johnson, and T. J. Martinez (Annual Reviews, Palo Alto, 2019), pp. 21. 136 T. R. Nelson et al., Chem. Rev. 120 (2020) 2215. 137 H. D. Meyer, and W. H. Miller, J. Chem. Phys. 70 (1979) 3214. 138 J. M. Cohen, and D. A. Micha, J. Chem. Phys. 97 (1992) 1038. 157 139 J. C. Tully, and R. K. Preston, J. Chem. Phys. 55 (1971) 562. 140 S. J. Cotton, and W. H. Miller, J. Phys. Chem. A 117 (2013) 7190. 141 B. J. Schwartz, and P. J. Rossky, J. Chem. Phys. 101 (1994) 6902. 142 E. R. Bittner, and P. J. Rossky, J. Chem. Phys. 103 (1995) 8130. 143 O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 107 (1997) 5863. 144 M. Ben-Nun, and T. J. Martinez, Adv. Chem. Phys. 121 (2002) 439. 145 I. Burghardt, H. D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 111 (1999) 2927. 146 G. A. Worth, and I. Burghardt, Chem. Phys. Lett. 368 (2003) 502. 147 G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127 (2004) 307. 148 R. Lambert, and N. Makri, J. Chem. Phys. 137, 22a552 (2012) 149 T. Banerjee, and N. Makri, J. Phys. Chem. B 117 (2013) 13357. 150 N. Makri, Int. J. Quantum Chem. 115 (2015) 1209. 151 P. L. Walters, and N. Makri, J. Phys. Chem. Lett. 6 (2015) 4959. 152 J. Y. Fang, and S. Hammes-Schiffer, J. Chem. Phys. 110 (1999) 11166. 153 M. D. Hack, and D. G. Truhlar, J. Chem. Phys. 114 (2001) 9305. 154 C. Zhu et al., J. Chem. Phys. 121 (2004) 7658. 155 C. Y. Zhu, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 120 (2004) 5543. 156 A. W. Jasper et al., Acc. Chem. Res. 39 (2006) 101. 157 J. E. Subotnik, and N. Shenvi, J. Chem. Phys. 134, 024105 (2011) 158 J. E. Subotnik, W. J. Ouyang, and B. R. Landry, J. Chem. Phys. 139, 214107 (2013) 159 L. J. Wang, A. Akimov, and O. V. Prezhdo, J. Phys. Chem. Lett. 7 (2016) 2100. 160 D. A. Fedorov, and B. G. Levine, J. Phys. Chem. Lett. 10 (2019) 4542. 161 W. C. Swope et al., The Journal of Chemical Physics 76 (1982) 637. 162 S. Blanes, F. Casas, and A. Murua, The Journal of Chemical Physics 124 (2006) 234105. 163 S. Hammes-Schiffer, and J. C. Tully, J. Chem. Phys. 101 (1994) 4657. 158 164 A. W. Jasper, M. D. Hack, and D. G. Truhlar, J. Chem. Phys. 115 (2001) 1804. 165 A. Jain, and J. E. Subotnik, J. Chem. Phys. 143, 134107 (2015) 166 W. J. Dou, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. 142, 084110 (2015) 167 W. J. Dou, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. 142, 234106 (2015) 168 W. J. Ouyang, W. J. Dou, and J. E. Subotnik, J. Chem. Phys. 142, 084109 (2015) 169 P. L. Walters, and N. Makri, J. Chem. Phys. 144, 044108 (2016) 8. 170 C. F. Craig, W. R. Duncan, and O. V. Prezhdo, Phys. Rev. Lett. 95 (2005) 163001. 171 B. Mignolet, B. F. E. Curchod, and T. J. Martinez, J. Chem. Phys. 145 (2016) 191104. 172 D. A. Fedorov, and B. G. Levine, J. Phys. Chem. Lett. (2019) 4542. 173 O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 107 (1997) 825. 174 D. V. Makhov et al., J. Chem. Phys. 141 (2014) 054110. 175 D. V. Makhov, T. J. Martinez, and D. V. Shalashilin, Faraday Discuss. 194 (2016) 81. 176 V. M. Freixas et al., Phys. Chem. Chem. Phys. 20 (2018) 17762. 177 M. P. Esch, and B. G. Levine, The Journal of Chemical Physics 152 (2020) 234105. 178 J. E. Subotnik, and N. Shenvi, J. Chem. Phys. 134, 024105 (2011) 19. 179 J. E. Subotnik, W. J. Ouyang, and B. R. Landry, J. Chem. Phys. 139, 214107 (2013) 16. 180 B. Smith, and A. V. Akimov, J. Chem. Phys. 151, 124107 (2019) 16. 181 D. Neuhauser, J. Chem. Phys. 100 (1994) 5076. 182 M. R. Wall, and D. Neuhauser, J. Chem. Phys. 102 (1995) 8011. 183 A. Jain, and J. E. Subotnik, J. Chem. Phys. 143, 134107 (2015) 13. 184 T. Nelson et al., J. Chem. Phys. 138, 224111 (2013) 13. 185 J. Kang, and L. W. Wang, Phys. Rev. B 99, 224303 (2019) 8. 186 A. E. Sifain et al., J. Chem. Phys. 150, 194104 (2019) 8. 187 X. Li et al., J. Chem. Phys. 123 (2005) 084106. 188 C. M. Isborn, X. Li, and J. C. Tully, J. Chem. Phys. 126 (2007) 134307. 159 189 W. T. Peng, B. S. Fales, and B. G. Levine, J Chem Theory Comput 14 (2018) 4129. 160