ADDRESSING WAVE FUNCTION DISCONTINUITIES AT CONICAL INTERSECTIONS AND NOVEL CHARGE TRANSFER PROCESSES IN NANOMATERIALS  By  Garrett A. Meek            A DISSERTATION  Submitted to Michigan State University in partial fulfillment of the requirements for the degree of  ChemistryDoctor of Philosophy  2016     ABSTRACT ADDRESSING WAVE FUNCTION DISCONTINUITIES AT CONICAL  INTERSECTIONS AND NOVEL CHARGE TRANSFER PROCESSES IN NANOMATERIALS  By  Garrett A. Meek  The following dissertation describes work completed in fulfillment of the requirements for the degree, doctor of philosophy, in Chemistry from Michigan State University.  This work addresses complications that arise when the adiabatic representation of electronic states is used during nonadiabatic molecular dynamics.  Specifically, this work provides dynamical solutions for discontinuities in the electronic wave functions at conical intersections.  Additionally, static quantum chemical approaches are used to investigate a novel charge transport mechanism in graphitic carbon nitrides, a class of metal-free materials capable of photocatalytic water splitting.   iii                          William Shakespeare, Henry VI    iv  ACKNOWLEDGEMENTS   First and foremost, I want to thank my adviser, Professor Benjamin G. Levine.  sincerity and approachability were what initially drew me to his research group.  I have since learned that these are some of the attributes that make him such a great scientist, and I endeavor to develop these attributes for myself.  Ben taught me how to think clearly and communicate effectively, through both his own example and dedicated advising.  These lessons have enriched my personal life as much as the way that I conduct scientific research.   Thanks also to fellow graduate students Yinan Shu and Scott Fales, the other members of whose ideas and contributions were helpful at many stages of this work.  Thank you Professor Piotr Piecuch, for your thoughtful advising and dedicated efforts as my second reader.  Finally, I want to thank Danielle Lepar.  I drew upon many people for support while finishing this, but Danielle most of all.       v  TABLE OF CONTENTS  LIST OF TABLES viii LIST OF FIGURES ix CHAPTER 1: INTRODUCTION 1 1.1 Motivation 1 1.2 Description of Chapters 4 REFERENCES 6  CHAPTER 2: EVALUATION OF THE TIME-DERIVATIVE COUPLING FOR ACCURATE ELECTRONIC STATE TRANSITION PROBABILITIES FROM NUMERICAL SIMULATIONS 8 2.1 Introduction 8 2.2 The norm preserving interpolation (NPI) method 11 2.3 Methods 13 2.4 Results 14 2.5 Conclusions 19 APPENDICES 20 APPENDIX A: Integrated expressions for the time-derivative coupling 21 APPENDIX B: Numerical Integration Details 23 APPENDIX C: Three state Landau-Zener model results 25 REFERENCES 29  CHAPTER 3: ACCURATE AND EFFICIENT EVALUATION OF TRANSITION PROBABILITIES AT UNAVOIDED CROSSINGS IN AB INITIO MULTIPLE SPAWNING 34 3.1 Introduction 34 3.2 Incorporation of NPI into AIMS 38 3.3 Methods 42 3.4 AIMS-NPI near a conical intersection 44 3.5 AIMS-NPI at a trivial unavoided crossing 46 3.6 Discussion 49 3.7 Time averaging of the TDC 53 3.8 Conclusions 56 APPENDIX 57 REFERENCES 62  CHAPTER 4: WAVE FUNCTION CONTINUITY AND THE DIAGONAL BORN-OPPENHEIMER CORRECTION AT CONICAL INTERSECTIONS 69 4.1 Introduction 69 4.2 Proof of non-integrability of the DBOC at a CI in the adiabatic representation 71 4.3 Discussion 73 4.4 Conclusions 78 vi  APPENDICES 79 APPENDIX A: Proof that a Continuous Molecular Wave Function can be Decomposed into Continuous Nuclear and Electronic Factors 80 APPENDIX B: Illustrative Example of a Continuous Molecular Wave Function which can be Represented as a Sum of Discontinuous Born-Oppenheimer Vibronic States 81 REFERENCES 85  CHAPTER 5: THE BEST OF BOTH REPSDIABATIZED GAUSSIANS ON ADIABATIC SURFACES 91 5.1 Introduction 91 5.2 Electronic representations for nonadiabatic dynamics 92 5.3 Diabatized Gaussians on adiabatic surfaces (DGAS) 96 5.3.1 Review of the AIMS method 96 5.3.2 Generalization of AIMS 100 5.3.3 The Diabatized Gaussians on Adiabatic Surfaces Wave Function 101 5.3.4 Calculation of Derivative Couplings in DGAS 104 5.3.5 Implementation of DBOC in Adiabatic AIMS 108 5.3.6 Computational methods and systems of interest 111 5.4 The effect of enforcing continuity at a conical intersection 114 5.5 Analysis of errors in DBOC gradients 117 5.6 The effect of enforcing continuity at a trivial crossing 120 5.7 Conclusions 125 APPENDICES 128 APPENDIX A: Codable expressions for NPI Time-Derivative Couplings in DGAS 129 APPENDIX B: Mathematical expressions for the first derivative coupling in DGAS 130 APPENDIX C: First and second derivative couplings for DGAS simulations of ethylene 132 REFERENCES 134  CHAPTER 6: POLARONIC RELAXATION BY THREE-ELECTRON BOND FORMATION IN GRAPHITIC CARBON NITRIDES 139 6.1 Introduction 139 6.2 Methods 142 6.2.1 Electronic structure calculations on molecular models 142 6.2.2 Analysis of nuclear relaxation 144 6.2.3 Periodic DFT calculations 145 6.3 Three-electron bond formation in melam/melem dimer cations 146 6.4 Three-electron bond formation in larger model cations 151 6.5 Error assessment for DFT energies 155 6.6 Charge carrier mobility 156 6.7 Lone pair interactions in extended systems 158 6.8 Conclusions 160 APPENDICES 162 APPENDIX A: Supplementary content 163 APPENDIX B: Molecular geometries for models of g-CN 175 REFERENCES 188 vii  CHAPTER 7: BIPOLARONIC RELAXATION IN GRAPHITIC CARBON        NITRIDES  195 7.1 Introduction 195 7.2 Methods 197 7.2.1 Electronic structure calculations 197 7.2.2 Decomposition of the dicationic relaxation energy 198 7.2.3 Nuclear relaxation and charge transport 200 7.3 Two-electron bond formation and bipolaronic relaxation 203 7.4 N1-N1 bond formation and aromatic stabilization 208 7.5 Polaronic and bipolaronic hopping rates 211 7.6 Discussion of g-CN charge transport incorporating recent experimental work 214 7.7 Conclusions 217 APPENDICES 218 APPENDIX A: Supplementary content 219 APPENDIX B: Molecular geometries for dicationic models of g-CN 223 REFERENCES 230    viii  LIST OF TABLES  Table 6.1 Cationic nuclear relaxation data for models of g-CN 147 Table 6.2 Inner sphere component of the polaronic reorganization energy 157 Table A6.1 Vertical and adiabatic ionization potentials for models of g-CN 163 Table A6.2 Neutral melam vibrational frequencies 167 Table A6.3 Cationic melam vibrational frequencies 169 Table A6.4 Neutral melem dimer vibrational frequencies 171 Table A6.5 Cationic melem dimer vibrational frequencies 173 Table 7.1 Structural and bond order data for models of g-CN 205 Table 7.2 Nuclear relaxation energies for models of g-CN with solvent dielectric 207 Table 7.3 Quantities related to the aromatic stabilization energy 209 Table 7.4 Inter-monomeric dihedral angles for models of g-CN 210 Table 7.5 Polaronic self-exchange reaction reorganization energies 211 Table 7.6 Free energy of self-exchange for a polaron in solvated g-CN 213 Table A7.1 Dicationic melam vibrational frequencies 219 Table A7.2 Dicationic melem dimer vibrational frequencies 221    ix  LIST OF FIGURES  Figure 2.1 Scenarios for integration of the time-derivative coupling 9 Figure 2.2 NPI population transfer error and time step 15 Figure 2.3 NPI population transfer error and offset time step 15 Figure 2.4 NPI population transfer error and diabatic coupling strength 18 Figure A2.1 NPI population transfer error with a focus on short time steps 26 Figure A2.2 NPI population transfer error with a focus on small coupling strengths 26 Figure A2.3 Adiabatic PESs for a three state Landau-Zener (LZ) model 27 Figure A2.4 NPI population transfer in a weakly coupled three state LZ model 27 Figure A2.5 NPI population transfer in a strongly coupled three state LZ model 28 Figure 3.1 Illustration of a conical intersection and a trivially unavoided crossing 36 Figure 3.2 AIMS-NPI test systems 40 Figure 3.3 AIMS-NPI population transfer in propene 45 Figure 3.4 AIMS-NPI anionic monomer charge in a long range separated dimer 46 Figure 3.5 AIMS-NPI population transfer in a long range separated dimer 47 Figure 3.6 AIMS-NPI long range charge transfer error and time step size 49 Figure 3.7 NPI and time-averaging of the time-derivative coupling 55 Figure A3.1 AIMS-NPI results for individual simulations of propene 59 Figure A3.2 AIMS results for individual simulations of propene 59 Figure A3.3 AIMS results (without adaptive integration) for simulations of propene 60 Figure A3.4 AIMS-NPI results for trifluoromethyl dimer (0.12 fs time step) 60 Figure A3.5 AIMS-NPI results for trifluoromethyl dimer (0.24 fs time step) 61 x  Figure A3.6 AIMS-NPI results for trifluoromethyl dimer (0.48 fs time step) 61 Figure 5.1 DGAS test systems 113 Figure 5.2 DGAS population transfer in ethylene 115 Figure 5.3 Ethylene pyramidalized conical intersection PES 115 Figure 5.4 DBOC in adiabatic+DBOC simulations of ethylene 116 Figure 5.5 Energy gap in DGAS and adiabatic+DBOC simulations of ethylene 116 Figure 5.6 DBOC and gradient error in the two-state approximation for ethylene 119 Figure 5.7 DGAS and adiabatic AIMS population transfer in trifluoromethyl dimer 121 Figure 5.8 Long range separated ethylene dimer cation simulation results 123 Figure 5.9 Long range separated ethylene dimer cation simulation results 124 Figure A5.1 First derivative couplings in DGAS simulations of ethylene 132 Figure A5.2 Second derivative couplings in DGAS simulations of ethylene 133 Figure 6.1 Phases of graphitic carbon nitrides (g-CN) 140 Figure 6.2 Molecular models of g-CN 141 Figure 6.3 Cationic melam and melem dimer PESs 148 Figure 6.4 Minimum energy structures for neutral and cationic models of g-CN 149 Figure 6.5 Molecular orbital diagram for Pauli-repulsive n-type orbitals in g-CN 151 Figure 6.6 Cationic relaxation energies and nitrogen atoms in models of g-CN 154 Figure 6.7 Band structures for the 2-D and polymeric phases of g-CN 159 Figure 6.8 Kohn-Sham orbitals for the valence and conduction bands of g-CN 160 Figure A6.1 Neutral melam PES 163 Figure A6.2 Polymeric melem trimer cationic LUMO (at neutral minimum) 164 Figure A6.3 Polymeric melem trimer cationic LUMO (at cationic minimum) 164 xi  Figure A6.4 2-D melem trimer cationic LUMO (at neutral minimum) 164 Figure A6.5 2-D melem trimer cationic LUMO (at cationic minimum) 165 Figure A6.6 (2-D melem hexamer cationic LUMO (at neutral minimum) 165 Figure A6.7 (2-D melem hexamer cationic LUMO (at cationic minimum) 166 Figure A6.8 Melam neutral minimum vibrational frequencies 168 Figure A6.9 Melam cationic minimum vibrational frequencies 170 Figure A6.10 Melem dimer neutral minimum vibrational frequencies 172 Figure A6.11 Melem dimer cationic minimum vibrational frequencies 174 Figure 7.1 Hypothesized mechanism for bipolaron formation in g-CN 195 Figure 7.2 Models and atomic labels for the analysis of aromatic stabilization 200 Figure 7.3 Dicationic Kohn-Sham molecular orbitals for all models of g-CN 204 Figure 7.4 Dicationic melam and melem dimer PESs 207 Figure 7.5 Dicationic melam DFT and coupled cluster theory PESs 208 Figure 7.6 Adiabatic ionization potentials and dielectric for models of g-CN 214 Figure A7.1 Melam dicationic vibrational frequencies 220 Figure A7.2 Melem dimer dicationic vibrational frequencies 222   1  CHAPTER 1: INTRODUCTION  1.1 Motivation Increasing energy demands and environmental concerns over fossil fuel use have heightened interest in sources of energy that are clean and fully renewable, such as sunlight. This interest has led to a search for materials which efficiently convert sunlight into usable energy.    Some of the simplest and cheapest materials with this capability are silicon-based solar cells.  Due in large part to low manufacturing costs, silicon solar cells have dominated the market for commercial photovoltaics, however typical silicon solar cells have low efficiencies.1  In many cases these low efficiencies are a consequence of structural impurities.2-3  Structural impurities decrease the efficiency of a device by providing sites for nonradiative decay4-5, in which the energy from absorbed light is lost through motion of the nuclei.  For example, recent work suggests that silicon nanoparticles may exhibit diminished photoluminescence as a result of bending and subsequent nonradiative decay at Si=O surface defects.6  Identifying structural defects that facilitate nonradiefficiency for a desired application. Nonradiative decay occurs when the nuclei relax to a structure where two or more electronic states have the same energy at a given nuclear configuration.  Computational modeling of potential energy surfaces, which express the electronic state energies as a function of the nuclear coordinates, has revealed that these degenerate electronic states intersect with a conical topology.  As such, these intersections have been termed conical intersections (CIs).  Recent work has shown that CIs may facilitate nonradiative decay6, suggesting that modeling molecular motion near CIs is a useful approach for understanding defects that limit material efficiencies.  One powerful tool 2  for identifying CIs is molecular simulation, in which the motions of the nuclei are modelled dynamically.  In most cases these simulations are performed on a single electronic state, so that one assumes the nuclear and electronic motions do not influence one another.  In order to study nonradiative decay, however, one applies nonadiabatic simulation methods, where the term ctronic and nuclear degrees of freedom are coupled to one another. Nonadiabatic approaches generally allow dynamical simulation near CIs through the solution of the time-dependent Schrodinger equation.  These simulations are challenging, however, due to mathematical complications that arise when electronic states intersect.  Established representations for the electronic wave function either have singularities7-8 or are not uniquely defined at CIs.  These representations are the adiabatic and quasi-diabatic representations, respectively.   Singularities in the adiabatic representation result from discontinuities in the electronic wave functions with respect to nuclear motion.  One can avoid these complications through the use of quasi-diabatic representations, but these representations are system-specific and not straight-forward to obtain.  Thus, when performing nonadiabatic simulations one must choose between two difficult representations for the electronic wave functions. In the present work we resolve considerable challenges that arise from dynamical modeling in the adiabatic representation near CIs.  We choose to focus on this representation since these electronic states are obtained straight-forwardly via the solution of the time-independent Schrodinger equation.  Unfortunately, adiabatic electronic wave functions are discontinuous with respect to both time and nuclear motion at CIs.  We will thus begin this work by presenting an approach for dealing with discontinuities in the electronic wave functions with respect to time in dynamical simulations.  This approach is known as the norm-preserving interpolation (NPI) 3  method.  We will then turn our attention to discontinuities with respect to nuclear motion.  First, we demonstrate that singularities arising from these discontinuities are non-integrable in the adiabatic representation.  Second, we present a new variant of the ab initio multiple spawning adiabatic surfaces of the electronic wave function (at CIs) with respect to nuclear motion, however, the wave function is still discontinuous with respect to time.  The NPI approach allows for an accurate treatment of this remaining discontinuity.  In tandem, NPI and DGAS ensure continuity of the molecular wave function in nonadiabatic simulations which use the adiabatic representation. Continuity in the molecular wave function improves the accuracy of dynamical simulations, thereby enhancing our ability to accurately model the process of nonradiative decay. While molecular simulations provide useful information about the dynamics of a molecule, time-independent electronic structure calculations for larger systems also provide a valuable means of studying efficiency-limiting processes in materials.  In contrast to small molecules, in which absorption of light results in an electronically excited state, absorption of light in a material y become uncorrelated and separate from one another.  If charge carriers become trapped, at a defect, for example, they may facilitate nonradiative recombination, in which the electron and hole recombine.  Trapped charge carriers are thus expected to limit efficiency in photovoltaic devices.  In catalytic materials, on the other hand, trapped charge carriers may play a catalytic role.  The desired application thus determines ideal charge carrier behavior.  This behavior is most commonly assessed through experimental observations of charge carrier lifetimes.  Experimental approaches 4  are often ill-equipped, however, to observe isolated surface features on a material.  Time-independent electronic structure calculations can provide this missing microscopic picture of charge carrier behavior. Herein we investigate the charge carrier behavior of a class of photocatalytic materials known as graphitic carbon nitrides9.  Graphitic carbon nitrides (g-CN) are unique in that they are fully organic and exhibit a host of desirable properties for photocatalysis, in which the energy from light absorption is used to catalyze a chemical reaction.  Unfortunately, the efficiency of g-CN is currently lower than other materials which contain precious metals.  A better understanding of efficiency-limiting processes in g-CN would thus be beneficial for the optimization of new fully organic devices with higher photocatalytic efficiencies.  The size of suitable molecular models for g-CN, which are up to two nanometers in diameter, is prohibitive for the application of common electronic structure methods to quantify the energy on a standard central processing unit (CPU).  We overcome this obstacle through the application of graphical processing unit (GPU)-accelerated electronic structure calculations, which offer a considerable speed-up over CPU-based calculations. In review, the present work has two primary foci: (1) to present new methods for the accurate simulation of nonradiative decay in the adiabatic representation, and (2) to use electronic structure calculations to understand the charge carrier behavior in g-CN.  1.2 Description of Chapters Chapters 2 and 3 discuss the NPI approach for handling discontinuities in the electronic wave function with respect to time.  We show that NPI improves the accuracy of computed population transfer rates in nonadiabatic simulations, particularly in systems which exhibit weakly 5  coupled electronic states.  In Chapter 4 we address the non-integrability of singularities at CIs in the adiabatic representation.  We show that these singularities are a consequence of discontinuities in the adiabatic wave functions.  In Chapter 5 we introduce the DGAS variant of AIMS, which removes these discontinuities.   We quantify the simulation error resulting from discontinuous adiabatic wave functions through comparison of the rates of nonradiative decay in DGAS and adiabatic AIMS simulations.  We demonstrate DGAS improves accuracy and enhances our ability to study more challenging problems in photochemistry. In Chapters 6 and 7 we investigate the transport of a single hole and two free holes in molecular models of g-CN.  We show that free holes may localize between nitrogen lone-pairs which are otherwise repulsive.  This localization results in the formation of a new three-electron bond when a single hole is localized, or a two-electron bond when two holes are localized.  It is expected that this mechanism for trapping free holes will influence the photocatalytic behavior of g-CN, most likely by decreasing the photocatalytic efficiency. This work addresses a variety of challenges in the computational modeling of materials for solar energy conversion.  It is unified by the goal of developing and applying new theoretical approaches that enhance our understanding of structural features that limit device efficiency.   6             REFERENCES    7  REFERENCES   (1) Shah, A.; Torres, P.; Tscharner, R.; Wyrsch, N.; Keppner, H. Photovoltaic technology: The case for thin-film solar cells. Science 1999, 285, 692-698.   (2) Galland, C.; Ghosh, Y.; Steinbruck, A.; Sykora, M.; Hollingsworth, J. A.; Klimov, V. I.; Htoon, H. Two types of luminescence blinking revealed by spectroelectrochemistry of single quantum dots. Nature 2011, 479, 203-U275.   (3) Kamat, P. V. Boosting the Efficiency of Quantum Dot Sensitized Solar Cells through Modulation of Interfacial Charge Transfer. Accounts of Chemical Research 2012, 45, 1906-1915.   (4) Hall, R. N. ELECTRON-HOLE RECOMBINATION IN GERMANIUM. Physical Review 1952, 87, 387-387.   (5) Shockley, W.; Read, W. T. STATISTICS OF THE RECOMBINATIONS OF HOLES AND ELECTRONS. Physical Review 1952, 87, 835-842.   (6) Shu, Y. N.; Fales, B. S.; Levine, B. G. Defect-Induced Conical Intersections Promote Nonradiative Recombination. Nano Letters 2015, 15, 6247-6253.   (7) Saxe, P.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS FOR MCSCF/CL WAVE-FUNCTIONS .4. 2ND DERIVATIVE TERMS USING ANALYTIC GRADIENT METHODS. Journal of Chemical Physics 1987, 86, 321-328.   (8) Gherib, R.; Ryabinkin, I. G.; Izmaylov, A. F. Why Do Mixed Quantum-Classical Methods Describe Short-Time Dynamics through Conical Intersections So Well? Analysis of Geometric Phase Effects. Journal of Chemical Theory and Computation 2015, 11, 1375-1382.   (9) Wang, X.; Blechert, S.; Antonietti, M. Polymeric Graphitic Carbon Nitride for Heterogeneous Photocatalysis. Acs Catalysis 2012, 2, 1596-1606.    8  CHAPTER 2: EVALUATION OF THE TIME-DERIVATIVE COUPLING FOR ACCURATE ELECTRONIC STATE TRANSITION PROBABILITIES FROM NUMERICAL SIMULATIONS  This chapter is taken in its entirety from: Evaluation of the time-derivative coupling for accurate electronic state transition probabilities from numerical simulations. Phys. Chem. Lett., 2014, 5, 2351-2356.  2.1 Introduction  The accurate theoretical modeling of molecular dynamics on multiple electronic states has allowed the elucidation of the ultrafast motions following the electronic excitation of molecules and materials.   Classical-trajectory-1-16 and wavepacket-based17-22 approaches allow the simulation of such phenomena on accurate potential energy surfaces (PESs) determined by ab initio calculations performed on-the-fly.  Both of these classes of methods generally rely on the numerical integration of the time-dependent Schrodinger equation (TDSE), which in turn requires the discretization of time into steps of finite size.  Numerical schemes for the integration of differential equations like the TDSE are well known,23 but such approaches only provide accurate results when the terms in the equation vary on timescales shorter than the time step.  Unfortunately, the simulation of nonadiabatic dynamics in the adiabatic representation requires the treatment of population transfer through conical intersections, points of degeneracy between electronic states.24 Such intersections introduce spikes in the time-derivative coupling (TDC), defined ,  (2.1) where  is the electronic wave function of adiabatic state j.  These spikes may be infinitesimally narrow in time and infinitely large in magnitude, and their accurate integration is essential to a correct prediction of the probability of population transfer.25-27    Working in the diabatic 9  representation can reduce these issues, but no unique diabatic representation exists, and the accuracy of simulations performed in the diabatic representation can be lower than those based on the adiabatic representation.7,28  Figure 2.1 Scenarios for integration of the time-derivative coupling Cartoon representation of the cases where the simulation steps over (left) or steps on (right) the maximum of the TDC.  Dotted lines indicate calculations of the TDC.  Approaches to numerically integrating the TDSE in the adiabatic basis can be split into two broad categories: schemes based on analytic calculations of the nonadiabatic coupling matrix element (NACME) vector, and those based on the numerical differentiation of the wave function in time.  The NACME , (2.2) where r and R represent the electronic and nuclear coordinates, respectivelycan now be computed for many approximate electronic structure methods using analytical gradient techniques.29-37  Treating nuclear motion classically, the TDC matrix element can then be computed according to 10  , (2.3) where  represents the time derivative of R.  Though efficient and widely used, this approach suffers when the TDSE is numerically integrated with a finite time step, as illustrated in Figure 2.1.  When  spikes on a timescale shorter than the integration time stepe.g. at a conical intersection or weakly avoided crossingthe simulation may step on the spike, resulting in a coupling which is far larger than the average coupling over the time step, or step over the spike, resulting in an erroneously small coupling.  Such errors in the determination of the coupling translate to errors in the population transfer predicted by integration of the TDSE.  These errors may be especially large in cases where systematically weak coupling between diabats leads to (N-1)-dimensional trivial unavoided crossings, such as in long-range energy transfer,25-27 long-range charge transfer,38 and intersystem crossing resulting from weak spin-orbit interactions.  Therefore, though analytical derivative techniques provide the exact TDC at discrete times, discretization of time itself is a coarse approximation which can lead to very large errors.  Adaptive time step integration (as implemented in ab initio multiple spawning; AIMS) alleviates this problem,38,39 but at a significant computational cost.  An appealing alternative is to use information about the electronic wave function at the beginning () and end () of the time step to approximate the TDC.  In such approaches, the true change in the electronic wave functions over the entire time step (rather than the derivative at a single time) is resolved by computing overlap integrals between the adiabatic wave functions at times  and .  Hammes-Schiffer and Tully (HST) proposed such an overlap-based approach,40 approximating the TDC according to . (2.4) 11  Though generally thought of as an approximation to the analytical scheme, the HST method has been demonstrated to predict a similar probability for population transfer upon averaging over large ensembles of surface hopping trajectories.41,42  The HST and other overlap-based approaches have the additional advantage that they can be applied in conjunction with electronic structure methods for which analytic NACMEs are not implemented,43 and in other cases can be less expensive to compute than such analytic calculations.41,42  In a similar spirit, several overlap-based schemes have been developed for directly determining the hopping probability for surface hopping simulations, bypassing the TDC altogether.10,25-27  These approaches have been shown to address the trivial unavoided crossing problem, but because they do not provide an explicit definition for the TDC they cannot be applied in wavepacket-based dynamical schemes.  The purpose of the present work is to develop an approach for approximation of the TDC that is accurate with standard simulation time steps in wavepacket-based dynamical schemes.    2.2 The norm preserving interpolation (NPI) method  Given the advantages of overlap-based approaches, we propose a new such approach, the norm-preserving interpolation (NPI) method.  Because NPI provides a recipe for the explicit calculation of the TDC it can be applied in conjunction with both wavepacket-based and surface hopping methods.  We will analyze the accuracy of simulations utilizing NPI compared to those based on analytical derivative and HST approaches by applying these methods to integrate the TDSE for the adiabatized Landau-Zener (LZ) model,44,45 for which the exact solution is known. In the NPI approach, we approximate the electronic wave function of adiabatic state  at time, , on the interval  as the product of a time-dependent transformation matrix, , with the wave function at the beginning of the time step, , 12  . (2.5) The diagonal and off-diagonal elements of  are respectively defined  (2.6) and . (2.7) with  . (2.8) Note that in the complete Hilbert space  is a unitary transformation for all , and thus Eq. (2.5) represents an interpolation of the adiabatic electronic wave function between its known values at the beginning and end of the time step,  and , that maintains the normalization of the wave function at all times, . We can thus approximate the TDC between states k and j within each time step as . (2.9) Averaging (2.9) over the interval  yields the NPI approximation to the TDC at the center of the time step, . (2.10) In a sense, this is not an approximation to the TDC at a specific point in time so much as a value of the TDC that represents the full change in the wave function over the time step.  See appendices 13  for integrated expressions and for comments about proper numerical implementation.  Though the expressions are somewhat more complicated than Eq. (2.4), they require essentially the same information from the electronic structure calculation: the overlaps of the adiabatic wave functions at the beginning of each time step with those at the end.  2.3 Methods  The Landau-Zener44,45 two-state model has been chosen as a test case for the performance of the NPI method because the exact solution for the population transfer probability is known, thus allowing for direct quantification of errors.  In the Landau-Zener two-state model the time-dependent energies of diabatic electronic states  and  are expressed as  and , respectively, where  is the slope of the potential energy surfaces.  These states are subject to a time-independent coupling, .  A set of time-dependent adiabatic states are defined by diagonalizing the resulting 2x2 Hamiltonian.  The exact probability of population transfer between adiabats is known to be .  All numerical simulations of population transfer performed here have been compared with this exact solution to determine the population transfer probability error, . The TDSE is numerically integrated to determine the population transfer probability.  Three approximations to the TDC have been tested.  The analytical scheme corresponds to the exact value of the TDC defined by Eq. (2.1) at each time step.  The NPI and HST schemes are computed according to Eqs. (2.10) and (2.4), respectively.  Integration of the TDSE has been performed using the fourth-order Runge-Kutta method, with a multiple time step algorithm similar 14  to that employed in the software implementation of AIMS.22  Additional details of the numerical integration scheme are presented in Appendix A.  2.4 Results In Figure 2.2 the error in the computed population transfer probability is plotted as a function of the length of the time step for the passage through a weakly avoided crossing .  These parameters were drastudy of the non-radiative decay of silicon epoxide defects, where similar slopes and several nonadiabatic events involving a gap of 10 meV or smaller were observed.46  Note that a similar size coupling was employed in a recent study of the behavior of surface hopping at trivial unavoided crossings as well.27   The exact LZ population transfer in this case is 0.995.  For time steps in the range typically chosen for realistic dynamical simulations (0.25-0.5 fs) the NPI method results in a smaller error than the HST and analytical schemes.  At a time step of 0.5 fs, the error in the computed population transfer probability from the NPI-based simulations has a negligible value of -0.001.  Application of the HST scheme yields a larger magnitude error of -0.093, while the analytical method results in an error of very large magnitude: -0.457.  At longer steps, NPI begins to yield noticeable errors, which we attribute to an inaccurate treatment of dephasing.  When a very short time stt < 0.15 fs) the error in all methods is small (), but it is interesting to note that the NPI scheme still yields the smallest error in these cases (Figure A2.1).   15   Figure 2.2 NPI population transfer error and time step The population transfer probability error,, for each method of evaluating the time-derivative coupling in the steps-over case is presented as a function of the length of the time step.   Figure 2.3 NPI population transfer error and offset time step The population transfer error is plotted as a function of the offset of the time step from the maximum of the TDC.  The symmetrical steps-over case falls at 0.25 fs, while the steps-on case falls at 0.0 and 0.5 fs. 16  It is important to note that in all of the simulations reported in Figure 2.2 the initial time was set such that the maximum of the TDC falls at the center of a time step, whereas the calculation of the potential, wave function, and (for the analytical method) the analytically computed TDC are performed at the beginning and end of the time step.  In more graphical language, the simulations in Figure 2.2 correspond to the case where the discretized trajectories step over the intersection.  In this case (weakly avoided crossing, stepping over the intersection), HST and NPI have the advantage of fully resolving the change in the electronic wave function, and thus exhibit smaller errors.  In contrast, the analytical simulations do not resolve the spike in the TDC, resulting in large errors.  Of course whether the discretized trajectory steps on or steps over the intersection in a dynamical simulation is determined by chance, and thus we consider how the population transfer error depends on the offset of the maximum in the TDC relative to the beginning of the time step.  In Figure 2.3, we show the error in the simulated population transfer probability for the same LZ model as a function of this offset.  All simulations were conducted with a 0.5 fs time step, so a 0.25 fs offset corresponds to the case where the trajectory symmetrically steps over the spike in the coupling (Figure 2.1 left), 0.0 and 0.5 fs offsets correspond to the steps-on case (Figure 2.1 right), and the other offsets correspond to the continuum of possibilities between these two extremes.   The NPI scheme yields population transfer probability errors between -0.001 and 0.000 for all offsets.  In contrast, the error in the population transfer varies wildly when analytical TDCs are employed, changing from -0.865 in the steps-on case, to nearly zero at a 0.10 fs offset, to -0.457 in the previously discussed steps-on case.  This is because the analytically calculated TDC is far too large at the steps-on point, and far too small at the steps-over point.  One might expect the error 17  to average out over many simulations, but this is not the case for a weakly avoided crossing.  At a weakly avoided crossing, the probability of population transfer is near unity.  Any error, whether it be that the TDC is too large or too small, results in a reduction of the predicted population transfer probability.  Thus, at weakly avoided crossings, the analytic scheme will systematically underpredict the probability of nonadiabatic population transfer.  In the context of realistic simulations, this could lead to, for example, the prediction of unphysical long-range charge and energy transfer events where diabatic motion through an intersection would be physically correct. Again, the HST scheme results in smaller population transfer errors than the analytical approach, but the errors rise as large as -0.093 at the steps-over point.  The errors in the HST approach can be understood by demonstrating how the HST definition of the TDC (Eq. (2.4)) can be derived from the same unitary formalism used to derive NPI.  The HST expression arises from the approximation .  (2.11) This expression differs from the NPI definition of the TDC (Eq. (2.10)) in that only the ket wave function is interpolated in time, while the bra wave function is frozen at  and  in the first and second terms, respectively.  These definitions are essentially equal when  is very small, but in the case that  approaches 1, the definitions of the TDC in Eqs. (2.4) and (2.10) approach  and , respectively, and only integration of the latter value results in the prediction of the correct LZ transition probability for very weakly avoided crossings. 18   We continue by investigating the utility of NPI at more strongly avoided crossings.  Figure 2.4 presents the error in the population transfer as a function of the coupling between diabats, .  In all cases a time step of 0.5 fs is used,  is 0.1 eV/fs, and the offset is chosen such that the trajectory steps over the intersection symmetrically.  For all values of  the use of NPI results in the smallest population transfer errors, with all errors between -0.007 and 0.000.  For larger values of  (strongly avoided crossings) HST and the analytical method both provide suitable accuracy, but as the coupling strength decreases these approaches yield significant errors.  Similar results are obtained when we consider the steps-on case, though HST yields considerably smaller errors in this case (see Figure A2.2 for details).  Figure 2.4 NPI population transfer error and diabatic coupling strength The population transfer error is plotted as a function of the coupling strength between the diabatic states, .    19  2.5 Conclusions  Thus, by using a time-dependent unitary transformation to interpolate the adiabatic electronic wave function across each time step and computing the TDC between these interpolated wave functions, NPI is able to provide TDCs which result in the accurate integration of the TDSE even at very weakly avoided crossings.  When applied in numerical simulations based on the adiabatized LZ model, the NPI approach results in population transfer errors of ~0.001 relative to the exact LZ solutions over a wide range of time steps and LZ model parameters, whereas application of the HST and analytical approaches results in much larger errors which depend strongly on both the choice of time step and the model parameters.  In addition, the population transfer probability at weakly avoided crossings predicted by NPI-based simulations does not depend strongly on the offset of the spike in the TDC relative to the beginning of the time step, whereas the HST and analytical schemes yield errors which are strongly dependent on this uncontrollable factor.  For the majority of offsets and realistic time steps, simulations based on HST yield smaller population transfer errors at weakly avoided crossings than do those based on analytical calculation of the TDC.  Lastly, NPI is trivially implemented into surface hopping and wavepacket-based simulation schemes, and does not require the analytical calculation of NACMEs.  Thus NPI can be employed with a wide range of electronic structure methods for which such analytic calculations have not yet been implemented.   20             APPENDICES    21  APPENDIX A: Integrated expressions for the time-derivative coupling The below assumes real wave functions, as are normally provided by electronic structure calculations.  For the sake of brevity, we define . (2.12) The integrated expression for the NPI TDC, , (2.13) can be shown to be  (2.14) where , (2.15) , (2.16) , (2.17) , (2.18) and , (2.19) 22  except in the case of a two-dimensional Hilbert space, in which .  In Eq. (2.19),  indexes the unit vector parallel to the component of  that is orthogonal to both  and :   The  term arises because the operator  spans the full Hilbert space, not just the space spanned by  and .  It can be shown that the elements  and  are   (2.20) and . (2.21)  In the special case that  is zero, we note that , and thus there is no need to compute . When numerically implementing Eq. (2.14) one must be careful of the denominators in the terms A, B, C, and D, which are all of the form .  When x approaches zero, the limiting expression  (2.22) is applied.  Similarly, the denominator in term E can approach zero.  In this case, we apply the limiting expression . (2.23)   23  APPENDIX B: Numerical Integration Details  In the LZ model, the electronic states depend explicitly on time, and thus there are no classical degrees of freedom to integrate.  Integration of the TDSE, , (2.24) has been performed using a multiple time step algorithm based on the fourth-order Runge-Kutta (RK) integrator, developed for use with ab initio multiple spawning.47  This approach was chosen to most closely model the application of NPI in realistic simulations.  In this algorithm, separate large and small time steps employed.  The energy, wave function, and (if used) analytic TDC are computed at the beginning of every large time step (referred to as  in the main text).  This larger time step is divided into smaller time steps on which the RK integrator is applied.  These RK time steps are adjusted adaptively.  In the context of an ab initio molecular dynamics run, this approach ensures the conservation of the norm of the wave function without requiring a large number of electronic structure calculations. Because the energy, wave function, and analytic TDC are only computed at the beginning of each large time step, it is necessary to define their values on the shorter RK time steps.  Both numerical time-derivative coupling methods discussed in the main text (HST and NPI) provide constant valued nonadiabatic couplings to be applied over the length of the long time step: . (2.25) For the analytical method, in which the TDC is evaluated exactly at the beginning of each time step, , the TDC is applied during the time-step as follows: . (2.26) 24  The same approach is taken for the energy: . (2.27)  Initially, all population is placed in adiabatic state .  The population of adiabatic state at the end of the simulation is taken as the population transfer probability.  All simulations have been run for 400 fs, with the crossing region occurring at a time of 200 fs.  This time has been chosen in conjunction with the slope for the diabatic state energies, 0.1 eV/fs in all cases, to ensure that the interaction between states is negligible at both the beginning and end of the simulation.    25  APPENDIX C: Three state Landau-Zener model results  Figures A1.3-A1.5 show results pertaining to tests of the NPI/HST TDC with three interacting states.  Figure A2.3 shows the adiabatic potential energy surfaces for one example three state LZ model of this study, in which the crossings between three states are trivially unavoided.  The energies of the diabatic states in this model are expressed as: . (2.28) We note that in this model the diabatic energies are set such that diabat  does not intersect with diabats and precisely at an integer value of time. Thus the offset of the spike in the TDC from the beginning of the time-step in this model is uncontrolled.  Initially all population is placed in adiabatic state .  The population of adiabatic state at the end of the simulation is plotted as a function of the length of the time step in Figure A2.4 for the case of trivial unavoided crossings, in which the coupling strength, , between all states is set to 0.01 eV.  The oscillatory behavior of HST in this case is likely due to variations in the offset of the spike in the TDC from the beginning of the time step, as previously discussed.  Figure A2.5 plots the final occupation probability of adiabatic state  as a function of the length of the time step for a strong coupling case, in which the coupling strength, , between all states is set to 0.25 eV.  All simulations were run for 400 fs.  This time has again been chosen in conjunction with the slope for the diabatic state energies, to ensure that the interaction between states is negligible at both the beginning and end of the simulation.     26   Figure A2.1 NPI population transfer error with a focus on short time steps The population transfer probability error, , for each method of evaluating the nonadiabatic coupling is presented as a function of the length of the time step.  This plot is similar to Figure 2.1, but focuses on time steps of less than 0.20 fs.   Figure A2.2 NPI population transfer error with a focus on small coupling strengths The population transfer error is plotted as a function of the coupling strength, , in the case that the f the TDC.  The inset shows the same results for the analytical method but is focused on smaller coupling strengths.  A time step of 0.50 fs is used and the slope of the potential energy surfaces, , is 0.1 eV/fs.   27   Figure A2.3 Adiabatic PESs for a three state Landau-Zener (LZ) model The adiabatic potential energy surfaces are shown for the three state LZ model discussed here.  The coupling strength,, has been set to 0.01 eV.  The parameter corresponding to the slope of the diabatic potential energy surfaces, , is set to 0.1 eV/fs.  These adiabatic energies correspond to the model for which occupation probabilities are presented in Figure A2.4.   Figure A2.4 NPI population transfer in a weakly coupled three state LZ model The final occupation probability for adiabatic state , , is shown for each numerical method of evaluating the TDC in the three state LZ model shown in Figure A2.3.  The occupation probability is plotted versus the length of the time step for the simulation.  These results correspond to the case of trivially unavoided crossings between adiabats, in which the coupling strength, , between all states has been set to 0.01 eV. 28   Figure A2.5 NPI population transfer in a strongly coupled three state LZ model The final occupation probability for adiabatic state , , is shown for each numerical method of evaluating the TDC in the three state LZ model.  The final occupation probability is plotted versus the length of the time step.  In these simulations the coupling strength, , between all states has been set to 0.25 eV.   29             REFERENCES    30  REFERENCES   (1) Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions - Reaction of H+ with D2. J. Chem. Phys. 1971, 55, 562-572.   (2) Micha, D. A. A Self-Consistent Eikonal Treatment of Electronic-Transitions in Molecular-Collisions. J. Chem. Phys. 1983, 78, 7138-7145.   (3) Tully, J. C. Molecular-Dynamics with Electronic-Transitions. J. Chem. Phys. 1990, 93, 1061-1071.   (4) Bittner, E. R.; Rossky, P. J. Quantum Decoherence in Mixed Quantum-Classical Systems - Nonadiabatic Processes. J. Chem. Phys. 1995, 103, 8130-8143.   (5) Coker, D. F.; Xiao, L. Methods for Molecular-Dynamics with Nonadiabatic Transitions. J. Chem. Phys. 1995, 102, 496-510.   (6) Martens, C. C.; Fang, J. Y. Semiclassical-Limit Molecular Dynamics on Multiple Electronic Surfaces. J. Chem. Phys. 1997, 106, 4918-4930.   (7) Prezhdo, O. V.; Rossky, P. J. Mean-Field Molecular Dynamics with Surface Hopping. J. Chem. Phys. 1997, 107, 825-834.   (8) Tully, J. C. Mixed Quantum-Classical Dynamics. Faraday Discuss. 1998, 110, 407-419.   (9) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919-8929.   (10) Granucci, G.; Persico, M.; Toniolo, A. Direct Semiclassical Simulation of Photochemical Processes with Semiempirical Wave Functions. J. Chem. Phys. 2001, 114, 10608-10615.   (11) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-Field Dynamics with Stochastic Decoherence (MF-SD): A New Algorithm for Nonadiabatic Mixed Quantum/Classical Molecular-Dynamics Simulations with Nuclear-Induced Decoherence. J. Chem. Phys. 2005, 123, 234106.   (12) Craig, C. F.; Duncan, W. R.; Prezhdo, O. V. Trajectory Surface Hopping in the Time-Dependent Kohn-Sham Approach for Electron-Nuclear Dynamics. Phys. Rev. Lett. 2005, 95, 163001.   (13) Li, X. S.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106.  31   (14) Fischer, S. A.; Chapman, C. T.; Li, X. S. Surface Hopping with Ehrenfest Excited Potential. J. Chem. Phys. 2011, 135, 144102.   (15) Shenvi, N.; Subotnik, J. E.; Yang, W. Simultaneous-Trajectory Surface Hopping: A Parameter-Free Algorithm for Implementing Decoherence in Nonadiabatic Dynamics. J. Chem. Phys. 2011, 134, 144102.   (16) Richter, M.; Marquetand, P.; Gonzalez-Vazquez, J.; Sola, I.; Gonzalez, L. SHARC: Ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. J. Chem. Theory Comput. 2011, 7, 1253-1258.   (17) Burghardt, I.; Meyer, H. D.; Cederbaum, L. S. Approaches to the Approximate Treatment of Complex Molecular Systems by the Multiconfiguration Time-Dependent Hartree Method. J. Chem. Phys. 1999, 111, 2927-2939.   (18) Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. J. Phys. Chem. A 2000, 104, 5161-5175.   (19) Shalashilin, D. V.; Child, M. S. Time Dependent Quantum Propagation in Phase Space. J. Chem. Phys. 2000, 113, 10028-10036.   (20) Worth, G. A.; Robb, M. A.; Burghardt, I. A Novel Algorithm for Non-Adiabatic Direct Dynamics Using Variational Gaussian Wavepackets. Faraday Discuss. 2004, 127, 307-323.   (21) Chen, X.; Batista, V. S. Matching-Pursuit/Split-Operator-Fourier-Transform Simulations of Excited-State Nonadiabatic Quantum Dynamics in Pyrazine. J. Chem. Phys. 2006, 125, 124313.   (22) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of Ab Initio Multiple Spawning in the MOLPRO Quantum Chemistry Package. Chem. Phys. 2008, 347, 3-16.   (23) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in Fortran 77: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992.   (24) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annu. Rev. Phys. Chem. 2012, 63, 325-352.   (25) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of Unavoided Crossings in Nonadiabatic Photoexcited Dynamics Involving Multiple Electronic States in Polyatomic Conjugated Molecules. J. Chem. Phys. 2012, 137, 014512.  (26) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts Due to Trivial Unavoided Crossings in the Modeling of Photoinduced Energy Transfer Dynamics in Extended Conjugated Molecules. Chem. Phys. Lett. 2013, 590, 208-213. 32    (27) Wang, L. J.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. 2014, 5, 713-719.   (28) Neria, E.; Nitzan, A. Semiclassical Evaluation of Nonadiabatic Rates in Condensed Phases. J. Chem. Phys. 1993, 99, 1109-1123.   (29) Lengsfield, B. H.; Saxe, P.; Yarkony, D. R. On the Evaluation of Nonadiabatic Coupling Matrix-Elements Using SA-MCSCF/CI Wave-Functions and Analytic Gradient Methods .1. J. Chem. Phys. 1984, 81, 4549-4553.   (30) Lengsfield, B. H.; Yarkony, D. R. Nonadiabatic Interactions between Potential-Energy Surfaces - Theory and Applications. Adv. Chem. Phys. 1992, 82, 1-71.   (31) Chernyak, V.; Mukamel, S. Density-Matrix Representation of Nonadiabatic Couplings in Time-Dependent Density Functional (TDDFT) Theories. J. Chem. Phys. 2000, 112, 3572-3579.   (32) Ichino, T.; Gauss, J.; Stanton, J. F. Quasidiabatic States Described by Coupled-Cluster Theory. J. Chem. Phys. 2009, 130, 174105.   (33) Tavernelli, I.; Curchod, B. F. E.; Laktionov, A.; Rothlisberger, U. Nonadiabatic Coupling Vectors for Excited States within Time-Dependent Density Functional Theory in the Tamm-Dancoff Approximation and Beyond. J. Chem. Phys. 2010, 133, 194104.   (34) Nelson, T.; Fernandez-Alberti, S.; Chernyak, V.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics Modeling of Photoinduced Dynamics in Conjugated Molecules. J. Phys. Chem. B 2011, 115, 5402-5414.   (35) Fatehi, S.; Alguire, E.; Shao, Y. H.; Subotnik, J. E. Analytic Derivative Couplings between Configuration-Interaction-Singles States with Built-in Electron-Translation Factors for Translational Invariance. J. Chem. Phys. 2011, 135, 234105.   (36) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martinez, T. J. Role of Rydberg States in the Photochemical Dynamics of Ethylene. J. Phys. Chem. A 2012, 116, 2808-2818.   (37) Fatehi, S.; Alguire, E.; Subotnik, J. E. Derivative Couplings and Analytic Gradients for Diabatic States, with an Implementation for Boys-Localized Configuration-Interaction Singles. J. Chem. Phys. 2013, 139, 124112.   (38) Wang, L. J.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-Like Transport in Organic Crystals. J. Phys. Chem. Lett. 2013, 4, 1888-1894.  (39) Virshup, A. M.; Levine, B. G.; Martinez, T. J. Steric and Electrostatic Effects on Photoisomerization Dynamics Using QM/MM Ab Initio Multiple Spawning Theor. Chem. Acc. 2014, under review. 33    (40) Hammes-Schiffer, S.; Tully, J. C. Proton-Transfer in Solution - Molecular-Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657-4667.   (41) Pittner, J.; Lischka, H.; Barbatti, M. Optimization of Mixed Quantum-Classical Dynamics: Time-Derivative Coupling Terms and Selected Couplings. Chem. Phys. 2009, 356, 147-152.   (42) Fabiano, E.; Keal, T. W.; Thiel, W. Implementation of Surface Hopping Molecular Dynamics Using Semiempirical Methods. Chem. Phys. 2008, 349, 334-347.   (43) Tao, H. L.; Levine, B. G.; Martinez, T. J. Ab Initio Multiple Spawning Dynamics Using Multi-State Second-Order Perturbation Theory. J. Phys. Chem. A 2009, 113, 13656-13662.   (44) Zener, C. Non-Adiabatic Crossing of Energy Levels. Proc. R. Soc. Lond. A 1932, 137, 696-702.   (45) Landau, L. Zur Theorie Der Energieubertragung.  II. Phys. Z. Sowjetunion 1932, 2, 46-51.   (46) Shu, Y.; Levine, B. G. Communication: Non-Radiative Recombination Via Conical Intersection at a Semiconductor Defect. J. Chem. Phys. 2013, 139, 081102.  (47) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of ab initio multiple spawning in the MOLPRO quantum chemistry package. Chemical Physics 2008, 347, 3-16.   34  CHAPTER 3: ACCURATE AND EFFICIENT EVALUATION OF TRANSITION PROBABILITIES AT UNAVOIDED CROSSINGS IN AB INITIO MULTIPLE SPAWNING  This chapter is taken in its entirety from: Accurate and Efficient Evaluation of Transition Probabilities at Unavoided Crossings in Ab Initio Multiple SpawningChem. Phys. 2015, 460, 117-124.  3.1 Introduction  Molecular dynamics (MD) simulations with electronic state transitions provide valuable insights into the behavior of photoexcited molecules and materials, including ultrafast processes such as non-radiative decay.  Leading approaches include classical-trajectory 1-22 methods, an -hopping, and fully quantum mechanical 23-32 schemes, including wavepacket-based approaches such as ab initio multiple spawning (AIMS).  Each of these approaches involves integration of the time-dependent Schrödinger equation (TDSE) with a Hamiltonian that includes non-Born-Oppenheimer terms.  Ab initio (AI) MD methodologies, which solve for the electronic structure of the system on-the-fly, can be particularly useful in describing such photochemical processes because they do not require assumptions about the potential energy surface (PES) or the character of the electronic states involved, but computational expense limits the system sizes, time scales, and accuracies accessible via such simulations.    The TDSE is most often solved by numerical integration, which requires the discretization of time into finite steps.  For the sake of numerical accuracy, it is necessary that the chosen time step is shorter than the timescale on which the Hamiltonian matrix elements vary.  However, at each time step computationally expensive electronic structure calculations are performed, and thus efficiency demands the choice of the longest time step possible.  When working in the adiabatic 35  representation, this dilemma can become intractable because the time-derivative coupling (TDC, eqn. (2.1)), may vary rapidly within a single simulation time step.  In modern AIMD simulations the TDC is often computed from the nonadiabatic coupling matrix element (NACME) vector (eqn. (2.2)) and the nuclear velocity via the chain rule (eqn. (2.3)).  NACMEs can now be efficiently calculated at many quantum chemical levels of theory using analytical gradient techniques.33-41    The TDC exhibits infinitesimally narrow spikes of infinite magnitude at conical intersections (CIs), points of degeneracy between PESs known to facilitate non-radiative decay,42-46 and thus simple numerical integration of the TDSE through a CI requires an infinitesimal time step, which is obviously computationally intractable.  In practice, however, CIs are (N-2)-dimensional features of the PES (where N is the total number of nuclear degrees of freedom), so the probability of a particular trajectory hitting one directly (Figure 3.1, arrow 1) is exceedingly low.  Most trajectories will instead pass around the intersection, leading to a more slowly varying TDC (Figure 3.1, arrow 2).  More troublesome are trivial unavoided crossings (TUCs), effectively (N-1)-dimensional intersections which arise in the adiabatic representation when the coupling between two diabatic electronic states is negligible (e.g. when modeling long range charge or energy transfer or intersystem crossing in the limit of weak spin-orbit coupling).  In such cases, direct passage of a trajectory through the intersection seam (Figure 3.1, arrow 3) is unavoidable, and adiabatic motion through the intersection (Figure 3.1, arrow 4) is impossible.  Passage through the TUC necessarily involves an infinitesimally narrow spike in the TDC, which a simulation based on finite time steps is unable to resolve, leading to inaccurate predictions of the probabilities of electronic state transitions.  36   Figure 3.1 Illustration of a conical intersection and a trivially unavoided crossing Illustrations of potential energy surfaces near A) a CI and B) a TUC.  The three green arrows illustrate possible nuclear motion near each crossing: 1) diabatic passage through a CI, 2) adiabatic passage around a CI, and 3) diabatic passage through a TUC.  The red arrow (4) illustrates unphysical adiabatic passage through a TUC, which could manifest itself as e.g. unphysical long range charge transfer in a real simulation.  Several strategies can be proposed to circumvent numerical issues associated with the integration of the TDSE through regions where the TDC spikes sharply.  Working in the diabatic representation eliminates the spike entirely, but the integration of the relatively small, temporally delocalized coupling between diabatic states can also lead to inaccurate prediction of transition probabilities.9,47  In the case of CIs, adaptive time step integration effectively circumvents numerical problems in all but the worst cases, but at a significant computational cost.48,49  In the context of trajectory surface hopping simulations, formulations for the hopping probability have been proposed which accurately predict transition probabilities resulting from sharply spiked 37  TDCs.12,50-53  However, formulations for the hopping probability are not directly applicable in fully quantum mechanical simulation methods such as AIMS.  An alternative to the above approaches is to eschew the analytical calculation of the NACME at a finite number of time steps altogether.  Approaches in which the TDC is computed directly from the overlaps of the adiabatic wave functions at the beginning of each time step with those at the end of each step have long been used when analytical NACME calculations are not available.4,54,55  An often overlooked advantage of these overlap-based approaches is that they guarantee that the full change in the adiabatic wave function over each time step is resolved, even if the instantaneous TDC is varying on timescales much shorter than this step.   In Chapter 2 we developed an overlap-based approach derived from a continuous interpolation of the adiabatic electronic wave function over the time step, the norm-preserving interpolation (NPI) scheme.55  In the present chapter we report the implementation of the NPI scheme into the AIMS nonadiabatic molecular dynamics method and assess its performance by application to two different chemical systems:  propylene, which exhibits a CI, and a long range charge transfer couple, which exhibits a TUC.  In section 3.2 we discuss the implementation of NPI in the AIMS method, and detail the two cases we will use to test the combined AIMS-NPI method.  In sections 3.4-3.7 we present analysis of the population transfer probabilities computed using the analytical and NPI schemes in these cases and discuss the primary sources of error in transition probabilities computed using AIMS-NPI.  Finally, in section 3.8 we draw conclusions and discuss future prospects.  38  3.2 Incorporation of NPI into AIMS The AIMS method is a hierarchy of approximations that allows the full solution of the TDSE, including all electronic and nuclear degrees of freedom explicitly.  The total wave function is defined as a sum over Born-Oppenheimer states   (3.1) where I indexes electronic states, is a time-dependent nuclear wave function associated with state I, and  is the nuclear-position-dependent electronic wave function associated with that state.  In the present work the basis of adiabatic electronic states is used.  Each nuclear wave function, , is expanded in a basis of frozen Gaussian56 trajectory basis functions (TBFs)  (3.2) where  indexes a TBF and  is a TBF which is parameterized by the average position, ,  the average momentum, , the nuclear phase, , and time-.  These TBFs are defined  (3.3) where  indexes nuclear degrees of freedom.  The average position and momentum of each TBF is propagated on the PES of electronic state I according to classical equations of motion and the phase is propagated semiclassically.  The time-dependent expansion coefficients, , are determined by integration of the TDSE. 39  In order to ensure that nonadiabatic events are modeled accurately, the nuclear basis is adaptively expanded via the creation (spawning) of new basis functions.  When a TBF encounters a region of the PES where the TDC between two electronic states is large in magnitude, such as near a CI or TUC, a new TBF is created on the coupled electronic state to accept population. In the most widely used implementation of AIMS the TDC is calculated from analytically-computed NACMEs, and accurate integration of the TDSE near CIs is achieved via an adaptive integrator.49  In this work, we have replaced the analytically computed TDC with the NPI approximation (Eq. (2.10)) in order to achieve accurate integration without the need for computationally expensive adaptive integration.  In doing so, we were required to modify two elements of the AIMS algorithm: the numerical integration procedure for the TDSE and the spawning algorithm.   In the standard integration scheme based on analytical NACMEs, the TDC values are computed at the beginning and end of a time step,  and , respectively.  However, when employing NPI the TDC is approximated at the midpoint of the time step, .  The standard numerical integrator employed in AIMS (described in detail in ref 29) has been modified to be consistent with this difference.  Additionally, the standard integration scheme makes use of adaptive integration to ensure that narrow spikes in the TDC are resolved, as described in ref 49.  In short, the overlap between the adiabatic electronic wave functions associated with each electronic state at the beginning and end of each time step is tracked.  When this overlap falls below a user-defined threshold (generally chosen to be ) for any state, the TDC is necessarily large, and thus the time step is reduced to ensure accurate integration.  This adaptation procedure is not required when NPI is used, reducing the computational cost of NPI calculations significantly. 40  The spawning procedure has also been modified.  In a typical AIMS simulation based on analytical calculations of the NACME, the spawning procedure proceeds as follows.  For each TBF, the NACME vector between the present state and all other states is computed at every time step.  When the norm of this vector rises above a user-specified threshold (the spawning threshold), a spawning event is triggered.  The simulation is stopped, and this initial TBF (hereafter known as the parent) is tentatively propagated forward in time until the norm of the NACME reaches a local maximum in time.  At this point, a new (child) basis function is spawned.  The basis function is created with the same average position, , as the parent.  In order to ensure long-time energy conservation, the momentum of the child is scaled in the NACME vector direction such that the parent and child have the same total energy.  This procedure can be justified by semiclassical arguments.57  The child TBF is then propagated backward to the point in time at which the spawning procedure was initiated and added to the nuclear basis set, and the simulation is continued.  For more detailed discussions of the spawning procedure, see references 58 and 59.  Figure 3.2 AIMS-NPI test systems  41  The fact that the NACME vector and its norm are never computed explicitly in the NPI approximation requires us to eliminate the dependence of the spawning procedure on these quantities.  Thus, when employing NPI, rather than tracking the norm of the NACME as the quantity which triggers the spawning event and determines the position of the child basis function, the absolute value of the TDC itself is used.  Similarly, with NPI it is no longer possible to adjust the momentum of the child TBF in the NACME direction, so the momentum itself is used.  Specifically, the momentum of the child is adjusted according to   (3.4) where , , and  are the average momentum, classical kinetic energy, and classical potential energy of TBF I, respectively (where classical indicates that the kinetic and potential energies are computed at  , rather than as the expectation value over the entire TBF.)  Alternative procedures for momentum adjustment that do not require the NACME vector include scaling in the direction of the energy difference gradient vector 60 or an optimal spawning algorithm.59  So as to most directly compare the NPI and analytical schemes, in the present work all calculations use this modified momentum scaling procedure, even those that are based on analytical calculations of the NACME.   It is noteworthy that the two most common approaches to scaling the momentum in the absence of the NACME each have a significant drawback.  Scaling in the momentum direction is not size-consistent, though the momentum necessarily overlaps strongly with the preferred NACME direction when the TDC is large.  On the other hand, scaling in the gradient difference direction is size-consistent, but in the linear region near the intersection the difference gradient is orthogonal to the NACME.  Though we take a different approach in the present work, when NPI 42  is used in conjunction with an electronic structure method for which analytical NACMEs are available, it is likely advantageous to compute the full NACME during spawns in order to perform the momentum scaling procedure.  This approach would achieve the best of both approaches: the accurate transition probabilities of NPI and correct momentum scaling during nonadiabatic events.  3.3 Methods  The AIMS-NPI approach will be applied to two test systems: one which exhibits a typical (N-2)-dimensional CI and one which exhibits an effectively (N-1)-dimensional TUC.  These systems are illustrated in Figure 3.2.  As an example of a system exhibiting a CI, we will model the non-radiative decay of propene, which relaxes via a pyramidalized CI similar to well-known intersections in other small unsaturated molecules such as ethylene.25  Three different schemes for computing the TDC have been tested: NPI, analytical NACMEs with adaptive integration (the r), and analytical this work adaptive integration is never used in conjunction with NPI.  A time step of 0.48 fs (20.0 a.u.) was employed to integrate for a total simulation time of 1.0 ps.  Within the analytical scheme, the time step is allowed to adapt down to 0.0024 fs (0.1 a.u.).  For each TDC scheme, we have run 50 simulations with different initial conditions (positions and momenta) sampled from the ground state vibrational Wigner distribution and placed on the lowest singlet excited (S1) electronic state.  The same set of 50 initial conditions was used for each scheme.  In all simulations of propene, we have computed the electronic wave function and PES at the state averaged (SA-) CASSCF level of theory, using an active space of two electrons in two orbitals, averaging over three electronic states (SA-3-CAS(2,2)), and using the 6-31G* basis set. 43   In order to assess the accuracy of the NPI scheme in a system exhibiting a TUC, we have run AIMS simulations of a dimer composed of a trifluoromethyl (CF3) radical and a trifluoromethyl anion (CF3-) separated by a distance of 1000 Å.  Here we consider only the NPI scheme and the analytical scheme.  In this system time steps of 0.12, 0.24, and 0.48 fs have been employed with the NPI scheme to investigate the scaling of error with time step.  With the analytical scheme an initial time step of 0.48 fs was used, and the time step is allowed to adapt down to 0.00024 fs.  All simulations were run for 1.0 ps.  As with propene we have averaged the results of 50 simulations, with initial conditions sampled from the ground state vibrational Wigner distribution.  Again the same set of 50 initial conditions was used for each scheme, and all population was initially placed in the S1 state.  We compute the PES at the SA-CASSCF level with three active electrons in two active orbitals and average over 2 electronic states (SA-2-CAS(3,2)).  The 6-31G* basis set is used.   The two active orbitals for this dimer are the nonbonding valence orbitals of each CF3 monomer. Thus, the two electronic states which arise from the SA-2-CAS(3/2) calculation differ by a long range charge transfer excitation, moving the excess electron from one CF3 monomer to the other.  Therefore, upon excitation the S1 state corresponds to a charge transfer state, and nuclear relaxation leads to a crossing between the S1 and S0 states.  Due to the large distance between monomers, the coupling between diabats is essentially zero.  Thus, this crossing is a TUC, and diabatic passage is expected.  Inaccurate integration of the TDSE, however, can lead to unphysical long range charge transfer, which is the result of adiabatic passage through the TUC (Figure 3.1, arrow 4).  Such unphysical behavior can be easily tracked by following the expectation value of the total Mulliken charge of a single monomer as a function of time.  As such, we will track this value as computed in the zeroth order saddle point approximation. 44   All calculations were performed with the FMSMolPro software package,29 which interfaces AIMS dynamics with the MolPro electronic structure software package.61-63  Further discussion of simulation parameters can be found in the appendix.   3.4 AIMS-NPI near a conical intersection  For propene, which exhibits a CI as noted above, we compare three schemes for computing the TDC: NPI, analytical NACMEs, and analytical NACMEs without adaptive integration.  Figure 3.3 shows the S0 population as a function of time averaged over 50 AIMS simulations using each of these three schemes.  The results from the 50 individual simulations are plotted in Figures A3.1-A3.3.  Agreement between NPI and the analytical scheme is excellent.  The difference between the two population curves is on average 0.6%, and the maximum difference at a particular time is 2.7%.  Comparing the analytical scheme with and without adaptive integration, the time-averaged error is 1.6%, with a maximum error of 10.5%.  Thus, NPI does not require adaptive integration to achieve accuracy comparable to that of the adaptively integrated analytical scheme, and thus is able to provide high accuracy at a low computational cost.  45   Figure 3.3 AIMS-NPI population transfer in propene The S0 population averaged over 50 AIMS simulations of propene is plotted as a function of time.  Results using NPI and the analytical scheme (with and without adaptive integration) for evaluation of the time-derivative coupling are shown in red, blue, and green, respectively.   Interestingly, the analytical simulations spawn slightly less often than the NPI simulations.  On average, each initial basis function in the analytical simulations spawns 4.4 child trajectories, while the initial basis functions in the NPI simulations each spawn 5.8 child trajectories on average.  This discrepancy, which likely arises because NPI better resolves spikes in the TDC, is not reflected in the computed population transfer probabilities.  However, as will be discussed below this difference can become extreme in the limit of a TUC.    46  3.5 AIMS-NPI at a trivial unavoided crossing Here we investigate the ability of AIMS-NPI to describe diabatic passage through a TUC, which we will characterize by following changes in the total charge of a single CF3 monomer (termed monomer 1) as a function of simulation time.  At the 50 randomly sampled initial positions used in our simulations, monomer 1 has an average charge of -1.0000 in the populated S1 state and a charge of 0.0000 in S0.  A large deviation in the charge from -1.0 would indicate unphysical adiabatic passage through the TUC between these states, and thus inaccurate integration of the TDSE.   The average charges as a function of time computed via AIMS-NPI and AIMS with analytical NACMEs are presented in Figure 3.4.  The NPI results were computed with a fixed time step of 0.12 fs, while adaptive integration (with time steps noted above) was used for the analytical simulations.      Figure 3.4 AIMS-NPI anionic monomer charge in a long range separated dimer The incoherently averaged total Mulliken charge of a single CF3 monomer for 50 simulations of the trifluoromethyl dimer is plotted as a function of time after excitation to S1.  Results computed with the analytical and NPI schemes are shown in blue and red, respectively.  47   Figure 3.5 AIMS-NPI population transfer in a long range separated dimer The S0 population averaged over 50 AIMS simulations of the trifluoromethyl dimer is plotted as a function of time.  Results using NPI and the analytical scheme to compute the TDC are shown in red and blue, respectively.  The NPI approach, which does not rely upon the value of the NACME computed at discrete times, predicts physically correct diabatic population transfer whenever a TUC is encountered, resulting in the fact that the charge remains essentially constant for the entire duration of the simulation.  The average charge on monomer 1 is -0.9993, indicating that very little unphysical long range charge transfer occurs.  In stark contrast, the analytical scheme results in dramatic oscillations in the charge, indicating frequent unphysical long range change transfer events due to inaccurate integration of the TDSE at TUCs.  The oscillations arise because the spike in the TDC is never resolved, and thus zero population decay to the ground state is observed.  Instead, the excited state population passes adiabatically back and forth between the two oppositely polarized charge transfer states (Figure 3.1, arrow 4), resulting in rapid and unphysical long range charge transfer events.  This assessment is confirmed by investigation of the adiabatic state populations 48  (Figure 3.5).  No decay to S0 is observed in the calculations based on analytical NACMEs.  In fact, the NACME-based simulations do not even spawn, because the spike in the NACME is too narrow to be resolved. When NPI is used with a time step greater than 0.12 fs the propensity for unphysical charge transfer increases.  Figure 3.6 shows the average charge on monomer 1 for 50 simulations of the trifluoromethyl dimer computed using the NPI scheme with different time steps.  As above, any significant change in the charge indicates unphysical charge transfer.  With a time step of 0.24 fs the average charge deviates from -1.0 by 0.008 (compared to 0.0007 for a 0.12 fs step).  When the time step is increased to 0.48 fs the average error increases to 0.075. (Plots of the results obtained for the 50 individual simulations with these three time steps are presented in Figures A3.4-A3.6).  These results suggest that choosing a time step larger than ~0.24 fs may be inadvisable in the present AIMS implementation of the NPI scheme in systems where TUCs are likely to be encountered.  In future work it would likely be advantageous to devise an adaptive integration scheme for use with NPI.    49   Figure 3.6 AIMS-NPI long range charge transfer error and time step size The average charge of a single CF3 monomer is plotted as a function of the simulation time for AIMS-NPI using various time steps.    3.6 Discussion In the NPI approximation the adiabatic wave functions are assumed to rotate smoothly and linearly over the duration of each time step.  In the case of a TUC or a direct hit on a CI, however, the true duration of the spike in the TDC is much shorter than the time step.  The result of this is that NPI overestimates dephasing effects, which translates to errors in the population transfer probability.  In our previous study of the accuracy of NPI for the numerical integration of the adiabatized Landau-Zener Hamiltonian, errors of 0.7% were attributed to dephasing when a time step of 0.48 fs was employed. 55  Thus, it is surprising that 7.5% of the population undergoes long range charge transfer in the trifluoromethyl dimer when a time step of 0.48 fs is used in the present simulations, and it seems unlikely that this stems entirely from the overestimation of dephasing effects.  We attribute this error, which increases with the size of the time step, to several sources, which we discuss in this subsection. 50  Population transfer between two trajectories in an AIMS simulation is governed, through the TDSE, by the off-diagonal elements of the Hamiltonian matrix .  (3.5) The exact evaluation of this integral would require knowledge of the electronic structure at all possible nuclear geometries, which is obviously unfeasible.  Thus in AIMS it is standard to compute the coupling in the zeroth order saddle-point approximation (SPA) .  (3.6) where  is the vector of nuclear coordinates at the maximum of the nuclear overlap density, .  An integral computed in the SPA requires only a single electronic structure calculation at .   As we have noted in our past study, when  a system will exhibit perfect diabatic passage of population from one adiabatic state to the other in a single time step in the absence of significant dephasing or decoherence effects.  It is a notable feature of the NPI approach that this is the value for the TDC given by Eq. (2.10) when two adiabatic states experience a complete exchange in character in a single time step, as at a TUC.  Thus, in our study we can assess the source of errors in the computed population transfer by comparing the computed value of  at a TUC to .  In our simulations of the trifluoromethyl dimer, the value of  upon crossing a TUC is never in error by more than 0.001%, but unphysical charge transfer can still occur when the nuclear overlap, , deviates from unity at the point of maximum coupling.   51  Such deviations can occur as a result of the fact that new TBFs are spawned only at discrete simulation time steps.  As discussed previously, when a child TBF is spawned the momentum is scaled such that the classical energy of the parent and child trajectories are identical.  In the ideal case that a spawning event takes place at the exact intersection geometry (where the energy gap is precisely zero) no adjustment of the momentum is required.  This leads to perfect nuclear overlap.  In practice, however, the point of zero energy gap is unlikely to fall exactly on a time step, and thus even at a TUC, some momentum adjustment occurs upon spawning, and in this case the nuclear overlap term in Eq. (3.6) deviates from unity.  This deviation will, in turn, result in an incorrect prediction of the population transfer probability. Analysis of the TBF overlaps in the trifluoromethyl dimer at the points of maximum coupling suggests that non-unity TBF overlap becomes a significant source of error as the time step is increased.  When a time step of 0.48 fs is used TBF overlaps as small as 0.76 are observed at the points of maximum overlap between parent and child, and the average deviation from unity TBF overlap at these points is 0.056.  The attribution of the error observed in Figure 3.6 to deviations from unity overlap is also supported by the observed scaling of the error with time step.  Assuming that spawning occurs in the linear region around a CI or TUC, one would expect the average momentum difference between parent and child to scale linearly with time step.  The deviation in the nuclear overlap in Eq. (3.6) from unity will scale quadratically with the momentum difference, thus the error in the coupling itself, , will scale quadratically with time step.  The deviation of the population transfer probability from unity will scale quadratically with the error in the coupling, therefore the overall error in population transfer probability will scale as , which is roughly what is observed in Figure 3.6.  It is noteworthy that within the SPA the overlap integrals between trajectories are computed exactly, and thus the error discussed here is not 52  attributable to the SPA itself, but would remain even if Eq. (3.5) was evaluated exactly.  Instead, the error is attributable to the choice of a finite nuclear basis and the discretization of time. Though minimal in the cases presented here, one additional source of error was observed in our work with the AIMS-NPI scheme near TUCs.  The AIMS TBFs are Born-Oppenheimer states, and thus can couple via the Born-Oppenheimer electronic Hamiltonian to other basis functions on the same adiabatic state.  When two TBFs come near each other, population will transfer between them, even if their centroids fall on opposite sides of a TUC.  Thus, it is possible to see unphysical long range charge transfer via an adiabatic mechanism.  This is an artifact of our finite nuclear basis and the choice of the adiabatic representation for our expansion of the AIMS wave function (Eq. (3.1)), and may lead to significant errors in systems which cross TUCs more frequently than the systems studied here.  In the complete basis limit one would expect this error to disappear.   Note that AIMS is exact quantum dynamics in the limit of 1) a complete nuclear basis, 2) exact electronic structure, 3) an infinitesimal time step, and 4) exact nuclear integral evaluation.  It is worth considering whether this property is maintained in the NPI approximation.  On one hand, in the limit of infinitesimal time step NPI reduces to the exact TDC.  However, there is no obvious recipe to apply NPI outside the SPA.  Thus NPI is not consistent with exact integral evaluation and is therefore an additional approximation which must be eliminated to achieve numerical exactness. It is also interesting to consider whether NPI can be of use in Ehrenfest and Ehrenfest-like dynamical simulations.13,30-32,64  Unlike spawning and surface hopping methods, the full NACME vector is required by the Ehrenfest method, because the NACME vector corresponds to the off-diagonal component of the force.  As such, we see no easy way to replace the calculation of the 53  NACME by the NPI approximation in this context.  However, in the case of very narrowly avoided crossings, NPI allows a more accurate prediction of electronic state transition probabilities than does the NACME, and thus a hybrid approach, utilizing NPI to solve the TDSE and analytical NACME vectors to compute the force, would capture the best features of both methods with no significant increase in computational cost.  This is exactly analogous to the hybrid approach for momentum scaling proposed in Section 3.2.  3.7 Time averaging of the TDC Here we investigate the role that time averaging of the TDC plays in the improved accuracy of NPI compared with the analytical approach.  We do so by numerical integration of the adiabatized Landau-Zener (LZ) model,65,66 building off of our previous work on NPI.55  In this previous study, we demonstrated that integration of the NPI approximation to the TDC provides a more accurate prediction of the population transfer probability (compared to the exact LZ solution) than does the integration of the analytically computed TDC at a set of finite time steps.  Here, in order to elucidate the fact that this is accomplished by averaging the TDC over the time step, and thus eliminating error due to discretization of time, we apply a new numerical scheme for integrating the TDSE of the adiabatized LZ model.  Specifically, we approximate the TDC at the middle of each time step as the average TDC over that time step . (3.7)  Note that in Eq. (3.7) we average over the analytical TDC computed from the exact LZ adiabatic wave functions (which are known at all times in the analytical LZ model), rather than over TDC computed from the approximated NPI wave functions (as in Eq. (2.10)).  Thus, this new 54  approximation allows us to consider separately the effect of interpolation and that of time averaging.  We refer to the scheme in which the TDC is computed according to Eq. (3.7) as  Our methodology is reported in detail in ref. 55, but here we briefly summarize.  The LZ model is adiabatized and numerically integrated using the fourth order Runge-Kutta method and a multiple time step algorithm similar to that used in AIMS.29  The error in population transfer probability relative to the exact LZ solution is computed as a function of time.  In the present work, the coupling between diabats, , and the slope of the diabats, , in the LZ model are set to 10 meV and 0.1 eV/fs, respectively.  These values correspond to a realistic but very weakly avoided crossing.  The error in the population transfer probability computed using the averaged analytical scheme as a function of time step is presented in Figure 3.7, along with those computed utilizing the analytical TDC without time averaging and the time-averaged NPI TDC.  The latter two sets of data are reproduced from ref. 55 for comparison.  Note that the error arising from the averaged analytical approach is much smaller than that computed without time averaging.  In fact, the time-averaged analytical approach results in very small errors comparable to those produced with NPI, suggesting that the success of NPI arises from the elimination of error due to discretization of time.  55   Figure 3.7 NPI and time-averaging of the time-derivative coupling The error in the numerically computed LZ transition probability as a function of time step.  Errors in simulations based on the averaged analytical approach described in the text are shown in green.  Errors from simulations based on NPI (blue) and analytical TDCs computed without time averaging (red) are reproduced from ref. 55 for comparison.   Note that in real AIMS simulations based on ab initio PESs the averaged analytical scheme cannot be practically implemented because it would require knowledge of the electronic wave function at all times.  Thus, interpolation is necessary to allow time averaging, and it is this time averaging which eliminates errors in population transfer probabilities due to discretization of time.    56  3.8 Conclusions In this work, we analyzed the behavior of the NPI scheme for computing time-derivative couplings used in conjunction with the AIMS dynamical method.  NPI offers a significant advantage over the more widely used approach of computing the TDC from the nonadiabatic coupling vector.  In a sense, in the analytical approach, the TDC is known exactly, but only at discrete points in time, while in NPI, the TDC is computed approximately as a continuous function of time.  By doing so, NPI provides more accurate predictions of population transfer probabilities than the analytical scheme when the TDC varies rapidly without the need to resort to adaptive integration.  In our studies, a time step of 0.12 fs was sufficient to accurately predict population transfer probabilities in the limit of a TUC, where the analytical TDC exhibits an infinitesimally narrow spike and thus is not numerically integrable at any finite time step.  Errors in the population transfer probability at time steps greater than 0.12 fs are attributable to the use of a finite time step and nuclear basis in the spawning procedure itself, and point to directions for future development.     57             APPENDIX   58   The initial conditions of our simulations were randomly sampled form the ground state vibrational Wigner distribution computed in the harmonic approximation at the B3LYP/6-31G** level of theory.  In the case of NPI and the analytical scheme, respectively, the spawning thresholds were chosen by identifying a baseline value for the TDC or norm of the NACME.  In the case of NPI, the spawning threshold was set to a value of 1.542x10-3 Hartree.  In all simulations employing the analytical scheme for evaluation of the time-derivative coupling, the spawning threshold was set to 3.0 Bohr-1.  When the TDC spikes above the spawning threshold, a child trajectory is spawned only if the population of the parent is above a threshold; this threshold was set to 0.1 for all propene simulations, and 0.05 for all trifluoromethyl dimer simulations.  Additionally, the spawning of trajectories is aborted if the overlap of the newly created trajectory with any other trajectory is greater than 0.6 in the case of propene, and 0.8 in the case of trifluoromethyl dimer.  Gaussian widths of 4.70, 22.7, and 8.50 Bohr-2 are associated with hydrogen, carbon, and fluorine degrees of freedom.     59   Figure A3.1 AIMS-NPI results for individual simulations of propene The ground state (S0) population of propene is plotted for 50 AIMS simulations (gray) employing the NPI approach for evaluation of the TDC.  The average is shown in red.   Figure A3.2 AIMS results for individual simulations of propene The ground state (S0) population of propene is plotted for 50 AIMS simulations (gray) employing the analytical approach for evaluation of the TDC.  The average is shown in red.   60   Figure A3.3 AIMS results (without adaptive integration) for simulations of propene The ground state (S0) population of propene is plotted for 50 AIMS simulations (gray) employing the analytical approach for evaluation of the TDC without adaptive integration.  The average is shown in red.    Figure A3.4 AIMS-NPI results for trifluoromethyl dimer (0.12 fs time step) The charge of monomer 1 of the trifluoromethyl dimer as a function of time is plotted for 50 AIMS simulations (gray) employing the NPI scheme for evaluation of the TDC with a time step of 0.12 fs. The average of all 50 simulations is shown in red.  61   Figure A3.5 AIMS-NPI results for trifluoromethyl dimer (0.24 fs time step) The charge of monomer 1 of the trifluoromethyl dimer as a function of time is plotted for 50 AIMS simulations (gray) employing the NPI scheme for evaluation of the TDC with a time step of 0.24 fs. The average of all 50 simulations is shown in red.  Figure A3.6 AIMS-NPI results for trifluoromethyl dimer (0.48 fs time step) The charge of monomer 1 of the trifluoromethyl dimer as a function of time is plotted for 50 AIMS simulations (gray) employing the NPI scheme for evaluation of the TDC with a time step of 0.48 fs. The average of all 50 simulations is shown in red.  62             REFERENCES    63  REFERENCES   (1) Tully, J. C.; Preston, R. K. TRAJECTORY SURFACE HOPPING APPROACH TO NONADIABATIC MOLECULAR COLLISIONS - REACTION OF H+ WITH D2. Journal of Chemical Physics 1971, 55, 562-572.   (2) Micha, D. A. A SELF-CONSISTENT EIKONAL TREATMENT OF ELECTRONIC-TRANSITIONS IN MOLECULAR-COLLISIONS. Journal of Chemical Physics 1983, 78, 7138-7145.   (3) Tully, J. C. MOLECULAR-DYNAMICS WITH ELECTRONIC-TRANSITIONS. Journal of Chemical Physics 1990, 93, 1061-1071.   (4) Hammes-Schiffer, S.; Tully, J. C. PROTON-TRANSFER IN SOLUTION - MOLECULAR-DYNAMICS WITH QUANTUM TRANSITIONS. Journal of Chemical Physics 1994, 101, 4657-4667.   (5) Bittner, E. R.; Rossky, P. J. QUANTUM DECOHERENCE IN MIXED QUANTUM-CLASSICAL SYSTEMS - NONADIABATIC PROCESSES. Journal of Chemical Physics 1995, 103, 8130-8143.   (6) Coker, D. F.; Xiao, L. METHODS FOR MOLECULAR-DYNAMICS WITH NONADIABATIC TRANSITIONS. Journal of Chemical Physics 1995, 102, 496-510.   (7) Martens, C. C.; Fang, J. Y. Semiclassical-limit molecular dynamics on multiple electronic surfaces. Journal of Chemical Physics 1997, 106, 4918-4930.   (8) Muller, U.; Stock, G. Surface-hopping modeling of photoinduced relaxation dynamics on coupled potential-energy surfaces. Journal of Chemical Physics 1997, 107, 6230-6245.   (9) Prezhdo, O. V.; Rossky, P. J. Mean-field molecular dynamics with surface hopping. Journal of Chemical Physics 1997, 107, 825-834.   (10) Tully, J. C. Mixed quantum-classical dynamics. Faraday Discussions 1998, 110, 407-419.   (11) Kapral, R.; Ciccotti, G. Mixed quantum-classical dynamics. Journal of Chemical Physics 1999, 110, 8919-8929.   (12) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. Journal of Chemical Physics 2001, 114, 10608-10615.  64   (13) Doltsinis, N. L.; Marx, D. Nonadiabatic Car-Parrinello molecular dynamics. Phys. Rev. Lett. 2002, 88.   (14) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. Mean-field dynamics with stochastic decoherence (MF-SD): A new algorithm for nonadiabatic mixed quantum/classical molecular-dynamics simulations with nuclear-induced decoherence. Journal of Chemical Physics 2005, 123.   (15) Craig, C. F.; Duncan, W. R.; Prezhdo, O. V. Trajectory surface hopping in the time-dependent Kohn-Sham approach for electron-nuclear dynamics. Physical Review Letters 2005, 95.   (16) Li, X. S.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab initio Ehrenfest dynamics. Journal of Chemical Physics 2005, 123.   (17) Jasper, A. W.; Nangia, S.; Zhu, C. Y.; Truhlar, D. G. Non-Born-Oppenheimer molecular dynamics. Acc. Chem. Res. 2006, 39, 101-108.   (18) Fischer, S. A.; Chapman, C. T.; Li, X. S. Surface hopping with Ehrenfest excited potential. Journal of Chemical Physics 2011, 135.   (19) Shenvi, N.; Subotnik, J. E.; Yang, W. Simultaneous-trajectory surface hopping: A parameter-free algorithm for implementing decoherence in nonadiabatic dynamics. Journal of Chemical Physics 2011, 134.   (20) Richter, M.; Marquetand, P.; Gonzalez-Vazquez, J.; Sola, I.; Gonzalez, L. SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. Journal of Chemical Theory and Computation 2011, 7, 1253-1258.   (21) Huo, P. F.; Coker, D. F. Consistent schemes for non-adiabatic dynamics derived from partial linearized density matrix propagation. J. Chem. Phys. 2012, 137.   (22) White, A. J.; Gorshkov, V. N.; Wang, R. X.; Tretiak, S.; Mozyrsky, D. Semiclassical Monte Carlo: A first principles approach to non-adiabatic molecular dynamics. J. Chem. Phys. 2014, 141.   (23) Koppel, H.; Domcke, W.; Cederbaum, L. S. MULTIMODE MOLECULAR-DYNAMICS BEYOND THE BORN-OPPENHEIMER APPROXIMATION. Advances in Chemical Physics 1984, 57, 59-246.   (24) Burghardt, I.; Meyer, H. D.; Cederbaum, L. S. Approaches to the approximate treatment of complex molecular systems by the multiconfiguration time-dependent Hartree method. J. Chem. Phys. 1999, 111, 2927-2939.   (25) Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics. Journal of Physical Chemistry A 2000, 104, 5161-5175. 65    (26) Shalashilin, D. V.; Child, M. S. Time dependent quantum propagation in phase space. Journal of Chemical Physics 2000, 113, 10028-10036.  (27) Worth, G. A.; Robb, M. A.; Burghardt, I. A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets. Faraday Discussions 2004, 127, 307-323.   (28) Chen, X.; Batista, V. S. Matching-pursuit/split-operator-Fourier-transform simulations of excited-state nonadiabatic quantum dynamics in pyrazine. Journal of Chemical Physics 2006, 125.   (29) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of ab initio multiple spawning in the MOLPRO quantum chemistry package. Chemical Physics 2008, 347, 3-16.   (30) Worth, G. A.; Robb, M. A.; Lasorne, B. Solving the time-dependent Schrodinger equation for nuclear motion in one step: direct dynamics of non-adiabatic systems. Mol. Phys. 2008, 106, 2077-2091.   (31) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics. J. Chem. Phys. 2014, 141.   (32) Makhov, D. V.; Saita, K.; Martinez, T. J.; Shalashilin, D. V. Ab initio multiple cloning simulations of pyrrole photodissociation: TKER spectra and velocity map imaging. Phys. Chem. Chem. Phys. 2015, 17, 3316-3325.   (33) Lengsfield, B. H.; Saxe, P.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS USING SA-MCSCF/CL WAVE-FUNCTIONS AND ANALYTIC GRADIENT METHODS .1. Journal of Chemical Physics 1984, 81, 4549-4553.   (34) Lengsfield, B. H.; Yarkony, D. R. NONADIABATIC INTERACTIONS BETWEEN POTENTIAL-ENERGY SURFACES - THEORY AND APPLICATIONS. Advances in Chemical Physics 1992, 82, 1-71.   (35) Chernyak, V.; Mukamel, S. Density-matrix representation of nonadiabatic couplings in time-dependent density functional (TDDFT) theories. Journal of Chemical Physics 2000, 112, 3572-3579.   (36) Ichino, T.; Gauss, J.; Stanton, J. F. Quasidiabatic states described by coupled-cluster theory. Journal of Chemical Physics 2009, 130.   (37) Tavernelli, I.; Curchod, B. F. E.; Laktionov, A.; Rothlisberger, U. Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm-Dancoff approximation and beyond. Journal of Chemical Physics 2010, 133.  66   (38) Nelson, T.; Fernandez-Alberti, S.; Chernyak, V.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics Modeling of Photoinduced Dynamics in Conjugated Molecules. Journal of Physical Chemistry B 2011, 115, 5402-5414.   (39) Fatehi, S.; Alguire, E.; Shao, Y. H.; Subotnik, J. E. Analytic derivative couplings between configuration-interaction-singles states with built-in electron-translation factors for translational invariance. Journal of Chemical Physics 2011, 135.   (40) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martinez, T. J. Role of Rydberg States in the Photochemical Dynamics of Ethylene. Journal of Physical Chemistry A 2012, 116, 2808-2818.   (41) Fatehi, S.; Alguire, E.; Subotnik, J. E. Derivative couplings and analytic gradients for diabatic states, with an implementation for Boys-localized configuration-interaction singles. Journal of Chemical Physics 2013, 139.   (42) Yarkony, D. R. Diabolical conical intersections. Reviews of Modern Physics 1996, 68, 985-1013.   (43) Bernardi, F.; Olivucci, M.; Robb, M. A. Potential energy surface crossings in organic photochemistry. Chemical Society Reviews 1996, 25, 321-&.   (44) Worth, G. A.; Cederbaum, L. S. Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection. Annu. Rev. Phys. Chem. 2004, 55, 127-158.   (45) Levine, B. G.; Martinez, T. J. Isomerization through conical intersections. Annu. Rev. Phys. Chem. 2007, 58, 613-634.   (46) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annual Review of Physical Chemistry 2012, 63, 325-352.   (47) Neria, E.; Nitzan, A. SEMICLASSICAL EVALUATION OF NONADIABATIC RATES IN CONDENSED PHASES. Journal of Chemical Physics 1993, 99, 1109-1123.   (48) Wang, L. J.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals. Journal of Physical Chemistry Letters 2013, 4, 1888-1894.   (49) Virshup, A. M.; Levine, B. G.; Martinez, T. J. Steric and electrostatic effects on photoisomerization dynamics using QM/MM ab initio multiple spawning. Theoretical Chemistry Accounts 2014, 133, 1506.   (50) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. Journal of Chemical Physics 2012, 137. 67    (51) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. Chemical Physics Letters 2013, 590, 208-213.   (52) Wang, L. J.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. 2014, 5, 713-719.   (53) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials. Acc. Chem. Res. 2014, 47, 1155-1164.   (54) Pittner, J.; Lischka, H.; Barbatti, M. Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings. Chemical Physics 2009, 356, 147-152.   (55) Meek, G. A.; Levine, B. G. Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations. Journal of Physical Chemistry Letters 2014, 5, 2351-2356.   (56) Heller, E. J. TIME-DEPENDENT APPROACH TO SEMICLASSICAL DYNAMICS. Journal of Chemical Physics 1975, 62, 1544-1555.   (57) Herman, M. F. NONADIABATIC SEMICLASSICAL SCATTERING .1. ANALYSIS OF GENERALIZED SURFACE HOPPING PROCEDURES. J. Chem. Phys. 1984, 81, 754-763.   (58) Ben-Nun, M.; Martinez, T. J. Ab initio quantum molecular dynamics. Adv. Chem. Phys. 2002, 121, 439-512.   (59) Yang, S.; Martinez, T. J.: Nonclassical Phase Space Jumps and Optimal Spawning. In Advances in the Theory of Atomic and Molecular Systems: Dynamics, Spectroscopy, Clusters, and Nanostructures; Piecuch, P., Maruani, J., Delgado-Barrio, G., Wilson, S., Eds.; Progress in Theoretical Chemistry and Physics 20, 2009; pp 35-45.   (60) Blais, N. C.; Truhlar, D. G. TRAJECTORY-SURFACE-HOPPING STUDY OF NA(3P2P) +H2- NA(3S2S)+H2(V',J',THETA). J. Chem. Phys. 1983, 79, 1334-1342.   (61) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schuetz, M. Molpro: a general-purpose quantum chemistry program package. Wiley Interdisciplinary Reviews-Computational Molecular Science 2012, 2, 242-253.   (62) Knowles, P. J.; Werner, H. J. AN EFFICIENT 2ND-ORDER MC SCF METHOD FOR LONG CONFIGURATION EXPANSIONS. Chem. Phys. Lett. 1985, 115, 259-267.  68   (63) Werner, H. J.; Knowles, P. J. A 2ND ORDER MULTICONFIGURATION SCF PROCEDURE WITH OPTIMUM CONVERGENCE. J. Chem. Phys. 1985, 82, 5053-5063.   (64) Mendive-Tapia, D.; Lasorne, B.; Worth, G. A.; Bearpark, M. J.; Robb, M. A. Controlling the mechanism of fulvene S-1/S-0 decay: switching off the stepwise population transfer. Phys. Chem. Chem. Phys. 2010, 12, 15725-15733.   (65) Zener, C. Non-adiabatic crossing of energy levels. Proceedings of the Royal Society of London Series a-Containing Papers of a Mathematical and Physical Character 1932, 137.   (66) Landau, L. Zur Theorie der Energieubertragung.  II. Phys. Z. Sowjetunion 1932, 2, 46-51.    69  CHAPTER 4: WAVE FUNCTION CONTINUITY AND THE DIAGONAL BORN-OPPENHEIMER CORRECTION AT CONICAL INTERSECTIONS  This chapter is taken in its entirety from: Wave function continuity and the diagonal Born-Oppenheimer correction at conical intersectionsphysics, 2016, 144, 184109.  4.1 Introduction  Nonadiabatic molecular dynamics simulations are typically performed in either the adiabatic or diabatic representation of electronic states.  When the nuclear-electronic time-dependent Schrodinger equation is solved exactly, the results are independent of the choice of representation, but when approximate methods are used, the results can depend strongly on this choice.1,2  Advantages of the adiabatic representation include its uniqueness and the locality of the coupling between electronic states, which simplifies the accurate treatment of coherent nuclear processes.  Working in the adiabatic representation is not without its challenges, however.  The first-derivative nonadiabatic coupling between electronic states becomes singular at conical intersections, points in nuclear configuration space at which adiabatic electronic states become degenerate.3-8   This difficulty is often circumvented by working in the diabatic representation,9-24 though graceful handling of these singularities within the adiabatic representation can also achieve accurate predictions of population transfer probabilities.25-29   A second-derivative nonadiabatic term, the diagonal of which is referred to as the diagonal Born-Oppenheimer correction (DBOC), also arises.  Though methods to compute the DBOC have been widely applied in high-accuracy quantum chemical calculations,30-44 with rare exceptions45-48 the DBOC is neglected in nonadiabatic molecular dynamics simulations.  On its surface this approximation appears to be questionable, however, as the DBOC is singular at conical intersections.33  Intuitively, one might expect inclusion of the DBOC to significantly reduce the 70  rate of nonadiabatic transitions by repelling population away from the conical intersection seam, for example.  Recently it has been argued that neglect of the DBOC in surface hopping simulations is justified by the fact that, in practice, errors introduced by this approximation are cancelled by errors associated with neglecting geometric phase effects.47,49  Perhaps most troubling is that the energy singularities associated with the DBOC do not arise naturally in the diabatic representation.  Should these effects not be independent of representation? In this work, we demonstrate that the singularity in the DBOC at a conical intersection is an artifact of the representation of the wave function and therefore is not reflected in the exact quantum dynamics of a real physical system.  First we review the definition of the DBOC in the context of the Born-Oppenheimer (BO) approximation.  Consider the total molecular wave function in the BO approximation, , (4.1) where r and R represent the electronic and nuclear coordinates, respectively, t is time, and  and  are the respective nuclear and electronic wave functions associated with adiabatic electronic state, I to a wave function that describes both nuclear and electronic motion.)  The DBOC, defined , (4.2) arises from the operation of the nuclear kinetic energy operator on the R-dependent electronic wave function of a particular adiabatic electronic state.  Here  indexes nuclear degrees of freedom, M is the mass associated with that degree of freedom, and the subscript r indicates that integration is performed over electronic degrees of freedom only.  Atomic units are used throughout this work.  The DBOC is typically a small correction to the BO potential energy surface (PES) which is 71  computed only when highly accurate energetics are required, but as noted above this is not the case at the geometry of a conical intersection, RCI.  At the intersection point,  is discontinuous with respect to R, and this discontinuity results in a singularity in .  4.2 Proof of non-integrability of the DBOC at a CI in the adiabatic representation  We will now prove that, unlike the singularity in the first-derivative coupling, the singularity in the second-derivative term is not integrable.  Consider a two-dimensional system described in polar coordinates, R and , with a conical intersection between adiabatic states I and J at the origin, R=0.  For convenience we define a normalized branching angle coordinate, .  Thus  is in units of distance, the -domain is , and the volume element for integration is simply .   The well-known Berry phase condition states that the integral of the derivative coupling projected onto a closed path encircling the intersection will have a value of .50  A circle of radius R centered at the origin is an example of such a path, thus for any R  (4.3) where we arbitrarily choose positive , and  is the -component of the first-derivative coupling vector, . (4.4) By insertion of a resolution of the identity into the expression for the -component of the DBOC of state I, one can show that . (4.5) 72  (For conciseness we do not include fundamental constants and masses and we assume the electronic wave functions to be real.)  The average value of  on some closed circular path, , is bounded from below by the square of the average of  along that same path, . (4.6) Dividing both sides of Eq. (4.3) R, yields , (4.7) therefore  (4.8) and . (4.9) One can demonstrate that the -component of the first-derivative coupling vector is integrable at the intersection by integrating over some neighborhood of radius , . (4.10) By similarly integrating over the lower bound to the -term of the second-derivative nonadiabatic term one can demonstrate that the DBOC is not integrable on the same domain: . (4.11) 73  This last integral diverges, therefore the DBOC is not integrable over domains including the intersection.  4.3 Discussion  As proven in Section 4.2, the total DBOC contribution to the energy of the BO molecular wave function, , will diverge whenever  is non-zero at the intersection.  (For brevity, throughout this work we refer to population that is located at RCI and corresponds to one of the intersecting adiabadifferently, because of the discontinuity in  at RCI there can be no continuous molecular wave function, , for which  is non-zero at RCI.  Thus, as previously demonstrated by several authors, the probability of the molecule being at a conical intersection is zero in all well-behaved molecular wave functions in the BO approximation.51,52  (In fact, past proofs that the DBOC is integrable46,51 were based on the assumption of a BO wave function that decays rapidly to zero at the intersection, and thus do not contradict our proof above, which makes no such assumption.)  Now consider the exact molecular wave function represented as a sum over BO vibronic states:  . (4.12) We emphasize that in this work we use the term BO vibronic state to refer to the product of a nuclear wave function with an adiabatic electronic wave function.  If each BO vibronic state is taken to be a well-behaved wave function, then the total molecular wave function would have the same properties as  above, i.e. the probability of population being at a conical intersection is zero.  There is clearly no physical basis for this constraint, however, when we 74  consider other representations of the exact molecular wave function.  For example, consider a similar sum over vibronic wave functions in a diabatic electronic basis, with the diabatic states chosen to be smoothly varying at the intersection point.  Such states would not lead to analogous singularities in second-derivative nonadiabatic terms because their electronic wave functions do not have analogous discontinuities.  Thus, in the diabatic representation it is trivial to construct a well-behaved molecular wave function with density at the conical intersection.  The same lack of singularities at the intersection point can be seen in the exact decomposition of the time-dependent molecular wave function,53-55 . (4.13) It can be proven that any well-behaved molecular wave function can be factored such that both  and  are continuous functions of R (see Appendix A).  Thus  is certain to be continuous regardless of whether a conical intersection is present or not.  Therefore, as is also apparent in the diabatic representation of the exact wave function described above, the exact molecular wave function may have density at a conical intersection (i.e.  may be non-zero even when  projects onto the intersecting adiabatic electronic states).  We also note that there are no discontinuities in  to give rise to singularities in the second-derivative nonadiabatic energy of ,  , (4.14) which we now call the exact second-derivative energy correction (ESDEC) to distinguish it from the DBOC, which corresponds to a particular adiabatic electronic state.  (Notice that the ESDEC is time-dependent due to the time-dependence of .  An identical term arises in the 75  exact time-dependent PES described by Abedi, Maitra, and Gross.53,54)  The ESDEC provides a concise description of second-derivative nonadiabatic effects and thus is a useful point of reference for further discussion below.   Consider now the apparent inconsistency that arises from the facts established above: 1) any continuous molecular wave function in the BO approximation (i.e. restricted to a single adiabatic electronic state) has no density at a conical intersection, 2) the exact molecular wave function may be expanded as a sum over BO vibronic wave functions, and 3) the exact molecular wave function may have density at a conical intersection.  We reconcile these facts by recognizing that there is no requirement that the individual terms in Eq. (4.12) be well-behaved, only that the total molecular wave function be so.  Thus, for the three points above to hold true, the BO vibronic wave functions of the intersecting states in the summation in Eq. (4.12) must be discontinuous at a conical intersection point if the probability of being at that intersection is non-zero.  As such, whenever density is present at a conical intersection, the DBOC matrix elements, , diverge as described above.  This divergence is an inconvenient consequence of the discontinuity of the BO vibronic wave functions.  While discontinuity of the intersecting BO vibronic states at conical intersections is a necessary condition for a continuous molecular wave function, it is obviously not sufficient, so it is worthwhile to consider the constraints on  and  near conical intersections.  Because the adiabatic electronic wave functions of the intersecting states,  and , are discontinuous at RCI, the adiabatic nuclear wave functions,  and , must be constructed such that the following condition in satisfied:  is independent of the direction in which the limit is taken.  (The intersection may be approached from any direction 76  in the branching plane.)  In some special cases, this condition may be satisfied when  and  are continuous functions of R, but there is no a priori reason to expect them to be continuous, and in the most general case  and  will be discontinuous at RCI.  (See Appendix B for an illustrative example.)  Thus, the choice of a wave function ansatz that requires  and  to be continuous is highly restrictive, and most often continuity of the total molecular wave function will be achieved via compensating discontinuities in , , , and .    To summarize, we have demonstrated that the singularities in the DBOC at conical intersections arise due to the fact that the BO vibronic wave functions summed over in Eq. (4.12) necessarily exhibit discontinuities that are not present in the well-behaved total molecular wave function.  These singularities are not integrable, and therefore result in undefined DBOC matrix elements.  In contrast, we have demonstrated that the exact electronic wave function must be continuous, and therefore the ESDEC does not exhibit analogous singularities at conical intersections.  Thus, undefined DBOC matrix elements at conical intersections are artifacts that reflect discontinuities in the chosen basis functions, and do not indicate that the true probability of a molecule being at a conical intersection is zero.   In other words, difficulties arise because one is representing a differentiable function as a sum of two (or more) non-differentiable functions, and then attempting to compute the derivative of the total function as the sum of the derivatives of the non-differentiable ones.  We now discuss the implications of the above arguments for nonadiabatic molecular dynamics simulations.  For all of the reasons described above, inclusion of the DBOC into fully quantum mechanical nonadiabatic molecular dynamics simulations is discouraged.  However, our 77  arguments say nothing about the expected magnitude of , so complete neglect of second-derivative nonadiabatic terms may still be ill-advised.  The exact R-dependent electronic wave function, , will depend on the specific dynamics of the system being studied, therefore  is dependent on those dynamics as well.  One can imagine that complete neglect of the ESDEC may be a good approximation in some cases but not in others.  More work is needed to investigate the magnitude of the ESDEC in realistic systems.  In the present communication we only argue that the singularities in the DBOC at conical intersections are not physically meaningful.  Similar arguments suggest that inclusion of the DBOC in mixed quantum-classical simulations of dynamics at conical intersections is not advisable, as has been demonstrated above for fully quantum simulations.  In the classical case, the singularity in the DBOC arises only when taking the semiclassical limit of a discontinuousand thus unphysicalmolecular wave function.  Additionally, inclusion of the DBOC in mixed quantum-classical simulations ensures that no population reaches the conical intersection seam itself, though such a constraint has no physical basis.  Thus, as with the fully quantum case, inclusion of the DBOC in mixed quantum-classical simulations is discouraged.  Because  cannot be easily determined from a single classical trajectory, there is no immediately apparent recipe to incorporate the ESDEC into widely used methods like surface hopping.  However, Gross and coworkers have recently reported a promising mixed quantum-classical approach which is based on the exact factorization of the molecular wave function and naturally includes the ESDEC.48   78  4.4 Conclusions  The above arguments suggest new avenues for research.  We hope that they will stimulate further investigations into the magnitude of the ESDEC in real systems and the role that it may or may not play in nonadiabatic dynamics near conical intersections.  Additionally, though the above arguments suggest that the adiabatic representation is highly problematic when second-derivative nonadiabatic matrix elements are incorporated, its uniqueness and the locality of the associated interstate couplings still strongly recommend it.  Further work may suggest means to achieve a best-of-both-worlds approach, which maintains these advantages to a great degree while circumventing its difficulties.     79             APPENDICES   80  APPENDIX A: Proof that a Continuous Molecular Wave Function can be Decomposed into Continuous Nuclear and Electronic Factors  Abedi, Maitra, and Gross have shown that the exact molecular wave functions can be factored according to56 . (4.15) As a well-behaved function,  must be finite-valued, continuous, and differentiable with respect to R.  Thus, the derivative of  with respect to R, , (4.16) must be finite-valued.  The factors,  and , must be finite-valued at all positions in order for the molecular wave function, , to be finite-valued.  Clearly  must be finite-valued in order for the nuclear kinetic energy to remain finite, thus  is continuous.  With , , and  all being finite-valued,  must also be finite-valued in order for  to be finite-valued.  Thus  is a continuous function of R.   81  APPENDIX B: Illustrative Example of a Continuous Molecular Wave Function which can be Represented as a Sum of Discontinuous Born-Oppenheimer Vibronic States  Consider a molecule with one nuclear degree of freedom near an unavoided crossing, which occurs at R=0.  The adiabatic electronic wave functions are defined  (4.17)  (4.18) where  and  are diabatic electronic wave functions that change continuously with R.  We define these states such that state 1 is higher in energy than state 0, except at R=0, where they are degenerate.  State 1 thus connects diabatically to state 0 at the intersection point and vice versa, and the adiabatic electronic wave functions change discontinuously at R=0.  Now consider the well-behaved total nuclear wave function, .  We can split this wave function into a term that is only non-zero when  and one which is only non-zero on the remainder of the domain, , (4.19) where  (4.20)  (4.21) If  is non-zero at the origin, then both  and  are discontinuous there.  We now define the total molecular wave function to be . 82   (4.22) Noting that  and  are orthogonal but not normalized, normalization requires that . (4.23) Continuity requires that the limits of the total molecular wave function as it approaches R=0 from the left or right are equivalent, . (4.24) Applying Eq. (4.24) to the definition of the total molecular wave function in Eq. (4.22) yields  (4.25) and  . (4.26) Therefore, there are three constraints on the four coefficients defining the total molecular wave function, so assigning a value to any one coefficient defines the rest.  For example, , , , and  satisfy all three constraints.  This wave function can be represented as a sum over BO vibronic states as  (4.27) with  (4.28) and . (4.29) These BO nuclear wave functions are discontinuous, as are the BO electronic wave functions defined in Eqs. (4.17) and (4.18).  The product BO vibronic functions are similarly discontinuous.  83  However, the total molecular wave function is continuous, as are the exact nuclear and electronic wave functions determined by factorization of this total wave function,  (4.30) where  (4.31) and . (4.32) Thus, a well behaved total molecular wave function can arise from summing over BO vibronic functions that are not, themselves, continuous.  Requiring the BO nuclear wave functions to be continuous introduces two additional constraints  (4.33) and . (4.34) Though a continuum of sets of coefficients meet the conditions for a well-behaved molecular wave function (Eqs. (4.23), (4.25), and (4.26)), only one of those also meets the conditions for continuity of the BO nuclear wave functions: , , , and , where  is a real angle that defines an arbitrary complex phase.  Representation of the 0 case as a sum over BO vibronic states (Eq. (4.27)) yields , (4.35) while exact factorization yields  (4.36) 84  and . (4.37) Thus, for this example, the space of possible continuous total molecular wave functions is strongly restricted when the continuity of the BO nuclear wave functions at a conical intersection is imposed as a constraint.    85             REFERENCES    86  REFERENCES   (1) Neria, E.; Nitzan, A. SEMICLASSICAL EVALUATION OF NONADIABATIC RATES IN CONDENSED PHASES. Journal of Chemical Physics 1993, 99, 1109-1123.   (2) Prezhdo, O. V.; Rossky, P. J. Mean-field molecular dynamics with surface hopping. Journal of Chemical Physics 1997, 107, 825-834.   (3) Yarkony, D. R. Diabolical conical intersections. Reviews of Modern Physics 1996, 68, 985-1013.   (4) Bernardi, F.; Olivucci, M.; Robb, M. A. Potential energy surface crossings in organic photochemistry. Chemical Society Reviews 1996, 25, 321-&.   (5) Yarkony, D. R. Conical intersections: The new conventional wisdom. J. Phys. Chem. A 2001, 105, 6277-6293.   (6) Worth, G. A.; Cederbaum, L. S. Beyond Born-Oppenheimer: Molecular dynamics through a conical intersection. Annu. Rev. Phys. Chem. 2004, 55, 127-158.   (7) Levine, B. G.; Martinez, T. J. Isomerization through conical intersections. Annu. Rev. Phys. Chem. 2007, 58, 613-634.   (8) Domcke, W.; Yarkony, D. R. Role of Conical Intersections in Molecular Spectroscopy and Photoinduced Chemical Dynamics. Annual Review of Physical Chemistry 2012, 63, 325-352.   (9) Baer, M. ADIABATIC AND DIABATIC REPRESENTATIONS FOR ATOM-DIATOM COLLISIONS - TREATMENT OF 3-DIMENSIONAL CASE. Chem. Phys. 1976, 15, 49-57.   (10) Macias, A.; Riera, A. CALCULATION OF DIABATIC STATES FROM MOLECULAR-PROPERTIES. Journal of Physics B-Atomic Molecular and Optical Physics 1978, 11, L489-L492.   (11) Mead, C. A.; Truhlar, D. G. CONDITIONS FOR THE DEFINITION OF A STRICTLY DIABATIC ELECTRONIC BASIS FOR MOLECULAR-SYSTEMS. Journal of Chemical Physics 1982, 77, 6090-6098.   (12) Pacher, T.; Cederbaum, L. S.; Koppel, H. APPROXIMATELY DIABATIC STATES FROM BLOCK DIAGONALIZATION OF THE ELECTRONIC HAMILTONIAN. J. Chem. Phys. 1988, 89, 7367-7381.   (13) Cave, R. J.; Newton, M. D. Generalization of the Mulliken-Hush treatment for the calculation of electron transfer matrix elements. Chem. Phys. Lett. 1996, 249, 15-19. 87   (14) Atchity, G. J.; Ruedenberg, K. Determination of diabatic states through enforcement of configurational uniformity. Theor. Chem. Acc. 1997, 97, 47-58.   (15) Dobbyn, A. J.; Knowles, P. J. A comparative study of methods for describing non-adiabatic coupling: diabatic representation of the (1)Sigma(+)/(1)Pi HOH and HHO conical intersections. Mol. Phys. 1997, 91, 1107-1123.   (16) Yarkony, D. R. On the adiabatic to diabatic states transformation near intersections of conical intersections. Journal of Chemical Physics 2000, 112, 2111-2120.   (17) Koppel, H.; Gronki, J.; Mahapatra, S. Construction scheme for regularized diabatic states. J. Chem. Phys. 2001, 115, 2377-2388.   (18) Nakamura, H.; Truhlar, D. G. The direct calculation of diabatic states based on configurational uniformity. J. Chem. Phys. 2001, 115, 10353-10372.   (19) Nakamura, H.; Truhlar, D. G. Direct diabatization of electronic states by the fourfold way. II. Dynamical correlation and rearrangement processes. J. Chem. Phys. 2002, 117, 5576-5593.   (20) Mota, V. C.; Varandas, A. J. C. HN2((2)A ') electronic manifold. II. Ab initio based double-sheeted DMBE potential energy surface via a global diabatization angle. J. Phys. Chem. A 2008, 112, 3768-3786.   (21) Zhu, X.; Yarkony, D. R. Toward eliminating the electronic structure bottleneck in nonadiabatic dynamics on the fly: An algorithm to fit nonlocal, quasidiabatic, coupled electronic state Hamiltonians based on ab initio electronic structure data. J. Chem. Phys. 2010, 132.   (22) Evenhuis, C.; Martinez, T. J. A scheme to interpolate potential energy surfaces and derivative coupling vectors without performing a global diabatization. J. Chem. Phys. 2011, 135.   (23) Sirjoosingh, A.; Hammes-Schiffer, S. Proton-Coupled Electron Transfer versus Hydrogen Atom Transfer: Generation of Charge-Localized Diabatic States. J. Phys. Chem. A 2011, 115, 2367-2377.   (24) Subotnik, J. E.; Alguire, E. C.; Ou, Q.; Landry, B. R.; Fatehi, S. The Requisite Electronic Structure Theory To Describe Photoexcited Nonadiabatic Dynamics: Nonadiabatic Derivative Couplings and Diabatic Electronic Couplings. Acc, Chem. Res. 2015, 48, 1340-1350.   (25) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. J. Chem. Phys. 2012, 137.   (26) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. Chemical Physics Letters 2013, 590, 208-213. 88   (27) Wang, L. J.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. Lett. 2014, 5, 713-719.   (28) Meek, G. A.; Levine, B. G. Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations. Journal of Physical Chemistry Letters 2014, 5, 2351-2356.   (29) Meek, G. A.; Levine, B. G. Accurate and efficient evaluation of transition probabilities at unavoided crossings in ab initio multiple spawning. Chem. Phys. 2015, 460, 117-124.   (30) Garrett, B. C.; Truhlar, D. G. NUCLEAR-MOTION CORRECTIONS TO BORN-OPPENHEIMER BARRIER HEIGHTS FOR CHEMICAL-REACTIONS. J. Chem. Phys. 1985, 82, 4543-4547.   (31) Lengsfield, B. H.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS FOR MCSCF CL WAVE-FUNCTIONS USING ANALYTIC DERIVATIVE METHODS .3. 2ND DERIVATIVE TERMS. Journal of Chemical Physics 1986, 84, 348-353.   (32) Handy, N. C.; Yamaguchi, Y.; Schaefer, H. F. THE DIAGONAL CORRECTION TO THE BORN-OPPENHEIMER APPROXIMATION - ITS EFFECT ON THE SINGLET-TRIPLET SPLITTING OF CH2 AND OTHER MOLECULAR EFFECTS. J. Chem. Phys. 1986, 84, 4481-4484.   (33) Saxe, P.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS FOR MCSCF/CL WAVE-FUNCTIONS .4. 2ND DERIVATIVE TERMS USING ANALYTIC GRADIENT METHODS. Journal of Chemical Physics 1987, 86, 321-328.   (34) Piecuch, P.; Li, X. Z.; Paldus, J. AN AB-INITIO DETERMINATION OF (1)A(1)-B-3(1) ENERGY-GAP IN CH2 USING ORTHOGONALLY SPIN-ADAPTED STATE-UNIVERSAL AND STATE-SPECIFIC COUPLED-CLUSTER METHODS. Chem. Phys. Lett. 1994, 230, 377-386.   (35) Wolniewicz, L.; Dressler, K. ADIABATIC POTENTIAL CURVES AND NONADIABATIC COUPLING FUNCTIONS FOR THE 1ST 5 EXCITED 1-SIGMA-G+ STATES OF THE HYDROGEN MOLECULE. J. Chem. Phys. 1994, 100, 444-451.   (36) Mahapatra, S.; Koppel, H.; Cederbaum, L. S. Reactive scattering dynamics on conically intersecting potential energy surfaces: The H+H-2 exchange reaction. J. Phys. Chem. A 2001, 105, 2321-2329.   (37) Schwenke, D. W. Beyond the potential energy surface: Ab initio corrections to the Born-Oppenheimer approximation for H2O. J. Phys. Chem. A 2001, 105, 2352-2360. 89   (38) Mielke, S. L.; Peterson, K. A.; Schwenke, D. W.; Garrett, B. C.; Truhlar, D. G.; Michael, J. V.; Su, M. C.; Sutherland, J. W. H+H-2 thermal reaction: A convergence of theory and experiment. Phys. Rev. Lett. 2003, 91.   (39) Valeev, E. F.; Sherrill, C. D. The diagonal Born-Oppenheimer correction beyond the Hartree-Fock approximation. J. Chem. Phys. 2003, 118, 3921-3927.   (40) Schuurman, M. S.; Muir, S. R.; Allen, W. D.; Schaefer, H. F. Toward subchemical accuracy in computational thermochemistry: Focal point analysis of the heat of formation of NCO and H,N,C,O isomers. J. Chem. Phys. 2004, 120, 11586-11599.   (41) Tajti, A.; Szalay, P. G.; Csaszar, A. G.; Kallay, M.; Gauss, J.; Valeev, E. F.; Flowers, B. A.; Vazquez, J.; Stanton, J. F. HEAT: High accuracy extrapolated ab initio thermochemistry. J. Chem. Phys. 2004, 121, 11599-11613.   (42) Gauss, J.; Tajti, A.; Kallay, M.; Stanton, J. F.; Szalay, P. G. Analytic calculation of the diagonal Born-Oppenheimer correction within configuration-interaction and coupled-cluster theory. J. Chem. Phys. 2006, 125.   (43) Czako, G.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. Accurate ab initio potential energy surface, dynamics, and thermochemistry of the F+CH4 -> HF+CH3 reaction. J. Chem. Phys. 2009, 130.   (44) Hu, C. P.; Sugino, O.; Watanabe, K. Second-order nonadiabatic couplings from time-dependent density functional theory: Evaluation in the immediate vicinity of Jahn-Teller/Renner-Teller intersections. J. Chem. Phys. 2011, 135.   (45) Morelli, J.; Hammes-Schiffer, S. Surface hopping and fully quantum dynamical wavepacket propagation on multiple coupled adiabatic potential surfaces for proton transfer reactions. Chemical Physics Letters 1997, 269, 161-170.   (46) Yarkony, D. R. Nuclear dynamics near conical intersections in the adiabatic representation: I. The effects of local topography on interstate transitions. Journal of Chemical Physics 2001, 114, 2601-2613.   (47) Gherib, R.; Ryabinkin, I. G.; Izmaylov, A. F. Why Do Mixed Quantum-Classical Methods Describe Short-Time Dynamics through Conical Intersections So Well? Analysis of Geometric Phase Effects. J. Chem. Theory Comput. 2015, 11, 1375-1382.   (48) Abedi, A.; Agostini, F.; Gross, E. K. U. Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion. Epl 2014, 106.   (49) Ryabinkin, I. G.; Joubert-Doriol, L.; Izmaylov, A. F. When do we need to account for the geometric phase in excited state dynamics? J. Chem. Phys. 2014, 140.  90   (50) Berry, M. V. QUANTAL PHASE-FACTORS ACCOMPANYING ADIABATIC CHANGES. Proc. R. Soc. London Ser. A-Math. Phys. Eng. Sci. 1984, 392, 45-57.   (51) Mead, C. A. ELECTRONIC HAMILTONIAN, WAVE-FUNCTIONS, AND ENERGIES, AND DERIVATIVE COUPLING BETWEEN BORN-OPPENHEIMER STATES IN THE VICINITY OF A CONICAL INTERSECTION. J. Chem. Phys. 1983, 78, 807-814.   (52) Varandas, A. J. C.; Xu, Z. R. On the behavior of single-surface nuclear wavefunctions in the vicinity of the conical intersection for an X-3 system. Chem. Phys. Lett. 2000, 316, 248-256.   (53) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. Phys. Rev. Lett. 2010, 105.   (54) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction. J. Chem. Phys. 2012, 137.   (55) Cederbaum, L. S. The exact molecular wavefunction as a product of an electronic and a nuclear wavefunction. J. Chem. Phys. 2013, 138.   (56) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. Phys. Rev. Lett. 2010, 105.     91  CHAPTER 5: THE BEST OF BOTH REPSDIABATIZED GAUSSIANS ON ADIABATIC SURFACES  5.1 Introduction Nonadiabatic dynamics near conical intersections (CIs) are modeled using electronic representations which have well-known tradeoffs.  In Chapter 4 we discussed many of these tradeoffs, with a focus on discontinuities in the adiabatic electronic wave functions near CIs.  We concluded that both the adiabatic and quasi-diabatic electronic representations present significant challenges for nonadiabatic dynamics.  In the present chapter we define a variant of the ab initio multiple spawning method which addresses these challenges and is continuous at CIs, known as diabatized Gaussians on adiabatic surfaces (DGAS).   We begin by reviewing the desirable properties and challenges of the adiabatic and quasi-diabatic representations.  After discussing these challenges, we introduce DGAS and discuss its implementation in the ab initio multiple spawning method for nonadiabatic molecular dynamics.  We analyze the performance of DGAS near conical intersections and trivially unavoided crossings between electronic states through comparison with simulations in the adiabatic representation.  We focus our analysis on changes in the rate of nonradiative decay when a continuous wave function is used.   One of the most significant consequences of discontinuities in the adiabatic representation is the occurrence of large spikes in the diagonal Born-Oppenheimer correction (DBOC) to the adiabatic potential energy surface, however the gradient of the DBOC is too computationally costly to evaluate analytically during dynamical simulations.  We will present and implement an approach for approximate evaluation of the gradient of the DBOC during on-the-fly dynamical 92  simulations.  This approach uses the branching coordinates for a conical intersection to construct a model potential energy surface for analytical differentiation.  We will analyze the error incurred with this approach through comparison of approximate gradients with numerical gradients in the region of a CI, where the DBOC is expected to be largest.  We conclude by reviewing advantages of the DGAS electronic representation.  5.2 Electronic representations for nonadiabatic dynamics  The two most commonly adopted representations for the electronic wave function are the adiabatic and quasi-diabatic representations.  A less widely adopted but very interesting representation arises from the exact factorization of the total molecular wave function as the product of a single nuclear wave function with a single electronic wave function.1-3  The adiabatic representation offers important advantages, particularly in the context of ab initio molecular dynamics simulations.  The adiabatic electronic states are those that are obtained straightforwardly by diagonalization of the electronic Hamiltonian.  However, expansion of the total molecular wave function as a sum over vibronic wave functions in the adiabatic electronic basis, , (5.1) is not without problems.  (Here  is the adiabatic electronic wave function of state J,  is the associated nuclear wave function, r and R are the electronic and nuclear coordinates, respectively, and t is time.)  Because the adiabatic electronic wave functions change discontinuously in R at conical intersection points, matrix elements involving the derivative of the adiabatic electronic wave functions with respect to R become singular.  As a consequence, the first-derivative coupling between intersecting states, 93  , (5.2) and the diagonal elements of the second-derivative coupling, known as the diagonal Born-Oppenheimer correction (DBOC), , (5.3) are singular at intersection points.  The singularities in the first-derivative couplings are integrable, and various strategies for handling them gracefully in numerical simulations of nonadiabatic dynamics are gaining popularity.4-9  However, the singularities in the second derivative couplings are not integrable, and therefore more difficult to incorporate in practice.    In fact, singularities in the DBOC are reflections of the discontinuity of the adiabatic vibronic basis functions, , in which the molecular wave function is expanded.   We have recently shown that representation of the majority of well-behaved molecular wave functions via the expansion in Eq. (5.1)  requires the adiabatic nuclear wave function, , to be discontinuous at the intersection points;10 continuity of the total molecular wave function is achieved by cancellation of compensating discontinuities in the nuclear and electronic wave functions of the intersecting states.  Upon such compensation and only upon such compensation, the singular DBOC contributions cancel, resulting in a finite molecular energy.  Managing these discontinuities and singular matrix elements is impractical in the context of a real calculation.  This formal analysis of the continuity of the total molecular wave function suggests that more accurate results are obtained when the DBOC is neglected in certain classes of approximate quantum dynamical simulations and mixed quantum classical simulations.10  After all, per the postulates of quantum mechanics, the total molecular wave function must be continuous, and the singular DBOC arises from discontinuities in the vibronic basis.  Izmaylov and coworkers have 94  also argued for the neglect of the DBOC in many cases on the grounds that a) errors due to the neglect of the DBOC are canceled by errors introduced by neglect of geometric phase, and b) in many cases simulations of dynamics at conical intersections which include the DBOC show greater errors than those that neglect it.11-13  Despite the mounting evidence against inclusion of the DBOC in most practical approximate simulations of nonadiabatic dynamics, it is unsettling to approximate divergent matrix elements as zero.  In order to circumvent the difficulties of the adiabatic representation, many adopt a quasi-diabatic electronic representation.14-24  Though a completely diabatic basisi.e. an electronic representation in which the derivative couplings are exactly zerocannot be constructed without introducing unphysical features to the electronic wave function, a quasi-diabatic representation eliminates the singularities at conical intersections, leaving behind a small non-removable nonadiabatic coupling.16  Unfortunately, there are infinitely many possible quasi-diabatic representations, leaving the user a difficult choice.  The optimal recipe for a well-behaved global quasi-diabatic representation for a particular chemical problem may not be obvious, even when one has significant prior knowledge of the photodynamics of the system in question.  A poorly chosen quasi-diabatic representation may introduce energy discontinuities and singularities in the nonadiabatic couplings at geometries other than the conical intersection seam, for example.25  Beyond their non-uniqueness, quasi-diabatic states are often coupled over a broader swath of nuclear configuration space than adiabatic states, complicating the simulation of electronic transitions.  Thus, it would be ideal to develop an approach to modeling nonadiabatic dynamics that eliminates the disadvantages of the adiabatic representationnamely the discontinuities in the vibronic basis functions and their consequenceswhile maintaining its advantagesthe relative 95  locality of the coupling between electronic states and the absence of the need to choose a quasi-diabatic basis that is globally well-behaved.  With these goals in mind, we present a new quantum dynamical ansatz: diabatized Gaussians on adiabatic surfaces (DGAS).  DGAS is a variant of the ab initio multiple spawning (AIMS) method for modeling nonadiabatic molecular dynamics26,27 that allows us to simultaneously achieve the best features of both electronic representations.  In this work we will primarily focus on one of the most troubling features of the expansion of the total molecular wave function in the adiabatic representation: the divergent DBOC at conical intersections.  To this end we compare simulation results obtained using three variants of AIMS:  the newly introduced DGAS approach, traditional AIMS in the adiabatic representation (adiabatic AIMS) in which the DBOC is neglected, and a variant of adiabatic AIMS including the DBOC (adiabatic AIMS+DBOC).  In Section 5.3 we will introduce the DGAS approach, and discuss the details of its implementation in AIMS.  We will also discuss the implementation of the adiabatic AIMS+DBOC method, as well as the simulation parameters used and the systems studied in the present work.  In Section 5.4 we will apply traditional adiabatic AIMS, DGAS, and adiabatic AIMS+DBOC to model the photochemistry of a simple molecule, ethene (ethylene), in order to explore the role that wave function continuity and the DBOC play in approximate simulations of nonadiabatic dynamics near conical intersections.  In Section 5.5 we quantify the errors in the approximations made in the adiabatic AIMS+DBOC approach.  In Section 5.6 we apply adiabatic and DGAS AIMS to a long-range separated dimer in order to explore the role that wave function continuity plays in dynamics near trivial crossings. In Section 5.7 we will draw conclusions and discuss future prospects for the DGAS method.  96  5.3 Diabatized Gaussians on adiabatic surfaces (DGAS)  In this section we will introduce the DGAS variant of the AIMS method.  In subsection 5.3.1 we will review the AIMS method as it is traditionally presented.  In subsection 5.3.2 we will introduce a generalized definition of AIMS which encompasses both the widely used AIMS method and the novel DGAS approach introduced in this work.  In subsections 5.3.3-5.3.4 we will introduce the DGAS formalism and discuss some important details of its implementation.  5.3.1 Review of the AIMS method   AIMS is a hierarchy of approximations that allows the solution of the full time-dependent Schrodinger equation (TDSE) for the wave function of an electronically excited molecule, including all electronic and nuclear coordinates explicitly.  The AIMS molecular wave function is typically written as an expansion in Born-Oppenheimer vibronic states, , (5.4) where  may be chosen to be either an adiabatic or quasi-diabatic electronic wave function.  Most often the adiabatic basis is chosen for the reasons described in Section 5.1, giving the expansion in Eq. (5.1).  Hereafter we will refer to this variant of AIMS as adiabatic AIMS.  This expansion is exact, of course; approximations are made in both the definition of  and .  The electronic wave functions, , are computed on the fly via ab initio electronic structure calculations.  Each nuclear wave function is expanded as a linear combination of frozen Gaussian trajectory basis functions (TBFs), 97   (5.5) with  (5.6) and . (5.7) Here  is the time-dependent number of TBFs on electronic state J,  j indexes TBFs,  is the total number of nuclear degrees of freedom, and  indexes such degrees of freedom.  The user-chosen Gaussian widths, momentum vectors of each TBF are represented by  and , respectively, and are propagated via the classical equations of motion  (5.8) and  . (5.9) Here  is the mass associated with nuclear degree of freedom , and , (5.10) where  is the electronic Hamiltonian.  The time-dependent phase,  , is propagated semiclassically 98  . (5.11)  The time-dependent expansion coefficients, , are solved by integration of the TDSE, . (5.12) Here S is the overlap matrix,  , (5.13) H is the full molecular Hamiltonian matrix,  , (5.14) and  is the matrix representation of the right-acting time derivative operator, . (5.15) Note that the second equivalences in equations (5.13) and (5.15) arise from the orthonormality of the electronic basis, whether it be adiabatic or quasi-diabatic.   The Hamiltonian matrix elements can be split into a nuclear kinetic energy term and a nuclear potential energy term arising from the electronic Hamiltonian,  (5.16) where  is the nuclear kinetic energy operator.  Integration over electronic degrees of freedom gives . (5.17) Because TBFs are non-zero for all R, the exact evaluation of matrix elements like that in Eq. (5.17) requires knowledge of the electronic wave functions of states I and J for all R.  Such infinite 99  knowledge is obviously not possible in practical applications.  One can take advantage of the locality of the TBFs, however, to approximate the matrix elements using only local information about the electronic structure.  Specifically, matrix elements over the Hamiltonian are typically computed in the zeroth-order saddle-point approximation (SPA),  (5.18) where . (5.19) The centroid position,  , is the maximum of the absolute overlap density, .  It is trivially defined   (5.20) in the common case that .  As can be seen, computing  in the zeroth-order SPA requires knowledge of the value of , and therefore the electronic structure, at only a single point in nuclear configuration space: .  In AIMS, the SPA is employed to approximate integrals over all operators of the form  (e.g. nonadiabatic coupling matrix elements (NACME) or other electronic observable).  In practice, it is obviously advantageous to minimize the number of TBFs in an AIMS simulation so as to reduce computational cost.  With this goal in mind, the initial nuclear basis is typically chosen to be small (often a single TBF), and then is adaptively expanded during the course of the simulation.  This expansion occurs when the norm of the NACME vector between two electronic states rises above a user defined threshold, through the creation of a new TBF.   This new TBF is created (spawned) in order to accept population transfer from a TBF on a different 100  electronic state.  In this way, TBFs are only added to the nuclear basis set when the algorithm determines that they are required for accurate solution of the TDSE.    Further details on the AIMS spawning algorithm, as well as simplified versions of many of the matrix elements discussed in this section, are provided in Chapter 3.    5.3.2 Generalization of AIMS   In order to introduce the DGAS approximation we first generalize the AIMS formalism.  Rather than associate each basis function with a set of global adiabatic or quasi-diabatic electronic states, as expressed in Eqs. (5.4) and (5.5), we associate a unique, R- and t-dependent electronic wave function, , with each TBF.  The AIMS wave function is then defined . (5.21) The widely used adiabatic AIMS ansatz is a special case of this generalized form of AIMS in which the electronic wave function is taken to be the time-independent adiabatic electronic wave function . (5.22) This generalization of AIMS also encompasses the recently developed ab initio multiple cloning algorithm.28  The generalized AIMS formalism defined in Eq. (5.21) is similar to that defined in Eqs. (5.4) and (5.5) in most respects.  The classical and semiclassical propagation, as defined in Eqs. (5.8)-(5.11) is identical.  The basis is still to be expanded adaptively in the spirit of the spawning procedure described above.  The TDSE (Eq. (5.12)) has the same form, and the matrix elements may still be approximated in the SPA.  What differs is the definitions of the overlap and derivative 101  matrix elements in Eqs. (5.13) and (5.15); the second equality in these equations does not hold because the electronic wave functions may not form an orthonormal set.  Furthermore, calculation of the elements of  is complicated by the time-dependence of the new electronic wave function, . (5.23)  5.3.3 The Diabatized Gaussians on Adiabatic Surfaces Wave Function   The generalized AIMS ansatz presented in the previous subsection gives us the freedom to eschew the primary difficulties associated with the expansion of the total molecular wave function as a sum over adiabatic vibronic wave functions without resorting to a globally quasi-diabatic electronic basis.  To this end, we introduce the diabatized Gaussians on adiabatic surfaces (DGAS) ansatz.  In DGAS, like in adiabatic AIMS, each TBF is associated with an adiabatic electronic state, thus we will index each nuclear wave function with both j and J:   However, the electronic wave function, , associated with  is not the adiabatic electronic wave function.  Instead it is a time-dependent quasi-diabatic electronic wave function, . (5.24) We will call  the DGAS electronic wave function.    The DGAS electronic wave function is anchored to the adiabatic electronic wave function at the time-dependent average position, , according to . (5.25) 102  At all other R, the DGAS electronic wave function does not correspond to a single adiabatic electronic state but instead is defined via a truncated expansion over adiabatic electronic wave functions,  . (5.26) The coefficients in this expansion, which we term the DGAS coefficients, are defined to maximize the overlap of the electronic wave function at all points R with the adiabatic electronic wave function at  of the associated TBF,  . (5.27) In short,  is a quasi-diabatic wave function chosen to be equal to the adiabatic electronic wave function at  and to vary continuously in the surrounding region of nuclear configuration space.  When reasonable simulation parameters are chosen, this DGAS electronic wave function will not exhibit discontinuities in R at conical intersections.  Discontinuities in the adiabatic electronic wave functions, , will be canceled by compensating discontinuities in the R-dependent DGAS coefficients.  Thus, singularities in the first and second derivative coupling matrix elements due to conical intersections do not arise in the DGAS basis.  Discontinuities at conical intersections are the primary disadvantage of the adiabatic electronic basis, and their absence is the primary advantage of the DGAS ansatz.  (A secondary advantage of the DGAS ansatz is that  naturally carries with it geometric phase, which is practically challenging to include in simulations in the adiabatic representation.) 103   The truncation of the summation in Eq. (5.26) to the lowest  adiabatic states appears on its face to be an approximation, but in fact truncation is required to achieve an accurate wave function.  In the limit that  is increased to include all possible many electron wave functions, the DGAS electronic wave functions become truly diabatic, with  no longer depending on R.  In this extreme, electron nuclear correlation is completely absent; the cusps in the electronic wave function will not follow the nuclei.  The complete absence of electron-nuclear correlation would result in unphysically large energies, thus truncation is essential.  Though  is continuous in R at conical intersections, the time dependence of the DGAS electronic wave function introduces a new complexity.  As seen in Eq. (5.23), the time-derivative of the electronic wave function contributes to .  Application of the zeroth order SPA to the term containing this derivative yields . (5.28) When the centroid of a TBF passes directly over a conical intersection, the adiabatic electronic wave function at  changes discontinuously in time.  Because the DGAS electronic wave function is anchored to this adiabatic electronic wave function at , it will also change discontinuously in time.  In this case, the matrix elements may become singular.  Thus, in moving from adiabatic AIMS to DGAS we have replaced a basis that is discontinuous in R with one that is discontinuous in t.  This discontinuity in t is advantageous relative to a discontinuity in R, however, because a discontinuity in t does not translate into divergent energies on the diagonal of the Hamiltonian.  In other words, discontinuities in t do not result in divergent DBOC matrix elements.   104   Nonetheless, the  matrix element in Eq. (5.28) is almost identical to the first-derivative coupling matrix elements that arise from the nuclear kinetic energy operator in adiabatic AIMS, and thus will exhibit an integrable singularity when one of the interacting TBFs passes directly over a conical intersection.  With these singularities in mind, we have recently developed an efficient approach to integrate the TDSE based on a norm-preserving interpolation (NPI) of the adiabatic electronic wave functions.  Integrating the TDSE with NPI ensures the prediction of accurate population transfer probabilities even in the presence of such singularities.8,9  In subsection 5.3.4 we describe the adaptation of this approach to the DGAS ansatz.  5.3.4 Calculation of Derivative Couplings in DGAS   Rapid changes in the adiabatic electronic wave function can result in sharp spikes in the time-derivative nonadiabatic coupling (as reflected in the  matrix) near conical intersections.  Such spikes become singularities if the centroid of one of the interacting TBFs passes directly over the intersection itself.  The integration of the TDSE in the presence of such spikes is numerically challenging; when implemented naively, errors in predicted population transfer probabilities may be large when the coupling varies on a time scale shorter than the time step.  Such errors can be reduced by application of an adaptive or multiple time step integrator in most cases,29,30 however this comes at a high computational cost, and in the worst cases such a short time step is required that the calculation becomes intractable.  As discussed in Chapters 2 and 3, NPI addresses this problem by interpolating the electronic wave function between time steps.  With this interpolated electronic wave function one can compute the time-derivative coupling as a continuous function of time.   105   Here we describe a modification of our original formulation of NPI designed to compute the nonadiabatic coupling term in Eq. (5.28).  Consider a single simulation time step from time  to .  The energy and electronic wave function are known at these two times, but no intermediate time.  Here we consider configuration interaction (CI) wave functions, but this formalism could trivially be generalized to other wave functions.  The CI vector of adiabatic state J at position  will be represented .  The molecular orbital basis at time  are rotated to maximize the overlap with those at  (diabatized), and we then approximate the overlap between orbital m at time  and orbital n at time  as .  The electronic overlap matrix elements between the adiabatic electronic wave functions of states I and J at times and , respectively, are thus well approximated as the inner products of CI vectors, (5.29) The CI vector of the DGAS electronic wave function at the centroid geometry can therefore be defined in terms of the DGAS coefficients, . (5.30)  Now we define the interpolated DGAS electronic wave function for all times, , between  and , , (5.31) where  , (5.32) , (5.33) 106  , (5.34) , (5.35) and .  (5.36) In short, Eqs. (5.31)-(5.36) define a linear rotation of the CI vector from its value at time  to that at .  With the CI vectors now defined as continuous functions of time, the time-derivative coupling between DGAS electronic wave functions can be computed as a continuous function of time . (5.37) The derivative coupling at the midpoint of the time step is approximated as the average of this continuously-defined coupling over the time step .(5.38) In our DGAS simulations, it is this average value that is employed in Eq. (5.28).  As with our previous formulation of NPI, the advantage of this approach is that the entire change in the 107  step is shorter than the timescale for decoherence and dephasing.8  This adaptation of NPI shares the simplicity of implementation of the version reported in our previous publications.8,9  The differentiation and integration in Eq. (5.38) are performed analytically leading to a relatively simple set of expressions for the time-derivative coupling.  These expressions, which are reported in Appendix A, require minimal input from electronic structure calculations: CI vectors of the adiabatic electronic states at , , , , , and .  As in adiabatic AIMS, nonadiabatic coupling terms also arise from the nuclear kinetic energy operator, , in DGAS.  However, due to the quasi-diabatic nature of the DGAS electronic wave functions, the singularities in the DBOC and the first-derivative couplings that exist at conical intersections in adiabatic AIMS are eliminated.  What is left behind is the non-removable portion of the nonadiabatic couplings.  This non-removable portion is small in general,31 and as such, we approximate the nonadiabatic coupling terms arising from the Hamiltonian to zero,   (5.39) and . (5.40) In order to test the validity of this approximation we will evaluate the left-hand side of eqn. (5.39) as 108  .  (5.41)  Here the subscript r is used to denote integration over the electronic coordinates.  We approximate the left-hand side of eqn. (5.40) as  , (5.42) for select simulations.  This evaluation is discussed in more detail in Appendix B.  5.3.5 Implementation of DBOC in Adiabatic AIMS   Separate from DGAS, in this work we introduce a second modification to the AIMS method: the addition of approximate DBOC terms to the PESs employed in adiabatic AIMS.  We introduce this method to quantify the unphysical effects of the wave function discontinuities that arise at conical intersections when the total molecular wave function is expanded as a sum over vibronic states in the adiabatic representation.  The only change to traditional adiabatic AIMS is that we add the DBOC to the PES of each state, i.e. we replace the Born-Oppenheimer surfaces with the Born-Huang surfaces,  . (5.43) Each term in the summation can be expanded as a sum over all electronic states,  . (5.44) 109  Here  is the energy of adiabatic electronic state I.  Despite the fact that analytic techniques for computing second derivative coupling terms are known32,33 and the role of the DBOC has been explored in molecular dynamics simulations based on analytic PESs on several occasions,12,34-38 implementation of the DBOC in the context of ab initio molecular dynamics simulations (AIMD) is impractical for several reasons: a) the summation over all electronic states in Eq. (5.44) is intractable, and b) the analytic differentiation of the Born-Huang PES requires analytic differentiation of this expression, which in turn requires computationally expensive second derivative terms.  Here we introduce approximations to circumvent these difficulties.  First, we truncate the summation in Eq. (5.44) to include only those states that we directly simulate in the dynamics.  This truncation is justified by the fact that only those terms which involve intersecting states I and J become singular, thus these terms dominate the DBOC near conical intersections.  To circumvent the need to compute second derivative terms, we approximate the gradient of the DBOC.  Specifically, to approximate the gradient at position, , we build an analytical model cone approximation to the PES in the region surrounding that point.  The model cone PESs are the adiabatic surfaces obtained by diagonalization of a two-state Hamiltonian in which all matrix elements are Taylor expanded to first order in R,  (5.45) where , (5.46) , (5.47) and 110  . (5.48) At , the model Hamiltonian and exact Hamiltonian predict exactly the same contribution to the DBOC from a given adiabatic electronic state, I, . (5.49) Here  and  are respectively the eigenstates and eigenvalues of .  Obviously the model cone also reproduces the energies and energy gradients of the Born-Oppenheimer PESs exactly at .  Unlike the exact electronic Hamiltonian, however, the model Hamiltonian is easily differentiated to arbitrary order.  Thus, in our simulations we approximate the gradient of the DBOC at  as the derivative of the DBOC computed from the model cone expanded around  . (5.50) We quantify the error in this approximation in Section 5.5.  Even after addition of the DBOC, the electronic Hamiltonian matrix elements, , (5.51) remain finite in the SPA so long as  does not fall exactly on the intersection seam.  In contrast, if  is evaluated exactly, the integral would diverge if an intersection involving state J exists at any R, regardless of , because the singularity in the DBOC is non-integrable and the nuclear basis functions are nonzero at all R.  Thus, in the SPA the divergent DBOC integrals are 111  approximated as finite.  This is arguably a harsh approximation, but still less harsh than the standard practice of complete neglect of the DBOC.   5.3.6 Computational methods and systems of interest   As discussed in the introduction, expansion of the total molecular wave function as a sum over vibronic states in the adiabatic representation is problematic; discontinuities in these vibronic wave functions result in non-integrable singularities in the DBOC.  Most often simulations based on such an expansion neglect the divergent DBOC terms (as in adiabatic AIMS), though through careful approximation one can include them in a practical manner (as in the new adiabatic AIMS+DBOC approach described in subsection 5.3.5).  More ideal is to eliminate the discontinuities altogether, as we do in the DGAS approach (subsection 5.3.3).  In this work we compare DGAS, adiabatic AIMS, and adiabatic AIMS+DBOC by application to a molecule with well-studied photochemistry: ethene.26,39-41  There are two goals of this comparison.  First, in comparing DGAS and adiabatic AIMS+DBOC we wish to quantify the effect of discontinuities in the vibronic basis and the associated energy spikes on nonadiabatic dynamics.  Second, by comparison of DGAS and adiabatic AIMS (without the DBOC) we investigate two different methods of eliminating the singular DBOC energy terms; in DGAS the wave function is carefully constructed to eliminate these singularities, while in adiabatic AIMS the vibronic basis remains discontinuous, but the singular DBOC terms are simply neglected.  We also quantify the effects of discontinuities in the vibronic basis near trivially unavoided crossings (TUCs), a variant of a conical intersection with dimensionality .  TUCs occur when two or more states with zero or effectively zero diabatic coupling intersect.  As discussed in Chapters 2 and 3, TUCs are expected in systems with long-range charge transfer states, as well as 112  in systems that have intersecting electronic states with different spin-symmetry and little or no spin-orbit coupling.  In the present work we model dynamics near trivial crossings using two systems (Figure 5.1).  The long-range separated trifluoromethyl and ethylene dimers were chosen in order to enable an analysis of the accuracy of the DGAS approximation for TUCs, since the correct dynamics in these systems are intuitively obvious.  These dimers have different dynamics as a consequence of the different masses of the terminal atoms (fluorine vs. hydrogen).  The trifluoromethyl dimer was used to demonstrate the accuracy of dynamics with standard Gaussian widths, while the ethylene dimer was used to demonstrate the dynamics of a system with a TUC near the Franck-Condon (equilibrium) region of the PES.  These systems are described further in the discussion of results in Section 5.6.  In all simulations the PES was computed at the state-averaged complete active space self-consistent field (SA-CASSCF) level of theory.  In ethylene the state average includes three states, and an active space of two electrons in two orbitals was used (SA-3-CAS(2/2)).  In the trifluoromethyl and ethylene dimers the state average includes two states, and an active space of three electrons in two orbitals (SA-2-CAS(3/2)) was used.  The 6-31G** basis was used.  All simulations of ethylene were initiated with a single TBF on the S1 state.  All simulations in the trifluoromethyl dimer were initiated with a single TBF on the D1 state, while simulations of the ethylene dimer were initiated with two TBFs on the D0 state, as will be described further in Section 5.6.  The average position and momentum vectors for the initial TBFs in simulations of ethylene and trifluoromethyl dimer were drawn at random from the ground state vibrational Wigner distribution computed in the harmonic approximation at the B3LYP/6-31G** level of theory.  The same set of 50 initial conditions was used in all three sets of ethylene simulations, and in both the adiabatic AIMS and DGAS AIMS simulations for the dimers.  The classical and quantum degrees 113  of freedom were integrated using adaptive integrators (velocity Verlet and 2nd order Magnus expansion, respectively) with a minimum time step of 0.1 au (0.002 fs).  Time-derivative couplings were computed using NPI in all simulations.  The DGAS simulations use NPI as presented in subsection 5.3.4, while the two variants of adiabatic AIMS use NPI as presented in our past work.9  All electronic structure calculations were performed with the MolPro electronic structure software package.42-44  All AIMS calculations (including DGAS) were performed with the FMSMolPro package.45  Figure 5.1 DGAS test systems Systems referenced in the present study.  The trifluoromethyl (A) and ethylene (B) dimers are used to investigate the quantitative effects of discontinuities in the adiabatic electronic wave functions at TUCs.    114  5.4 The effect of enforcing continuity at a conical intersection  The average S0 population of ethylene after excitation to S1 is plotted as a function of the simulation time in Figure 5.2.  Adiabatic+DBOC simulations predict a considerably slower decay to the ground state than DGAS.  Note that the total population on S0 never reaches 1.0 because simulations are stopped when 90% of the population reaches the ground state.)  That inclusion of the DBOC slows population transfer is intuitive.  The DBOC energy is always positive, and as such the spikes in the DBOC at conical intersection points are purely repulsive.  Figure 5.3 shows the Born-Oppenheimer (adiabatic) and Born-Huang (adiabatic+DBOC) PESs for ethylene in the region of the pyramidalized conical intersection that is accessed in these simulations.  It is clear from the Born-Huang surface that the DBOC causes an infinite spike at the conical intersection geometry.  As such, inclusion of the DBOC reduces the probability that a given TBF passes close to the conical intersection seam, where the coupling between adiabatic states is strongest.  We emphasize that these spikes are a consequence of the fact that the adiabatic vibronic wave functions exhibit discontinuities.  They do not exist in the DGAS approximation, which is designed to provide a well-behaved molecular wave function at conical intersection points.  Analysis of the DBOC observed in our dynamical simulations shows that this intuitive picture of repulsion at the DBOC is born out in practice.  Figure 5.4 presents the DBOC experienced by the 50 initially excited TBFs as a function of time.  Spikes in the DBOC above 0.25 eV are regularly seen, with the largest spike close to 1.5 eV.  Clearly significant repulsive forces are attributable to the DBOC.  Figure 5.5 shows the energy gap between S1 and S0 for each of the 50 initially excited TBFs as a function of time.  Smaller energy gaps are observed in the DGAS simulations than in the adiabatic AIMS+DBOC simulations, demonstrating that the DGAS TBFs pass closer to the intersection seam. 115   Figure 5.2 DGAS population transfer in ethylene The coherently averaged S0 population for 50 SA-3-CAS(2,2)/6-31G** AIMS simulations of ethylene is plotted as a function of the simulation time for the adiabatic, adiabatic+DBOC, and DGAS approaches.   Figure 5.3 Ethylene pyramidalized conical intersection PES The SA-3-CAS(2,2)/6-31G** Born-Oppenheimer (A) and Born-Huang (B) PESs are shown in the region of the pyramidalized conical intersection for ethylene.  The Born-Huang PESs are calculated by adding the DBOC to the Born-Oppenheimer energies.  The DBOC is calculated as discussed in Section 5.5. 116   Figure 5.4 DBOC in adiabatic+DBOC simulations of ethylene The DBOC for the first TBF in 50 adiabatic+DBOC simulations of ethylene is shown as a function of the simulation time.   Figure 5.5 Energy gap in DGAS and adiabatic+DBOC simulations of ethylene The log10 of the S0-S1 energy gap (in eV) is plotted for the first TBF in 50 DGAS (A) and adiabatic+DBOC (B) simulations of ethylene, with the log10 of the average gap for all TBFs shown in red.   117   Now we consider the adiabatic AIMS results without the DBOC.  As noted above, adiabatic AIMS is based on the expansion of the total molecular wave function as a sum over vibronic states in the adiabatic representation, and therefore the individual vibronic states exhibit discontinuities at conical intersection points.  The DBOC (which is therefore singular at intersection points) is neglected, however.  It is noteworthy that the S0 population curves for DGAS and adiabatic AIMS (Figure 5.2) are effectively identical.  Thus, the rate of decay is similar regardless of whether the elimination of singularities in the DBOC is achieved by careful construction of a continuous molecular wave function (DGAS) or by simple neglect (adiabatic AIMS).  This supports the widespread practice of neglect of the DBOC in simulations in the adiabatic basis.  Given the classical nature of the TBF basis, this result also supports the neglect of the DBOC in mixed quantum classical simulations (MQC; e.g. surface hopping simulations).  5.5 Analysis of errors in DBOC gradients  Adiabatic AIMS+DBOC simulations are based on the approximation to the gradient of the DBOC defined in Eq. (5.50).  Here we analyze the error in these approximate gradients by comparison to numerical gradients of the DBOC.  Specifically, we have optimized the minimal energy conical intersection (MECI) of ethene at the CASSCF level of theory described in the previous section.  The difference gradient and nonadiabatic coupling vectors at the intersection point were computed, and a grid of geometries displaced from the MECI geometry in these two directions was computed.  The electronic structure was solved at the CASSCF level of theory at each grid point.  The DBOC is then computed at each point considering only the two intersecting states (i.e. the summation in Eq. (5.44) is truncated to two states).  The derivative of the DBOC was then computed analytically according to Eq. (5.50) and numerically by central difference with 118  a step size of 0.001 Å.  This numerical derivative is taken as the reference, and the error in the model cone derivative is computed relative to it.  The DBOC is plotted in Figure 5.6A.  Figure 5.6B shows the relative error in the gradient vector, defined , (5.52) where  is the numerical derivative of the DBOC with respect to R,  is the analytical derivative of the DBOC in the model cone approximation (the right hand side of Eq. (5.50)), and  is the gradient of the Born-Oppenheimer PES of the first excited state.  In other words, this is the norm of the error in the model cone gradient vector divided by the norm of the total force on S1.  Not surprisingly, this relative error is largest in the region where the DBOC itself is large, exceeding 1% of the total gradient norm in roughly the same region where the DBOC exceeds 10-0.5 (0.32) eV.  The maximum error is 6.4% (in the red region in the upper right quadrant of Figure 5.6B).  The error in the gradient results in non-conservation of the classical energy of the 50 initial TBFs in the adiabatic AIMS+DBOC simulations described in the previous section.  The mean absolute deviation (mean deviation) of the energy at the end of the simulations relative to that at the beginning of the simulations is 0.17 (0.06) eV.  Though not insignificant, this is a small fraction of the average total kinetic energy at the end of the simulations (2.99 eV) and thus we expect that it does not affect the conclusions of our study.  Note that these averages include 49 of the 50 initial TBFs.  The remaining TBF experienced a massive gain of energy (36.8 eV) due to an artifact of the CASSCF PES, and was therefore excluded as an outlier.  119   Figure 5.6 DBOC and gradient error in the two-state approximation for ethylene The log10 of the DBOC in eV (A) at points in a region surrounding the pyramidalized MECI of ethylene.  The DBOC is computed in the two-state approximation, as described in the text.  The log10 of the relative error in the DBOC gradient (defined in Eq. (5.52)) in the same region is shown in B.  Note that in both cases the color legend defines the upper bound to the value in that color region (e.g. in A, red indicates a DBOC below 103 eV, white indicates a DBOC below 100.5 eV, and so on).   120  5.6 The effect of enforcing continuity at a trivial crossing  Thus far we have compared AIMS simulations with the adiabatic, adiabatic+DBOC, and DGAS approaches in ethylene, which exhibits a CI.  We will now compare AIMS simulations with the adiabatic and DGAS wave functions near a trivially unavoided crossing (TUC).  We begin by noting an important difference between the simulations we discuss in Section 5.6 and the simulations discussed in Section 5.4: the application of a different spawning algorithm.  The spawning algorithm is responsible for determining how new TBFs are added to the nuclear basis when a CI or TUC is encountered.  In the standard AIMS implementation, as applied in Section 5.4, new TBFs are created at the simulation time when the time-derivative coupling (eqn. (5.38)) between two TBFs is largest.  This is a sensible algorithm for simulating dynamics near CIs, which occupy a significant region of the PES, but is problematic for the study of trivial crossings, which are infinitesimally narrow in R.  We addressed this issue at length when discussing simulations of the trifluoromethyl dimer in Chapter 3.  We showed that in order to completely eliminate long-range charge transfer in the AIMS approach near TUCs it is necessary to create new TBFs where the energy gap between the crossing adiabatic states is smallest.   In order to eliminate long-range charge transfer that arises from the standard spawning algorithm, the simulations in this section were completed using a line search of the intersection region when a TBF experiences a spike in the time-derivative coupling.  The purpose of this line search is to identify the simulation time at which the energy gap between the intersecting electronic states is smallest.  For example, if TBF i, on adiabatic state I, experiences a spike in the time-derivative coupling with state J at time t0, when its velocity is , this line search minimizes the following function with respect to :  . (5.53) 121  Once this function is minimized, TBF i is propagated forward by  and a new TBF is spawned at simulation time . Unphysical long-range transfer of charge may be quantified by plotting the average Mulliken charge of the anion in the trifluoromethyl dimer as a function of the simulation time.  As a consequence of the distance between the molecules any loss of charge from the anion is unphysical.  Figure 5.7 shows that the average Mulliken charge is exactly correct in all simulations with both the adiabatic and DGAS representations for the electronic wave function.  Perfect agreement between the adiabatic and DGAS representations suggests that the discontinuity in the adiabatic basis has no effect in standard dynamical simulations near TUCs.  Figure 5.7 DGAS and adiabatic AIMS population transfer in trifluoromethyl dimer The average Mulliken charge of the anion in 50 AIMS simulations of the trifluoromethyl dimer is shown.  There is no error in either the adiabatic or DGAS representations, so that they are directly on top of one another.    122  The simulation results presented in Figure 5.7 are influenced by another parameter which we now discuss further.  In order to limit the size of the nuclear basis, TBFs that occupy the ground electronic state and do not interact with other TBFs for a user-determined amount of time are removed.  In the simulations discussed here, non-interacting TBFs were removed from the basis after 5 fs.  One might expect to observe somewhat different dynamical behavior if these TBFs were not removed from the nuclear basis.  Unfortunately, such a simulation would be challenging to execute since the number of electronic structure calculations quickly becomes intractable.  In simulations of the trifluoromethyl dimer (Figure 5.1A) the simulation becomes intractable when the number of TBFs reaches 20.  On average, this occurs within the first 50 fs of the simulation time.  Furthermore, linear dependence, which is caused by large overlaps between two TBFs, complicates the solution of eqn. (5.12).  We thus face a choice between a small nuclear basis, in which the TBFs are largely non-interactive, or a large nuclear basis in which more TBFs interact but we cannot solve the TDSE.  This is, in effect, a choice between modeling and not modeling coherence of the TBFs. In order to model coherent TBFs while managing the computational cost, we model the nonadiabatic dynamics of a cationic ethylene dimer (Figure 5.1B).  These simulations were prepared with two TBFs located at the degenerate D0 minima.  These minima differ with respect to which monomer is cationic.  As each TBF is located at a minimum on the PES, they are effectively frozen, with no momentum in any nuclear direction.  This system is depicted in Figure 5.8. 123   Figure 5.8 Long range separated ethylene dimer cation simulation results The initial conditions for the ethylene dimer simulation described in the present work are shown via a one-pts cationic.  Thick red and blue lines depict intersecting adiabatic state PESs at a TUC.  Thinner lines depict TBFs at each of the D0 minima.  Note that the TBFs interact when confined to these minima (labeled by points on the PES).  Figure 5.9 shows the average Mulliken charge of the cation in the ethylene dimer as a function of the simulation time using the adiabatic and DGAS electronic representations.  As time progresses, the simulation in the adiabatic representation gives an oscillatory value for the Mulliken charge of the monomer that was initially cationic.  This is unphysical, since the monomers are separated by a distance of 1000 Å, and stems from population exchange between the TBFs at each D0 minimum.  The same simulation in the DGAS representation exhibits no population exchange.  This is a consequence of the fact that the DGAS electronic wave function is continuous with respect to displacements in the nuclear coordinates at TUCs, while the adiabatic wave function is discontinuous. We can better understand the performance of the DGAS and adiabatic representations by reviewing the range of possible values for the wave function overlaps in each approach.  The two 124  TBFs in the ethylene dimer model discussed here occupy the same adiabatic state.  They may then exchange population in the adiabatic representation through the overlap matrix since, by definition,  . (5.54) By contrast, in DGAS the electronic wave function overlaps may have the following range of values (when the wave function is real):  , (5.55) where i and j are different TBFs.  This flexibility is essential to guaranteeing continuity in the molecular wave function near CIs and TUCs while maintaining a continuous definition for the nuclear wave function.  Disagreement between the DGAS and adiabatic simulation data in Figure 5.9 is a result of uncompensated discontinuities in the adiabatic wave functions at a TUC.     Figure 5.9 Long range separated ethylene dimer cation simulation results The Mulliken charge of the cationic monomer in AIMS simulations of an ethylene dimer.   Results obtained using the DGAS (blue) and adiabatic (black) representations for the electronic wave function are shown.  125   Figure 5.9 suggests that the DGAS representation is superior to the adiabatic representation for dynamics in which TBFs are located near a TUC, particularly if this TUC is near the Franck-Condon (equilibrium) region of the PES.  We demonstrated that this superiority makes DGAS advantageous for the study of long-range charge transfer dimers.  We will now discuss applications in which these advantages would be useful.   Donor-acceptor architectures are now quite common in organic solar devices.  The simulation of non-radiative decay in these devices is likely to involve trivial crossings between wave functions that are localized on spatially separated molecules.  We would expect the DGAS representation to enhance our understanding of different mechanisms for nonradiative decay in these highly complex systems.   In the field of inorganic photovoltaics, the introduction of transition metals often results in a manifold of electronic states with different spin-symmetry and spin-orbit couplings.  We expect the advantages of the DGAS representation to become apparent when simulating the dynamics of these systems as well, since TUCs are expected between electronic states with weak spin-orbit coupling.  5.7 Conclusions  In this work we have introduced DGAS, a new approximate quantum dynamical ansatz that simultaneously achieves the advantages of expanding in the adiabatic and quasi-diabatic bases.  Like the quasi-diabatic representation, A) the total molecular wave is expanded in a basis of vibronic states that is well behaved at conical intersections, and B) all matrix elements are finite.  However, similar to the adiabatic representation: C) there is no need to choose a globally quasi-diabatic representation of the electronic wave function, and D) the couplings between electronic states spike intermittently and are negligibly small the remainder of the time. 126   For comparison we introduced a variant of the AIMS method in which the kinetic energy associated with discontinuities in the vibronic basis is explicitly included.  This kinetic energy takes the form of infinite spikes in the Born-Huang surfaces at conical intersections.  In practice these spikes repel population from the conical intersection seam, slowing population transfer to the ground state significantly in simulations of photoexcited ethene.  Elimination of these spikes either by rigorously eliminating wave function discontinuities (DGAS) or by simple neglect (adiabatic AIMS) resulted in similarly fast lifetimes.  In order to quantify the effect of discontinuities in the adiabatic surfaces on the dynamical rate of nonradiative decay near TUCs we simulated the dynamics of a trifluoromethyl dimer, and found that both the adiabatic and DGAS representations give exact results when an appropriate spawning algorithm is applied.  Simulations of the ethylene dimer, however, suggest that the DGAS representation is more accurate in certain long range charge transfer systems.  This accuracy is derived from the fact that the electronic wave functions are continuous with respect to nuclear motion, and suggests that the DGAS representation may be advantageous for the simulation of interfaces that are important in emerging solar energy technologies.  In closing, we highlight the key idea behind DGAS, which we hope may stimulate developments in other areas:  Many quantum dynamical ansatzes based on the adiabatic representation are designed such that the molecular wave function itself and/or the basis in which that wave function is expanded are discontinuous in R, but vary continuously in time.  In contrast, the DGAS vibronic basis is chosen to vary continuously in R, but may change discontinuously in time.  In this way, all singularities in the derivative coupling terms are eliminated from the Hamiltonian.  Instead, singular terms that are akin to first-derivative couplings arise on the time-derivative side of the TDSE.  Because these singularities are integrable and occur only in off-127  diagonal matrix elements, they are far easier to deal with in practical simulations than the divergent DBOC matrix elements that arise in the adiabatic representation.      128             APPENDICES      129  APPENDIX A: Codable expressions for NPI Time-Derivative Couplings in DGAS  The NPI derivative couplings are defined in eqn. (5.38),  . (A5.1) By inserting the definition of  (Eqs. (5.30)-(5.33)) into Eq. (A5.1), analytically differentiating  with respect to , and then analytically integrating over , one obtains the following codable expressions for the time-derivative coupling.  (5.57)  (5.58)  (5.59)  (5.60)  (5.61)  (5.62)  (5.63)  (5.64) 130  APPENDIX B: Mathematical expressions for the first derivative coupling in DGAS   In the expressions that follow, unless otherwise noted, integration is over the electronic degrees of freedom.  The DGAS first derivative coupling (FDC) is defined as (eqn. (5.41)): .   (5.65) Analytical evaluation of the nonadiabatic coupling vector is straight-forward, and may be done trivially for small molecules.  The only term in eqn. (5.65) that remains to be defined is .  To simplify the following expressions, we define two new quantities, SI and:    (5.66)  . (5.67) Evaluation of  is not straight-forward when .  This can be overcome through insertion of resolution of the identity  , (5.68) in eqn. (5.67) as follows:  .   (5.69)  In terms of SI, the DGAS coefficient is defined as 131   . (5.70) Using the chain rule, and introducing a new variable, B, whose definition is  , (5.71) we may express  as:  . (5.72) Using the product rule:  . (5.73) Using the chain rule again for the terms in eqn. (5.73), we obtain  , (5.74) and  . (5.75) Using eqns. (5.73)-(5.75), eqn. (5.72) becomes  .    (5.76) One would think that  in the denominator would present issues when this coefficient has a value of zero, however, the Kronecker delta in eqn. (5.65) eliminates this possibility.  132  Substituting eqn. (5.76) into eqn. (5.65), we obtain a solvable expression for the DGAS first derivative coupling. APPENDIX C: First and second derivative couplings for DGAS simulations of ethylene   The FDC between the first and second TBFs in 50 DGAS simulations of ethylene are plotted as a function of the simulation time in Figure A5.1.  This data shows that the FDC with the DGAS representation is negligible for simulations near a CI, justifying its omission during a standard implementation.  The approximate second derivative coupling (SDC) is shown in Figure A5.2.  Calculated as the square of the FDC, which is already negligible, the SDC is insignificant.  This justifies the omission of the SDC during dynamical simulations which apply the DGAS representation.  Figure A5.1 First derivative couplings in DGAS simulations of ethylene The FDC between the first parent and child basis functions is plotted as a function of the simulation time for 50 DGAS simulations of ethylene.  Note that the y-axis scale is in eV, and that their small values suggest that ignoring this term is a good approximation. 133   Figure A5.2 Second derivative couplings in DGAS simulations of ethylene The approximate SDC between the first parent and child basis functions is plotted as a function of the simulation time for 50 DGAS simulations of ethylene.  Note that the y-axis scale is in eV, and that their small values suggest that ignoring this term is a good approximation.   134             REFERENCES    135  REFERENCES   (1) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Exact Factorization of the Time-Dependent Electron-Nuclear Wave Function. Phys. Rev. Lett. 2010, 105.   (2) Abedi, A.; Maitra, N. T.; Gross, E. K. U. Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction. J. Chem. Phys. 2012, 137.   (3) Cederbaum, L. S. The exact molecular wavefunction as a product of an electronic and a nuclear wavefunction. J. Chem. Phys. 2013, 138.   (4) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. Journal of Chemical Physics 2001, 114, 10608-10615.   (5) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. J. Chem. Phys. 2012, 137.   (6) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. Chemical Physics Letters 2013, 590, 208-213.   (7) Wang, L. J.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. Lett. 2014, 5, 713-719.   (8) Meek, G. A.; Levine, B. G. Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations. Journal of Physical Chemistry Letters 2014, 5, 2351-2356.   (9) Meek, G. A.; Levine, B. G. Accurate and efficient evaluation of transition probabilities at unavoided crossings in ab initio multiple spawning. Chem. Phys. 2015, 460, 117-124.   (10) Meek, G. A.; Levine, B. G. Wave function continuity and the diagonal Born-Oppenheimer correction at conical intersections. J. Chem. Phys. 2016, 144.   (11) Ryabinkin, I. G.; Joubert-Doriol, L.; Izmaylov, A. F. When do we need to account for the geometric phase in excited state dynamics? J. Chem. Phys. 2014, 140.   (12) Gherib, R.; Ryabinkin, I. G.; Izmaylov, A. F. Why Do Mixed Quantum-Classical Methods Describe Short-Time Dynamics through Conical Intersections So Well? Analysis of Geometric Phase Effects. J. Chem. Theory Comput. 2015, 11, 1375-1382.  136   (13) Gherib, R.; Ye, L. Y.; Ryabinkin, I. G.; Izmaylov, A. F. On the inclusion of the diagonal Born-Oppenheimer correction in surface hopping methods. J. Chem. Phys. 2016, 144.   (14) Baer, M. ADIABATIC AND DIABATIC REPRESENTATIONS FOR ATOM-DIATOM COLLISIONS - TREATMENT OF 3-DIMENSIONAL CASE. Chem. Phys. 1976, 15, 49-57.   (15) Macias, A.; Riera, A. CALCULATION OF DIABATIC STATES FROM MOLECULAR-PROPERTIES. Journal of Physics B-Atomic Molecular and Optical Physics 1978, 11, L489-L492.   (16) Mead, C. A.; Truhlar, D. G. CONDITIONS FOR THE DEFINITION OF A STRICTLY DIABATIC ELECTRONIC BASIS FOR MOLECULAR-SYSTEMS. Journal of Chemical Physics 1982, 77, 6090-6098.   (17) Pacher, T.; Cederbaum, L. S.; Koppel, H. APPROXIMATELY DIABATIC STATES FROM BLOCK DIAGONALIZATION OF THE ELECTRONIC HAMILTONIAN. J. Chem. Phys. 1988, 89, 7367-7381.   (18) Cave, R. J.; Newton, M. D. Generalization of the Mulliken-Hush treatment for the calculation of electron transfer matrix elements. Chem. Phys. Lett. 1996, 249, 15-19.   (19) Atchity, G. J.; Ruedenberg, K. Determination of diabatic states through enforcement of configurational uniformity. Theor. Chem. Acc. 1997, 97, 47-58.   (20) Yarkony, D. R. On the adiabatic to diabatic states transformation near intersections of conical intersections. Journal of Chemical Physics 2000, 112, 2111-2120.   (21) Mota, V. C.; Varandas, A. J. C. HN2((2)A ') electronic manifold. II. Ab initio based double-sheeted DMBE potential energy surface via a global diabatization angle. J. Phys. Chem. A 2008, 112, 3768-3786.   (22) Evenhuis, C.; Martinez, T. J. A scheme to interpolate potential energy surfaces and derivative coupling vectors without performing a global diabatization. J. Chem. Phys. 2011, 135.   (23) Zhu, X. L.; Yarkony, D. R. Fitting coupled potential energy surfaces for large systems: Method and construction of a 3-state representation for phenol photodissociation in the full 33 internal degrees of freedom using multireference configuration interaction determined data. J. Chem. Phys. 2014, 140.   (24) Hoyer, C. E.; Xu, X. F.; Ma, D. X.; Gagliardi, L.; Truhlar, D. G. Diabatization based on the dipole and quadrupole: The DQ method. J. Chem. Phys. 2014, 141.   (25) Zhu, X. L.; Yarkony, D. R. On the Construction of Property Based Diabatizations: Diabolical Singular Points. J. Phys. Chem. A 2015, 119, 12383-12391. 137   (26) Ben-Nun, M.; Quenneville, J.; Martinez, T. J. Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics. Journal of Physical Chemistry A 2000, 104, 5161-5175.   (27) Ben-Nun, M.; Martinez, T. J. Ab initio quantum molecular dynamics. Advances in Chemical Physics, Volume 121 2002, 121, 439-512.   (28) Makhov, D. V.; Glover, W. J.; Martinez, T. J.; Shalashilin, D. V. Ab initio multiple cloning algorithm for quantum nonadiabatic molecular dynamics. J. Chem. Phys. 2014, 141.   (29) Wang, L. J.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals. Journal of Physical Chemistry Letters 2013, 4, 1888-1894.   (30) Virshup, A. M.; Levine, B. G.; Martinez, T. J. Steric and electrostatic effects on photoisomerization dynamics using QM/MM ab initio multiple spawning. Theoretical Chemistry Accounts 2014, 133, 1506.   (31) Kendrick, B. K.; Mead, C. A.; Truhlar, D. G. Properties of nonadiabatic couplings and the generalized Born-Oppenheimer approximation. Chemical Physics 2002, 277, 31-41.   (32) Lengsfield, B. H.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS FOR MCSCF CL WAVE-FUNCTIONS USING ANALYTIC DERIVATIVE METHODS .3. 2ND DERIVATIVE TERMS. Journal of Chemical Physics 1986, 84, 348-353.   (33) Saxe, P.; Yarkony, D. R. ON THE EVALUATION OF NONADIABATIC COUPLING MATRIX-ELEMENTS FOR MCSCF/CL WAVE-FUNCTIONS .4. 2ND DERIVATIVE TERMS USING ANALYTIC GRADIENT METHODS. Journal of Chemical Physics 1987, 86, 321-328.   (34) Morelli, J.; Hammes-Schiffer, S. Surface hopping and fully quantum dynamical wavepacket propagation on multiple coupled adiabatic potential surfaces for proton transfer reactions. Chemical Physics Letters 1997, 269, 161-170.   (35) Yarkony, D. R. Nuclear dynamics near conical intersections in the adiabatic representation: I. The effects of local topography on interstate transitions. Journal of Chemical Physics 2001, 114, 2601-2613.   (36) Shenvi, N. Phase-space surface hopping: Nonadiabatic dynamics in a superadiabatic basis. J. Chem. Phys. 2009, 130.   (37) Abedi, A.; Agostini, F.; Gross, E. K. U. Mixed quantum-classical dynamics from the exact decomposition of electron-nuclear motion. Epl 2014, 106. 138   (38) McKemmish, L. K.; McKenzie, R. H.; Hush, N. S.; Reimers, J. R. Electron-vibration entanglement in the Born-Oppenheimer description of chemical reactions and spectroscopy. Physical Chemistry Chemical Physics 2015, 17, 24666-24682.   (39) Michl, J.; Bonacic-Koutecky, V.: Electronic Aspects of Organic Photochemistry; Wiley: New York, 1990.   (40) Levine, B. G.; Martinez, T. J. Isomerization through conical intersections. Annu. Rev. Phys. Chem. 2007, 58, 613-634.   (41) Mori, T.; Glover, W. J.; Schuurman, M. S.; Martinez, T. J. Role of Rydberg States in the Photochemical Dynamics of Ethylene. Journal of Physical Chemistry A 2012, 116, 2808-2818.   (42) Werner, H. J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schutz, M. Molpro: a general-purpose quantum chemistry program package. WIREs Comput. Mol. Sci. 2012, 2, 242-253.   (43) Knowles, P. J.; Werner, H. J. AN EFFICIENT 2ND-ORDER MC SCF METHOD FOR LONG CONFIGURATION EXPANSIONS. Chem. Phys. Lett. 1985, 115, 259-267.   (44) Werner, H. J.; Knowles, P. J. A 2ND ORDER MULTICONFIGURATION SCF PROCEDURE WITH OPTIMUM CONVERGENCE. J. Chem. Phys. 1985, 82, 5053-5063.   (45) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of ab initio multiple spawning in the MOLPRO quantum chemistry package. Chemical Physics 2008, 347, 3-16.    139  CHAPTER 6: POLARONIC RELAXATION BY THREE-ELECTRON BOND FORMATION IN GRAPHITIC CARBON NITRIDES  This chapter is taken in its entirety from: Meek, G. A.; Baczewski, A. D.; Little, D. J.; Levine, B. Polaronic Relaxation by Three-Electron Bond Formation in Graphitic Carbon NitridesPhys. Chem. C, 2014, 118, 4023.  6.1 Introduction  Carbon nitrides (C3N4) have attracted wide attention for their applications as organic photocatalysts capable of driving water splitting (including both the oxidative and reductive steps),1-3 oxidative organic reactions,4-7 and photodegradation of pollutants.8,9  As a metal-free material, C3N4 offers many advantages over traditional inorganic catalysts, including tunable properties afforded by the range of organic synthetic techniques and compatibility with functional groups that interact strongly with metals.  Selectivities of 99% have been reported for the photocatalytic oxidative coupling of amines6 and oxidation of benzyl alcohol to benzaldehyde4 on C3N4, and water splitting quantum efficiencies as high as 4.8% have also been reported.2 Though carbon nitride can exist as phases where the constituent atoms take on either graphite-like sp2 (graphitic, g-C3N4) or diamond-like sp3 -C3N4) character,10 x-ray diffraction has established that the graphitic form is most often observed under typical synthetic conditions.1  Studies employing density functional theory11,12 (DFT) have established the graphitic phase to be the more thermodynamically stable form under ambient conditions, as well.13-16  Graphitic carbon nitride can exist either as a fully-condensed two-dimensional (2-D) graphitic sheet (Figure 6.1A) 17-19  Both phases are believed to be composed of repeating melem (2,5,8-triamino-1,3,4,6,7,9,9b-heptaazaphenalene; Figure 6.2A) monomers, while phases with identical stoichiometry constructed from smaller melamine (1,3,5-triazine-2,4,6-triamine) units are believed to be less stable.20,21  Many previous theoretical investigations 140  have been devoted to the determination of the most stable structure of the material,13,14,16,21-23 with recent studies suggesting that the Pauli repulsion arising from the close contact between lone pair electrons on the nitrogen atoms of adjacent monomer units (hereafter referred to as position 1 nitrogens, N1; see Figure 6.2) force the fully-condensed 2-D structure to buckle out-of-plane and the polymeric form to curve within the plane.20-22  A three-dimensional corrugated structure that reduces these unfavorable interactions has also been proposed and demonstrated to be consistent with experimental x-ray diffraction data.21  Given that experimental studies report electronic, optical, and catalytic properties that depend strongly on the conditions of synthesis,19,24 it seems likely that experimentally realized g-C3N4 is a heterogeneous material comprising polymeric, 2-D, and other phases.    Figure 6.1 Phases of graphitic carbon nitrides (g-CN) Chemical structures of the fully condensed 2-D (A) and polymeric (B) phases of graphitic carbon nitride. 141   Figure 6.2 Molecular models of g-CN Chemical structures of molecular models of g-C3N4.  Position 1 nitrogen atoms, whose lone pair electrons interact with neighboring monomers, have    The electronic band structures of idealized 2-D and polymeric g-C3N4 have been elucidated by past DFT calculations.  Interestingly, they are quite different, with an indirect band gap reported for the 2-D phase25,26 and a direct gap reported for the polymer.1  In both phases the conduction band minimum is comprised of carbon-centered pz orbitals (where the z-axis is perpendicular to the plane of the material).1,20  In the polymer, the valence band maximum (VBM) is composed of nitrogen-centered pz orbitals,1 while the VBM of the 2-D material is composed of nitrogen px, py, and s orbitals (in the plane of the graphitic sheet).20  The electronic structure of melem is also well 142  known and, similar to polymeric g-C3N4, is characterized by a nitrogen-centered pz highest occupied molecular orbital (HOMO), a carbon-centered pz lowest unoccupied molecular orbital ese orbitals.27,28  Polaronic relaxation, the stabilization of charge carriers by a local nuclear distortion, reduces the mobility of charge carriers in organic semiconductors, and strong polaronic relaxation can lead to self-trapping.29  Though many theoretical studies have focused on determining the structural morphology and bulk electronic structure of g-C3N4, the possibility that such relaxation slows the transport of charge has yet to be addressed.  Given that recent experiments suggest low mobility as a factor limiting the photocatalytic efficiency in g-C3N4,24 it is certainly important to investigate polaronic relaxation in these materials, and as such, we present a combined DFT and coupled cluster study of polaronic relaxation in positively charged models of carbon nitrides.  We also investigate the role that interactions between lone pairs of neighboring monomers play in determining the electronic properties of g-C3N4.      6.2 Methods 6.2.1 Electronic structure calculations on molecular models   Herein, the behavior of holes in g-C3N4 is analyzed through the study of cations of various molecular models of the bulk material, with polaronic relaxation in g-C3N4 characterized through geometry optimization of these cations.  Because the interaction between position 1 nitrogen lone pairs on adjacent melem units is important in determining the electronic structure of the material, a set of molecular models of g-C3N4 were designed which each include at least one such N1-N1 contact (Figure 6.2B-F).  The melamine dimer, or melam (Figure 6.2B), was studied due to its small size, making it a convenient molecule for the application of ab initio electronic structure 143  methodologies including a highly accurate treatment of electron correlation.  The melem dimer, which consists of two linked melem units (Figure 6.2C), and the polymeric melem trimer, consisting of three melem units in a linear chain (Figure 6.2D), serve as larger models of polymeric g-C3N4. The 2-D melem trimer (Figure 6.2E), which contains three melem units in a triangular orientation, and hexamer (Figure 6.2F), which contains six melem units in a similar triangular structure, serve as models for 2-D g-C3N4.  All dangling bonds left by truncation of the extended structures are capped with hydrogen atoms.  All cation and neutral structures were optimized with the TeraChem30-33 graphics-processing-unit-accelerated quantum chemistry software package using unrestricted (U-) DFT for the cation and restricted DFT for the neutral, with the CAM-B3LYP34 functional and the cc-pVDZ basis set35 employed in both cases.  All optimizations were performed without assumption of symmetry or planarity of the molecular models.  The CAM-B3LYP functional was chosen due to its known accuracy for the prediction of the ionization potentials and polarizabil-conjugated systems.34  All cation doublet Kohn-Sham wave-functions have  values smaller than 0.78 in the results reported here.  Throughout this work the lowest unoccupied Kohn-Sham molecular orbital of the minority spin is reported as the hole orbital created by ionization.  Reported bond orders between pairs of atoms were determined from the off-diagonal elements of the Kohn-Sham density matrix.  Vibrational frequencies and infrared intensities were calculated with the GAMESS electronic structure package36,37 and are reported in Supporting Information for the neutral (cation) minimum structures of melam and melem dimer at the (U-)CAM-B3LYP/cc-pVDZ level of theory.  To assess the errors in our DFT study, the two hole, one particle variant of the ionization potential equation of motion coupled cluster theory with single and double excitations,38-41 IP-144  EOM-CCSD(2h,1p)/cc-pVDZ (henceforth referred to as IP-EOM), as implemented42,43 in GAMESS was applied to the DFT-optimized structures of the smallest model system, melam.   6.2.2 Analysis of nuclear relaxation  In this work, we define the cation relaxation energy, , as the difference between the energy of the cation at the neutral minimum structure, , and the cation energy at the cation minimum structure, , . (6.1) The corresponding definition for the neutral relaxation energy, , is . (6.2) As noted above, an estimate of the error in the U-CAM-B3LYP cation relaxation energy was obtained by computing the relaxation energy of the smallest molecular model at the IP-EOM level of theory.  The difference in relaxation energy between IP-EOM and U-CAM-B3LYP in melam can then be applied to define a corrected relaxation energy for larger systems, , (6.3) with the IP-EOM relaxation energy computed at the DFT-optimized structures.  Corrected cation relaxation energies are reported in conjunction with unmodified U-CAM-B3LYP cation relaxation energies.  Marcus theory of electron transfer provides a simple framework for analysis of the impact of nuclear relaxation on charge carrier mobility.29,44  When the driving force for electron transfer, 145  , is zero, as is the case for the hole self-exchange responsible for charge transport in g-C3N4, the electron-transfer rate-constant () is then related to the reorganization energy, , by . (6.4) In the present work we focus on the inner sphere contribution to the reorganization energy in order to address the impact of nuclear relaxation on hole mobility, with the inner sphere reorganization energy defined as . (6.5) This reorganization energy was determined for all molecular models of g-C3N4.   6.2.3 Periodic DFT calculations  Density functional theory calculations of idealized planar 2-D and polymeric phases of g-C3N4 were performed to correlate changes in N1-N1 distance with the band structure in the extended system.  Between planes, 14 Å of vacuum is utilized to avoid interaction between neighboring layers.  The gross physics reported (e.g., band gaps, orbital character, N1-N1 distances) should not be dramatically different from those of the layered bulk structures.    These periodic calculations were carried out using the Vienna Ab Initio Simulation Package (VASP v 5.3.3),45,46 which uses a plane wave basis with the electron-ion interaction treated efficiently using the projector augmented-wave (PAW) method.47,48  For all calculations, the plane wave cutoff and number of k-points used in sampling the Brillouin Zone (BZ) were adjusted to achieve convergence in the total energy to within 1 meV/atom.  To this end, a 500 eV cutoff and distributions of 5 and 6 points in the irreducible BZ (25 and 12 in the first BZ) were used for the 2-D and polymeric systems, respectively.  The HSE06 hybrid functional49-52 was used to treat exchange-correlation, given its success in remediating issues with the underestimation of 146  band gaps by semilocal functionals. Additionally, in an attempt to capture the inter-chain hydrogen bonding in the polymeric structure, the dispersion-corrected DFT-D2 method was applied.53 6.3 Three-electron bond formation in melam/melem dimer cations  The depth and character of positively charged polarons in g-C3N4 has been assessed by studying nuclear relaxation in the cations of the molecular models described above.  The PESs of melam and melem dimer as a function of N1-N1 distance, with all other nuclear coordinates relaxed at the DFT level of theory, are presented in Figure 6.3A and B, respectively.  Note that in both cases the cation minimum structure occurs at a very short N1-N1 distance (2.16 and 2.13 Å for melam and the melem dimer, respectively) compared to the N1-N1 distance in the neutral minimum structure (2.89 and 2.85 Å; see Figure 6.4A and B for optimized neutral and cationic structures).  This relaxation is very energetically favorable in both systems, as indicated by the total cation relaxation energies (Table 6.1), which are 649 meV and 422 meV, respectively.  It is noteworthy that neither melam nor melem dimer have a planar structure in the cationic or neutral state, with out-of-plane distortion of the rings relative to one another present in all structures.  (Optimized geometries, absolute energies, and ionization potentials are presented in Supporting Information.) The hole orbitals at selected points along the relaxed PESs are presented as insets in Figure 6.3A and B, with similar features observed in each.  At larger N1-N1 distance corresponding to the neutral minimum structure, ionization occurs from a nitrogen-n -bonding (n) character in melam.  As the distance between position 1 nitrogen atoms decreases towards the cationic minimum there is a change in the character of the hole orbitals; they lose their pair n orbitals.  For conciseness, we term the hole left behind after removal of an electron from an 147   orbitals seen at short N1-N1 distances correspond to out-of-phase (OOP) superpositions of lone pair orbitals on the N1 atoms of neighboring monomers.  Such orbitals are hereafter referred to as nOOP orbitals (see Figure 6.5 top center).    Table 6.1 Cationic nuclear relaxation data for models of g-CN Cation nuclear relaxation data for all model systems, including uncorrected and corrected cation relaxation energies.   Cation N1-N1 distance relaxation (Å) Cation minimum N1-N1 bond order  (meV)  (meV) Melam  0.73 0.32 649 531 Melem dimer 0.72 0.34 422 304 2-D melem trimer 0.47 0.34 362 244 Polymeric melem trimer 0.72 0.33 330 212 2-D melem hexamer 0.49 0.31 276 158  148   Figure 6.3 Cationic melam and melem dimer PESs Relaxed cation PESs for melam (A) at the U-CAM-B3LYP and IP-EOM-CCSD(2h,1p) levels of theory and melem dimer (B; U-CAM-B3LYP only) as a function of N1-N1 distance.   The bond order for both molecules as a function of the same distance is presented in C.  The U-CAM-B3LYP hole orbitals at the N1-N1 distance are shown as insets in A and B. 149   Figure 6.4 Minimum energy structures for neutral and cationic models of g-CN Minimum energy structures for the neutral (left) and cation (right) electronic states of the molecular models.  Shown on each structure are the distances between position 1 nitrogen atoms.  Between each pair of structures is the change in energy and N1-N1 distance associated with relaxation from the neutral to the cation minimum.    150  The change in character of the hole upon shortening of the N1-N1 distance can be intuitively understood in terms of the interaction of molecular orbitals on neighboring monomer units, as illustrated in the molecular orbital diagram in Figure 6.5.  In isolated melem, the highest occupied preferred.  However, when the N1-N1 distance is shortened, strong splitting of interacting N1 lone pair orbitals into an in-phase superposition of lone pair orbitals (nIP) and the nOOP orbital is unfavorable nOOP orbital becoming the HOMO.  At short N1-N1 distance both melam and melem dimer cations possess a single electron occupying the nOOP orbital and two electrons in the corresponding nIP orbital.  Noting the bonding character of the nIP orbital and the antibonding character of the nOOP orbital, this electron configuration indicates the presence of a two-centered, three-electron bond between position 1 nitrogen atoms.  The formation of this bond can also be observed in the computed bond order, which rises from less than 0.10 at N1-N1 distances similar to the neutral minimum structure, to 0.32 at the cation minimum (Figure 6.3C; compare to an ideal value of 0.5 for a three-electron bond).  Two-centered, three-electron N-N bonds have been observed in other heteroaromatic compounds,54-56 but to the knowledge of the authors this is the first prediction of polaronic behavior due to intramolecular three-electron bond formation.  Both PESs exhibit a double-well-like shape, indicating the crossing of states corresponding to ionization out of an n -well-like character is more pronounced in melem dimer than in melam, and the cation relaxation energy is smaller (422 meV vs. 649 meV).  hole in the larger system relative to the n hole, the energy of which is less sensitive to the system size. 151   Figure 6.5 Molecular orbital diagram for Pauli-repulsive n-type orbitals in g-CN Model molecular orbital diagrams for melem dimer at different N1-N1 distances.   At shorter distances, splitting of the in-phase and out-of-phase position 1 nitrogen lone pair orbitals is considerably  to n character.  The diagram is plotted using orbital eigenvalues obtained from CAM-B3LYP/cc-pVDZ calculations of neutral melem and melem dimer molecules.  Shown on the right are pictoral examples of the dimer molecular orbitals corresponding to the in-phase and out-of-phase superpositions of the highest occupied molecular orbitals of melem at an N1-N1 distance of 2.89 Å.  The correspondence of energy levels (left) to orbitals (right) is indicated by color.  6.4 Three-electron bond formation in larger model cations Larger models of g-C3N4 were investigated to address several questions:  In a system containing multiple N1-N1 contacts, does the positive charge preferentially localize at a single such contact?  Is three-electron bond formation possible in the more structurally constrained 2-D phase of g-C3N4, or only in the more flexible polymeric phase?  If formation is possible in both phases, how does the cation relaxation energy depend on the phase?  Indeed, the hole is found to localize at a single N1-N1 contact in the polymeric melem trimer cation, which contains two such contacts.  The trimer exhibits similar nuclear relaxation to that of the melem dimer (Figure 6.4C), with a decrease in N1-N1 distance of 0.72 Å and N1-N1 bond order of 0.33 at the cation minimum, signaling the formation of a three-electron bond at a single N1-N1 152  contact.  The second N1-N1 pair on the polymeric trimer exhibits a negligible decrease in distance (0.01 Å) at the optimized cation structure relative to the neutral.   Relaxation from the neutral mthe melem dimer.  (Images of hole orbitals are presented in Figures A5.2 and A5.3.)  The cation relaxation energy for the polymeric melem trimer, 330 meV (Table 6.1), is 92 meV less than that of the melem dimer, despite an identical change in N1-N1 distance and a similar N1-N1 bond order predicted for larger N1-N1 distances by delocalization of the hole onto the third melem unit.  Three-electron bond formation is also observed in models of the more structurally constrained 2-D phase of g-C3N4.  The N1-N1 relaxation in the 2-D melem trimer cation is again confined to a single N1-N1 pair, with a decrease in N1-N1 distance of 0.47 Å (Figure 6.4D) from the neutral to the cation minimum structure.  The cation relaxation energy for the 2-D melem trimer is 362 meV and the cation minimum structure has an N1-N1 bond order of 0.34.    It is instructive to compare cation relaxation energies and geometric distortions associated with three-electron bond formation in the 2-D and polymeric melem trimers.  Though the N1-N1 distance decrease in the polymeric trimer is greater than that in the 2-D trimer by 0.25 Å and the cation minimum N1-N1 bond orders are almost equivalent, the relaxation energy of the polymeric trimer is actually 32 meV less than that of the 2-D trimer.  This may be attributed to the fact that the 2-D trimer has a buckled structure at the neutral minimum geometry (Figure 6.4D), caused by Pauli repulsion of N1 lone pairs (consistent with DFT calculations showing buckling of extended 2-D sheets20-22).    The formation of a cation alleviates some buckling of the structure as removal of a lone pair electron replaces Pauli repulsion at one N1-N1 contact with a weak covalent bond.  153  Thus the reduction of strain due to lone pair repulsion in the neutral minimum structure of the 2-D trimer may explain the greater cation relaxation energy compared to the polymeric trimer.    The shorter N1-N1 distance in the neutral 2-D trimer also affects the character of the hole observed upon vertical ionization.  The equilibrium N1-N1 distance of the neutral 2-D trimer (2.64 this N1-N1 distance the splitting of nIP and nOOP orbitals, in combination with the buckled shape of the molecule, results in a hole olocalized nOOP form as was observed in the polymeric models.  (Images of hole orbitals are presented in Figures A5.4 and A5.5.)  These results suggest that the 2-D phase of the material may favor n holes even before nuclear relaxation, while the polymeric phase requires such relaxation for the n hole to be favored.  As will be discussed in Section 6.7, this is consistent with the in periodic DFT calculations of fully condensed 2-D g-C3N4.20  The 2-D melem hexamer exhibits a buckled neutral minimum structure similar to that of the 2-D melem trimer.  Much like the 2-D melem trimer the hexamer character at the neutral minimum structure and n character at the cation minimum.  (Images of hole orbitals are presented in Figures A5.6 and A5.7.)  In relaxing to the cation minimum the N1-N1 distance of a single N1-N1 pair decreases by 0.49 Å (Figure 6.5E), with an energy decrease of 276 meV.  This cation relaxation energy is 86 meV lower than that of the 2-D trimer, despite a similar decrease in N1-N1 distance.    Though these larger model calculations suggest that three-electron bond formation occurs even in large models, there is a significant decline in the cation relaxation energy with increasing system size, as seen in Figure 6.6.  The reason for the further decrease in the relaxation energy of 154  the 2-D hexamer relative to the 2-D trimer is not completely clear.  The fact that three-electron bond formation reduces buckling in the surrounding trimer may play a role.  In the hexamer, this trimer is bound to additional melem units which buckle out-of-plane as well, and it is possible that the planarization of one unit strains the surrounding non-planar structure.  In fact, it is observed that the total of all N1-N1 distances in the melem hexamer excluding the one participating in the three-electron bond is decreased by 0.10 Å upon relaxation to the cation minimum, suggesting the existence of such strain.  Figure 6.6 Cationic relaxation energies and nitrogen atoms in models of g-CN Cation relaxation energies in the model systems as a function of the number of nitrogen atoms in the molecule.   Though our study has been limited to the polymeric and 2-D phases of g-C3N4, three-dimensional structures with reduced N1-N1 repulsion have been suggested in previous theoretical studies.21  Though the N1-N1 lone pair repulsion is reduced in the proposed 3-D structures, it is noteworthy that close contact between N1-N1 pairs still exists in this phase, and thus three-electron bond formation is likely possible. 155   The role of the reported polaronic relaxation in determining photocatalytic efficiency may be complex.  The suggestion that polaronic relaxation involves the formation of an n hole is interesting in light of a recent experimental study reporting that photocatalytic efficiency for spectrum of g-C3N4 materials synthesized under different conditions.24  The authors suggest that a trap state for either electrons or holes may be responsible for this behavior, and the trapping of -electron bonds provides a plausible explanation for this observation.  On the other hand, though polarons clearly slow charge transport, it is possible that the localization of charge in strained three-electron bonds on the surface of the material may favor efficient photocatalysis of oxidative reactions by providing a reactive and highly localized source of positive charge.  6.5 Error assessment for DFT energies  To assess the error in the U-CAM-B3LYP cation relaxation energies described above, IP-EOM-CCSD(2h,1p) calculations were performed at the DFT-optimized structures of the smallest model, melam, to provide a more accurate description of the relaxation process.  Note that the IP-EOM PES for relaxation on the cationic surface in Figure 6.3A is nearly identical to the U-CAM-B3LYP curve.  In addition, IP-EOM predicts a cation relaxation energy of 531 meV, compared with 649 meV obtained with U-CAM-B3LYP.  The small (118 meV) difference stems largely from disagreement in the ionization potential at the neutral minimum geometry, where U-CAM-B3LYP predicts an ionization potential of 8.179 eV, while IP-EOM predicts 8.080 eV, for a difference of 99 meV (all calculated IPs are available in Table A6.1).  A U-CAM-B3LYP treatment of melam thus overestimates the ionization potential in this odd-electron intramolecular bonded 156  system when compared with IP-EOM, as has been seen in previous DFT studies of intermolecular three-electron bonds.57-59  The difference between the cation relaxation energies obtained with each method, +118 meV, is used as an estimate for the error in the U-CAM-B3LYP cation relaxation energies of larger model systems and corrected relaxation energies are reported in conjunction with uncorrected U-CAM-B3LYP energies in Table 6.1.  The impact of this IP-EOM correction on the cation relaxation energies is significant in the case of the polymeric trimer and the 2-D melem trimer/hexamer, with the corrected cation relaxation energy of the 2-D hexamer being only 158 meV.  However, bond formation is still favorable in all models.  6.6 Charge carrier mobility  Calculation of the inner sphere reorganization energy provides an estimate of the mobility of holes in g-C3N4 following nuclear relaxation and three-electron bond formation.  Table 6.2 lists the inner sphere reorganization energies for all model systems, as well as their cation and neutral relaxation energy contributions.   The cation relaxation energies of these models, as previously discussed, decrease as the system size increases.  The neutral relaxation energy, on the other hand, correlates with the phase of the material2-D or polymeric.  The reorganization energies of the 2-D phase models are more than 0.5 eV lower than those in the polymeric models, with the 2-D phase models exhibiting neutral relaxation energies of 1.13 and 1.05 eV for the trimer and hexamer, respectively, while the polymeric models exhibit neutral relaxation energies of 1.65 eV for the melem dimer and 1.68 eV for the polymeric melem trimer.    Strong Pauli repulsion forces between lone pairs in the neutral molecules at the cation minimum structures are responsible for the larger energy associated with relaxation of the neutral molecule relative to the cation, which can be seen in the increased slope of the neutral melem PES 157  at short N1-N1 distance (see Figure A6.1). The differences between the neutral relaxation energies of the 2-D and polymeric models can be rationalized in terms of the N1-N1 distance.  In both polymeric phase models (melem dimer and polymeric trimer) the N1-N1 distances in the cation and neutral minimum structures differ by 0.72 Å (Figure 6.4, Table 6.1).  In the 2-D models the relaxation of the N1-N1 distance is reduced0.47 Å for the trimer and 0.49 Å for the hexamerthus less of the lone pair repulsion is relieved in these systems.  Table 6.2 Inner sphere component of the polaronic reorganization energy Values for the inner sphere reorganization energy, as well as its neutral and cation relaxation energy contributions, in all model systems.   cat (eV) neut (eV)  Melam 0.649 1.502 2.151 Melem dimer 0.422 1.649 2.071 2-D melem trimer 0.362 1.125 1.487 Polymeric melem trimer 0.330 1.676 2.006 2-D melem hexamer 0.276 1.054 1.329   The inner sphere reorganization energies of all models are large (Table 6.2) and suggest that transport may be slow in both the 2-D and polymeric materials.  The reorganization energies of the polymeric systems are similar to those predicted for iron oxide (2.1 eV),60 another photocatalytic material with performance limited by a very low hole collection distance.61  This result suggests a larger hole mobility in more condensed g-C3N4 materials, and also suggests that doping schemes which disrupt the N1-N1 contacts may increase charge carrier mobility and therefore photocatalytic efficiency.  Though far from the only possible explanation, this is consistent with the observation of increased photocatalytic efficiency in carbon- and sulfur-doped g-C3N4 materials.62,63 158   We report inner sphere reorganization energies here to allow the rough estimation of hole transport rates, however the grossly anharmonic shape of the melem dimer cation PES (Figure 6.3B) suggests that the use of Marcus theory for quantitative estimation of the transfer rate may be inadvisable.  6.7 Lone pair interactions in extended systems  Though there are earlier reports of the electronic band structure of extended graphitic carbon nitrides,1,20,25,26 here we discuss the role that N1-N1 lone pair interactions play in its determination.  This has been investigated by application of periodic DFT to idealized planar 2-D and polymeric g-C3N4, with results presented in Figures 6.7 and 8.  The polymeric and 2-D structures differ in that their N1-N1 distances are predicted to be 2.70 Å and 2.43 Å, respectively, and as suggested by the previously discussed model studies, the chemical character of the valence band maximum correlates with this distance.  Figures 6.7C and D report the atomic character resolved band structures of the two phases of g-C3N4z) and x, and py) character (where z is the direction perpendicular to the graphitic plane).  Consistent with previous theoretical reports,1,20 the valence band maxima of the 2-D and polymeric structures  character, respectively (Figure 6.7).  Further, looking directly at the Kohn-Sham orbitals pictured in Figure 6.8, it is evident that the valence band maximum of the 2-D system corresponds to an out-of-phase (antibonding) superposition of nitrogen lone pair orbitals, comparable to the nOOP orbital described above.  This change in valence band character between the polymeric and 2-D phases can be explained in terms of the increased splitting between n orbitals with decreasing N1-N1 distance.  Aside from the character of the near-gap orbitals, there is a 0.86 eV difference in the band gap between the 2-D (2.78 eV) and polymeric (3.64 eV) phases 159  which are indirect and direct, respectively. While it has been previously reported for the 2-D phase,20,25,26 it is worth remarking that the DFT gap of the polymeric phase is opened considerably (~1.0 eV relative to Reference 1) through the use of a hybrid functional.   Figure 6.7 Band structures for the 2-D and polymeric phases of g-CN The 2-D (A) and polymeric (B) planar bulk structures.  The 2-D lattice constant is 7.08 Å, the polymeric lattice constants are 17.10 Å and 12.80 Å.  Atomic character resolved Kohn-Sham band structures of the 2-D (C) and polymeric (D) periodic systems are also shown.    160   Figure 6.8 Kohn-Sham orbitals for the valence and conduction bands of g-CN Real part of the Kohn-Sham orbitals at the -point for the 2-D valence (A) and conduction bands (B) as well as the polymeric valence (C) and conduction bands (D) in each periodic system.   6.8 Conclusions  Taken together, these model molecular and periodic DFT studies of g-C3N4 suggest several important features of hole transport related to the interaction between lone pairs on adjacent melem monomers.  The periodic DFT calculations suggest that the character of the valence band maximum in idealized planar 2-D and polymeric g-C3N4 are correlated with the distance between N1 upon vertical ionization, while the planar 2-D material exhibits an n hole corresponding to an out-of-phase superposition of nitrogen lone pair orbitals.  The model calculations suggest that, upon creation of a hole, strong nuclear relaxation occurs along the N1-N1 distance coordinate.  Regardless of the phase of the material, the formation of a three-electron bond between N1 atoms on adjacent monomers stabilizes the cation significantly.  The Marcus reorganization energy for hole self-exchange between N1-N1 contacts 161  in polymeric g-C3N4 is predicted to be approximately 2.1 eV, while it is significantly smaller (~1.4 eV) in the 2-D phase due to the more constrained structure of the fully condensed material.  These large reorganization energies suggest slow transport in both phases, which may be alleviated by modification of the material (for example by substitutional doping of the position 1 nitrogen atoms62,63).  In addition to the role relaxation plays in slowing charge transport, the strained and local nature of the three-electron bond suggests that this may be the active species in oxidative reactions catalyzed by g-C3N4.    162             APPENDICES    163  APPENDIX A: Supplementary content   Figure A6.1 Neutral melam PES   Table A6.1 Vertical and adiabatic ionization potentials for models of g-CN    Vertical Ionization Potential (eV) Adiabatic Ionization Potential (eV) Melamine dimer 8.179 7.530 Melem dimer 8.004 7.582 Polymeric melem trimer 7.909 7.578 2-D melem trimer 7.787 7.425 2-D melem hexamer 7.662 7.386   164   Figure A6.2 Polymeric melem trimer cationic LUMO (at neutral minimum)   Figure A6.3 Polymeric melem trimer cationic LUMO (at cationic minimum)  Figure A6.4 2-D melem trimer cationic LUMO (at neutral minimum)  165   Figure A6.5 2-D melem trimer cationic LUMO at cationic minimum    Figure A6.6 (2-D melem hexamer cationic LUMO at neutral minimum 166   Figure A6.7 (2-D melem hexamer cationic LUMO (at cationic minimum)   167  Table A6.2 Neutral melam vibrational frequencies  Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 21.74 0.00363 1000.98 0.00970 39.57 0.06073 1011.31 0.02138 79.84 0.00454 1015.15 0.06627 163.85 0.02063 1029.16 0.81220 187.85 0.01275 1078.02 0.00443 190.98 0.00442 1103.66 0.16557 194.56 0.29539 1149.98 0.01065 199.07 0.11663 1181.23 0.06639 224.94 0.06284 1187.32 0.06904 255.76 0.10550 1219.43 1.22349 257.68 2.51809 1317.71 2.13015 273.69 4.16530 1341.20 0.09345 285.65 5.27366 1392.92 4.34015 305.15 11.24085 1479.83 2.81761 309.20 2.23226 1485.83 29.99354 340.59 0.27255 1497.72 0.28346 364.18 0.01773 1507.01 0.11881 475.42 0.15021 1511.53 1.94709 510.04 0.00127 1560.64 45.64420 546.56 0.09057 1615.41 1.62765 548.44 0.27260 1616.19 13.81858 557.00 0.00590 1646.66 18.19809 569.49 0.04684 1661.86 9.91719 586.61 0.06521 1673.21 11.15572 594.08 0.20577 1684.20 3.15297 618.21 0.39302 1700.13 0.04076 639.03 0.04912 1703.03 29.34736 701.73 1.19530 3616.31 1.83705 758.19 0.06051 3616.74 2.99781 758.85 0.24789 3617.82 1.36597 765.79 0.09908 3618.92 1.55038 766.14 0.01627 3637.65 1.54808 770.38 0.00904 3755.80 1.18527 795.37 0.27192 3756.16 1.23568 852.81 0.01148 3759.40 2.78654 854.36 2.66032 3760.24 0.06095 168   Figure A6.8 Melam neutral minimum vibrational frequencies   169  Table A6.3 Cationic melam vibrational frequencies  Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 40.08 0.06760 964.09 0.91750 47.27 0.00473 993.23 0.09330 79.13 0.00894 1007.1 0.23077 128.79 0.00831 1008.01 0.60855 156.76 0.15335 1059.43 0.22823 181.43 0.00906 1088.17 0.05664 198.46 0.50981 1175.11 0.01301 201.76 0.02238 1176.12 0.01331 232.08 0.07373 1194.82 0.52863 281.68 0.25354 1204.28 0.00467 284.03 0.14126 1362.09 0.14797 338.41 1.61224 1376.36 0.02084 355.85 3.38656 1387.56 0.00772 363.8 0.25608 1437.66 1.84117 401.57 2.51550 1515.66 0.84816 404.37 6.46164 1523.58 2.97140 419.71 3.30948 1540.06 1.71404 496.5 0.30778 1540.62 0.02432 527.1 0.76881 1543.49 0.08814 544.41 1.77449 1592.16 20.55611 554.92 0.58065 1633.64 0.09302 560.59 0.20661 1642.89 4.00481 583.32 0.98741 1659.57 22.66425 596.65 1.68706 1683.37 60.13834 603.12 0.17130 1707.78 15.99771 615.18 0.00247 1718.41 3.71891 615.87 0.00836 1723.45 14.74345 630.82 0.10076 3600.16 0.92132 722.62 0.00177 3600.91 7.24440 727.87 0.75581 3605.39 12.56668 739.06 0.01322 3607.08 0.32886 740.32 1.11628 3639.38 3.70504 751.26 0.26318 3735.5 3.72736 794.27 0.21346 3736.63 0.05623 817.38 0.01643 3747.18 2.13963 820.37 2.49736 3747.43 3.05878 170   Figure A6.9 Melam cationic minimum vibrational frequencies   171  Table A6.4 Neutral melem dimer vibrational frequencies Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 18.95 0.03303 649.20 0.00084 1335.48 0.70916 40.92 0.01132 659.10 0.09602 1336.99 24.70108 75.11 0.00211 695.33 0.00093 1345.94 0.60251 95.38 0.00650 725.80 0.62182 1380.62 6.14004 99.05 0.00074 741.27 0.00750 1408.27 0.24554 101.79 0.17186 745.08 0.01405 1453.29 23.07768 104.50 0.00508 746.36 0.00724 1491.54 56.88236 145.27 0.20199 751.85 0.19840 1497.48 4.86312 166.30 0.07044 760.79 0.01469 1508.29 2.54074 182.16 0.01638 763.12 0.01117 1515.40 100.03272 223.73 0.06500 764.65 0.00077 1561.53 30.50813 225.31 0.19090 772.51 0.95138 1564.49 0.66912 227.90 2.29348 777.57 0.00497 1576.95 0.48485 234.89 1.75951 794.19 0.00477 1590.71 1.74551 235.60 1.38815 798.23 0.00433 1618.37 16.85308 239.29 5.21557 798.46 0.00430 1621.60 1.62274 245.40 0.00777 802.57 0.08950 1634.32 31.86102 250.24 0.00287 804.66 0.06798 1644.10 81.51322 251.30 0.00024 845.11 0.10084 1659.60 0.18144 265.30 6.15931 854.71 0.03481 1665.50 5.00914 267.13 0.77381 856.10 3.43631 1698.45 2.78729 292.81 0.08222 884.33 0.15338 1700.19 1.32223 298.39 0.73555 1008.51 0.00125 1736.57 19.82190 299.83 0.01741 1029.68 0.00143 1744.22 17.88364 360.99 0.05817 1031.78 0.02164 1746.53 48.60511 381.66 0.00240 1051.69 0.05295 1754.08 5.89168 454.91 0.00015 1071.43 0.06096 3619.18 2.53073 456.50 0.00150 1075.27 0.08621 3619.73 9.04812 470.87 0.02914 1080.61 0.38719 3620.86 4.94841 489.68 0.01133 1084.88 0.99850 3621.23 1.73667 533.99 0.04011 1176.77 0.00339 3622.54 2.04733 566.31 0.00184 1191.65 1.77676 3769.02 1.99356 590.44 0.01631 1201.51 0.00691 3770.03 2.00471 591.99 0.01100 1203.32 0.09658 3770.57 0.87775 602.79 0.15665 1209.21 0.08934 3770.73 2.68380 610.78 0.00084 1220.57 0.07206 613.44 0.00023 1255.73 4.52440 616.74 0.00105 1261.72 0.39457 619.28 0.00005 1314.04 6.12683 172   Figure A6.10 Melem dimer neutral minimum vibrational frequencies   173  Table A6.5 Cationic melem dimer vibrational frequencies  Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 18.23 0.00030 642.50 0.00164 1316.55 5.27361 20.27 0.02151 645.70 0.00402 1342.33 0.36828 41.27 0.00132 654.56 0.47542 1347.93 0.42934 65.94 0.00056 656.98 0.00707 1370.13 1.06069 86.06 0.07800 689.92 0.00380 1385.23 0.91744 101.49 0.04182 716.02 0.00316 1430.04 0.28132 101.72 0.00263 719.64 0.72133 1455.31 13.57349 102.67 0.35986 722.04 0.00236 1469.58 10.54304 153.41 0.08760 730.58 0.57605 1514.32 3.91612 172.60 0.04656 748.42 0.41011 1518.72 4.46637 198.63 0.00712 748.67 0.05340 1524.90 2.32739 210.99 0.01149 754.19 0.00365 1553.68 1.75628 221.74 0.36825 760.76 0.52685 1571.70 1.27637 237.14 0.00090 770.23 0.00300 1572.25 3.60928 243.66 0.02021 773.89 0.90894 1588.47 67.44387 244.03 0.35882 784.03 0.00189 1602.74 13.05370 246.37 0.17476 794.10 0.94541 1618.61 9.89786 260.28 0.00371 797.22 0.00115 1657.27 110.37583 293.97 0.06885 797.86 0.01207 1666.79 25.04645 300.27 0.10309 813.38 0.08627 1674.80 2.44874 305.48 0.10534 838.81 0.07735 1683.71 41.22708 375.70 0.02984 839.39 3.26887 1719.79 32.70653 386.11 0.02949 880.62 0.32556 1723.47 1.66919 398.14 0.00304 1011.38 0.03601 1727.67 0.62832 410.20 7.60588 1015.55 0.10729 1758.44 1.14641 441.99 3.85491 1029.00 0.01802 1760.83 19.24052 444.15 0.62468 1030.32 0.01931 1769.47 49.85151 444.62 5.98246 1072.48 0.47710 3598.06 13.46771 458.69 0.00447 1075.29 0.34587 3598.55 4.69208 464.29 0.01351 1077.44 0.16748 3604.39 0.89295 493.85 0.04244 1078.51 0.39553 3606.35 9.22201 543.10 0.02963 1168.12 0.53976 3627.04 4.25060 554.23 1.09595 1178.18 0.02554 3742.68 2.35925 583.92 0.01908 1183.83 0.04523 3742.89 2.73678 588.82 0.00886 1184.81 0.00978 3749.94 4.67540 592.15 0.80479 1223.55 0.37311 3751.90 0.23494 629.61 0.00788 1231.95 0.01059 633.17 1.88978 1289.37 0.23295 634.45 0.30689 1309.11 0.01062  174   Figure A6.11 Melem dimer cationic minimum vibrational frequencies    175  APPENDIX B: Molecular geometries for models of g-CN                                              176                                                177                                                178                                                179                                                180                                                181                                                182                                                183                                                184                                                185                                                186                                                187                                                188             REFERENCES    189  REFERENCES   (1) Wang, X.; Maeda, K.; Thomas, A.; Takanabe, K.; Xin, G.; Carlsson, J. M.; Domen, K.; Antonietti, M. A metal-free polymeric photocatalyst for hydrogen production from water under visible light. Nature Materials 2009, 8, 76-80.   (2) Li, Q. Y.; Yue, B.; Iwai, H.; Kako, T.; Ye, J. H. Carbon Nitride Polymers Sensitized with N-Doped Tantalic Acid for Visible Light-Induced Photocatalytic Hydrogen Evolution. J. Phys. Chem. C 2010, 114, 4100-4105.   (3) Wang, X.; Blechert, S.; Antonietti, M. Polymeric Graphitic Carbon Nitride for Heterogeneous Photocatalysis. Acs Catalysis 2012, 2, 1596-1606.   (4) Su, F. Z.; Mathew, S. C.; Lipner, G.; Fu, X. Z.; Antonietti, M.; Blechert, S.; Wang, X. C. mpg-C3N4-Catalyzed Selective Oxidation of Alcohols Using O-2 and Visible Light. J. Am. Chem. Soc. 2010, 132, 16299-16301.   (5) Wang, Y.; Zhang, J. S.; Wang, X. C.; Antonietti, M.; Li, H. R. Boron- and Fluorine-Containing Mesoporous Carbon Nitride Polymers: Metal-Free Catalysts for Cyclohexane Oxidation. Angewandte Chemie-International Edition 2010, 49, 3356-3359.   (6) Su, F. Z.; Mathew, S. C.; Mohlmann, L.; Antonietti, M.; Wang, X. C.; Blechert, S. Aerobic Oxidative Coupling of Amines by Carbon Nitride Photocatalysis with Visible Light. Angew. Chem. Int. Ed. 2011, 50, 657-660.   (7) Wang, Y.; Li, H. R.; Yao, J.; Wang, X. C.; Antonietti, M. Synthesis of boron doped polymeric carbon nitride solids and their use as metal-free catalysts for aliphatic C-H bond oxidation. Chemical Science 2011, 2, 446-450.   (8) Yan, S. C.; Li, Z. S.; Zou, Z. G. Photodegradation Performance of g-C3N4 Fabricated by Directly Heating Melamine. Langmuir 2009, 25, 10397-10401.   (9) Yan, S. C.; Li, Z. S.; Zou, Z. G. Photodegradation of Rhodamine B and Methyl Orange over Boron-Doped g-C3N4 under Visible Light Irradiation. Langmuir 2010, 26, 3894-3901.   (10) Teter, D. M.; Hemley, R. J. Low-compressibility carbon nitrides. Science 1996, 271, 53-55.   (11) Hohenberg, P.; Kohn, W. INHOMOGENEOUS ELECTRON GAS. Phys. Rev. B 1964, 136, B864-&.   (12) Kohn, W.; Sham, L. J. SELF-CONSISTENT EQUATIONS INCLUDING EXCHANGE AND CORRELATION EFFECTS. Physical Review 1965, 140, 1133-&. 190   (13) Liu, A. Y.; Wentzcovitch, R. M. STABILITY OF CARBON NITRIDE SOLIDS. Physical Review B 1994, 50, 10362-10365.   (14) Ortega, J.; Sankey, O. F. RELATIVE STABILITY OF HEXAGONAL AND PLANAR STRUCTURES OF HYPOTHETICAL C3N4 SOLIDS. Physical Review B 1995, 51, 2624-2627.   (15) Nesting, D. C.; Badding, J. V. High-pressure synthesis of sp(2)-bonded carbon nitrides. Chemistry of Materials 1996, 8, 1535-1539.   (16) Lowther, J. E. Relative stability of some possible phases of graphitic carbon nitride. Physical Review B 1999, 59, 11683-11686.   (17) Lotsch, B. V.; Doeblinger, M.; Sehnert, J.; Seyfarth, L.; Senker, J.; Oeckler, O.; Schnick, W. Unmasking melon by a complementary approach employing electron diffraction, solid-state NMR spectroscopy, and theoretical calculations-structural characterization of a carbon nitride polymer. Chemistry-a European Journal 2007, 13, 4969-4980.   (18) Bojdys, M. J.; Mueller, J.-O.; Antonietti, M.; Thomas, A. Ionothermal Synthesis of Crystalline, Condensed, Graphitic Carbon Nitride. Chemistry-a European Journal 2008, 14, 8177-8182.   (19) Wang, Y.; Wang, X.; Antonietti, M. Polymeric Graphitic Carbon Nitride as a Heterogeneous Organocatalyst: From Photochemistry to Multipurpose Catalysis to Sustainable Chemistry. Angewandte Chemie-International Edition 2012, 51, 68-89.   (20) Deifallah, M.; McMillan, P. F.; Cora, F. Electronic and structural properties of two-dimensional carbon nitride graphenes. J. Phys. Chem. C 2008, 112, 5447-5453.   (21) Gracia, J.; Kroll, P. Corrugated layered heptazine-based carbon nitride: the lowest energy modifications of C3N4 ground state. Journal of Materials Chemistry 2009, 19, 3013-3019.   (22) Sehnert, J.; Baerwinkel, K.; Senker, J. Ab initio calculation of solid-state NMR spectra for different triazine and heptazine based structure proposals of g-C3N4. Journal of Physical Chemistry B 2007, 111, 10671-10680.   (23) Thomas, A.; Fischer, A.; Goettmann, F.; Antonietti, M.; Muller, J. O.; Schlogl, R.; Carlsson, J. M. Graphitic carbon nitride materials: variation of structure and morphology and their use as metal-free catalysts. Journal of Materials Chemistry 2008, 18, 4893-4908.   (24) Jorge, A. B.; Martin, D. J.; Dhanoa, M. T. S.; Rahman, A. S.; Makwana, N.; Tang, J. W.; Sella, A.; Cora, F.; Firth, S.; Darr, J. A.; McMillan, P. F. H-2 and O-2 Evolution from Water Half-Splitting Reactions by Graphitic Carbon Nitride Materials. J. Phys. Chem. C 2013, 117, 7178-7185.  191   (25) Du, A. J.; Sanvito, S.; Li, Z.; Wang, D. W.; Jiao, Y.; Liao, T.; Sun, Q.; Ng, Y. H.; Zhu, Z. H.; Amal, R.; Smith, S. C. Hybrid Graphene and Graphitic Carbon Nitride Nanocomposite: Gap Opening, Electron-Hole Puddle, Interfacial Charge Transfer, and Enhanced Visible Light Response. J. Am. Chem. Soc. 2012, 134, 4393-4397.   (26) Wei, W.; Jacob, T. Strong excitonic effects in the optical properties of graphitic carbon nitride g-C3N4 from first principles. Phys. Rev. B 2013, 87.   (27) Jurgens, B.; Irran, E.; Senker, J.; Kroll, P.; Muller, H.; Schnick, W. Melem (2,5,8-triamino-tri-s-triazine), an important intermediate during condensation of melamine rings to graphitic carbon nitride: Synthesis, structure determination by X-ray powder diffractometry, solid-state NMR, and theoretical studies. Journal of the American Chemical Society 2003, 125, 10288-10300.   (28) Zheng, W.; Wong, N.-B.; Li, W.-K.; Tian, A. Absorption spectra of tri-s-triazines: time dependent density functional theory calculations. New Journal of Chemistry 2006, 30, 1307-1318.   (29) Coropceanu, V.; Cornil, J.; da Silva, D. A.; Olivier, Y.; Silbey, R.; Bredas, J. L. Charge transport in organic semiconductors. Chem. Rev. 2007, 107, 926-952.   (30) Ufimtsev, I. S.; Martinez, T. J. Quantum chemistry on graphical processing units. 1. Strategies for two-electron integral evaluation. J. Chem. Theory Comput. 2008, 4, 222-231.   (31) Ufimtsev, I. S.; Martinez, T. J. Quantum Chemistry on Graphical Processing Units. 2. Direct Self-Consistent-Field Implementation. J. Chem. Theory Comput. 2009, 5, 1004-1015.   (32) Ufimtsev, I. S.; Martinez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients, Geometry Optimization, and First Principles Molecular Dynamics. Journal of Chemical Theory and Computation 2009, 5, 2619-2628.   (33) Kaestner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. Journal of Physical Chemistry A 2009, 113, 11856-11865.   (34) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chemical Physics Letters 2004, 393, 51-57.   (35) Dunning, T. H. GAUSSIAN-BASIS SETS FOR USE IN CORRELATED MOLECULAR CALCULATIONS .1. THE ATOMS BORON THROUGH NEON AND HYDROGEN. Journal of Chemical Physics 1989, 90, 1007-1023.   (36) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; 192  Montgomery, J. A. GENERAL ATOMIC AND MOLECULAR ELECTRONIC-STRUCTURE SYSTEM. Journal of Computational Chemistry 1993, 14, 1347-1363.   (37) Gordon, M. S.; Schmidt, M. W. Advances in electronic structure theory: GAMESS a decade later. Theory and Applications of Computational Chemistry: the First Forty Years 2005, 1167-1189.   (38) Nooijen, M.; Snijders, J. G. COUPLED CLUSTER APPROACH TO THE SINGLE-PARTICLE GREEN-FUNCTION. Int. J. Quantum Chem 1992, 55-83.   (39) Nooijen, M.; Snijders, J. G. COUPLED-CLUSTER GREEN-FUNCTION METHOD - WORKING EQUATIONS AND APPLICATIONS. Int. J. Quantum Chem 1993, 48, 15-48.   (40) Stanton, J. F.; Gauss, J. ANALYTIC ENERGY DERIVATIVES FOR IONIZED STATES DESCRIBED BY THE EQUATION-OF-MOTION COUPLED-CLUSTER METHOD. Journal of Chemical Physics 1994, 101, 8938-8944.   (41) Bartlett, R. J.; Stanton, J. F.: Applications of Post-HartreeFock Methods: A Tutorial. In Rev. Comput. Chem.; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH: New York, 1994; Vol. 5; pp 65-169.   (42) Gour, J. R.; Piecuch, P.; Wloch, M. Active-space equation-of-motion coupled-cluster methods for excited states of radicals and other open-shell systems: EA-EOMCCSDt and IP-EOMCCSDt. Journal of Chemical Physics 2005, 123.   (43) Gour, J. R.; Piecuch, P. Efficient formulation and computer implementation of the active-space electron-attached and ionized equation-of-motion coupled-cluster methods. Journal of Chemical Physics 2006, 125.   (44) Marcus, R. A. ON THE THEORY OF OXIDATION-REDUCTION REACTIONS INVOLVING ELECTRON TRANSFER .1. Journal of Chemical Physics 1956, 24, 966-978.   (45) Kresse, G.; Furthmuller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science 1996, 6, 15-50.   (46) Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Physical Review B 1996, 54, 11169-11186.   (47) Blochl, P. E. PROJECTOR AUGMENTED-WAVE METHOD. Physical Review B 1994, 50, 17953-17979.   (48) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Physical Review B 1999, 59, 1758-1775.  193   (49) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207-8215.   (50) Heyd, J.; Scuseria, G. E. Efficient hybrid density functional calculations in solids: Assessment of the Heyd-Scuseria-Ernzerhof screened Coulomb hybrid functional. Journal of Chemical Physics 2004, 121, 1187-1192.   (51) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential (vol 118, pg 8207, 2003). J. Chem. Phys. 2006, 124.   (52) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 2006, 125.   (53) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787-1799.   (54) Alder, R. W.; Arrowsmith, R. J.; Casson, A.; Sessions, R. B.; Heilbronner, E.; Kovac, B.; Huber, H.; Taagepera, M. PROTON AFFINITIES AND IONIZATION ENERGIES OF BICYCLIC AMINES AND DIAMINES - THE EFFECTS OF RING STRAIN AND OF 3-ELECTRON SIGMA-BONDING. Journal of the American Chemical Society 1981, 103, 6137-6142.   (55) Alder, R. W. MEDIUM-RING BICYCLIC COMPOUNDS AND INTRABRIDGEHEAD CHEMISTRY. Accounts of Chemical Research 1983, 16, 321-327.   (56) Alder, R. W. INTRABRIDGEHEAD CHEMISTRY. Tetrahedron 1990, 46, 683-713.   (57) Braida, B.; Hiberty, P. C.; Savin, A. A systematic failing of current density functionals: Overestimation of two-center three-electron bonding energies. Journal of Physical Chemistry A 1998, 102, 7872-7877.   (58) Zhang, Y. K.; Yang, W. T. A challenge for density functionals: Self-interaction error increases for systems with a noninteger number of electrons. Journal of Chemical Physics 1998, 109, 2604-2608.   (59) Grafenstein, J.; Kraka, E.; Cremer, D. Effect of the self-interaction error for three-electron bonds: On the development of new exchange-correlation functionals. Physical Chemistry Chemical Physics 2004, 6, 1096-1112.   (60) Rosso, K. M.; Dupuis, M. Reorganization energy associated with small polaron mobility in iron oxide. J. Chem. Phys. 2004, 120, 7050-7054.   (61) Hamann, T. W. Splitting water with rust: hematite photoelectrochemistry. Dalton Trans. 2012, 41, 7830-7834. 194   (62) Zhang, J.; Chen, X.; Takanabe, K.; Maeda, K.; Domen, K.; Epping, J. D.; Fu, X.; Antonietti, M.; Wang, X. Synthesis of a Carbon Nitride Structure for Visible-Light Catalysis by Copolymerization. Angewandte Chemie-International Edition 2010, 49, 441-444.   (63) Zhang, J. S.; Sun, J. H.; Maeda, K.; Domen, K.; Liu, P.; Antonietti, M.; Fu, X. Z.; Wang, X. C. Sulfur-mediated synthesis of carbon nitride: Band-gap engineering and improved functions for photocatalysis. Energy Environ. Sci. 2011, 4, 675-678.     195  CHAPTER 7: BIPOLARONIC RELAXATION IN GRAPHITIC CARBON NITRIDES  7.1 Introduction In Chapter 6 we discussed the impact of N1-N1 contacts on the electronic structure of graphitic carbon nitrides (g-CN).  We showed that N1-N1 contacts influence the nature of the band gap, with a direct band gap in the polymeric phase when these contacts are weak (longer interatomic distance), and an indirect band gap in the 2-D phase when they are more pronounced (shorter interatomic distance).  This suggests that N1-N1 contacts may play an important role in the photocatalysis of g-CN, since the nature of the band gap affects both the absorption and rate of radiative recombination in a material. We showed that Pauli repulsion at N1-N1 contacts is alleviated through the removal of an electron, allowing relaxation in the N1-N1 distance and formation of a three-electron bond.  This polaronic relaxation is expected to impact charge transport by trapping free holes (see Figure 7.1 for illustration).  The purpose of the present work is to investigate the possibility of a similar mechanism for bipolaronic relaxation with the formation of a two-electron bond.   Figure 7.1 Hypothesized mechanism for bipolaron formation in g-CN The hypothesized mechanism for three-electron (polaronic) and two-electron (bipolaronic) bond formation in g-CN.  Note the aromaticity of the newly-formed five-membered ring in the bipolaron.    196  Bipolaronic relaxation, whereby two of the same charge carriers localize at a single nuclear lattice distortion, will only occur if the lattice distortion releases enough energy to overcome the Coulombic repulsion between them.  The structures of some materials enable bipolaronic relaxation through a chemically intuitive mechanism.  One such material is doped polythiophene1, in which the removal of two electrons is followed by a change from aromatic to quin-bonding in the thiophene bridge2.  A similar bipolaronic relaxation mechanism has been identified in general polyparaphenylenes (PPP)3.  In both cases bipolaronic relaxation results in diminished charge transport.  Identification of a similar mechanism in g-CN would improve our understanding of what structural properties affect charge transport, and how the mechanism of charge transport impacts photocatalytic efficiency. The mechanism of charge transport has important implications for photocatalysis.  Previous work has suggested that, when g-CN is used to catalyze the water-splitting reaction, nitrogen atoms serve as catalytic sites for water oxidation while carbon atoms serve as catalytic sites for hydrogen reduction4.  Polaronic and bipolaronic relaxation may impact catalysis by trapping holes at lattice distortions on the surface, where they are accessible to substrate molecules, or by trapping holes in the material bulk, where they are inaccessible to substrate.  Computational modeling of g-CN under different dielectric conditions would provide insight into the preferred mode of nuclear relaxation (polaronic or bipolaronic) at the surface and in the material bulk.  We begin the present work with an electronic structure investigation of dicationic molecular models of g-CN, and then compare the structural properties of dicationic and cationic models (Chapter 6) in order to distinguish the properties of a single polaron and a bipolaron. We demonstrate that two-electron bond formation and aromatic stabilization increase the favorability of bipolaron formation from two polarons.  Marcus theory reorganization energies are then used 197  to discuss the relative lifetimes of polarons and bipolarons.  Finally, we discuss consistencies between conclusions in the present work and recent experimental studies of g-CN.  7.2 Methods 7.2.1 Electronic structure calculations  The molecules shown in Figure 6.2 were used to model g-CN, in correspondence with the work in Chapter 6.  All structures were optimized using density functional theory (DFT) in TeraChem5-8, a graphical processing unit-accelerated electronic structure package, and the cc-pVDZ basis set9.  The CAM-B3LYP10 functional was used in order to improve the description of long--systems of the models.  All singlet electronic state energies were calculated using a restricted wave function.  Doublet and triplet electronic state energies were calculated using an unrestricted wave function.  In the present work the term t spin-state is the the singlet spin-state is lowest in energy.  Potential energy surfaces were calculated as a function of the N1-N1 distance  for select models using constrained optimizations of this internal coordinate. In order to estimate the error in DFT energies of larger models we have performed coupled cluster (CC) calculations with single and double excitations (CCSD), as well as the triples-corrected left-eigenstate completely renormalized variant of CCSD, CR-CC(2,3). These calculations were completed for melam (Figure 1B) using the GAMESS11 implementation12,13 and CAM-B3LYP optimized structures.  The CR-CC(2,3) method was chosen for its accurate treatment of bond breaking14, a consequence of a non-iterative triples correction to the energy.  198  The CC second ionization potential was calculated using the one-hole, two-particle variant of the electron-attachment equation of motion coupled cluster theory with single and double excitations15-19, EA-EOM-CCSD(1h,2p)/cc-pVDZ (referred to in the present study as EA-EOM), as implemented20,21 in GAMESS.  This calculation is performed using the dicationic closed shell wave function as a reference, in order to determine the ionization potential of the cation.  Solvent effects were modeled implicitly through the polarizable continuum model22,23 (PCM).  Unless otherwise noted, results reported in this work were obtained using PCM.  7.2.2 Decomposition of the dicationic relaxation energy  The hypothesized mechanism for bipolaronic N1-N1 bond formation is shown in Figure 7.1.  This mechanism illustrates that relaxation in the N1-N1 distance can cause a free hole to become trapped, forming a new three-electron bond.  If this process occurs a second time, the two-electron bond that is formed may result in a new heteroaromatic ring.  The aromatic stabilization energy (ASE) of this ring is expected to play an important role in counteracting unfavorable Coulombic repulsion between polarons on adjacent monomeric units.  In an effort to quantify the energy released by bond formation and aromatic stabilization we introduce a quantity we call the nuclear relaxation energy, .  Here Z and  denote different charge states for the molecule.  The relaxation energy is defined as  .  (7.1) where RZ and R denote the minimum energy nuclear coordinates when the charge is Z and Z, respectively.  In some cases, the name of the molecule is also listed with the charge: . For example, the change in energy when dicationic melam (Figure 1B) relaxes from the cationic to the dicationic minimum is defined as 199   . (7.2) that the molecular charge has changed in addition to the nuclear geometry: .  The first adiabatic ionization potential, for example, is defined as  . (7.3) The dicationic relaxation energy encompasses both the bond formation and aromatic stabilization energies.  An approach for approximation of one of these quantities will thus allow estimation of the other.  The 1,3,5-triazine dimer (Figure 7.2) allows calculation of the N1-N1 dicationic two-electron bond energy without energetic contributions from aromatic stabilization.  This was achieved through an initial optimization of individual cationic triazine monomers, followed by optimizations of the cationic and dicationic dimers.  The triazine dimer dihedral angle, (Figure 7.2A), was constrained to values of either 0º or 90º, while the N1-N1 distances  at the cationic and dicationic minima were constrained to their values in melam.  Dihedral angles of 0º and 90º were used in order to assess the impamonomers.  A corresponding dihedral angle, , was defined for select atoms that occur -monomeric dihedral -CN.   Once the most appropriate triazine dimer dihedral angle was identified, the dicationic melam relaxation energy () was used to approximate the energy released by aromatic stabilization as  . (7.4) 200  Eqn. (7.4) These quantities are reported for CR-CC(2,3) and CAM-B3LYP, and were calculated using the basis sets and software packages described in section 7.2.1.  Figure 7.2 Models and atomic labels for the analysis of aromatic stabilization The 1,3,5-triazine dimer is shown (A), where select carbon and nitrogen atoms are labeled with subscripts used to define the dihedral angle, .  Also shown (B) is a schematic close-up of the atoms used to calculate the inter-monomeric dihedral angle,   7.2.3 Nuclear relaxation and charge transport  Due to the importance of N1-N1 Pauli repulsion and the hypothesized bond formation mechanism, properties of the model systems in Figure 1 were calculated as a function of the N1-N1 distance through constrained optimizations.  The N1-N1 bond order was determined from the off-diagonal elements of the density matrix at each of these structures.  Vibrational frequencies for select molecular models are available in Appendix A. The error in the DFT nuclear relaxation energies for all systems, , was estimated by taking the difference between the DFT and CC energies in melam: 201   . (7.5) A CC-corrected DFT nuclear relaxation energy was then calculated for molecule X as:  . (7.6) Note that different CC methods were used, depending upon the charge of the molecule.  When Z=0 the CCSD energy was used, when Z=+1 the EA-EOM energy was used, and when Z=+2 the CR-CC(2,3) energy was used.  The CR-CC(2,3) method includes a triples correction that is absent in the EA-EOM(1h,2p) method.  This results in a more accurate treatment of electron correlation in the dicationic molecular models, in which bond formation is hypothesized.  Marcus theory was applied in order to analyze the relative hopping rates of polarons and bipolarons in g-CN, where the rate constant for hopping may be expressed as:  . (7.7) Here  is a pre-exponential factor, and  , (7.8) where  is the reorganization energy.  Treating polaron hopping as a self-exchange reaction, in which , the rate constant may be expressed as  . (7.9) Forming a bipolaron from two polarons is not a self-exchange reaction.  We may express this reaction, for molecule X, as:  . (7.10) The adiabatic ionization potentials allow quantification of  for eqn. (7.10) as  . (7.11) 202  Eqn. (7.11) enables analysis of the stability of polarons and bipolarons in g-CN from a thermodynamic perspective, while insertion of eqn. (7.11) into the Marcus theory expression for the rate constant allows determination of relative hopping rates for polarons and bipolarons. Calculation of  is beyond the scope of the present work. Reorganization energies for a reaction are defined via the nuclear relaxation energies.  The reorganization energy for the self-exchange of the neutral and polaronic states is defined as  , (7.12) and the reorganization energy for the self-exchange of the neutral and bipolaronic states is defined as  . (7.13) Eqn. (7.13) however, it is also possible for the minimum energy dicationic state to be a polaron pair (triplet electronic state).  This reorganization energy is instead expressed as , which has the same definition as   (eqn. (7.13)), with the important difference that the energies of triplet electronic spin states were used during the calculation.   First and second adiabatic ionization potentials (IPs) were calculated using a range of dielectric constants.  Trends in the IPs as a function of the dielectric may provide some insight into the preferred relaxation mechanism (polaronic or bipolaronic) in g-CN at the surface and in the bulk.  We have found no dielectric constant for g-CN reported in the literature.  The dielectric range reported for isolated carbon nitride thin films, which are distinct from but structurally related to g-CN, is ~2-4.24  Consequently we assume a dielectric of 3.0 for g-CN in the present work.  203  Pure water has a dielectric of 80.1.  The dielectric constant for water was used to describe the -CN in solution.  7.3 Two-electron bond formation and bipolaronic relaxation The properties of bipolaronic g-CN are best understood in the context of our previous work on cationic molecular models (Chapter 6).  This work demonstrated the importance of N1-N1 Pauli repulsion.  Analysis of the hole orbitals (orbital from which an electron is removed) for the systems in Figure 6.2 showed that at the neutral minimum, ionization occurs from a nitrogen-molecular orbital.  As the distance between position 1 nitrogen atoms decreases towards the and become nitrogen lone pair (n-type) orbitals.  More specifically, these n-type molecular orbitals are composed from an out-of-phase superposition of N1 lone pair orbitals from neighboring monomers.  This orbital character is reflective of significant Pauli repulsion between N1 lone pairs.   Removal of a second electron occurs from the same lone pair orbital at the cationic minimum.  Figure 7.3 shows the dicationic lowest unoccupied Kohn-Sham molecular orbitals (LUMO) at the cationic and dicationic minimum structures, where the dicationic LUMO may be considered a hole orbital for the dication.  All results reported in the current section were obtained from gas phase calculations (without application of PCM).  The dicationic LUMO orbitals at the cationic minima exhibit the same out-of-phase n-type orbitals seen in cationic models (see the cationic minima orbitals in Figure 7.3 for an example).  At the dicaorbital for all models except the 2-D melem hexamer, with character that reflects aromatic stabilization across the newly formed ring.  This change in orbital character is a result of bond formation. 204   Figure 7.3 Dicationic Kohn-Sham molecular orbitals for all models of g-CN Dicationic lowest unoccupied Kohn-Sham molecular orbitals are shown for melam (A), melem dimer (B), polymeric melem trimer (C), 2-D melem trimer (D), and 2-D melem hexamer (E) at their cationic and dicationic minimum structures.  Nitrogen atoms are blue, carbon atoms are gray and hydrogen atoms are white.  All models exhibit a large change in N1-N1 distance when relaxing from the cationic to the dicationic minimum structure.  These distances are reported in Table 7.1 (for N1-N1 pairs that experience significant nuclear relaxation).  The cationic N1-N1 distances in each model are nearly identical, with values of 2.1-2.2 Å; the dicationic N1-N1 distances are similar in all cases except the 2-D melem hexamer, with lengths of 1.39-1.42 Å.   This is similar to the equilibrium length of a two-electron bond between two nitrogen atoms (1.45 Å), and supports the hypothesis that bond formation results from dicationic nuclear relaxation.     205  Table 7.1 Structural and bond order data for models of g-CN Neutral, cationic and dicationic minimum N1-N1 distances  and bond orders are reported for each molecular model.  Note that two distances and bond orders are reported in the dicationic 2-D melem hexamer, which exhibits a polaron pair.   neutral (Å) cationic (Å) dicationic (Å) cationic bond order dicationic bond order melam 2.89 2.16 1.40 0.32 0.96 melem dimer 2.85 2.13 1.39 0.34 1.00 polymeric melem trimer 2.85 2.13 1.40 0.34 0.99 2-D melem trimer 2.64 2.17 1.42 0.33 0.99 2-D melem hexamer 2.69 2.20 2.21/2.21 0.31 0.32/0.33  molecular orbital localized at a single N1-N1 pair, the 2-D melem hexamer has an n-type dicationic LUMO localized at two different N1-N1 pairs on opposite sides of the molecule.  This is a reflection of the fact that the minimum energy spin-state for the dicationic 2-D melem hexamer is a triplet, suggesting that a polaron pair is favored over a bipolaron in the 2-D phase.  This will be discussed further in Section 7.6.  Comparison of the 2-D melem hexamer LUMO at the cationic and dicationic minima shows that the first polaron migrates from the middle of the structure to an edge N1-N1 pair in order to accommodate a second polaron.  This is likely caused by steric strain.  Strain is diminished in the polymeric phase, however, which may provide a partial explanation for the thermodynamic favorability of polaron pair formation in the 2-D melem hexamer. The polymeric models (and even the 2-D melem trimer) all favor bipolaron formation. Having discussed the unique properties of the 2-D melem hexamer, we will now return our focus to a more detailed analysis of the smaller molecular models.  Potential energy surfaces (PESs) for dicationic melam, melem dimer, and 2-D melem trimer (Figure 7.4) demonstrate a rapid 206  decrease in energy with decreasing N1-N1 distance.  Dicationic melam (Figure 7.4A) relaxes in energy by 3.01 eV when the N1-N1 distance decreases from 2.16 to 1.40 Å, while the melem dimer and 2-D melem trimer relax by 3.09 and 1.92 eV, respectively, with similar changes in N1-N1 distance.  Figure 7.4B shows that this relaxation in energy is accompanied by an increase in the N1-N1 bond order.  At the dicationic minimum geometry the N1-N1 bond order of all models except the 2-D melem hexamer is approximately 1.0 (Table 7.1); this is the bond order for a two-electron bond.  Thus, both the bond lengths and bond orders observed in these systems suggest that a new two-electron bond is formed at the dicationic minimum.  The 2-D melem hexamer is the exception, possessing two N1-N1 pairs with bond orders of 0.324 and 0.325. As described in section 7.2.3, we assessed the error in DFT through the application of CC theories to melam.  Figure 7.5 shows PESs for dicationic melam as a function of N1-N1 distance using CAM-B3LYP, CCSD, and CR-CC(2,3).  At the cationic minimum distance  there is notable disagreement between the relative CR-CC(2,3) and CAM-B3LYP/CCSD dicationic energies.  This disagreement is a result of N1-N1 bond formation.  The dicationic energy at the cationic minimum structure is the energy of the broken bond; CAM-B3LYP and CCSD are not well-suited to describe bond-breaking, while the triples-correction in CR-CC(2,3) enables an accurate treatment of this process.  Consequently, the vertical ionization potential of the cation (at R+1) is 0.41 eV greater with CAM-B3LYP than with CR-CC(2,3).  The differences between the DFT and CC nuclear relaxation energies that contribute to the neutral/polaronic self-exchange reaction were calculated for melam.  It was found that  and . This is a small difference, and suggests that DFT is a suitable method to calculate thermodynamic quantities for the self-exchange of a neutral and polaronic molecule. 207   Figure 7.4 Dicationic melam and melem dimer PESs CAM-B3LYP dicationic potential energy surfaces (A) are shown for melam, melem dimer, and the 2-D melem trimer over a range of N1-N1 distances.  The N1-N1 bond order is shown (B) over the same range of distances.  Table 7.2 Nuclear relaxation energies for models of g-CN with solvent dielectric Nuclear relaxation energies at the surface  are reported in eV.       melam 1.51 0.66 4.76 7.51 melem dimer 1.58 0.23 3.94 7.22 polymeric melem trimer 1.63 0.13 3.66 6.21 2-D melem trimer 0.99 0.39 2.21 6.98 2-D melem hexamer 1.04 0.13 0.59 1.57 208    Figure 7.5 Dicationic melam DFT and coupled cluster theory PESs Dicationic potential energy surfaces are shown for melam using CAM-B3LYP, CCSD, and CR-CC(2,3) over a range of N1-N1 distances.  7.4 N1-N1 bond formation and aromatic stabilization The enthalpy for a single bond between two nitrogen atoms is 1.67 eV25.  By comparison, the melam dicationic relaxation energy is 2.60 eV (Table 7.2).  As discussed previously, aromatic stabilization is expected to contribute to the dicationic relaxation energy, since this relaxation centered on N1 atoms from neighboring monomers.  The dicationic minimum LUMOs in Figure 7.3 (for all molecules except the 2-D melem hexamer) de Analysis of the dicationic relaxation energies in the triazine dimer and melam, as described in subsection 7.2.2, allows estimation of the aromatic stabilization energy ().  The triazine dimer N1-N1 distance was constrained to the values at the cationic and dicationic minima in melam (2.16 and 1.40 Å, respectively).  Table 7.3 lists  (eqn. (7.4)) for CAM-B3LYP and CR-209  CC(2,3), with inter-monomeric dihedral angles  of 0º and 90º.  These -system interactions, which are maximized when  and minimized when .  When  the CAM-B3LYP dicationic relaxation energy is 0.4 eV larger than the same CR-CC(2,3) energy in the triazine dimer.  In melam the CAM-B3LYP dicationic relaxation energy is 0.4 eV smaller than the CR-CC(2,3) energy.  In total, the CAM-B3LYP  is 0.8 eV lower than in CR-CC(2,3).  The same difference is observed when the dihedral angle is 90º.  These calculations suggest that  may range between 1.10 and 1.54 eV for bipolaronic relaxation at N1-N1 pairs.    Table 7.3 Quantities related to the aromatic stabilization energy Dicationic relaxation energies (gas phase) are reported for melam ( ) and the triazine dimer  with different values for the dihedral angle , as depicted in Figure 7.2A.  These quantities are used to calculate the aromatic stabilization energy, . All quantities are reported for the CAM-B3LYP and CR-CC(2,3) methods.    (triazine dimer)     (eV)   (eV)   (eV) CAM-B3LYP 0º 1.86 2.60 0.74 CAM-B3LYP 90º 2.28 2.60 0.32 CR-CC(2,3) 0º 1.47 3.01 1.54 CR-CC(2,3) 90º 1.91 3.01 1.10     210  Table 7.4 Inter-monomeric dihedral angles for models of g-CN The inter-monomeric dihedral angle, , as depicted in Figure 7.2B.  Note that in the 2-D melem hexamer two inter-monomeric dihedral angles are reported since a polaron pair is observed.     melam 2.3º melem dimer 4.0º polymeric melem trimer 4.2º 2-D melem trimer 0.0º 2-D melem hexamer 10.6º/11.5 º  Our approach for estimating  makes some significant approximations and merits further discussion.  First and foremost, we assume that the aromatic stabilization and bond formation energies are the only significant contributors to .  The approximate nature of this assumption is apparent when we consider that some nuclear relaxation will counteract aromatic stabilization or N1-N1 bond formation.  Out-of-plane nuclear relaxation is observed in the structure of the dicationic melem dimer, for example, as shown in Figure 7.3B.  At R+2 the interacting terminal amino groups have relaxed out-of-plane in order to minimize steric hindrance.  This out-of-plane relaxation will likely offset some aromatic stabilization and reduce the dicationic relaxation energy.  We may assess the extent of out-of-plane relaxation using the inter-monomeric dihedral angle, , as reported in Table 7.4.  All dihedral angles are small, justifying the use of a triazine dimer dihedral angle of 0º to estimate the energy released by aromatic stabilization.  When the dihedral angle is 0º, CR-CC(2,3) calculations suggest that the aromatic stabilization energy may be as large as the bond formation energy (1.5 eV).  Thus, despite the approximate nature of our approach,  provides valuable insight into the large dicationic relaxation energy that is observed in these model systems. 211  7.5 Polaronic and bipolaronic hopping rates  Relative polaronic and bipolaronic hopping rates were estimated using the neutral/polaronic , neutral/bipolaronic , and neutral/polaron pair  self-exchange reorganization energies, as defined in eqns. (7.12) and (7.13).  Table 7.5 lists these quantities for the models in the present work.  The neutral/bipolaronic reorganization energy is at least 7 eV larger than the neutral/polaronic reorganization energy.  This large difference suggests that, after N1-N1 bond formation, bipolarons are kinetically trapped for a much longer time than polarons.  The consistent difference between the neutral/polaronic and neutral/bipolaronic reorganization energies suggests that this result does not depend on the phase of the material.  Overall, this result suggests that bipolaron formation is favorable at the surface in g-CN.  Higher accuracy electronic structure theory, as well as explicit solvent modeling, would provide a more quantitative comparison of polaronic and bipolaronic reorganization energies at the surface of g-CN.  Table 7.5 Polaronic self-exchange reaction reorganization energies Reorganization energies for the self-exchange reactions of neutral/polaronic models  and neutral/(bipolaronic/polaron pair) models  are reported for surface g-CN .    (surface)  (surface) melam 2.17 12.27  melem dimer 1.81 11.16 polymeric melem trimer 1.77 9.87 2-D melem trimer 1.38 9.19  212  The reorganization energies enable analysis of polaronic and bipolaronic hopping rates in the framework of Marcus theory.  Determination of rate constants is beyond the scope of the present study due to the complexity of calculating pre-exponential factors, which are expected to differ significantly for polaronic and bipolaronic/polaron pair hopping.  Comparison of the neutral/polaronic and neutral/bipolaronic reorganization energies (Table 7.5), however, provides valuable insight into relative hopping rates.  This comparison is justified by the fact that the reorganization energies are different by more than 7 eV in all cases.  This suggests that the bipolaronic hopping rate should be orders of magnitude slower than the polaronic hopping rate in g-CN.  Both hopping rates, however, should be exceedingly slow, indicating that polarons and bipolarons may have a long lifetime in g-CN. Discussion of relative hopping rates in the heterogeneous reaction that converts two polaronic molecules to a bipolaron and a neutral molecule requires . Table 7.6 lists the values for the adiabatic ionization potentials which contribute to , as well as .  The negative sign of   in all molecular models suggests that the formation of a bipolaron/polaron pair and neutral molecule from two polarons is thermodynamically favorable at the surface of g-CN (both polymeric and 2-D phase).    213  Table 7.6 Free energy of self-exchange for a polaron in solvated g-CN , and quantities relevant to the calculation of  for the reaction defined in eqn. (7.10).   is calculated as defined in eqn. (7.11)).      (surface)  (surface)  (surface) melam 6.20 10.89 -1.51 melem dimer 6.46 11.29 -1.63 polymeric melem trimer 6.58 11.52 -1.64 2-D melem trimer 6.32 12.46 -0.84 2-D melem hexamer 6.65 13.05 -0.25   Figure 7.6 shows the first and second adiabatic IPs as a function of dielectric constant for all systems except the polymeric melem trimer; graphical results are omitted for the polymeric trimer since they are visually indistinguishable from the melem dimer.  In melam, the second IP is lower than the first IP above a dielectric of 2.5.  In the melem dimer the second IP is lower than the first IP above a dielectric of 2.0.  In the 2-D melem trimer and hexamer these crossover dielectric values are 14.0 and 7.0, respectively.  These results further suggest that bipolaron formation is favored over polaron formation for a broad range of dielectrics, and is competitive with polaron formation at low dielectrics which include the estimated dielectric for g-CN. 214   Figure 7.6 Adiabatic ionization potentials and dielectric for models of g-CN The first and second adiabatic ionization potentials (IPs) are shown for melam (A), melem dimer (B), 2-D melem trimer (C), and the 2-D melem hexamer (D).  Inset with each figure is a closer look at the range of dielectric values for which the first and second IPs are similar.  Note the difference in the scales of the y-axes for the polymeric (melam and melem dimer) and 2-D (trimer and hexamer) ionization potentials.  7.6 Discussion of g-CN charge transport incorporating recent experimental work We now address how our interpretation of these results is influenced by recent experimental work on g-CN.  To summarize, thus far we have demonstrated that: 1) dicationic models of g-CN exhibit two-electron bond formation between nitrogen atoms on adjacent monomeric units 215  2) Aromatic stabilization effectively doubles the energy released as this bond is formed (3.0 eV with ASE, ~1.5 eV without ASE) 3) This molecular mechanism is expected to mirror a bipolaronic relaxation mechanism in g-CN 4) Comparison of the reorganization energies for the self-exchange reactions, in which a cationic/dicationic molecule is formed from a neutral molecule, suggest that both polarons and bipolarons have a long lifetime in g-CN, though the bipolaronic lifetime is expected to be orders of magnitude longer than the polaronic lifetime 5) Comparison of free energy changes for the formation of a dicationic and neutral molecule from two cationic molecules suggests that bipolarons are favored in both the polymeric and 2-D phases at the surface of g-CN 6) First and second adiabatic ionization potentials further suggest that bipolaron formation is favored over a single polaron in the polymeric phase, and less favorable than a single polaron in the 2-D phase Few experimental studies have directly addressed the catalytic mechanism of g-CN, however Wang, et al.4, have proposed that nitrogen atoms act as catalytic sites for the water oxidation reaction and carbon atoms for the hydrogen reduction reaction.  With respect to charge transport in g-CN, Merschjann, et al.26, recently showed that g-CN exhibits low intra-planar charge carrier mobilities. We will discuss the ways in which these results pertain to the conclusions above.  The mechanism of bipolaronic relaxation involves localization of two holes between two nitrogen atoms.  Wang, et al., propose that nitrogen atoms act as catalytic sites for water oxidation.  If this proposal is correct, our results suggest that bipolaronic relaxation may have a direct impact on photocatalysis through the localization of two positive charges at a single catalytic site.  The 216  localization of multiple holes is expected to be catalytically beneficial for water splitting, in particular, which is a four electron process.  Further work with explicit solvent modeling would provide insight into the catalytic role of bipolaronic relaxation in g-CN.   Merschjann, et al., observe low intra-planar charge carrier mobilities in g-CN.  Large neutral/polaronic and neutral/bipolaronic reorganization energies are consistent with this observation.  Calculation of a rate constant for polaron/bipolaron hopping is beyond the scope of the present work, but would allow quantitative comparison of theoretically and experimentally obtained hopping rates.  This comparison would aid investigation of the dominant relaxation mechanism (polaronic/bipolaronic) in g-CN. We conclude by discussing the ramifications of these results for our understanding of the photocatalytic mechanism in g-CN.  The change in the free energy for the reaction that converts two polaronic molecules to a polaronic and neutral molecule is significant (>0.25 eV) at the surface, suggesting that bipolaron formation is favored.  Polarons/bipolarons on the material surface may serve a catalytic role, while those trapped in the bulk are inaccessible to solvent.  The polaronic/bipolaronic relaxation mechanism identified in this work thus may provide a partial explanation for low photocatalytic efficiencies in g-CN: polarons/bipolarons trapped in the bulk are inaccessible to solvent and unable to catalyze the water oxidation reaction.  Further work with layered models and, as previously mentioned, explicit solvent would allow exploration of this hypothesis.    217  7.7 Conclusions    The present work complements existing work on g-CN by providing two new and important insights into the charge transport and photocatalytic mechanism.  First, calculations on dicationic molecular models suggest that a two-electron bond between nitrogen atoms may be formed with bipolaronic relaxation.  This two-electron bond releases 3.0 eV, with an estimated 1.5 eV released from aromatic stabilization.  Second, large reorganization energies suggest that bipolarons are favored over polarons at the surface of g-CN.  Furthermore, large reorganization energies suggest low charge carrier mobilities, in agreement with experimental observations. The low mobilities of polarons and bipolarons indicate that they are unlikely to reach the surface of the material quickly if they are generated in the bulk.    Further work is needed to determine whether polaronic/bipolaronic relaxation in the bulk is directly related to low photocatalytic efficiencies in g-CN.    218             APPENDICES    219  APPENDIX A: Supplementary content   Table A7.1 Dicationic melam vibrational frequencies (R-CAM-B3LYP/cc-pVDZ)  Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 9.19 0.00737 710.90 1.39309 1173.17 0.96667 8.02 0.00515 716.82 0.21927 1177.61 0.04654 7.62 0.00736 738.48 0.09995 1191.63 0.20362 0.02 0.00002 751.66 0.03971 1317.25 0.13037 0.01 0.00424 793.91 0.00050 1347.64 0.31516 0.03 0.00000 808.55 2.57160 1377.4 0.0405 47.55 0.06836 716.82 0.21927 1397.31 3.17671 71.52 0.00125 738.48 0.09995 1425.82 0.06472 137.50 0.95472 751.66 0.03971 1499.45 0.00365 160.19 0.01081 793.91 0.00050 1514.88 6.89531 191.73 0.00895 808.55 2.57160 1527.09 9.85712 197.54 1.45121 820.92 0.23388 1548.87 0.09134 203.83 1.00260 980.06 0.03033 1553.97 1.34925 206.14 0.07591 984.97 0.21095 1576.2 0.35776 281.53 0.09757 997.29 0.63652 1606.36 0.21586 310.60 0.35001 1046.04 0.00471 1695.77 3.61166 342.19 1.60326 1057.26 0.52783 1712.98 25.452 354.76 0.00303 1134.61 0.57833 1734.35 64.23654 370.14 0.31729 1173.17 0.96667 1758.12 0.24213 413.96 0.15428 1177.61 0.04654 1813.21 34.75014 521.57 0.00453 1191.63 0.20362 1838.56 3.17029 554.12 1.40582 1317.25 0.13037 3542.8 9.64145 560.21 0.21539 1347.64 0.31516 3544.03 9.08755 577.71 0.07448 1377.4 0.0405 3548.89 14.00762 593.02 0.65288 1397.31 3.17671 3550.96 0.98207 598.18 0.80291 1425.82 0.06472 3575.23 6.93299 611.91 1.86313 1499.45 0.00365 3675.53 0.10524 619.42 0.49745 1514.88 6.89531 3677.18 5.95486 626.64 0.73718 820.92 0.23388 3680.4 4.2232 629.10 11.08892 980.06 0.03033 3680.68 4.79556 642.89 0.13133 984.97 0.21095 670.50 0.05578 997.29 0.63652 670.78 0.19905 1046.04 0.00471 684.08 0.00916 1057.26 0.52783 691.75 0.12089 1134.61 0.57833 220   Figure A7.1 Melam dicationic vibrational frequencies Frequencies and IR intensities are plotted for melam at the dicationic minimum geometry.    221  Table A7.2 Dicationic melem dimer vibrational frequencies (R-CAM-B3LYP/cc-pVDZ)  Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu Frequency (cm-1) IR Intensity (Debye/Å)2/amu 8.05 0.00081 598.35 2.45510 1222.13 4.79583 5.53 0.00027 598.48 0.10004 1268.39 1.84035 5.22 0.00281 620.74 0.30039 1278.07 0.15703 0.13 0.00000 648.06 0.50435 1345.97 0.65481 0.02 0.00145 650.3 0.29464 1350.37 1.43649 0.02 0.00092 659.57 0.02241 1378.34 5.00702 28.85 0.00993 669.59 0.05917 1382.59 1.52929 35.00 0.00091 671.02 0.08949 1391.55 1.01731 68.71 0.00197 673.51 0.01364 1405.88 0.70395 94.91 0.00928 680.88 0.0143 1436.09 2.82459 100.11 0.02388 690.24 0.07277 1456.64 1.34146 103.30 0.50008 711.07 0.01136 1460.77 2.35395 104.97 0.04382 711.72 0.29267 1500.28 13.97976 162.29 0.08564 731.33 0.00551 1510.94 8.80751 176.91 0.00099 738.5 0.00118 1516.24 1.64081 189.43 0.00478 743.33 0.31003 1528.69 3.79874 225.86 0.39585 750.5 0.15132 1529.42 3.58612 229.75 0.08731 759.65 0.05632 1583.9 13.84566 241.23 0.03085 766.26 0.05284 1587.33 0.3207 249.40 0.00741 782.36 0.00094 1594.49 7.48006 249.70 0.00728 795.85 0.01455 1602.7 16.30599 265.31 0.00070 795.97 0.13015 1707.21 100.2219 286.80 0.25941 806.79 0.08133 1711.89 25.15662 297.23 0.14585 815.41 0.29693 1737.95 61.7428 306.69 0.08965 827.54 0.1127 1744.96 8.3276 324.92 0.09726 828.16 3.19966 1745.82 0.7937 361.69 0.05025 832.25 0.60241 1788.57 12.27852 407.57 0.00212 923.46 0.21327 1789.12 19.79158 411.69 0.02921 983.62 0.10155 1803.03 30.75576 455.09 0.03089 1004.9 0.08467 1813.79 11.19506 461.77 0.02037 1025.54 0.07023 1839.57 0.86421 471.37 0.00230 1040.53 0.97696 3563.9 20.60676 507.88 4.83627 1070.78 0.83775 3564.61 3.88755 521.86 0.25444 1077.11 0.45599 3576.15 4.69818 525.84 1.54517 1078.68 0.43609 3576.6 10.5406 545.53 1.16721 1129.39 0.10766 3584.19 6.41002 568.78 0.00333 1150.5 0.03269 3702.57 3.25605 577.80 3.92596 1160.96 0.08828 3703 3.69177 580.56 4.27805 1166.09 2.10043 3718.21 1.72324 586.98 0.14692 1188.2 0.00767 3718.68 4.5802 222   Figure A7.2 Melem dimer dicationic vibrational frequencies Frequencies and IR intensities are plotted for melem dimer at the dicationic minimum geometry.    223  APPENDIX B: Molecular geometries for dicationic models of g-CN  Melam dicationic minimum 26       N   -3.7751132251    0.2379852250    1.6639163610      H   -4.6472172337    0.4489169587    1.1783227879      C   -2.5897564024   -0.0338045097    1.0451776165      N   -2.3252078840   -0.0017747268   -0.2032202740      C   -0.9903189306   -0.1609689073   -0.4929006329      N   -0.0064811910   -0.0649454301    0.4297371558      C   -0.3180402741   -0.0364609342    1.6939624270      N    0.5983253883    0.1438002932    2.6406461759      H    0.3635190517    0.4162699080    3.5890494182      H    1.5449531274    0.3386995374    2.3167907923      N   -0.6633858468   -0.3453135193   -1.7468478769      H    0.3148432006   -0.4287904078   -2.0118840802      H   -1.3846589914   -0.3730190117   -2.4631191165      N   -1.6543067447   -0.2948699992    2.0502725496      N   -2.2947534055   -0.1029745973    3.2832757795      C   -3.6358481702    0.1891205660    3.0203189490      N   -4.5407316582    0.3135404448    3.9120578460      C   -4.1042511332    0.0455139589    5.1882537636      N   -2.9344342259   -0.5738774865    5.4655386686      C   -2.0602584667   -0.7407326894    4.5150531705      N   -0.9310641323   -1.4166176966    4.7079106970      H   -0.8277698933   -1.8639226205    5.6180304133      H   -0.3698095068   -1.7780500141    3.9441556059      N   -4.9077243291    0.3472291883    6.1765816641      H   -4.6421183628    0.1306881623    7.1340468501      H   -5.8135106058    0.7688197345    5.9871589320   40       N    -0.0386338116    1.0005354282    9.9065717603      C    -1.0246075718    0.1700483897    9.6059396449      N    -1.2066658468   -0.2074529026    8.2458389888      C    -0.2592235825    0.2218321455    7.3418405013      N     0.6863948844    1.0349383786    7.6032002218      C     0.7363142825    1.4515256609    8.9295366717      C    -2.3658792530   -0.9353987364    7.8903281987 224       N    -2.6556802911   -1.0599982043    6.5545906415      C    -1.7308037640   -0.7431850474    5.7410389872      N    -0.4199898073   -0.3827940068    6.0910848843      N    -1.8379745525   -0.3140781443   10.5109000331      C    -2.8325373911   -1.1137187921   10.1011482101      N    -3.1559912308   -1.3884797035    8.7894036194      N    -1.8169334653   -0.6281577025    4.3823385930      C    -0.6433020106   -0.1711983028    3.8524735845      N    -0.3506168205   -0.0253136738    2.6234074928      C     0.9558760744    0.2825744173    2.3372596135      N     1.9088859337    0.1254626175    3.3704312242      C     1.5846773421   -0.1330657092    4.6846765131      N     0.2218156775    0.0576722232    4.9335623708      C     3.2923817030    0.1762218478    3.0405460689      N     3.6413848180    0.5364086853    1.8308748708      C     2.6692026087    0.8022137868    0.9473589592      N     1.3184668094    0.6231246707    1.1581356797      N     2.4158645075   -0.4654735207    5.5914572507      C     3.7405034207   -0.5093200871    5.1687199504      N     4.1774524837   -0.1469681259    3.9702999037      N     3.0075966584    1.2403015192   -0.2457110629      N     4.6115036090   -0.9323028698    6.0665965757      N    -3.6055716900   -1.6699693320   11.0082870852      N     1.6551223198    2.3581006438    9.2083694613      H     2.2528628521    2.7279109646    8.4796067545      H     1.7251872967    2.7123828309   10.1575208257      H     5.5911578633   -1.0011827658    5.8080931358      H     4.3038080563   -1.2248198780    6.9857876627      H     2.2900434166    1.4234246362   -0.9396963811      H     3.9886948912    1.3746563138   -0.4708689848      H    -2.6478031782   -0.8524156705    3.8358935824      H    -4.3811006999   -2.2579561165   10.7204072060      H    -3.4324125422   -1.4909088681   11.9927897011   58       N    -4.4566192927    2.9044870198    1.1823368333      C    -3.3935431060    3.4137340159    0.5842111581      N    -2.4298873606    2.5420615497    0.0342096480      C    -2.6101717736    1.1571354819    0.1454164048      N    -3.6647535856    0.6676758440    0.7368503676      C    -4.5515997818    1.5730605054    1.2292122109      N    -1.6691704326    0.3393090975   -0.3722667123 225       C    -0.6625274769    0.9524575288   -0.9432733091      N    -0.4058556571    2.2313973699   -1.1341565373      C    -1.3130105296    3.0829356802   -0.6180217639      N    -3.2146775891    4.7176889825    0.4931582821      C    -2.1131365334    5.1422135192   -0.1377665601      N    -1.1519456631    4.3728741309   -0.7120708401      N     0.3750942304    0.1530237214   -1.5078629996      C     0.7417261018   -1.0979659635   -1.2744701070      N     1.8584773194   -1.4729593716   -2.0051615433      C     2.2486786721   -2.6839783856   -1.9031273846      N     1.5457709854   -3.6051066889   -1.0915912790      C     0.5135473612   -3.1063799419   -0.2448922133      N     0.1225418892   -1.8662168201   -0.3746145023      N     3.4140022341   -3.1144590596   -2.4941962336      C     3.6370201020   -4.3636134326   -2.4638117971      N     2.7547319905   -5.3482402286   -1.9863095622      C     1.8638905911   -4.9392784234   -0.9858393779      N    -0.0079390094   -3.9260279526    0.6584673829      C     0.4540393117   -5.1634775963    0.7408384192      N     1.3891833871   -5.7390236379   -0.1142782497      N     4.7884369996   -5.0126671577   -2.8189099012      N    -0.0041855070   -5.9568126962    1.6948157955      N    -1.9287836486    6.4573122116   -0.2155223035      N    -5.6288471843    1.0642230103    1.8256247736      H     0.9967365408    0.6896577573   -2.1119122453      H     5.6047349556   -4.5565990969   -3.2238627345      H    -1.1083341049    6.8185089645   -0.6833539926      H    -2.6137028465    7.0770437257    0.1974738453      H    -0.6789303920   -5.5889327065    2.3579215298      H    -5.7416337522    0.0615325149    1.8781768882      H     0.3408927452   -6.9036835144    1.7864268300      H    -6.3229206755    1.6931819760    2.2078288079      C     4.7084663818   -6.3443104186   -2.5369085258      N     3.4556494898   -6.5550689836   -1.9432930164      N     5.5685055418   -7.2458032889   -2.7987033514      C     5.1850999361   -8.5418421908   -2.5720216559      N     3.8118074232   -8.7928824835   -2.3449325099      C     2.8686430387   -7.8041160249   -2.1594931275      N     1.6080745218   -7.9920592886   -2.1397203097      C     3.3475379174  -10.1376284216   -2.3718574752      N     6.0245789797   -9.5064904905   -2.6404662473      N     2.0429413366  -10.3595150206   -2.3330276392      C     1.2213723704   -9.3196393728   -2.2774807219      N     4.2239774630  -11.1090722089   -2.4356019117      C     5.5216962281  -10.7820228550   -2.5102999393      N    -0.0814529800   -9.5299950744   -2.3333466969 226       H    -0.4233146133  -10.4792647793   -2.4467273300      H    -0.7296493022   -8.7528814614   -2.2996582831      N     6.4150035592  -11.7481436357   -2.4943470691      H     6.1088550025  -12.7137257613   -2.4264659840      H     7.4012260765  -11.5248159729   -2.5795591644   54       N    -3.4082153310    1.1773234250    9.7622327616      C    -2.1999841194    0.8906520558    9.3660279537      N    -2.0019512143   -0.1943718625    8.4931559633      C    -3.1046939894   -0.7667775289    7.8395308894      N    -4.3107395048   -0.4879991216    8.2471681414      C    -4.4099402124    0.4121066594    9.2586537157      C    -0.7256954345   -0.6692967204    8.2425480467      N    -0.5211877835   -1.5876124187    7.3125540914      C    -1.5670309826   -1.7425096718    6.5183750404      N    -2.8400456985   -1.4994103783    6.7401838584      N    -1.1194735145    1.6296510849    9.6902049537      C     0.0332147665    1.1046569241    9.3507530913      N     0.3007344778   -0.0991435793    8.8410122244      N    -1.2549704575   -2.0788510017    5.1844943799      C    -0.3118423723   -1.4175157659    4.4674102578      N    -0.6336734000   -1.2210925508    3.2033872271      C     0.1600798930   -0.4806698132    2.4429265943      N     1.3888637381   -0.0228630077    2.9861318034      C     1.6491214742   -0.3092196787    4.3227998918      N     0.8584849149   -0.9692224898    5.0718759161      C     2.2973407915    0.5833762319    2.0896507245      N     1.9437306860    0.9353183424    0.9092872157      C     0.6620411059    0.6271064255    0.5269943914      N    -0.1824715861   -0.1442617469    1.2265145854      N     2.8883781514    0.2539944049    4.7554028911      C     3.8293933701    0.3930344997    3.7115589735      N     3.5961332580    0.6527289405    2.4889663952      N     0.2714471235    1.0934295436   -0.6390065592      N     5.0535822804    0.1184776582    4.2229204171      C     4.9749945362   -0.1893312801    5.5408899041      N     5.9769119244   -0.4253170398    6.2879854512      C     5.7591941175   -0.3866346159    7.6300600686      N     4.5249241637    0.1403593348    8.0645616120      C     3.4224039299    0.3903599667    7.2411949826      N     3.6200194318   -0.1174130506    5.9185870360      C     4.4641240814    0.5925405470    9.4041767912      N     3.3727469833    1.2456555028    9.8127684644 227       C     2.3593231841    1.3357204957    8.9953853941      N     2.3696845015    0.9734059258    7.6447007471      N     5.4330326621    0.3262276882   10.2351450679      C     6.4685729025   -0.3987508375    9.7766640680      N     6.6904332629   -0.6927467739    8.4586536693      N     1.2004309988    1.8581483608    9.5053158854      N     7.3683199574   -0.8195174549   10.6374363724      N    -5.6142050391    0.6009538853    9.7712392194      H     0.9092017416    1.6463579728   -1.2025082254      H    -0.6561743432    0.8663840758   -0.9840732497      H    -2.0814462712   -2.2792423630    4.6190535010      H     5.9208079585    0.1437567401    3.6881292222      H     7.2614995043   -0.5985338456   11.6228690637      H     8.1773109761   -1.3385546798   10.3111744861      H    -6.4000942827    0.0748588548    9.4077148757      H    -5.7417550567    1.2865190666   10.5063238423      2-D melem hexamer dicationic minimum 100       N    -3.6911544307   -0.7154949453   10.4291905582      C    -2.5008563228   -0.7743197132    9.8752523304      N    -2.3458341893   -0.3445375906    8.5615999924      C    -3.4730491163   -0.0935385582    7.7719060374      N    -4.6246891670    0.1311532144    8.3932481855      C    -4.6070107318   -0.0654255064    9.7125837111      C    -1.0738023663   -0.1616541566    8.0305511867      N    -0.0324375872   -0.2115302504    8.8563230044      C    -0.2766922339   -0.8036421923   10.0111984734      N    -1.4241126037   -1.1945812673   10.5342236903      N    -3.3167996138   -0.0954511473    6.4533239832      C    -2.0545623627   -0.1811167533    6.0297485503      N    -0.9416688837    0.0147908295    6.7349957714      N     0.8702929628   -1.0594553928   10.8048792099      C     2.0671698100   -1.4325445826   10.1671566897      N     1.9304137496   -2.0765494103    8.9822164804      C     3.0381327832   -2.4317119221    8.3982052176      N     4.2571779980   -2.3401897106    9.0689390218      C     4.2991316483   -1.6933893840   10.3066652339      N     3.1816817812   -1.1291398175   10.7787235712      N     3.0955943991   -2.9544630926    7.1536803854      C     4.2968998614   -3.2851951580    6.6433032855 228       N     5.4362623172   -3.3584331204    7.2933679680      C     5.4440424582   -2.8970106370    8.5538179312      N     5.4051444509   -1.6789444784   11.0008135937      C     6.4200085847   -2.4133224169   10.5435336548      N     6.5130807718   -2.9512898205    9.3036929224      N     4.2953222430   -3.5955927324    5.3160521545      C     3.2383257430   -3.2581390314    4.5146591052      N     2.2104163855   -2.6807219040    5.1500478633      C     1.1489169022   -2.2341810048    4.4455026228      N     1.2146297470   -2.4331548154    3.0681492116      C     2.3054635962   -3.0869501509    2.4582258134      N     3.3354643612   -3.4863838560    3.2300947397      N     2.3124032936   -3.2895479833    1.1755299714      C     1.2349278217   -2.8473086981    0.4758063069      N     0.1640441160   -2.2017783923    0.9762692703      C     0.1438508555   -1.9848689304    2.2663364326      N     0.1361073629   -1.6442812331    5.0001344852      C    -0.8147366743   -1.1782920458    4.1576232959      N    -0.8666034315   -1.3324137516    2.8478322573      N     1.2489142853   -3.0744187986   -0.8275474020      N    -1.8855930718   -0.5089198793    4.7006731260      N     7.4599704039   -2.6977066997   11.3690425222      N    -5.6668648233    0.4674023330   10.4160534234      H     8.1832196021   -3.2926987090   10.9698801876      H     5.1564510657   -3.9194983147    4.8826110356      H    -6.4931850129    0.7316288253    9.8848077059      H     2.0389960767   -3.5535686467   -1.2424987503      H     0.4716470753   -2.7604584593   -1.3964123852      H    -2.6930285494   -0.4086981494    4.0900688274      C     0.7867815333   -0.9198700801   12.2036131644      N    -0.1584652920   -0.0573206175   12.6486508324      N     1.6162196526   -1.6277715569   12.9254372254      C     1.6819109843   -1.3376934414   14.2318829638      C    -0.2519887907    0.0659387672   13.9405770371      N     0.7071894127   -0.4960529279   14.7833276417      N    -1.2303151920    0.7674255523   14.5544158029      N     2.6334022804   -1.8153440755   14.9882186703      C     0.7191511050   -0.2588272919   16.1707087755      C     2.6833560309   -1.3827988291   16.2551694963      C    -1.2473082350    0.8252971756   15.8978942022      N    -0.3047327281    0.4153479108   16.7156726327      N     1.7148195442   -0.6946968292   16.8971238984      N    -2.3622002116    1.3904340689   16.4415087558      N     3.7792415483   -1.6792887707   16.9913071176      C    -3.4971127425    1.5698888787   15.6994697485      H    -2.4084923801    1.5221004888   17.4488874661 229       N    -4.5686146354    2.0160230429   16.3034228503      N    -3.3905707654    1.2252531078   14.4094785119      C    -4.4803888740    1.2561472448   13.6126381728      N    -4.4535043473    0.8873619814   12.3699272622      C    -5.6409695623    0.8988472596   11.7207924535      N    -6.7995634638    1.3307705722   12.1820547422      C    -6.8330580153    1.7658319621   13.4446356291      N    -5.6542169676    1.7017376592   14.2170896523      C    -5.6910165006    2.1155806925   15.5647117461      N    -6.7923977814    2.5723858070   16.0792291704      N    -7.9331587804    2.2332051334   13.9765716279      C    -7.8776331208    2.6262585255   15.2635610215      N    -8.9877510571    3.1106942912   15.7959117791      H    -8.9870412137    3.4133471668   16.7625540749      H    -9.8270466945    3.1689842255   15.2315541357      C     5.0193315118   -2.0713692858   16.4581081485      N     5.7683983885   -2.8286229545   17.2347911628      N     5.2740554588   -1.6000135158   15.2452833983      C     6.3103235286   -2.1432419091   14.6303792151      N     6.5282192241   -1.9071000420   13.3483222721      C     7.4431117557   -2.6690827197   12.7745872803      N     7.1444015666   -3.0122285519   15.3252578321      C     6.9122229302   -3.2761376000   16.6854864361      N     8.3339693805   -3.4694163654   13.3389196065      C     8.2369518727   -3.6081194310   14.6705933275      N     9.1002114427   -4.3084762216   15.3622324184      N     7.7815834221   -3.9765197412   17.3667524902      C     8.8609543908   -4.4282614043   16.6870812607      N     9.7738326452   -5.0804303230   17.3968066483      H     9.6327360361   -5.2066018164   18.3909721341      H    10.5944038150   -5.4484932396   16.9327174288      H     3.7028385497   -1.5944589638   18.0015152183    230             REFERENCES    231  REFERENCES   (1) Bredas, J. L.; Themans, B.; Fripiat, J. G.; Andre, J. M.; Chance, R. R. HIGHLY CONDUCTING POLYPARAPHENYLENE, POLYPYRROLE, AND POLYTHIOPHENE CHAINS - AN ABINITIO STUDY OF THE GEOMETRY AND ELECTRONIC-STRUCTURE MODIFICATIONS UPON DOPING. Physical Review B 1984, 29, 6761-6773.   (2) Bredas, J. L.; Scott, J. C.; Yakushi, K.; Street, G. B. POLARONS AND BIPOLARONS IN POLYPYRROLE - EVOLUTION OF THE BAND-STRUCTURE AND OPTICAL-SPECTRUM UPON DOPING. Physical Review B 1984, 30, 1023-1025.   (3) Heeger, A. J. Semiconducting and metallic polymers: The fourth generation of polymeric materials (Nobel lecture). Angew. Chem. Int. Ed. 2001, 40, 2591-2611.   (4) Wang, X.; Maeda, K.; Thomas, A.; Takanabe, K.; Xin, G.; Carlsson, J. M.; Domen, K.; Antonietti, M. A metal-free polymeric photocatalyst for hydrogen production from water under visible light. Nature Materials 2009, 8, 76-80.   (5) Ufimtsev, I. S.; Martinez, T. J. Quantum chemistry on graphical processing units. 1. Strategies for two-electron integral evaluation. J. Chem. Theory Comput. 2008, 4, 222-231.   (6) Ufimtsev, I. S.; Martinez, T. J. Quantum Chemistry on Graphical Processing Units. 2. Direct Self-Consistent-Field Implementation. J. Chem. Theory Comput. 2009, 5, 1004-1015.   (7) Ufimtsev, I. S.; Martinez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients, Geometry Optimization, and First Principles Molecular Dynamics. Journal of Chemical Theory and Computation 2009, 5, 2619-2628.   (8) Kaestner, J.; Carr, J. M.; Keal, T. W.; Thiel, W.; Wander, A.; Sherwood, P. DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations. Journal of Physical Chemistry A 2009, 113, 11856-11865.   (9) Dunning, T. H. GAUSSIAN-BASIS SETS FOR USE IN CORRELATED MOLECULAR CALCULATIONS .1. THE ATOMS BORON THROUGH NEON AND HYDROGEN. Journal of Chemical Physics 1989, 90, 1007-1023.   (10) Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP). Chemical Physics Letters 2004, 393, 51-57.   (11) Gordon, M. S.; Schmidt, M. W. Advances in electronic structure theory: GAMESS a decade later. Theory and Applications of Computational Chemistry: the First Forty Years 2005, 1167-1189. 232   (12) Piecuch, P.; Kucharski, S. A.; Kowalski, K.; Musial, M. Efficient computer implementation of the renormalized coupled-cluster methods: The R-CCSD T , R-CCSD(T), CR-CCSD T , and CR-CCSD(T) approaches. Computer Physics Communications 2002, 149, 71-96.   (13) Piecuch, P.; Gour, J. R.; Wloch, M. Left-Eigenstate Completely Renormalized Equation-of-Motion Coupled-Cluster Methods: Review of Key Concepts, Extension to Excited States of Open-Shell Systems, and Comparison With Electron-Attached and Ionized Approaches. International Journal of Quantum Chemistry 2009, 109, 3268-3304.   (14) Piecuch, P.; Wloch, M. Renormalized coupled-cluster methods exploiting left eigenstates of the similarity-transformed Hamiltonian. Journal of Chemical Physics 2005, 123.   (15) Nooijen, M.; Snijders, J. G. COUPLED CLUSTER APPROACH TO THE SINGLE-PARTICLE GREEN-FUNCTION. Int. J. Quantum Chem 1992, 55-83.   (16) Nooijen, M.; Snijders, J. G. COUPLED-CLUSTER GREEN-FUNCTION METHOD - WORKING EQUATIONS AND APPLICATIONS. Int. J. Quantum Chem 1993, 48, 15-48.   (17) Nooijen, M.; Bartlett, R. J. EQUATION-OF-MOTION COUPLED-CLUSTER METHOD FOR ELECTRON-ATTACHMENT. Journal of Chemical Physics 1995, 102, 3629-3647.   (18) Hirata, S.; Nooijen, M.; Bartlett, R. J. High-order determinantal equation-of-motion coupled-cluster calculations for ionized and electron-attached states. Chemical Physics Letters 2000, 328, 459-468.   (19) Bartlett, R. J.; Stanton, J. F.: Applications of Post-HartreeFock Methods: A Tutorial. In Rev. Comput. Chem.; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH: New York, 1994; Vol. 5; pp 65-169.   (20) Gour, J. R.; Piecuch, P.; Wloch, M. Active-space equation-of-motion coupled-cluster methods for excited states of radicals and other open-shell systems: EA-EOMCCSDt and IP-EOMCCSDt. Journal of Chemical Physics 2005, 123.   (21) Gour, J. R.; Piecuch, P. Efficient formulation and computer implementation of the active-space electron-attached and ionized equation-of-motion coupled-cluster methods. Journal of Chemical Physics 2006, 125.   (22) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Ab initio study of solvated molecules: A new implementation of the polarizable continuum model. Chemical Physics Letters 1996, 255, 327-335.   (23) Barone, V.; Cossi, M.; Tomasi, J. A new definition of cavities for the computation of solvation free energies by the polarizable continuum model. Journal of Chemical Physics 1997, 107, 3210-3221. 233   (24) Aono, M.; Nitta, S. High resistivity and low dielectric constant amorphous carbon nitride films: application to low-k materials for ULSI. Diamond and Related Materials 2002, 11, 1219-1222.   (25) Sanderson, R. T.: Chemical Bond and Bond Energy; Academic Press Inc., 1976.   (26) Merschjann, C.; Tschierlei, S.; Tyborski, T.; Kailasam, K.; Orthmann, S.; Hollmann, D.; Schedel-Niedrig, T.; Thomas, A.; Lochbrunner, S. Complementing Graphenes: 1D Interplanar Charge Transport in Polymeric Graphitic Carbon Nitrides. Advanced Materials 2015, 27, 7993-7999.