Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00037520/00001
## Material Information- Title:
- Analysis of the Hâ‚‚âº with Hâ‚‚ reaction using electron nuclear dynamics
- Creator:
- Oreiro, Juan
- Publication Date:
- 1997
- Language:
- English
- Physical Description:
- xiv, 90 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Charge transfer ( jstor )
Chemicals ( jstor ) Dihedral angle ( jstor ) Electronics ( jstor ) Ions ( jstor ) Mathematical independent variables ( jstor ) Physics ( jstor ) Projectiles ( jstor ) Reactants ( jstor ) Trajectories ( jstor ) Chemistry thesis, Ph.D ( lcsh ) Dissertations, Academic -- Chemistry -- UF ( lcsh ) Hydrogen ( lcsh ) Hydrogen ions ( lcsh ) Quantum chemistry ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Summary:
- Electron nuclear dynamics theory.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 1997.
- Bibliography:
- Includes bibliographical references (leaves 82-89).
- General Note:
- Typescript.
- General Note:
- In formulas in title "2"s are subscript "2"s; "+" is superscript "+"
- General Note:
- Vita.
- Statement of Responsibility:
- by Juan Oreiro.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 028637110 ( ALEPH )
38855423 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

ANALYSIS OF THE H2 WITH H2 REACTION USING ELECTRON NUCLEAR DYNAMICS By JUAN OREIRO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1997 In memory of Maria Augusta. ACKNOWLEDGMENTS I would like to thank my advisor, Professor Yngve 6hrn, for his always encouraging, warm, and inspiring attitude, not only in science but also in a personal level. I also wish to acknowledge the support, enthusiasm, and guidance of Dr. Erik Deumens who is always in a good mood despite my repeated attempts to bring chaos onto his neat computer network. Of course, this work would not be possible without the help and the many discussions with the rest of the 6hrnDeumens group. I thank, in particular, Augie Diz, Ricardo Longo, and Hugh Taylor for their invaluable help and exchange of ideas. I couldn't forget to mention Monique Chacon-Taylor, who made things make more sense with her "Margarita nights," albeit I couldn't remember much in the morning. Finally, I would like to extend my gratitude to the other QTP members and personnel for creating a both stimulating and fun environment to do research. iii TABLE OF CONTENTS page ACKNOWLEDGMENTS .......................................... iii LIST OF TABLES ........................................... vi LIST OF FIGURES ................................ viii ABSTRACT .................................................xiii 1 INTRODUCTION ............................................. 1 2 THEORETICAL, TIME-DEPENDENT APPROACHES TO MOLECULAR DYNAMICS ......................................*......... 11 Classical Approach . ............. ............ - 11 Semiclassical Approaches...... ......... ............ 13 Quantum Mechanical Approaches......................... 15 3 THE ELECTRON-NUCLEAR DYNAMICS (END) FORMALISM ........... 17 Method.... ............... ........................ 17 Electronic and Nuclear Descriptions ....o..................18 The END Equations of Motion.............................. 19 4 BIBLIOGRAPHIC REVIEW ..... .......................25 Collision Induced Dissociation .................... . 25 Charge Transfer ................. ............. 28 5 BASIS SETS AND ENERGETICS ................................. 31 Basis and H 31 STO-3G Basis ......................................... 35 Dunning's DZ Basis ..................................... 36 Dunning's pVDZ Basis ...... .......................36 Energetics............................................. . 37 6 CROSS-SECTIONS .......................................... 44 The Collinear Approximation ............................. 45 Collision Induced Dissociation ....................... 48 Charge Transfer ...................................... 52 Full Reaction Space ..................................... 58 Collision Induced Dissociation ...................... 59 Charge Transfer ................................. ..... 65 Comparison with Molecular Dynamics ...................... 70 Final Probabilities ................................. 71 pVDZ basis ........................................... 76 7 CONCLUSIONS .............................................. 79 LIST OF REFERENCES ......................................... 82 BIOGRAPHICAL SKETCH ........................................ 90 LIST OF TABLES Table page Table 5-1: Contraction coefficients (c) and exponents (a) for the STO-3G basis ....................................... 36 Table 5-2: Contraction coefficients (c) and exponents (a) for Dunnin's pVDZ basis, with scaling factor 1.44 .......... 36 Table 5-3: Total energies, in atomic units, for the ground states .................................................. 37 Table 5-4: Bond lengths for the equilibrium geometries in the STO-3G basis, values in atomic units ................... 38 Table 5-5: Bond lengths for the equilibrium geometries in the DZ basis, values in atomic units ....................... 38 Table 5-6: Bond lengths for the equilibrium geometries in the pVDZ basis, values in atomic units ..................... 38 Table 5-7: Energetics in kcal/mol. Energy required to break H4+ into H2+ + H2, H3 + H, and total energy released by CID, respectively ........................................... 39 Table 6-1: CID cros-sections, in A2, for 0.519 eV (CM), DZ basis ................................................. 51 Table 6-2: CT cros-sections, in A2, for 0.519 eV (CM), DZ basis .................................................. 58 Table 6-3: CID cros-sections, in A2, for 4 eV (CM), STO-3G basis ................................................... 64 Table 6-4: CT cros-sections, in A2, for 4 eV (CM), STO-3G basis .................................................. 70 Table 6-5: Probabilities for the final electronic configurations with initial conditions 606060 (4 eV, STO3G). The shaded box displays the highest mixed character and a comparison between END and dynamics on a surface is provided for this case. (a, b are the impact parameters; CT: charge transfer; NR: non-reactive) ................. 72 vii LIST OF FIGURES Figure page Figure 5-1: Available channels and cross-sections (A2) for energies 0.1 to 105 eV (lab) ........................... 33 Figure 5-2: Cross-sections (A2) for the range of energies studied, 1-8 eV (lab) ................................. 34 Figure 5-3: Energetics for the STO-3G, DZ, and pVDZ bases sets ................................................... 40 Figure 5-4: UHF ground-state surfaces for the STO-3G, DZ, and pVDZ bases ............................................ 41 Figure 5-5: Bond lengths diagram for the H3 and H4+ species in the D3h and C2, geometries, respectively ................. 42 Figure 6-1: Degrees of freedom for the H2+ + H2 reaction: a and 0 are the projectile and target tilt angles, respectively; y is the dihedral angle; a and b are the impact parameters in the x and y directions, respectively.45 Figure 6-2: Collinear series, 0.519 eV, DZ basis, impact parameter 5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) .................................................... 49 Figure 6-3: Collinear series, 0.519 eV, DZ basis, impact parameter 3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom viii panel: total Mulliken populations vs. time (all data in a. u.) ... .................. .................................. 50 Figure 6-4: Collinear series, 0.519 eV, DZ basis, impact parameter 0.5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) ........................................................ 53 Figure 6-5: Collinear series, 0.519 eV, DZ basis; the projectile (H2 ) total Mulliken population, as a function of time and impact parameter (all data in a. u.) ....... 55 Figure 6-6: Collinear series, 0.519 eV, DZ basis; the projectile (H2 ) bond length, as a function of time and impact parameter (all data in a. u.) ................... 56 Figure 6-7: Collinear series, 0.519 eV, DZ basis, impact parameter 0.3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) .................................................... 57 Figure 6-8: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 60 Figure 6-9: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 61 Figure 6-10: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 62 Figure 6-11: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees ...... 63 Figure 6-12: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 66 Figure 6-13: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 67 Figure 6-14: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees ....... 68 Figure 6-15: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees ...... 69 Figure 6-16: 606060 END series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) ............... 74 Figure 6-17: 606060 DRC series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) ............... 75 Figure 6-18: 000000 END series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) .................. 77 Figure 6-19: 000000 DRC series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) ............... 78 xii Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ANALYSIS OF THE H2 WITH H2 REACTION USING ELECTRON NUCLEAR DYNAMICS By Juan Oreiro December, 1997 Chairman: N. Yngve 6hrn Major Department: Chemistry The END formalism addresses the solution of the timedependent Schr6dinger equation, treating both electrons and nuclei simultaneously. It differs from the other traditional fully quantum mechanical time-dependent methods in that it does not require a potential energy surface (PES) to carry the nuclear motion. The interaction between electronic and nuclear motion is, therefore, obtained in a transparent way, not relying on PES gradients to obtain the coupling between electrons and nuclei. We analyze the H2+ + H2 reaction at energies below 4 eV using different approximations and basis sets. Other than the choice of initial conditions, form of wave function, and basis set, no constraints are imposed on the system evolution. The nuclei are treated in the classical limit, and the electronic part is described by a single determinantal, unrestricted wave function. We obtain properties, such as xiii Mulliken populations, transition probabilites, and crosssections, from the resulting trajectories. These results are compared with other current theoretical approaches, and with experimental values. The relevance of the electron-nuclear coupling is estimated by comparing the END trajectories with molecular dynamics calculations for selected initial conditions in different basis sets. xiv CHAPTER 1 INTRODUCTION Hydrogen is the most abundant element in the universe. Whenever matter begins to aggregate, the first molecules to form are H2 and carbon, nitrogen, and oxygen hydrides. From dense interstellar clouds to the atmosphere of the Jovian planets, these molecules play a central role, and ultimately determine, the chemistry of the system. On Earth, most of the hydrogen has escaped from the atmosphere. Its reactions do not affect us in a direct way. That did not prevent an early interest in hydrogen chemistry from developing. In particular, the reaction H2+ + H2 has attracted a crowd of experimental and theoretical researchers alike. This enthusiasm is justified. The reaction is rich in features, yet simple enough to allow an understanding of its driving principles. A great number of chemical processes is displayed by H2+ + H2. Long-range Coulomb forces, selective vibrational enhancement, multi-channel competition at low energies, are among the many features present. The ability to predict the outcome of a reaction, based on its initial conditions, constitutes the very core of chemistry. In this context, the meaning of "outcome" has undergone an update in the last few years with the advent of feasible time-dependent methods on the theoretical side and femtosecond spectroscopic techniques in the lab. Again, the H2+ + H2 processes provide an excellent opportunity to put those new tools to the test. At collision energies of about 8 eV (lab frame) or less, only two competing channels are present. The first one, commonly referred to as collision induced dissociation (CID), yields H3+ + H as the product. The other channel involves an electron transfer from H2 to H2+, and is usually called the charge transfer (CT) reaction. Above 8 eV, a multitude of channels is accessible but either CT or CID dominate in the range from 0 to 104 eV. Minor competition is offered by proton and electron production starting at around 8 and 100 eV, respectively, and the atomic Hydrogen Lyman and Balmer series at even higher energies. No experimental or theoretical data is available for energies above 104 eV and results in the lower end are scarce. Each of the species involved in this reaction, either as a reactant, product, or intermediate, has a history of its own. The H2 weak quadrupole transitions, initially measured by Herzberg in the infrared range[l], have become an important tool in identifying astronomical objects[2]. The hydrogen molecule has been widely studied by both theoretical and experimental methods. Its ground state vibrational and rotational spectra were analyzed in detail[3] and data is available for even highly excited Rydberg states[4]. The H2+ ion is equally well known, being a common example in most introductory Quantum Chemistry texts. J. J. Thomson was the first to observe the existence of H3+ in 1912 using a mass spectrometer coupled to his "positive rays" experiment[5]. Its preponderance over H2+ was soon explained by the highly efficient and exothermic CID reaction, which converts the latter to the triatomic ion[6]. Since its discovery, H3+ has puzzled chemists. For many years its structure was totally unknown, leading Eyring to refer to the problem as the "scandal of modern Chemistry[7]". Bohr was the first to publish a theoretical paper about the ion[8]. In the mid and late 30's Eyring, Hirschfelder[9], and colaborators suggested that H3+ exists in a triangular structure. However, it was not until 1973 that an accurate value for its ground state energy was made available[10]. It took a few more years and the, then new, technique of Coulomb explosion to establish its ground state as belonging to a C2v point group[11]. A benchmark calculation employing an extensive Gaussian basis and a Hylleraas-type CI wave function was able to improve on these results by a few tenths of milihartree[12]. In recent years, the interest in H3+ has been revived. This triatomic has been found to be the main protonator agent in, and ultimately responsible for, all known interstellar chemistry[13]. A cascade of reactions in a highly ionizing environment, where free electrons abound, form most of the common hydrides: H3+ + X -+ XH+ + H2 XH+ + H2 -4 H2X+ + H H2X+ + H2 -+ H3X+ + H H3X+ + e- -4 H2X + H -4 XH + H2 where X could be, for instance, 0 or C, and e- is a free electron. It was also suggested that the abundance of H3+ could provide a measure of cosmic ray flux in an interstellar medium[14]. Due to its high rate of formation in hydrogen plasmas[15], its suppression is desirable in nuclear fusion[16]. It is surprising that for a molecule considered the simplest polyatomic and with only one barely stable electronic excited state[17], its spectrum was completely unknown until recently. Only in 1980 the first 15 spectral lines belonging to the V2 fundamental band were identified[18]. There is a wide selection of high-quality PES's available for the H3 system, starting with the one obtained by Carney and Porter in 1974[19] and subsequently used to calculate the vibrational modes of H3+[20]. Their work has been further improved and it is now possible to find even full CI surfaces, such as the Meyer-Botschwina-Burton (MBB) obtained with a (10s,4p,2d) Gaussian basis[21]. Very accurate fittings have been obtained for those PES's. There are no detectable discrepancies between experimental and theoretical spectroscopic results based on those analytical expansions. In fact, these calculations have been used to predict rovibrational lines beyond the limits of the available experimental apparatus[22, 23]. These highly accurate potentials have made H3+ a favorite target for studies. Another participating species, H4+, has been no less controversial. Deemed as non-stable for many years, it was only observed experimentally in 1984(24], after detailed ab initio calculations suggested the structure was viable[2527]. Now, it is known that H4+ is not only stable, but its ground state PES is capable of holding several vibrational levels. Its lifetime is of the order of microseconds. A favorite topic of discussion seems to be whether H4+ can be obtained from the H2 + H2+ reaction. According to current consensus, it can only be generated from breaking heavier clusters. This is based on experimental observations and no reasonable theoretical explanation seems to exist that can support this belief. If H4+ were to be formed, due to its shallow ground state minimum, it would likely be found at low collision energies. No reliable calculations have been performed for collisions energies below 2 eV (lab frame). The difficulties associated with measuring a low-stability charged species might account for its absence in the list of products. This ion exists in a supermolecule-like ground state. The H3+ ion, in an almost C2v geometry, is orbited by a H atom which lies far away aligned to one of the triangle vertices. Conventional detection techniques, such as those used in time of flight (TOF) measurements, would allow smaller charged ions to interfere, a problem known to exist for H2+ and H3+. That would be less likely to occur when H4+ is the result of splitting a heavier cluster, and might explain why it has been detected that way. The possibility of breakage into H3+ + H should not be discarded either, due to the relatively long times required for some detection methods. Unlike the situation for H3+, there are not that many potential energy surfaces available for the H4+ system. The higher number of geometrical conformations to be spanned makes the work of mapping the reaction space both tedious and costly. Accurate fitting is also significantly more complex. The situation is even worse for excited electronic states, required to describe electron transfer properly. That fact might explain the dominant position of the CID channel as compared to CT, which has been borne out in a number of published studies. The inclusion of even simple excited states, such as the ones obtained via self-consistent field (SCF) or singles and doubles configuration interaction (SDCI,) seems to provide good results. It allows for the correct behavior in the CT channel, i.e., an increase of the cross section with the impact energy, and a decrease with the impact parameter[28]. Although models developed within stationary-state frameworks have dominated Theoretical Chemistry in the past, frequently yielding excellent results, a full understanding of a chemical reaction encompasses its evolution in time. A dynamical process, with spectroscopic and structural properties evolving continuously, needs a time-dependent treatment. As successful as traditional approaches to molecular dynamics have been, these did not possess the tools to allow an insight into the reaction mechanisms in terms of basic principles. They have frequently resorted to extraneous assumptions in order to achieve a simple theoretical structure. The inception of time-dependent methods is not recent, but until the last decade they could not be classified as working tools. Both computational and theoretical problems contributed to this delay. For similar reasons, Chemistry has been depicted in terms of mostly adiabatic processes. Both the convenience and the success of the Born-Oppenheimer approximation have led to the belief that the coupling between electronic and nuclear motion was either small enough to be negligible or confined to short segments of a reaction. As a consequence, the Potential Energy Surface picture flourished and it is almost impossible not to rely on it when discussing chemical processes. That well established view has been challenged more intensely in the last few years. A particular example being the study of charge transfer. The refinement of existing experimental techniques and the introduction of new ones made more apparent the deficiency of older models. The increased computer power allowed the development of new ones. The Electron-Nuclear Dynamics (END)[29] theoretical framework, goes beyond the above restrictions. It allows for a full and instantaneous interaction between electrons and nuclei, while developing approximations to the time-dependent Schr6dinger equation. The END equations have been coded into the software package ENDyne and will be discussed in some detail later. The analysis of the resulting evolution is done by the program Evolve, which was revised for this study. In its present incarnation, END employs a single determinant built from spin-unrestricted orbitals with complex coefficients to describe the electrons. The nuclei are considered in the limit of narrow wave-packets. The system evolution proceeds without the need for a PES. Not only does that save computer time but also eliminates inaccuracies due to the fitting procedure. It is a one-step, general process. Traditional methods involve three steps: calculation of the PES, fitting of the surface, and carrying out the dynamics. (The second of them frequently is property or application oriented.) With END it is possible to calculate systems for which no PES exist and retrieve most properties, as a function of time, from the resulting wave function. Despite being investigated by an unusually high number of experimental techniques and theoretical approaches, the reaction of H2+ with H2 is still a challenge. The fact that the system is the simplest one involving a charged and a neutral molecule could not explain alone the intense interest it has motivated. The existence of a high quality PES, and analytical fittings for those, has certainly helped to attract part of the crowd of dynamicists eager to test their approaches in a reaction with more than three nuclei. Reliable PES for systems with four or more atoms are still an exception at this time. Since most of the available approaches for molecular dynamics rely on one or more PES's, systems for which such surfaces exist are always a preferred target. The vast number of experimental results available for this reaction has also contributed to intensify the efforts of studying it using a variety of theoretical techniques. A fully quantum-mechanical, three-dimensional study was published only recently on the subject[30]. The authors calculated both the reactive and the charge transfer cross sections for two values of collisional energy (0.252 eV and 0.519 eV). They used only 14 randomly selected initial geometries due to the high computational demand of the problem. Their results agree remarkably well with experimental measurements for the reactive cross section, but not for the charge transfer (CT) one. Several reasons come to mind to explain the failure in predicting the CT cross section. It might be an intrinsic limitation of their approach (Infinite Order Sudden Approximation - IOSA), an inaccuracy due to the rather limited sampling of the reaction 10 space, a basis set artifact, a feature of the potential energy surface or fitting employed, or a combination of any of those and other factors. Also, even for the reactive cross section, there are many details that are not clear or could be explored more deeply. In order to do that, we decided to carry out calculations on the H2+ with H2 reaction using the END approach. CHAPTER 2 THEORETICAL, TIME-DEPENDENT APPROACHES TO MOLECULAR DYNAMICS Classical Approach The classical approach to molecular dynamics consists of two steps: the determination of the potential energy surface (PES) and the numerical tracing of the collision trajectories on that surface, following classical equations of motion. The latter is a rather trivial procedure, especially with the fast computers available now. The generation of the PES and its fitting is the crucial part of this method. Despite the fact that only one PES is required, the accurate representation of the potential topology is still a demanding task for most systems containing more than three atoms. In addition, processes involving multiple PES's are not adequately treated within this framework. There is a great number of processes that are not adiabatic, among them, charge-transfer, electronic energy transfer, collisional quenching, spin-forbidden reactions, and radiationless transitions. To accurately depict such phenomena one has to include not only all the PES's involved, but also the couplings between those surfaces[31]. Also, processes where quantum effects are noticeable cannot be correctly described, those include systems for which tunneling and zero-point energy are important. The usual way to obtain the potential is through ab initio calculations, although semiempirical and even empirical methods have also been employed[32]. The discrete set of energies must then be converted to a more suitable analytical form, in order to allow for the rapid evaluation of the derivatives of the potential energy. However, the large number of nuclear configurations involved in a molecular collision precludes the calculation of a complete surface, even for systems for which accurate energies can be obtained. The fitting must be based on a limited set of geometries, usually chosen by convenience or symmetry considerations. According to McLaughlin and Thompson[33]: "-The calculation of surfaces remains somewhat of a stumbling block to a completely ab initio approach to classical dynamical studies of chemical dynamics." The fitting of PES's has become a favorite topic of interest by itself and some examples can be found in the study on 2D cubic spline interpolation by Sathyamurthy, Kellerhals, and Raff[34], the scattering claculations for He + H2+ by Sathyamurthy, Rangarajan, and Raff[35], or the H + H2 reaction using classical trajectories by Gray and Wright[36]. Semiclassical APProaches This class of approaches comprises methods that employ simple approximate procedures to extend the validity of the classical models. They provide a better description of the molecular dynamics in the sense that one is no longer limited to purely adiabatic effects. In their most common implementation, two PES's are used. The coupling between the two surfaces is normally introduced in an ad hoc way. The evolution is carried out until a certain point of the trajectory is reached, where diabatic effects are believed important. At that point, a somewhat arbitrary procedure is used to trigger the surface hopping. How well the technique describes the system depends on how accurately the transition point is chosen and how good the guess for the coupling between the PES's is. The same limitations concerning the calculation and fitting of the potential surfaces apply, with the extra complication that now more than one surface is needed and they will not be necessarily of the same quality. One popular implementation of the semiclassical methodology is due to Miller and George[37] and it's aimed at describing electronic transitions in low energy atomic and molecular collisions, or involving quantized degrees of freedom. Their framework is based on the integration of Hamilton's equations using complex-valued time. Both coordinates and momenta are allowed to take on complex values and the transition happens when a trajectory crosses a branch point (corresponding to an adiabatic surface intersection). The approach relies on the Feynman propagator and involves the construction of complex-valued time contours to find a hopping point. It is assumed that at low energy the dynamics is essentially a classical single-surface phenomena, with transitions being isolated events that happen in localized points of space and time. Another common technique was devised by Tully and Preston[38] and called the "surface hopping" method. It differs from the Miller and George approach in that the momentum is not required to be a smooth function of time. Rather, the momentum changes abruptly when an avoided crossing is reached, to allow for the propagation on the different PES. The main advantage in this case is that one does not have to deal with the complex-valued trajectories. Both techniques have been shown to be equivalent, provided there are only two surfaces involved in the crossing, the splitting of the surfaces may be expressed as a function of a single variable, and the trajectories have no turning points close to the intersection ("high energy limit")[39]. There's a great number of papers applying semiclassical methods, particularly surface hopping and its variations, to light systems. The semiclassical approach is sometimes called quasiclassical in these studies[40-45]. Quantum Mechanical Approaches This class of methods, unlike the previous ones, rely on fully quantum mechanical equations of motion. Coupling between different surfaces arise naturally and do not depend on predefined transition schemes. Presently, there are two main approaches in the literature for quantum mechanical evolutions on potential surfaces. The first one was proposed by M. Baer and coworkers and it is referred to simply as the quantum mechanical (QM) method. It's based on what the authors call "the negative imaginary decoupling potentials" (NIPD) to convert multi-arrangement channels into inelastic singlearrangement systems, in order to reduce the configuration space. The inelastic collision is then treated by means of the infinite order sudden approximation (IOSA). That method has been applied to a number of collinear and three-dimensional systems[46, 47]. Baer and Ng[30] have published the only paper where the full dimensionality of the H2+ + H2 reaction was taken into account in a quantum mechanical level, where they have used 28 randomly selected configurations (initial relative molecular orientations). The second approach in use is usually referred to as exact quantum mechanical methods[48, 49]. They are "exact" in the sense that the error could, in principle, be made as small as required. This is in practice frustrated by computational constraints. In order to represent the system wave function, a basis of a numerical grid and weights is required. The smaller the interval between the grid points, the better the accuracy achieved. The total size of the grid is also important, since boundary effects may compromise the results. It is frequently not possible to generate the grid with dimensions large enough to describe the reaction going on. Wave packets can approach the rim of the grid and impart non physical effects to the dynamics. To minimize that, one has to apply "cushion" potentials to simulate a natural exit for the reactants. Exact methods give excellent results when they can be applied. Currently, the larger systems studied with this technique consist of six atoms. CHAPTER 3 THE ELECTRON-NUCLEAR DYNAMICS (END) FORMALISM Method The END formalism addresses the solution of the timedependent Schr6dinger equation, treating both electrons and nuclei simultaneously. It differs from other traditional fully quantum mechanical time-dependent methods in that it does not require a PES to carry the nuclear motion. The interaction between electronic and nuclear motion is therefore obtained in a transparent way, not relying on PES gradients to obtain the coupling between electrons and nuclei. There are several advantages in not depending on PES. First, there is no need to calculate such surfaces per se. For mid-sized and large systems, obtaining a good PES is a major task. With several degrees of freedom, it can easily be the most time-consuming part of the calculation. Commonly, the groups who work on dynamics are not the ones who generate the PES. That creates two problems. First, the dynamicists don't have control over the PES quality. Second, they are limited to whatever data is available. The actual calculation of the surface is just one of the steps, however. The data is given in a set of points that must be interpolated and extrapolated to provide the necessary surface. The quality of the interpolation is critical to obtain a good representation of the dynamical event. So, another source of errors is introduced, over which one has little control. Besides, a chemical process frequently requires the use of several surfaces, corresponding to different electronic states. That is especially true at higher energies. Finding one highquality PES is difficult, finding several of them for the same system is almost a unique event. END is based on the Time-Dependent Variational principle (TDVP)[50], with a wave function parametrized in terms of Coherent States[51, 52]. The main advantage of using Coherent States is that one gets simplified equations, non-redundant parameters in the definition of the electronic state, and it is easier to interpret the results. This kind of parametrization also guarantees that all wavefunctions of a particular family can be reached, should the dynamics so require. The END formalism has been implemented in the program ENDyne[29] for a number of wave functions, with some restrictions and approximations which will be discussed below. Electronic and Nuclear Descriptions As it is now, ENDyne is limited to treat nuclei as classical particles. This corresponds to the limit of infinitely narrow gaussian wave packets and defines their time-dependent variables as their positions and momenta and retains all nonadiabatic coupling terms. The electronic part is expressed in terms of a truncated Gaussian basis and is restricted to single-determinantal wave functions. The formalism for a multi-determinantal electronic description is developed and under implementation[53]. That means that we can set the initial conditions for an evolution as a set of generalized coordinates for the nuclei and an electronic state, the latter expressed in terms of complex spin-orbital coefficients. The only additional data needed are the mass of the nuclei and their atomic number. There are three different choices for the electronic wave function presently available: a restricted singledeterminantal wavefunction, with either pure-spin orbitals or mixed-spin orbitals. A two-level model, for both classical quantum (Gaussian) nuclei is available. In this work we use the unrestricted description with the pure-spin option. The END Equations of Motion Here we don't give a detailed derivation of the END equations but just show the structure of them and where they come from. We are going to outline the expressions built in the ENDyne code. Instead of considering the most general case, we will limit ourselves to classical nuclei and assume that the electronic part of the wave function can be described by a complex, spin-unrestricted determinant. A rigorous and more general account of the formalism behind the method can be found in the references listed previously. In our case, the atomic spin orbitals (AOs) are given by a set of Gaussians, ()lzm ex4 -b( -1] (3.1) with an appropriate spin part, a or 0, attached to them. We can define a general, complex spin orbital for the electrons as VO lil Ui(ci ,(3.2) where the total number of AOs, K, consists half of a and half of P orbitals. For a Coherent State wave function ), defined by a set of parameters 1}, we can write the first variation of the quantum mechanical action functional, A, as M= { ( - (I -( ) )'dt (3.3) and the solutions we look for are the ones for which 6A= 0, according to the TDVP (principle of least action). By introducing the functions 4 )=[Hk)/[I) and for the total energy and overlap, respectively, and setting h =1, after several algebraic steps we can rewrite the first variation of the action functional as .5A~f _2_In_ dE fr+ .(,d 2 InS f }ldt =0, (3.4) dCdf C'-pdIdf where the denotes differentiation of the complex parameter with respect to time. The solution of the TDVP expression will then be the dynamical equations given by where we define the general element of the hermitian matrix C as C. =d2 In There are a few ways one can choose the electronic basis dependence on the nuclear motion. Three of the most useful choices are called the nonfollowing basis (NFB), the static following basis (SFB), and the dynamic following basis (DFB). In the limit of a complete basis set, the three types are equivalent. However, actual calculations are carried out in truncated basis sets, with a finite number of orbitals. The three bases will not necessarily yield the same results under those circumstances. A NFB is one for which the basis grid is fixed in time and space, the electronic orbitals do not move at all during the system evolution. The equations of motion assume their simplest form in this case. This is the ideal choice for carrying the algebraic manipulations and interpreting the equations, if one can assume completeness. A SFB is one for which the basis are centered on the average position of the nuclear wave packets, or the nuclear coordinates in the limit of infinitely narrow wave packets. The AOs do not carry an explicit term containing the orbital velocity in these cases, they just follow the nuclei. A DFB is one for which the AOs are defined in such a way that they have an intrinsic velocity and their positions change with time. The velocities are embodied in factors with the properties of a plane wave and are usually referred to as electron translation factors (ETFs)[54]. Those terms are known to play an important r6le in high energy collisions. The general form of an AO (omitting spin), including the ETF would be u xyz ex)2 +imev.0 .. (3.6) Saying that the three bases above are equivalent in the completeness limit amounts to stating that their equations yield the same dynamics. Their metrics are the same. They can be therefore related by a transformation that does not change the metric. Such transformation is referred to as a symplectic one (or canonical, when the metric only contains zeroes and ones). At this point, ENDyne does not include ETFs in its equations, but a SFB representation is used. Ordinary chemistry takes place at relatively low energies. We therefore don't expect ETFs contributing significantly to the study proposed here. The omission of such factors have a negligible effect in the final results. (Their effect is taken into account approximately by the basis chosen, since the MO coefficients are complex.) For NFB, the equations of motion with real nuclear coordinates, R, and momenta, P, become (z are the Thouless parameters for the electrons),[55] yielding the individual equations, E c / 4 - (3.7) dR Pk =-VRkE, and (3.9) P Rk= k 3.10) Assuming a complete basis, we can obtain the equivalent result for a SFB, after a symplectic transformation, as iCI dEf, -0 -iCR o " ' C -' R CR 3 .R1 where I is the identity matrix, denotes a Thouless parameter in the SFB and the following definitions were made, CX, = dkdXk' R'=RP'=P and cxdyd2 =nS(-2,RPiRP . (3.13) dX~dY, 'RP= The individual equations in vector form for this case are therefore, dE C! + CR R, = -i d (3.14) k + 2ImCi j-C,,, R, =-VRE (3.15) k (3.16) Mk Those are the equations coded in the ENDyne program package, where a summation over the nuclear coordinates label, 1, is implied. CR and C. are simply the so called nonadiabatic correction terms, that in the END theory appear in the dynamical metric. CHAPTER 4 BIBLIOGRAPHIC REVIEW Collision Induced Dissociation The reaction H2+ + H2 -- H3+ + H has been extensively studied over the years. Its importance in interstellar[56J and atmospheric chemistry[57], as well as its suitability for rigorous theoretical treatment, account for part of the interest. Collision induced dissociation (CID) between ion and molecules are part of a class of reactions generally referred to as collisional energy transfer. However, they differ from the most general case in two important aspects: first, the amount of energy transferred is large enough to comprise a significant part of the translational energy; second, the reaction mechanism may be complex, involving several electronic surfaces and competing reaction channels[58]. Experimentally the H4+ system has been studied by crossed beam[59, 60], mass spectrometric[61, 62], merged beam[63, 64], photoionization[65], ICR[66], and ion beam-gas cell[67] techniques. These experiments showed the reaction to be exoergic by 1.7 eV and, contrary to initial beliefs, proceeding through a direct mechanism[63, 59, 60] where no barrier is involved[63, 64] and confirming the instability of H4+ with respect to the H3 + H products[68]. In addition, its total cross-section was found to be very large at thermal energies, and to fall off rapidly with increasing collision energy[63, 64]. In principle, there would be two ways for the collision induced dissociation to proceed. The H3+ product could be formed by either proton or by a H atom transfer. The first mechanism would involve the breakage of the H2, ion (requiring approximately 2.6 eV,) while the second would entail the rupture of the H2 molecule (4.5 eV.) Energetics alone would make it seem more likely to have a predominance of the proton transfer. Trajectory surface hopping calculations have shown this to be true[69], the H2+ rupture not only dominates but it seems to be the exclusive mechanism for the reaction. Attempts to confirm this experimentally, by using isotopic substitution, have yielded a very different picture. It was found that both mechanisms occur, and have comparable crosssections, apparently due to intensive charge transfer between the reagents[63, 69, 60, 59]. Photoionization experiments have shown that at collision energies below 1 eV (center of mass,) the H3 formation can be inhibited by vibrationally exciting the H2 ion[70, 65]. At higher energies, the opposite seems to hold. The reason for this behavior has been explained in terms of the interaction of several potential energy surfaces[60, 59, 71], putting a rigorous treatment of the problem out of reach of most current theoretical methods. Such interaction appears to be the underlying cause for the competing mechanisms found in experiments involving isotopic substitution. The (H2 + D2) ground PES could be entered by two alternative paths, one belonging to H2 + D2, the other to D2 + H2. At infinite separation each channel would have its own PES, but as they approach, making the ion and the neutral bond equal, an avoided crossing puts them on the same ground PES. The barrier between the two PES decreases until it becomes smaller than the reagent zero-point energy, which happens at 8-10 Bohr. At this point charge hopping between the two channels sets in[72, 59, 60, 27], accounting for the formation of both D2H and H2D . The multiple surface nature of the H4 system also further complicates the comparison between experimental and theoretical results. Most experiments involving electron bombardment are likely to excite the H2 ion, populating higher vibrational states. Such experiments include crossed and merged beam studies[63, 64], and their results should be considered as involving the all the accessible potential energy surfaces[73]. The CID channel was also explored with the molecular beam technique, using the enhanced four-photon ionization to prepare the H24 ion in selected vibrational states distributions, for collision energies ranging from 1.5 to 5.3 eV (c m.)[74] The resulting data included the state-selected integral and differential cross sections. On a different experiment, using the merging beams technique, Lees and Rol[64], obtained the reaction cross sections for several isotopic substitutions, at energies below 10 eV. Charge Transfer When a charge transfer (CT) process involving an ion and a molecule is exothermic, it is usually the dominant product channel for the reaction. Due to the large CT cross-sections, these reactions play a major role in both relaxation processes and reaction kinetics in ionized gases. The understanding of this class of processes is of special interest to gas discharges, lasers, flames, and controlled thermonuclear fusion. Although present in low concentration in flames, the ions are responsible for the thermal charge exchange reactions which create free radical chain carriers. For collisions between molecular ions and molecules, where spatial orientation provides an extra challenge in generating the interaction potentials, the availability of reliable potential energy surfaces has been one of the major stumbling blocks for theoretical treatments aiming at obtaining accurate cross-sections. For the H4+ system, a number of PES has been made available, mostly for the ground state[75, 27, 71, 72, 28, 25]. On the experimental side, the challenges have not been smaller. Most of the studies done involve the measurement of cross sections as a function of the kinetic energy, and preparing the H2 ion by electron ionization[76, 77, 61, 7886]. The H2 ions formed by electron ionization usually have a broad distribution of vibrational states. Experimental results seem to indicate that a strong dependence of the cross section on the vibrational population of H2' is to be expected[76, 77]. This finding, with the difficulty in obtaining reliable PES, has hindered the efforts to obtain a sound comparison between experiment and theory. For that reason, in recent years, the preferred empirical approach to prepare the system in well characterized distribution of states has become photoionization[87]. It is worth pointing out that most of the available reliable experimental results concentrate on higher collision energies, in particular above 200 eV[88]. One explanation might be the suggested poor resolution of the time-of-flight (TOF) analysis at low energies. Due to inelastic charge transfer channels, the slow H2+ products can be scattered by more than 10' away from the initial neutral H2 beam direction[89]. The collection efficiency in the TOF analysis for product ions scattered at wide angles is poor. In addition, inelastic charge transfer can also lead the arrival of the H2+ product ions to be spread in time. As a result, for the reaction H2+ + H2, the TOF does not provide a reliable detection method at energies below 8 eV (center of mass), since it cannot distinguish between the H3+ and the H2+ ions. To circumvent those problems, Liao, Liao, and Ng[87], develop a ion-molecule reaction apparatus which combined the crossed ion-neutral beam method, high resolution photoionization mass spectroscopy, and charge transfer detection. Using this technique, they obtained the total charge transfer cross-section, as a function of the vibrational state of H2* (for states up to 0'= 4), at energies ranging from 8 to 200 eV (center of mass.) It was also found, that the low-lying rotational states of the H2 ion had no discernible impact in the outcome of the reaction, at least when the ion was in its ground vibrational state and the collision energy was below 4 eV (c.m.) A strong dependency on the vibrational state was observed, with its maximum around 2 eV. It however decreases as the collision energy is reduced, contradicting the theoretical results trajectory surface hopping[90]. In a later experiment, with improved resolution, Liao and Ng[91] extended those results finding a good agreement with the semiclassical energy conserving trajectory calculations[92] at 16 eV, but with significant deviation at lower energies. Other experiments include the use of pulsed synchroton radiaton and photoelectron-photoion coincidence to select the internal state of the H2 ion, enabling the measurement of relative cross sections for states with H2+( ï¿½= 0-10) at 16 eV. The coupling among several vibrational states of the reactants and products was also observed for the CT channel in the energy range 8-400 eV, using the crossed beam technique[93]. CHAPTER 5 BASIS SETS AND ENERGETICS BasLis ndH4 + As mentioned previously, only the two channels, CID and CT, are present at the energy range studied, in addition to the non-reactive one. A compilation of all the available theoretical and experimental data to date on the H2 + H2 reaction was provided by Phelps[94], and is shown in figure 5-1. A detail of the region analyzed here is presented in figure 5-2. Before we present the results, we will discuss briefly the basis sets employed and our reasons to choose them. Due to the number of degrees of freedom involved, and the intensive computational nature of the END, our main concern is to find a basis that allows us to obtain a qualitative picture of the reaction, while minimizing the computer demand. Unlike the calculation of stationary states, in dynamics there is no common "feeling" of which basis will yield the desired results. That is especially true for END. In order to quantify how well each of the basis sets performs, we ran several calculations using an SCF wave function to minimize the energy of the projectile, the target, and some relevant intermediates. It's known that H4+ is a stable molecule with a lifetime of about 10 ps, capable of supporting several vibrational levels[24]. It's also known that the SCF wave function can correctly describe that stability, yielding a ground state of C2v symmetry[26]. Although it's generaly accepted that H4+ cannot be formed by collisions between smaller molecules, but instead should arise only from breaking heavier fragments, an accurate depiction of the reaction energetics demands that at least the relative energies between the species involved are reasonably well reproduced. We have chosen three different sets. Since our goal is to achieve a qualitative picture, we used the system energetics as a guide to initially gauge the performance of each. Here, we compare the results for the STO-3G basis, a double-zeta basis[95], and Dunning's pVDZ basis[96]. In a related calculation, the stability and geometry of H4+ was analyzed by Wright and Borkman[25]. Since they do not make clear what scaling factor they have used in their double-zeta case, if any, we have chosen to carry out the tests with the factor suggested by Dunning (1.44). To study the collisional process we use a factor equal to 1.54 for both the STO-3G and the DZ basis. The total energy at key points, involving intermediates, entrance and exit channel products was used to compare their effectiveness. I v -13+ +H 0 CT A VE * H*+H t' Lyman Balmer , e" I 0.01 0.001 0.0001. E Ib (ev) Figure 5-1: Available channels and cross-sections (A2) for energies 0.1 to 105 eV (lab) Elb'(eV) = 2 EO 4 5 6 7 8 -0- H3++H -0- CT --A- VE --- H++H * Lyman * Balmer -4- e" Figure 5-2: Cross-sections (A2) studied, 1-8 eV (lab) for the range of energies 2 3 I 2oI N I5 I 1M I 10 I 3 IV .111 III I HVII~ _ III N II The three sets analyzed were the STO-3G, and Dunning's DZ, and pVDZ. Both the number of trajectories to sample the reaction space and the time necessary to simulate a single collision allowing for product separation are high. Since the END effort scales as N4, where N is the number of basis functions, the size of the basis is a major concern. In practical terms, we have limited this study to the DZ basis and the pVDZ is used only as a reference for selected initial conditions. The PVDZ basis is our benchmark set, being very successful in previous similar calculations on smaller systems. It also allows us to see how the introduction of polarization functions affect the results. A complete calculation of the reaction using this basis is not feasible due to the computer load imposed, requiring nearly 40 times the effort of a minimal basis. The DZ basis falls in between the previous two and gives an estimate of the importance of diffuse functions. STO-3G Basis This is a standard minimal basis, where a linear combination of three Gaussians has been fit to a single Slater-type function. The is orbital on each center is then represented by one of the resulting functions. The optimal fit for a Hydrogen atom, according this procedure yields the basis listed in table 5-1. Table 5-1: Contraction coefficients (c) and exponents (a) for the STO-3G basis. c a is orbital 0.44463 0.109818 0.535328 0.405771 0.154329 2.22766 The fitted exponents above were then modified to account for the tighter Is orbitals in the reactants as compared to the atom. The scaling factor, 1.5376, was obtained from the Slater exponent, 1.24. Dunnina's DZ Basisf95] The DZ basis has been used successfully by Wright and Borkman to study the stability of Hn+ clusters, with n as high as 9. The resulting energies and geometries were comparable to a larger [5s,2p] basis. Dunning's OVDZ Basisf96] Table 5-2: Contraction coefficients (c) and exponents (a) for Dunnin's pVDZ basis, with scaling factor 1.44. c a is orbital 0.01969 18.7344 0.13798 2.82528 0.47815 0.64022 2s orbital 0.50124 0.17568 2p orbital 1.00000 0.72700 In order to ascertain the performance of the above sets, we assume the reaction proceeds through a simple path going from the reactants, to the intermediate, to the products. The differences in energy between these points is compared in the three bases. The pVDZ results are considered as the reference values. A geometry optimization, using a UHF wave function, was performed for each of the species involved. The calculations were not constrained to any particular point group, with the exception of the ones for the H4+ ion. In this case, the optimization was limited to the C2. symmetry to accelerate convergence. This should not impose a major constraint, since the H4 ground state is known to belong to the C2, geometry. In all sets, H3+ converged to the D3h point group, as expected. Table 5-3 contains the resulting total energies. Table 5-3: Total energies, in atomic units, for the ground states STO-3G DZ pVDZ H -0.466582 -0.497637 -0.499278 H2 -0.582697 -0.586387 -0.600286 H2 -1.117506 -1.126658 -1.128746 H3 -1.246860 -1.275824 -1.294005 H4+ -1.743480 -1.777434 -1.799292 Table 5-4: Bond lengths for the equilibrium geometries in the STO-3G basis, values in atomic units. R R" Rff H2 + 2.00421 H2 1.34592 - H3+(D3h) 1.82718 1.82718 1.82718 H4+(C,2 ) 1.89523 2.24201 1.57174 Table 5-5: Bond lengths for the equilibrium geometries in the DZ basis, values in atomic units. R R" R"f H2 + 1.98507 H2 1.38142 - H3+(D3h) 1.60570 1.60570 1.60570 H4+(C2v) 3.30625 1.64697 1.57097 Table 5-6: Bond lengths for the equilibrium geometries in the pVDZ basis, values in atomic units. R R' R" H2 1.98053 H2 1.41343 - H3+(D3h) 1.68039 1.68039 1.68039 H4+(C2) 2.096235 1.74943 1.62159 In tables 5-4, 5-5, and 5-6 we list the equilibrium geometries for the above species. They essentially agree with values listed in the literature and confirm H4' as a stable species. In table 5-8 we show the energies required to break H4' into H2' + H2 (either CT or non-reactive), H3+ + H, and the total amount of energy released by CID. The energy barriers for the three sets is shown in figure 5-3 as the minimal energy path for the exit channel of a bond-breaking reaction (H2' + H2 -4 H4+ -+ H3+ + H) proceeding through a C2v H4+ intermediate, and the reactants UHF surfaces in figure 5-4. Table 5-7: Energetics in kcal/mol. Energy required to break H4 into H2 + H2, H3 + H, and total energy released by CID, respectively. H2 +H2 H3 + AE STO-3G 27.16 18.85 8.31 DZ 40.40 2.49 37.91 pVDZ 42.84 2.36 40.48 40 H2*+H2 H4+ (C2) H+ H3 (DM) -1060 I -1070-1080. _ -10e0-_ _ _ _ P-1100, -1110 -1120 S DZ pVDZ -1130 Figure 5-3: Energetics for the STO-3G, DZ, and pVDZ bases sets. Some conclusions can be drawn from the above numbers. The stability of the H4' species is a function of the basis size. As the basis become larger, going from the STO-3G set to the pVDZ, H4+ becomes more stable with respect to the entrance and the CT channels. The opposite happens for the CID channel. Increasing the basis makes it easier to have a collision resulting in dissociation. H2 (STO-3G) H -2 (STO-3G) HA2 (DZ) & H2' (DZ) o H2 (PVDZ) * H2+(PVDZ) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Internuclear Distance (Angstron) Figure 5-4: UHF ground-state surfaces for the STO-3G, DZ, and pVDZ bases 41 UHF PES for the Reactants in Their Ground State Figure 5-5: Bond lengths diagram for the H3 and H4 species in the D3h and C2, geometries, respectively As a consequence, in the STO-3G basis, it is difficult to break the H44 and, when it happens, the channels H2+ + H2 and H3 + H will be more balanced than in larger basis. Between CT and CID, the later has a better chance. The DZ and pVDZ basis seem to perform about the same, at least in terms of energy differences. It is interesting to note that the equilibrium geometry of H4 in the STO-3G basis resembles the two reactants aligned perpendicularly. For the other two basis, the geometry is similar to an H3+ with an H atom aligned to one of the vertexes and orbiting far away from the H3+ core. That weak bond would make it more feasible to attain the CID channel. L-L3 (D3h) R R R ~H4+ (C,2) It is clear that increasing the basis set reduces the H4+ stability. The percent relative difference between the DZ basis and the pVDZ one is just 5%. The DZ basis overshoots the exit barrier by about 5.5% of the pVDZ value, while the minimal basis requires over 7.5 times the energy when compared to the DZ number. Our results disagree with the findings of Wright and Borkman in that they report an increase in stability with the size of the basis. Their double-zeta result gives a bond-breaking energy of approximately 1.3 kcal/mol and a value of 3.2 kcal/mol for a [5s,2p] basis. That's too much of a discrepancy to be accounted for by differences in the scaling factor. According to our numbers, the DZ should perform comparably to a pVDZ basis, at least in energetic terms. The minimal basis, while displaying the correct qualitative features of the larger sets, might lead to extremely stable H4+ intermediates, especially at lower values of collision energy. CHAPTER 6 CROSS-SECTIONS The full treatment of the H2 + H2 collision involves six internal nuclear degrees of freedom, corresponding to a choice of five parameters for the initial conditions. Consider the reactants in the lab frame, and take the z-axis as being the direction of the reaction. A single trajectory can be described by two tilt angles, a dihedral angle, and two impact parameters. The tilt angles refer to the inclination of each of the reactant bonds with respect to the z-axis. The dihedral angle is defined by planes formed by each of the reactant bonds and the z-axis. The two impact parameters are obtained from the distance between the center of the reactants and the z-axis. We show a diagram of the above in figure 6-1. In order to obtain the total cross-sections, it is necessary to run several trajectories. Each trajectory is defined by a choice of initial conditions, which include the assignment of a value for each of the degrees of freedom. With five degrees of freedom, the number of possibilities increases quickly, making it important to devise a way to sample the reaction space without incurring in a prohibitively large calculation. For example, a relatively coarse grid with 6 points per degree of freedom, would involve 7,776 trajectories. To decrease the mesh size to half of that would require 248,832 trajectories. Below, we explore a few possibilities to keep the calculation manageable. %II x . 4 / Figure 6-1: Degrees of freedom for the H2 + H2 reaction: a and 1 are the projectile and target tilt angles, respectively; y is the dihedral angle; a and b are the impact parameters in the x and y directions, respectively. The Collinear Approximation The H2 molecule, and to some extent the H2 ion, possess an almost radial potential. This spherically symmetric potential suggests that the reactants could, in principle, be treated as atoms, which would reduce the problem to just one dimension. In fact, that idea has been used successfully in the past at energies of the order of 300 eV[97]. Although, it is more likely that the asymmetry of the potential plays a ~.4- Y 'K I I I / I ~JL4 I I I I I L I r8le at low-energies, there is a possibility that this approximation could yield acceptable results concerning the total cross-sections. Since the calculation of the total cross-section involves averaging over all the degrees of freedom, any asymmetry could be averaged as well. To test that possibility, we performed calculations at very low energies, 0.519 eV, assuming both reactants as radially symmetric, and taking into account just one of the degrees of freedom. The reduction in the number of trajectories, by considering just the impact parameter, is especially welcome at very low energies. To avoid unnecessary oscillations in the electronic populations, due to the overlap of the reactants' orbitals, the projectile and the target must be kept at a reasonably large separation. The higher the frequency of the oscillations, the harder the integrator has to work to keep track of the evolution, leading to extended computation time. This problem is somewhat alleviated for small basis sets, where the overlap is only significant at close range. For the H2' + H2 reaction, the total elimination of the interaction between the fragments is not possible, regardless of the amount of separation. The Coulomb interaction, due to the charged species, is always present, even at asymptotic distances. In this series of calculations, we use an initial separation of 20 a.u. between the centers of mass of the reactants. Since we are assuming radial symmetry for both reactants, we must choose an initial spatial orientation for them which will be the same for every trajectory. It is important to notice that although we choose the initial orientation, the reactants are free to reorient as the reaction proceeds. This dynamically changing evolution is significantly less restrictive than some previous treatments of this system. A case in point is the full three-dimensional IOSA methodology[30]. Despite taking into account all the available degrees of freedom, the IOSA evolution makes no allowance for any re-orientation of the reaction participants[98]. For all practical purposes, the reactants are frozen in their initial orientations throughout the reaction. In our understanding, this is far more restrictive than a choice of fewer degrees of freedom, in particular for charge-transfer processes. In our calculations we choose to put both projectile and target initially to lie parallel to the reaction coordinate. For that reason, we name this series the collinear approximation from now on. As mentioned previously, the STO-3G basis yields exceptionally stable H4+ intermediates. This is especially critical at low-energies, where not enough translational energy from the projectile might be transferred into an appropriate vibrational mode of the intermediate to allow for a bond breaking. In our experience, this excessive stability renders the STO-3G basis unusable for energies bellow 1 eV. The majority of the trajectories under these conditions yield the H4 ion as a stable species even after 40,000 a.u. of time. For that reason, we choose instead to make use of a Dunning double-zeta set[95] in this series. Collision Induced Dissociation In table 6-1 we show the total cross-section for the CID channel, i.e., the one leading to H3 + H products. We have found this to be the prevalent channel for impact parameters above 0.6 a.u., and up to 5 a.u. An example of one such trajectories is given in figure 6-2. The first box shows the nuclear positions in the plane of the collision (xy plane), as the reaction proceeds. The bond lengths for both projectile and target are shown in the second box. From that, we can easily see the H2 bond breaking, as its nuclei separate by up to 8 a.u. when the reaction reaches 8,000 a.u. of time. In the third box, we can see the changes in the Mulliken populations. The formation of the H3+ ion is clear above 5,000 a.u. of time. It is also apparent that the nucleus leaving as the H atom belonged originally to the target (H2), and was the one facing the projectile. This is not always the case, and in figure 6-3 we show a trajectory, with impact parameter of 3 a.u., for which the leaving nucleus originally belonged to the projectile. I _ H2(1) - H2+(2) - H 2(1) 1 H2(2) - 6-, 45 3 0 (-0 , -1 1 I I 1 1 1 1 7 1 7 1 " ï¿½ï¿½ ï¿½ 1" " ï¿½ 1 1 1 1 J -5 0 5 10 15 20 25 30 35 X Position (a.u.) I --H2+ --H2 31 2.1 0 1000 2000 3000 4000 5000 6000 7000 8000 Time (a.u.) a "0 1. 0.8 0,6. 'S 0.4" cc 0.2 C 0 a) 0 - Nudeus(l) - Nucleus(2) - Nudeus(3) - Nudleus(4) ] 1000 2000 3000 4000 5000 6000 7000 8000 Time (a.u.) Figure 6-2: Collinear series, 0.519 eV, DZ basis, impact parameter 5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 1 I - H2+(1) - H2+(2) - H2(1) - H2(2) 15 ( X Postion (au.) 25 30 35 40 I - IV - H2 - Nudous(l) - Nuclsus(2) - NucuW3) Nudeus(4) 4000 (au. Mime (a.u.) edoo 7000 8000 Figure 6-3: Collinear series, 0.519 eV, DZ basis, impact parameter 3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) The total-cross section for the radial approximation is given bya=2yrJJPb)bdb, where P(b) is the probability, and b the impact parameter. Therefore, the contribution of a single trajectory increases with the square of its impact parameter. 20 15 5 -5. ,12. 1 0.6 .4= 100 2000 3000 5 0 5 1b We can project the final wave function to obtain the probabilities P(b), as described earlier. We have found the main channel for each of trajectory to be totally predominant, with percentages higher than 99.9999%. In practice, P(b) is either one or zero for any choice of impact parameter. We also observed that there is no significant perturbation of the electronic populations for impact parameters above 6 a.u. and we take this value as the outer boundary for the reactive region. The dominance of the CID in the outer portion of the reaction space accounts for its high cross-section, and seems to be in accordance with both the IOSA results (shown in table 6-1 as QM), and with the ones from semiclassical theory (SCL). All the three methods, END, QM[30], and SCL[45], seem to be in good agreement with each other and with the experimental results (Exp.)[99]. It is particularly notable that END shows a percentual relative deviation of just about 7.7%, with respect to experimental number. Since the CID trajectories dominate the space close to the outer reactive shell, it is less likely that the choice of initial orientation for the reactants has any major impact on the final cross-sections. Table 6-1: CID cros-sections, in A2, for 0.519 eV (CM), DZ basis CI N MSCL EXP.24.0 19 30ï¿½1 26ï¿½1 Charge Transfer There are no available experimental results for the charge-transfer channel at this energy. In table 6-2 we list the IOSA, semiclassical[69], and END cross-sections. This channel was found to be the dominant one for impact parameters from slightly above 0.3 a.u. to just less than 0.6 a.u. A representation of one such trajectory is given in figure 6-4. We can see, in the first panel, the products as two distinct vibrationally excited diatomic fragments. The second panel, shows the reactants exchanging bond lengths, indicating an electron transfer at around 4,000 a.u. of time. Notice the higher frequency and lower amplitude of the H2 product. In the final panel, the Mulliken populations of the reactants are reversed, showing that the projectile(H2 ) acquired one electron from the target. - H2(1) - H21(2) - H2(1) - H2(2) I " " . I " . " " . I I I I . 0 5 10 15 20 X Position (a.u.) 25 3D 35 Time (a.u.) Figure 6-4: Collinear series, 0.519 eV, DZ basis, impact parameter 0.5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) The CT process is extremely localized in terms of initial conditions and limited to small impact parameters. As a consequence, its contribution for the total cross-section is minimal. An extrapolation for the CT channel carried by Phelps[94], and shown in Figures 5-1 and 5-2, seems to 2.5 2 1.5 tw 1 *05 -1.5 -2 I% -% ] indicate that the correct value for this cross-section should be close to 1.5 A2. This is far less than what is obtained by SCL and, specially, the IOSA methods. Still, it is higher than the END result. One possibility for this discrepancy might lie in the choice of initial conditions. By initially putting both reactants parallel to the reaction coordinate, we are effectively decreasing their level of interaction for small impact parameters. The combined half-bond lengths of the reactants would yield an extra 1.7 a.u. to the reaction shell, if aligned perpendicularly to the reaction coordinate. That would have a negligible effect on the CID channel, but could contain the major portion of the CT shell. Studying the effects of the orientation on the crosssections is beyond the scope of the collinear approximation and we leave it for a more detailed calculation, in the following sections. As an illustration of the whole CT region, we present the projectile Mulliken population as a function of the impact parameter and evolution time, in figure 6-5. The projectile bond length is shown in a similar graph in figure 6-6. We can see in the lower right quadrant the increase in population (figure 6-5), and the decrease in bond lenght(figure 6-6). A I P1 00 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 6-5: Collinear series, 0.519 eV, DZ basis; the projectile (H24) total Mulliken population, as a function of time and impact parameter (all data in a. u.) 000 0.5 1 1.5 2 2.5 3 35 rr A -" I Figure 6-6: Collinear series, 0.519 eV, DZ basis; the projectile (H2 ) bond length, as a function of time and impact parameter (all data in a. u.) For the purpose of comparison, figure 6-7 depicts a nonreactive trajectory for an impact parameter of 0.3 a.u. This is the lower boundary of the CT channel and we can see its similarities with the collision shown in figure 6-4. The main differences occur in the region where the H4+ is formed, 'I FTJ I Rn around 4,000 a.u. of time. The non-reactive channel displays a temporary electron transfer, clearly visible in the Mulliken populations. As with the CT collision, the products leave in an excited vibrational state. I - H 2+(1) - K2+(2) - H2() - H2(2) ........... (u.... . X Position (a~u.) 0 3..... VVVV Time (a.u.) Time (a.u.) Figure 6-7: Collinear series, 0.519 eV, DZ basis, impact parameter 0.3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 3, 2.5, 1.54 051 -.50 >- -14 - .. . . . . . M m . . I . II I I 1H2+ H2I Table 6-2: CT cros-sections, in A2, for 0.519 eV (CM), DZ basis CT END QM SCL 0.20 14 9ï¿½1 Full Reaction Space At higher energies, the excessive stability of H4 in the STO-3G basis is less of a problem. We learnt from practice that above 2 eV (CM) H4+ does not show up as a stable species, instead breaking up to yield a CT, a CID, or a NR product. On the other hand, from figure 5-2, we see that up to 4 eV (CM) it is feasible to carry on calculations without worrying about additional channels. At this collision energy, the evolutions are fast enough to allow us to be more comprehensive in the depth of treatment. With that in mind, we obtain the total cross-sections taking the full dimensionality of the reaction into account. We employ the STO-3G basis, and a choice of initial conditions that span all the five degrees of freedom albeit in a limited way. These limitations include restricting the tilt angles to just three values, 30, 60, and 90 degrees. We also limit the dihedral angle to 30, 60, and 120 degrees. The impact parameters range from 0 to 6 a.u., in intervals of 1 or 2 a.u. Further, we assume the reaction to be symmetric with respect to negative and positive impact parameters, decreasing the number of required trajectories to about onefourth. As it occurred with the 0.519 eV series, the probabilities for a given channel are sufficiently close to either 100% or zero for almost all the trajectories. Collision Induced Dissociation In figures 6-8 and 6-9, we show the reaction profiles for CID, in a variety of angles and impact parameters. Figure 6-8 groups the profiles for a dihedral angle of 30 degrees, and figure 6-9 for 60 degrees. It would seem likely by comparing those two groups of trajectories that the dihedral angle is playing little or no r6le in the outcome of the reaction. In figures 6-10 and 6-11, we present a more detailed view for the same channel but with a different choice of angles. The influence of the dihedral angle is now clear. These cases exemplify the difficulty of sampling the reaction space with a small number of trajectories. While figures 6-8 and 6-9 seem to point at some kind of symmetry that could be explored in order to reduce the number of required trajectories, it is not clear what the rule is. To determine it, a larger sampling of the space would be needed, which would defeat its purpose. -6 06030 Series 0 2- 4 6- - - - - - -6 -4 -2 0 2 4 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 Impact Parameter b (a.u.) A-. '2 2. Impact Parameter b (a.u.) 1 0.90.80.70.6050.40.30.20.1 0 .09030 Series -4. -- -2 0 2- 6- -6 -4 -2 0 2 4 6 Impact Parameter b (a.u.) Figure 6-8: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. J % to 1 0.90.80.70.60.50.40.30.20.1 0 -6 606060 Series -4 -- --2 - - - - - 62 61 1 i ll , -6 -4 -2 0 2 4 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 --I I I I I -4 -2 0 2 4 Impact Parameter b (a.u.) j 0 0 U. 1 0.90.80.70.60.50.40.30.20.1 0 I 1 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 -6 909060 Series -4- -2-- 2- 4,-- 6 -6 -4 -2 0 2 4 6 Impact Parameter b (a.u.) Figure 6-9: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. 1 0.90.80.70.60.50.40.30.20.1 0 -606060 Series ' -3.- - - - - - - - - -5. ---4-321 1234 6 -2.-----------1.------------2-----------6 -5 -4-3-2 -10 1 2 3 4 5 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 W' 1 I 1I 306060 Series .-5 1I1 11l1lEIElE l lE I El as E I- 0--- 13- -- - -- --6-5-4-3-2-1 0 1 2 3 4 5 6 Impact Parameter b (a.u.) I. CT 1 0.90.80.70.60.50.40.30.20.1 0 603060 Series -5 -4 1 " o 1 2 3 -6-54-3-2-1 01 3 4 5 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 2f llRfl Arian -6-------------2-- -------------1--------------3- - 4- 2--- 5--- -6-5-4-3-2-1 0 1 2 3 4 5 6 Impact Parameter b (a.u.) Figure 6-10: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. 1 0.90.80.70.60.50.40.30.20.1 0 6060120 Series -5. --4 -3 -2 0 14 2 5. -ii-- - 3 4 5 6 -6-5-4-3-2-1 0 1 23 45 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40-30,20.1 0 MOtMEMELLA I I 1 0.90.80.70.60.50.40.3020.1 0 6030120 Series -5, -- --4 -3, --- --2 2 3. 4 5 -6-5-4-3-2-1 0 1 2 3 4 5 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 I Impact Parameter b (a.u.) Impact Parameter b (a.u.) Figure 6-11: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees. a .1-, I It is also striking to note the difference in contrast to the reaction at low energies, where the CID channel was overwhelmingly prevalent. At 4 eV this is hardly the case, with only about 1% of the trajectories yielding H3' + H. Another important distinction is the location of these CID trajectories, all of them lying close to the small impact parameters region. This is in accordance to the what is expected, and shown in figure 5-2. It is encouraging to see that even with a minimal basis set the correct behavior of the system prevails. Table 6-3 compares the END total cross-section with both the trajectory-surface hopping(TSH)[42] method and the experimental value[58]. The discrepancy between the END result and the other two is within the precision of the grid of points used in the calculation. With so few reactive trajectories and a relatively coarse grid, it is reasonable to expect either missing a significant portion of the product, or "over-integrating" the outside boundary of the reaction shell. Table 6-3: CID cros-sections, in A2, for 4 eV (CM), STO-3G basis CID END TSH Exp. 0.9 1.1 0.6 Charge Transfer As we did for the CID channel, figures 6-12 and 6-13 display the reaction profiles for 60 and 90 degrees tilts, and 30 and 60 degrees dihedral. Once again, for this choice of angles, there is very little influence of the dihedral angle. With a finer grid and a different choice of angles, in figures 6-14 and 6-15, the reaction takes a different look. It shows a high asymmetry and a strong dependency on the choice of initial orientation. The most important feature observed is, perhaps, the dominance of the CT channel at both large and small impact parameters. A marked difference in behavior when compared to energies below 1.5 eV. This holds true for any choice of initial conditions tested and we take it as a sign that the integration will give a good representation for this basis set, despite the strong asymmetry and sensitivity to the initial conditions. 1 0.90.80.70.60.50.40.30.20.1 0 MMMt ! I 1I -6 606030 Series ~ 2 -6 4 -2 0 2 4 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.3020.1 0 | II :j -4. 0 -. Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 I Iz A-4. -2' 1 S2. ~4. 6. Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 ,M I I 1I -6 909030 Series -2- - m4 -2m 2 4 6 Impact Parameter b (a.u.) Figure 6-12: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. 1 0.90.80.70.60.50.40.30.20.1 0 -4 0'- I. 6 -4 -2 0 2 4 6 Impact Parameter b (a, u.) 1 0.90.80.70.60.50.40.30.20.1 0 -4 -2 0 2 4 Impact Parameter b (a.u.) Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 W17 I I I -4 -2 0 2 4 Impact Parameter b (a.u.) Figure 6-13: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. J1 0.90.80.70.60.50.40.30.20.1 0 - 606060 -6 1 0.90.80.70.60.50.40.30.20.1 0 1 V' I= I -6-5-4-3-2-1 0 1 23 45 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 I- IW- - 1-- Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.3020.1 0 603060 Series -5 -4-- -2 1 4 5 6 -6-5-4-3-2-1 0 1 2 3 4 5 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 -- - It 303060 Series -6 -5-i "2 1 2 3 4 5- 6- -6-5-4-3-2-1 0 1 2 3 45 6 Impact Parameter b (a.u.) Figure 6-14: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. 101 0.90.80.70.60.50.40.30.20.1 0 1, 'I AI 6060120 Series -5 -,4 -4i I I I1I I I 1 II .I C (. I-0 5 - -6-5-4-3-2-1 0 1 2 3 4 5 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 71 I= 3060120 Series ,I I I I I I I I I I I T -3 -2 6 -6--4--2-56 ImatPaaee-b(~. 6030120 Series -4' - -1 0 1 24 - - -- 3,5 - -- - -+ " 0 " r -I"1U I 0 -+ U 0 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 3030120 Series _,; I I 1 ,1 1 1 1 ! ! : :,-4. 2- Ne I I I I I I ! 1 I 5-I1 11111-11M .I I I III I I I I -6-5-4 -3-2-1 0 1 2 3 4 5 Impact Parameter b (a.u.) Figure 6-15: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees. I.......... I 411! III I! .1 S1 0.90.80.70.60.50.40.3020.1 0 i I I I 4- In table 6-4, we compare the total cross-sections as obtained from END, TSH, and experiment. Despite the small number of trajectories relative to the 40,000 in the TSH, END does remarkably well with a relative difference with respect to the experiment of just above 10%, compared to 103% for TSH. Table 6-4: CT cros-sections, in A2, for 4 eV (CM), STO-3G basis CT END TSH Exp. 0 10.5 19.3 9.5 Comparison with Molecular Dynamics Considering the integer probabilities we observed in the low-energy region, the question arises as to whether the electron-nuclei coupling is relevant for this system. And, if so, to what extent and in what form. To address this, we performed calculations using the molecular dynamics (MD) approach. In particular, we have conducted the calculations using the dynamic reaction coordinate (DRC) methodology[100-102], as implemented in the software package GAMESS[103]. The DRC is a classical trajectory method, based on quantum chemical potential energy surfaces, which can be either ab initio or semi-empirical. The initial conditions for a DRC trajectory are similar to the ones for END. One must provide the kinetic energy for the reactants, as well as the direction of the initial velocity. We have chosen those to match the values of selected END trajectories. The spatial distribution of the reactants and their separation were made to agree as well. There is a choice of wave-functions to describe the PES on which the DRC trajectory will evolve. We have selected the unrestricted Hartree-Fock (UHF) in the following examples. Final Probabilities Unlike END, the MD trajectories are constrained a single potential energy surface. As mentioned previously, almost all of the END evolutions for this system yielded 100% probability for the dominant channel. This is not an indication that these trajectories are limited to a single energy surface. However, it is more likely that we find the participation of higher lying surfaces for trajectories with split character (both channels present.) One can measure the probability for each channel by projecting the final wave-function, representing the products, on the known wave-function for the desired fragment. In table 6-5, we shown the probabilities for a series of calculations, with initial conditions of 60 degrees for both tilt angles, and the dihedral angle, collision energy 4 eV, and a STO-3G basis. This series provided some of the largest deviations from unity probability, of all the series studied. It is clear that even for this set of calculations the deviation is minor. The trajectory with the most pronounced split character is highlighted, and corresponds to the one with impact parameters 0 and 3 a. u., for the z and y directions, respectively. For this instance, the final probabilities are 8% for charge-transfer and 92% for the non-reactive channel. Table 6-5: Probabilities for the final electronic configurations with initial conditions 606060 (4 eV, STO-3G). The shaded box displays the highest mixed character and a comparison between END and dynamics on a surface is provided for this case. (a, b are the impact parameters; CT: charge transfer; NR: non-reactive) a b PC PM a b PC PM (a.u.) (a.u.) (a.u.) (a.u.) 0 0 1.00 0.00 2 5 0.00 1.00 0 1 0.00 1.00 3 0 0.00 1.00 0 2 0.00 1.00 3 1 0.00 1.00 3 0.08 0.92 3 2 0.00 1.00 0 4 0.00 1.00 3 3 0.00 1.00 0 5 0.00 1.00 3 5 0.00 1.00 0 6 0.00 1.00 4 0 0.98 0.02 1 1 0.01 0.99 4 1 0.00 1.00 1 2 0.00 1.00 4 2 0.00 1.00 1 3 0.00 1.00 4 3 0.00 1.00 1 4 0.00 1.00 4 4 0.00 1.00 1 5 0.00 1.00 5 0 0.00 1.00 2 0 0.00 1.00 5 1 0.00 1.00 2 1 0.00 1.00 5 2 0.00 1.00 2 2 0.00 1.00 5 3 0.00 1.00 2 3 0.00 1.00 6 0 0.00 1.00 2 4 0.00 1.00 In figure 6-16, we show the END evolution for the trajectory discussed. The first panel, represents the nuclei trails for the entire collision. The representation has been shifted to the center of mass to allow an easier comparison with the MD results. The projectile (H2+) is initially situated at the origin, while the target sits at 15 a. u. As the reaction proceeds, both come into contact. They acquire some vibrational energy, which can be seen in the undulating motion of the nuclei (notice that H2 mostly rotates around its center,) as the products leave the interaction region. The second panel shows the vibrational excitation more clearly, representing the bond lengths of the fragments. Most of the vibrational energy for this collision is transferred to the H2 ion. The bond lengths also suggest there was no charge transfer, since they oscillate around their initial values for both fragments. This is confirmed in the third panel, were we show the total Mulliken population for each center, as a function of the evolution time. I - H2+ - - H2 I ime (au.) I - (1 - [2] - 13] - [41 IZF 1500 2000 2500 ime (a.u.) 3 I....I.... 3000 3500 4000 Figure 6-16: 606060 END series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) Figure 6-17, displays the same trajectory above, this time obtained through the MD approach. A quick comparison shows the qualitative features of the END evolution are all present. The fragments follow roughly the same path, H2+ carries most of the final vibrational energy, and as before there is no charge transfer. There are however noticeable + 120.2 0 i 0 500 1000 changes in the amplitudes of vibration for both reactants. Such differences would not account for a change in the total cross sections, since the products would still be the same, but could certainly affect the distribution of final states. I - H2+ - H2 Time (a.u.) - NUCeU1) - Nudeue(2) _ Nucleus(3) - de 4) I 08 0.02 9. 01" 9 0 500 1000 1500 2D0 2500 3000 3500 4000 : Tlrne (a.u,) Figure 6-17: 606060 DRC series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) A more sensitive way to gauge the difference between END and the MD, for this system, is perhaps a comparison involving a collision for which the prevalent channel results in charge transfer. As in the previous section, we choose a collision energy of 4 eV but, instead of the STO-3G basis, we use the more sensitive pVDZ set. For this example, we set the reactants initially lying parallel to the reaction coordinate, reducing the problem by one impact parameter, and making it simpler to interpret. Figures 6-18 and 6-19 show an example of such trajectories for END and MD, respectively. Once again, we shifted the nuclei trajectories to the center of mass frame. It is apparent that both methods give comparable results. The main discrepancy being the vibrational states of the products, which show not only different amplitudes but also higher frequencies for END. I - H2(1) - H+(2) - H2(1) - H2(2) , 10 X Position (a.u.) Time (a.u.) -Nudeua() - NUcleus(2) -Nudeus3) - thice*s4) 0 500 1000 1500 2000 Time (a u.) 1~~~~ 250 I 2, 0 I W I W Figure 6-18: 000000 END series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) - 2 . -A . . . . . . . . , ' "] + 1., 048. 0.24 02 I H- H2+ - H2 I 15 2D 25 - H2'(1) - 12+(2) - H2(1) - -t,(2) j -8.1 -10 -5 0 5 10 X Position (a.u.) 15 2D 25 I - H2+ - H I 2.6 2A. S 22 2 '1.2 ,812 I . . . .5 0 I I * I I I " " I I ' I " I . . . . I . . . ' . . . ' . . . 0 500 1000 1500 2000 2500 3000 3500 4000 Tme (a.u.) I - Nucleue(I) - Nucdeue(2) - Nuleus(3) - Nucleus(4) I Time (au.) Figure 6-19: 000000 DRC series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) .2 1 o8 OO , '" """ "" I + a 0. 0. 90 8 CL 0 9 .Iw CHAPTER 7 CONCLUSIONS The H2 + H2 system has long been a favorite topic for both theoretical and empirical studies. In recent years, this interest has been revived, as the reaction not only displays both electron transfer and chemical exchange at energies low enough to be of chemical application, but also serves as a prototype for more complex ion-molecule interactions. We have analyzed the reaction by using the electron nuclear dynamics (END) theory using a single determinantal wave function for the electrons and a classical description for the nuclei. The END results compare favorably with the existing "state of the art" approaches in both semi-classical and quantum methods. We have found that the system energetics and dynamics can be qualitatively described by using small basis sets, even a minimal basis, provided the collision energy is large enough to allow for the dissociation of the H4+ intermediate. For energies below 2 eV, at least a DZ set is required. Also, the introduction of polarization functions, as in the case of the pVDZ set seems not to improve on the DZ results in terms of energetics. Our calculations at 0.519 eV indicate that the reactants can be considered as radially symmetrical without too much error. This approximation yields specially good results for the dominant channel (CID.) Since the chemical exchange region dominates the outer shell of the reaction space, it is unlikely the total cross-section for this channel could be significantly affected by the choice of initial orientations. We obtain a charge transfer cross-section considerably lower than other theoretical approaches for the collinear approximation. One possible explanation for that is: since the CT region predominates only at small impact parameters, its area is comparable to the size of the bond lengths. In other words, the collinear approximation might be showing its shortcomings, by reducing the effective reaction shell. We should mention that both the semiclassical approach and the IOSA method grossly overestimate the cross-section reported by Phelps, which in itself seems to be extrapolated higher than what it should be. The calculations exploring the full dimensionality of the system, at 4 eV, yield excellent results for both channels. They compare favorably to the semiclassical calculations, if the experimental result is taken as a reference. It is somewhat surprising, how much the reaction profiles change with the choice of initial conditions for both channels. That seems to contradict the results based on the collinear approximation and, perhaps, point to an overall averaging over the spatial orientations. We have found that the dynamics of the system is in almost its entirety deterministic concerning the products, for the range of energies studied, in the sense that only one channel is accessed for a given choice of initial conditions. This fact, and the range of low-energies we considered, appears to indicate the reaction proceeds mainly on a single PES. A comparison with the molecular dynamics approach for some of the most critical cases supports that and does not seem to be dependent on the size of the basis set, at least as far as the total cross-sections are involved. At higher energies, as well as for properties dependent on the vibrational state of the products, this might not be the case. LIST OF REFERENCES 1. Herzberg, G., Nature 163, 170 (1949). 2. Joseph, R.D., Wright, G.S., & Wade, R., Nature 311, 132 (1984). 3. Kolos, W. & Wolniewicz, L., The Journal of Chemical Physics 49(1), 404-410 (1968). 4. Dehmer, P.M. & Chupka, W.A., Journal of Physical Chemistry 99(6), 1686-1699 (1995). 5. Thomson, J.J., Philos. Mag. 24, 209 (1912). 6. Dempster, A.J., Philos. Mag. 31, 438 (1916). 7. Acufta, M.H., Neubauer, F.M., & Ness, N.F., J. Geophys. Res. 86, 8513 (1981). 8. Bohr, N., Medd. K. Vet.-Akad: S. Nobelinstitut 5, 1 (1919). 9. Hirschfelder, J.O., Journal of Chemical Physics 6(795)(1938). 10. Salmon, L. & Poshusta, R.D., J. Chem. Phys. 59, 3497 (1973). 11. Gaillard, M.J., Gemmell, D.S., Goldring, G., Levine, I., Pietsch, W.J., Poizat, J.C., Ratkowski, A.J., Remillieux, J., Vager, Z., & Zabransky, B.J., Phys. Rev. A 17, 1797 (1978). 12. Frye, D., Preiskorn, A., Lie, G.C., & Clementi, E., J. Chem. Phys. 92, 4948 (1990). 13. Geballe, T.R. & Oka, T., Astrophys. J. 342, 855 (1989). 14. Dalgarno, A., TERRESTRIAL AND EXTRATERRESTRIAL H3+, in ADVANCES IN ATOMIC, MOLECULAR, AND OPTICAL PHYSICS. 1994, Academic Press, Inc.: p. 57-68. 15. Saporoschenko, M., Phys. Rev. A 139, 349 (1965). 16. Lie, G.C. & Frye, D., J. Chem. Phys. 96(9), 6784-6790 (1992). 17. Frye, D., Preiskorn, A., Lie, G.C., & Clementi, E. Modern Techniques in Computational Chemistry. in ESCOM. 1990. Leiden: ESCOM. 18. Oka, T., Phys. Rev. Lett. 45, 531 (1980). 19. Carney, G.D. & Porter, R.N., J. Chem. Phys. 64, 4251 (1974). 20. Carney, G.D. & Porter, R.N., J. Chem. Phys. 65, 3547 (1976). 21. Meyer, E., Botschwina, P., & Burton, P., J. Chem. Phys. 84, 891 (1986). 22. Miller, S. & Tennyson, J., J. Mol. Spectrosc. 128, 530 (1988). 23. Miller, S. & Tennyson, J., J. Mol. Spectrosc. 136, 223 (1989). 24. Kirchner, N.J., Gilbert, J.R., & Bowers, M.T., Chemical Physics Letters 106(1,2), 7-12 (1984). 25. Borkman, R.F. & Cobb, M., J. Chem. Phys. 74, 2920 (1981). 26. Wright, L.R. & Borkman, R.F., The Journal of Chemical Physics 77(4), 1938-1941 (1982). 27. Poshusta, R.D. & Zetik, D.F., The Journal of Chemical Physics 58(1), 118-125 (1973). 28. Borkman, R.F. & Cobb, M., The Journal of Chemical Physics 74(5), 2920-2926 (1981). 29. Deumens, E., Diz, A., Taylor, H., & bhrn, Y., The Journal of Chemical Physics 96(9), 6820-6833 (1992). 30. Baer, M. & Ng, C.Y., The Journal of Chemical Physics 93(11), 7787-7799 (1990). 31. Tully, J.C., The Journal of Chemical Physics 59(9), 5122-5134 (1973). 32. Karplus, M., Porter, R.N., & Sharma, R.D., The Journal of Chemical Physics 59, 4393 (1965). 33. McLaughlin, D.R. & Thompson, D.L., Journal of Chemical Physics 59, 4393 (1973). 34. Sathyamurthy, N., Kellerhals, G.E., & Raff, L.M., The Journal of Chemical Physics 64(5), 2259-2261 (1976). 35. Sathyamurthy, N., Rangarajan, R., & Raff, L.M., Journal of Chemical Physics 64, 4606 (1976). 36. Gray, S.K. & Wright, J.S., Journal of Chemical Physics 66, 2867 (1977). 37. Miller, W.H. & George, T.F., The Journal of Chemical Physics 56(11), 5637-5652 (1972). 38. Tully, J.C. & Preston, R.K., Journal of Chemical Physics 55, 562 (1971). 39. Stine, J.R. & Muckerman, J.T., Chemical Physics Letters 44(1), 46-49 (1976). 40. Preston, R.K., Sloane, C., & Miller, W.H., Journal of Chemical Physics 60, 4961 (1974). 41. Lin, Y.-W., George, T.F., & Morokuma, K., Chemical Physics Letters , 49 (1975). 42. Eaker, C.W. & Schatz, G.C., The Journal of Chemical Physics 89(11), 6713-6718 (1988). 43. Schatz, G.C., Badenhoop, J.K., & Eaker, C.W., International Journal of Quantum Chemistry XXXI, 57-63 (1987). 44. Eaker, C.W. & Schatz, G.C., The Journal of Physical Chemistry 89, 2612 (1985). 45. Eaker, C.W. & Schatz, G.C., Chemical Physics Letters 127(4), 343-346 (1986). 46. Neuhauser, D. & Baer, M., Journal of Chemical Physics 93, 2862 (1989). 47. Neuhauser, D. & Baer, M., Journal of Chemical Physics 94, 185 (1990). 48. Kosloff, R., Journal of Chemical Physics 92, 2087 (1988). 49. Kosloff, R. & Kosloff, D., Journal of Chemical Physics 94, 185 (1990). 50. Kramer, P. & Saraceno, M., Geometry of the TimeDependent Variational Principle in Quantum Mechanics. 1981, New York: Springer. 51. Klauder, J.R. & Skagerstam, B.-S., COHERENT STATES APPLICATIONS IN PHYSICS AND MATHEjATICAL PHYSICS. ist ed. 1985, Singapore: World Scientific Publishing Co. Pte. Ltd. 911. 52. Perelomov, A.M., Generalized Coherent States and Their Applications. 18t ed. Texts and Monographs in Physics, eds. W. Beiglb6ck, et al. 1986, Berlin, Germany: Springer-Verlag. 320. 53. Deumens, E., 6hrn, Y., & Weiner, B., Journal of Mathematical Physics 32, 1166 (1991). 54. Delos, J.B., Review of Modern Physics 58, 287 (1981). 55. Thouless, D.J., Nuclear Physics 21, 225 (1960). 56. Herbst, E. & Klemperer, W., Astrophys. J. 183, 117 (1973). 57. Prasad, S.S. & Tan, A., Geophys. Res. Lett. 1, 337 (1979). 58. Guyon, P.M., Baer, T., Cole, S.K., & Govers, T.R., Chemical Physics 119, 145-158 (1988). 59. Krenos, J.R., LLehman, K.K., Tully, J.C., Hierl, P.M., & Smith, G.P., Chem. Phys. 16, 109 (1976). 60. Hierl, P.M. & Herman, Z., Chem. Phys. 50, 249 (1980). 61. Vance, D.W. & Bailey, T.L., The Journal of CHemical Physics 44(2), 486-493 (1966). 62. Futrell, J.H. & Abramson, F.P., Adv. Chem. Ser. 58, 123 (1966). 63. Douglass, C.H., McClure, D.J., & Gentry, W.R., The Journal of Chemical Physics 67(11), 4931-4940 (1977). 64. Lees, A.B. & Rol, P.K., The Journal of Chemical Physics 61(11), 4444-4449 (1974). |

