UNPHYSICAL BEHAVIOUR OF THE HINDE POTENTIAL ENERGY SURFACE: IMPLICATIONS FOR THE CALCULATION OF DISSOCIATION RATE COEFFICIENTS OF H2 + H2 by Abid Afzal B.Sc., Government College University Faisalabad, 2018 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN CHEMISTRY UNIVERSITY OF NORTHERN BRITISH COLUMBIA April 2024 © Abid Afzal, 2024 UNIVERSITY OF NORTHERN BRITISH COLUMBIA PARTIAL COPYRIGHT LICENCE I hereby grant the University of Northern British Columbia the right to lend my project/thesis/dissertation to users of the library or to other libraries. Furthermore, I grant the University of Northern British Columbia library the right to make single copies only of my project/thesis/dissertation for users of the library or in response to a request from other libraries, on their behalf or one of their users. Permission for extensive copying of this project/thesis/dissertation for scholarly purposes may be granted by me or by a member of the university designated by me. It is understood that copying or publication of this thesis/dissertation for financial gain shall not be allowed without my written permission. Title of Project/Thesis/Dissertation: Unphysical Behaviour of the Hinde Potential Energy Surface: Implications for the Calculation of Dissociation Rate Coefficients of H2 + H2 Author: Abid Afzal Signature: Date: April 2024 ii Abstract Molecular hydrogen is the dominant molecular species present in the interstellar medium and has an important role in the cooling of shocks that are associated with star formation. The two mechanisms of cooling are collisional excitation followed by quadrupole emission and collisional dissociation. Modelling the role of dissociation in this cooling needs detailed information on the state specific dissociation rate coefficients. The initial goal of this research was to compare the trajectory outcomes on the Hinde potential energy surface (PES) with those on the BMKP2 PES to assess whether it is required to do extensive and more expensive calculations to determine state specific rate coefficients for dissociation of H2 + H2 to supersede those previously determined with the BMKP2 surface. A phenomenon of double dissociation was unexpectedly identified within the Hinde PES, despite the absence of sufficient energy for such an occurrence. These results prompted a comprehensive analysis of the Hinde PES, which in turn involved an exploration of the regions that exhibit unphysical behavior. This detailed examination unveiled problematic aspects of the potential energy surface. As a result of this, it has been determined that the Hinde PES is unsuitable for calculating dissociation rate coefficients for H2 + H2 . iii TABLE OF CONTENTS Abstract iii Table of Contents iv List of Tables v List of Figures vi Acknowledgements ix 1 Introduction 1.1 Potential energy surfaces . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Computational methods in molecular collision studies . . . . . . . . 1 6 8 2 Literature survey of potential energy surfaces for H2 + H2 2.1 Recent potential energy surfaces . . . . . . . . . . . . . . . . . . . . . 2.2 Calculations on potential energy surfaces . . . . . . . . . . . . . . . 11 16 19 3 Methodology 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Broad survey of initial states . . . . . . . . . . . . . . . . . . . 3.2 Selection of initial conditions . . . . . . . . . . . . . . . . . . . . . . . 3.3 Integration of the trajectory . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Analysis of trajectory outcomes . . . . . . . . . . . . . . . . . . . . . 3.4.1 Calculation of cross-sections . . . . . . . . . . . . . . . . . . . . 22 22 22 23 26 28 29 4 Results and Discussion 4.0.1 Numerical Integrator Check . . . . . . . . . . . . . . . . . . . . 4.1 Failure of the Hinde PES when R is small . . . . . . . . . . . . . . . . 4.1.1 Potential plots for the Hinde PES when R is small . . . . . . . 4.2 Assessment and Validation of PES . . . . . . . . . . . . . . . . . . . . 4.2.1 PES Paranoia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 35 37 38 50 51 5 Conclusion and Future Directions 5.1 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Concluding statement . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 58 59 iv Bibliography 60 v LIST OF TABLES 3.1 4.1 4.2 4.3 Selected vibrational v and rotational j quantum states of the H2 molecule are presented with corresponding internal energies. Values are expressed both in Hartree (Eh ) and electron volts (eV) . . . Summary of the collision outcomes for H2 + H2 molecular interactions for the lowest impact parameter range. The data reflects the number of trajectories computed for each set of vibrational (v) and rotational (j) quantum states, including reactive (nr ), reactive dissociative (nrd ), dissociative (nd ) and double dissociative (ndd ) events. The table displays the total internal energy in Eh . Results from the Hinde PES are shown in bold. . . . . . . . . . . . . . . . . . Analysis of molecular geometries and the calculated angles. This table presents a comparison of 24 permutations of a molecular arrangement evaluated using the paranoia program, showing the calculated distances (in a0 ), angles (in radians) and the potential energy values (in Eh ) obtained through both the Hinde and BMKP2 PESs for each permutation. The angles α, β, and γ are as given in Figure 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The table represents the 24 permutations of second arrangement of nuclei. For these 24 permutations, measured interatomic distances are in a0 , and angular relationships in radians. Angles α and β make the same angle of 1.57080 radian for all 24 conformations. Hinde and BMKP2 potential values are in Eh . . . . . . . . . . . . . . vi 24 32 55 56 LIST OF FIGURES 1.1 1.2 1.3 3.1 3.2 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 Pillars of creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H2 internal energy levels in mEh . . . . . . . . . . . . . . . . . . . . . Initial survey comparison between the Hinde and BMKP2 potential for the study of H2 + H2 . . . . . . . . . . . . . . . . . . . . . . . . Selected vibrational v and rotational j quantum states of the H2 molecule. The energy is in mEh and is on the vertical axis and horizontal axis represents the rotational quantum number (j). The dotted line is the minimum energy required for dissociation and the highlighted parts are the selected v and j states [1]. . . . . . . . . . . Coordinate system used by Aguado et al. for the study of two interacting H2 molecules. . . . . . . . . . . . . . . . . . . . . . . . . . . 2 5 10 24 25 The apparent double dissociation cross sections at various translational energy levels for (5,24) (v, j) state. The vertical axis is cross section in a20 and horizontal axis is the translational energy in Eh . The arrowhead represents the translational energy for H2 (0,0) and H2 (5,24) having enough to achieve double dissociation. . . . . . . . . . . . . 34 The apparent double dissociation cross sections at various translational energy levels for (12,0) (v, j) state. The vertical axis is cross section in a20 and horizontal axis is the translational energy in Eh . The arrowhead represents the translational energy for H2 (0,0) and H2 (12,0) having enough to achieve double dissociation. . . . . . . . . . . . . 34 Relative error (∆E/Etot ) against time-step . . . . . . . . . . . . . . . 35 Average deviation against time-step . . . . . . . . . . . . . . . . . . . 36 Average of the absolute deviation against time-step . . . . . . . . . 37 Comparison between the Hinde PES and the BMKP2 PES for collinear conformation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Illustration of the Hinde PES and the BMKP2 PES in terms of their collinear conformation diagonal cuts. . . . . . . . . . . . . . . . . . . 41 Comparison between the Hinde PES and the BMKP2 PES for tetrahedral conformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Illustration of the Hinde PES and the BMKP2 PES in terms of their tetrahedral conformation diagonal cuts. . . . . . . . . . . . . . . . . 45 vii 4.10 4.11 4.12 4.13 Comparison between the Hinde PES and the BMKP2 PES in terms of trapezoidal conformation. . . . . . . . . . . . . . . . . . . . . . . . Illustration of the Hinde PES and the BMKP2 PES in terms of their trapezoidal conformation diagonal cuts. . . . . . . . . . . . . . . . . Illustration of one of the 24 atomic permutations studied by the Paranoia program, showcasing interatomic distances in atomic units (a0 ). The program evaluates the potential energy surface across all permutations to ensure comprehensive testing and validation . . . This figure represents one of the twenty-four permutations studied by the Paranoia program, the two pairs of atoms reach positions corresponding to the minima of the potential, which are proximal to the depth of potential. . . . . . . . . . . . . . . . . . . . . . . . . . viii 47 48 52 53 Acknowledgements First and foremost, I extend my deepest gratitude to my supervisor, Dr. Margot Mandy, for her exceptional guidance, patience, and unwavering support throughout the duration of this research. Dr. Mandy has been more than a supervisor to me; she has been a mentor, a guiding light, and a constant source of encouragement in the face of challenges. Her readiness to assist at any moment, her insightful feedback, and her ability to inspire and motivate have been pivotal to my journey. The depth of her commitment to my academic and personal growth has left an indelible mark on my professional development, and for that, I am eternally grateful. I would also like to express my sincere appreciation to my committee members, Dr. Ian Hartley and Dr. Andy Wan, for their valuable insights, constructive criticism, and thoughtful guidance. Their expertise and feedback have significantly contributed to the quality of this thesis. A special word of thanks goes to my family, who have provided me with endless love, support, and encouragement. To my father, Muhammad Afzal, and my mother, Rashida Naseem Zar, for their sacrifices and unwavering belief in my abilities; to my sister, Zunaira Khalid, for her constant encouragement and being my sounding board; and to my brother, Zain Ali Afzal, for his humor and lightheartedness that provided much-needed relief during stressful times. Your faith in me has been my strength and motivation. My journey would not have been the same without the support of my friends, whose companionship and understanding have been a source of joy and comfort. To all of you, I owe a debt of gratitude. This thesis stands as a testament to the collective effort and support of each individual mentioned and many who have contributed silently. I am profoundly grateful for the privilege of this experience and the opportunity to grow under such guidance and support. ix Chapter 1 Introduction The interstellar medium (ISM) constitutes a dynamic and complex environment characterized by immense turbulent structures, radiative emission, and shockfronts for stellar formation and termination. Recent research has revealed an expanding list of more than 250 distinct molecules detected within the interstellar medium [1] [2]. These molecules encompass various chemical compositions and are vital indicators of the diverse chemical reactions and physical processes occurring within interstellar space. This evolving understanding underscores the interplay among chemistry, physics, and astrophysical phenomena that shape the ISM’s composition and behavior. Hydrogen is the most abundant nucleus in the universe [1], and H2 is the most abundant molecule in the universe [1] [3] [4] [5]. Hydrogen nuclei constitute around ninety percent of the universe in terms of number density or approximately seventy-five percent of the mass of universe [6]. H2 , because of its dissociation into H atoms, can act as an essential chemical intermediate in the formation of larger species [7]. Hydrogen gas in the ISM is present in the cold clouds, also called molecular clouds, which contain much of the mass of ISM. Star formation occurs within these molecular clouds, and cooling of shockfronts plays a vital role in this 1 process. Figure 1.1 illustrates a region of star formation, as evidenced in the Pillars of Creation image captured by NASA’s James Webb Space Telescope. Figure 1.1: NASA’s James Webb Space Telescope has captured a lush, highly detailed landscape of the iconic Pillars of Creation, where new stars are forming within dense clouds of gas and dust (https://www.nasa.gov/universe/nasaswebb-takes-star-filled-portrait-of-pillars-of-creation/). It has been proposed that interstellar dust grains can act as a catalyst in the formation of molecular hydrogen [3] [8] [9]. The formation of molecular hydrogen is possible because hydrogen atoms can efficiently adsorb and diffuse over these surfaces. An H-atom adsorbs to a dust grain surface either due to van der Waals forces, which is known as physisorption, or by the formation of chemical bonds, which is known as chemisorption [8] [9]. Baryonic matter is the term used to describe matter composed of baryons, which are subatomic particles such as protons 2 and neutrons. This type of matter makes up the ordinary material in the universe, including stars, planets, and living organisms, and is what we interact with in our daily lives. A significant fraction of non-stellar baryonic matter in spiral galaxies is in H2 . The formation process can be considered in a series of steps: an H-atom adsorbs on the surface of a dust grain, encounters another H-atom, diffuses across the surface, and finally, a H2 molecule desorbs [9]. The formation of the molecule releases a binding energy of up to 4.5 eV. In regions of the ISM, particularly those characterized by lower dust densities or intense ultraviolet (UV) radiation emanating from nearby stars, the conditions are suitable for gas phase H2 formation. It is motivated by the powerful radiation from nearby stars, particularly in the form of UV photons. These high energy photons can ionize hydrogen atoms (H), creating free electrons and protons. The presence of these free electrons instigates the process of gas phase H2 formation [10]. This molecular formation process can occur through radiative association and associative detachment. Radiative association involves the interaction of hydrogen atoms and free electrons. During this process, a photon is emitted as excess energy is released. This emission typically falls within the ultraviolet or visible range. Equation 1.1 describes a step where an electron is attached to a hydrogen atom, creating a negatively charged hydrogen ion and emitting a photon. H + e− → H− + hν (1.1) In associative detachment, H2 molecules can be formed when an H-atom attaches to a negatively charged hydrogen ion, effectively detaching it from its electron. Equation 1.2 illustrates the formation of an H2 molecule by combining a hydride 3 ion with a hydrogen atom and releasing a free electron. H− + H → H2 + e− (1.2) The three-body recombination of H atoms is possible, but its rate is too slow to be responsible for the substantial amount of H2 present when the densities of atoms are below 108 cm−3 [1] [10]. Although dust grains typically provide surfaces on which hydrogen atoms can adhere and efficiently transform into H2 , in regions where dust densities are relatively low or where the influence of intense UV radiation is prevalent, gas phase H2 formation takes precedence. Understanding the prevalence and efficacy of the gas phase mechanism in the early universe is crucial for modelling the formation of primordial stars. Later-generation stars rely heavily on H2 formed on dust grains. This highlights the distinct nature of primordial star formation and its significance in the broader narrative of cosmic evolution. Molecular hydrogen is vital to astrophysics [11] because it is the dominant molecular species present in the ISM [12]. The outcome of a collision involving H2 depends on the nature of the collision partner. H2 is a closed-shell molecule with its own internal structure. In addition to having few electrons H2 , being a light molecule, has widely spaced energy levels (shown in Figure 1.2) [1]. When a low energy H2 molecule collides with another H2 molecule, the interaction is often repulsive. In such cases, the targeted H2 molecule may experience rotational excitation. Conversely, when the targeted H2 molecule collides with a partner at higher internal energy levels, it primarily undergoes interconversion between rotational and vibrational forms of its internal energy [11]. It is crucial to consider the role of the transitions among the rotational and vibrational (v, j) energy states of H2 in its ground electronic state (shown in Figure 1.2) and their cross-sections and rate coefficients for inelastic, exchange, and dissociative processes. 4 Figure 1.2: The (v, j) states and corresponding energy levels. The energy is in mEh and is on the vertical axis and horizontal axis represents the rotational quantum number j. The dotted line is the minimum energy required for dissociation [1]. The cooling effect of H2 involves the two key mechanisms which are dissociation and radiative cooling [1]. Radiative cooling involves the emission of photons as H2 molecules transition from higher energy states to lower ones. However, excited states are generally short-lived, and H2 molecules tend to rapidly return to their lower energy states by releasing the excess energy in the form of infrared photons. This emission of photons carries away energy from the shocked gas, leading to a cooling effect. The photon is permanently lost rather then being reabsorbed because of the low density. Dissociative cooling becomes important in density regimes where the reverse process has a negligible probability of occurring. This leads to permanently removing the bond energy of the molecule from the bulk system. Even at low densities, dissociative cooling can rival radiative cooling in importance, and its significance grows with increasing density [13]. Dissociation produces two hydrogen atoms from a hydrogen molecule, which can subsequently collide with other molecules. 5 An H-atom is a more efficient collider than a H2 from which it originates. This efficiency in collisions can be attributed to the smaller size of atom and its ability to interact more directly and effectively with other molecules or atoms. The reservoir of mass from which new stars are formed is mostly H2 in the interstellar clouds. The formation of a star by the gravitational collapse of an interstellar cloud from a typical density of 103 H2 cm−3 and dimension of about a parsec must increase the density by approximately 21 orders of magnitude, compressing much of the mass into the size scale of stellar dimensions, ∼ 109 m (3.24 × 10−8 parsec). The importance of cooling can be seen from the expression for the Jeans mass [14]. (︃ 5 MJ = 2 × 10 Mo T 100 )︃ 3 2 1 n− 2 (1.3) Equation 1.3 represents the baryonic mass above which a self gravitating cloud is unstable to collapse, where T is the gas temperature (K), n is the number density cm−3 , Mo is cloud mass and MJ is Jeans mass. All things made of ordinary atomic matter are referred to as baryonic by astronomers, who disregard the presence of electrons (which are fermions). 1.1 Potential energy surfaces As the lightest molecule with few electrons, the interaction of H2 with other H2 molecules or H-atoms serves as an ideal candidate for the development of accurate potential energy surfaces (PES). A PES describes the relationship between the arrangement of atoms or particles and the electronic energy associated with that arrangement. The PES concept finds application in chemistry and physics, especially in the theoretical sub-branches of these subjects. Theories involving H2 are also easiest to refine due to its simplicity, high symmetry, well-understood physics, abundance of experimental data, and its status as a benchmark system 6 for testing computational methods. This can help develop techniques to construct potential energy surfaces for larger and more complicated molecules [7]. There have been many potential energy surfaces developed to study collisions between H2 and H2 [15] [16], but this research focuses on using the Hinde potential energy surface [17]. At the beginning of this work, the Hinde PES was believed to be the best available potential for H2 + H2 study. The collisional outcomes from the Hinde PES are subsequently contrasted with those arising from another established global PES called the BMKP2 PES. A ’global’ PES, such as the BMKP2 PES, encompasses a comprehensive range of molecular conformations and energy states, offering a universal representation of potential energy for a given molecular system over various geometries and conditions. This global nature allows for a broad application in simulating molecular collisions across different energy regimes and physical scenarios. An initial survey of some potential plots of the Hinde PES [17] with BMKP2 PES [18] was carried out to observe the differences between these PESs. Cuts of the BMKP2 and Hinde potential energy surface for ”T” and collinear conformations can be seen in Figure 1.3. The figure also includes a difference plot for each of the conformations. The differences in the PESs support the hypothesis that there will be differences in the dissociative behaviour. In the top row plots in Figure 1.3, a repulsive interaction is observed when H2 approaches the target molecule, which tends to suppress dissociation. The BMKP2 plots show ”bends around the corner” representing attractive interactions promoting dissociation and energy transfer. 7 1.2 Computational methods in molecular collision studies In studying molecular collisions, particularly those involving two hydrogen molecules, researchers employ various computational methods to navigate the intricate dynamics of these interactions. Quantum mechanical methods, such as Full Configuration Interaction (FCI), Density Functional Theory (DFT), and wave function methods like Hartree-Fock, are foundational in this endeavour. These methods are most effective at low to intermediate energy ranges where quantum effects are significant. They become computationally impractical at higher energy ranges or for larger systems due to the increase in computational complexity [19]. Semiclassical methods, such as the Time-Dependent Schrödinger Equation (TDSE) and Wigner Phase Space Methods, present a viable alternative by blending quantum and classical approaches. They treat some degrees of freedom classically while handling others quantum mechanically, thus capturing essential quantum features while significantly reducing computational demands. Semiclassical methods are suitable for intermediate energy ranges where some quantum effects must be considered, but a full quantum treatment is unnecessary. These methods may lose accuracy at very high energy levels or for large systems where quantum effects become less pronounced. Given these constraints, the Quasi-Classical Trajectory (QCT) method emerges as a practical choice, particularly in the context of this study. The QCT method, treats the atoms as particles from a classical perspective. It utilizes Monte Carlo sampling for setting initial conditions, accounting for total energy, particle positions, velocities, and other relevant parameters. Trajectories are propagated using classical mechanics, specifically Newton’s equations, to calculate subsequent positions and velocities. The QCT method becomes advantageous for extensive col- 8 lision simulations where complete quantum calculations are impractical, allowing for large-scale studies. Though it overlooks quantum effects like interference and tunnelling, the QCT method offers a viable approach, combining computational efficiency with meaningful insights into molecular dynamics [20]. Its strategic employment is crucial for understanding reaction dynamics comprehensively across various energy ranges when full quantum or semiclassical methods are not feasible. The significance of this study lies in its comprehensive comparison of collision outcomes derived from two distinct potential energy surfaces for H2 + H2 . Through this comparison, a clear comprehension of both commonalities and disparities between these surfaces is achieved. The disparity between the two potential energy surfaces becomes evident upon analyzing the results. Initially, the objective of this research was to determine whether additional calculations on the Hinde potential energy surface were necessary. However, the preliminary findings indicate that trajectory outcomes from the Hinde PES are unphysical. An in-depth investigation of the unviable Hinde PES was required, which led to an exploration of the reasons that render studying H2 + H2 dissociation on the Hinde PES . 9 "' "' ... ... N N ------------------------4 2 --------------- ---------- 8 4 2 R 8 R "' ... ... I\ N ------------ - - - - - - - - - I...._ - '/'( - - N ----------------2 R "' ... II II II ... N 4 8 ,,, N 4 R 8 Ill I Ill I Ill I Ill I Ill I Ill I Ill I Ill I Ill I "' I 2 4 R 8 0 0 2 R Figure 1.3: Comparison between Hinde (first row) and BMKP2 (second row) surface plots in terms of T conformation (first column) and collinear conformation (second column). The third row shows the difference between both plots. r1 (a0 ) is the internuclear separation of one molecule and R is the distance between centers of mass of two molecules. The bond length of second molecule is fixed at the inner turning point of the v = 12 state. The energy is in hartree (Eh ) [12]. 10 Chapter 2 Literature survey of potential energy surfaces for H2 + H2 Many potential energy surfaces for H2 + H2 have been developed. This study is focused on the use of the Hinde PES and comparing its results with the BMKP2 PES. Research on the interaction potential between H2 molecules began in the 1970s and has continued since then. In 1974, George Zarur and Herschel Rabitz [21] explored molecular collision phenomena where both colliding partners were molecules. They noted that most prior studies were limited to atom-diatom systems due to the complexity and computational demands of non-perturbative calculations in molecule-molecule collisions. To overcome these challenges, they developed an approach which reduced the dimensionality of the problem by assuming rigid rotors. This simplified the computational complexity and allowed the assessment of rotational and vibrational transitions. They calculated crosssections for twelve rovibrational states and found the quantum statistical effects to be significant. They developed an analytical potential by averaging over angular momentum (M). 11 The twelve rovibrational states were ordered by increasing energy (v,j) = (0,0), (0,2), (0,4), (0,6), (0,8), (1,0), (1,2), (1,4), (0,10), (1,6), (1,8), and (2,0) respectively. They calculated integral and differential cross-sections. Their research indicated that the effective potential method could be helpful in various ways, including fitting rates or cross sections and dealing with rotational and vibrational-rotational inelastic collisions. It also showed that quantum statistical effects played a significant role in differential cross-sections. The study laid the groundwork for further exploration of molecular collision behaviour and its experimental examination. In 1980, Brown and Silver extensively investigated the collision dynamics of the H2 + D2 system, employing quasiclassical dynamics to explore reactivity across different energy regimes. This research, building upon two earlier studies that used distinct potential energy surfaces and electronic energy levels [22] [23], focused on how variations in PES influenced reaction probabilities and energy transfer within this molecular system. The first study explored energy transfer and reactivity through double exchange reactions (H2 + D2 → 2HD) at a total system energy of 159 millihartree. The second part of their study utilized a valence bond model surface, incorporating a semi-empirical evaluation of integrals, including manyatom contributions. Specifically, Brown and Silver examined vibrational states in the range of v ′ < 12, uncovering that, in specific scenarios, one product molecule would be in a low vibrational state (v ′ < 2). At the same time, the other was highly excited, with 10 < v ′ ⩽ 12 [24]. Their findings highlighted common trends, such as the critical role of initial vibrational energy distributions in dictating the outcomes of collisions and the significant gains in rotational energy. In non-reactive collisions, the study noted predominant translational and vibrational energy transfers. A notable divergence was observed when comparing the outcomes from valence bond model surfaces with those from London-type surfaces, particularly in double exchange reactions. This contrast brought to light the sensitivity of collision 12 dynamics to the specific properties of the PES. Furthermore, the research of Brown and Silver raised questions about lowenergy reaction pathways on certain surfaces, emphasizing the impact of surface parameterization on accurately predicting reaction barriers. This comprehensive study substantially enhanced the understanding of collision dynamics in the H2 + D2 system, particularly under reactive conditions, contributing significantly to molecular dynamics. They also pointed to the need for an improved PES. In 1988, Schwenke et al. [25] focused on calculating recombination rate coefficients for hydrogen containing compounds in the gas phase. This research was motivated by the need for accurate rate coefficients as inputs for modeling hydrogen combustion processes. The study initially focused on the process where hydrogen atoms were combined to produce H2 while an abundance of thermal molecular hydrogen gas (H2 ) was present. This system was chosen for its simplicity and because experimental measurements were available for comparison. The study employed the resonance complex theory, a model that considers the quantum mechanical nature of meta-stable states involved in recombination, and combined it with the quasiclassical trajectory method. The relationship between resonance complex theory (RCT) and a PES lies in their combined use to study chemical reactions and processes involving resonance behaviour. RCT models the complex electronic structure and transitions between multiple electronic states in such cases. Each electronic state is associated with its distinct PES, which may intersect with the PES of other electronic states. This intersection allows for a comprehensive examination of the transitions among various electronic states. At the same time, the PES provides the necessary information about how the energy changes as the molecular geometry evolves during the reaction. The focus of their study was not on electronically excited states of molecular hydrogen but rather on the resonant metastable states, utilizing the ground electronic state of molecular 13 hydrogen. These metastable states play a pivotal role in understanding the energy transfer mechanisms within molecular hydrogen. The PES often serves as a foundation for understanding the potential energy barriers, reaction pathways, and energy changes during resonance-driven reactions. To improve the precision of the computations, a new potential energy function was constructed based on ab initio electronic structure calculations. Despite these efforts, the agreement between the theoretical results and experimental data was poor, with the theoretical values being approximately half of the recommended experimental values over a range of temperatures from 300 to 5000 K [26]. This discrepancy suggested that quantum effects and other factors needed further investigation to provide more reliable predictions for recombination rate constants in hydrogen hydrogen recombination processes. In approximately 1991, Boothroyd and colleagues conducted calculations involving ab initio energies for 6046 different conformations of H4 in the ground electronic state. They employed a multiple reference (single and) double excitation conformation interaction (MRD-CI) program for this purpose [16]. Through their analysis, both systematic and random errors were assessed, revealing a root mean square (rms) error size of 0.6 mEh . This led to a cumulative rms error of around 0.9 mEh (equivalent to 0.55 kcal/mol) in the final ab initio energy values. Notably, existing analytic potential energy surfaces (PESs) for H4 were demonstrated to accurately represent pairs of H2 molecules only when the distance between them was greater than roughly 3 a0 . The conclusion by this study highlighted the necessity for a more comprehensive analytic PES for H4 , as the existing ones were confined to low energy values, had limited practicality and were not feasible to study dissociation. Based on the data generated by Boothroyd et al. [16], Aguado and collaborators [15] fitted a potential energy surface for H4 . This surface marked the first 14 global PES to encompass the full range of conformations necessary to study of four hydrogen atoms. The Aguado PES employed functional forms that were especially convenient for quantum calculations. The resulting global H4 PES achieved through this study exhibited complete symmetry concerning permutations of hydrogen atoms. The significance of both the Boothroyd and Aguado surfaces persists in ongoing research [27], underscoring their enduring importance for global potentials for H2 + H2 . Flower and Roueff [28] employed the potential fitted by Schwenke et al. [25] to examine the rovibrational relaxation of H2 by H2 using quantal techniques. However, these techniques could not encompass all degrees of freedom due to their limitations. A comparison with previous semi-classical calculations revealed the unreliability of semi-classical methods at low temperatures. Moreover, their quantal results underestimated experimental measurements of the vibrational relaxation rate coefficient at intermediate temperatures (around 500 K). This indicated potential shortcomings in the theoretical model. To enhance accuracy, Flower and Roueff suggested further research on the H2 + H2 interaction potential, especially its dependence on vibrational coordinates. Diep and Johnson [29] undertook the task of constructing a potential energy surface (PES) for the H2 + H2 system, with a particular focus on the (H2 )2 dimer. This endeavor employed the highly accurate coupled-cluster theory with singles, doubles, and perturbational triples excitations (CCSD(T)) for PES calculations. Acknowledging the advantages and limitations of CC methods, like those employed by Diep and Johnson, is crucial. CC methods, especially the CCSD(T) variant used in their research, are renowned for their high accuracy in small to moderatelysized molecules with well-defined ground states. They offer a systematic approach to improving accuracy through higher-level excitations. However, these methods have significant computational demands, limiting their use for larger systems. A 15 notable limitation is their reliance on a single-reference ground state, rendering them less effective for highly correlated or multi-reference systems. Additionally, the Born-Oppenheimer approximation integral to CC methods may not suit studies involving non-adiabatic processes where electronic and nuclear motions are coupled. This necessitates careful consideration of the methodological choices in computational studies of molecular interactions. To model the (H2 )2 dimer, Diep and Johnson extrapolated their calculations to the complete basis set limit using a rigid monomer model, producing a PES characterized by 37 angular conformations with notable anisotropy. Specifically, they calculated the potential energy surface for a vibrationally-averaged rigid monomer bond length model of the (H2 )2 dimer, with the bond length of H2 + H2 held constant at 1.449 a0 , consistent with the rigid rotor approximation. This PES for rigid-monomer H2 + H2 interactions is described in terms of the center of mass intermolecular distance (R) and three angles (θ and ϕ). In this study, angular increments of 30 degrees were chosen for each of the three angles, resulting in 37 unique conformations after accounting for symmetry. 2.1 Recent potential energy surfaces The potential developed by Boothroyd et al. often referred to as BMKP2 [18], put emphasis on constructing a global H4 PES. This PES is of great importance in collision dynamics and serves as a valuable test case for understanding molecular interactions. The aim of the study was to attain a highly precise H4 PES, achieving chemical accuracy with precision better than one millihartree, which is vital for comprehending reaction dynamics. The approach involves extensive computations of ab initio H4 energies, encompassing both ground and excited electronic states, to improve the coverage of H4 conformational space to include regions im- 16 portant for dissociation. The resulting enhanced H4 PES offers significant advantages. It is particularly valuable for studies important to astrophysical applications, especially in examining interactions between H2 molecules in low density environments like giant molecular clouds. These regions experience heating from shockfronts, leading to various molecular excitation and collision induced dissociation processes. These phenomena are challenging to investigate through laboratory experiments due to the substantial mean free paths of molecules . Therefore, computer simulations are essential for exploring the physics and chemistry of such regions. This PES represents a substantial advancement over prior H4 potential energy surfaces. The BMKP2 PES, designed for molecular dynamics computations on a fulldimensional surface, was found to have higher anisotropy, resulting in higher values of vibrational relaxation rate coefficients, than the Hinde PES. Moreover, the BMKP2 PES was developed with analytical derivatives, unlike the Hinde PES which relies on numerical derivatives. Due to its lower computational demands, these characteristics make BMKP2 a preferred choice for quasi-classical trajectory calculations in H2 + H2 reactions. However, it is crucial to acknowledge that this computational advantage may sometimes result in discrepancies compared to experimental data [30], raising concerns about its accuracy and reliability for real world predictions. The primary objective of the study done by Hinde was to construct an accurate potential energy surface for the (H2 )2 dimer. This dimer involves two hydrogen molecules interacting, and understanding its behavior is crucial in various contexts, including molecular dynamics, astrophysics and pressure-broadening in spectroscopy. To create this PES, Hinde employed Gaussian to calculate the interaction energy between two H2 molecules (H2 –H2 interaction energy). This calculation involved a sophisticated treatment of electron-electron interactions known as 17 the coupled-cluster method with single and double excitations and a perturbative treatment of triple excitations, abbreviated as CCSD(T) [31]. To ensure the accuracy of the calculations, a restricted Hartree–Fock approach was utilized for the reference wave function, and it was confirmed that this selection remained stable for the studied H2 bond lengths. A basis set called aug-cc-pVQZ was also used for the four hydrogen atoms involved in the calculations [32]. The aug-cc-pVQZ basis set plays a pivotal role in ensuring the accuracy of the potential energy surface calculations particularly in regions where accurate modeling of electron distributions is crucial. This basis set is a comprehensive Gaussian set characterized by its ’augmented’ feature. It includes additional diffuse functions essential for accurately modelling the electron distributions in molecules, particularly in systems with non-bonding electron pairs or anions. The ’cc-pVQZ’ part stands for correlation-consistent polarized valence quadruple zeta (ζ), indicating a high-quality set with multiple (ζ) functions representing various angular momentum states and polarization functions to capture the directional nature of chemical bonds and electron interactions. While this basis set provides enhanced accuracy, especially for small to mediumsized molecules, it has the drawback of increased computational demand. This higher computational demand is a trade-off for the detailed and precise electronic structure calculations it facilitates. Nevertheless, the choice of the aug-cc-pVQZ basis set underscores the commitment of this study to precision in modelling the electronic structure of the H2 + H2 system, ensuring that the calculated PES is as accurate as possible within the constraints of computational resources. The constructed PES was then tested by comparing its predictions to experimental data, particularly regarding infrared and Raman transitions. Remarkably, the empirical adjustments made to the PES to better align with experimental results led to excellent agreement. This means that the calculated energies for transitions involving the (H2 )2 dimer closely matched what was observed in real world experiments [30]. 18 Furthermore, the study did not stop at reproducing existing experimental data, it also predicted the energies of additional transitions that had yet to be observed. This predictive capability opens up possibilities for future experiments to validate the PES accuracy or identify areas where further refinements are needed. In summary, this study significantly advanced our understanding of the (H2 )2 dimers behavior by providing a highly accurate potential energy surface. The fundamental drawback of utilizing this surface is the extra computational resources the calculations need since they employ numerical derivatives rather than the analytical derivatives used for the BMKP2 surface [1]. 2.2 Calculations on potential energy surfaces The H2 + H2 molecular system is of interest in chemical physics due to its fundamental role in understanding molecular interactions. The research by Balakrishnan [33] centred on conducting quantum dynamics calculations to analyze the transfer of rotational and vibrational energy during collisions involving two paraH2 molecules. The initial rotational state was set at j = 0. Balakrishnan spanned collision energies from the ultracold limit to thermal energies. Specifically, the research examined the vibrational relaxation of H2 (v = 1, j = 0) through collisions with H2 (v = 0, j = 0), as well as rotational excitations in collisions between ground state H2 molecules. Vibrational and rotational transition rate coefficients were determined from quantum mechanical calculations on these surfaces. These coefficients were then compared with experimental observations and calculations on the Diep and Johnson potential [29], unveiling distinct patterns and PES-independent characteristics. The BMKP2 PES exhibited higher anisotropy, which is the dependence of potential energy on the molecular orientations during collisions. The heightened 19 anisotropy observed in the BMKP2 PES is the consequence of how interaction potential was fitted. Specifically, to achieve accurate results, the BMKP2 PES necessitated the inclusion of terms up to λ1 = λ2 = 8 and 10, a requirement that underscores its significantly greater anisotropic nature ( The symbol λ represents the angular momentum quantum number used to expand the interaction potential. Specifically, in describing molecular interactions and potential energy surfaces like BMKP2, λ is used to quantify the angular dependence of the interaction potential). This significant anisotropy led to larger relaxation cross sections, nearly an order of magnitude greater than those observed experimentally at energies below 0.1 eV. Such anisotropy indicates a heightened sensitivity to molecular orientation, suggesting an overestimation of the available pathways for vibrational energy redistribution, thereby enhancing the vibrational relaxation rates. Experimental data shows that the actual molecular interactions are less orientation-dependent than suggested by the BMKP2 PES [30] [34]. In contrast, the Hinde PES demonstrates a lesser degree of anisotropy and gives results which are aligned more closely with experimental results [34]. The Hinde PES appears more accurate and reliable surface for studying the H2 + H2 collisions, especially within lower internal energy states. Consequently, the assessment of anisotropy present in a PES is a critical to the accuracy of simulated collision dynamics and their correspondence with actual experimental data. Bohr et al. [35] assessed how well the Centrifugal-Sudden (CS) approximation works when applied to collisions between H2 molecules at low energies (below 100 K). The CS approximation simplifies rotational dynamics during molecular collisions by assuming a ’frozen’ or ’sudden’ state of rotational motion. It reduces computational complexity and makes it suitable for low-energy collision scenarios. This quantum mechanical approach focuses on rotational and vibrational states, efficiently predicting outcomes in molecular collisions, especially where rotational 20 energy is significant. One of the key strengths of the CS approximation lies in its computational efficiency and effectiveness, where it provides a practical tool for studying molecular interactions with substantial rotational energy. However, its suitability diminishes in high-energy collisions and in systems with complex interactions, where rotational states may change significantly during the collision. Bohr compared the results from the CS approximation to those obtained using the highly precise Coupled Cluster (CC) method, and both methods were used on a complicated potential energy surface that accurately describes the interactions between the molecules [35] [17]. The CS method accurately predicted what happens when two molecules scatter off each other, undergo pure vibrations, or experience single-molecule rotations. Single-molecule rotation refers to the rotational motion of an individual molecule around its center of mass. This motion is distinct from the translational and vibrational movements of the molecule and is quantized, meaning it occurs at discrete energy levels. Overall, research by Bohr showed that the CS approximation is practical for studying similar molecular collisions with significant rotational energy the results of which could have applications in astrophysical models. Sultanov and Guster [36] considered rotational energy transfer for HD + H2 , H2 + H2 and H + H2 collisions in an attempt to estimate the HD cooling equation. This study was done on four potential energy surfaces, including the Diep and Johnson potential and BMKP2 potential. Thermal rate coefficients were then found from the previously mentioned molecular systems and the cooling function was developed. 21 Chapter 3 Methodology To address the computational challenges associated with the size of system, and energy range, this research uses quasi-classical trajectory (QCT) calculations. The decision to employ QCT is grounded in the impracticality of conducting quantum calculations due to the substantial increase in the number of required wave functions, a challenge that escalates with the fifth power of the total energy of the system. The computational constraints make a purely quantum methodology intractable for the objectives of this study. This study adopts methodology based on that used by Mandy et al. in 1998 [11], with differences in the potential energy surface and type of derivatives. 3.1 Overview 3.1.1 Broad survey of initial states In our preliminary examination of H2 + H2 molecular collision dynamics, a broad range of initial vibrational and rotational states were selected. This range is represented in Table 3.1 and Figure 3.1. For the purpose of this section one H2 molecule is designated as the collider with an associated kinetic energy vector, and the other 22 as target. In our exploration of H2 molecular collision dynamics, the quantum state identified by v = 0 and j= 0 delineates the ground rovibrational state of the hydrogen molecule in its ground electronic state where the molecule possesses the lowest possible vibrational and rotational energies. Three intermediate states were carefully chosen roughly halfway to the dissociation limit. One state is characterized by a near equal partition of energy between vibrational and rotational forms, another is dominated by vibrational energy, and the third state as primarily rotational energy, except for the zero-point vibrational energy inherent to the system. Finally, three states were chosen close to the dissociation limit. This selection facilitates a detailed examination of the energy distribution within the molecule effects the collisional outcomes. It enables an assessment of how the molecule responds when internal energy is variably allocated between vibration and rotation yet remains below the dissociation threshold. States that are close to the dissociation limit with rotational energy can be quasibound. These quasi-bound states, which have internal energy above the dissociation energy, occupy a special place in molecular dynamics because of dissociative tunnelling. This selection of states will allow us to survey the effects of the differences between the two PESs. 3.2 Selection of initial conditions The selection of initial conditions for our exploration of molecular collision dynamics within the H4 system are based on the coordinate system presented by Aguado et al. as illustrated in Figure 3.2. 23 v j 0 0 2 5 0 5 12 0 20 14 0 32 24 0 Eh 0.0099 0.0974 0.0909 0.0942 0.1764 0.1755 0.1662 eV 0.2693 2.6503 2.4735 2.5633 4.8000 4.7755 4.522 Table 3.1: Selected vibrational v and rotational j quantum states of the H2 molecule are presented with corresponding internal energies. Values are expressed both in Hartree (Eh ) and electron volts (eV) 0 10 20 30 Figure 3.1: Selected vibrational v and rotational j quantum states of the H2 molecule. The energy is in mEh and is on the vertical axis and horizontal axis represents the rotational quantum number (j). The dotted line is the minimum energy required for dissociation and the highlighted parts are the selected v and j states [1]. 24 ,. •.,/ , •· - -- - ···· ·· ··---· - - · - - - - A Figure 3.2: Coordinate system used by Aguado et al. for the study of two interacting H2 molecules. Three distances and three angles are selected as independent coordinates named r1 , r2 , R, α, β, and γ, r1 being the distance between atoms A and B; r2 being the distance between atoms C and D; R being the distance between the center of mass of the AB and CD atomic molecules [15]. 25 The initial conditions for each H2 molecule involves Monte Carlo selection of angular momentum based on the rotational quantum number j. Angular momentum of first molecule determines α, the angular momentum of second molecule determines β, and γ, as depicted in Figure 3.2. The initial orientation of one molecule out of this plane defines the γ angle, completing the spatial orientation for the molecular system. The vibrational phase is determined for each molecule using its effective potential for each diatomic. The vibrational phase selection uses an acceptance-rejection algorithm within the inverse quadrature framework, resulting in a precise sampling of r1 and r2 for each molecule. The distance R is the separation between the centers of mass of the two interacting molecules. Initially, R is set to a value where the molecular interaction is negligible, ensuring an accurate starting point for the simulation. The selection of the impact parameter specifies the direction of relative velocity. A stratified sampling technique is employed for its selection to optimize computational resources, defining an inner and outer limit for the impact parameter range and selecting uniformly from the annular area. Counts of inelastic events are collected for the sub-strata to determine the need to explore further annuli of of impact parameter. This systematic approach allows for a focused exploration of collision dynamics, concentrating computational efforts on regions where significant interactions are likely. 3.3 Integration of the trajectory In this section, we delve into the numerical integration of the equations governing the molecular system dynamics under study. The Hamiltonian H that governs the dynamics of the H2 + H2 system incorporating the Hinde potential can be ex- 26 pressed as: H = T + VHinde (R, r1 , r2 , θ1 , θ2 , ϕ) Here, T represents the kinetic energy term, including monomers’ translational and rotational kinetic energies. VHinde denotes the Hinde potential energy surface, which is a function of the intermolecular distance R, the intramolecular bond lengths r1 and r2 , and the angular coordinates θ1 , θ2 , and ϕ, encapsulating the full range of molecular orientations and interactions in the six-dimensional space. A set of trajectories (a group of individual simulations, each representing a possible path that molecules may follow given a set of initial conditions) commences with an initial kinetic energy set at 0.4780 Eh , reducing decrementally, first by 0.0080 Eh steps and then, as the molecules near the dissociation threshold, in intervals of 0.0020 Eh . The translational energy is decreased until no inelastic events are observed. For the integration process, a coordinate system centered around the center of mass of the interacting molecules is employed. This system adopts cartesian coordinates, explicitly accounting for the momenta coordinates of three atoms, with the fourth atom’s coordinates inferred through momentum conservation principles and the positioning of the center of mass. Using a variable step Runge-Kutta method, the equations of motion are integrated, monitoring energy conservation and momentum throughout. To gather meaningful results, a substantial number of trajectories, at least 80,000, were simulated for each translational energy and impact parameter range. This approach thoroughly explores collision behaviour and outcomes across a range of conditions. DVERK, a double-precision subroutine for solving systems of first-order differential equations, has been used to determine the trajectories. The DVERK algorithm, which offers adaptive step-size capabilities based on Runge-Kutta methods, 27 is widely used in scientific and engineering applications for its efficiency in solving ordinary differential equations. Each step requires the calling of potential energy function, which determines the forces acting on the atoms. The BMKP2 PES utilizes the same cartesian coordinates as the integration coordinates as it only requires a single call for the potential and its analytical derivatives. The Hinde PES employs the same set of coordinates as Aguado et al. [15] for molecular interactions within the H4 system. These coordinates include three intermolecular distances, which are r1 , r2 , and R and three angular coordinates, α, β, and γ that define the orientation of the molecules in space. The use of this potential is computationally demanding, as it requires numerically calculating six gradients at each integration step, leading to thirteen calls to the potential energy function. Each call requires a transformation of the coordinates from the integration coordinates to the Hinde PES coordinates and back again. Such computational demands have significant implications for the efficiency of simulations. The numerical approach of the Hinde PES, while providing the precision needed in force calculations without explicit analytical expressions, leads to increased computational resource requirements. This difference is a consideration in choosing between these two PESs, particularly in large-scale simulations where computational efficiency is paramount. 3.4 Analysis of trajectory outcomes A trajectory is considered complete when, out of the six interatomic distances in a four-atom system, four surpass a specified threshold. This threshold ensures that two pairs of atoms have achieved adequate spatial separation, such that the forces between the molecules are neglegible. 28 Once a trajectory is complete, the analysis focuses on the two smallest atomic distances. While the closest pair of atoms may form a molecule, this is not a guaranteed outcome. The other two atoms, by implication, may constitute another pair. This phase of the study involves a detailed final energy analysis. For each diatomic pair, the rotational motion about their respective centers of mass and their relative motion along the axis that connects these centers are considered. This allows the determination of an effective potential for the pair of atoms. The angular momentum of each diatomic pair around their center of mass is evaluated, resulting in continuously valued rotational quantum numbers. These rotational quantum numbers, in turn, inform the calculation of the vibrational quantum numbers by inverse quadrature on the associated effective potential. This process is repeated for the second pair of atoms if energy transfer is under consideration. 3.4.1 Calculation of cross-sections The cross-sections for energy transfer can be calculated using the bin histogram method. This technique categorizes the continuously valued QCT quantum numbers into discrete intervals or bins. It is necessary to consider the potential systematic issues with binning, such as the choice of bin size and boundaries, which can affect the accuracy and resolution of the data. Since dissociation is under consideration binning is not necessary in order to obtain counts of dissociative events. In alignment with the methodology established by Mandy [11] [12], the probabilistic assessment of collision outcomes within this study was conducted through ni , wherein i indexes the individual impact parameter strata. the equation Pi = -N i The variable ni represents the count of specific events within a stratum, while Ni denotes the total number of trajectories computed for that stratum. When a batch of trajectories has been converged with respect to impact parameter the integral cross sections, σ, can be determined by using following formula 29 ⌜ ⃓ Pi (1−Pi ) ⃓∑︂ + N12 + N13 ∑︂ Ni ⃓ i i A2i [︃(︂ σ ± dσ = Pi Ai ± ⃓ )︂ (︂ )︂2 ]︃ ⎷ 3 2 i i 1 + Ni 1 + Ni (3.1) where i denotes the specific impact parameter stratum, Pi represents the ratio of ni to Ni , signifying the probability of the desired process taking place within that particular stratum. Here, ni signifies the count of trajectories that result in the desired process, Ni corresponds to the total count of trajectories executed within stratum i, Ai denotes the area of this stratum, and σ represents the cross section. 30 Chapter 4 Results and Discussion The preliminary survey of the Hinde and BMKP2 PESs showed notable differences in potential plots. This initial observation motivated the current study, particularly since the Hinde PES was previously believed to be the best option for studying the collisions of H2 molecules. This context sets the stage for Table 4.1, which showcases a subset of the number of counts for the lowest impact parameter range from a simulation involving 80,000 trajectories. Each trajectory in this dataset represents a unique H2 collision. Surprisingly, the analysis revealed an identical number of reactive (nr ) and reactive dissociative (nr,d ) events, hinting at a direct link between reactive collisions and dissociation with the Hinde PES. Moreover, the data presents an anomalously high count of double dissociative trajectories (ndd ), a phenomenon occurring despite the apparent insufficiency of total energy for such events. This unexpected finding underscores the need to delve deeper into the discrepancies between the Hinde and BMKP2 PESs, especially in their representation of high-energy molecular interactions. 31 v1 j1 Collider 0 0 0 0 Eh nr,d nd ndd 0 0 0 0 0 0 0 2 0 0 0 0 0 0 2 14 H2(0,0) 0.25 0 2 14 H2(0,0) 0.25 19 0 1 0 0 0 0 5 5 3 3 0 0 3 0 H2(0,0) 0.25 0 H2(0,0) 0.25 10 0 20 H2(0,0) 0.25 0 20 H2(0,0) 0.25 0 0 12 0 12 0 nr H2(0,0) 0.25 3 H2(0,0) 0.25 45 H2(0,0) 0.25 112 112 685 643 H2(0,0) 0.25 296 24 452 0 5 24 H2(0,0) 0.25 79 79 369 327 5 24 H2(0,0) 0.25 106 47 388 0 0 32 H2(0,0) 0.25 0 32 H2(0,0) 0.25 0 2 0 2 2 1 0 0 Table 4.1: Summary of the collision outcomes for H2 + H2 molecular interactions for the lowest impact parameter range. The data reflects the number of trajectories computed for each set of vibrational (v) and rotational (j) quantum states, including reactive (nr ), reactive dissociative (nrd ), dissociative (nd ) and double dissociative (ndd ) events. The table displays the total internal energy in Eh . Results from the Hinde PES are shown in bold. 32 This study brings to light significant concerns regarding ability of the Hinde PES to model dissociation processes in molecular hydrogen interactions. A considerable deviation has been observed, with dissociation events occurring at energies below the required threshold. This observation suggests a deficiency in the Hinde PES, particularly in high-energy regimes where dissociation events can occur. Further evidence of this issue is presented in Figure 4.1 and Figure 4.2, which displays the apparent double dissociation cross-sections for specific states of the hydrogen molecule, namely (5,24) and (12,00). The cross-sections, measured in a20 , reveals a stark contrast in outcomes between the BMKP2 and Hinde PES. Notably, the Hinde PES forecasts a significantly higher probability of double dissociation at elevated energy levels. However, the total energy present appears insufficient to support such doubly dissociative events. This unexpected result raises two pivotal questions that this study aims to resolve. Firstly, could an issue with the integrator lead to erroneous predictions of such high probabilities for double dissociation? It is essential to consider whether the numerical method employed for integration could introduce errors, especially at the higher energy calculations where the discrepancies are most apparent. Secondly, a detailed look at the Hinde PES itself is needed. Especially regarding its ability to simulate high-energy dissociative phenomena in molecular hydrogen interactions accurately. Addressing these questions will help assess why the Hinde PES is giving such unphysical results. 33 0 - r--..--..---.--~--r-~--r----r---.---.---.-----.-----.----.---,----,-, BMKP2 ........ Hinde 0 ..............................-.,........................--..ii:::=...:i.J 0 0.1 0.3 0.2 Figure 4.1: The apparent double dissociation cross sections at various translational energy levels for (5,24) (v, j) state. The vertical axis is cross section in a20 and horizontal axis is the translational energy in Eh . The arrowhead represents the translational energy for H2 (0,0) and H2 (5,24) having enough to achieve double dissociation. I{) 0 N 0 0 -....... 0 BMKP2 Hinde ...,_._.............i.............¥.....~------i:::l!:::............J 0 0 .1 0.2 0.3 Figure 4.2: The apparent double dissociation cross sections at various translational energy levels for (12,0) (v, j) state. The vertical axis is cross section in a20 and horizontal axis is the translational energy in Eh . The arrowhead represents the translational energy for H2 (0,0) and H2 (12,0) having enough to achieve double dissociation. 34 4.0.1 Numerical Integrator Check Previously, rigorous testing was undertaken to ensure that the integration conditions were not the source of computational inaccuracies. This diligence is reflected in the accompanying Figure 4.3, which plots the relative conservation of energy against the time-step duration employed in our simulations. The precision with which energy conservation was achieved, maintained at the machine’s accuracy limit of 1 part per million (ppm), suggests the absence of significant numerical errors in the integration. Aside from small fluctuations, the satisfactory energy conservation reinforces the Runge-Kutta method’s suitability for numerical integration of H2 + H2 trajectories. tO I 0 X N tO I 0 0 tO I 0 tO I 0 ..... .__..............................____.____.___.____.____.____.___.__.__._....._....._.....__.....__.....___....___. ~o I 100 200 300 400 Figure 4.3: This plot of relative error (∆E/Etot ) against time-step demonstrates the stability of the numerical integration method employed. The plots provided comprehensively analyse the behaviour of the system under simulation. Figure 4.3 showcases the deviation of a single trajectory, illustrat- 35 ing the dynamic nature of the simulation at an individual level. Figure 4.4 expands upon this by presenting a statistical average of deviations across multiple trajectories, offering a macroscopic view of the stability of the system over time. This plot elucidates the overall consistency and reliability of the simulation by averaging out the individual variations. Figure 4.5 furthers this analysis by considering the average of absolute values of the deviations. This plot provides a critical measure of the total deviation experienced by the system, thereby offering an integrated check on the precision of the simulation. Such a measure is particularly appropriate for assessing the effectiveness of numerical integrators, as it captures the full scope of the deviation. The collective insights from these plots affirm the efficacy of the numerical integrator within the established tolerances. The numerical integrator check has shown us that the integrator is working as it should be within the tolerance, meaning the issue lies within the potential. 0 ,..._ I ....0X N ,-..I I ....0X -.:t- ,-..I I ....0X lO ,-..I I 0 X IX) ID I I 0 10 100 200 300 400 Figure 4.4: This plot illustrates the change in average deviation over thirty trajectories against time-step, providing insight into the stability and convergence of the system over the simulation period. 36 co I 0 X N co I 0 r- x I() co I 0 ,.._ I 0 r- x I() 0 ...................___.___.___.____.___.___.___._...................................................._...._...__........... 0 100 200 300 400 Figure 4.5: This figure traces the average of the absolute deviation over thirty trajectories against time-steps, clearly showing how closely the simulation follows the expected trajectory. 4.1 Failure of the Hinde PES when R is small A particular focus is required on conformations involving small separation between the two molecules in stretched states, as these conditions are crucial for understanding collision induced dissociation. For this purpose, an in-depth exploration of the Hinde PES was conducted, explicitly considering instances of small intermolecular separation R and assessing the complete range of potential values for r1 and r2 , the separations within each diatomic molecule. This investigation is imperative for assessing whether the PES used in molecular dynamics simulations accurately captures both the physical reality of molecular interactions. 37 4.1.1 Potential plots for the Hinde PES when R is small A survey of the Hinde Potential PES was conducted to elucidate its unphysical behaviour, using three different geometries for small separation of the two molecules. The first one is the collinear arrangement. This geometry is characterized by the internuclear separations r1 and r2 within the two diatomic molecules, respectively, while R remains a fixed value representing the distance between the centers of mass of the two molecules. The plots examined, as depicted in the provided Figure 4.6 effectively map the internuclear separations r2 and r1 , with r2 plotted along the vertical axis and r1 along the horizontal axis. This arrangement allows for an comparison of the of the BMKP2 and Hinde PES. In this series of plots, the first column illustrates the BMKP2 PES, the second column portrays the Hinde PES and the third column, shows the differences in the two potential energy surfaces which is the value of the BMKP2 PES subtracted from the value of Hinde PES. All measurements and comparisons in these plots are made in atomic units. 38 R 0 0 0 co co co <.O <.O <.O '