# Full text of "Exploring a new time-dependent method for molecular quantum dynamics"

## See other formats

EXPLORING A NEW TIME-DEPENDENT METHOD FOR MOLECULAR QUANTUM DYNAMICS By RICARDO LUIZ LONGO 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 To Dalva, Adhemar, Renata, Milton, and Ivani ACKNOWLEDGMENTS I wish to acknowledge my advisor Professor Yngve Ohm not only for his counsel, but mainly for what he stands for. Also, I would like to thank Dr. Erik Deumens and their group, namely, Augie Diz, Juan Oreiro, and Hugh Taylor for their support, encouragement, and fruitful discussions. I thank Michael Zemer and his group for their friendship in the beginning of my time in Lei 212 and the QTP personnel. Finally, I could not have so much fun without some of my good friends, like Marshall and Genny, Rajiv and Pria, Sylvio, Osiel, Caco, Luciana, Andres, Victor, Monique, Bill, Dave, Mark, Charlie. lU TABLE OF CONTENTS ACKNOWLEDGMENTS . . LIST OF TABLES LIST OF FIGURES ABSTRACT CHAPTERS 1 INTRODUCTION . 2 TIME-DEPENDENT METHODS 5 PES Based Time-Dependent Methods 6 PES Free Time-Dependent Methods 9 Close-coupling Methods 10 Molecular Dynamics Methods 10 Time- Dependent Hartree-Fock Methods 11 3 THE THEORETICAL APPROACH 14 The END Formalism 15 Coherent States 15 The Time-Dependent Variational Principle 20 The END-SD-FGWP Method 22 The Equations of Motion 22 The NDDO Approximation 31 The END-SD(NDDO)-FGWP Method 34 AMI Model 37 Implementation 40 4 H+ + H, He, AND H 2 CHARGE TRANSFER COLLISION . . . The H-^ + H Collision Total Cross Sections Differential Cross Sections Trajectory Analysis Transfer Probability Reduced Differential Cross Section The H^ + He Collision Trajectory Analysis The Deflection Function and the Interaction Potential Classical Differential Cross Section Semiclassical Elastic Differential Cross Section Experimental and Theoretical Elastic Differential Cross Sections The H^ + H 2 Collision Trajectory Analysis Vibrational Analysis Transfer Probability Deflection Function Differential Cross Sections Summary and Conclusions . 44 . 46 . 46 . 50 . 52 . 54 . 59 . 63 . 64 . 67 . 73 . 78 . 85 . 87 . 88 . 92 . 95 . 96 . 97 . 99 5 INTRAMOLECULAR CHARGE TRANSFER DYNAMICS 101 Electron Transfer Formalisms 101 Structure, Energetics, and Electron Density of LiCNLi 104 Structure, Energetics, and Electronic Population of LiHLi 110 Dynamics of Electron Transfer in LiHLi 112 Summary and Conclusion 117 6 CONCLUSIONS 119 REFERENCES 122 BIOGRAPHICAL SKETCH . 129 LIST OF TABLES Table 4-1: Contraction coefficients c and exponents a for the basis sets used in this Table 4-2: Ionization potentials, electron affinities and polarizabilities of H, He and Table 4-3: Total cross section for + H collision at 500 eV using pVDZ basis set. . 48 Table 4-4: Total transfer cross sections for colliding with a H (x 10"^^ cm^). . 49 Table 4-5: Experimental and theoretical rainbow angles for H'*' + He system. ... 86 Table 4-6: Vibrational frequencies (cm“^) as function of the impact parameter (a.u.) and molecular orientation (a°, 0°) 95 Table 5-7: Minima for the linear structure of LiCNLi at various levels of theory. Bond distance in pm 105 Table 5-8: Relative energies of minima I and II and of transition state I-II (C*). Energies in kcal/mol 106 Table 5-9: Mulliken population (a.u.) and dipole moments (D) of structures I and II, transition state I-II and global minimum at UHF/3-21G level. . . 107 Table 5-10: Mulliken population, structure, and relative energy of the LiHLi molecule at SCFAJHF level 113 VI LIST OF FIGURES Figure 3-1: Figure 4-2: Figure 4-3: Figure 4-4: Figure 4-5: Figure 4-6: Figure 4-7: Figure 4-8: Functional diagram of the ENDyne program 42 Weighted transition probabilities for total electron transfer at 500 eV as a function of impact parameter from ENDyne. All data in atomic units. . 49 Scattering angle as a function of impact parameter using ENDyne, bare nuclei, and hydrogen atom potential for the + H system. Energy: 500 eV. Basis set: pVDZ. Scattering angle in degrees and impact parameter in a.u.. Full line: ENDyne, dotted line: bare nuclei, and dashed line: hydrogen atom potential 53 Scattering angle as a function of impact parameter using ENDyne, bare nuclei, and hydrogen atom potential for H'^ + H system. Small impact parameter region. Energy: 500 eV. Basis set: pVDZ. Scattering angle in degrees and impact parameter in a.u.. Full line: ENDyne, dotted line: bare nuclei, and dashed line: hydrogen atom potential 54 Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 250 eV. Basis set: pVDZ. Experimental angular resolution ±0.6°, ’+’ ENDyne, from Ref. , and ’x’ from Ref. ... 56 Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 410 eV. Basis set: pVDZ. Experimental angular resolution ±0-2° — ±0.6°, ’+’ ENDyne, from Ref. , and ’x’ from Ref 57 Transfer probability versus the scattering angle (in degrees) for the + H system. Energy: 500 eV. Basis set: pVDZ. ’+’ ENDyne 57 Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 700 eV. Basis set: pVDZ. Experimental angular resolution ±0.02° — ±0.6°, ’+’ ENDyne and from Ref • • Vll 58 Figure 4-9: Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 1(KX) eV. Basis set: pVDZ. Experimental angular resolution ±0.07° — ±0.6°, ’-i-’ ENDyne and from Ref. 58 Figure 4-10: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the -f- H system. Energy: 250 eV. Basis set: pVDZ, ENDyne transfer, ’+’ ENDyne elastic, ’ X ’ transfer and ’o’ elastic from Ref. 60 Figure 4-11: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the H'*' -i- H system. Energy: 410 eV. Basis set: pVDZ, ’*’ ENDyne transfer, ’-i-’ ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. 61 Figure 4-12: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the -(- H system. Energy: 5(X) eV. Basis set: pVDZ, ’*’ ENDyne transfer, ENDyne elastic, ’ X ’ transfer and ’o’ elastic from Ref. 61 Figure 4-13: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the -(- H system. Energy: 700 eV. Basis set: pVDZ, ’*’ ENDyne transfer, ’-i-’ ENDyne elastic, ’ X ’ transfer and ’o’ elastic from Ref. 62 Figure 4-14: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the + H system. Energy: 1(X)0 eV. Basis set: pVDZ. ’*’ ENDyne transfer, ’+’ ENDyne elastic, ’ X ’ transfer and ’o’ elastic from Ref. 62 Figure 4-15: Dynamical trajectories of being scattered by He. ENDyne results at energy of 50.0 eV with a pVDZ basis set 65 Figure 4-16: Dynamical trajectories for scattering angle of 2.65°. ENDyne results for H"'’ + He system at 50.0 eV with a pVDZ basis set. Impact parameters: solid line ( — ) = 1.12 a.u., dashed line ( — ) = 1.60 a.u., dotted line (• • •) = 2.10 a.u 66 VUJ Figure 4-17: Motion of the atomic target (He) during a collision with H"*" at 50.0 eV as computed by ENDyne with a pVDZ basis set. The initial coordinate of the He atom is (15.0, 0.0). Impact parameters: solid line ( — ^) = 1.12 a.u., dashed line ( — ) = 1.60 a.u., dotted line (• • •) = 2.10 a.u 67 Figure 4-18: Deflection function for H"^ + He. ENDyne results at 50.0 eV with a pVDZ basis set 68 Figure 4-19: Potential energy curve obtained from the inversion of the deflection function generated by ENDyne for H"^ + He at 50.0 eV with a pVDZ basis set 70 Figure 4-20: Potential energy curve obtained from the direct time evolution. ENDyne results for H"^ + He at 50.0 eV with a pVDZ basis set 71 Figure 4-21: Potential energy curves for H"^ + He obtained from a CIS calculation with a pVDZ basis set 72 Figure 4-22: Deflection function for H'^ + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: I) diamonds ’o’, II) triangles ’a’, and. III) open circles ’o’ 74 Figure 4-23: Elastic differential cross section for H"^ + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: I) diamonds ’o’, II) triangles ’a’, and, in) open circles ’o’ 75 Figure 4-24: Charge transfer differential cross section for H^ + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: I) diamonds ’o’, II) triangles ’a’, and, IE) open circles ’o’ 76 Figure 4-25: Elastic and charge transfer differential cross sections for H'*' + He. ENDyne results at 500 eV with a pVDZ basis set 77 Figure 4-26: Elastic and charge transfer differential cross sections for H"*^ + He. ENDyne results at 1500 eV with a pVDZ basis set 77 Figure 4-27 : Elastic and charge transfer differential cross sections for H"^ + He. ENDyne results at 5000 eV with a pVDZ basis set 78 IX Figure 4-28; Elastic differential cross section for H+ + He. Energy: 5000 eV. Basis set: pVDZ. Solid line corresponds to the semiclassical corrections to the ENDyne results. Open circles are experimental data 87 Figure 4-29: Dynamical trajectories for H'*’ + H2 collision from ENDyne calculation for (0°, 0°) orientation at 500 eV with a pVDZ basis set 89 Figure 4-30: Dynamical trajectories of the target Hj. ENDyne results for (0°, 0°) orientation at 500 eV with a pVDZ basis set 90 Figure 4-3 1 : Dynamical trajectories for + H2 collision from ENDyne calculation for (90°, 0°) orientation at 500 eV with a pVDZ basis set 91 Figure 4-32: Dynamical trajectories of the target H2. ENDyne results for (90°, 0°) orientation at 500 eV with a pVDZ basis set 91 Figure 4-33: Interatomic distances of the target H2 as function of time. ENDyne results for (0°, 0°) and (90°, 0°) orientations at 500 eV with a pVDZ basis set. (1 a.u. = 0.024195 fs) 93 Figure 4-34: Mbrational spectrum of the target H2. ENDyne results for (0°, 0°) orientation with impact parameter of 1.25 a.u. Energy: 500 eV. Basis set: pVDZ 94 Figure 4-35: Probability for charge transfer in -1- H2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Energy: 5(X) eV. Basis set: pVDZ. . 96 Figure 4-36: Deflection functions for + H2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Energy: 5(X) eV. Basis set: pVDZ 97 Figure 4-37 : Elastic differential cross section for + H2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Experimental data from reference . Energy: 5(X) eV. Basis set: pVDZ 98 Figure 4-38: Charge transfer differential cross section for H"*" -1- H2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Experimental data from reference . Energy: 500 eV. Basis set: pVDZ 98 Figure 5-39: Diabatic and adiabatic potential curves for normal electron transfer. . 102 Figure 5-40: Transition state for structures I and II and global minimum structures. Bond distances in pm and angles in degrees 105 Figure 5-41: Alpha and beta isodensity for structure I at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: Li-C = 197.98 pm (= 1.980 A = 3.741 a.u.), C-N = 115.02 pm (= 1.150 A = 2.174 a.u.), Li-N = 192.59 pm (= 1.926 A = 3.639 a.u.). SCF energy: -106.63313290 a.u. = -2901.6368 eV 108 Figure 5-42: Alpha and beta isodensity for structure II at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: Li-C = 214.10 pm (= 2.141 A = 4.046 a.u.), C-N = 115.77 pm (= 1.158 A = 2.188 a.u.), Li-N = 178.79 pm (= 1.788 A = 3.379 a.u.). SCF energy: -106.63909258 a.u. = -2901.7989 eV 108 Figure 5-43: Alpha and beta isodensity of the transition state I-II structure (Cs) at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: Li-C = 208.75 pm (= 2.088 A = 3.945 a.u.), C-N = 115.83 pm (= 1.158 A = 2.189 a.u.), Li-N = 181.89 pm (= 1.819 A = 3.437 a.u.), Li-C-N = 148.49°, C-N-Li = 138.61°. SCF energy: -106.63244750 a.u. = -2901.6180 eV 109 Figure 5-44: Alpha and beta isodensity of the global minimum structure (C 2 v) at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: C-N = 117.46 pm (= 1.175 A = 2.220 a.u.), Li-N = 189.22 pm (= 1.892 A = 3.576 a.u.), C-N-Li = 137.12°. SCF energy: -106.64697039 a.u. = -2902.0134 eV 109 Figure 5-45: Mulliken charge differences at SCFAJHF level in the Li(l)-H-Li(2) molecule as function of relative bond distances (ri - r 2 ) Ill Figure 5-46: Potential energy curves of the Li(l)-H-Li(2) molecule at SCFAJHF level 112 Figure 5-47: Bond distances, ri = distance Li(l)-H and T 2 = distance H-Li(2) in a.u. as a function of time. ENDyne calculation with a H/3— 21G and Li/3-2 1-Kj basis set 114 Figure 5-48: Atomic Mulliken population of Li(l) and Li(2) as a function of time. ENDyne calculation with a H/3-21G and Li/3-2 1-K3 basis set 115 Figure 5-49: Atomic Mulliken charge differences q[Li(l)]-q[Li(2)] as a function of time. ENDyne calculation with a H/3-21G and Li/3-21+G basis set. 116 Figure 5-50: Alpha and beta spin atomic Mulliken populations of Li(l) and Li(2) in the 15% structure as a function of time. ENDyne calculation with a H/3-21G and Li/3-2 1-K5 basis set 117 XU 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 EXPLORING A NEW TIME-DEPENDENT METHOD FOR MOLECULAR QUANTUM DYNAMICS By RICARDO LUIZ LONGO May, 1993 Chairman: N. Yngve Ohm Major Department: Chemistry The electron-nuclear dynamics in chemical processes is described by a method founded on the Time-Dependent Variational Principle (TDVP). In addition, to avoid redundancies in the parametrization of the electronic wave function, coherent states (CS) are used. The equations resulting from this treatment enable us to describe the dynamics in general molecular system, since it does not presume a priori knowledge of the potential energy surface of the molecular system. The theory is called Electron Nuclear Dynamics (END). In order to make this time-dependent method based upon TDVP and CS practical and still realistic, we treat the nuclei in the limit of narrow Gaussian wavepackets (which corresponds to classical nuclei) and use a single determinant to describe the electrons. We name this approach END-SD-FGWP. We use the END-SD-FGWP approach to study the dynamics of charge transfer in • • • XUl the systems Li^-H— Li ^ Li— H-Li^ and Li+-CN-Li ^ Li-CN-Li-". The collision of a proton with H, He, and H2 is also investigated by the END-SD-FGWP approach. Several properties of these collisions are compared directly with the experimental results. Since we are interested in the application of the END formalism to large systems we propose one more approximation, namely, the neglect of the differential diatomic over- lap (NDDO), and introduce the END-SD(NDDO)-FGWP approach. The MNDO/AMl parametrization is used for the NDDO approach. The computational savings are twofold: a) the number of integrals computed are reduced by several orders of magnitude com- pared to the ab initio approach, b) the equations of motion are simplified, and c) the high frequency motions of the core electrons is removed by an effective core. XIV CHAPTER 1 INTRODUCTION Time-dependent methods for solving the Schrodinger equation were introduced in quantum mechanics since its beginning. For instance, already in the 1930s Dirac sug- gested the time-dependent Hartree-Fock method [1]. However, due to the theoretical and computational complexity of time-dependent methods, their development and ap- plication to molecular systems have been delayed by several decades in comparison to time-independent methods. Despite these difficulties, time-dependent methods for molec- ular quantum dynamics have gained in importance in the last two decades. Certainly, the interest in time-dependent solutions of the Schrddinger equation is very much due to the development in recent years of novel experimental techniques that make it possible to measure quantum state-specific properties of rather complicated systems and to probe the dynamics of chemical reactions at the molecular level for very short times. However, a solid theoretical understanding of non-stationary processes is highly desirable, and due to the development of new computational techniques and powerful computer hardware, assumptions, hypotheses and theories describing these processes can now be explored. The time-dependent approach to quantum dynamics offers a number of advantages over time-independent methods. For example, in a reactive scattering process, plots of the particle density as it evolves throughout the collision are easily interpreted using simple concepts of the underlying physics. These plots can be used to improve one’s understanding of mechanisms that are important during a reaction or in other physical- 1 2 chemical processes. From a theoretical point of view, the time-dependent approach makes it straightforward to include the effect of intense laser fields or external time-dependent heat baths, where the dynamics is treated approximately by classical mechanics, statistical or stochastic methods. It is also well adapted to improvement in the dynamics for those degrees of freedom that are described classically, by the wavepacket propagation, as well as those given by explicit quantum mechanical treatments. As a result, the END and other time-dependent approaches offer great flexibility for the choice of various approximations and is also capable of a dynamically accurate description. It is also possible to use within these approaches some of the experience acquired for decades with time-independent methods to provide a practical and realistic description of the chemical system while at the same time keeping an accurate description of the dynamics. It also allows the blending of classical, wavepacket, quantum operator, state expansion and numerical techniques necessary for solving complex dynamical problems. We present a time-dependent method for solving the Schrddinger equation and explore its capabilities in describing charge exchange collisions, molecular vibrations, and electron transfer dynamics. The method presented here attempts to describe the simultaneous electron and nuclear dynamics. It has been developed during the last half decade and has the very nice feature of not requiring potential energy surfaces of the molecular system. It allows both electrons and nuclei to be treated quantum mechanically. The dynamical equations of motion are obtained using the time-dependent variational principle (TDVP) and the wavefunctions are parametrized using coherent states (CS). The combination of TDVP/CS gives rise to a new time-dependent approach to molecular dynamics called 3 Electron Nuclear Dynamics (END). The END formalism provides dynamical equations that look simple and are easy to interpret. The computer implementation of this approach is not trivial, and in order to make it practical some approximations are necessary. Based upon the experience acquired with time-independent methods, the nuclei are treated classically or, more appropriately, in the limit of narrow Gaussian wavepackets, and the electrons are described by a single determinant. The theory for quantum nuclei and multi-determinantal description of the electrons has already been developed [2], but its implementation is not yet done. So we think it would be interesting and fruitful to explore in some detail the single determinant/classical nuclei model. Certainly, this simplest model for the END formalism may be inadequate for describing some problems, mainly the ones involving many product channels (multi-reference character) and tunneling dominant processes [3]. However, since the END formalism does not require the knowledge of the potential energy surface(s) of the molecular system it can be applied to a large variety of systems and dynamical problems in an unified manner. For instance, it can be used to study scattering problems (total, differential and state to state reactive cross sections), charge transfer dynamics (rates and mechanisms), molecular vibrations (interpretation of transition state spectroscopy experiments and prediction of reaction rates). It permits us to follow chemical reactions, predict and interpret Raman spectra, and many other dynamical processes. Future implementation of time-dependent field in the molecular Hamiltonian would make the END formalism able to describe photodissociation processes, laser chemistry, dynamical nuclear magnetic resonance (NMR), etc. It should be noted however, that a fully ab initio treatment of large molecular systems using 4 the END/single determinant/classical nuclei approach is not yet practical. Therefore, in addition to exploring the ab initio version for describing reactive collisions, molecular vibrations, and electron transfer in small systems, we propose to employ the neglect of the differential diatomic overlap (NDDO) approximation [4]. It is our expectation that the NDDO approximation on top of the END/single determinant/classical nuclei model will decrease the computational demand in many ways. a) The number of molecular integrals and their gradients involved would be several orders of magnitude smaller than the ab initio version. b) Since it is commonplace to use effective potential for the core electrons it should, in addition to decreasing the number of degree of freedom (differential equations needed to be solved) also eliminate the high frequency motions of the core electrons. This will allow the integrator of the differential equation to take larger steps decreasing by several orders of magnitude the time to solve the differential equations. c) In general, minimum basis sets are the working frame for the NDDO approximation, which decreases significantly the number of equations of motion. The outline of this work is as follows. Chapter 2 describes some of the available time-dependent methods for molecular dynamics. The END formalism is introduced in Chapter 3 with the NDDO approximation and concluded with the END-NDDO equations. In Chapter 4 we present the results, analysis and discussion of the reactive scattering of on H, He, and H 2 . Chapter 5 is devoted to the study of the intramolecular charge transfer systems y-CN-Li ^ Li-CN-Li-" and y-H— Li ^ Li— H-Li-". The conclusions are presented in Chapter 6. CHAPTER 2 TIME-DEPENDENT METHODS The existing time-dependent methods that are relevant to our work are described in this chapter. The aim is to provide a general view of these methods with emphasis on performance, limitations, and comparisons among them. It seems that the best way to give an overview of these time-dependent methods is to provide some kind of classification that compare their main assumptions. In order to provide such a classification, the concept of potential energy (hyper)surface (PES) is fundamental. For practical purposes a PES is a function that gives the value of the electronic energy of a group of N interacting atoms with respect to their 37V — 6 relative position coordinates. For each electronic state of the molecular system we have one PES which contains information about the isolated reactant and product species, their long-range interactions, their molecular deformations which lead to the breaking and forming of chemical bonds, the geometries and properties of transient reaction intermediates, and energy barriers for going from one of those intermediates to another. As we can see, the task involved in obtaining such a (hyper)surface for polyatomic molecules is tremendous. Highly correlated electronic structure calculations with very large basis sets are necessary in order to provide an accurate description of the different regions of the PES. Once the PES has been mapped, a fitting procedure will be employed to give the functional form of the PES. This fitting procedure is very time consuming and tedious, and it is still more an art than a scientific procedure. The functional form of a PES can be quite complicated. As evidence for that 5 6 we can mention that even for a triatomic molecule like H3, an attempt to use the well developed algebra program (REDUCE) to evaluate the first, second, third, and fourth derivatives analytically failed due to the algebraic complexity [5]. We classify the time- dependent methods into two categories a) PES based methods, and b) PES free methods. PES Based Time-Dependent Methods Once the potential energy surface for the molecular system has been determined, there are basically three approaches to perform dynamics on this surface, namely, a) the classical trajectory method, b) the semiclassical method, and c) the quantum time- dependent method. The given order a), b), and c) is related to the increasing theoretical and implementation difficulties, computational demand, and accuracy. The widely used classical trajectory approach [6-8] employs Newton’s equation of motion to treat the dynamics of the nuclei on the PES. This approach is easily implemented and the interpretation of the results are simple. However, the computational demand can be large for polyatomic systems, and a good statistical analysis program is necessary. A detailed classical trajectory simulation requires an extensive sampling in the phase space (50 000 to 100 000 trajectories for a triatomic system) for each initial condition (collision energy, vibrational and rotational initial states). This computational demand is not the only drawback of the classical method. It neglects quantum effects like tunneling, interference, penetration (generalized tunneling), zero-point motion, resonance behavior, etc, which can be fundamental for chemical encounters. These deficiencies represent an inherent obstacle for the classical trajectory approach to provide qualitative or quantitative results for chemical processes. Despite these problems, the classical 7 trajectory approach has been applied to a widely variety of dynamical problems including scattering of atoms or molecules by solid surfaces [9] and bimolecular reactions [8], In most of these applications (semi)empirical potential energy surfaces are used, such as DIM (Diatomic-in-Molecules), MIM (Molecule-in-Molecules), LEPS (London-Eyring- Polany-Sato), etc [10, 11]. Semiclassical methods have been designed to overcome some of the flaws in the classical description. They assume that the quantum effects can be treated as corrections to the classical motion, based on the fact that the nuclei are relatively heavy compared to the electronic mass. The semiclassical approaches have, in principle, the advantage of incorporating some of the important quantum effects while maintaining the simple implementation. However, the generalization of the semiclassical methods is not simple, and, as result, each system requires a special implementation, and the validity of the approximations is often not clear. The perhaps most developed semiclassical dynamical method has been proposed by Heller [12—15] and employs the Gaussian wavepacket approximation. That is, the quantum wavepacket is constrained to be of a Gaussian form. This semiclassical method has been widely used in many applications with some success [12-15]. Another well-developed semiclassical method that uses the so called eikonal approach to dynamics has been proposed by Micha [16]. Comparisons of the results provided by these two semiclassical methods, namely the Gaussian wavepacket and the eikonal approach, have shown that they differ only in some minor details [17, 18]. These semiclassical methods have been applied to describe photodissociation processes, reactive scattering and bimolecular reactions [15]. 8 Among the various quantum time-dependent methods on PES two have been widely used, namely, the ’exact’ time-dependent quantum mechanical and the quantum time- dependent self-consistent field (TDSCF) methods [19, 20]. Even the so called ’exact’ time-dependent methods do use some approximations, which lead to errors that can, in principle, be made as small as desired. The ’exact’ time-dependent methods can be classified according to the representation of the wavefunction and the algorithm for the time evolution. We distinguish two approaches to wavefunction representation, namely, discretization on a grid of points, sometimes referred to as DVR (discrete variable representation), and expansion in a basis set. In general, the numerical cost of the discretization and basis set expansion representations scales as 0{V^) and 0{N^), respectively, where V is the volume of phase space containing the system and N is the number of expansion functions. There are basically four algorithms in use for the time evolution: SOD (second-order differencing), SO (split operator), SIL (short-time iterative Lanczos), and CE (Chebychev expansion). For a detailed presentation and comparisons of the different ’exact’ time-dependent methods, see references [20, 21]. Applications of the ’exact’ time-dependent methods have been restricted to at most three degrees of freedom and include mainly photodissociation processes, absorption and photodetachment spectra, reactive scattering (atom-diatom), and unimolecular reactions [22]. The TDSCF method assumes variables separation such that each degree of freedom is treated in the mean field of the remaining degrees of freedom. In the static SCF this approximation, also called the independent particle model, leads to the lack of correlation between the degrees of freedom of the system. However, in the case of TDSCF method energy is allowed to flow 9 between degrees of freedom leading to some correlation among them. As a result, the TDSCF scheme is expected to be accurate for calculations of such averaged properties as lifetimes, mean fragment energies, energy flows, etc. In fact, the dissociation dynamics of the l 2 He van der Waals complex has been examined by ’exact’ and TDSCF methods. It has been shown that for this system the TDSCF scheme gives results in excellent agreement with ’exact’ method when the calculated properties are obtained during a period of time much larger than the vibrational periods and so long as the excitation level is not too high. The TDSCF approach also reproduces well the dephasing events in the dynamics and gives very good results for the autocorrelation function [19, 20]. We can summarize this section by saying that the greatest drawback of these methods (classical, semiclassical, ’exact’, TDSCF) is the need for one or more potential energy surfaces. In the case of classical methods the lack of any quantum effects can cause problems, the semiclassical approaches are too system specific, and the ’exact’ schemes are very numerically intensive and are restricted to problems with few degrees of freedom. No doubt these methods have their advantages and their contribution to the understanding of time-dependent processes must be emphasized. PES Free Time-Dependent Methods We focus in this section on time-dependent methods that do not require potential energy surfaces. Among these methods we intend to discuss the time-dependent methods that fall into these three classes: a) close-coupling, b) molecular dynamics, and c) time- dependent Hartree-Fock. 10 Close-coupling Methods The close-coupling methods [23-25] consist in projecting the time-dependent Schrodinger equation on a basis. The electronic states are described by a Cl (config- uration interaction) expansion where the time-dependence is put into the Cl coefficients. The nuclear dynamics are not treated explicitly and are constrained to follow a priori trajectories, such as linear or Coulombic trajectories. The description of the electronic states by a Cl wavefunction imposes a heavy computational demand restricting the application of the close-coupling methods to simple one electron (or pseudo one-electron) systems. Recent extensions to two-electron systems seems not to be very practical [26]. The prescribed trajectories for the nuclei limit the application of the close-coupling methods to high energy (above 1.0 keV) collisions. It also makes the application of these methods for the calculation of differential cross sections other than the smallest scattering angles quite doubtful. Molecular Dynamics Methods Molecular dynamics (MD) methods assume that the atomic dynamics is governed by the laws of classical mechanics. That is, MD methods neglect quantum effects on the nuclear motion and the electrons are considered as being always in the ground- state corresponding to the instantaneous nuclear configuration. As a result the explicit knowledge of the interatomic potential is required. Obtaining this interatomic potential is 11 a difficult many-body problem and different ways of calculating this potential distinguish the different MD methods. The most widely used MD methods are the ones that employ (semi)empirical force fields for describing the interatomic potential [27]. These MD/force-field methods can be combined with perturbation methods to provide thermodynamical properties [28], in addition to structural (geometric and conformational) information. They can be applied to a great variety of systems including large biological molecules [29]. Some of the other MD formulations evolve computing the interatomic potential by solving the Schrddinger equation in some approximate manner (e.g. ab initio or semiempirical Hartree-Fock (HF), density functional (DF) theory, etc) for each nuclear configuration during the time evolution [30, 31]. Time-Dependent Hartree-Fock Methods Despite the fact that the time-dependent Hartree-Fock (TDHF) method was suggested by Dirac already in 1930s [1] only recently has this method been developed and applied to reactive and inelastic scattering [32-35] and molecular problems [36, 37]. In TDHF the quantum dynamics of the electrons represented by a single Slater determinant is coupled to a classical description of the nuclear dynamics in a self- consistent fashion. The differences between some of the TDHF implementations are related to the way the Hartree-Fock equation is approximately solved and to the way that the translation of the basis set with the nuclei is handled. This is related to the electron translation factor (ETF) problem. 12 The NDDO approximation has been used to simplify the TDHF equations and applied in the calculation of molecular properties, such as atomic charges, structure, dipole moments, etc, during time evolution of LiH, H 2 O, and CH 2 O molecules [36]. The density functional theory (DFT) using the local density approximation (LDA) has been employed to compute the Hartree-Fock potential for the TDHF method. This approach has been used in the dynamics of Na 4 cluster for the singlet and triplet electronic states [37]. Both approaches, namely, TDHF-NDDO and TDHF-DFT, do not consider the effects of the translation of the basis set with the nuclear motion upon the electron dynamics. This means that the electron translation factors are not properly taken into account by these approaches, neither are the electron-nuclear couplings. The proper derivation of TDHF equations with the inclusion of the electron translation factors has been done using the density matrix instead of the orbital coefficients as the time -dependent quantities [34]. However, no applications of this method have yet been made, and only after the neglect of the electron translation factors (RTF’s) was this approach applied to the H'*’ + H collision [34]. It is expected that in a short time applications including the ETF’s will appear and we shall be able to see the magnitude of the errors caused by the neglect of these effects. In summary: 1 . PES based time-dependent methods have the drawback of requiring a priori knowl- edge of the interparticle potentials. This requirement limits the applicability of these methods to small systems (just a few degrees of freedom) and low energy processes (high energy implies in a large number of excited PES’s), 13 2. PES free time-dependent methods are very promising, but general computer codes, extensive testing of the approximations, and a solid foundation of the formal devel- opments are still required. As a conclusion to this chapter we can say that there are still a lot of formal and numerical investigations to be performed on time-dependent methods for molecular quantum dynamics. In the next few chapters we shall be able to contribute to the investigation of a new time-dependent method. CHAPTER 3 THE THEORETICAL APPROACH In the last chapter we expressed our concern for potential energy surface (PES) based time-dependent methods being limited to systems where these PES’s are known. In this chapter we present the equations for a recently proposed time-dependent approach to molecular dynamics called Electron Nuclear Dynamics (END). The time-dependent variational principal (TDVP) is used to obtain an approximate time evolution of the molecular system. In order to obtain the equations of motion it is necessary to decide how the electrons and nuclei will be treated. From the experience with time-independent methods and some of the time-dependent methods, like TDSCF and TDHF, it seems that treating the electrons at the single determinant level with classical nuclei should provide a very realistic approach for a wide variety of chemical processes. Since the classical nuclei approximation can be obtained as the narrow width limit of a frozen Gaussian wave packet treatment, we denote the method resulting from these approximations as END-SD-FGWP. That is, END stands for Electron Nuclear Dynamics from the approximate dynamical equations obtained from the TDVP, SD is the Single Determinantal description of the electrons, and FGWP is the classical, or narrow width limit of Frozen Gaussian Wave Packet representation of the nuclei. The coherent state (CS) for the electrons is the Thouless single determinant representation, which gives a suitable parameterization. In order to make the END-SD-FGWP approach suitable for studying dynamics of 14 15 large molecular systems, we employ the neglect of diatomic differential overlap (NDDO) approximation. We have divided this chapter into five sections, namely, the END formalism, the END-SD-FGWP method, the NDDO approximation, the END-SD(NDDO)-FGWP model, and implementation. In the first section the CS, the TDVP, and the dynamical equations are introduced. In the NDDO section we have a subsection where we present the Austin model 1 (AMI) realization of the NDDO approximation. The END Formalism We present a formalism that treats the electron nuclear dynamics without the need of a priori knowledge of potential energy surfaces or interatomic potentials. In other words, in this formalism only the Hamiltonian (e.g., Coulombic) and the initial conditions (initial nuclear positions and momenta, for instance) are needed for the calculation to proceed. This formalism relies upon the time-dependent variational principal (TDVP) of a quantum state describing the system through a group theoretical coherent state (CS). In the next two subsections we present the dynamical equations for a molecular system derived using the TDVP and the CS representation. The actual derivation has been presented elsewhere [38, 39] and will not be repeated here. Coherent States We now present the parametrization of the wavefunction in terms of coherent states (CS). A coherent state [38, 40] of a molecular wave function |C) is a continuous function of the labeled set of electronic and nuclear parameters (. This set of parameters has a 16 positive measure d( such that the resolution of the identity j ICHCWC = 1 (3.1) holds. This choice of representing the wave function is very important since it removes any parameter redundancies and leads to a simple interpretation of the dynamical equations. The removal of the redundant parameters is essential to avoid complications with the phase space that could lead to artificial singularities and difficulties with integration. A coherent state representation also leads to a natural division of the parameters into (gen- eralized) canonical coordinates and conjugate momenta. Thus, they form a generalized phase space with the expectation value of the energy being the Hamiltonian function for the time evolution of the parameters. As a result, the equations for the CS parameters form a classical Hamiltonian system facilitating their physical interpretation. The generalization to a multi-determinantal description would follow the same general approach and has been done elsewhere [2] . We are interested in those coherent states that are related to compact Lie groups and to a single determinant. Using the group theoretical approach we can eliminate the redundancy of the single determinant representation. Let us assume that we have an orthonormal spin orbital basis set made up of K vectors and N particles which will be described by a single determinant. This determinant can be found as a unitary transformation acting on the orbital expansion coefficients of some reference single determinant |^o)- We define this reference state 17 l^o) = n^^K^) (3.2) /=1 where b] (6,) are the fermion creation (annihilation) operators of the spin orbital basis used. The reference state is a lowest weight state for the irreducible representation [l'^0(^“^)] of the group U{K) with generators bjbj because the weight or number operators b]b{ have eigenvalues 1 for * = 1, . . . , iV and 0 for i = iV + 1, . . . , ii'. As we can see we have divided our set of K spin orbitals into two subsets N and K — N. We now denote any quantity q as: 1. q if q refers to the set of all K spin orbitals; 2. q* if q refers to the subset of N spin orbitals; and 3. 9° if q refers to the subset of K - N spin orbitals. In the context of a determinantal state this division of the set of K spin orbitals into the two subsets: N and K — N takes the meaning of occupation, where the subset N is called occupied and K — N is the unoccupied or virtual subset. In the general case, however, these subsets N and K — N will be bullet and open circle, respectively. Also, the following notation will be used when dealing with matrices: / “11 ••• ^IN “l(7V+l) ••• ^'iK \ u = u ^N1 (AT+1)1 u n K\ u u n NN u It KN u N{N^\) (TV-fl)A^ ^(Ar-hi)(AT+i) u A'(iV+l) u NK u u KK (3.3) ^NxN If" ^K-NxN ^NxK-N \ _ ^K-NxK-N • • • / 18 We now assume that there is an t/ matrix which transforms the spin orbital basis to the occupied and virtual spin orbitals, such that, all redundancies due to unitary transformations among the occupied (virtual) spin orbitals are removed, that is. ( 3 . 4 ) Then, the coherent state representation of a determinantal state takes the form: N 1=1 = n ( E + E 1 1'"“^) i=l \i=l ji=7V+l N / N K N = n I E 1 + E E f'yikVtr' I f'l* I \o“c) j=N+l k=l i=l \/=l ( 3 . 5 ) “HK+ E i=l V >=A^+lib=l =a n('+ E i=i \ y=jv+iJb=i / 1=1 19 and since we work with an unnormalized coherent state we have that N K =n n t=l j=N-\-l (3.6) N K 1=1 j=N+l where ^ The unnormalized determinantal wavefunction can be expressed in terms of the orthonormal spin orbital basis set and the dynamical variables {zji,Zji} in the following way where the dynamical spin orbitals K Xi = V’i + (1 < * < N) (3.8) j=N+i are nonorthogonal, but linearly independent. The corresponding unoccupied dynamical orbitals may be chosen as and although mutually nonorthogonal they are orthogonal to the occupied space. \z) = det{x,} (3.7) N (3.9) 20 The Time-Dependent Variational Principle We now give an outline of how the time-dependent variational principle (TDVP) can be used to provide us with an approximate time evolution of the quantum mechanical system. There are several forms of the TDVP [1, 41, 42], but they are all equivalent when the trial wave function is described by complex parameters and is analytic in this parameter space [43]. The TDVP approach consists in minimizing the variations of the quantum mechanical action functional A, that is. f[(Cl|(IO) dt «ci)ic) - (ci^io Mcicr' = 0 (3.10) where \Q is quantum state of the system labeled with a set of electronic and nuclear parameters ( and H is the Hamiltonian of the system. The TDVP approach leads then to the following equation for the time-dependent state |C) (3.11) This is the time -dependent Schrodinger equation, provided that \6() is allowed to vary over all Hilbert space and the right-hand side vanishes. It can be shown [38] the right- hand side of equation (3.11) vanishes when we consider explicitly the overall phase factor of the time-dependent state |(). 21 The dynamical equations resulting from this approach are [38, 39] E = 0 dE d(o (3.12) which in matrix form becomes C 0 0 -C* c c* / dE ( w \ M V ac (3.13) and the elements of the metric matrix defined as ^ a^InS ^010 — p,. (3.14) dQ9C0 where the energy of the system is given by the expectation value of the Hamiltonian H as follows E{C\0 = {C\H\0 (CIO (3.15) and the norm is 5(C*,C) = (CIC). (3.16) It also should be mentioned that in this formulation the global phase factor exp (^ 7 ) of the wavefunction l^') = exp (n) 10 (3.17) has been factorized from the evolving state and can be obtained by a simple quadrature integration of the following differential equation 7 4E(<. a a,n5(c,o_^.ai^ \ CK dCa d\nSiC,0 dCa da (3.18) -^(C*,C) 22 during the evolution. This global phase should be important in applications where the time-correlation functions need to be computed. The END-SD-FGWP Method We have so far obtained dynamical equations that describe the dynamics of a general quantum molecular system. As mentioned before, in order to obtain equations of motion for the electron nuclear dynamics that are practical and realistic, some approximations are needed. We choose a model for which (1) the electrons are described by a single Slater determinant, and (2) since we would like to make use of the available computer codes for molecular integrals we use a truncated spin orbital basis without electron translation factors (ETF’s). In other words, the basis set is independent of the nuclear velocity (or momentum), and (3) use a classical description of the nuclei. In mathematical terms, this means that when employing the TDVP approach the limit of narrow Gaussian wave packets is taken. As a result, the nuclear parameters become the positions and momenta of the nuclei. TTiese approximations allow us to obtain the equations of motion for END-SD-FGWP method. The Equations of Motion Taking the limit of narrow Gaussian wave packets or the classical description of the nuclei permit us to define the dynamic variable for the nuclei as [38, 39] Zk = Rk + iPk (3.19) 23 where R and P are the nuclear coordinates and momenta, respectively. Since we have assumed that the dependence of the basis on the nuclear parameters is only through the nuclear positions, the equations of motion for the electron and nuclei become [38, 39] / iC 0 iCR 0 ^ (dEldz*\ 0 j -iC* 0 • ♦ dE/dz iC R -iCl Crr -I k dEjdR 0 0 I 0 ) \p) \ dEjdP J where d‘^\nS{z\R\P\z,R,P) dz*dXk R'=R,P'=P and CxkYi = -2Im d^\nS{z\R',P\z,R,P) dX'dVi R'=R,P'z=P with X and Y standing for R or P, with the norm S being (3.20) (3.21) (3.22) 5 = S{z'\R',P',z,R,P) = {z'\R!,P'\z,R,P) = det(/* + z'^z) (3.23) In order to implement the END-SD-FGWP formalism developed so far we need to specify how to construct the spin orbital basis. Usually, in electronic structure calculation, a set of K localized atomic orbitals are used to define the spin orbital basis. These atomic orbitals are given as linear contractions of Slater (STO) or Gaussian (GTO) type orbitals and are, in general, centered on the nuclei. These orbitals are also, in general, non- orthonormal. 24 Let VL be a matrix that transforms the non-orthonormal localized atomic orbital basis to the orthonormal spin orbital basis tl^ = <f>W (3.24) f W* W' \ (-/-• = I,) It should be noted that the bullet and open circle notation lose their meaning of occupied and virtual in the atomic basis. However, we are still using occupied and virtual for bullet and open circle even in the atomic basis representation. Using the fact that the only transformation that will have any effect in a determinantal state is the one that mixes unoccupied orbitals into the occupied orbitals, we can have that W" = 0. This means that the reference spin orbitals = (f>*W' + <I>°W° (3.25) are constructed by first orthonormalizing the bullet (occupied) atomic orbitals among themselves, and then orthonormalizing the open circle (unoccupied) atomic orbitals to the bullet (occupied) space and among themselves. This has the effect that the occupied (bullet) space is the same whether defined in the atomic or the spin orbital basis. That is the reason why we keep using occupied and virtual for the bullet and open circle spaces even in the atomic basis representation. We have then created a recipe of how to construct the reference state from the atomic orbitals basis. The construction of the spin orbitals involves what is called, in electronic structure theory, the MO (molecular 25 orbital) transformation. This is a very demanding step in terms of computational effort, since it scales as > K^, and the atomic integral calculation scales as w K^. As a result, despite the fact that the equations for the metric and dynamical variables get a little bit more complicated, it is of great computational advantage to work in the atomic basis. The spin orbitals that are delocalized over the entire molecular system are also called molecular orbitals. We denote the set of parameters of the coherent state with respect to the reference (spin orbital, or molecular) basis as and with respect to the atomic basis. The following transformations exist 2“* = -I- ( 3 . 26 ) Z^nol ^ and also the basis field operators transform as b = wU ( 3 . 27 ) where {a, a’} are the field operators in the atomic basis representation. The coherent state is thus equal to N K t=l >=^ + 1 N K = q;“* n «=1 j=N+l ( 3 . 28 ) at\at = 1 ^.-) 26 and as before, we work with the unnormalized coherent state k“r=nK+ E “pji »=i \ Note that because the anticommutation relations among are not canonical^ it is not possible to write the coherent state as an exponential with the atomic parameters and However, it is still possible to express the dynamical orbitals K Xi^^i+ 4>jZji, (l<i<N) (3.31) j=N+l ^ I vac) (3.29) N Xj = - XZ 4-V’.^ (A" + 1 < i < K) 1=1 in the same simple form K Xi = ^i+ X] (1 < * < N) j=N+l (3.32) (3.33) where N Xj = 4>}~Y. {N + l<j< K) 1=1 V = - (A* -I- A'z)"^ = - (a* + A'^) (a' + (3.34) (3.35) 1. the creation and annihilation field operators in the spin orbital basis are canonical since [6, 6^].^. = 1, however, in the atomic basis [a, = A the anticommutation is not canonical, since the atomic overlap matrix A has elements A.J = {4>i\<l>j) = J 4>i{x)4>j{x)dx and {<?!>, A' being non-orthononmal makes A different of the unit matrix (3.30) 27 with the overlap being given by S = S{z'*,E',P\z,R,P) = {z'\R',P'\z,R,P) = det(A* + a' z + z'^ a" + z'^A°z) ( 3 . 36 ) We now define the following matrices A* A' /• A° ^ = A* + z^A'^ + a' z + z^A°z ( 3 . 37 ) and A° = (t. =A» + ..A' + AV + rAV ( 3 . 38 ) which naturally occur in many expressions. Note that A* and A° are dependent upon {z,z^} and respectively. The metric C can now be expressed in the atomic orbital basis as C.. = dHnS - = -(a-^(a' + 2^A°))^.^(a-^(a' + 2+A°)) nk ( 3 . 39 ) c...= d'^lnS = ((aV + A->)a»-> (vA' + A->))^^{A->)^„. and Cz'z' dHnS The mixed derivatives dz*V r involved in Cr are d = (a%^ + A°)a°-^(u ( 3 . 40 ) ^ + + <3.41) ( 3 . 42 ) 28 And the double derivative with respect to the nuclear coordinates is given by Crr = V g^\n{z' ,Fi\z,R%i^^^R>^R = Tr A— ^(/* ( 3 . 43 ) /• •— 1 / /• z where P-'M) We have introduced the notation of a vertical bar to indicate that the derivative of the overlap matrix is to be taken with respect to the R dependence of one side only. Now that explicit expressions for the metric have been given, we turn to the expression for the energy derivatives appearing in the right hand side of the equation of motion. In order to obtain explicit expressions for the derivatives of the energy it is appropriate to define the one-particle density matrix (also called 1-density or charge- bond-order matrix) as follows r = (^J)A-'(/* 2'). which was obtain from the kernel Tji{z*,z) = {z\zy {z\blbj\z) and has the familiar block form ( 3 . 45 ) ( 3 . 46 ) ( 3 . 47 ) 29 The energy derivatives then become dE{z\z) + /°)(/i + Tr(Va6;„6r)J = ^ (a'^i;+ + A°) ( u P )F ^ ^ A*" where we recognize the Fock operator F = h + Tr(yoj,;aj,r)^ = h-\- Tr([(aa|66) — (a6|fea)]r)^ . The energy derivatives with respect to the nuclear coordinates are ki 1=1 (3.48) (3.49) F{z*,R',P\z,R,P) I l^ZkZie^(^Rk-Ri) / \ / \ = -5 E + Tr(v^.Ar) - Tr(v^-,Arftr) (3.50) + 5Tr(Tr(v^^Ki;,»r) j) - Tr('Tr(v^^Ary,i.,tr) J and since the basis set is independent of the velocity (or momentum) of the nuclei the energy derivatives with respect to the momenta become VpE{z\R',P\z,R,P) = fL R’=R,P'=P Mk (3.51) We now have explicit equations for all terms in the equations of motion / iC 0 iCR i 1 ( OF !dz*\ 0 • -iC* -iC*p 0 z* dEldz iCn -iCp Crr -7 R dEjdR \ 0 0 I 0 \pj \ dEIdP / (3.52) 30 and we can try to simplify the main equations in such a way that they become efficient to implement in a computer program. Starting with the equation + iCRR dE dz* ( 3 . 53 ) we can use the expressions for Cg>z, Cr, and dEjdz* to further reduce it to /• 2: = ( -2 1° )F /• ( 3 . 54 ) The right hand side has a very simple interpretation: the Fock operator acts on the occupied states and is then projected on the virtual space. Thus only the projection on the virtual states gives rise to a change in the 2 coefficients. This is also the way that this particular equation of motion is implemented since it does not involve any matrix inversion, which would be necessary if the metric did not have this particular structure. In the equation of motion that involves the derivatives of the energy with respect to the nuclear coordinates, E{z\H',P\z,R,P) I = 4E ^ ZtZ,e‘‘{Rt - Rl) /=! \Rk - Ri\^ + Tr(v^-/r)-Tr(v^^Ar/.r) (3.55) + ^Tv[Tr(v^Vab-,abT)r^ - Tv ( tv ( v ^ATVab.,abr) T 31 we can see that it will be necessary to compute the derivative of the 1-density matrix r= (^*j(^A* + A'z + z^A'^ + z^A°zy\r z^) _M 1a— !(/. /• Z V(/* \i' ^') with respect to the nuclear coordinates, which can obtained as follows = ; r )(V^.A-)(/' r\..- z (/• .')(v^.a)(^Qa-‘(/' ^') (3.56) (3.57) = r(v^,A)r. As a result, in order to compute the energy gradient with respect to the nuclear coordinates all we need are the gradients of the one- and two-electron integrals (which are available from many quantum chemical program package) in addition to the derivatives of the atomic overlap matrix. The NDDO Approximation The NDDO (neglect of diatomic differential overlap) approach employs the following approximations [4]: 32 1. Only the valence electrons are considered explicitly in the calculations. The nucleus and inner shell electrons are replaced by a fixed core function. This is the so called core approximation [44]. 2. A minimum basis set is used. This basis set is constructed with Slater type orbitals (STO’s). The advantages of using STO’s will be clarified later. The STO’s are products of a radial function and a normalized real spherical harmonics with quantum numbers n, /, m, Xnlm = Rnl{r)Yi^{e,(p), (or ,'\n+l/2 (3.58) R ,(r) = -n-1 -Cnir ' [(2n)!]V2 ' • Here Cn/ is the orbital exponent of the Slater atomic orbital. 3. The overlap integrals are neglected, except when they appear in the one-electron resonance integrals (one-electron two-center integrals). As a result, the NDDO approximation consists of applying the following relation (1) ^ SabX^{1)x^{1) . (3.59) Since we are treating only valence electrons and using STO’s we have that (3.60) the atomic overlap matrix is approximated by the identity matrix. 4. Two-electron integrals then become ^^C<Bc(#^(l)tf(2)|^^(l)^f(2)) (3.61) 33 We have chosen the NDDO model [45, 4] since it is considered the lowest level at which there exist a basis set for which the zero differential overlap (ZDO) [46] approximation is valid, or almost valid [47, 48]. The NDDO model is rotationally invariant so that it is not necessary to perform any spherical averaging of the atomic orbitals with angular momentum larger than zero. As a result, unlike the other ZDO models, such as CNDO (complete neglect of differential overlap) and INDO (intermediate neglect of differential overlap), the NDDO model includes orbital anisotropies. In addition, in the HF-NDDO method the computational dominant step is the diagonalization of the Fock matrix. We now present some justifications for the approximations adopted and their advan- tages: (1) The core approximation in the NDDO model is reasonable because the inner electrons are tighdy bound and are, therefore, unlikely to be significantly perturbed by changes in the valence shell. Consequently, most chemical processes can be accurately described by valence electrons only. Also, there are solid theoretical foundations for the core approximation in terms of an effective Hamiltonian [49]. (2) The use of STO’s as building blocks for the molecular orbitals ensures that only a few of these STO’s are necessary to give a good description of neutral molecules or positive ions. This justifies the approximation of minimal basis set. Only in the cases of anions or hypervalent compounds must the atomic basis set be more flexible to accommodate the extra charge density or extra bond(s). Adding diffuse or polarization atomic orbitals is one way of solving this problem. Another is to allow the orbital 34 exponents to vary with the atomic charge. This provides an appropriate description for molecular anions [50], and also allows us to keep the minimum basis set formalism. The use of STO’s is useful in the gradient calculation. That is, the differentiation of an STO function with respect to the nuclear coordinate does not yield higher angular momentum functions. In other words, once the basis set is defined in terms of STO’s no higher angular momentum functions are needed in the gradients of the molecular integrals. (3) The neglect of overlap is consistent with the approximation employed in the two- electron integrals. However, keeping the overlap in the one-electron resonance integrals (one-electron two-center integrals) is important since these integrals are responsible for the bonding in molecules. That is, the resonance integrals = (x„(l)l(-iv?-^:?i)|x.(l)), (3-62) which involve an overlap distribution are not neglected. (4) It can be shown [51, 52] that the approximation of the two-electron integrals by the neglect of diatomic differential orbital is correct up to second order in the overlap expansion. It is also possible to evaluate the three- and four-center integrals in terms of overlap and two-center integrals [52]. However, if we parametrize the two-electron integrals with respect to experimental data the errors caused by neglecting the three- and four-center integrals can be minimized. The END-SD(NDDO)-FGWP Method Since the dynamical equations are dependent upon the norm we start by applying 35 the NDDO approximation to S = det(A* + A' z + z'^A" -f z'^A^z). The NDDO approximation implies that A* = r, A' = A" = o, A° = r so that the overlap takes the same form as in the orthonormal basis, S = det(7* + z'^z) as do the A matrices A* = /• + zh A° = 7° + zzK Also, the density matrix is simplified to r= 1^*1 A— ^(7* z^) = (^'’yr + z'z)-'{r zt) (3.63) (3.64) (3.65) (3.66) (3.67) and u = -(A'^ + A°2)(A* + A'z)-^ = -z = -(A* + z^A'^)-^(A' + z^A®) = -z+ As a result, the relevant terms for the dynamical equations become Metric d'^lTiS -((r + z*z)-'z^).,{{r + z<z)-'z\, (3.68) (3.69) 36 (3.70) Im ^ki (3.71) d = (/° + ^zt)-l(_^ /'>)V^-^A|i: (/* + zM K\-i Crr = V ^\n{z' ,Ii\z,R)\^'^^^Il-^R (3.72) = Tr (/• + 2^2)-^(/* V^AI ^ -(r + z'z)-'(r z^)v^A\(''){r + z'z)-\r z «) v ^ a /^‘ Energy Gradients dE{z\z) 8z* = ((/° + 22')-‘(-z r)F{’'\r + zU)-' ki V^E{z\B!,P\z,R,P)\^ "* \R'=R,P>= ^ZkZ,e^(^Rk-Ri) 1 =—y 2 ^ 1=1 l^k - RiP + Tr(Vjj,Ar)-Tr( V^AThr) + 5Tr('Tr(v^^Ki..ir) - Tr(Tr(v^^ArK.,,.ir) F (3.73) (3.74) (3.75) Vpy{z\R',P\z,R,P) = A R'=R,P'=P Mk (3.76) 37 where the Fock operator must also be considered under the NDDO approximation f = A + Tr(v;^S'>°r)^ (3.77) = h + TT{[{aAaA\bBbB) ~ {aAbA\bBaB)]^)a • Caution must be exercised when considering the derivatives of the atomic orbital overlap A with respect to the nuclear position. If we make the NDDO approximation for the overlap, that is, A = /, then take the derivative we obtain no contribution. As a result, both Cr and Crr vanish and we do not have any coupling between the electrons and the nuclei. This is not consistent with the ab initio END equations, since the NDDO approximation should not affect the proper coupling between electrons and nuclei. Also, making the NDDO approximation before taking the derivatives is inconsistent with the analytical gradient method. That is, the contribution of the two-center one-electron (resonance) integrals to the energy gradient would then vanish, which is not true when finite difference method is employed. Consequently, we should first take the derivative of the non-approximated atomic orbital overlap and then perform the NDDO approximation on the resulting expression. AMI Model The AMI (Austin Model 1) approach to the NDDO approximation [53] has been extensively tested and has been shown to reproduce geometry, heats of formation, dipole moments, and ionization potentials very well [53-55]. In addition, the overall perfor- mance of the AMI for computing molecular properties (excitation energies, vibrational frequencies and intensities) and quantities relevant for chemical reactions (enthalpy of 38 reaction and activation energies) is also satisfactory [55]. It, thus, seems appropriate to use the AMI realizations of the NDDO method in our END-SD(NDDO)-FGWP model. Results for dynamical properties like intra- and intermolecular charge transfer, transition state spectroscopy, photodissociation, etc, can be used to check the performance of the END- AMI method. In case this approach fails a different realization of the NDDO method can implemented where the parametrization should also include experimental dynamical properties in the data base. Before we present the AMI model we need to consider first its Fock matrix elements [56, 53, 55]: 1 . Diagonal terms R = u, (tA/iA Bi^A (3.78) B^A Xb 2. Off-diagonal terms on the same atom Bi^A + l:P^^Al'A[^y'Al'A\^^AVA) - {tiAiiA\vAt'A)\ (3.79) + '^'Y^P\B<TBit^AVA\^B(^B)- B^A Xb ^b 39 3. General terms on different atoms I'A <TB where <^b) are atomic orbitals centered in A {B) and P^x is the bond order matrix. In order to complete the definition of the AMI model we need expressions or numbers for each one of the terms that go into the elements of the Fock matrix. That is, !• UfiAHA one-electron one-center terms and are parametrized from spectroscopic data for valence states of atom A and its ions; 2. Za is the atomic number of atom A minus the number of core electrons; 3. The nuclear repulsion energy is given by where af, bf, and cf are parameters; 4. The resonance integral is made proportional to the atomic orbital overlap. The proportionality constant is the average of the parameters; Eab = ZaZb{sasa\sbsb) 4 (3.81) 5. The repulsion integrals are expanded in terms of multipole Mi^ and are approximated 40 by the following expression 00 OO imtn / l =0 / 2— 0 where, Rij = ^AB + + vf,D^ (3.83) with tj = ...,±2,±1,0 depending upon the order of the dipole, p and D are derived parameters, that is, D = /^(C/i) and p = p{U,i„,D). In the case of a minimum sp basis set, there are 22 distinct repulsion integrals. The inclusion of d orbitals increase this number to 450, which make this type of approximation for the repulsion integrals cumbersome to be extended beyond sp basis. The number of parameters in the AMI model ranges from 13 to 16. Around 150 molecules were included in the reference data base used to obtain those parameters. These numbers can give an idea of how difficult is to obtain parameters for a semiempirical method. Implementation The END-SD-FGWP method has been implemented in the ENDyne program package. This program is capable of performing simultaneous optimization of the electronic and nuclear (geometry, for classical nuclei) wave functions. ENDyne also can perform the time evolution of a molecular system state. The input consists of the initial states 41 (electronic and nuclear), number and type of particles (electrons, nuclei, nuclear masses, atomic number), initial and final time of the evolution, the (Gaussian) basis set for each nucleus, and the type of differential equation solver. The program works in the laboratory frame and uses atomic units. As output, ENDyne generates a file, called restart which contains the history (the 2 parameters and phase) of the system evolving in time. Another program, called EVOLVE, is used to extract, plot, and manipulate relevant information like nuclear position, electronic and nuclear momenta, phase, electron density and atomic population, etc, from the restart file. A functional diagram of the ENDyne program can be seen in the Figure 3-1. 42 Figure 3-1: Functional diagram of the ENDyne program. The main program endyne calls getinp and decides the structure and type of calcula- tion. After that, it makes a dummy run and determines exactly how much memory will be needed for the calculation. Then, the rundyn is called and it takes control of the ENDyne calculation. It calls dynini to initialize the calculation for the models (quartic, Morse, and oscillator potentials) and theories (Hartree-Fock, AGP, RPA). It then decides if the 43 calculation is an optimization (frprnm, djpmin), optimization (Jrprmn, dfpmin) followed by evolution {deint, drint, odeint), only evolution(dc, driveb, odeint), or exact oscillator (xctosc). Dynfun is then called and computes the gradient of the energy with respect to the wave function parameters in the case of optimization or the derivatives of all wave function parameters with respect to time for an evolution calculation. The appropriate forces (gradients) for the semiempirical or ab initio model is then called. In the case of the Hartree-Fock theory (END-SG-FGWP) the atomic integrals and overlap and their derivatives are computed every time the force routine {hffrc) is called. The onedrv and twoint routines are drivers for one- and two-electron integral calculations, with averll, averll, and disint being the interfaces with the ABACUS integral package. In the case of AMI model the structure is the same, except that the drivers onedrv and twoint calls some other appropriate interfaces, instead of calling averll, averll, and disint. In addition, the implementation of the AMI method, necessitates the development of equations and source code for computing analytical derivatives of the integrals (repulsion and overlap) and of the energy. The molecular orbital programs AMPAC 2.1 and MOPAC 6.0 based upon the AMI model do not perform analytical derivatives at the SCF (self-consistent field) level. However, analytical derivatives for similar expressions of the repulsion integrals have already been performed [57], which insures that the AMI model implementation is viable. CHAPTER 4 H+ + H, He, AND H 2 CHARGE TRANSFER COLLISION In this Chapter we present results for the electron-nuclear dynamics that takes place when H"^ collides with H, He, or H 2 . Such a study consists of computing the probability of the charge transfer, energy transfer, etc as a function of the impact parameter and, in the case of the H 2 target, also molecular orientation. Plots of the transfer probability, vibrational excitation, and cross sections versus the scattering angle are presented and comparisons with other theoretical approaches as well as experimental data discussed. This chapter has the following subsections: the H'*' + H, + He, and + H 2 results, with a fourth subsection containing general conclusions. Before presenting the results some general considerations are in order. i) All the calculations were performed with the ENDyne program package using the following computers: IBM RS/6000-550, Sun 4/490, Sun 690MP, and Cray Y-MP/432. ii) An Is2s2p basis set was used for both H and He. This basis sets, namely pVDZ for H: [4slp]/(2slp) [58] and 6-31G* for He: [4slp]/(2slp) [59] are listed in Table 4-1. In the pVDZ basis set the s-exponents were scaled by 1.2^ = 1.44, as was done in the original DZP basis set by Dunning [60]. iii) Some of physical properties for the atoms and molecules involved in this work are presented in Table 4-2. iv) Due to the physics of the problem (a proton hitting atomic and molecular targets) the 44 45 initial wave function consisted of a totally localized electron density on the H, He, and H 2 targets. In other word, the starting wave functions are not the ground states of the complex species (H-H)"^, (H-He)"^, or (H-H 2 )‘^. It should be mentioned, however, that the atomic and molecular targets are in their electronic (and vibrational) ground states. The difference in energy associated with these initial wave functions and the ground state energy of the (super)molecule ion species is due to the localization of the electronic density around the targets. Table 4-1: Contraction coefficients c and exponents a for the basis sets used in this work. H/pVDZ“ He/6-3 c a c a Is orbital 0.01969 18.7344 0.02311 38.4216 0.13798 2.82528 0.15468 5.77803 0.47815 0.64022 0.46930 1.24177 2s orbital 0.50124 0.17568 0.29796 1.00000 2p orbital 1.00000 0.72700 1.10000 1.00000 a: Ref. [58], b: Ref. [59]. Table 4-2: Ionization potentials, electron affinities and polarizabilities of H, He and H 2 . Ionization Potential (eV) Electron Affitinity (eV) Polarizability (lO'^"* cm^) H* 13.598 0.754 0.6668 He® 24.587 not stable 0.2050 H 2 ® 15.427 — 0.803 a: Ref. [61]. 46 The + H Collision The study of H'*' + H collisions is of special interest because of its great simplicity but rich physics. This has inspired much theoretical and experimental work involving charge transfer (total and differential cross sections), excitation, alignment measurements, etc. As a result, H'*’ + H collisions provide a good test for any time-dependent method that intends to describe electron and nuclear dynamics. It should be noted that since this is a one electron system the effects of electronic correlation are absent. Consequently, the same methods should be tested for systems with more than one electron, since like in stationary quantum chemistry, the electronic correlation effects can be essential when describing molecular properties. Most of the theoretical calculations on the H'*’ + H system are done by close- coupling methods [25, 24]. Since these methods use prescribed trajectories they are most commonly applied to collision energies above 1 keV in the laboratory frame. One of the investigations in this chapter includes the comparison of dynamical trajectories obtained with the END-SD-FGWP method with prescribed Coulomb trajectories [34]. Below 1 keV there are just a few theoretical studies where the PSS (perturbed stationary states) method is used to separate the nuclear and electronic degrees of freedom [62]. Total Cross Sections The scattering of a proton by a hydrogen atom has already been studied by the END-SD-FGWP method [38]. In that work an Is2s2p basis set was also employed but each function was expressed in terms of six primitive Gaussian functions. Since in this 47 work we are using a smaller number of primitives we decided to compute the total cross sections at several energies and compare them with experimental and previous END work. This comparison should provide us with some idea of how the basis set affects the total cross section. We have used the following initial conditions: i) the hydrogen atom is in its ground state with the following coordinates (20.0, 0.0, 0.0) in a.u. with respect to the laboratory frame; ii) the proton is placed at the (0.0, b, 0.0) coordinates, with b (the impact parameter) varying from 0.0 to 15.0 a.u.; iii) the hydrogen atom has zero momentum; iv) the proton has initial momentum corresponding to the energy of the collision; v) the time evolution of the system is carried out until the nuclei are separated by more than 1(X) a.u. The transition probability is obtained by projecting the asymptotic state of the proton on the eigenstates of the moving proton for elastic processes or hydrogen atom for charge exchange [38]. Total cross sections of the scattering of H'*’ on H are obtained by numerically integrating [63, 64] CX) <7{E) = 27t / bP{b;E)db (4.1) 0 for several energies E (0.5, 1.5, 5.0 keV, e.g.), where b is the impact parameter and P{b] E) the transfer probability as a function of P(6; E) for a given E. The upper integration limit is chosen to be between 11.0 and 12.0 a.u.. This choice is reasonable, since the transfer probability in this interval is smaller than 10~^. For a collision energy of 0.5 keV, we obtain the results summarized in the Table 4-3, where also the total cross section is shown as calculated using different integration 48 steps (increments in the impact parameter). As we can notice from Table 4-3, the best compromise between intervals in the impact parameter and the cross section accuracy is the last row, where we have an increment of 0.10 a.u. in the oscillatory region (0.0 - 6.0 a.u.) and 0.50 a.u. in the monotonically decreasing region (6.0 - 11.0 a.u.), see Figure 4-2 for visualization. Table 4-3: Total cross section for + H collision at 500 eV using pVDZ basis set. Interval (a.u.) Increment in b (a.u.) Interval (a.u.) Increment in b (a.u.) Cross Section ([a.u.]2) 0.0 - 11.0 0.05 64.0583 O • 1 O • o 0.10 64.0583 0.0 - 11.0 0.15 64.0671 0.0 - 6.0 0.05 6.0 - 11.0 0.15 64.0589 0.0 - 6.0 0.10 6.0 - 11.0 0.50 64.0521 49 Impact Parameter (a.u.) Figure 4-2: Weighted transition probabilities for total electron transfer at 500 eV as a function of impact parameter from ENDyne. All data in atomic units. Table 4-4 summarizes the total cross section results obtained by the END-SD-FGWP method using two different sets of primitive functions in the basis sets. Table 4-4: Total transfer cross sections for H* colliding with a H (xlO"*^ cm^). Total Transfer Cross-section. Collision END-SD-FGWP Theory Experiment^*’ Energy (eV) pVDZ basis STO-6G basis 10 35.93 36.37 37.0+5.7 100 24.31 25.60 23.7+3.5 500 17.94 19.44 18.9+3.2 1000 16.55 16.78 16.3±2.9 a: Newman et ai. Ref. [65], b: Gealy and Van Zyl, Ref. [66]. 50 As we can see from Table 4-4 the basis set do not have a large effect in the calculated integral cross section for electron transfer in the + H collision. In addition, as noted in early studies the END-SD-FGWP method provides excellent results for charge transfer cross sections when compared with experimental data. Differential Cross Sections Since the integral cross section involves an integration much of the detailed behavior of a charge exchange collision is masked. Such details can be obtained from state-to- state cross section, differential cross section, and transfer probability as a function of the scattering angle and consequently should be the basis for theoretical model analysis. We present results for the transfer probability and the reduced differential cross section as a function of the scattering angle for the H'*' colliding with H at several energies (0.25, 0.41, 0.50, 0.70, 1.0 keV). The expression for the differential cross section (aj) at a given collision energy E is [63, 64] CTd = da{e) b P{b) dO, sin 9 d6 IE (4.2) where b, 6, P{b), dil = sin 0 dO d(f> are the impact parameter, the scattering angle, the transfer probability, and the solid angle, respectively. The modulus of db/d6 is used to insure that the differential cross section is always a positive quantity. In addition, since the scattering angle by definition always lies between zero and tt, we have that sin 0 > 0. It should be noted that the differential cross section depends upon the derivative and is not integrated over any variable so that more details of the collision process are revealed. 51 We should now make the connection between the classical deflection function 0 and the scattering angle d. The classical deflection function 0 is given by [63, 64] where V(r) is the central potential, E the collision energy, and the distance of closest approach a is obtained by solving the equation (4.4) From the expression for the deflection function we can see that 0 may have any value between — oo and x, and as we mentioned before 0 < 0 < tt. As a result, the two angles have to be connected by Q = TjO — 2nv where t; = ±1 and n is a positive integer (or zero). For a given 0, the numbers t) and n are chosen so that 6 lies between zero and tt. We can also show that de dd de de db ~"^db ^ db — db (4.6) and consequently. (4.7) 52 Trajectory Analysis As mentioned before time-dependent methods using prescribed trajectories are not appropriate for calculating the differential cross sections since they involve the derivative of the deflection function 0 with respect to the impact parameter b. The END-SD-FGWP method with its dynamical trajectories should be more suitable. We present a comparison of 0 versus 6 using bare nuclei and a hydrogen atom potential. The central potential in the case of bare nuclei has the following form r) = z= r~^ (4.8) r for two protons. The resulting analytical solution for the relation between impact parameter and scattering angle is the so-called Rutherford scattering formula This is certainly not a realistic interaction potential between a proton and a hydrogen atom. Instead, using the exact wave function for the ground state of a hydrogen atom we obtain the following potential for the electrostatic interaction between a proton and a hydrogen atom in the ground state [67] I/(r) = (l-|-r-^)e-2". (4.10) There is no analytical solution of the classical deflection function for this potential, so we have to solve it numerically and the results are plotted in Figure 4-3. 53 Figure 4-3: Scattering angle as a function of impact parameter using ENDyne, bare nuclei, and hydrogen atom potential for the + H system. Eneigy: 500 eV. Basis set: pVDZ. Scattering angle in degrees and impact parameter in a.u.. Full line: ENDyne, dotted line: bare nuclei, and dashed line: hydrogen atom potential. A more detailed illustration of the differences between these three trajectories can be seen in Figure 4-4. It is clear from Figure 4-4 that the errors due to prescribed trajectories can be quite large for small impact parameter. We also notice that the dynamical trajectory lies between the bare nuclear and the hydrogen atom potentials. This is due to the fact that in the H'*’ + H system charge transfer is expected so the interaction potential should be represented by a combination of the bare nuclei and the atomic potential. The contribution of each extreme potential (bare nuclei and atomic) is related to the amount of charge transfer during the collision. Consequently, the use of prescribed trajectories in calculating properties that depends upon the scattering angle becomes doubtful for energies below 1 keV and for small impact parameters. 54 Figure 4-4: Scattering angle as a function of impact parameter using ENDyne, bare nuclei, and hydrogen atom potential for + H system. Small impact parameter region. Energy: 500 eV. Basis set: pVDZ. Scattering angle in degrees and impact parameter in a.u.. Full line: ENDyne, dotted line: bare nuclei, and dashed line: hydrogen atom potential. Transfer Probability As we can see from the expression for the differential cross section, equation (4.7), the transfer probability is an important contributing factor for this quantity. We then present in this subsection END-SD-FGWP results for this probability as a funciton of the scattering angle and compare them with the experimental data. Before performing these comparisons we present a brief discussion of the experi- mental procedure and try to identify the main sources of experimental uncertainties. The experimental procedure consists of [68, 69] a mass analyzed H"'' beam issued from a single discharge ion source crossed with a thermal H beam. Ions and atoms scattered 55 at a given angle are selected and analyzed. The overall energy resolution varies from 0.8 eV up to 2 eV in the 250 eV — 2000 eV energy range. The uncertainties affecting the experimental results arise from i) the finite angular resolution AO (± 0.07° for small angles up to ± 0.2° for large angles), ii) the errors in the dissociation fraction determina- tion (H/H 2 ratio in the collision region), and iii) the accuracy of the relative calibrations of the cross sections for the different processes (elastic, charge transfer, and excitation). Because the rapid variation of the transfer probability and the differential cross section with the angle, more particles arise from the region 6 — AO than from the region 0 -f- A^. Thus the center angle is slightly larger than the one corresponding to the most prob- able scattering. However, this error never exceeds 0.1° [68, 69]. As we shall see in the comparison between the experimental data and the END-SD-FGWP results the finite angular resolution can account for most of the damping in the oscillatory behavior of the observed versus calculated transfer probability and differential cross sections. Problems with the measured direction of the incident beam with respect to the detection system can be eliminated by studying the scattering on both sides of the incident beam direction and choosing the zero index from the symmetry of the data curves. In addition to these uncertainties, the experimental results are normalized to the theoretical results [70] for the elastic differential cross section at small angles. From the discussion of the experimental uncertainties it seems clear that the measured transfer probability versus the scattering angle is less susceptible to errors than the differential cross section and they are independent of scaling factor. Thus in Figures 4-5-4-9 we present the results for the transfer probability calculated by the ENDyne 56 program and compare with the experimental data. The theoretical procedure is the same as described in the calculation of the integral cross section. The only additional information needed is the asymptotic momenta of the proton, which are used to compute the scattering angle, that is Px Px 1 / .1 1 1 \ COS0 = = . , = = (4.11) Ptot y/Px+Py y/1 -f {PylpxY where px and py are the x- and y- components of the asymptotic proton momentum. Scattering Angle (deg) Figure 4-5: Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the H* + H system. Energy: 250 eV. Basis set: pVDZ. Experimental angtilar resolution ±0.6°, ’+’ ENDyne, from Ref. [68], and ’x’ from Ref. [69]. 57 Figure 4-6: Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 410 eV. Basis set: pVDZ. Experimental angular resolution ±0.2° — ±0.6°, '+' ENDyne, firom Ref. [68], and ’x’ from Ref. [69]. Figure 4-7: Transfer probability versus the scattering angle (in degrees) for the H* + H system. Energy: 500 eV. Basis set: pVDZ. ’+’ ENDyne. 58 Figure 4-8: Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the + H system. Energy: 700 eV. Basis set: pVDZ. Experimental angular resolution ±0.02° — ±0.6°, ’+’ ENDyne and from Ref. [68]. Figure 4-9: Transfer probability versus the scattering angle (in degrees). Comparison between the experimental and ENDyne results for the H"^ + H system. Energy: 1000 eV. Basis set: pVDZ. Experimental angular resolution ±0.07° — ±0.6°, ’+’ ENDyne and from Ref. [68]. 59 Both the qualitative and quantitative features (position of maxima and minima) of the experimental data are very well reproduced by the END-SD-FGWP method The quantitative agreement can be made much better if we apply an angular resolution window to the theoretical data. This would simulate the damping in the experimental data and explain why the oscillations in the transfer probability do not extend form zero to unity. There is still some discrepancies, mainly at high energies and large angles, between our theoretical treatment and the experimental data. We believe that the extra features in the theoretical results are real and were not picked up by the experiment due to the poor angular resolution. It also should be mentioned that the semiclassical treatment [70] using a three-term molecular basis provides results in perfect agreement with the experimental data. This is puzzling since the experimental data have a damping of the oscillation due to the finite angular resolution and the theoretical work does not mention how this damping is handled. Because of this the quantitative agreement between the semiclassical method and the experimental results seems to be fortuitous or the true meaning of the theoretical data is difficult to understand. We should add that this damping is also not observed in other theoretical treatments [34]. Reduced Differential Cross Section The experimental results available for the H'*’ + H collision are often presented as ’reduced’ differential cross sections (p) which are expressed as [68] P be P{b) db (4.12) 60 The reason for this is due to the rapid variation of the differential cross sections with the scattering angle. So, in the following Figures (4-10-4-14) we present both elastic and charge transfer reduced cross sections calculated using the ENDyne program and compare them with the experimental data. 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Scattering Angle, 0 (deg) Figure 4-10: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the H"^ -h H system. Eneigy: 250 eV. Basis set: pVDZ, ENDyne transfer, ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. [68]. 61 Figure 4-11: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the -i- H system. Energy: 410 eV. Basis set: pVDZ, ENDyne transfer, ’+’ ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. [68]. Scattering Angle, 0 (deg) Figure 4-12: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the H system. Energy: 500 eV. Basis set: pVDZ, ’*’ ENDyne transfer, ’+’ ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. [68]. 62 Figure 4-13: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the + H system. Eneigy: 700 eV. Basis set: pVDZ, ENDyne transfer, ’+’ ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. [68]. Scattering Angle, 0 (deg) Figure 4-14: Reduced differential cross sections versus the scattering angle. Experimental and the ENDyne results for the H* + H system. Energy: 1000 eV. Basis set: pVDZ. ’*’ ENDyne transfer, ’+’ ENDyne elastic, ’x’ transfer and ’o’ elastic from Ref. [68]. 63 The same comments and conclusions about the transfer probability hold for the reduced differential cross sections. We should only add that the agreement between the END-SD-FGWP and experimental results for the reduced cross sections is worse than for the transfer probability. The reduced differential cross sections are explicitly dependent upon the scattering angle. The experimental data were normalized and fitted to theoretical results, which are based upon prescribed trajectories. We can then argue that the experimental results for the reduced cross sections have been downgraded in quality by fitting to unsuitable theoretical results. The H~^ + He Collision The behavior of proton-helium scattering is different from that of the proton-hydrogen collisions. The -f He system exhibits electron-electron interaction, the charge transfer is non-resonant, and the interaction potential has repulsive and attractive parts. These facts make the scattering of a proton by a helium atom an important system to study with any new time-dependent method that has demonstrated the ability to describe well the -t- H collision. Since the early 1930’s the proton-helium scattering has been studied by theoretical methods [71]. The accumulated theoretical data falls mostly in the energy range above 10 keV and are related mainly to total (integral) cross sections. Experimental and theoretical studies of differential cross sections for + He are more scarce. However, it is possible to find several experimental studies for energies down to 5 keV and scattering angles greater than 1° [72]. These experiments are probing the repulsive part of the interaction potential and no structures are then observed. Recent studies of angular differential 64 scattering at keV energies in ion-atom collisions have focused on very small angles, in general, below 1° [73]. These studies are motivated by the highly forward-peaked character of the cross sections and by the location of the classical rainbow angle within the 0-1° range (at keV energies). The fact that structures occur for very small angle has been the main reason for the delay in investigating this region, since high angular resolution is needed. From the theoretical point of view, the oscillatory patterns seen in elastic and inelastic scattering have been studied extensively using semiclassical scattering theories [63, 74, 75]. Trajectory Analysis Some of the dynamical trajectories obtained from an ENDyne calculation of a proton being scattered by a helium atom, with collision energy of 50.0 eV is presented in Figure 4-15. In general, these trajectories should be expected in any ion-atom collision non- resonant scattering. 65 Figure 4-15: Dynamical trajectories of being scattered by He. ENDyne results at energy of 50.0 eV with a pVDZ basis set. The appearance of classical rainbow scattering implies that the classical deflection function has at least one extremum. This means that there are at least two impact parameters that will yield the same scattering angle. In the case of a proton colliding with a helium atom we show in Figure 4-16 an example of three dynamical trajectories that yield the same scattering angle. 66 Figure 4-16: Dynamical trajectories for scattering angle of 2.65°. ENDyne results for + He system at 50.0 eV with a pVDZ basis set. Impact parameters: solid line ( — ) =1.12 a.u., dashed line ( — ) = 1.60 a.u., dotted line (• • •) = 2.10 a.u. We should point out that apparently the target (He) moves only in the perpendicular direction (y-axis) to the scattering axis (x-axis). However, this is due to the scale of the graphic. In fact, the atomic target is first attracted to the incoming ionic projectile and then repelled, as can be seen in Figure 4-17 for three selected impact parameters. 67 Figure 4-17: Motion of the atomic target (He) during a collision with at 50.0 eV as computed by ENDyne with a pVDZ basis set. The initial coordinate of the He atom is (15.0, 0.0). Impact parameters: solid line ( — ) = 1.12 a.u., dashed line ( — ) = 1.60 a.u., dotted line (•••) = 2.10 a.u. The Deflection Function and the Interaction Potential In general, the deflection function of an ion-atom system has at least one extreme point. The exception is the resonant collision of a proton with a hydrogen atom. This behavior of the deflection function is due to the interaction potential which for the non- resonant case has both repulsive and attractive parts. We present a detailed analysis of the deflection function for the + He system at 50.0 keV. At this energy, the elastic scattering is the only important process. Consequently, the classical elastic scattering theory can be used to invert the deflection function to yield the intermolecular potential. This potential is compared to the dynamic potential and to the Bom-Oppenheimer potential energy surface. However, before this comparison is made, we show the deflection function in Figure 4-18. 68 Figure 4-18: Deflection function for + He. ENDyne results at 50.0 eV with a pVDZ basis set. From the expression for the classical deflection function, equation 4,3, it is possible to devise a procedure that takes monotonic parts of the deflection function and through an inversion yields the potential function. This inversion procedure consists of the following steps [76, 63]: 1. define some distance of closest approach s. 2. evaluate numerically the integral I{s) = 7T 1 I ^0 0([s^ + \/ x'^ 1 / y/y^S^ -I- 1/2 (4.13) 0 3. compute r from r = (4.14) 69 4. determine a point on the potential energy curve (r, V'(r)) from the above expression for r and from the following expression for y(r) V(r) = E{1 - e"2^W). (4.15) 5. determine values of V (r) at increasing separations r by gradually increasing the value of the lower limit s from which 0(6) is integrated. In order for this inversion procedure to be valid, the following conditions must be met. (i) Only the monotonic branches of the deflection function must used. In the case of a proton colliding with helium, the deflection function must be divided into three parts to be inverted, namely, {0, 6o), {6q, K), and {6r, oo}, and (ii) E > K(ro), that is, the potential can only be evaluated up to the classical turning point ro at the energy E, in addition (iii) only data that do not involve orbiting collisions may be inverted uniquely, that is, ,,/ X rdV(r) When this procedure is applied to the deflection function. Figure 4-18, generated by ENDyne for + He at 50.0 eV, the potential curve in Figure 4-19 is obtained. 70 Figure 4-19: Potential energy curve obtained from the inversion of the deflection function generated by ENDyne for + He at 50.0 eV with a pVDZ basis set. During an evolution calculation the ENDyne program creates a file, called restart file, that contains the history of the system evolving in time. We then use another program, called EVOLVE, to extract relevant information about the evolution. At each instant in time during the evolution the position of each particle and the potential energy of the system are recorded. As a result, we can compute the dynamical potential, that is, the potential energy at each instant in time with its corresponding interatomic distance. This dynamical potential can be obtained from data before and after the collision, and also from several impact parameters. The data from both before and after the collision of proton -i- helium for impact parameters ranging from 0.01 to 1.0 a.u., yielded the same (up to six decimal places) potential curve. One of this curves is presented in Figure 4-20. 71 Figure 4-20: Potential energy curve obtained from the direct time evolution. ENDyne results for -I- He at 50.0 eV with a pVDZ basis set. We also used an ab initio CIS (configuration interaction with single excitations) to obtain the potential energy curves for the ground state and several excited states. The results for + He are presented in Figure 4-21. 72 Figure 4-21: Potential energy curves for + He obtained from a CIS calculation with a pVDZ basis set. Attempts to plot the ground state of a CIS calculation, together with the inverted deflection function and the dynamical potential show that they all lie on top of each other. This is a very reassuring result, since consistency between three different approach has been obtained. The reason for this agreement is due to the large separation between the potential curves as can be seen in Figure 4-21. This means that the dynamics will remain on only one surface and the non-Bom-Oppenheimer effects are negligible. In addition, unlike the H'*' -h H system, the energy to localize the electron around the helium atom in the incoming asymptotic channel is very small and the initial state basically is the ground state of the (H-t-He)"^ system. 73 Classical Differential Cross Section Recent improvements in data acquisition, analysis and detection techniques, error handling, and beam collimation have enabled experimentalist to explore regions of very small scattering angle [77, 78, 73]. The angular resolution achieved by these improvements is now good enough to show the oscillatory behavior of the ion-atom differential cross sections. Oscillations in the differential cross sections can be interpreted as due to interferences between different trajectories leading to the same scattering angle. Since our treatment of the nuclear motion is classical, we should not expect the appearance of these oscillations in the classical differential cross sections. The only structure appearing in the classical treatment is due to the oscillations in the probabilities for each process. In this subsection, we present the theoretical results obtained by the END-SD-FGWG method. The semiclassical corrections to our classical dynamics and comparisons with experimental cross sections are presented in the next subsection. We present the classical differential cross section for H'*’ + He at 500 eV in some detail, and just a summary of the results for 1500 and 5000 eV. Since the differential cross section depends upon the deflection function, we start by showing the three regions of the deflection function in Figure 4-22. The regions are, I) {0, 6o } where 6 q corresponds to the value of the impact parameter where the deflection function vanishes, the so-called glory angle, 0(6 q), II) ( 6o , 6r } where W defines the rainbow angle, that is, the extremum of the deflection function 74 and, in) {br,oo} the attractive region. (4.17) Figure 4-22: Deflection function for + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: I) diamonds ’o’, II) triangles ’A’, and. III) open circles ’o’. As stated before, from the shape of the deflection function in Figure 4-22, we can see that more than one impact parameter leads to the same scattering angle 0. Thus, the expression for the classical differential cross section becomes t bi P{bi) sin 0 dQ I dh l6=6i (4.18) 75 i.e., a sum of contributions from these different trajectories. The contributions of each one of the three regions for the elastic and for the charge transfer differential cross sections are shown in detail in Figures 4-23-4-24. Figure 4-23: Elastic differential cross section for H"" + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: I) diamonds ’o’, II) triangles ’A’, and, III) open circles ’o’. 76 Figure 4-24: Charge transfer differential cross section for + He. ENDyne results at 500 eV with a pVDZ basis set. Regions: 1) diamonds ’o’, II) triangles ’A’, and, III) open circles ’o’. The elastic and charge transfer cross sections for + He at 500 eV, 1500 eV, and 5000 eV are shown in Figures 4-25-4-27. 77 Figure 4-25: Elastic and charge transfer differential cross sections for + He. ENDyne results at 500 eV with a pVDZ basis set. Figure 4-26: Elastic and charge transfer differential cross sections for H* + He. ENDyne results at 1500 eV with a pVDZ basis set. 78 Figure 4-27: Elastic and charge transfer differential cross sections for + He. ENDyne results at 5000 eV with a pVDZ basis set. As we can see from Figures 4-25—4-27 the classical differential cross section becomes singular at the rainbow angle, and tends to infinity when the scattering angle goes to zero. We also should notice the difference in orders of magnitude between the elastic and the charge transfer processes. Semiclassical Elastic Differential Cross Section As mentioned before, the resolution of the experimental differential cross sections is high enough to allow us to observe quantum effects. Consequently, the results presented for the classical differential cross sections cannot directly be compared with the experimental data, except for the value of the rainbow angle. In order to remedy this problem we use semiclassical corrections [63, 75, 74, 79, 80] to our classical elastic differential cross sections. 79 We start by reviewing one of the old stablished techniques to solve the elastic scattering problem using time-independent methods . The elastic scattering by a spherical potential V{r) can be described by the wave function ^(f) which is an eigenfunction of the Schrodinger equation _ ^y(r) -f ^{f) = 0 (4.19) with p = hk, k^ = (4.20) where p is the momentum. The elastic differential cross section can be expressed as follows ^ = \Am where A{0) is called the scattering amplitude. (4.21) In order to obtain the scattering amplitude A{0) we expand the solution of the SchrOdinger equation for the scattering problem in terms of Legendre polynomials . OO ^(^ = - Citpi{r)Pi{cos 6) (4.22) /=o where C/ are the amplitudes of the partial waves ipi{r), Pi{cos0) are the Legendre polynomials, and / is the angular momentum quantum number. This is the so called partial wave expansion, and it yields the following set of radial equations for V’/(r) - ?^V(r) - h H V’/(r) = 0 . (4.23) Using the normalization condition for ^'(f) and the orthogonality of the Legendre polynomials, we find that .. OO ^ /=0 (4.24) 80 where the phase shift is defined as 6i = lim {Phase tpi{r) — Phase ‘4’i{r,V{r) = 0)} (4.25) Consequently, the crux of quantum mechanical scattering the calculation of the phase shift. Although the expression for the scattering amplitude, equation (4.24), is exact, the series converge extremely slowly. Several approximate techniques have been developed that give analytical expressions for the scattering amplitude. For instance, employing the WKB semiclassical approximation, the following ex- pression for the phase shift can be obtained [75], OC (4.26) where (4.27) and (4.28) Differentiating the WKB phase shift 6i with respect to / we obtain OO (4.29) To which when compared with the classical deflection function (a = tq) OO OO 2J / drr-^FI~^^^{r) (4.30) 81 suggests that (4.31) This shows that the semiclassical phase shift is connected to the classical deflection function of the same angular momentum. Now that we have an analytical expression for the semiclassical phase shift we can use approximations of the Legendre polynomials Pi{cos0) to obtain expressions for the scattering amplitude. We analyze three cases, namely, a) angles not too small and not close to the rainbow angle Or, b) very small angles, and c) angles very near to Or. (a) {kb)~^ < 0 < TT {4 and 0 / Or- This range covers scattering at most angles, excluding very small angles and near the rainbow angle Or. In this range, the Legendre polynomials P/(cos 0) are approximated by 2 cos [(/ -t- 1/2)0 — 7r/4] + 0(/ ^) (4.32) and the WKB scattering amplitude can be simplified to (4.33) 0 with (4.34) Because of the small value of h these integrals may be evaluated by the method of stationary phase. After some algebraic manipulations, the semiclassical (sc) elastic 82 differential cross section is given by -SC = da 2 cos(a, - oy) + E(‘'5')i (4.35) • • X<] where (4)i = - bi sin^ d© db (4.36) b=bi is the classical cross sections originating from each impact parameter that gives rise to the same scattering angle. Assuming that -tt < 0 < x, the phase difference can be expressed as Aq^ = a,- — Qj bi = k J db 0 ( 6 ) - ke{r]ibi - rjjbj) + ^ [( 77 , • - 27 /^;) - {T)j - 2 t)s.)] where rjs is the algebraic sign of the slope of the deflection function (4.37) (4.38) and , + 1 if 0(6,) > 0 ^ . (4.39) ^ - 1 if 0(6,) < 0 We find, therefore, that unless the deflection function 0 is a monotonic function of the impact parameter 6 and is confined to values between — x < 0 < x (so that for a given 6 there is but one value of 6), we do not get the classical result for the cross section, no matter how short the wavelength. Instead there are additional interference terms. However, for truly macroscopic scattering the beam of projectiles is not exactly monoenergertic. Thus, since the phase angles a,- contain the microscopic variable k, or 83 the wavelength, a spread in the energy, gives averaged phase angles yielding the purely classical result. (b) Very small angles: 9 < {kb)~^ ~ 0 and 9 ^ 9^. This range of scattering angle is called forward glory scattering. One of the classical trajectories which contributes to the forward glory corresponds to a large impact parameter (6i in Figure 4-18). The other is the glory trajectory at small b, corresponding the coalescence of 62 and 63 trajectories in Figure 4-18. The asymptotic expression for the Legendre polynomials Pi{cos9) at very small angles is 0 As a result, the scattering amplitude for very small angles can be expressed as P/(cos 9) ~ (cos 9yj() [(/ -|- 1 /2) sin 9\ ~ (cos 9)^Jq [(/ -|- 1 /2)0] (4.40) where the Bessel’s integral is (4.41) (4.42) We expand about the glory value, (0(6o) = 0, e{i) = e{io) + -^ , , (/-/o) + ... ai /=/o (4.43) and obtain the scattering amplitude db 6=60 (4.44) 84 TTius, the elastic differential cross section for angles close to the forward glory angle is given by de dh b=bo [Jo(^^osin 0)]^. (4.45) (c) Near the rainbow angle: B k Or- The deflection function has an extremum at 0r, so we expand 0 about it 0 = 0 r + 1 SQ db^ b=br so that. c c 1 ^ / T T ^ 1 d^0 h ^ h , + 56,(7 - 7,) + b=br (J-Jr)' The scattering amplitude close to the rainbow angle is given by A{B) ~ , 1 2irkb sin $ nl/2 ,«[a+»/fc6r(Sr-^)] Ai(.) where 7 T Q = 26 b - kbrQr + TJT 4 and the Airy function is defined as 00 Ai( 2 ) = — J ducos{zu + 0 with the argument given by 2 = 2^l^k^l^1)rn{0T - B) <Pe -1/3 b=br Vt = sgn <i 62 b=br (4.46) (4.47) (4.48) (4.49) (4.50) (4.51) 85 The value of the elastic differential cross section near the rainbow angle can then be written in the form 25/Vibi/36r sin 6 (Pe d9 - 2/3 |Ai(2))^ (4.52) The expression in equation (4.52) shows the complex structure of the elastic differential cross section in the region of Br when compared to the classical differential cross section. That is, the classical singularity at ^ has been smoothed out and replaced by an exponential damping 2tX B > Or and by an oscillatory feature at ^ This oscillating behavior arises from the Airy function for ^ < 0 and the maxima in the differential cross section, called rainbow maxima, correspond to the extrema in Ai( 2 ), that is, 2 = -1.0188, -3.2482, -4.8201, -6.1633, etc. Since, in general, 0(6 « 6r) < 0 and d‘^Q/db'^\i,=i^ > 0 we have that 2 = -{Br - B) 1 d'^Q 2P d9 - 1/3 h=hr = {B- Br)q-^l^ (4.53) 1 <PQ 9 = 2it2 dip b=br Thus, the main peak occurs below Br, where 2 = -1.0188, i.e. Omax =Br- I.OISS?^/^ (4 54) Subsidiary supernumerary rainbows occur at lower angles (2 = -3.2482, -4.8201, -6.1633, etc). Experimental and Theoretical Elastic Differential Cross Sections The value of the rainbow angle, Br, can be determined experimentally as the position of the inflection in Od{B) beyond the rainbow maximum. For the scattering of a proton 86 with a helium atom the experimental determination of the rainbow angle is compared with the theoretical result in Table 4-5. The agreement between the ENDyne results and the experimental rainbow angles is excellent. Table 4-5: Experimental and theoretical rainbow angles for + He system. Energy (eV) Rainbow angle (degrees) Impact parameter (a.u.) ENDyne Experimental® ENDyne 50.0 — 2.963 1.826 500 0.32 0.3015 1.778 1500 0.11 0.1013 1.772 5000 0.03 0.0302 1.772 a: Ref. [73]. We are working in the angular range of 0.01° — 1.0° and 0.04° < {kb)~^ < 0.29° for H'*’ + He at 500 eV. Consequently, we are in the regime of very small angles, ^ ~ 0 and should use equation (4.45) for the differential cross section far from the rainbow angle. The semiclassical corrections to the classical elastic differential cross section obtained by the END-SD-FGWP method is shown in Figure 4-28. 87 Figure 4-28: Elastic differential cross section for + He. Energy: 5000 eV. Basis set: pVDZ. Solid line corresponds to the semiclassical corrections to the ENDyne results. Open circles are experimental data [73]. We have presented the results for 500 eV only since at this energy the elastic process is several orders of magnitude larger that the charge transfer. Thus, the semiclassical corrections to the classical elastic differential cross section should be good, as can be seen from the comparison between the theoretical and the experimental results. It remains to develop an appropriate formalism within END that yield these quantum interferences for inelastic processes. The H 2 Collision We have seen so far, that the END formalism based upon a single determinant and classical nuclei gives very good results for resonant and non-resonant ion-atom scattering. We now examine the performance of this method for ion-molecule collisions, in particular + H 2 . In these type of collisions not only do the electrons play an important role 88 (transfer and excitation) but so do the nuclei, since vibrational and rotational excitations take place. There is a surprising lack (virtually none) of rigorous theoretical investigations on dynamics which occur in ion-molecule collisions in low- to intermediate-energy regime (50 eW < E < 50 keV). There are several reason for this lack of theoretical treatment, but mainly that, it is quite a complex problem to obtain reasonably accurate adiabatic potentials and eigenstates as functions of intemuclear distances and molecular orientations. In addition, for an ion-molecule system, the number of internal degrees of freedom that need a proper dynamical treatment increases dramatically. The END formalism in this context gains a lot of importance since it does not need adiabatic potentials and treats properly all degrees of freedom. In the present formulation, namely, END-SD-FGWP, we have not incorporated the ETF’s (electron translation factors), so it is limited to low-energies (E < I keV) applications. The END formalism is also capable of treating dynamics at very low-energies (E < 10 eV), which will be presented in Chapter 5. In this subsection we present the collision of a proton with a hydrogen molecule in its electronic and vibrational ground state at 500 eV. We examine some of the extreme molecular orientations, namely, (0°, 0°) and (90°, 0°). This should be representative of the effects of orientation of the target molecule on the charge-transfer mechanism, cross sections (charge transfer), and energy transfer mechanisms. Trajectory Analysis The nuclear motion of an ion-molecule collision is quite a bit more complicated (and more interesting) than the ion-atom scattering showed in the last two subsections. In an 89 ion-molecule collision one would expect that vibrational and rotational excitations occur. These excitations are shown graphically via the evolution of the nuclear coordinates in time. Figure 4-29 shows the head-on collision (0°, 0° orientation) for three impact parameters that lead to the same scattering angle (0.29°). We should clarify the notation for the molecular orientation (a, /3). The first angle, a, is the angle between the scattering direction and the molecule axis, and is the dihedral angle between the scattering plane and the molecular plane. Figure 4-29 represents the nuclear motions in the laboratory frame where the proton trajectories are the lines (almost straight lines) and the motion of the hydrogen atoms in the H 2 molecule are the oscillatory lines. A more detail representation of these motions, namely, the target atom coordinates, are plotted in Figure 4-30. Figure 4-29: Dynamical trajectories for + H 2 collision from ENDyne calculation for (0°, 0°) orientation at 500 eV with a pVDZ basis set. 90 Figure 4-30: Dynamical trajectories of the taiget H 2 . ENDyne results for (0°, 0°) orientation at 500 eV with a pVDZ basis set. The nuclear motion in the laboratory frame for the (90°, 0°) orientation is plotted in Figure 4-31. These trajectories lead to the same scattering angle (0.29°) as in the case of (0°, 0°) orientation. As we can see, the trajectories for both molecular orientations are quite distinct, but a closer look shows the same qualitative behavior, i.e., at small impact parameters the target and projectile initially attract each other then the repulsive interaction overcomes this attraction and they move apart from each other. At large impact parameters they attract each other and move close together, but the target cannot keep up with the projectile velocity. In both cases the H 2 molecule starts vibrating and rotating after the collision. We should notice that even for the longest evolution times, namely, 1200 a.u. (= 30 fs) for (0.0°, 0.0°) and 800 a.u. (= 20 fs) for (90.0°, 0.0°), considered here, the rotational motion is almost absent. 91 X-Coordinate (a.u.) Figure 4-31: Dynamical trajectories for H"^ + H 2 collision from ENDyne calculation for (90°, 0°) orientation at 500 eV with a pVDZ basis set. X'Coordinate (a.u.) Figure 4-32: Dynamical trajectories of the target H 2 . ENDyne results for (90°, 0°) orientation at 5(X) eV with a pVDZ basis set. 92 Vibrational Analysis The trajectory analysis seems to show that the H 2 molecular target vibrates in quite different vibrational states (frequency and amplitude). We emphasize again that the trajectories shown in the last subsection lead to the same scattering angle (0.29°). The angular dependence of the vibrational populations p„ of H 2 molecules after a collision with H''^ have been extensively studied by theoretical and experimental methods for low- (E < 30 eV) and high-energies (E > 10 keV). However, due to recent experimental developments vibrationally resolved measurements [81], differential in scattering angle 6, have been made for the intermediate-energy range (100 eV < < 10 keV) and started challenging the theoreticians. Since we have a classical description of the nuclear motion, some approximate vibrational wavefunction would have to be determined to yield p„ as a function of 6. We have not fully analyzed this problem yet, so we present in this section the vibrational analysis of the classical motion of the nuclei. We find that the best way to obtain vibrational frequencies of the H 2 target is to use the Prony method [82] to analyze the interatomic distance of the H 2 molecule as function of time. Figure 4-33 shows the interatomic distances for all six impact parameters that lead to the same scattering angle (0.29°) for both (0°, 0°) and (90°, 0°) molecular orientations. 93 Figure 4-33: Interatomic distances of the target H 2 as function of time. ENDyne results for (0°, 0°) and (90°, 0°) orientations at 5(X) eV with a pVDZ basis set. (1 a.u. = 0.024195 fs). The Prony method is an alternative to the DFT (discrete Fourier transform) for spectral analysis and is widely used in signal processing. It is a technique for modeling sampled data as linear combination of (complex) exponentials with damping coefficients e and angular frequencies uj V i(n) = (4.55) i=\ where the number of terms p is called the order of the method and the equally spaced steps of the data sample is n{6t). Details of the algorithm and its implementation can be found in reference [82]. Once the parameters of the exponential model are determined, the frequency spectrum is obtained directly, using the infinite domain analytic Fourier 94 transform N-l x{k) = N{6t) n=0 (4.56) where N is the number of points in the data sample. In order to enhance the signal it is useful to use the relative spectral densities (RSD) as a function of the frequency u to obtain the spectrum. The RSD, expressed in decibel (dB), is defined as RSDHldbl = 101og.„(J^). (4.57) The vibration spectrum of H 2 for a specific trajectory (b = 1.25 a.u. and (a, = (0°, 0°)) is shown in Figure 4-34. The other trajectories present similar spectra. The largest order p of the method before linear dependency occurs was around 40 for all trajectories analyzed. Figure 4-34; Vibrational spectrum of the target H 2 . ENDyne results for (0°, 0°) orientation with impact parameter of 1.25 a.u. Eneigy: 500 eV. Basis set: pVDZ. 95 In Table 4-6 we present the vibrational frequencies for all trajectories shown in Figure 4-33 (impact parameters leading to the same scattering angle (0.29°) for both (0°, 0°) and (90°, 0°) molecular orientations). The differences in the vibrational states are basically due to differences in the energy transfer and the electronic state of the H 2 target after the collision. Table 4-6: Vibrational frequencies (cm *) as function of the impact parameter (a.u.) and molecular orientation (a°, 13 °). b (a.u.) 0 0 Vibrational Frequency (cm'^) 1.25 0,0 4441 1.59 0,0 4467 3.42 0,0 4612 1.70 90,0 4415 2.20 90,0 4308 3.50 90,0 4580 Transfer Probability Since the vibrational ftequencies are different for a given scattering angle, it suggests that the electronic states of the H 2 target are also different. The transfer probability reflects these differences and is plotted in Figure 4-35. 96 Impact Parameter (a.u.) Figure 4-35: Probability for charge transfer in + H 2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Eneigy: 500 eV. Basis set: pVDZ. It is interesting to note that even in these two extreme molecular orientations the transfer probability has similar qualitative behavior, however the quantitative differences are large enough to yield a variety of vibrational frequencies for the same scattering angle Deflection Function As we have seen before, the deflection function is an indicator of the appearance of the intermolecular potential. The adiabatic potential curves as a function of the - H 2 (center of mass) distance are highly dependent upon the molecular orientation for intermolecular distances smaller that 2.5 a.u. [83]. We should expect then some signif- icant differences between the deflection function for (0°, 0°) and (90°, 0°) orientations at small impact parameters. These deflection functions are shown in Figure 4-36. They 97 do confirm the findings of ab initio calculations for the dependence of the potential on the molecular orientation. Figure 4-36: Deflection functions for + H 2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Energy: 500 eV. Basis set: pVDZ. Differential Ooss Sections Since we have the probabilities and the deflection function we can compute the classical differential cross section for these two molecular orientation. The calculation of the differential cross sections is also motivated by the experimental measurements [84] of the transfer and elastic processes in + H 2 collision. Figures 4-37—4-38 show the elastic and charge transfer differential cross sections calculated by ENDyne for the two molecular orientations considered here and compared with the experimental data [84]. 98 Figure 4-37: Elastic differential cross section for + H 2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Experimental data from reference [84]. Energy: 500 eV. Basis set: pVDZ. Figure 4-38: Charge transfer differential cross section for + H 2 collision. ENDyne results for (0°, 0°) and (90°, 0°) orientations. Experimental data from reference [84]. Energy: 500 eV. Basis set: pVDZ. 99 The calculated classical elastic differential cross sections for both molecular orien- tations agree quite well with the experimental results. From the inflection in the exper- imental data for the elastic cross section we estimate that the rainbow angle should lie between 4.0° and 4.7°, which agrees quite well with the ’averaged’ theoretical rainbow angle from the molecular orientations considered here. The agreement with the experi- mental data for the charge transfer cross section is not as good as for the elastic one. The reason perhaps is that we are comparing classical cross sections with the experimental results. The latter, by nature, includes the quantum effects due to interference between different trajectories leading to the same scattering angle. In the case of a molecular target, these effects should be even greater than for atomic target, since we have also interference between trajectories from different molecular orientations. Summary and Conclusions Both quantum and classical dynamical quantities calculated by the END-SD-FGWP method is very satisfactory. Some of the quantum dynamical quantities like transfer probability and interparticle interactions are compared with experimental results and it shows an excellent agreement. So do some classical quantities like rainbow angles. However, some work has to be done to develop appropriate techniques and formalisms to correct classical results like vibrational populations and quantum effects in differential cross sections. The systems treated here + H, He, and H 2 present a wide variety of situations. For instance, we have resonant (H'*' + H), near-resonant (H"^ - 1 - H 2 ), and non-resonant (H'*’ 100 + He) situations, and the END-SD-FGWP method represented them quite well. Also, we have shown that the molecular orientation plays an important role in charge transfer processes as do the quantum effects due to contributions from different trajectories. CHAPTER 5 INTRAMOLECULAR CHARGE TRANSFER DYNAMICS The dynamics of inter- and intramolecular charge transfer has been a very important subject in chemical-physics sciences [85-89]. The main reason is that charge transfer is the essence of some most fundamental biological and biochemical processes. Also, charge transfer is important in many fields such as atomic physics, plasma physics, astrophysics, semiconductor physics, organic and inorganic chemistry, and many others. In this chapter we study the following charge transfer systems Li — H-Li'*' ^ Li"^- H — Li and Li-CN-Li"^ ^ Li'^-CN-Li. These molecular systems represent a good prototype for theoretical studies of charge transfer dynamics. They are small enough to be treated at some high level ab initio time-dependent and time-independent methods. In addition, the Li-CN-Li system can be seen as exhibiting charge transfer between two mixed valence metal atoms mediated by a bridge. Electron Transfer Formalisms Due to its conceptual simplicity the formalism for electron transfer proposed by Marcus [90-92, 89] in the late 50’s has drawn a lot of attention from both theoreticians and experimentalists. This theory has been extensively reviewed, revised, and extended. However, two basic assumptions are still being employed, namely, a) that there is a reaction coordinate that takes the reactants to the products via nuclear motion, and b) there is a coupling (// 12 ) between two electronic states (donor and acceptor). A typical electron transfer reaction coordinate is shown in Figure 5-39. 101 102 Figure 5-39: Diabatic and adiabatic potential curves for normal electron transfer. Non-dynamical type approaches for electron transfer assume that the transfer rate k is usually given by an expression based on the Fermi golden rule [92, 89, 93] 2ir »: = y|/ri2p(FC). (5.1) The quantity (FC) is the Franck-Condon factor and is related to the vibrational spectrum of the donor-acceptor system and its surroundings, if any. Analytical expressions for (FC) derived from classical and semiclassical theories exist and are discussed in references [92, 89]. The matrix element H \2 represents the electronic coupling of the donor and acceptor (or reactant and product). This is the quantity that has been challenging theoreticians and experimentalists. There are several approaches to compute the electron transfer matrix element. Usually, for large systems such as proteins, the transfer dynamics is treated as a one-electron problem [87, 93]. However, for small molecular systems many-electrons 103 ab initio and semiempirical theories have been used [89]. In this case, the most common approach to calculate H \2 is to compute the diabatic states which, in general, are non- orthogonal and use them as zeroth order states in a perturbational approach [94, 89] or use them in a two-dimensional non-orthogonal Cl (configuration interaction) problem [95, 96]. More specifically, the broken space symmetry wave functions, |1) and |2), are used in the following expression for the coupling element [97, 94, 96, 89] 1 1 Si2 {\\H\2)-Sn (l|/f|l)-H(2|i/l2) (5.2) 5i2 = (1|2). Since the nuclear motion is responsible for electron transfer [97, 89], some attempts have been made to treat the motion of the nuclei explicitly. Most of these dynamical approaches are also two-state models with a coupling H\ 2 . One of these approaches uses an ab initio method to calculate the coupling H 12 at each nuclear configuration as the system is evolving in time [98, 99]. Some simpler approaches use a model potential function for the electronic state and the nuclear motion is treated explicitly by the time- dependent variational principle [100] or by the semiclassical approach of Nikitin [101]. All these approaches for electron transfer need to compute the coupling matrix element. In order to obtain H 12 it is necessary to know first what is the reaction coordinate, and that might be very difficult to determine for large systems. We use the ENDyne program based on the END-SD-FGWP formalism to study intramolecular electron transfer. The main advantage is that we do not need to compute 104 the coupling matrix element, and to initialize the electron transfer we just distort the molecule and let the system evolve in time. Structure, Energetics, and Electron Density of LiCNLi Before presenting the results for LiCNLi we should specify that all results presented in this subsection were obtained using the 3-21G basis set [102]. In addition, we perform time-independent ab initio self-consistent field (SCF) calculations to determine potential surfaces and electronic densities of the LiCNLi and LiHLi molecules. We do not describe the details of the Hartree-Fock method and the methods that introduce electronic correlation, since we use standard approaches widely tested and reviewed in the literature [103-105]. The potential energy surface of the LiCNLi molecule has been recently studied where the basis set and electronic correlation effects were studied in some detail [106]. We augmented those results by adding data from a UHF/3-21G calculation. The linear structure of LiCNLi has two minima at all levels of theory explored. Table 5-7 presents the results for these two linear structures at several levels of theory. As we can see the general geometrical structure of linear LiCNLi is quite independent of basis set and electronic correlation corrections. 105 Table 5-7: Minima for the linear structure of LiCNLi at various levels of theory. Bond distance in pm. Theory Level Structure P Structure IF 1 LiC CN NLi LiC CN NLi UHF/3-21G 198.0 115.0 192.6 214.1 115.8 178.8 UHF/6-31G* 198.3 114.5 198.8 216.6 115.3 182.1 UHF/6-31+G* 196.5 114.6 197.9 215.9 115.4 180.7 UMP2/6-31G* 196.0 118.4 200.1 213.1 118.8 183.6 CASS/3-21G 197.0 117.4 194.3 212.0 117.7 180.5 a: All results from Ref. [106] except UHF/3-21G. These two linear structures are local minima in the potential energy surface. The proposed global minimum has a C 2 v structure at UHF/3-21G level with the lithium atoms bonded to the nitrogen. Another interesting fact about the potential energy surface of this system is that the transition state between the linear structures (I and II) has Cg symmetry. The structure of this transition state and global minimum at the UHF/3-21G level is presented in Figure 5-40. Transition Slate I-D Figure 5-40: Transition state for structures I and II and global minimum structures. Bond distances in pm and angles in degrees. The energetics of the extrema of the potential energy surface of LiCNLi is indicated 106 in Table 5-8. As we can see, the qualitative and quantitative features of the potential energy surface of LiCNLi are quite basis set independent, and only the relative stability of the linear structures depends upon the theoretical level of electronic correlation. Table 5-8; Relative energies of minima I and II and of transition state I-II (Cs). Eneigies in kcal/mol. Theory Level® Structure I Structure II Transition State Global Minimum UHF/3-21G 8.68 4.94 9.11 0.0 UHF/6-31+G* 9.78 7.00 11.58 0.0” UMP2/6-31G* 7.75 9.00 9.15 4D O • o a: AU results from Ref. [106] except UHF/3-21G, b: The global minimum has a Cj structure. We think that the LiCNLi molecule is a good prototype system for charge transfer dynamics. We present Mulliken population analysis, dipole moments, and isodensity surfaces for all extrema in the potential energy surface of LiCNLi. Table 5-9 contains charges and dipole moments for these structures. As we see from Table 5-9 almost a complete electron transfer takes place when going from structure I to II and vice versa. Since the atomic charge is not an observable, we plot in Figures 5-41-5-44 the electron isodensity surfaces of these structures for both alpha and beta spins, and confirm what has been stated using ordinary Mulliken analysis. We should mention that the assignment of alpha or beta electrons is arbitrary, and we choose to have ten alpha electrons and 9 beta electrons. As a result, assuming that the spin contamination is small, the unpaired electron is alpha and this explains why the alpha electron density is asymmetric. 107 Table 5-9: MuUiken population (a.u.) and dipx)le moments (D) of structures I and II, transition state I-II and global minimum at UHF/3-21G level. Structure I Structure II Transition State Global Minimum Li(C) 0.704 -0.168 0.017 0.283 C 0.061 0.213 0.138 0.256 N -0.599 -0.761 -0.674 -0.821 Li(N) -0.167 0.716 0.518 0.283 Dipole (D) -16.74 15.68 -6.16 -0.54 alpha density beta density Figure 5-41: Alpha and beta isodensity for structure I at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A3|. Structure: Li-C = 197.98 pm (= 1.980 A = 3.741 a.u.), C-N = 115.02 pm (= 1.150 A = 2.174 a.u.), Li-N = 192.59 pm (= 1.926 A = 3.639 a.u.). SCF energy: -106.63313290 a.u. = -2901.6368 eV alpha density beta density Figure 5-42: Alpha and beta isodensity for structure II at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A3|. Structure: Li-C = 214.10 pm (= 2.141 A = 4.046 a.u.), C-N = 115.77 pm (= 1.158 A = 2.188 a.u.), Li-N = 178.79 pm (= 1.788 A = 3.379 a.u.). SCF energy: -106.63909258 a.u. = -2901.7989 eV ti cn 109 alpha density beta density Figure 5-43: Alpha and beta isodensity of the transition state I-II structure (Cs) at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: Li-C = 208.75 pm (= 2.088 A 3.945 a.u.), C-N = 115.83 pm (= 1.158 A = 2.189 a.u.), Li-N = 181.89 pm (= 1.819 A = .437 a.u.), Li-C-N = 148.49°, C-N-Li = 138.61°. SCF energy: -106.63244750 a.u. = -2901.6180eV alpha density beta density Figure 5-44: Alpha and beta isodensity of the global minimum structure (C2v) at UHF/3-21G level. Isodensity is 0.003 a.u. (1 a.u. = 6.748 e/A^). Structure: C-N = 117.46 pm (= 1.175 A = 2.220 a.u.), Li-N = 189.22 pm (= 1.892 A = 3.576 a.u.), C-N-Li = 137.12°. SCF energy: -106.64697039 a.u. = -2902.0134 eV no The choice of an isodensity value of 0.003 a.u. (1 a.u. = I/oq = 6.748 e/A^) is based on the fact that the isodensity in the range of 0.002-0.003 a.u. yields molecular dimensions that agree with kinetic theory data, so the expected shape of the molecule can be shown [107]. Figures 5-41 and 5-42 corroborate the finding of the Mulliken population analysis and the dipole moment calculations for structures I and II. We are unable to find any transition state that is linear. Consequently, the reaction coordinate for the electron transfer is non-trivial and application of a two-state type model is not obvious. As a result, we decided not to perform any dynamical calculation and instead to explore the intramolecular electron transfer in the LiHLi molecule. The behavior of the LiHLi molecular wave function and potential surface is common among the electron transfer systems and dynamical results from the END-SD-FGWP formalism can be compared with other theoretical approaches. Structure, Energetics, and Electronic Population of LiHLi Before presenting the results for LiHLi we should specify that all results presented in this subsection were obtained using the 3-21G basis set for H [102] and 3-214G for Li [108]. In addition, we perform time-independent ab initio self-consistent field (SCF) calculations to determine potential surfaces and electronic densities of the LiHLi and LiCNLi molecules. The linear structure of the open shell LiHLi molecule at the SCFAJHF level has a broken symmetry wave function which is strongly (charge) localized. This can be seen in Figure 5-45, where we plot the difference of Mulliken charges on the lithium Ill atoms versus the relative Li-H bond distances, and the potential energy as function of the reaction coordinate in Figure 5-46. Figure 5-45: Mulliken charge differences at SCFAJHF level in the Li(l)-H-Li(2) molecule as function of relative bond distances (rj - T 2 ). 112 Figure 5-46: Potential energy curves of the Li(l)-H-Li(2) molecule at SCF/UHF level. As we can see from Figures 5-45-5-46, the broken symmetry electronic wave function (lower surface) exhibits an abrupt change at the symmetric structure (ri - r 2 = 0, where ri = distance Li(l)-H and r 2 = distance H-Li(2)), This leads to an unphysical discontinuous behavior of the charge density. The wave function describing the upper surface does not show this behavior perhaps because it has less spin contamination. We should mention that the two asymmetric linear structure minima are in fact saddle points, since they have one imaginary vibrational frequency leading to the bend structure. Nevertheless, the linear structures are a good candidate for charge transfer studies since the system is small enough to allow us to perform some computational experiments within the ENDyne program. Dynamics of Electron Transfer in LiHLi We present the dynamical results for the intramolecular electron transfer in LiHLi 113 obtained using the ENDyne program within the END-SD-FGWP formalism. The procedure utilized here is the following, (1) make a distortion of the molecular geometry (contract and/or stretch bonds) relative to the linear minima structures; (2) for this non-equilibrium geometry perform a SCF/UHF calculation and use the converged molecular orbitals as the initial z parameters; (3) let the system evolve in time and analyze the molecular properties as function of time. We have performed time-dependent calculations for four different initial conditions (geometry and electronic wave function), namely, asymmetric stretch and contraction of the Li-H bonds by 5.6%, 10%, and 15%, and an almost symmetric structure (Iri - T 2 \ = 0.004 a.u.). Table 5-10 presents the structure and energies of these initial geometries for evolution. Table 5-10: MuUiken population, structure, and relative energy of the LiHLi molecule at SCF/UHF level. Structure MuUiken population Bond distances (a.u.) Relative energy (kcal/mol) Li(l) Li(2) Li(l)-H Li(2)-H Minimum 2.387 3.243 3.0904 3.4559 0.0 5.6% 2.478 3.195 2.9170 3.6491 0.97 10% 2.498 3.186 2.7814 3.8015 3.2 15% 2.522 3.178 2.6110 3.9743 8.03 SYM 2.436 3.219 3.2549 3.2586 0.87 The results of the time evolution of these structures (Table 5-10) are presented in 114 Figures 5-47-5-50. The changes of the Li-H bond distances as functions of time are plotted in Figure 5-47. Figure 5-47: Bond distances, ri = distance Li(l)-H and tz = distance H-Li(2) in a.u. as a function of time. ENDyne calculation with a H/3-21G and Li/3-21+G basis set. Figure 5-48 has the results for the Mulliken population of the lithium atoms during the evolution, and the differences in the Mulliken atomic charges are plotted in Figure 5-49. 115 Figure 5-48: Atomic Mulliken population of Li(l) and Li(2) as a function of time. ENDyne calculation with a H/3-21G and Li/3-21+G basis set. 116 Figure 5-49: Atomic Mulliken charge differences q[Li(l)]-q[Li(2)] as a function of time. ENDyne calculation with a H/3-21G and Li/3-2 1+G basis set. Comparing the atomic charges as functions of time in Figure 5-49 with the charges in terms of the reaction coordinate in Figure 5-45 we see that the dynamical treatment seems to better describe the charge transfer since it does not exhibits unphysical discontinuity like the time-independent calculation. 117 Figure 5-50: Alpha and beta spin atomic Mulliken populations of Li(l) and Li(2) in the 15% stnicture as a function of time. ENDyne calculation with a H/3-21G and Li/3-2 1+G basis set. Figure 5-50 represents a more detailed Mulliken population for the 15% structure. As we can see the beta spin populations remain quite constant in time and all the electron transfer occur between the alpha spin electronic densities. The reason for this is that, as we said before, the assignment of alpha or beta spin for the electrons is arbitrary, and as in the case of LiCNLi, we have assigned alpha spin to the unpaired electron. Summary and Conclusion In conclusion of this chapter we again point out that the application of the END-SD- FGWP method to electron transfer does not depend upon any coupling matrix element neither is an a priori knowledge of the reaction coordinate necessary. In addition the ENDyne results for intramolecular charge transfer in the LiHLi molecules are smooth in 118 both time and nuclear displacement, unlike the SCF/UHF results that show an unphysical discontinuity at the symmetric structure. CHAPTER 6 CONCLUSIONS The exploration of the capabilities of the END-SD-FGWP method show that it is successful in providing an appropriate description of electron nuclear dynamics. The main advantages of the END formalism are that it (i) does not need the knowledge of potential energy surfaces, (ii) treats appropriately the coupling between electrons and nuclei, and (iii) does not require a priori knowledge of the reaction coordinate neither the coupling matrix element for calculating electron transfer. The approximations employed, namely a single determinantal wave function for the electrons and a classical description of the nuclei seem to be excellent starting points for generating qualitative and in some cases even quantitative results in molecular dynamics. For instance, in the case of the intramolecular electron transfer in LiHLi, the Mulliken charges generated by the ab initio SCFAJHF method show an unphysical discontinuous behavior at the linear symmetric structure. The dynamical calculations, however, do not exhibit this behavior. In fact, the Mulliken charges, calculated by the ENDyne program based upon the END-SG-FGWP formalism, change with a continuous behavior. Although, these are successful results, there is still room for further developments, from both formal and computational. Inclusion of electron correlation with an multireference wave function [2] would bring the END formalism to a level where, predictive, quantitative results could be obtained for most modest size molecular system. As we have seen from the scattering results of the 119 120 + He and + H 2 systems the quantum interferences between different trajectories are important for an appropriate description of the differential cross sections. We present semiclassical corrections to this problem, for the elastic differential cross section, which seem to work well. However, for inelastic scattering the semiclassical corrections cannot be straightforwardly employed. Thus, an appropriate approach within the END-SD- FGWP formalism needs to be developed in order to take into account these quantum interference effects. The proper way of treating this problem would be the extension to include a quantum description of the nuclei. In other words, if we avoid the FGWP (classical) approximation we could solve the problem of quantum interferences in the differential cross section calculations. This quantum treatment of the nuclei would also provide the appropriate approach for processes such as nuclear tunneling, multi-channel reactions, and zero point energy effects. The ab initio implementation of the END-SD-FGWP method requires the evaluation of all integrals and integral derivatives for each step in the differential equation solver. This makes the ab initio implementation computationally demanding, and limits the applications to small and medium sized molecules. We have presented a simplification of the END-SD-FGWP method that uses the NDDO approximation with effective cores. An appropriate parametrization of the NDDO, such as provided by the AMI and PM3 methods has proved useful and accurate in describing stationary molecular properties. It is our expectation that the same quality of description will also be adequate in the case of time-dependent methods. From a computational point of view the savings will be twofold: (i) a decrease by several orders of magnitude of the number of molecular integrals, (ii) 121 the elimination of high frequency motions (core electrons) allowing the integrator to take larger time steps, and (iii) a decrease of the number of dynamical equations. REFERENCES 1. P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930). 2. E. Deumens, Y. Ohm, and B. Weiner, J. Math. Phys. 32, 1166 (1991). 3. A. D. Hammerich, R. Kosloff, and M. A. Ratner, Chem. Phys. Lett. 171, 97 (1990). 4. J. A. Pople and D. L. Beveridge, Approximate Molecular Orbital Theory, McGraw- Hill, New York, 1970. 5. M. J. Cohen, N. C. Handy, R. Hernandez, and W. H. Miller, Chem. Phys. Lett. 192, 407 (1992). 6. M. Karplus, R. Porter, and R. D. Sharma, J. Chem. Phys. 49, 3259 (1965). 7. J. D. Doll, Chem. Phys. 3, 257 (1974). 8. L. M. Raff and D. L. Thompson, The classical trajectory approach to reactive scattering, in The Theory of Chemical Reaction Dynamics, edited by M. Baer, pages 1-121, CRC Press, Boca Raton, FL., 1985. 9. T. F. George, K. T. Lee, W. C. Murphy, M. Hutchinson, and H. W. Lee, Theory of reactions at a solid surface, in The Theory of Chemical Reaction Dynamics, edited by M. Baer, pages 139-169, CRC Press, Boca Raton, FL., 1985. 10. J. C. Tully, Adv. Chem. Phys. 42, 63 (1980). 11. B. J. Kuntz, Semiempirical potential energy surfaces, in The Theory of Chemical Reaction Dynamics, edited by M. Baer, pages 71-90, CRC Press, Boca Raton, FL., 1985. 12. E. J. Heller, J. Chem. Phys. 62, 1544 (1975). 13. E. J. Heller, Chem. Phys. Lett 34, 321 (1975). 122 123 14. E. J. Heller, J. Chem. Phys. 65, 4979 (1976). 15. E. J. Heller, Acc. Chem. Res. 14, 386 (1981). 16. D. A. Micha, J. Chem. Phys. 78, 7138 (1983). 17. C. D. Stodden and D. A. Micha, Int. J. Quantum Chem.. S21, 239 (1987). 18. S.-Y. Lee and E. J. Heller, J. Chem. Phys. 76, 3035 (1982). 19. R. H. Bisseling et al., J. Chem. Phys. 87, 2760 (1987). 20. R. Kosloff, J. Phys. Chem. 92, 2087 (1988). 21. C. LeForestier et al., J. Comput. Physics 94, 59 (1991). 22. K. C. Kulander, Comput. Phys. Commun. 63, 1 (1991). 23. M. Kimura and W. R. Thorson, Phys. Rev. A 24, 1780 (1981). 24. W. Fritsch and C. D. Lin, Phys. Rep. 202, 1 (1991). 25. M. Kimura and N. F. Lane, The low-energy, heavy-particle collisions - a close- coupling treatment, in Advances in Atomic, Molecular and Optical Physics, edited by D. Bates and B. Bederson, page 79, Academic Press, New York, 1990. 26. M. Kimura, Phys. Rev. A 44, R5339 (1991). 27. M. Karplus and J. A. McCammon, Annu. Rev. Biochem. 53, 263 (1983). 28. P. A. Kollman and K. M. Merz, Acc. Chem. Res. 23, 246 (1990). 29. W. A. Eaton and A. Szabo, Chem. Phys. 158, 191 (1991). 30. M. Parrinello, First-principles molecular dynamics, in Modern Techniques in Computational Chemistry: MOTECC-90, edited by E. dementi, pages 731-743, ESCOM, Leiden, The Netherlands, 1990. 31. P. A. Bash, M. Field, and M. Karplus, J. Amer. Chem. Soc. 109, 8092 (1987). 124 32. K. R. S. Devi and S. E. Koonin, Phys. Rev. Lett. 47, 27 (1981). 33. K. C. Kulander, K. R. S. Devi, and S. E. Koonin, Phys. Rev. A 25, 2968 (1982). 34. K. Runge, D. A. Micha, and E. Q. Feng, Int. J. Quantum Chem. S24, 781 (1990). 35. D. H. Tiszauer and K. C. Kulander, Comput. Phys. Commun. 63, 351 (1991). 36. M. J. Field, J. Chem. Phys. 96, 4583 (1992). 37. B. Hartke and E. A. Carter, Chem. Phys. Lett. 189, 358 (1992). 38. E. Deumens, A. Diz, H. Taylor, and Y. Ohm, J. Chem. Phys. 96, 6820 (1992). 39. A. Diz, Electron Nuclear Dynamics: A Theoretical Treatment Using Coherent States and the Time-Dependent Variational Principle., PhD thesis. University of Florida, Gainesville, FL, 1992. 40. J. R. Klauder and B.-S. Skagerstam, Coherent States, Applications in Physics and Mathematical Physics, World Scientific, Singapore, 1985. 41. A. D. McLachlan and M. A. Ball, Rev. Mod. Phys. 36, 844 (1964). 42. P. Kramer and M. Saraceno, Geometry of the Time-Dependent Variational Principle in Quantum Mechanics, Springer, New York, 1981. 43. J. Broeckhove, L. Lathouwers, E. Kesteloot, and P. Van Leuven, Chem. Phys. Lett. 149, 547 (1988). 44. M. C. Zemer, Mol. Phys. 23, 963 (1972). 45. J. A. Pople, D. P. Santry, and G. A. Segal, J. Chem. Phys. S43, 129 (1965). 46. R. G. Parr, J. Chem. Phys. 20, 239 (1952). 47. K. R. Roby, Chem. Phys. Utt. 11, 6 (1971). 48. K. R. Roby, Chem. Phys. Lett. 12, 579 (1972). 125 49. K. F. Freed, J. Chem. Phys. 60, 1765 (1974). 50. R. Iffert and K. Jug, Theor. Chim. Acta 72, 373 (1987). 51. K. Riidenberg, J. Chem. Phys. 19, 1433 (1951). 52. P. O. Lowdin, J. Chem. Phys. 21, 374 (1953). 53. M. S. J. Dewar, E. G. Zoebisch, E. F. Healy, and J. J. P. Stewart, J. Am. Chem. Soc. 107, 3902 (1985). 54. W. Thiel, Tetrahedron 44, 7393 (1988). 55. J. J. P. Stewart, Semiempirical molecular orbital methods, in Reviews in Computa- tional Chemistry, edited by K. B. Lipkowitz and D. B. Boyd, pages 45-81, VCH Publishers, New York, 1990. 56. M. S. J. Dewar and W. Thiel, J. Am. Chem. Soc. 99, 4899 (1977). 57. P. Hoggan and D. Rinaldi, Theor. Chim. Acta 72, 47 (1987). 58. T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 59. W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. 56, 2257 (1972). 60. T. H. Dunning, Jr., J. Chem. Phys. 53, 2823 (1970). 61. R. C. Weast, Handboodk of Chemistry and Physics, CRC Press, Boca Raton, FL, 67 edition, 1986. 62. J. P. Davis and W. R. Thorson, Can. J. Phys. 56, 996 (1978). 63. G. C. Maitland, M. Rigby, E. B. Smith, and W. A. Wakeham, Intermolecular Forces. Their Origin and Determination, Clarendon Press, Oxford, 1981. 64. J. N. Murrell and S. D. Bosanac, Chem. Soc. Rev. 95, 21 (1992). 65. J. H. Newman et al., Phys. Rev. A 25, 2976 (1982). 126 66. M. W. Gealy and B. Van Zyl, Phys. Rev. A 36, 3091 (1987). 67. J. D. Jackson, Classical Electrodynamics, John Wiley, New York, 2nd edition, 1975. 68. J. C. Houver, J. Fayeton, and M. Barat, J. Phys. B: At. Mol. Phys. 7, 1358 (1974). 69. H. F. Helbig and E. Everhart, Phys. Rev. 140, 715 (1965). 70. C. Gaussorgues, C. L. Sech, F. Masnou-Seeuws, R. McCarroll, and A. Riera, J. Phys. B: At. Mol. Phys. 8, 253 (1975). 71. M. Kimura, Phys. Rev. A 31, 2158 (1985). 72. R. L. Fitzwilson and E. W. Thomas, Phys. Rev. A 6, 1054 (1972). 73. L. K. Johnson et al., Phys. Rev. A 40, 3626 (1989). 74. R. M. C. McDowell and J. P. Coleman, Introduction to the Theory of I on- Atom Collisions, North-Holland, Amsterdam, 1970. 75. R. G. Newton, Scattering Theory of Waves and Particles, Springer- Verlag, New York, 1982. 76. U. Buck, J. Chem. Phys. 54, 1923 (1971). 77. D. E. Nitz, G. S. Gao, L. K. Johnson, K. A. Smith, and R. F. Stebbings, Phys. Rev. A 35, 4541 (1987). 78. G. S. Gao, L. K. Johnson, J. H. Newman, K. A. Smith, and R. F. Stebbings, Phys. Rev. A 38, 2789 (1988). 79. K. W. Ford and J. A. Wheeler, Ann. Phys. 7, 259 (1959). 80. M. V. Berry, Proc. Phys. Soc. (London) 89, 479 (1966). 81. D. Dhuicq and C. Benoit, J. Phys. B: At. Mol. Opt. Phys. 24, 3588 (1991). 82. S. L. Marple, Jr., Digital Spectral Analysis, Prentice-Hall, Englewood Cliffs, New Jersey, 1987. 127 83. M. Kimura, Phys. Rev. A 32, 802 (1985). 84. R. S. Gao et al., Phys. Rev. A 44, 5599 (1991). 85. B. Chance et al.. Tunneling in Biological Systems, Academic Press, New York, 1979. 86. D. C. De Vault, Quantum-Mechanical Tunneling in Biological Systems, Cambridge University, New York, 2nd edition, 1984. 87. D. N. Beratan, J. N. Onuchic, J. N. Betts, B. E. Bowler, and H. B. Gray, J. Am. Chem. Soc. 112, 7915 (1990). 88. J. J. Hopfield, J. N. Onuchic, and D. N. Beratan, Science 241, 817 (1988). 89. M. D. Newton, Chem. Rev. 91, 767 (1991). 90. R. A. Marcus, J. Chem. Phys. 24, 966 (1956). 91. R. A. Marcus, J. Chem. Phys. 46, 679 (1965). 92. R. A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 265 (1985). 93. J. W. Evenson and M. Karplus, J. Chem. Phys. 96, 5272 (1992). 94. K. Ohta, G. L. Closs, K. Morokuma, and N. J. Green, J. Am. Chem. Soc. 108, 1319 (1986). 95. C. F. Jackels and E. R. Davidson, J. Chem. Phys. 64, 2908 (1976), Non-orthogonal Cl. 96. A. Broo and S. Larsson, Chem. Phys. 148, 103 (1990). 97. M. D. Newton, Int. J. Quantum Chem. S14, 363 (1980). 98. K. V. Mikkelsen, E. Dalgaard, and P. Swanstrom, J. Phys. Chem. 91, 3081 (1987). 99. K. V. Mikkelsen and M. A. Ratner, J. Phys. Chem. 93, 1759 (1989). 100. E. Deumens, Y. Ohm, and L. Lathouwers, Int. J. Quant. Chem. S21, 321 (1987). 128 101. S. Larsson, Theor. Chim. Acta 60, 111 (1981). 102. J. S. Binkley, J. A. Pople, and W. J. Hehre, J. Am. Chem. Soc. 102, 939 (1980). 103. A. Szabo and N. S. Ostlund, Modern Quantum Chemistry. Introduction to Advanced Electronic Structure Theory, McGraw-Hill, New York, 1989. 104. P. Jorgensen and J. Simons, Second Quantization-Based Methods in Quantum Chemistry, Academic Press, New York, 1981. 105. R. McWeeny, Methods of Molecular Quantum Mechanics, Academic Press, New York, 2nd edition, 1992. 106. A. E. Dorigo and P. von R. Schleyer, Chem. Phys. Lett. 186, 363 (1991). 107. R. F. W. Bader, W. H. Henmeker, and P. E. Cade, J. Chem. Phys. 46, 3341 (1967). 108. T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. v. R. Schleyer, J. Comput. Chem. 4, 294 (1983). BIOGRAPHICAL SKETCH Ricardo Longo was bom in Sao Paulo, Brazil on June 11, 1964, to Dalva and Adhemar Longo. He has a sister Renata and a brother Milton, both younger than he. As a kid, in Pen^polis, he used to mix chemicals and make experiments. A chemist was being bom. However, after reading a Portuguese translation of Gamow’s ’The Great Physicists, from Galileo to Einstein’ he was very enthusiastic about physics, so he became a physical-chemist. In August, 1982 he enrolled in chemistry at Federal University of Sao Carlos. During his chemistry major in the Department of Chemistry at Sao Carlos, he received several scholarships from the Brazilian government for scientific research and teaching. As a result, he worked in organic synthesis, electrochemistry, high vacuum pressure gauges, and quantum chemistry. He also started a physics major in the Physics Department at University of Sao Paulo, at Sao Carlos campus. He did not have time to graduate from the latter. After getting his BA in chemistry he almost gave up quantum chemistry, when he had the opportunity to perform his master degree’s research at Fundamental Chemistry Department at Federal University of Pernambuco, in Recife. There he met Michael Zemer and made up his mind to come to Gainesville to get his degree in the Quantum Theory Project. So, after getting his master’s degree on May, 1988 he was granted a scholarship from the Coordenadoria de Aperfei 9 oamento de Pessoal de Ensino Superior (CAPES) and came to Gator’s Country to start working on his doctoral degree in the Fall of 1988. He hopes to start a research group in Brazil, perhaps in a good university, and make some contributions to science and society. 129 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. N. Yngve 0hm, Chairman Professor ot 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. Michael C. Zemer 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. // I /I . . If . William Weltner, Jr. 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. 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 Doctor of Philosophy. Professor of Chemical Engineering 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. May 1993 Dean, Graduate School