Full Text |

PAGE 1 ANALYSIS OF THE H 2 + WITH H 2 REACTION USING ELECTRON NUCLEAR DYNAMICS By JUAN OREIRO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1997 PAGE 2 In memory of Maria Augusta. PAGE 3 ACKNOWLEDGMENTS I would like to thank my advisor. Professor Yngve Ohrn, for his always encouraging, warm, and inspiring attitude, not only in science but also in a personal level. I also wish to acknowledge the support, enthusiasm, and guidance of Dr. Erik Deumens who is always in a good mood despite my repeated attempts to bring chaos onto his neat computer network. Of course, this work would not be possible without the help and the many discussions with the rest of the OhrnDeumens group. I thank, in particular, Augie Diz, Ricardo Longo, and Hugh Taylor for their invaluable help and exchange of ideas. I couldn't forget to mention Monique Chacon-Taylor, who made things make more sense with her "Margarita nights," albeit I couldn't remember much in the morning. Finally, I would like to extend my gratitude to the other QTP members and personnel for creating a both stimulating and fun environment to do research. iii PAGE 4 TABLE OF CONTENTS page ACKNOWLEDGMENTS iii LIST OF TABLES vi LIST OF FIGURES viii ABSTRACT xiii 1 INTRODUCTION 1 2 THEORETICAL, TIME-DEPENDENT APPROACHES TO MOLECULAR DYNAMICS 11 Classical Approach 11 Semiclassical Approaches 13 Quantum Mechanical Approaches 15 3 THE ELECTRON-NUCLEAR DYNAMICS (END) FORMALISM 17 Method 17 Electronic and Nuclear Descriptions 18 The END Equations of Motion 19 4 BIBLIOGRAPHIC REVIEW 25 Collision Induced Dissociation 25 Charge Transfer 28 5 BASIS SETS AND ENERGETICS 31 Basis and H 4 + 31 STO-3G Basis 35 Dunning's DZ Basis 36 Dunning's pVDZ Basis 36 Energetics 37 6 CROSS-SECTIONS 44 The Collinear Approximation 45 IV PAGE 5 Collision Induced Dissociation 48 Charge Transfer 52 Full Reaction Space 58 Collision Induced Dissociation 59 Charge Transfer 65 Comparison with Molecular Dynamics 70 Final Probabilities 71 pVDZ basis 76 7 CONCLUSIONS 79 LIST OF REFERENCES 82 BIOGRAPHICAL SKETCH 90 v PAGE 6 LIST OF TABLES Table page Table 5-1: Contraction coefficients (c) and exponents (a) for the ST0-3G basis 36 Table 5-2: Contraction coefficients (c) and exponents (a) for Dunnin's pVDZ basis, with scaling factor 1.44 36 Table 5-3: Total energies, in atomic units, for the ground states 37 Table 5-4: Bond lengths for the equilibrium geometries in the ST0-3G basis, values in atomic units 38 Table 5-5: Bond lengths for the equilibrium geometries in the DZ basis, values in atomic units 38 Table 5-6: Bond lengths for the equilibrium geometries in the pVDZ basis, values in atomic units 38 Table 5-7: Energetics in kcal/mol. Energy required to break H 4 + into H 2 + + H 2 , H 3 + + H, and total energy released by CID, respectively 39 Table 6-1: CID cros-sections, in A 2 , for 0.519 eV (CM), DZ basis Table 6-2: CT cros-sections, in A 2 , for 0.519 eV (CM), DZ basis Table 6-3: CID cros-sections, in A 2 , for 4 eV (CM), ST0-3G basis 64 vi PAGE 7 Table 6-4: CT cros -sections , in A 2 , for 4 eV (CM), ST0-3G basis 70 Table 6-5: Probabilities for the final electronic configurations with initial conditions 606060 (4 eV, STO3G) . The shaded box displays the highest mixed character and a comparison between END and dynamics on a surface is provided for this case, (a, b are the impact parameters; CT: charge transfer; NR: non-reactive) 72 vxi PAGE 8 LIST OF FIGURES Figure page Figure 5-1: Available channels and cross-sections (A 2 ) for energies 0.1 to 10 5 eV (lab) 33 Figure 5-2: Cross-sections (A 2 ) for the range of energies studied, 1-8 eV (lab) 34 Figure 5-3: Energetics for the STO-3G, DZ, and pVDZ bases sets 40 Figure 5-4: UHF ground-state surfaces for the STO-3G, DZ, and pVDZ bases 41 Figure 5-5: Bond lengths diagram for the H 3 + and H 4 + species in the D 3h and C 2v geometries, respectively 42 Figure 6-1: Degrees of freedom for the H 2 + + H 2 reaction: a and (3 are the projectile and target tilt angles, respectively; y is the dihedral angle; a and b are the impact parameters in the x and y directions, respectively. 45 Figure 6-2: Collinear series, 0.519 eV, DZ basis, impact parameter 5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 49 Figure 6-3: Collinear series, 0.519 eV, DZ basis, impact parameter 3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom viii PAGE 9 panel: total Mulliken populations vs. time (all data in a. u. ) 50 Figure 6-4: Collinear series, 0.519 eV, DZ basis, impact parameter 0.5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u. ) 53 Figure 6-5: Collinear series, 0.519 eV, DZ basis; the projectile (H 2 + ) total Mulliken population, as a function of time and impact parameter (all data in a. u.) 55 Figure 6-6: Collinear series, 0.519 eV, DZ basis; the projectile (H 2 + ) bond length, as a function of time and impact parameter (all data in a. u.) 56 Figure 6-7: Collinear series, 0.519 eV, DZ basis, impact parameter 0.3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 57 Figure 6-8: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 60 Figure 6-9: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact ix PAGE 10 parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 61 Figure 6-10: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 62 Figure 6-11: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees 63 Figure 6-12: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u .)7 each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 66 Figure 6-13: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 67 x PAGE 11 Figure 6-14: Charge transfer profiles for collision energy 4 eV and ST0-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees 68 Figure 6-15: Charge transfer profiles for collision energy 4 eV and ST0-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees 69 Figure 6-16: 606060 END series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively ) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 74 Figure 6-17: 606060 DRC series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 75 Figure 6-18: 000000 END series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively); top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 77 xi PAGE 12 Figure 6-19: 000000 DRC series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively ) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) 78 Xll PAGE 13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy ANALYSIS OF THE H 2 + WITH H 2 REACTION USING ELECTRON NUCLEAR DYNAMICS By Juan Oreiro December, 1997 Chairman: N. Yngve Ohrn Major Department: Chemistry The END formalism addresses the solution of the timedependent Schrodinger equation, treating both electrons and nuclei simultaneously. It differs from the other traditional fully quantum mechanical time-dependent methods in that it does not require a potential energy surface (PES) to carry the nuclear motion. The interaction between electronic and nuclear motion is, therefore, obtained in a transparent way, not relying on PES gradients to obtain the coupling between electrons and nuclei. We analyze the H 2 + + H 2 reaction at energies below 4 eV using different approximations and basis sets. Other than the choice of initial conditions, form of wave function, and basis set, no constraints are imposed on the system evolution. The nuclei are treated in the classical limit, and the electronic part is described by a single determinantal , unrestricted wave function. We obtain properties, such as xiii PAGE 14 Mulliken populations, transition probabilites , and crosssections, from the resulting trajectories. These results are compared with other current theoretical approaches, and with experimental values. The relevance of the electron-nuclear coupling is estimated by comparing the END trajectories with molecular dynamics calculations for selected initial conditions in different basis sets. xiv PAGE 15 CHAPTER 1 INTRODUCTION Hydrogen is the most abundant element in the universe. Whenever matter begins to aggregate, the first molecules to form are H2 and carbon, nitrogen, and oxygen hydrides. From dense interstellar clouds to the atmosphere of the Jovian planets, these molecules play a central role, and ultimately determine, the chemistry of the system. On Earth, most of the hydrogen has escaped from the atmosphere. Its reactions do not affect us in a direct way. That did not prevent an early interest in hydrogen chemistry from developing. In particular, the reaction H2 + + H2 has attracted a crowd of experimental and theoretical researchers alike. This enthusiasm is justified. The reaction is rich in features, yet simple enough to allow an understanding of its driving principles. A great number of chemical processes is displayed by H2 + + H2. Long-range Coulomb forces, selective vibrational enhancement, multi-channel competition at low energies, are among the many features present. The ability to predict the outcome of a reaction, based on its initial conditions, constitutes the very core of chemistry. In this context, the meaning of "outcome" has undergone an update in the last few years with the advent of 1 PAGE 16 2 feasible time-dependent methods on the theoretical side and femtosecond spectroscopic techniques in the lab. Again, the H 2 + + H 2 processes provide an excellent opportunity to put those new tools to the test. At collision energies of about 8 eV (lab frame) or less, only two competing channels are present. The first one, commonly referred to as collision induced dissociation (CID), yields H 3 + + H as the product. The other channel involves an electron transfer from H 2 to H 2 + , and is usually called the charge transfer (CT) reaction. Above 8 eV, a multitude of channels is accessible but either CT or CID dominate in the range from 0 to 10 4 eV. Minor competition is offered by proton and electron production starting at around 8 and 100 eV, respectively, and the atomic Hydrogen Lyman and Balmer series at even higher energies. No experimental or theoretical data is available for energies above 10 4 eV and results in the lower end are scarce. Each of the species involved in this reaction, either as a reactant, product, or intermediate, has a history of its own. The H 2 weak quadrupole transitions, initially measured by Herzberg in the infrared range[l], have become an important tool in identifying astronomical objects[2]. The hydrogen molecule has been widely studied by both theoretical and experimental methods. Its ground state vibrational and rotational spectra were analyzed in detail[3] and data is PAGE 17 3 available for even highly excited Rydberg states[4]. The H 2 + ion is equally well known, being a common example in most introductory Quantum Chemistry texts. J. J. Thomson was the first to observe the existence of H 3 + in 1912 using a mass spectrometer coupled to his "positive rays" experiment [ 5 ] . Its preponderance over H 2 + was soon explained by the highly efficient and exothermic CID reaction, which converts the latter to the triatomic ion[6]. Since its discovery, H 3 + has puzzled chemists. For many years its structure was totally unknown, leading Eyring to refer to the problem as the "scandal of modern Chemistry [ 7 ]" . Bohr was the first to publish a theoretical paper about the ion[8]. In the mid and late 30' s Eyring, Hirschfelder [ 9 ] , and colaborators suggested that H 3 + exists in a triangular structure. However, it was not until 1973 that an accurate value for its ground state energy was made available [ 10 ] . It took a few more years and the, then new, technique of Coulomb explosion to establish its ground state as belonging to a C 2 V point group[ 11], A benchmark calculation employing an extensive Gaussian basis and a Hylleraas-type Cl wave function was able to improve on these results by a few tenths of milihartree[ 12 ] . In recent years, the interest in H 3 + has been revived. This triatomic has been found to be the main protonator agent in, and ultimately responsible for, all known interstellar chemistry [ 13 ] . A cascade of reactions in a highly ionizing PAGE 18 4 environment, where free electrons abound, form most of the common hydrides: H 3 + + X -Â» XH + + H 2 XH + + H 2 H 2 X + + H H 2 X + + H 2 H 3 X + + H H 3 X + + e~ H 2 X + H XH + H 2 where X could be, for instance, 0 or C, and e~ is a free electron. It was also suggested that the abundance of H 3 + could provide a measure of cosmic ray flux in an interstellar medium[14]. Due to its high rate of formation in hydrogen plasmas [15], its suppression is desirable in nuclear fusion[ 16] . It is surprising that for a molecule considered the simplest polyatomic and with only one barely stable electronic excited state[17], its spectrum was completely unknown until recently. Only in 1980 the first 15 spectral lines belonging to the v 2 fundamental band were identified [ 18 ] . There is a wide selection of high-quality PES's available for the H 3 system, starting with the one obtained by Carney and Porter in 1974 [19] and subsequently used to calculate the vibrational modes of H 3 + [20]. Their work has been further improved and it is now possible to find even full Cl surfaces, such as the Meyer-Botschwina-Burton (MBB) obtained with a (10s,4p,2d) Gaussian basis[21]. Very accurate fittings have been obtained for those PES's. There are no detectable discrepancies between experimental and theoretical PAGE 19 5 spectroscopic results based on those analytical expansions. In fact, these calculations have been used to predict rovibrational lines beyond the limits of the available experimental apparatus [22, 23]. These highly accurate potentials have made H 3 + a favorite target for studies. Another participating species, H 4 + , has been no less controversial. Deemed as non-stable for many years, it was only observed experimentally in 1984 [24], after detailed ab initio calculations suggested the structure was viable[2527], Now, it is known that H 4 + is not only stable, but its ground state PES is capable of holding several vibrational levels. Its lifetime is of the order of microseconds. A favorite topic of discussion seems to be whether H 4 + can be obtained from the H 2 + H 2 + reaction. According to current consensus, it can only be generated from breaking heavier clusters. This is based on experimental observations and no reasonable theoretical explanation seems to exist that can support this belief. If H 4 + were to be formed, due to its shallow ground state minimum, it would likely be found at low collision energies. No reliable calculations have been performed for collisions energies below 2 eV (lab frame). The difficulties associated with measuring a low-stability charged species might account for its absence in the list of products. This ion exists in a supermolecule-like ground state. The H 3 + ion, in an almost C 2 V geometry, is orbited by a H atom which lies far away aligned to one of the triangle PAGE 20 6 vertices. Conventional detection techniques, such as those used in time of flight (TOF) measurements, would allow smaller charged ions to interfere, a problem known to exist for H2 + and H3 + . That would be less likely to occur when H4 + is the result of splitting a heavier cluster, and might explain why it has been detected that way. The possibility of breakage into H3 + + H should not be discarded either, due to the relatively long times required for some detection methods . Unlike the situation for H3 + , there are not that many potential energy surfaces available for the H4 + system. The higher number of geometrical conformations to be spanned makes the work of mapping the reaction space both tedious and costly. Accurate fitting is also significantly more complex. The situation is even worse for excited electronic states, required to describe electron transfer properly. That fact might explain the dominant position of the CID channel as compared to CT, which has been borne out in a number of published studies. The inclusion of even simple excited states, such as the ones obtained via self-consistent field (SCF) or singles and doubles configuration interaction (SDCI,) seems to provide good results. It allows for the correct behavior in the CT channel, i.e., an increase of the cross section with the impact energy, and a decrease with the impact parameter [ 28 ]. PAGE 21 7 Although models developed within stationary-state frameworks have dominated Theoretical Chemistry in the past, frequently yielding excellent results, a full understanding of a chemical reaction encompasses its evolution in time. A dynamical process, with spectroscopic and structural properties evolving continuously, needs a time-dependent treatment. As successful as traditional approaches to molecular dynamics have been, these did not possess the tools to allow an insight into the reaction mechanisms in terms of basic principles. They have frequently resorted to extraneous assumptions in order to achieve a simple theoretical structure. The inception of time-dependent methods is not recent, but until the last decade they could not be classified as working tools . Both computational and theoretical problems contributed to this delay. For similar reasons. Chemistry has been depicted in terms of mostly adiabatic processes. Both the convenience and the success of the Born-Oppenheimer approximation have led to the belief that the coupling between electronic and nuclear motion was either small enough to be negligible or confined to short segments of a reaction. As a consequence, the Potential Energy Surface picture flourished and it is almost impossible not to rely on it when discussing chemical processes. That well established view has been challenged more intensely in the last few years . A particular example being the study of charge transfer. The refinement of PAGE 22 8 existing experimental techniques and the introduction of new ones made more apparent the deficiency of older models. The increased computer power allowed the development of new ones. The Electron-Nuclear Dynamics (END) [29] theoretical framework, goes beyond the above restrictions. It allows for a full and instantaneous interaction between electrons and nuclei, while developing approximations to the time-dependent Schrodinger equation. The END equations have been coded into the software package ENDyne and will be discussed in some detail later. The analysis of the resulting evolution is done by the program Evolve, which was revised for this study. In its present incarnation, END employs a single determinant built from spin-unrestricted orbitals with complex coefficients to describe the electrons. The nuclei are considered in the limit of narrow wave-packets. The system evolution proceeds without the need for a PES. Not only does that save computer time but also eliminates inaccuracies due to the fitting procedure. It is a one-step, general process. Traditional methods involve three steps: calculation of the PES, fitting of the surface, and carrying out the dynamics. (The second of them frequently is property or application oriented.) With END it is possible to calculate systems for which no PES exist and retrieve most properties, as a function of time, from the resulting wave function. Despite being investigated by an unusually high number of experimental techniques and theoretical approaches, the PAGE 23 9 reaction of H 2 + with H 2 is still a challenge. The fact that the system is the simplest one involving a charged and a neutral molecule could not explain alone the intense interest it has motivated. The existence of a high quality PES, and analytical fittings for those, has certainly helped to attract part of the crowd of dynamicists eager to test their approaches in a reaction with more than three nuclei. Reliable PES for systems with four or more atoms are still an exception at this time. Since most of the available approaches for molecular dynamics rely on one or more PES's, systems for which such surfaces exist are always a preferred target. The vast number of experimental results available for this reaction has also contributed to intensify the efforts of studying it using a variety of theoretical techniques. A fully quantum-mechanical, three-dimensional study was published only recently on the subject [30], The authors calculated both the reactive and the charge transfer cross sections for two values of collisional energy (0.252 eV and 0.519 eV). They used only 14 randomly selected initial geometries due to the high computational demand of the problem. Their results agree remarkably well with experimental measurements for the reactive cross section, but not for the charge transfer (CT) one. Several reasons come to mind to explain the failure in predicting the CT cross section. It might be an intrinsic limitation of their approach (Infinite Order Sudden Approximation IOSA), an inaccuracy due to the rather limited sampling of the reaction PAGE 24 10 space, a basis set artifact, a feature of the potential energy surface or fitting employed, or a combination of any of those and other factors. Also, even for the reactive cross section, there are many details that are not clear or could be explored more deeply. In order to do that, we decided to carry out calculations on the H2 + with H2 reaction using the END approach . PAGE 25 CHAPTER 2 THEORETICAL, TIME -DEPENDENT APPROACHES TO MOLECULAR DYNAMICS Cla ss i cal APPf Qflch The classical approach to molecular dynamics consists of two steps : the determination of the potential energy surface (PES) and the numerical tracing of the collision trajectories on that surface, following classical equations of motion. The latter is a rather trivial procedure, especially with the fast computers available now. The generation of the PES and its fitting is the crucial part of this method. Despite the fact that only one PES is required, the accurate representation of the potential topology is still a demanding task for most systems containing more than three atoms. In addition, processes involving multiple PES's are not adequately treated within this framework. There is a great number of processes that are not adiabatic, among them, charge-transfer, electronic energy transfer, collisional quenching, spin-forbidden reactions, and radiationless transitions. To accurately depict such phenomena one has to include not only all the PES's involved, but also the couplings between those surfaces [ 31 ] . Also, processes where quantum effects are noticeable cannot be correctly described. 11 PAGE 26 12 those include systems for which tunneling and zero-point energy are important. The usual way to obtain the potential is through ab initio calculations, although semiempirical and even empirical methods have also been employed[ 32 ] . The discrete set of energies must then be converted to a more suitable analytical form, in order to allow for the rapid evaluation of the derivatives of the potential energy. However, the large number of nuclear configurations involved in a molecular collision precludes the calculation of a complete surface, even for systems for which accurate energies can be obtained. The fitting must be based on a limited set of geometries, usually chosen by convenience or symmetry considerations. According to McLaughlin and Thompson[ 33 ] : "...The calculation of surfaces remains somewhat of a stumbling block to a completely ab initio approach to classical dynamical studies of chemical dynamics." The fitting of PES's has become a favorite topic of interest by itself and some examples can be found in the study on 2D cubic spline interpolation by Sathyamurthy, Kellerhals, and Raff [34], the scattering claculations for He + H 2 + by Sathyamurthy, Rangarajan, and Raff [35], or the H + H 2 reaction using classical trajectories by Gray and Wright[36], PAGE 27 13 genticiaggica l A ppcgaghaa This class of approaches comprises methods that employ simple approximate procedures to extend the validity of the classical models. They provide a better description of the molecular dynamics in the sense that one is no longer limited to purely adiabatic effects. In their most common implementation, two PES's are used. The coupling between the two surfaces is normally introduced in an ad hoc way. The evolution is carried out until a certain point of the trajectory is reached, where diabatic effects are believed important. At that point, a somewhat arbitrary procedure is used to trigger the surface hopping. How well the technique describes the system depends on how accurately the transition point is chosen and how good the guess for the coupling between the PES's is. The same limitations concerning the calculation and fitting of the potential surfaces apply, with the extra complication that now more than one surface is needed and they will not be necessarily of the same quality. One popular implementation of the semiclassical methodology is due to Miller and George [37] and it's aimed at describing electronic transitions in low energy atomic and molecular collisions, or involving quantized degrees of freedom. Their framework is based on the integration of Hamilton's equations using complex-valued time. Both coordinates and momenta are allowed to take on complex values and the transition happens when a trajectory crosses a branch PAGE 28 14 point (corresponding to an adiabatic surface intersection). The approach relies on the Feynman propagator and involves the construction of complex-valued time contours to find a hopping point. It is assumed that at low energy the dynamics is essentially a classical single-surface phenomena, with transitions being isolated events that happen in localized points of space and time. Another common technique was devised by Tully and Preston [38] and called the "surface hopping" method. It differs from the Miller and George approach in that the momentum is not required to be a smooth function of time. Rather, the momentum changes abruptly when an avoided crossing is reached, to allow for the propagation on the different PES. The main advantage in this case is that one does not have to deal with the complex-valued trajectories. Both techniques have been shown to be equivalent, provided there are only two surfaces involved in the crossing, the splitting of the surfaces may be expressed as a function of a single variable, and the trajectories have no turning points close to the intersection ("high energy limit") [39]. There's a great number of papers applying semiclassical methods, particularly surface hopping and its variations, to light systems. The semiclassical approach is sometimes called quasiclassical in these studies [40-45], PAGE 29 15 Quantum Mechanical Approaches This class of methods, unlike the previous ones, rely on fully quantum mechanical equations of motion. Coupling between different surfaces arise naturally and do not depend on predefined transition schemes. Presently, there are two main approaches in the literature for quantum mechanical evolutions on potential surfaces. The first one was proposed by M. Baer and coworkers and it is referred to simply as the quantum mechanical (QM) method. It's based on what the authors call "the negative imaginary decoupling potentials" (NIPD) to convert multi-arrangement channels into inelastic singlearrangement systems, in order to reduce the configuration space. The inelastic collision is then treated by means of the infinite order sudden approximation (IOSA). That method has been applied to a number of collinear and three-dimensional systems [46, 47]. Baer and Ng[30] have published the only paper where the full dimensionality of the H 2 + + H 2 reaction was taken into account in a quantum mechanical level, where they have used 28 randomly selected configurations (initial relative molecular orientations). The second approach in use is usually referred to as exact quantum mechanical methods [48, 49], They are "exact" in the sense that the error could, in principle, be made as small as required. This is in practice frustrated by computational constraints. In order to represent the system PAGE 30 16 wave function, a basis of a numerical grid and weights is required. The smaller the interval between the grid points, the better the accuracy achieved. The total size of the grid is also important, since boundary effects may compromise the results. It is frequently not possible to generate the grid with dimensions large enough to describe the reaction going on. Wave packets can approach the rim of the grid and impart non physical effects to the dynamics. To minimize that, one has to apply "cushion" potentials to simulate a natural exit for the reactants. Exact methods give excellent results when they can be applied. Currently, the larger systems studied with this technique consist of six atoms. PAGE 31 CHAPTER 3 THE ELECTRON-NUCLEAR DYNAMICS (END) FORMALISM Method The END formalism addresses the solution of the timedependent Schrodinger equation, treating both electrons and nuclei simultaneously. It differs from other traditional fully quantum mechanical time-dependent methods in that it does not require a PES to carry the nuclear motion. The interaction between electronic and nuclear motion is therefore obtained in a transparent way, not relying on PES gradients to obtain the coupling between electrons and nuclei. There are several advantages in not depending on PES. First, there is no need to calculate such surfaces per se. For mid-sized and large systems, obtaining a good PES is a major task. With several degrees of freedom, it can easily be the most time-consuming part of the calculation. Commonly, the groups who work on dynamics are not the ones who generate the PES. That creates two problems. First, the dynamicists don't have control over the PES quality. Second, they are limited to whatever data is available. The actual calculation of the surface is just one of the steps, however. The data is given in a set of points that must be interpolated and 17 PAGE 32 18 extrapolated to provide the necessary surface. The quality of the interpolation is critical to obtain a good representation of the dynamical event. So, another source of errors is introduced, over which one has little control. Besides, a chemical process frequently requires the use of several surfaces, corresponding to different electronic states. That is especially true at higher energies. Finding one highquality PES is difficult, finding several of them for the same system is almost a unique event. END is based on the Time-Dependent Variational principle (TDVP)[50], with a wave function parametrized in terms of Coherent States [51, 52], The main advantage of using Coherent States is that one gets simplified equations, non-redundant parameters in the definition of the electronic state, and it is easier to interpret the results. This kind of parametrization also guarantees that all wavefunctions of a particular family can be reached, should the dynamics so require . The END formalism has been implemented in the program ENDyne[29] for a number of wave functions, with some restrictions and approximations which will be discussed below. Electronic and Nuclear Descriptions As it is now, ENDyne is limited to treat nuclei as classical particles. This corresponds to the limit of infinitely narrow gaussian wave packets and defines their PAGE 33 19 time-dependent variables as their positions and momenta and retains all nonadiabatic coupling terms. The electronic part is expressed in terms of a truncated Gaussian basis and is restricted to single-determinantal wave functions. The formalism for a multi-determinantal electronic description is developed and under implementation [ 53 ]. That means that we can set the initial conditions for an evolution as a set of generalized coordinates for the nuclei and an electronic state, the latter expressed in terms of complex spin-orbital coefficients. The only additional data needed are the mass of the nuclei and their atomic number. There are three different choices for the electronic wave function presently available: a restricted singledeterminantal wavef unction, with either pure-spin orbitals or mixed-spin orbitals. A two-level model, for both classical quantum (Gaussian) nuclei is available. In this work we use the unrestricted description with the pure-spin option. The END Equations of Motion Here we don't give a detailed derivation of the END equations but just show the structure of them and where they come from. We are going to outline the expressions built in the ENDyne code. Instead of considering the most general case, we will limit ourselves to classical nuclei and assume that the electronic part of the wave function can be described by a complex, spin-unrestricted determinant. A PAGE 34 20 rigorous and more general account of the formalism behind the method can be found in the references listed previously. In our case, the atomic spin orbitals (AOs) are given by a set of Gaussians, *t(r) = exp[-fc( r-r 0 f (3.1) with an appropriate spin part, a or (3, attached to them. We can define a general, complex spin orbital for the electrons as â€¢(r)-X^Â«,(r)c (r (3.2) where the total number of AOs, K, consists half of a and half of (3 orbitals. For a Coherent State wave function | Â£), defined by a set of parameters {Â£} , we can write the first variation of the quantum mechanical action functional, A, as (3.3) and the solutions we look for are the ones for which 5A= 0, according to the TDVP (principle of least action). By introducing the functions Â£(Â£*Â£) =(^\h\^)/(^) and sfcyMck) for the total energy and overlap, respectively, and setting h = 1, after several algebraic steps we can rewrite the first variation of the action functional as 8A = j x. xÂ„ B 2 In 5 BE V \ BOtr Bt; p ) p . B 2 In 5 Â‘xjc X, dt =0 , (3.4) KJ PAGE 35 21 where the Â£ denotes differentiation of the complex parameter with respect to time. The solution of the TDVP expression will then be the dynamical equations given by ( iC 0 0 -*Cf V c A? / % (3.5) where we define the general element of the hermitian matrix C f _<9 2 ln s/ as c <* Ax, Â• There are a few ways one can choose the electronic basis dependence on the nuclear motion. Three of the most useful choices are called the nonfollowing basis (NFB), the static following basis (SFB) , and the dynamic following basis (DFB ) . In the limit of a complete basis set, the three types are equivalent. However, actual calculations are carried out in truncated basis sets, with a finite number of orbitals. The three bases will not necessarily yield the same results under those circumstances. A NFB is one for which the basis grid is fixed in time and space, the electronic orbitals do not move at all during the system evolution. The equations of motion assume their simplest form in this case. This is the ideal choice for carrying the algebraic manipulations and interpreting the equations, if one can assume completeness. A SFB is one for which the basis are centered on the average position of the nuclear wave packets, or the nuclear coordinates in the limit of infinitely narrow wave packets. PAGE 36 22 The AOs do not carry an explicit term containing the orbital velocity in these cases, they just follow the nuclei. A DFB is one for which the AOs are defined in such a way that they have an intrinsic velocity and their positions change with time. The velocities are embodied in factors with the properties of a plane wave and are usually referred to as electron translation factors (ETFs)[54]. Those terms are known to play an important role in high energy collisions. The general form of an AO (omitting spin), including the ETF would be Saying that the three bases above are equivalent in the completeness limit amounts to stating that their equations yield the same dynamics. Their metrics are the same. They can be therefore related by a transformation that does not change the metric. Such transformation is referred to as a symplectic one (or canonical, when the metric only contains zeroes and ones ) . At this point, ENDyne does not include ETFs in its equations, but a SFB representation is used. Ordinary chemistry takes place at relatively low energies. We therefore don't expect ETFs contributing significantly to the study proposed here. The omission of such factors have a negligible effect in the final results. (Their effect is PAGE 37 23 taken into account approximately by the basis chosen, since the MO coefficients are complex.) For NFB, the equations of motion with real nuclear coordinates, R, and momenta, P, become (z are the Thouless parameters for the electrons ),[ 55 ] f iC 0 0 0 0 -K? 0 0 0 0 0 -/ v 0 0 / 0 yielding the individual equations, Yz] % * Z % R % %-J dE Cz = -ir~7 , dz P*=V R t Â£/ and R * = m ( 3 . 7 ) ( 3 . 8 ) ( 3 . 9 ) ( 3 . 10 ) Assuming a complete basis, we can obtain the equivalent result for a SFB, after a symplectic transformation, as r iC 0 iC R o' (i'l (dt/\ 0 -d Â—iCf R 0 -J* z % ic; -iCl Crr -/ R % l 0 0 I Â°, <%>) where I is the identity matrix , z denotes a Thouless parameter in the SFB and the following definitions were made, C x d 2 Ji^zJR'r) dzdX. R'=R,PÂ’=P ( 3 . 12 ) and PAGE 38 24 ( 3 . 13 ) The individual equations in vector form for this case are therefore * dE Cz + C R R, Â— Â—i + 2ImC^z -C RtR( R, =Â— V r E ( 3 . 15 ) ( 3 . 14 ) ( 3 . 16 ) Those are the equations coded in the ENDyne program package, where a summation over the nuclear coordinates label, 1, is implied. C R and C,^ are simply the so called nonadiabatic correction terms, that in the END theory appear in the dynamical metric. PAGE 39 CHAPTER 4 BIBLIOGRAPHIC REVIEW Collision Induced Dissociation The reaction H 2 + + H 2 Â— Â» H 3 + + H has been extensively studied over the years. Its importance in interstellar [ 56 ] and atmospheric chemistry [ 57 ] , as well as its suitability for rigorous theoretical treatment, account for part of the interest. Collision induced dissociation (CID) between ion and molecules are part of a class of reactions generally referred to as collisional energy transfer. However, they differ from the most general case in two important aspects: first, the amount of energy transferred is large enough to comprise a significant part of the translational energy; second, the reaction mechanism may be complex, involving several electronic surfaces and competing reaction channels [58] . Experimentally the H 4 + system has been studied by crossed beam[59, 60], mass spectrometric[61, 62], merged beam[63, 64], photoionization [ 65 ] , ICR[66], and ion beam-gas cell [67] techniques. These experiments showed the reaction to be exoergic by 1.7 eV and, contrary to initial beliefs, proceeding through a direct mechanism[63, 59, 60] where no barrier is involved[63, 64] and confirming the instability of 25 PAGE 40 26 H 4 + with respect to the H 3 + + H products [ 68 ] . In addition, its total cross-section was found to be very large at thermal energies, and to fall off rapidly with increasing collision energy [63, 64]. In principle, there would be two ways for the collision induced dissociation to proceed. The H 3 + product could be formed by either proton or by a H atom transfer. The first mechanism would involve the breakage of the H 2 + ion (requiring approximately 2.6 eV, ) while the second would entail the rupture of the H 2 molecule (4.5 eV. ) Energetics alone would make it seem more likely to have a predominance of the proton transfer. Trajectory surface hopping calculations have shown this to be true[69], the H 2 + rupture not only dominates but it seems to be the exclusive mechanism for the reaction. Attempts to confirm this experimentally, by using isotopic substitution, have yielded a very different picture. It was found that both mechanisms occur, and have comparable crosssections, apparently due to intensive charge transfer between the reagents [63, 69, 60, 59]. Photoionization experiments have shown that at collision energies below 1 eV (center of mass,) the H 3 + formation can be inhibited by vibrationally exciting the H 2 + ion[70, 65]. At higher energies, the opposite seems to hold. The reason for this behavior has been explained in terms of the interaction of several potential energy surfaces [60, 59, 71], putting a rigorous treatment of the problem out of reach of most current theoretical methods . Such interaction appears to be PAGE 41 27 the underlying cause for the competing mechanisms found in experiments involving isotopic substitution. The (H 2 + D 2 ) + ground PES could be entered by two alternative paths, one belonging to H 2 + + D 2 , the other to D 2 + + H 2 . At infinite separation each channel would have its own PES, but as they approach, making the ion and the neutral bond equal, an avoided crossing puts them on the same ground PES. The barrier between the two PES decreases until it becomes smaller than the reagent zero-point energy, which happens at 8-10 Bohr. At this point charge hopping between the two channels sets in [72, 59, 60, 27], accounting for the formation of both D 2 H + and H 2 D + . The multiple surface nature of the H 4 + system also further complicates the comparison between experimental and theoretical results. Most experiments involving electron bombardment are likely to excite the H 2 + ion, populating higher vibrational states. Such experiments include crossed and merged beam studies [63, 64], and their results should be considered as involving the all the accessible potential energy surfaces [ 73 ] . The CID channel was also explored with the molecular beam technique, using the enhanced four-photon ionization to prepare the H 2 + ion in selected vibrational states distributions, for collision energies ranging from 1.5 to 5.3 eV (c m. ) [ 74 ] The resulting data included the state-selected integral and differential cross sections. On a different experiment, using the merging beams technique, Lees and PAGE 42 28 Rol[64], obtained the reaction cross sections for several isotopic substitutions, at energies below 10 eV. Charge .Transfer When a charge transfer (CT) process involving an ion and a molecule is exothermic, it is usually the dominant product channel for the reaction. Due to the large CT cross-sections, these reactions play a major role in both relaxation processes and reaction kinetics in ionized gases. The understanding of this class of processes is of special interest to gas discharges, lasers, flames, and controlled thermonuclear fusion. Although present in low concentration in flames, the ions are responsible for the thermal charge exchange reactions which create free radical chain carriers . For collisions between molecular ions and molecules, where spatial orientation provides an extra challenge in generating the interaction potentials, the availability of reliable potential energy surfaces has been one of the major stumbling blocks for theoretical treatments aiming at obtaining accurate cross-sections. For the H 4 + system, a number of PES has been made available, mostly for the ground state [75, 27, 71, 72, 28, 25]. On the experimental side, the challenges have not been smaller. Most of the studies done involve the measurement of cross sections as a function of the kinetic energy, and preparing the H 2 + ion by electron ionization [ 76, 77, 61, 7886], The H 2 + ions formed by electron ionization usually have a PAGE 43 29 broad distribution of vibrational states. Experimental results seem to indicate that a strong dependence of the cross section on the vibrational population of H 2 + is to be expected[76, 77]. This finding, with the difficulty in obtaining reliable PES, has hindered the efforts to obtain a sound comparison between experiment and theory. For that reason, in recent years, the preferred empirical approach to prepare the system in well characterized distribution of states has become photoionization [ 87 ] . It is worth pointing out that most of the available reliable experimental results concentrate on higher collision energies, in particular above 200 eV[88]. One explanation might be the suggested poor resolution of the time-of-f light (TOF) analysis at low energies. Due to inelastic charge transfer channels, the slow H 2 + products can be scattered by more than 10Â° away from the initial neutral H 2 beam direction [ 89 ] . The collection efficiency in the TOF analysis for product ions scattered at wide angles is poor. In addition, inelastic charge transfer can also lead the arrival of the H 2 + product ions to be spread in time. As a result, for the reaction H 2 + + H 2 , the TOF does not provide a reliable detection method at energies below 8 eV (center of mass), since it cannot distinguish between the H 3 + and the H 2 + ions. To circumvent those problems, Liao, Liao, and Ng[87], develop a ion-molecule reaction apparatus which combined the crossed ion-neutral beam method, high resolution photoionization mass spectroscopy, and charge transfer PAGE 44 30 detection. Using this technique, they obtained the total charge transfer cross-section, as a function of the vibrational state of H 2 + (for states up to $ 0 '= 4), at energies ranging from 8 to 200 eV (center of mass.) It was also found, that the low-lying rotational states of the H 2 + ion had no discernible impact in the outcome of the reaction, at least when the ion was in its ground vibrational state and the collision energy was below 4 eV (c.m.) A strong dependency on the vibrational state was observed, with its maximum around 2 eV. It however decreases as the collision energy is reduced, contradicting the theoretical results trajectory surface hopping [90] . In a later experiment, with improved resolution, Liao and Ng[91] extended those results finding a good agreement with the semiclassical energy conserving trajectory calculations [92] at 16 eV, but with significant deviation at lower energies. Other experiments include the use of pulsed synchroton radiaton and photoelectron-photoion coincidence to select the internal state of the H 2 + ion, enabling the measurement of relative cross sections for states with H 2 + (fl'= 0-10) at 16 eV. The coupling among several vibrational states of the reactants and products was also observed for the CT channel in the energy range 8-400 eV, using the crossed beam technique [93] . PAGE 45 CHAPTER 5 BASIS SETS AND ENERGETICS Bagjg and H 4 + As mentioned previously, only the two channels, CID and CT, are present at the energy range studied, in addition to the non-reactive one. A compilation of all the available theoretical and experimental data to date on the H 2 + + H 2 reaction was provided by Phelps [94], and is shown in figure 5-1. A detail of the region analyzed here is presented in figure 5-2. Before we present the results, we will discuss briefly the basis sets employed and our reasons to choose them. Due to the number of degrees of freedom involved, and the intensive computational nature of the END, our main concern is to find a basis that allows us to obtain a qualitative picture of the reaction, while minimizing the computer demand. Unlike the calculation of stationary states, in dynamics there is no common "feeling" of which basis will yield the desired results. That is especially true for END. In order to quantify how well each of the basis sets performs, we ran several calculations using an SCF wave function to minimize the energy of the projectile, the 31 PAGE 46 32 target, and some relevant intermediates. It's known that H 4 + is a stable molecule with a lifetime of about 10 iis, capable of supporting several vibrational levels [24]. It's also known that the SCF wave function can correctly describe that stability, yielding a ground state of C 2 V symmetry[26] . Although it's generaly accepted that H 4 + cannot be formed by collisions between smaller molecules, but instead should arise only from breaking heavier fragments, an accurate depiction of the reaction energetics demands that at least the relative energies between the species involved are reasonably well reproduced. We have chosen three different sets. Since our goal is to achieve a qualitative picture, we used the system energetics as a guide to initially gauge the performance of each. Here, we compare the results for the STO-3G basis, a double-zeta basis [95], and Dunning's pVDZ basis [96]. In a related calculation, the stability and geometry of H 4 + was analyzed by Wright and Borkman[25]. Since they do not make clear what scaling factor they have used in their double-zeta case, if any, we have chosen to carry out the tests with the factor suggested by Dunning (1.44). To study the collisional process we use a factor equal to 1.54 for both the STO-3G and the DZ basis. The total energy at key points, involving intermediates, entrance and exit channel products was used to compare their effectiveness. PAGE 47 Cress Section (X 2 ) 33 v H 3 + + H o CT a VE h + +H d Lyman o Balmer * e 10000 PAGE 48 34 E |ab ion (eV) = 2E C(ri 1 2 3 4 5 6 7 8 H 3 + + H CT VE H + +H -aLyman -Â©Balmer e' Figure 5-2: Cross-sections (A 2 ) for the range of energies studied, 1-8 eV (lab) PAGE 49 35 The three sets analyzed were the ST0-3G, and Dunning's DZ, and pVDZ. Both the number of trajectories to sample the reaction space and the time necessary to simulate a single collision allowing for product separation are high. Since the END effort scales as N 4 , where N is the number of basis functions, the size of the basis is a major concern. In practical terms, we have limited this study to the DZ basis and the pVDZ is used only as a reference for selected initial conditions. The PVDZ basis is our benchmark set, being very successful in previous similar calculations on smaller systems. It also allows us to see how the introduction of polarization functions affect the results. A complete calculation of the reaction using this basis is not feasible due to the computer load imposed, requiring nearly 40 times the effort of a minimal basis. The DZ basis falls in between the previous two and gives an estimate of the importance of diffuse functions. ST0-3G Basis This is a standard minimal basis, where a linear combination of three Gaussians has been fit to a single Slater-type function. The Is orbital on each center is then represented by one of the resulting functions. The optimal fit for a Hydrogen atom, according this procedure yields the basis listed in table 5-1. PAGE 50 36 Table 5-1: Contraction coefficients (c) and exponents (a) for the ST0-3G basis. c a Is orbital 0.44463 0.109818 0.535328 0.405771 0.154329 2.22766 The fitted exponents above were then modified to account for the tighter Is orbitals in the reactants as compared to the atom. The scaling factor, 1.5376, was obtained from the Slater exponent, 1.24. Dunning f s DZ Basis r 951 The DZ basis has been used successfully by Wright and Borkman to study the stability of H n + clusters, with n as high as 9 . The resulting energies and geometries were comparable to a larger [5s,2p] basis. Dunning ' s pVDZ Basis r 961 Table 5-2: Contraction coefficients (c) and exponents (a) for Dunnin's pVDZ basis, with scaling factor 1.44. c a Is orbital 0.01969 18.7344 0.13798 2.82528 0.47815 0.64022 2s orbital 0.50124 0.17568 2p orbital 1.00000 0.72700 PAGE 51 37 In order to ascertain the performance of the above sets, we assume the reaction proceeds through a simple path going from the reactants, to the intermediate, to the products. The differences in energy between these points is compared in the three bases. The pVDZ results are considered as the reference values . A geometry optimization, using a UHF wave function, was performed for each of the species involved. The calculations were not constrained to any particular point group, with the exception of the ones for the H 4 + ion. In this case, the optimization was limited to the C 2v symmetry to accelerate convergence. This should not impose a major constraint, since the H 4 + ground state is known to belong to the C 2v geometry. In all sets, H 3 + converged to the D 3h point group, as expected. Table 5-3 contains the resulting total energies. Table 5-3: Total energies, in atomic units, for the ground states STO-3G DZ pVDZ H -0.466582 -0.497637 -0.499278 h 2 + -0.582697 -0.586387 -0.600286 h 2 -1.117506 -1.126658 -1.128746 h 3 + -1.246860 -1.275824 -1.294005 h 4 + -1.743480 -1.777434 -1.799292 PAGE 52 38 Table 5-4: Bond lengths for the equilibrium geometries in the ST0-3G basis, values in atomic units. R R' R" h 2 + 2.00421 h 2 1.34592 H 3 + (D 3h ) 1.82718 1.82718 1.82718 H4 + (C 2v ) 1.89523 2.24201 1.57174 Table 5-5: Bond lengths for the equilibrium geometries in the DZ basis, values in atomic units. R R' R" h 2 + 1.98507 h 2 1.38142 H 3 + (D 3h ) 1.60570 1.60570 1.60570 H 4 + (C 2v ) 3.30625 1.64697 1.57097 Table 5-6: Bond lengths for the equilibrium geometries in the pVDZ basis, values in atomic units. R R' R" h 2 + 1.98053 h 2 1.41343 H 3 + (D 3h ) 1.68039 1.68039 1.68039 H/(C 2v ) 2.096235 1.74943 1.62159 PAGE 53 39 In tables 5-4, 5-5, and 5-6 we list the equilibrium geometries for the above species. They essentially agree with values listed in the literature and confirm H 4 + as a stable species. In table 5-8 we show the energies required to break H 4 + into H 2 + + H 2 (either CT or non-reactive), H 3 + + H, and the total amount of energy released by CID. The energy barriers for the three sets is shown in figure 5-3 as the minimal energy path for the exit channel of a bond-breaking reaction (H 2 + + H 2 -Â» H 4 + Â— Â» H 3 + + H) proceeding through a C 2 V H 4 + intermediate, and the reactants UHF surfaces in figure 5-4. Table 5-7: Energetics in kcal/mol. Energy required to break H 4 + into H 2 + + H 2 , H 3 + + H, and total energy released by CID, respectively. h; + h 2 H* + H Ae ST0-3G 27.16 18.85 8.31 DZ 40.40 2.49 37.91 pVDZ 42.84 2.36 40.48 PAGE 54 40 Species H 4 + PAGE 55 Energy (Hartree) 41 UHF PES for the Reactants in Their Ground State o H 2 (STO-3G) Â• H 2 + (STO-3G) a Hg (DZ) a H, + (DZ) H 2 (PVDZ) H 2 + (PVDZ) Internuclear Distance (Angstron) Figure 5-4: UHF ground-state surfaces for the STO-3G, DZ, and pVDZ bases PAGE 56 42 As a consequence, in the ST0-3G basis, it is difficult to break the H 4 + and, when it happens, the channels H 2 + + H 2 and H 3 + + H will be more balanced than in larger basis. Between CT and CID, the later has a better chance. The DZ and pVDZ basis seem to perform about the same, at least in terms of energy differences. It is interesting to note that the equilibrium geometry of H 4 + in the STO-3G basis resembles the two reactants aligned perpendicularly. For the other two basis, the geometry is similar to an H 3 + with an H atom aligned to one of the vertexes and orbiting far away from the H 3 + core. That weak bond would make it more feasible to attain the CID channel. PAGE 57 43 It is clear that increasing the basis set reduces the H 4 + stability. The percent relative difference between the DZ basis and the pVDZ one is just 5%. The DZ basis overshoots the exit barrier by about 5.5% of the pVDZ value, while the minimal basis requires over 7.5 times the energy when compared to the DZ number. Our results disagree with the findings of Wright and Borkman in that they report an increase in stability with the size of the basis. Their double-zeta result gives a bond-breaking energy of approximately 1.3 kcal/mol and a value of 3.2 kcal/mol for a [5s,2p] basis. That's too much of a discrepancy to be accounted for by differences in the scaling factor. According to our numbers, the DZ should perform comparably to a pVDZ basis, at least in energetic terms. The minimal basis, while displaying the correct qualitative features of the larger sets, might lead to extremely stable H 4 + intermediates, especially at lower values of collision energy . PAGE 58 CHAPTER 6 CROSS-SECTIONS The full treatment of the H 2 + + H 2 collision involves six internal nuclear degrees of freedom, corresponding to a choice of five parameters for the initial conditions. Consider the reactants in the lab frame, and take the z-axis as being the direction of the reaction. A single trajectory can be described by two tilt angles, a dihedral angle, and two impact parameters. The tilt angles refer to the inclination of each of the reactant bonds with respect to the z-axis. The dihedral angle is defined by planes formed by each of the reactant bonds and the z-axis. The two impact parameters are obtained from the distance between the center of the reactants and the z-axis. We show a diagram of the above in figure 6-1. In order to obtain the total cross-sections, it is necessary to run several trajectories. Each trajectory is defined by a choice of initial conditions, which include the assignment of a value for each of the degrees of freedom. With five degrees of freedom, the number of possibilities increases quickly, making it important to devise a way to sample the reaction space without incurring in a prohibitively large calculation. For example, a relatively coarse grid with 6 points per degree of freedom, would 44 PAGE 59 45 involve 7,776 trajectories. To decrease the mesh size to half of that would require 248,832 trajectories. Below, we explore a few possibilities to keep the calculation manageable. Figure 6-1 s Degrees of freedom for the H 2 + + H 2 reaction: a and (3 are the projectile and target tilt angles, respectively; y is the dihedral angle; a and b are the impact parameters in the x and y directions, respectively. The Collinear Approximation The H 2 molecule, and to some extent the H 2 + ion, possess an almost radial potential. This spherically symmetric potential suggests that the reactants could, in principle, be treated as atoms, which would reduce the problem to just one dimension. In fact, that idea has been used successfully in the past at energies of the order of 300 eV[97]. Although, it is more likely that the asymmetry of the potential plays a PAGE 60 46 role at low-energies, there is a possibility that this approximation could yield acceptable results concerning the total cross-sections. Since the calculation of the total cross-section involves averaging over all the degrees of freedom, any asymmetry could be averaged as well. To test that possibility, we performed calculations at very low energies, 0.519 eV, assuming both reactants as radially symmetric, and taking into account just one of the degrees of freedom. The reduction in the number of trajectories, by considering just the impact parameter, is especially welcome at very low energies. To avoid unnecessary oscillations in the electronic populations, due to the overlap of the reactants' orbitals, the projectile and the target must be kept at a reasonably large separation. The higher the frequency of the oscillations, the harder the integrator has to work to keep track of the evolution, leading to extended computation time. This problem is somewhat alleviated for small basis sets, where the overlap is only significant at close range. For the H 2 + + H 2 reaction, the total elimination of the interaction between the fragments is not possible, regardless of the amount of separation. The Coulomb interaction, due to the charged species, is always present, even at asymptotic distances. In this series of calculations, we use an initial separation of 20 a.u. between the centers of mass of the reactants. PAGE 61 47 Since we are assuming radial symmetry for both reactants, we must choose an initial spatial orientation for them which will be the same for every trajectory. It is important to notice that although we choose the initial orientation, the reactants are free to reorient as the reaction proceeds. This dynamically changing evolution is significantly less restrictive than some previous treatments of this system. A case in point is the full three-dimensional I OS A methodology [ 30 ] . Despite taking into account all the available degrees of freedom, the IOSA evolution makes no allowance for any re-orientation of the reaction participants [ 98 ] . For all practical purposes, the reactants are frozen in their initial orientations throughout the reaction. In our understanding, this is far more restrictive than a choice of fewer degrees of freedom, in particular for charge-transfer processes. In our calculations we choose to put both projectile and target initially to lie parallel to the reaction coordinate. For that reason, we name this series the collinear approximation from now on. As mentioned previously, the STO-3G basis yields exceptionally stable H 4 + intermediates. This is especially critical at low-energies, where not enough translational energy from the projectile might be transferred into an appropriate vibrational mode of the intermediate to allow for a bond breaking. In our experience, this excessive stability renders the STO-3G basis unusable for energies bellow 1 eV. PAGE 62 48 The majority of the trajectories under these conditions yield the H 4 + ion as a stable species even after 40,000 a.u. of time. For that reason, we choose instead to make use of a Dunning double-zeta set [95] in this series. Collision Induced Dissociation In table 6-1 we show the total cross-section for the CID channel, i.e., the one leading to H 3 + + H products. We have found this to be the prevalent channel for impact parameters above 0.6 a.u., and up to 5 a.u. An example of one such trajectories is given in figure 6-2. The first box shows the nuclear positions in the plane of the collision (xy plane), as the reaction proceeds. The bond lengths for both projectile and target are shown in the second box. From that, we can easily see the H 2 bond breaking, as its nuclei separate by up to 8 a.u. when the reaction reaches 8,000 a.u. of time. In the third box, we can see the changes in the Mulliken populations. The formation of the H 3 + ion is clear above 5,000 a.u. of time. It is also apparent that the nucleus leaving as the H atom belonged originally to the target (H 2 ), and was the one facing the projectile. This is not always the case, and in figure 6-3 we show a trajectory, with impact parameter of 3 a.u., for which the leaving nucleus originally belonged to the projectile. PAGE 63 49 Â— h 2 *0) Â— h 2 V) Â— h 2 0> Â— Â” 2 < 2 > X Position (a.u.) Â— H 2 + Figure 6-2: Collinear series, 0.519 eV, DZ basis, impact parameter 5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) PAGE 64 50 H 2 + (1) H/(2) ^(1) H 2 (2) Â— V Â— H 2 Figure 6-3: Collinear series, 0.519 eV, DZ basis, impact parameter 3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mul liken populations vs. time (all data in a. u.) The total-cross section for the radial approximation is given by <7 = 2 nj P{b)bdb , where P(b) is the probability, and b the impact parameter. Therefore, the contribution of a single trajectory increases with the square of its impact parameter. PAGE 65 51 We can project the final wave function to obtain the probabilities P(b), as described earlier. We have found the main channel for each of trajectory to be totally predominant, with percentages higher than 99.9999%. In practice, P(b) is either one or zero for any choice of impact parameter . We also observed that there is no significant perturbation of the electronic populations for impact parameters above 6 a.u. and we take this value as the outer boundary for the reactive region. The dominance of the CID in the outer portion of the reaction space accounts for its high cross-section, and seems to be in accordance with both the IOSA results (shown in table 6-1 as QM) , and with the ones from semiclass ical theory (SCL). All the three methods, END, QM[30], and SCL[45], seem to be in good agreement with each other and with the experimental results (Exp.) [99]. It is particularly notable that END shows a percentual relative deviation of just about 7.7%, with respect to experimental number. Since the CID trajectories dominate the space close to the outer reactive shell, it is less likely that the choice Â°f initial orientation for the reactants has any major impact on the final cross-sections. PAGE 66 52 Table 6-1: CID cros-sections, in A 2 , for 0.519 eV (CM), DZ basis CID END QM SCL Exp. C 24.0 19 30Â±1 26Â±1 Charge Transfer There are no available experimental results for the charge-transfer channel at this energy. In table 6-2 we list the IOSA, semiclassical [ 69 ] , and END cross-sections. This channel was found to be the dominant one for impact parameters from slightly above 0.3 a.u. to just less than 0.6 a.u. A representation of one such trajectory is given in figure 6-4. We can see, in the first panel, the products as two distinct vibrationally excited diatomic fragments. The second panel, shows the reactants exchanging bond lengths, indicating an electron transfer at around 4,000 a.u. of time. Notice the higher frequency and lower amplitude of the H 2 product. In the final panel, the Mulliken populations of the reactants are reversed, showing that the projectile(H 2 + ) acquired one electron from the target. PAGE 67 53 H 2 + (1) R, + <2) H 2 (1) H 2 (2) Â— V H, Figure 6-4: Collinear series, 0.519 eV, DZ basis, impact parameter 0.5 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) The CT process is extremely localized in terms of initial conditions and limited to small impact parameters. As a consequence, its contribution for the total cross-section is minimal. An extrapolation for the CT channel carried by Phelps [94], and shown in Figures 5-1 and 5-2, seems to PAGE 68 54 indicate that the correct value for this cross-section should be close to 1.5 A 2 . This is far less than what is obtained by SCL and, specially, the IOSA methods. Still, it is higher than the END result. One possibility for this discrepancy might lie in the choice of initial conditions. By initially putting both reactants parallel to the reaction coordinate, we are effectively decreasing their level of interaction for small impact parameters. The combined half -bond lengths of the reactants would yield an extra 1.7 a.u. to the reaction shell, if aligned perpendicularly to the reaction coordinate. That would have a negligible effect on the CID channel, but could contain the major portion of the CT shell. Studying the effects of the orientation on the crosssections is beyond the scope of the collinear approximation and we leave it for a more detailed calculation, in the following sections. As an illustration of the whole CT region, we present the projectile Mulliken population as a function of the impact parameter and evolution time, in figure 6-5. The projectile bond length is shown in a similar graph in figure 6-6. We can see in the lower right quadrant the increase in population (figure 6-5), and the decrease in bond lenght ( f igure 6 -6 ) . PAGE 69 55 Figure 6-5: Collinear series, 0.519 eV, DZ basis; the projectile (H 2 + ) total Mulliken population, as a function of time and impact parameter (all data in a. u.) PAGE 70 56 Figure 6-6: Collinear series, 0.519 eV, DZ basis; the projectile (H 2 + ) bond length, as a function of time and impact parameter (all data in a. u.) For the purpose of comparison, figure 6-7 depicts a nonreactive trajectory for an impact parameter of 0.3 a.u. This is the lower boundary of the CT channel and we can see its similarities with the collision shown in figure 6-4. The main differences occur in the region where the H/ is formed. PAGE 71 57 around 4,000 a.u. of time. The non-reactive channel displays a temporary electron transfer, clearly visible in the Mulliken populations. As with the CT collision, the products leave in an excited vibrational state. H 2 *(1) H, + (2) H a (1) H z (2) Â— V Figure 6-7: Collinear series, 0.519 eV, DZ basis, impact parameter 0.3 a. u.; top panel: nuclei trajectories on the xy plane; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) PAGE 72 58 Table 6-2: CT cros-sections, in A 2 , for 0.519 eV (CM), DZ basis CT END QM SCL a 0.20 14 9Â±1 Full Reaction Space At higher energies, the excessive stability of H 4 + in the STO-3G basis is less of a problem. We learnt from practice that above 2 eV (CM) H 4 + does not show up as a stable species, instead breaking up to yield a CT, a CID, or a NR product. On the other hand, from figure 5-2, we see that up to 4 eV (CM) it is feasible to carry on calculations without worrying about additional channels. At this collision energy, the evolutions are fast enough to allow us to be more comprehensive in the depth of treatment. With that in mind, we obtain the total cross-sections taking the full dimensionality of the reaction into account. We employ the STO-3G basis, and a choice of initial conditions that span all the five degrees of freedom albeit in a limited way. These limitations include restricting the tilt angles to just three values, 30, 60, and 90 degrees. We also limit the dihedral angle to 30, 60, and 120 degrees. The impact parameters range from 0 to 6 a.u., in intervals of 1 or 2 a.u. Further, we assume the reaction to be symmetric PAGE 73 59 with respect to negative and positive impact parameters, decreasing the number of required trajectories to about onefourth . As it occurred with the 0.519 eV series, the probabilities for a given channel are sufficiently close to either 100% or zero for almost all the trajectories. Collision Induced Dissociation In figures 6-8 and 6-9, we show the reaction profiles for CID, in a variety of angles and impact parameters. Figure 6-8 groups the profiles for a dihedral angle of 30 degrees, and figure 6-9 for 60 degrees. It would seem likely by comparing those two groups of trajectories that the dihedral angle is playing little or no role in the outcome of the reaction. In figures 6-10 and 6-11, we present a more detailed view for the same channel but with a different choice of angles. The influence of the dihedral angle is now clear. These cases exemplify the difficulty of sampling the reaction space with a small number of trajectories. While figures 6-8 and 6-9 seem to point at some kind of symmetry that could be explored in order to reduce the number of required trajectories, it is not clear what the rule is. To determine it, a larger sampling of the space would be needed, which would defeat its purpose. PAGE 74 60 1 0.90.80.70.60.50.40.30.20.1 0 Â— ~L I 1 I II I 1 0.90.80.70.60.50.40.30.20.1 0 Â— r-| "I I I m w PAGE 75 61 1 0.90.80.70.60.50.40.30.20.1 0 Â—Â—I 1 0.90.80.70.60.50.40.30.20.1 0 Â— Â— 1 1 1 1 1 1 u Impact Parameter b (a.u.) Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 n=n _L o 1 0.90.80.70.60.50.40.30.20.1 0 Figure 6-9: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. PAGE 76 62 1 0.90.80.70.60.50.40.30.20.1 0 i 1 -| |~1 n 1 0.90.80.70.60.50.40.30.20.1 0 HLLJ i iÂ»i 601 50 60 Series -6-5-4-3-2 -1 0 1 23 45 6 Impact Parameter b (a.u.) 3 -4 It E i 3 a Q 6 -5-4-32 -1 01 23456603C I6C Series Â•-1 -6 ) \ -> 5 -: > . ) Â• ) * \ Â£ 6 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 Â— I I I I I I 1 0.90.80.70.60.50.40.30.20.1 0 Impact Parameter b (a.u.) Impact Parameter b (a.u.) Figure 6-10: Bond breaking profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. PAGE 77 63 1 0.90.80.70.60.50.40.30.20.1 0 i Â— Â— iiimi i' r i i i 1 0.90.80.70.60.50.40.30.20.1 0 ^ 1 1 1 1 1 U Impact Parameter b (a.u.) 60301 20 Series -5d -4to i HI -1Â£ 0 !5 1 (8 1 a 2XJ -6 -5-4 -3 -2-1 0 1 2 3 4 5 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 o Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 H IT -6 -5. 3 Â“4(C n tof U! -1I O' 30301 20 Series 2J-J C,, j L Â±Â± ' -6 -5-4-3 -2-1 01 23456 Impact Parameter b (a.u.) Figure 6-11: Bond breaking profiles for collision energy 4 eV and ST0-3G basis; the shaded area marks the probability for the CID channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzzz, where xx is the projectile tilt angle, yy the target tilt angle, and zzz the dihedral angle, all in degrees. PAGE 78 64 It is also striking to note the difference in contrast to the reaction at low energies, where the CID channel was overwhelmingly prevalent. At 4 eV this is hardly the case, with only about 1% of the trajectories yielding H 3 + + H. Another important distinction is the location of these CID trajectories, all of them lying close to the small impact parameters region. This is in accordance to the what is expected, and shown in figure 5-2. It is encouraging to see that even with a minimal basis set the correct behavior of the system prevails. Table 6-3 compares the END total cross-section with both the trajectory-surface hopping (TSH) [42] method and the experimental value[58]. The discrepancy between the END result and the other two is within the precision of the grid of points used in the calculation. With so few reactive trajectories and a relatively coarse grid, it is reasonable to expect either missing a significant portion of the product, or "over-integrating" the outside boundary of the reaction shell. Table 6-3: basis CID cros -sections, in A 2 , for 4 eV (CM), ST0-3G CID END TSH Exp. o 0.9 1.1 0.6 PAGE 79 65 C ha r ge Tr ansfer As we did for the CID channel, figures 6-12 and 6-13 display the reaction profiles for 60 and 90 degrees tilts, and 30 and 60 degrees dihedral. Once again, for this choice of angles, there is very little influence of the dihedral angle. With a finer grid and a different choice of angles, in figures 6-14 and 6-15, the reaction takes a different look. It shows a high asymmetry and a strong dependency on the choice of initial orientation. The most important feature observed is, perhaps, the dominance of the CT channel at both large and small impact parameters. A marked difference in behavior when compared to energies below 1.5 eV. This holds true for any choice of initial conditions tested and we take it as a sign that the integration will give a good representation for this basis set, despite the strong asymmetry and sensitivity to the initial conditions. PAGE 80 66 1 0.90.80.70.60.50.40.30.20.1 0 Â— Â— Â— rn i i 1 0.90.80.70.60.50.40.30.20.1 0 ^ Â» 1 1 1 1 1 1 Impact Parameter b (a.u.) Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 mi i i i Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 Impact Parameter b (a.u.) Figure 6-12: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel , as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. PAGE 81 67 1 0 . 90 . 80 . 70 . 60 . 50 . 40 . 30 . 20.1 0 Â—"I -i-|"T 1 Impact Parameter b (a.u.) 1 0 . 90 . 80 . 70 . 60 . 50 . 40 . 30 . 20.1 0 an t i Impact Parameter b (a.u.) 1 0 . 90 . 80 . 70 . 60 . 50 . 40 . 30 . 20.1 0 htt r 'I I II 609060 Series 3 4 . si ctf 2 . i Â£ 2o ns f 4 -6 -4 -2 0 2 4 Impact Parameter b (a.u.) 1 0 . 90 . 80 . 70 . 60 . 50 . 40 . 30 . 20.1 0 l QflQnfin Qorioc -6 4 -2 0 2 4 6 Impact Parameter b (a.u.) Figure 6-13: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel, as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. PAGE 82 68 1 0.90.80.70.60.50.40.30.20.1 0 I -iÂ—i r ;i" n 6 . -5 9 -4S.-3o % Â£ 0RJ I Â£ 2Â« 3F 4 5' 6o CD 50 60 Series | m 1 i m L * Â— i -6-5-4-3-2 -1 0 1 23456 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 gHT a 603060 Series -5d -4Â£-3rt o Hi -1E 0123456iÂ£ -6 -5 -4 -3 -2 -1 0 1 23456 Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 Em =T"~r~T i ii Impact Parameter b (a.u.) 1 0.90.80.70.60.50.40.30.20.1 0 JUD J L 6 C _ -4 q . j L j -i \ E | i 1 1 * i 3 6 -6-5-4-3 -2-1 0 1 2 3 4 5 Impact Parameter b (a.u.) Figure 6-14: Charge transfer profiles for collision energy 4 eV and STO-3G basis; the shaded area marks the probability for the CT channel , as a function of both impact parameters (in a. u.); each series is labeled as xxyyzz, where xx is the projectile tilt angle, yy the target tilt angle, and zz the dihedral angle, all in degrees. PAGE 83 69 1 0.90.80.70.60.50.40.30.20.1 0 I 3 Â«j v -(G to I PAGE 84 70 In table 6-4, we compare the total cross-sections as obtained from END, TSH, and experiment. Despite the small number of trajectories relative to the 40,000 in the TSH, END does remarkably well with a relative difference with respect to the experiment of just above 10%, compared to 103% for TSH. Table 6-4: CT cros-sections , in A 2 , for 4 eV (CM), ST0-3G basis CT END TSH Exp. O 10.5 19.3 9.5 Comparison with Molecular Dynamics Considering the integer probabilities we observed in the low-energy region, the question arises as to whether the electron-nuclei coupling is relevant for this system. And, if so, to what extent and in what form. To address this, we performed calculations using the molecular dynamics (MD) approach . In particular, we have conducted the calculations using the dynamic reaction coordinate (DRC ) methodology [ 100-102 ] , as implemented in the software package GAMESS[103]. The DRC is a classical trajectory method, based on quantum chemical potential energy surfaces, which can be either ab initio or semi-empirical . The initial conditions for a DRC trajectory are similar to the ones for END. One must provide the kinetic energy for PAGE 85 71 the reactants, as well as the direction of the initial velocity. We have chosen those to match the values of selected END trajectories. The spatial distribution of the reactants and their separation were made to agree as well. There is a choice of wave-functions to describe the PES on which the DRC trajectory will evolve. We have selected the unrestricted Hartree-Fock (UHF) in the following examples. Final Probabilities Unlike END, the MD trajectories are constrained a single potential energy surface. As mentioned previously, almost all of the END evolutions for this system yielded 100% probability for the dominant channel. This is not an indication that these trajectories are limited to a single energy surface. However, it is more likely that we find the participation of higher lying surfaces for trajectories with split character (both channels present.) One can measure the probability for each channel by projecting the final wave-function, representing the products, on the known wave-function for the desired fragment. In table 6Â—5, we shown the probabilities for a series of calculations, with initial conditions of 60 degrees for both tilt angles, and the dihedral angle, collision energy 4 eV, and a STO-3G basis. This series provided some of the largest deviations from unity probability, of all the series studied. It is clear that even for this set of calculations the deviation is minor. The trajectory with the PAGE 86 72 most pronounced split character is highlighted, and corresponds to the one with impact parameters 0 and 3 a . u . , for the z and y directions, respectively. For this instance, the final probabilities are 8% for charge-transfer and 92% for the non-reactive channel. Table 6-5: Probabilities for the final electronic configurations with initial conditions 606060 (4 eV, STO-3G). The shaded box displays the highest mixed character and a comparison between END and dynamics on a surface is provided for this case, (a, b are the impact parameters; CT: charge transfer; NR: non-reactive) a (a.u. ) b (a.u. ) P CT p * NR a (a.u. ) b (a.u. ) Per P(JR 0 0 1.00 0.00 2 5 0.00 1.00 0 1 0.00 1.00 3 0 0.00 1.00 0 2 0.00 1.00 3 1 0.00 1.00 1 o 3 0.08 0.92 3 2 0.00 1.00 0 4 0.00 1.00 3 3 0.00 1.00 0 5 0.00 1.00 3 5 0.00 1.00 0 6 0.00 1.00 4 0 0.98 0.02 1 1 0.01 0.99 4 1 0.00 1.00 1 2 0.00 1.00 4 2 0.00 1.00 1 3 0.00 1.00 4 3 0.00 1.00 1 4 0.00 1.00 4 4 0.00 1.00 1 5 0.00 1.00 5 0 0.00 1.00 2 0 0.00 1.00 5 1 0.00 1.00 2 1 0.00 1.00 5 2 0.00 1.00 2 2 0.00 1.00 5 3 0.00 1.00 2 3 0.00 1.00 6 0 0.00 1.00 2 4 0.00 1.00 In figure 6-16, we show the END evolution for the trajectory discussed. The first panel, represents the nuclei trails for the entire collision. The representation has been shifted to the center of mass to allow an easier comparison with the MD results. The projectile (H 2 + ) is initially PAGE 87 73 situated at the origin, while the target sits at 15 a. u. As the reaction proceeds, both come into contact. They acquire some vibrational energy, which can be seen in the undulating motion of the nuclei (notice that H 2 mostly rotates around its center,) as the products leave the interaction region. The second panel shows the vibrational excitation more clearly, representing the bond lengths of the fragments. Most of the vibrational energy for this collision is transferred to the H 2 + ion. The bond lengths also suggest there was no charge transfer, since they oscillate around their initial values for both fragments. This is confirmed in the third panel, were we show the total Mulliken population for each center, as a function of the evolution time. PAGE 88 74 Â— h 2 + Â— h 2 Figure 6-16: 606060 END series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively ) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) Figure 6-17, displays the same trajectory above, this time obtained through the MD approach. A quick comparison shows the qualitative features of the END evolution are all present. The fragments follow roughly the same path, H 2 + carries most of the final vibrational energy, and as before there is no charge transfer. There are however noticeable PAGE 89 75 changes in the amplitudes of vibration for both reactants. Such differences would not account for a change in the total cross sections, since the products would still be the same, but could certainly affect the distribution of final states. Figure 6-17: 606060 DRC series, 4 eV, STO-3G basis, impact parameters 3 and 0 a. u. (y and z axes, repectively) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) PAGE 90 76 EYM-basis A more sensitive way to gauge the difference between END and the MD, for this system, is perhaps a comparison involving a collision for which the prevalent channel results in charge transfer. As in the previous section, we choose a collision energy of 4 eV but, instead of the ST0-3G basis, we use the more sensitive pVDZ set. For this example, we set the reactants initially lying parallel to the reaction coordinate, reducing the problem by one impact parameter, and making it simpler to interpret. Figures 6-18 and 6-19 show an example of such trajectories for END and MD, respectively. Once again, we shifted the nuclei trajectories to the center of mass frame. It is apparent that both methods give comparable results . The main discrepancy being the vibrational states of the products, which show not only different amplitudes but also higher frequencies for END. PAGE 91 77 H 2 + (1) Hg + (2) n^ 1 ) Â— H,(2) Â— h 2 + Â— h 2 Figure 6-18: 000000 END series, 4 eV, pVDZ basis, impact parameters 0 and 4 a . u . ( y and z axes , repectively ) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) PAGE 92 78 Â— H 2 + (1) Â— H, + (2) Â— H 2 (1) _ K,(2) Â— H / Figure 6-19: 000000 DRC series, 4 eV, pVDZ basis, impact parameters 0 and 4 a. u. (y and z axes, repectively) ; top panel: nuclei trajectories of the coordinates; center panel: bond lengths vs. time; bottom panel: total Mulliken populations vs. time (all data in a. u.) PAGE 93 CHAPTER 7 CONCLUSIONS The H 2 + + H 2 system has long been a favorite topic for both theoretical and empirical studies. In recent years, this interest has been revived, as the reaction not only displays both electron transfer and chemical exchange at energies low enough to be of chemical application, but also serves as a prototype for more complex ion-molecule interactions. We have analyzed the reaction by using the electron nuclear dynamics (END) theory using a single determinantal wave function for the electrons and a classical description for the nuclei. The END results compare favorably with the existing "state of the art" approaches in both semi-classical and quantum methods. We have found that the system energetics and dynamics can be qualitatively described by using small basis sets, even a minimal basis, provided the collision energy is large enough to allow for the dissociation of the H 4 + intermediate. For energies below 2 eV, at least a DZ set is required. Also, the introduction of polarization functions, as in the case of the pVDZ set seems not to improve on the DZ results in terms of energetics. Our calculations at 0.519 eV indicate that the reactants can be considered as radially symmetrical without too much 79 PAGE 94 80 error. This approximation yields specially good results for the dominant channel (CID.) Since the chemical exchange region dominates the outer shell of the reaction space, it is unlikely the total cross-section for this channel could be significantly affected by the choice of initial orientations. We obtain a charge transfer cross-section considerably lower than other theoretical approaches for the collinear approximation. One possible explanation for that is: since the CT region predominates only at small impact parameters, its area is comparable to the size of the bond lengths. In other words, the collinear approximation might be showing its shortcomings, by reducing the effective reaction shell. We should mention that both the semiclassical approach and the IOSA method grossly overestimate the cross-section reported by Phelps, which in itself seems to be extrapolated higher than what it should be. The calculations exploring the full dimensionality of the system, at 4 eV, yield excellent results for both channels. They compare favorably to the semiclassical calculations, if the experimental result is taken as a reference. It is somewhat surprising, how much the reaction profiles change with the choice of initial conditions for both channels. That seems to contradict the results based on the collinear approximation and, perhaps, point to an overall averaging over the spatial orientations. PAGE 95 81 We have found that the dynamics of the system is in almost its entirety deterministic concerning the products, for the range of energies studied, in the sense that only one channel is accessed for a given choice of initial conditions. This fact, and the range of low-energies we considered, appears to indicate the reaction proceeds mainly on a single PES. A comparison with the molecular dynamics approach for some of the most critical cases supports that and does not seem to be dependent on the size of the basis set, at least as far as the total cross-sections are involved. At higher energies, as well as for properties dependent on the vibrational state of the products, this might not be the case. PAGE 96 LIST OF REFERENCES 1. Herzberg, G., Nature 163, 170 (1949). 2. Joseph, R.D., Wright, G.S., & Wade, R., Nature 311, 132 (1984) . 3. Kolos, W. & Wolniewicz, L., The Journal of Chemical Physics 49(1), 404-410 (1968). 4. Dehmer, P.M. & Chupka, W.A. , Journal of Physical Chemistry 99(6), 1686-1699 (1995). 5. Thomson, J.J., Philos. Mag. 24, 209 (1912). 6. Dempster, A.J., Philos. Mag. 31, 438 (1916). 7. Acuna, M.H., Neubauer, F.M., & Ness, N.F., J. Geophys. Res. 86 , 8513 (1981). 8. Bohr, N., Medd. K. Vet.-Akad: S. Nobelinstitut 5, 1 (1919) . 9. Hirschfelder, J.O., Journal of Chemical Physics 6 ( 795 ) ( 1938 ) . 10. Salmon, L. & Poshusta, R.D., J. Chem. Phys. 59, 3497 (1973) . 11. Gaillard, M.J., Gemmell, D.S., Goldring, G. , Levine, I., Pietsch, W.J., Poizat, J.C., Ratkowski, A.J., Remillieux, J., Vager, Z., & Zabransky, B.J., Phys. Rev. A 17, 1797 (1978). 12. Frye, D., Preiskorn, A., Lie, G.C., & Clementi, E., J. Chem. Phys. 92, 4948 (1990). 82 PAGE 97 83 13. Geballe, T.R. & Oka, T., Astrophys. J. 3 42, 855 (1989). 14. Dalgarno, A., TERRESTRIAL AND EXTRATERRESTRIAL H 3 + , in ADVANCES IN ATOMIC , MOLECULAR, AND OPTICAL PHYSICS. 1994, Academic Press, Inc.: p. 57-68. 15. Saporoschenko, M., Phys. Rev. A 139, 349 (1965). 16. Lie, G.C. & Frye, D., J. Chem. Phys. 96(9), 6784-6790 (1992). 17. Frye, D., Preiskorn, A., Lie, G.C., & Clementi, E. Modern Techniques in Computational Chemistry, in ESCOM. 1990. Leiden: ESCOM. 18. Oka, T., Phys. Rev. Lett. 45, 531 (1980). 19. Carney, G.D. & Porter, R.N., J. Chem. Phys. 64, 4251 (1974). 20. Carney, G.D. & Porter, R.N., J. Chem. Phys. 65, 3547 (1976) . 21. Meyer, E., Botschwina, P., & Burton, P., J. Chem. Phys. 84, 891 (1986). 22. Miller, S. & Tennyson, J., J. Mol. Spectrosc. 128, 530 (1988) . 23. Miller, S. & Tennyson, J., J. Mol. Spectrosc. 136, 223 (1989). 24. Kirchner, N.J., Gilbert, J.R., & Bowers, M.T., Chemical Physics Letters 106(1,2), 7-12 (1984). 25. Borkman, R.F. & Cobb, M. , J. Chem. Phys. 74, 2920 (1981). 26. Wright, L.R. & Borkman, R.F., The Journal of Chemical Physics 77(4), 1938-1941 (1982). PAGE 98 84 27. Poshusta, R.D. & Zetik, D.F., The Journal of Chemical Physics 58(1), 118-125 (1973). 28. Borkman, R.F. & Cobb, M., The Journal of Chemical Physics 74(5), 2920-2926 (1981). 29. Deumens, E., Diz, A., Taylor, H., & Ohrn, Y., The Journal of Chemical Physics 96 ( 9 ), 6820-6833 (1992). 30. Baer, M. & Ng, C.Y., The Journal of Chemical Physics 93(11), 7787-7799 (1990). 31. Tully, J.C., The Journal of Chemical Physics 59 ( 9 ), 5122-5134 (1973). 32. Karplus, M. , Porter, R.N., & Sharma, R.D., The Journal of Chemical Physics 59 , 4393 (1965). 33. McLaughlin, D.R. & Thompson, D.L., Journal of Chemical Physics 59, 4393 (1973). 34. Sathyamurthy, N., Kellerhals, G.E., & Raff, L.M., The Journal of Chemical Physics 64(5), 2259-2261 (1976). 35. Sathyamurthy, N., Rangarajan, R. , & Raff, L.M., Journal of Chemical Physics 64 , 4606 (1976). 36. Gray, S.K. & Wright, J.S., Journal of Chemical Physics 66, 2867 (1977). 37. Miller, W.H. & George, T.F., The Journal of Chemical Physics 56(11), 5637-5652 (1972). 38. Tully, J.C. & Preston, R.K., Journal of Chemical Physics 55 , 562 (1971). 39. Stine, J.R. & Muckerman, J.T., Chemical Physics Letters 44(1), 46-49 (1976). PAGE 99 85 40. Preston, R.K., Sloane, C., & Miller, W.H., Journal of Chemical Physics 60 , 4961 (1974). 41. Lin, Y.-W., George, T.F., & Morokuma, K., Chemical Physics Letters , 49 (1975). 42. Eaker, C.W. & Schatz, G.C., The Journal of Chemical Physics 89(11), 6713-6718 (1988). 43. Schatz, G.C., Badenhoop, J.K., & Eaker, C.W., International Journal of Quantum Chemistry XXXI, 57-63 (1987) . 44. Eaker, C.W. & Schatz, G.C., The Journal of Physical Chemistry 89, 2612 (1985). 45. Eaker, C.W. & Schatz, G.C., Chemical Physics Letters 127(4), 343-346 (1986). 46. Neuhauser, D. & Baer, M. , Journal of Chemical Physics 93, 2862 (1989). 47. Neuhauser, D. & Baer, M. , Journal of Chemical Physics 94, 185 (1990). 48. Kosloff, R., Journal of Chemical Physics 92, 2087 (1988) . 49. Kosloff, R. & Kosloff, D., Journal of Chemical Physics 94, 185 (1990). 50. Kramer, P. & Saraceno, M. , Geometry of the TimeDependent Variational Principle in Quantum Mechanics . 1981, New York: Springer. 51. Klauder, J.R. & Skagerstam, B.-S., COHERENT STATES APPLICATIONS IN PHYSICS AND MATHEMATICAL PHYSICS. 1 st PAGE 100 86 ed. 1985, Singapore: World Scientific Publishing Co. Pte. Ltd. 911. 52. Perelomov, A.M., Generalized Coherent States and Their Applications. 1 st ed. Texts and Monographs in Physics, eds. W. Beiglbock, et al. 1986, Berlin, Germany: Springer-Verlag. 320. 53. Deumens, E., Ohrn, Y., & Weiner, B., Journal of Mathematical Physics 32, 1166 (1991). 54. Delos, J.B., Review of Modern Physics 58, 287 (1981). 55. Thouless, D.J., Nuclear Physics 21, 225 (1960). 56. Herbst, E. & Klemperer, W. , Astrophys. J. 183, 117 (1973) . 57. Prasad, S.S. & Tan, A., Geophys. Res. Lett. 1, 337 (1979) . 58. Guyon, P.M., Baer, T., Cole, S.K., & Govers, T.R., Chemical Physics 119, 145-158 (1988). 59. Krenos, J.R., LLehman, K.K., Tully, J.C., Hierl, P.M., & Smith, G.P., Chem. Phys. 16, 109 (1976). 60. Hierl, P.M. & Herman, Z., Chem. Phys. 50, 249 (1980). 61. Vance, D.W. & Bailey, T.L., The Journal of CHemical Physics 44(2), 486-493 (1966). 62. Futrell, J.H. & Abramson, F.P., Adv. Chem. Ser. 58, 123 (1966) . 63. Douglass, C.H., McClure, D.J., & Gentry, W.R., The Journal of Chemical Physics 67(11), 4931-4940 (1977). 64. Lees, A.B. & Rol, P.K., The Journal of Chemical Physics 61(11), 4444-4449 (1974). PAGE 101 87 65. Chupka, W.A. & Russell, M.E., The Journal of Chemical Physics 49(12), 5426-5437 (1968). 66. Huntress, W.T., Elleman, D.D., & Bowers, M.T., J. Chem. Phys. 55, 5413 (1971). 67. Doverspike, L.D. & Champion, R.L., J. Chem. Phys. 46 , 4718 (1967). 68. Anderson, S.L., Hirooka, T., Tiedemann, P.W., Mahan, B.H., & Lee, Y.T., J. Chem. Phys. 73, 4479 (1980). 69. Muckerman, J.T., Theoretical Chemistry 6 , 1 (1981). 70. Koyono, I. & Tanaka, K., J. Chem. Phys. 72, 4858 (1980). 71. Stine, J.R. & Muckerman, J.T., The Journal of Chemical Physics 68(1), 185-194 (1978). 72. Cobb, M. , Moran, T.F., Borkman, R.F., & Childs, R., Chem. Phys. Lett. 57, 326 (1978). 73. Anderson, S.L., Houle, F.A., Gerlich, D., & Lee, Y.T., The Journal of Chemical Physics 75(5), 2153-2162 (1981). 74. Pollard, J.E., Johnson, L.K., Lichtin, D.A., & Cohen, R.B., The Journal of Chemical Physics 95(7), 4877 (1991) . 75. Polak, R., Chem. Phys. 16 , 353 (1976). 76. Latimer, C.J., Browning, R., & Gilbody, H.B., J. Phys. B 2, 1055 (1969). 77. Hayden, H.C. & Amme, R.C., Phys. Rev. 172, 104 (1968). 78. Cramer, W.H., J. Chem. Phys. 35, 836 (1961). 79. Rothwell, H.L., Zyl, B.V., & Ame, R.C., J. Chem. Phys. 61 , 3851 (1974). PAGE 102 88 80. Wolf, F., Ann. Phys. 29, 33 (1937). 81. Simons, J.H., Fontana, C.M., E. E. Muschlitz, J., & Jackson, S.R., J. Chem. Phys. 11, 312 (1943). 82. Moran, T.F. & Roberts, J.R., J. Chem. Phys. 49, 3411 (1968). 83. Stedeford, J.B.H. & Hasted, J.B., Proc. R. Soc. London Ser. A 227, 466 (1955). 84. Cramer, W.H. & Marcus, A.B., J. Chem. Phys. 32, 186 (1960) . 85. Hollricher, 0., Z. Phys. 187, 41 (1965). 86. Koopman, D.W., Phys. Rev. 154, 79 (1967). 87. Liao, C.-L., Liao, C.-X., & Ng, C.Y., The Journal of Chemical Physics 81(2), 5672-5691 (1984). 88. Chupka, W.A., , in Ion-Molecule Reactions , J.L. Franklin, Editor. 1972, Plenum: New York. p. 33. 89. Campbell, F.M., Browning, R., & Latimer, C.J., J. Phys. B 14, 3493 (1981). 90. Muckerman, J.T., , in Theoretics Chemistry, D. Henderson, Editor. 1981, Academic: New York. p. 1. 91. Liao, C.L. & Ng, C.Y., J. Chem. Phys. 84(1), 197-200 (1986) . 92. Lee, C.-Y. & DePristo, A.E., The Journal of Chemical Physics 80(3), 1116-1125 (1984). 93. Lee, C.-Y., Depristo, A.E., Liao, C.-L., & Ng, C.Y., Chemical Physics Letters 116(6), 534 (1985). 94. Phelps, A. V. , Journal of Physical Chemistry Reference Data 19(3), 653-674 (1990). PAGE 103 89 95. Dunning, T.H., Jr., The Journal of Chemical Physics 53(7), 2823-2833 (1970). 96. Dunning, T.H., Jr., The Journal of Chemical Physics 90(2), 1007-1023 (1989). 97. Giese, C.F. & II, W.B.M., The Journal of Chemical Physics 39 ( 3 ), 739-748 (1963). 98. Baer, M. & Nakamura, H., The Journal of Physical Chemistry 91 , 5503-5509 (1987). 99. Shao, J.D. & Ng, C.Y., The Journal of Chemical Physics 84 ( 8 ), 4317-4325 (1986). 100. Stewart, J.J.P., Davis, L.P., & Burggraf, L.W., J. Comput. Chem. 8, 1117-1123 (1987). 101. Maluendes, S.A. & Dupuis, M. , J. Chem. Phys. 93, 59025911 (1990). 102. Taketsugu, T. & Gordon, M.S., J. Phys. Chem. 99 , 84628471 (1995). 103. Schmidt, M.W., Baldridge, K.K., Boatz, J.A., Elbert, S.T., Gordon, M.S., Jensen, J.J., Koseki, S., Matsunaga, N., Nguyen, K.A., Su, S., Windus, T.L., Dupuis, M. , & Motgomery, J.A., J. Comput. Chem. 14 , 1347-1363 (1993). PAGE 104 BIOGRAPHICAL SKETCH Juan Oreiro was born in Rio de Janeiro, Brazil. He enrolled in Chemical Engineering at Federal University of Rio de Janeiro (EQ-UFRJ. ) As an undergraduate he received scholarships for scientific research in the fields of organic synthesis and analysis, natural products, and inorganic synthesis . After graduation he enrolled in a Master's degree in Chemistry, at the Institute of Chemistry-UFRJ. During that period he received a scholarship from the Braziliam government and started doing research in Quantum Chemistry. When visiting Gainesville, during the 1988 Winter Institute on Quantum Chemistry he met Professor Ohrn. Later, he came to the University of Florida and joined the OhrnDeumens group. 90 PAGE 105 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, a dissertation for the degree of Doctor /ot Philosophy. N. Yngve Professo as lrn. Chairman of Chemistry I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctoj^?of Philosophy. ett Rodney-Barytlet Graduate Research Professor of Chemistry I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctorvof^Philosophy . JohjF' 5abin Professor of Physics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, a dissertation for the degree of Dyoctor of Philosophy. as William Weltner, Professor of Chemistry I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy Professor of Mathematics PAGE 106 This dissertation was submitted to the Graduate Faculty of the Department of Chemistry in the College of Liberal Arts and Sciences and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December, 1997 Dean, Graduate School |