Nitramine Conformers and Detonation Mechanism Via Standard Coupled Cluster and Double-Electron Attached Coupled Cluster ...


Material Information

Nitramine Conformers and Detonation Mechanism Via Standard Coupled Cluster and Double-Electron Attached Coupled Cluster Theory
Physical Description:
1 online resource (121 p.)
Molt, Robert W
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Thesis/Dissertation Information

Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Committee Chair:
Committee Co-Chair:
Committee Members:


Subjects / Keywords:
Chemistry -- Dissertations, Academic -- UF
Chemistry thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation


The mechanism of detonation and conformational diversity of nitramine explosives have been explored using coupled-cluster techniques, old and new. Experimental efforts to unveil the mechanism have been contradictory for decades; we look to modern methods of calculating wave functions systematically for clarity. We have used CCSD(T) to describe the energies of the single-reference systems, a tried-and-true method known as the gold standard of quantum chemistry. Energetic minima along the nitramine potential energy surface have been thus characterized. Transition states are intrinsically multi-reference systems for which even CCSD(T) may be questionable. We have validated the reliability of the EOM-DEA-CCSD method for notoriously pathological molecules with the same static symmetry as most transition states in chemistry. EOM-DEA-CCSD performs spectacularly, and has thus been used to analyze the nitramine decomposition mechanism. We have achieved a characterization of the nitramine potential energy surface between different conformers and degradation paths among the family of nitramine explosive.
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility:
by Robert W Molt.
Thesis (Ph.D.)--University of Florida, 2013.
Electronic Access:

Record Information

Source Institution:
Rights Management:
Applicable rights reserved.
lcc - LD1780 2013
System ID:

This item is only available as the following downloads:

Full Text




2 2013 Robert William Molt Jr


3 To Dr. Caroline Joan Picart and Stephanie Frangos Wilson: A Mentor and a Mother who suppor ted me in the difficult times, without whom my thesis would not be possible


4 ACKNOWLEDGMENTS I offer my gratitude to the following individuals for their scientific support over the years. Dedra Demaree You were my firs t true role model of what a scientist should be. You taught me what it was to be a professional. You showed me what true dedication is to science. You believed in me when others did not. I would have quit science if not for you. Paul Oxley Everything I know about doing science and being a PI I learned from you. You were a terrific advisor; you knew when to support me, when to punish me, when to push me, when to help me. Michael Tate You are the best combination of scientist, statesman, p hilosopher, and fun friend yet known. Doing problems s ets with you was fun despite Goldstein, Jackson, Shankar, and Kittel. Rebecca Ansari The kindness you offered me is ineffable, and I owe you in life. You gave me the strength to remain scientificall y honest amidst scientists who cheated. Matthew Strasberg, To m Watson, and Alexandre Bazante You are the 3 smartest graduate students I have known. If you three choose it, you will have fantastic scientific careers. I admire you intellectually, to say nothing of what great guys you are. Victor Lotrich If not for you, I would have not lasted these past two years. You are a great scientist because you emphasize the simplicity of science with neither pretense nor arrogance. You are practical and kind. I hope to work with you again soon. Rodney Bartlett and Ajith Perera I thank you for choosing to take me on as your research student.


5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 LIST OF ABBREVIATIONS ................................ ................................ ........................... 10 ABSTRACT ................................ ................................ ................................ ................... 11 CHAPTER 1 NITRAMINES AND THE NEED FOR AB INITIO METHODOLOGIES .................... 12 RDX ................................ ................................ ................................ ........................ 14 HMX ................................ ................................ ................................ ........................ 16 CL 20 ................................ ................................ ................................ ...................... 17 DEA CCSD ................................ ................................ ................................ ............. 19 RDX Kinetics ................................ ................................ ................................ ........... 20 2 RDX CONFORMERS AND HOMOLYTIC N N DISSOCIATION ............................. 21 Relevant Literature ................................ ................................ ................................ 21 Computational Methodology ................................ ................................ ................... 25 Structural Analysis and Stationary States ................................ ............................... 26 Energetics ................................ ................................ ................................ ............... 29 Electronic Excitation Spectra ................................ ................................ .................. 33 3 HMX CONFORMERS AND CRYSTAL STRUCTURES ................................ .......... 36 Relevant Literature ................................ ................................ ................................ 36 Computationa l Methodology ................................ ................................ ................... 38 Searching Over the Potential Energy Surface ................................ ........................ 39 Resulting Conformers and Nomenclature ................................ ............................... 42 4 CL 20 CONFORMERS AND CRYSTAL STRUCTURES ................................ ........ 48 Relevant Literature ................................ ................................ ................................ 48 Methodology ................................ ................................ ................................ ........... 49 The Six Conformers ................................ ................................ ................................ 51 Energetics and Geometric Speculation Therein ................................ ...................... 56 Implications Toward Reactivity ................................ ................................ ............... 58


6 5 DEA CCSD CONQUERING PATHOLOGICAL MOLECULES ................................ 61 The Multi reference Problem ................................ ................................ .................. 61 Multi reference Coupled Cluster Theory in Context ................................ ................ 61 Theory ................................ ................................ ................................ ..................... 64 Methodology ................................ ................................ ................................ ........... 67 Methylene ................................ ................................ ................................ ............... 68 Trimethylenemthane ................................ ................................ ............................... 74 Benzynes ................................ ................................ ................................ ................ 81 Further Thoughts in Conc lusion ................................ ................................ .............. 90 6 MECHANISM OF RDX DECOMPOSITION ................................ ............................ 91 The Stage Has Been Set, Let Us Dance ................................ ................................ 91 Methodology ................................ ................................ ................................ ........... 91 Proposed Mechanisms ................................ ................................ ........................... 92 N N Homolysi s ................................ ................................ ................................ ........ 94 HONO Elimination ................................ ................................ ................................ 100 Triple Whammy ................................ ................................ ................................ ..... 103 Validation of Methodology on Transition States ................................ .................... 103 Further Work ................................ ................................ ................................ ......... 1 04 Concluding Remarks ................................ ................................ ............................ 106 LIST OF REFERENCES ................................ ................................ ............................. 107 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 121


7 LIST OF TABLES Table page 2 1 MP2/6 311++G(d,p) level of theory for Four Lowest Harmonic Vibrational Frequencies of RDX ................................ ................................ ........................... 27 2 2 Relative Energies of RDX Conformers Using Different Many Body Methods ..... 29 2 3 Relative CCSD(T) Energies of RDR Conformers ................................ ............... 32 2 4 Relative Energies of RDR Conformers Relative to RDX ................................ ..... 32 2 5 EE EOM CCSD and CIS Excitation Energies, eV ................................ .............. 34 3 1 Stationary States of the HMX Molecule ................................ .............................. 42 4 1 Energetics of the CL 20 Conformers ................................ ................................ .. 56 4 2 Bond Lengths of CL 20 Conformers, ................................ ................................ 58 4 3 Lowest 4 Harmonic Vibrational Frequencies Calcu lated using MBPT(2)/6 311G(d,p) ................................ ................................ ................................ ........... 59 5 1 Methylene Singlet Triplet Splittings, kcal/mol ................................ ..................... 70 5 2 Methylene Singlet/Triplet gap at CCSD(T)/6 311++G(2d,2p) Geometry. ........... 73 5 3 Geometry of 1 A 1 C 2v structure. ................................ ................................ ............ 76 5 4 Geometry of 1/3 B 2 C 2v / 3 A 2 D 3h structure ................................ ............................. 77 5 5 TMM Singlet Triplet Splittings. ................................ ................................ ............ 79 5 6 TMM Singlet Triplet Splittings at CCSD(T)/6 311++G(2d,2p) Geometries ......... 80 5 7 Geometric Parameters of Para Benzyne, Singlet followed by Triplet ................. 82 5 8 Geometries of Meta Benzyne ................................ ................................ ............. 86 5 9 Meta and Para Benzyne Splittings ................................ ................................ ...... 88 5 10 Meta and Para Benzyne Singlet/Triplet Splittings, kcal/mol. ............................... 88 5 11 HOMO/LUMO Gaps of Benzynes ................................ ................................ ....... 90 6 1 RDX and RDR Energies, kcal/mol ................................ ................................ ...... 95 6 2 DEA CCSD and CCSD(T) cc pVDZ Energies on Transition States ................. 104


8 LIST OF FIGURES Figure page 1 1 The Three Major Nitramines: RDX, HMX, and CL 20 ................................ ......... 12 2 1 Relative Energies of RDX Conformers, kcal/mol ................................ ................ 21 3 1 Boat Boat ................................ ................................ ................................ ........... 40 3 2 Chair ................................ ................................ ................................ ................... 40 3 3 Chair Chair ................................ ................................ ................................ ......... 41 3 4 Crown ................................ ................................ ................................ ................. 41 4 1 Side View of Equatorial/Axial CL 20 ................................ ................................ ... 51 4 2 Side View of Diaxial CL 20. ................................ ................................ ................ 51 4 3 ................................ ................................ ............ 52 4 4 ................................ ........ 52 4 5 equatorial. ................................ ................................ ................. 53 4 6 ................................ ................................ ........... 53 4 7 equatorial ................................ ................................ .................. 53 4 8 Transition State Conformer; diaxial ................................ ................................ .... 54 5 1 Reference Determinant and Achieved Static Correlation in the DEA CCSD Scheme. ................................ ................................ ................................ ............. 65 5 2 1 A 1 C 2v structure. ................................ ................................ ................................ 75 5 3 3 B 1 C 2v / 3 A 2 D 3h structure. ................................ ................................ ................... 77 5 4 Para Benzyne and our associated axis choice. ................................ .................. 82 5 5 Orbital diagram of Para Benzyne. ................................ ................................ ...... 83 5 6 Meta Benzyne. ................................ ................................ ................................ ... 85 6 1 N N Homoly sis Mechanism ................................ ................................ ................ 93 6 2 HONO Elimination Mechanism ................................ ................................ ........... 93


9 6 3 Triple Whammy Mechanism ................................ ................................ ............... 93 6 4 Secondary N N Homolysis ................................ ................................ ................. 96 6 5 CN Ring Scission from the Secondary N N Homolysis ................................ ....... 96 6 6 Primary CN Bond Cleavage of RDR ................................ ................................ ... 97 6 7 Second CN Cleavage of RDR ................................ ................................ ............ 98 6 8 CN Cleavage of Twice CN Cleaved RDR ................................ ........................... 98 6 9 I somerization of CN Cleaved RDR ................................ ................................ ..... 99 6 10 CN Cleavage of the Most Stable N N Homolysis Path ................................ ....... 99 6 11 Second HONO Elimination ................................ ................................ ............... 100 6 12 Third HONO Elimination ................................ ................................ ................... 101 6 13 Simultaneous Heterolytic Bond Cleavage ................................ ........................ 102 6 14 Simultaneous Heterolytic Bond Cleavage on Pre Aromatic Product ................ 102


10 LIST OF ABBREVIATIONS CCSD Coupled Cluster Singles and Doubles Refers to a choice in operator to sample excitation space of a many body problem, T=T 1 +T 2 CCSD(T) Coupled Cluster Singles and Doubles and perturbative Triples CCSD with the use of perturbation theory to approximate the effect of the T 3 operator; it thus samples excitation space in a manner in between CCSD and CCSDT DEA CCSD Double Electron Attached CCSD Ref ers to the use of the R operator to sample excitations of a coupled cluster wavefunction reference for two added particles to the reference.


11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Part ial Fulfillment of the Requirements for the Degree of Doctor of Philosophy NITRAMINE CONFORMERS AND DETONATION MECHANISM VIA STANDARD COUPLED CLUSTER AND DOUBLE ELECTRON ATTACHED COUPLED CLUSTER THEORY By Robert William Molt Jr December 2013 Chair: Rodney Bartlett Major: Chemistry The mechanism of detonation and conformational diversity of nitramine explosives have been explored using coupled cl uster techniques, old and new. Experimental efforts to unveil the mechanism have been contradictory for decades; we look to modern methods of calculating wavefunctions systematically for clarity. We have used CCSD(T) to describe the energies of the single reference systems, a tried and true method known as the gold standard of quantum chemistry. Energetic minima along the nitramine potential energy surface have been thus characterized. Transition states are intrinsically multi reference systems for which ev en CCSD(T) may be questionable. We have validated the reliability of the EOM DEA CCSD method for notoriously pathological molecules with the same static symmetry as most transition states in chemistry. EOM DEA CCSD performs spectacularly, and has thus been used to analyze the nitramine decomposition mechanism. We have achieved a characterization of the nitramine potential energy surface between different conforme rs and degradation paths for the family of nitramine explosive s


12 CHAPTER 1 NITRAMINES AND THE NE ED FOR AB INITIO METHODOLOGIES Nitramines are an exotic functional group (R 2 N NO 2 ) characteristic of a class of explosives. The three major nitramine molecules are RDX, HMX, and CL 20; they are featured in Figure 1 1 Note that the basic units for RDX and HMX are cyclohexane and cyclooctane, respectively Figure 1 1 The Three Major Nitramines: RDX, HMX, and CL 20 The RDX molecule is the basic constituent of C4 explosive, and CL 20 is one of the current champions in exergonic yield of chemical explosive s. It is always desirable, however, to have explosives which have adjustable shock sensitivity (readiness to explode upon mechanical stress) and a larger explosive yield than current explosives


13 The shock sensitivity is indeed a function of multiple vari ables from the molecular level to the solid state, but an important governing factor is the kinetic mechanism. The hardest steps will be the most energetically demanding; intramolecular forces are much greater than intermolecular. The larger scale factor s (surface area, crystal polymorphism, etc.) can be adjusted more malleably and are secondary effects; the intrinsic qualities begin with the chemistry. The mechanism for decomposition is unknown, despite 80 years of research. Establishing this mechanism of decomposition would allow us to adjust shock sensitivity to our purposes (methylate an electron deficient transition state to increase sensitivity, for example). Whether we intend to increase or decrease shock sensitivity depends on the specific milit ary situation in question; we merely offer the knowledge of how to adjust this. The explosive yield (defined in terms of either the negative enthalpy of combustion or the negative of Gibbs energy of combustion) is a more complex feature to optimize. The s ingle closest parameter to correlate explosive yield is crystal density (although it can mislead on occasion ). There is thus a great value in predicting new crystal structures, as this may lead to greater explosive yield. In organic crystals with weak i ntermolecular forces, the degree of perturbation in geometry between the gas phase conformer and the unit cell monomers is small. Consequently, in finding new conformers, we may find new crystal structures, and perhaps more efficient explosive forms. Expe rimental studies of the mechanism have been inconclusive; contradictory evidence is everywhere 1 21 In using methodologies as rigorous as the gold standard of quantum chemistry, coupled cluster singles, doubles, and perturbative triples,


14 CCSD(T) 123 125, 1 52 154 we can address definitively which mechanisms are energetically favored in the gas phase. The value of such methods lies in that they are universal descriptions of wavefunctions for which the system is single reference in nature. When the system i s NOT single reference in nature, as is the case for transition states, we will use the double electron attached coupled cluster singles and doubles (DEA CCSD) method. For all transition states in which there are only two degenerate states, DEA CCSD captu res the static correlation present precisely, and reli es on coupled cluster to recover the majority of the dynamic correlation. Additionally, the use of MBPT(2) (a subset of the CCSD equations, truncated for greater cost efficiency at an accuracy sacrific e) for geometries of gas phase conformers will allow us to identify new crystalline structures for the various nitramine units. RDX RDX (formally known as 1,3,5 trinitro 1,3,5 triazacyclohexane, informally as royal demolition explosive), as shown in Figure 1 1 is the subject of intensive research as a result of its uses as an explosive. Above and beyond its own particular uses, it is the smallest and simplest chemical structure that has the essential qualities of some more powerful explosives, such a s HMX (octahydro 1,3,5,7 tetranitro 1,3,5,7 tetraazocine, or hexanitro 2,4,6, 8,10,12 hexaazaisowurtzitane, or China Lake compound 20, also shown in Figure 1 1 ). Many investigations have been conducted conce rning the RDX structure, 1 9 and even more concerning the mechanism of decomposition 10 21 Knowledge of form and function (structure and exploding) go hand in hand. Our work here serves to improve the understanding of the RDX structure and, based on this, to understand the mechanism of decomposition. This will help to design more stable or more energetic materials.


15 The mechanism of decomposition of RDX has been studied experimentally for over 80 years. It is established that multiple decomposition pathways exist; however, the viability of the various pathways is still unclear. Some mechanisms exist in one phase, but not in others; that is, the gas phase mechanism involves NO 2 formation, but the condensed phase seemingly does not 16 It was proposed that th e mechanism begins, in some cases, via homolytic bond dissociation of the N N bond 12 It has been noted that only one nitro group dissociates in this manner per molecule of RDX 17 We are interested in the N N bond cleavage method, in which a radical NO 2 an d the corresponding RDX radical (RDR) is formed, as well as having refined structural information from which to begin a larger mechanistic study of various transition states involved in various proposed mechanisms. It should be noted that implicit in any experimental study is an interpretation of the data with respect to one particular molecular mechanism in mind. This has led to disagreements in the literature about whether some evidence uniquely supports a given mechanism. Additionally, the degree to whi ch experimental contention in interpreting the data. For example, one can study optical versus thermal experiments for the mechanistic trigger; one can study the initial kinetic step or the ki netics that dominate the pathway after the first few steps that change the temperature and the chemical species present. This level of interpretation is exceptionally difficult. Through trustworthy calculation methods, we can offer insight that experiments cannot: unambiguous interpretation of the data (because we have complete control of what is occurring) and the safety of not working with explosives.


16 HMX The HMX molecule is formally known as octahydro 1,3,5,7 tetranitro 1,3,5,7 tetraazocine It is a hete rocyclic cyclooctane, shown in Figure 1 1, along with related members of its chemical family, RDX and CL 20. The tell tale feature of all of these molecules is the nitramine group. This functional group is the seat of the basic explosive properties of this m olecule. Clearly, factors such as crystalline density, phonon magnitude, and many others contribute toward the overall explosive quality, but all of these properties are secondary. The nitramine groups are sine qua non for these explosives. As such, a clea r understanding of reaction kinetics at the chemical level is necessary to appreciate and improve upon the explosive utility of these molecules. The greatest gains in improving the material properties (explosive yield ) can be made in changing the most fund amental properties (chemical), rather than bulk phase. The HMX literature is quite extensive, spanning 60 years. There have been numerous kinetic, structural, explosive, and physical properties studies. Experiments done on the solid state have been careful to distinguish which crystalline polymorph was used; however, the same cannot always be said for the study of the gas phase with respect to conformers. Some computational kinetic studies do acknowledge that there is more than one fo rm of HMX, but choose only one or two conformers for study; 35,36 it would be most systematic and thorough first to consider our procedure below, it is hard to imagi ne other plausible gas phase conformers. This is the primary effort of this study. Once the conformers are known, one may study the kinetics more carefully, as well as consider new crystalline polymorphs. To study the reaction kinetics, it is critical that one know all of the germane conformers. If one wishes


17 to be accurate to within even an order of magnitude in the relative populations of conformers/products to reactants/any two comparative states, one must have an accuracy of about 1.4 kcal/mol according to the Boltzmann distribution. Consequently, one may not speak meaningfully about the kinetics of the process unless one has such accuracy. The energy differences between conformers of related nitramine molecules are easily smaller than this, suggesting t he importance of considering all plausible HMX conformers. 72,73 Additionally, there is a literature history of most gas phase conformers subject to the crystalline field. 1 Given that crystalline field effect s have been observed to be small in related nitramines, 72,73 we might expect HMX crystalline field effects to be modest, all the more speaking to the importance of understanding the gas phase conformers for the discover y of new possible HMX crystalline polymorphs. CL 20 CL 20 (shown in Figure 1 1 along with its chemical brethren) has been referred is statement inherently presages the importance of the study of th is molecule. 82,83 The relevant conformers of CL 20 are an important question to be studied. In chemical ly similar explosives, many conformers discovered in the gas phase were eventually discovered to be the basis for a crystalline polymorph. All 2 groups) compounds is inherently worthy of comparison. What is c hemically true of the chemical cousins of CL 20 may be true of it; however, others dispute the seeming similarities in detonation mechanism based on various kinetic studies. 94, 95 In particular, CL 20 has been studied far less than its RDX and HMX cousins. This is partly because its existence has not


18 be en known as long as RDX and HMX and partly to the fact that it is an even more dangerous explosive than its counterparts (and t hus undesirable experimentally). Moreover, there are restrictions o n inform ation as a classified weapon and on the computational side, its larger size (which inhibits easy computational study). Pace and co workers first studied the change in the concentration of nitro groups in the three major nitramines (RDX, HMX, and CL 20) through EPR experiments .96 They concluded that nitro groups become trapped and stabilized within the CL 20 crystal structure. Oxley and co workers s tudied thermal decomposition of CL 20 in dilute solution. 97 Their studies of dilute solutions had the advantage of inherent avoidance of autocatalys is and multiphase decomposition. Additionally, the samples could be kept isothermal. Indeed, one of the gre atest challenges in the study of RDX and HMX is the ability to generalize a mechanistic study of a single phase to a multiphase system, or the ability to discriminate reactions of one phase amidst data from a multiphase ows that the rate constant under the aforementioned conditions for CL 20 decomposition is about 3.5 times greater than for RDX and about 35 times greater than HMX. The activation energies observed, in each case, are between 40 and 50 kcal/mol. This is sugg estive of a common mechanism; the slight variations in activation energy could be a result of structural differences endemic to RDX, HMX, and CL 20, rather than a different mechanism. Other experimental studies, referenced above, 95 disagree; some evidence suggests that the RDX and HMX mechanisms differ significantly. In all instances, however, we emphasize that different experimental conditions may mean that the results do not contradict one another; the interpretation of the data differs widely in various cases. Ryzhkov and McBride also


19 performed EPR studies of CL 20 crystals and came to a similar conclusion as Oxley and co workers. 98 Rice, Thompson, and Sorescu performed constant particle number, pressure, temperature molecular dynamics calculations (NPT MD) on CL 20 in the crystalline phase using MD parameters associated with RDX. 99 This work demonstrated that the essential qualities of the basic crystal structures could be replicated for CL 20 from the basic physical chemistry of RDX, i.e., that CL 20 is essentially RDX, chemically (PT) curves of all of the nitramines, further validating their MD calculations. 100 However, DFT) cal culations on the solid state of nitramines showed drastic failures 101 in predicting lattice parameters This is a good example of the subtle nuance of computational modeling: success in modeling one property by no means implies successful prediction of oth er properties. Additionally, this is a good example of where common wisdom that KS DFT is better than molecular mechanics is far from true. This is a usual phenomenon in highly parametrized methodologies, such as many (but far from all KS DFT ) functional a nd MM calculations D EA CCSD Multi reference wavefunctions can be made single reference using electron attachment methodologies. By judi cious removal of a certain number of particles in the reference wavefunction, then addition of them during an EOM CCSD calculation, we may circumvent difficulties associated with multi reference wavefunctions. We demonstrate the accuracy of the method on notable multi reference problems: the singlet/triplet splittings of methylene, meta benzyne, para benzyne, and TMM, us ing basis set extrapolation methods and CCSD for dynamic correlation. We also


20 additionally present the geometries of the aforementioned challenging diradicals which are at least as accurate, if not more than, yet known geometries. Comparison is made to o ther multi reference methodologies. As a consequence, any system with the same essential static symmetry of the wavefunction will be well described via DEA CCSD (double electron attached CCSD) This is of particular use for transition states described by one leaving group, one attacking group, as such systems have the same degeneracy patterns as diradicals. RDX Kinetics The final act in appreciating nitramine chemistry is to examine the reaction kinetics of detonation of the simplest nitramine, RDX. This unto itself will allow for the knowledge of how to adjust the shock sensitivity as well as potentially increase explosive yield. We consider the most significant mechanisms in the literature as well as some newly proposed pathways original to this work. In calculating the activation energy barriers, we demonstrate that the putative mechanism of N N homolysis is, in fact, not the dominant mechanism. The lesser favored mechanism of HONO elimination from the RDX ring is also shown to be un tenable due to e xtremely stable aromatic intermediates. The surprise champion is in fact a simultaneous breaking of surprising and fascinating result; all previous experiments and ca lculations failed to find any evidence for this path 10 21 aside from the originator of the idea 15


21 C HAPTER 2 RDX CONFORMERS AND HOMOLYTIC N N DISSOCIATION Relevant Literature X ray crystallography and neutron diffraction have provided the geometries of crystal phases of RDX. 1,7, 113 Various Kohn Sham density functional theory (KS DFT) methods with semi empirical exchange correlation functional have been used to analyze specific proposed molecular geometries. 3 5,8,37,114 115 The resulting conformers are referred to in terms of the nomenclature associated with cyclohexane chemistry, that is, in terms of axial/equatorial chairs and twist/boat. 116 These structures are shown in Figure 1 Figure 2 1 Relative Energies of RDX Conformers, kcal/mol


22 Rice and co workers 3 used B3LYP/6 311 + G(d,p) to show that the form of RDX with two axial nitro groups and one equatorial (referred to as AAE) closely co rresponds crystal phase RDX. They also showed crystal phase and vapor phase RDX correspond to three axial nitro groups (AAA). However, the observed matching in geometries and IR spectra between AAA and the gas phase structures was found to be incidental by Torres and co workers; 114 the Raman spectra and depolarization of AAA gas phase electronic structure calculation versus experimental RDX gas phase do not match. Harris and co workers 4 studied five of the six conformers using B3LYP/6 31G(d) for geometries and frequencies. Wu and Fried 115 effectively tried a variety of KS DFT func tional using the cc use of any one functional. It also seems they did not state which conformer of RDX they were considering in their mechanistic study. Chakraborty and co workers 37 studied conformers using the same computational methodology as Harris and co workers. 4 Refinement of these geometries, frequencies, and energies is necessary to describe accurately the mechanism in question. As noted already for the Boltzmann distribution with an assumed degeneracy factor of 1, the ratio of products to reactants of even a 10:1 ratio corresponds to 1.4 kcal mol 1 energy difference. Thus, to estimate the viability of one pathway relative to another and be right within a factor of 10 for the populations of molecules, one requ ires accuracy beyond the capabilities of the afore mention ed computational studies. 117, 118 say, the accuracy necessary to describe chemical populations, one needs to use better many body methods with more functions i n the basis. Almost all of these prior


23 calculations use B3LYP, which is disconcerting in light of recent work by Shimojo and co workers 8 as well as by Byrd and co workers, 5 that shows that functionals not specifically parametrized for dispersion will not a ccurately predict crystal properties if dispersion is significant B3LYP does not have chemical accuracy when subject to systematic tests. 117 119 In more detail, work using Quantum Monte Carlo (QMC) and the ab initio dft method (not to be confused with DFT methods, emphasis on the capitalization) of Bartlett and co workers has shown that most all KS DFT potentials get the signs of exchange potential versus correlation potential wrong. 120 121 Additionally, semi empirical KS DFT studies are only reliable in s o much as one believes one has parametrized for the molecule in question. The lack of systematic correct ability of KS DFT and difficulty in predicting whether one is within the parametrization range of a functional led us to methods such as many body perturbation theory and coupled cluster theory. 1 22 123 The use of methods that are systematically improvable allows us to investigate RDX with more confidence. The basic omission of Hartree Fock methodology is, by definition, the correlation energy. Coupled cluster as a many body method in the form of CCSD(T) achieves greater than 99% of the correlation energy in single reference molecules. 124 125 As such, within the confines of eigenfunction expansion of basis vectors 110 111 and the degree to whic h the molecule is of single reference character CCSD(T) energies are rigorous; all that remains is to ensure that a sufficiently large basis set is used. Within the confines of single reference molecules, second order perturbation theory geometries have b een shown to be accurate to within +/ 0.15 pm for bond lengths using a triple 12 6 We thus rely on such methodology for our own


24 geometry optimization. We acknowledge that harmonic frequencies, however, remain a difficult issue to guarantee for any method short of CCSD(T). The closely related theory, CCSD[T], is generally trustworthy within 31 cm 1 MP2 to within 70 cm 1 assuming the molecule reference organic molecule. 108 We als o emphasize that the extent to which harmonic frequencies sufficiently account for the total frequencies is less than clear for an arbitrary molecule. It is nonetheless useful to calculate harmonic frequencies for the sake of being qualitatively internally consistent (QIC); that is, there is great confidence that the calculation will predict whether a point is a minimum or maximum accurately. The ability to perform such calculations as MP2 and CCSD(T) with triple basis sets is no small feat computational ly. Leaders in the computational studies of RDX, such as Betsy Rice, have noted as recently as 2007 the intractability of using CCSD(T) on RDX molecules. 9 However, the ACESIII software 102 was designed to handle large molecules with methods such as perturb ation theory and coupled cluster quickly. The previous studies that used KS DFT had to speculate about energy differences that are below the resolution of KS DFT with the assumed functionals based on generalized studies. We thus seek to use higher resolut ion methodology. The improved structures, frequencies, and energies allow us to interpret how the structural features of RDX relate to the decomposition mechanism. In particular, certain conformers may be more likely to be the starting point for certain me chanistic pathways. The improved geometries, frequencies, and energies found for RDX serve as a foundation to begin a mechanistic study in future work for transition states and activation barriers.


25 There has not been any study of which RDX structure corres ponds to the proposed homolytic bond cleavage. It is not obvious that there is or is not a dependence on whether the dissociating nitro group is an axial or equatorial nitro group. It is also not clear what the resultant geometry is of the newly formed radical ringed species in terms of equatorial versus axial positions of the remaining two nitro groups. Here, we address all these issues. Computational Methodology Starting from the KS DFT geometries of RDX of Rice and Chabalowski, we optimized the geomet ries at the MP2/6 311 ++ G(d,p) level. 106, 127 This choice in method was appropriate in light of the accuracies associated with MP2 with a triple set mentioned in the Introduction. 126 The calculation of harmonic frequencies was conducted at the same level to confirm minima rather than maxima for the stationary points found through optimization, as well as to analyze the vibrational qualities with respect to the reaction mechanism. Single point energy calculations then were conducted using CCSD(T) with both 6 311 ++ G(d,p) and cc pVTZ bases. 104 We opted for as many polarization functions as we could add, given the strong dependence of polarization functions with CCSD(T), and hence opted for 104 This is also important for the s ake of future work, in which basis set extrapolation methods may be used, since basis set extrapolation methods. For single point energy calculations, core molecular orbitals were dropped due to their irrelevance to the energy relative to the accuracy of our study. 128 All radical species used an unrestricted reference, as spin contamination is unlikely to be an issue for simple doublets using MP2 or CCSD(T). 124 Electronic excited states were calcu lated using equation of motion (EOM) coupled cluster for single and


26 double excitations (EOM CCSD). 129 We opted to use the 6 311 (2 + ,2 + )G(d,p) basis sets for excited states, based on the strong performance of this basis set relative to its cost (in contrast to, say, aug cc pVTZ). 130 All cores were also dropped in the excited state calculations. We note the recent work of Wiberg et al., who established that there is not a need for more polarization functions in the basis set of an EOM CCSD calculation to matc h experimental results. 131 The starting geometries for the RDR (RDX minus one nitro group, resultant from N N homolysis) molecules were obtained by removing the distinct NO 2 groups from the RDX optimized structures. These structures were then optimized at the MP2/6 311 ++ G(d,p) level, and we used only CCSD(T)/6 311 ++ G (d,p) for further single point energy computations. The rationale for this choice will be presented later. When a many body system is not dominated by a single reference, the T 2 amplitudes in t he coupled cluster equations will rise dramatically. In all calculations performed, we made sure to monitor the values of the T 2 amplitudes to make sure the molecules were dominated by a single reference. The T 2 amplitudes never rose to a troubling degree, a solid indication that we are within the range of validity of the method. Structural Analysis and Stationary States Rice and co workers noted that, based on the near indistinguishability of crystal geometries versus KS DFT molecular geometries, the cryst al field must exert very little influence on the structure. 3 Harris and Lammertsma 4 disagreed, based on their interpretation of the data. They stated that the crystal field effect must be quite strong to account for the variation in angles between the experimental crystal data and their gas phase data. +/ 0.01 of crystal structures; results for our bond lengths and bond angles are virtually s. On the


27 interp retation of negligibility of crystal field effects on geometry; given that our methods agree, and the success of our methods in generalized tests, there can be little doubt that this is the correct conclusion. Small angle variation can have negligible ener gy effects. This finding resolves the contradiction in the literature between Rice et al. and Harris et al. Our calculations on the harmonic vibrational frequencies are enlightening about the conformers. Previous work on harmonic frequencies has been done using B3LYP/6 31G(d) for AAA by Harris and Lammertsma and 6 311+ G(d,p) by Rice and Chabalowski. An inherent problem in using KS DFT is that numeric al integration is almost always required, unlike ab initio methods (excepting explicitly correlated methods). This is not a practical problem in energies or geometries, in general, but is a problem for harmonic frequencies if the frequencies are small. This is especially a concern given that both of the aforementioned studies used the default integration grid of Gaussian 94 132 which is acknowledged to have too few points by Gaussian Inc. 132, 133 W e list the smallest frequencies for each conformer in Table 2 1, as these are most relevant for the sake of establishing minima/maxima. Table 2 1 MP2/6 311++G(d,p) level of theory for Four Lowest Harmonic Vibrational Frequencies of RDX Conformer 1 2 3 4 AAE 32 64 65 71 Twist 62 65 75 92 AAA 20i 18 18 110 EEA 50 57 64 72 Boat 27i 41 56 75 EEE 57 57 66 84


28 Frequencie s are in wavenumbers. We report all values to the ones place without claiming we have accuracy in the harmonic modes to the ones place. Rather, we report these for the sake of being QIC. Absolute values for frequencies predicted for MP2 with a triple expected to be more accurate than those from B3LYP. 108 109 However, in light of the numeric integration grid issue of Gaussian 94, and the fact that ab initio methods do not suffer this problem, MP2 frequencies are qualitatively better in examining the sign of small frequencies. We thus confirm that all of the conformers, save the boat and AAA form,are minima, in agreement with Rice. 112 Note the similarity to cyclohexane chemistry: i n both cyclohexanes and in RDX, the boat form is a transition state, and all other for ms are stationary states. 116 We do not, however, make any firm statement concerning whether the AAA form is an energetic maximum or minimum. The AAA form we isolated was a maximum, based on one imaginary frequency. The frequency in question corresponds to a rotation of two nitro groups with respect to one another. It is possible that the AAA form may be essentially kept with a small rotation of nitro groups with respect to one another; resolving that issue would require careful analysis of the potential en ergy surface a question we leave open It is not crucial to this study, as will be shown in the next section on Energetics. The product of the N N bond homolysis is referred to as RDR. RDR has only previously been studied with very modest calculations, 37 insufficient to rely on with confidence. Moreover, this study only considered two forms of RDR, when there are, in fact, many conformers. We have systematically considered every RDR cyclohexane geometry with accuracy within the aforementioned lesser uncer tainties, by contrast. We


29 considered all eight distinct geometries that could emerge from removing each distinct nitro group from all RDX conformers as an initial guess. W e label these forms analogously with RDX: AA refers to a cyclohexane ring with two ax ial nitro groups, etc. There are three distinct twist cyclohexane forms that result. We note that one half of the geometries only converged to an rms force of 10 4 Hartree/Bohr The other one half of the conformers (AE, cis boat form, and two of the twist forms) converged to 3 x 10 4 with a maximum of 3.2 x 10 4 iterations with negligible change in the geometry. This is indicative of a very flat potential energy surface and suggests that these geometries are more than adequate to provide detailed energetic profiles for the RDR conformers. Vibrational analyses were not performed on the RDR conformers, however. Energetics In Table 2 2 we present the electronic energies of the RDX conformers from MP2/6 31 1++ G(d,p), CCSD(T )/6 311 ++ G(d,p), and CCSD(T)/cc pVTZ calculations Table 2 2 Relative Energies of RDX Conformers Using Different Many Body Methods In kcal/mol Conformer CCSD(T) a CCSD(T) b MP2 a MP2 b B3LYP a B3LYP b 6 311++G(d,p) AAE 0.0 0.0 0.0 0.0 0.0 0.0 AAA 1.8 1.5 1.2 0.9 1.6 0.6 EEE 7.7 7.3 8.3 7.8 5.7 4.6 AEE 3.0 2.6 3.2 2.9 1.7 Twist 0.4 0.7 0.2 0.4 0.0 Boat 4.4 4.0 4.1 3.8 2.5 We also note the electronic energies predicted by the largest basis set employed using B3LYP from Rice and Chabalowski. 3 The MP2 energies are presented gi ven that


30 they were available as a result of geometry optimizations and are an important check on B3LYP energies. By contrast, as presented in the Introduction, CCSD(T) offers a definitive answer for single reference molec ules in so much as the basis set is complete. To this end, we have calculated the energy using the smaller 6 311++ G( d,p) and the larger cc pVTZ with CCSD(T) to observe the energy differences. In each column, we arbitrarily select a different zero of energy based on the deepest bound state electronic energy conformer. We also consider the zero point energies as calculated using MP2/6 311 ++ G(d,p) for the total energy at 0 K. We emphasize the importance of studying the ZPEs, given that small differences in ZPE s can drastically influence the energy ordering when conformers are so close in electronic energy. Rice and Chabalowski 3 concluded that the energy differences were small and that it was unclear which conformer was the minimum in terms of the electronic ene rgies and zero point energies Given the uncertainty range commonly up to 5 kcal mol 1 for B3LYP 117, 119 the results were not capable of solving this issue We confirm some of their conclusions, in that the AAE and AAA forms are close in energy. However, we also note that the twist form is just as likely to be the minimum conformer as the AAE form (which incidentally agrees with calculations done by Harris and Lammertsma 4 ). The EEA form is also close in electronic energy. If one considers the effect of zer o point energies (ZPEs) at the MP2 level, one sees that all conformers have approximately the same ZPEs (a difference of at most 0.7 kcal mol 1 ). Although important to note for the sake of establishing the true minimum energy conformer, this much accuracy is beyond reasonable expectation of CCSD(T)/cc pVTZ and effectively


31 issues. We show the comparative energies of the various conformers in Figure 2 1 The accuracy of CCSD(T) energies at MP2 geometries cannot be rig orously exp ected to be accurate to within 1 kcal mol 1 for an arbitrary single reference molecule. 128 Future work done on a more refined geometry, perhaps CCSD(T) geometries, may be able to answer this more nearly definitively, especially if one incorporat es a basis set extrapolation method. Our results give confidence that four of the six conformers are competitive for the minimum and close in energy. However, the question of which is the absolute minimum is not relevant to the mechan istic question of deco mposition. W ith such small energy differences, variations in temperature and solution dynamics make it likely that all of these conformers will be present significantly, excepting those higher than 5 kcal mol 1 in energy. This finding agrees with Vladimiroff and Rice, who came to the same conclusion via electron diff raction experiments and by use of a crude B3LYP/6 31G(d) estimate of energies 112 It is also consistent with experimental observation of H NMR spectra, 134 which show only one signal, in dicating that the protons are either identical by symmetry or time averaged to be the same. Given the competitive energies of several conformers of nonequivalent hydrogens, it must be a time average. That there is a time average implies fast kinetics of in terconversion between various forms, indicating a low barrier. It is reasonable to propose beginning a mechanism with any one of the four lowest lying stationary states, given how floppy the conformers are. We list in Table 2 3 the energies of the various RDR conformers, with the converged geometries in bold.


32 Table 2 3 Relative CCSD(T) Energies of RDR Conformers Conformer Energy (kcal/mol) AA 1.0 AE 0.6 EE 3.8 cis boat 4.9 trans boat 0.7 twist 1 0.2 twist 2 1.3 twist 3 0.0 W e note that, again, the conformers are all reasonably close in energy. The conclusions we have reached concerning RDX, mechanistically, equally apply to RDR. Only the EE and cis boat forms are likely to be minor contributors. We did not consider the CCSD(T)/cc pVTZ single point calculations for the RDR conformers, in light of all previous evidence. Given how floppy the molecule is, and that CCSD(T) calculations cannot be trusted to be accurate to within 1 kcal mol 1 without basis set extrapolation for a single reference molecule, there is little value in using only one more set of polarization functions. In Table 2 4, we see the comparative energies of the RDR conformers with respect to the RDX conformers. Table 2 4 Relative Energies of RDR Conformers Relative to RDX RDR Energy (k cal/mol) RDX Energy (kcal/mol) AA 46.3 EEA 3.0 AE 45.9 AAA 1.8 EE 49.0 AAE 0.0 cis boat 50.1 Boat 4.4 trans boat 46.0 EEE 7.7 twist 1 45.4 Twist 0.4 twist 2 46.6 twist 3 45.3


33 We neglect the ZPEs of RDR in this table, because the ZPEs for all RDX/RDR forms are comparable and therefore do not change the ordering relative to the accuracy of the calculation method. We see that there is insensitivity to the cyclohexane conformation or equatorial/axial positions for RDR energies. We can see that the energy range for the N N bond dissociation energy is between 45 and 50 kcal mol 1 depending on the conformer. This value was calculated via summing the energies of an isolated nitro radical and the RDR conformer as compared to the AAE RDX conformer. We note that previous estimates, relying on B3LYP/6 311G(d,p) / / B3LYP/6 31G (d), predict only 30 kcal mol 1 37 This discrepancy highlights the inadequacy of parametrized mo dels, which cannot be expected to get the right answer for an arbitrary molecule outsi de the parametrized set. Electronic Excitation Spectra We present the CCSD EOM/6 311(2 + ,2 + )G(d,p) results in Table 2 5, as well as the CIS (CI with only single excitations in the CI framework) results (CIS is a biproduct of the EOM calculation and are reported for convenience). Calculations used 6 311(2+,2+)G(d,p). These oscillator strengths do not follow the Thomas Reiche Kuhn sum rule by virtue of the L=R t assumption discussed below; if values are not given for oscillator strengths, the matrix elem ent is too small to be considered. Note that the excitation energies are almost identical for all conformers. Hence, research on an RDX mechanism based on observing the electronic excitation spectrum would have difficulty in distinguishing which conformer is being observed. This is an example in which a computational study helps. The example is quite relevant, given the number of mechanistic studies done using electronic excitation energies for RDX. 18


34 Table 2 5 EE EOM CCSD and CIS Excitation Energies, eV Conformer AAA CIS EE EOM CCSD EOM Oscillator Strength 4.92 5.44 0.009 4.93 5.47 0.009 5.52 6.09 0.024 AAE 3.34 5.38 0.008 4.88 4.60 0.001 4.89 5.37 0.013 4.98 5.39 0.007 5.44 6.07 0.140 5.46 Boat 5.12 5.48 0.018 5.14 5.61 0.004 5.56 5.69 EEA 4.97 5.45 0.037 5.06 5.46 0.010 5.53 5.47 0.007 5.56 6.24 0.229 Twist 3.45 4.57 0.002 3.47 4.62 0.001 3.58 5.39 0.003 4.91 5.45 0.009 4.94 5.53 0.006 5.09 6.01 0.083 5.46 5.53 It must be mentioned that an approximate left hand eigenvector was used to determine the EOM CCSD density and, consequently, the oscillator strengths for each root. We used the nonorthogonal {R k } instead of the true {L k }. In the limit of using the full T operator, R k = L k such that this approximation also assesses the quality of the wave function. For a variety of test systems, this procedure turns out to be an excellent approximation. A study on three different molecules was performed (C 3v ammonia, C 1 ammonia, and hydrazine), takin g the X, Y, and Z components of the oscillator strength for each root using the aforementioned approximation. Only nonzero values were taken into consideration. The average signed error is 0.0016, and the average unsigned error


35 is only 0.0020 atomic units, with a standard deviation of 0.0077. Therefore, the oscillator strengths are indeed reliable, up to the level of theory implemented.


36 CHAPTER 3 HMX CONFORMERS AND CRYSTAL STRUCTURES Relevant Literature The literature of HMX is inherently also tied to the literature of RDX and CL 20; we shall endeavor to focus on what has been specifically proven for HMX and the most key points about RDX and CL 20 pertinent to HMX. Lea concluded that liquid phase RDX and HMX follow first order kinetics with activation ener gies of 47.5 and 52.7 kcal/mol, respectively 1 0 Cosgr ove and Owens disagreed with that analysis of Lea for RDX and instead thought that it was merely the vapor of the solid or liquid which began the reaction, and thus produces enough energy to begin a liq uid or solid phase reaction 11 Rauch and Fanelli, again working with RDX, believed that the mechanisms were fundamentally different between gas and liquid phases, above and beyond just differing activation energies 66 From these original papers, great co ntroversy arose as to what the mechanisms were, as noted in the above listed kinetics study references. Crystal structures of HMX were first reported in 1950 42 ; more crystalline polymorphs were discovered over time, thereby showing the immense diversity of structures possible from this flexible cyclooctane in the solid phase 23, 43 45, 64, 65 It has been observed that the crystal form has 4 coplanar carbon atoms with a C 2 axis; HMX has a C 2 axis; HMX has a center of inversion symmetry 45 The stabili ties of the crystalline forms at room temperature are the same as the order of the densities 6 5 In addition to collecting Raman spectra of the polymorphs, Brill and Goetz established that the crystalline form has a different ring structure comp ared to the forms. They believed that the form was more chair like and the are more boat like, as well as observing greater anharmonicity in the form 23


37 The first significant computational study on gas phase HMX found four conformers and estimated the energies via MP2/6 311G(d,p) 35 The first computational reference to the fact that there are indeed multiple boats or chairs in cyclooctanes 36 Addition ally, this study relied on energies from KS DFT functionals known to be less accurate than the energy differences between the two conformers in their study 117 121, 136 This approach is also problematic given that the study of nitrami ne geometries has sho wn that tho se functionals fail to p redict nitramine lattice parameters accurately 5 Such failure is a nother good example of how a parameterized method can fail outside of its parameterized scope, despite stellar success in other arenas 4 7 A similar kinetic study was conducted using many body methodologies again known to be insufficiently accurate for comp aring energy levels or kinetics 137 albeit thorough and effective in finding relevant geometric structures. There are many reactions in which smal l differences in geometry are the difference between exergonicity and endergonicity 116 The most thorough structural quantum chemistry study was performed by Lippert et al 138 That work closely examined the relationship between gas phase conformers and s olid state structures. It included careful examination of lengths and angles in crystals vs. KS DFT ge ometries. However, it considered only two types of boats and one type of chair. It is highly plausible that there are multiple matches to each crystal structure, given the sheer number of cyclooctane conformers. We have to walk before we can run. No one has the ability to figure out the mechanism rigorously in condensed phases. Interpreting the experiments is hard and different experiments contradict one another. What we can offer is an unambiguous


38 description of what is going on in the gas phase at STP. Once established, experim ental groups can try to see which methods match to this simpler system (gas, STP) to help sort out whose experiment is correct f or the condensed phase (assuming it starts from the condensed phase at all; some of the literature argues that vapor from the condensed phase IS the starting point of mechanism). Additionally, a computational estimate of condensed phase mechanism would ben efit from a comparison to a closely similar system for corroboration (the gas phase). Computational Methodology We need to approach chemical accurac y as best we can in this work. Refinement of our work for chemical accuracy will be the subject of future co mputational efforts. General purpose, non empirically parameterized functionals, such as PBE 143 LDA 144 VWN 145 146 typically do not offer accuracy as high as 3 kcal/mol for an arbitrary molecule, despite their more physically grounded quality c ompared to highly parameterized functionals 147 151 Nitramines are a somewhat exotic functional group, and thus we cannot be sure which functional to trust. Consequently, we opt for ab initio coupled cluster singles, doubles, and perturbative triples (CC SD(T)) 123 125 152 154 in a triple zeta basis set. Second order many body perturbation theory (MBPT(2)), or MP2 if one uses a canonical Hartree Fock reference 155 provides bond length geometries as accurate as 0.01 for single reference systems with a tr iple zeta basis set 1 26 One of the best software for massively parallelized calculations using MBPT(2) and CCSD(T) is ACES III 102 and was thus used. A basis set with some diffuse functions as needed, but not an excessive amount, is in the augmented Dunn ing style 104 156 We thus performed our geometry optimizations using MBPT(2)/6 311++G(d,p) 157, 158 as well as for the calculation of harmonic vibrational frequencies to establish energetic maxima vs.


39 minima. Single point energies were performed using CCSD(T)/cc pVTZ 104 We wish to perform energy extrapolation eventually, and thus choose the Dunning basis sets. However, given the aforementioned excess of diffuse functions present, we opt for the non augmented form. Searching Over the Potential Energy Surface Having stated how we wish to optimize and understand the energy differences, we must establish what to optimize and study. There has been great success in predicting the possible RDX conformers in terms of basic cyclohexane chemistry 3 72, 73 Ana logously, therefore look to cyclooctane chemistry for the basic HMX conformers. Previous work on cyclooctane chemistry shows that there are 10 conformers: boat chair, twist boat chair, crown, chair chair, twist chair chair, boat boat, twist boat, boat, ch air, and twist chair. This nomenclature of cyclooctanes derives from looking at each side of 6 carbons on either end and concluding whether it is a boat, chair, twist boat, or twist chair from cyclohexane chemistry. If the two sides have a perpendicular mirror plane and are either boat boat at either end or chair chair at either end we call it boat boat or chair chair. If such a mirror plane does not exist perpendicular to the ring, it is simply a chair or a boat. An exception to this is the crown conf ormer, which is a distinct concept of cyclooctane chemistry. Together, the two ends give the two part name (except the crown). For visual aid, in Figure s 3 1 through 3 4 we have examples of the HMX forms of boat boat, chair, chair chair, and crown.


40 F igure 3 1 Boat Boat Figure 3 2 Chair


41 Figure 3 3 Chair Chair Figure 3 4 Crow n We constructed each of these conformers out of the heterocyclic cyclooctane of HMX (accounting for all the additional conformers possible because nitrogens in the ring introduce more variations). Beyond this, we constructed every possible combination of axial/equatorial positions of the NO 2 coming off of the ring nitrogens for each of the ring conformers. From this large set of possible conformers, we calculated the op timal geometry via Newton Raphson optimization 107 Many optimizations went nowhere, with large RMS values persisting (greater than 1 mHart/Bohr persisting over at least 10 cycles, with no rea l downward trend in value). Tho se conformers were thus discarde d as un viable. The default convergence criterion was 0.1 mHart/Bohr. Many others simply converged toward one of the other stationary states found. The remaining conformers were subject to an energy Hessian with respect to normal coordinates to determin e if the conformer was a minimum or maximum. All of these were then subject


42 to a single point calculation of energy using CCSD(T)/cc pVTZ. We believe this technique to have been a rather exhaustive search of the PES of HMX. Resulting Conformers and Nomen clature There were 13 stationary states discovered from this analysis; they are given in Table 3 1. Table 3 1 Stationary States of the HMX Molecule Conformer Electronic E (kcal/mol) Smallest 2 Frequencies (cm 1 ) ZPE (kcal/mol) Electronic E + ZPE (kcal/mol) 1,7 diaxial 3,5 dipseudo boat chair 0 +28, +61 122.5 0 1,5 diaxial 3,7 diequatorial chair 0.7 +26, +61 122.3 0.5 2,4,8 tripseudo 6 axial twist boat chair 1.1 +21, +52 122.3 0.9 1,3,5,7 tetraaxial twist boat 1.2 +64, +73 122.4 1.2 2,4 diaxial 6,8 dipseudo boat chair 2.1 9, +50 122.2 1.8 2,8 dipseudo 4,6 diaxial boat chair 2.2 +16, +41 122.3 2.01 1,3 diaxial 5,7 diequatorial boat boat 2.3 +29, +37 122.2 2.1 1,5 diaxial 3,7 diequatorial crown 4.2 +39, +50 122.1 3.9 2,6,8 tripseudo 4 axial boat chair 5.7 43, 22 121.8 5.01 2,4,6,8 tetrapseudo twist boat 11.9 +26, +43 122.6 12.01 1,5 diaxial 3,7 diequatorial chair chair 13.5 9, 7 121.5 12.6 2,4 diequatorial 6,8 dipseudo twist boat 15.1 22, 9 122.3 15.01 CCSD( T)/cc pVTZ Energies, MBPT(2)/6 311++G(d,p) Harmonic Frequencies We shall explain our proposed nomenclature of these states. Clear nomenclature will help future efforts between solid state studies of crystals and gas phase conformers. The bridgehead atom


43 numbers to the ring atoms. We then classify each nitro group as either axial, equatorial, indicate a nitro group whose posi tion is halfway between being a full axial or a full equatorial. As is seen in Table 3 1, this is common, a result of the conformational flexibility the cyclooctane ring. Finally, we add the basic cyclooctane conformer label of the ring (crown, chair cha ir, etc.). The activation barriers in RDX are very small indeed, based on research of Vladimiroff and Rice 112 RDX chemistry is based on cyclohexane. Cyclohexanes are more rigid than cyclooctanes. Q ualitatively, we would thus expect that the barriers in HMX would also be small, if not smaller. The barriers to conformer interconversion are themselves not relevant to the detonation mechanism. The energetic barriers we calculate in a mechanistic study, however, may be off significantly if we start from the wrong conformer. For example, even within RDX conformers, N N homolysis varies by about 5 kcal/mol depending on the conformation 72 The N N homolysis would presumably be less sensitive, mechanisticall y speaking, compared to say the HONO elimination m echanism (which requires steric relationships to be just right). This highlights the importance of the right conformer as a starting point. The energy levels of HMX feature several very closely lying states beyond the resolution of CCSD(T)/cc pVTZ to guarantee. We are aware of this, and thus emphasize that the en ergetic ordering of states listed in Table 3 1 is NOT meant to be taken as the true ordering of states except when the differences are greate r tha n 2 kcal/mol Future work performing basis set extrapolation will hopefully resolve this issue, as well as taking into account the anharmonicity.


44 Four of these states are transition states / saddle points according to MBPT(2)/6 311++G(d,p) harmonic freque ncies. We list the first two harmonic vibrations in each case due to the presence of very small harmonic frequencies. Our Hessian calculations are done analytically, not requiring a grid as KS DFT generally does; this is a significant adv antage for small frequencies (the calculations do not suffer issues of numerical instability in the same manner). We emphasize that the absolute values of MBPT(2) frequencies are not generally more accurate in magnitude than previous studies done using B3LYP; only in cas es of numerical sensitivity to small values is MBPT(2) more reliable 108,109 However, small differences in energy curvature may or may not be accurately resolved with the MBPT(2) estimate of the correlation energy (or anything less than CCSD(T), for that matter). The PES curvature is quite small. Beyond just the basis set and correlation, anharmonic effects are important for such small frequencies. The presence of two negative frequencies makes it seem more plausible for a stationary state to be a genui ne saddle point or at least transition state, as is the case for 2,4 diequatorial 6,8 dipseudo twist boat and 2,6,8 tripseudo 4 axial boat chair. It is less clear whether 2,4 diaxial 6,8 dipseudo boat chair and 1,5 diaxial 3,7 diequatorial chair chair are genuine transition states/saddle points. This is especially true so far for the boat chair given that the next harmonic vibrational frequency is solidly positive (+50cm 1 ). Further refinement s of basis set, correlation, and anharmonicity are needed. Fo r the states, given our agnosticism about the signs of some o f the harmonic frequencies. However, clearly, if a geometry is a saddle point or energetic maximum, it is no t a true


45 One may notice that there are, in fact, only 12 states listed in Table 3 1, despite the claim of 13 distinct stationary states. There is no contradiction; rather, in our conformer searches, we came across an enantiomeric pair, and th us only list 2,8 dipseudo 4,6 diaxial b oat chair in Table 3 1's physical properties. The energy levels of the conf ormers show that most of them are energetically accessible at STP, excepting 1,5 diaxial 3,7 diequatorial chair chair, 2,6,8 tripseudo 4 axial boat chair, 2,4 diequatorial 6,8 dipseudo twist boat, and 2,4,6,8 tetrapseudo twist boat. As noted earlier, the procedure of CCSD(T) typically accounts for 99% of the correlation energies in a single reference system in a triple basis. T he only conseq uential approximation is thus the number of functions in the expansion of basis vectors 110,111 The use of our triple basis set makes it plausible that 1,5 diaxial 3,7 diequatorial chair might yet go lower in energy to be in the thermally populated range. This possibility is of importance, given that certain proposed mechanisms (such as HONO elimination) are likely to be sterically sensitive. This disproves the previous assertions that there are on ly 2 gas phase conformers 36 ; there are certainly more than two. We now proceed to correlate our gas phase conformers with known solid state results. Some conformers feature significant symmetry; this is important for the sake of matching to observed exper imental spectra. The 1,3,5,7 tetraaxial twist boat is approximately S 4 in point group. This is of particular interest, in that no HMX species with this symmetry has been observed experimentally; the nitramine orientation is unlike that of or 47 Moreover, the energy of this conformer is quite low with stable harmonic frequencies; one would expect its observation. It is conceivable that rapid


46 gas phase interconversion of different forms makes the observation of this distinct species difficult; suc h behavior has been shown for RDX 112 We nonetheless would be interested to see if experimental i dentification is yet possible. The 1,5 diaxial 3,7 diequatorial crown features C 2v symmetry. The literature has sta ted that both the and phases has C 2v s ymmetry 36,47, 51 137 when the gas phase conformer is optimized starting from either crystal monomer or crystal monomer. However, all of these references describe the and conformers as being types of boats. This is probably due to a fast and loose comparison to cyclohexane chemistry rather than cyclooctane chemistry. When one examines the gas phase version of the and crystal monomers, one sees the same nitramine orientation as is present in the 1,5 diaxial 3,7 diequatorial crown. The nomencla ture used in describing HMX should change to reflect this fact (cyclohexane chemistry should not be used to describe cyclooctane chemistry; this is a crown, not a boat). There is general agreement between the geometry of our proposed crown s tructure and t he crystalline and forms 51 We confirm the conclusion of Brand and coworkers 51 that the most logical geometric assignment of lengths/angles is such that the C 2v structure corresponds to the and crystalline monomers. W e thus do not lengthily repeat our geometric data confirmation. We do submit that our geometries are more trustworthy by nature of MBPT(2)/6 311++G(d,p) calculations compared to B3LYP/6 31G(d) calculations, but the chemistry conclusion is the same. Given that the crystalline forms of and distort monomers down to C 2 and C 1 respectively, we agree with the inference the crystal field effect is strong in these cases. The monomer from the is known to have C i symmetry, consistent with the 1,5 diaxia l 3,7 diequatorial chair form. T he bond lengths


47 also very closely match the crystalline geometry 43 This fact also speaks to the relatively weak intermole cular forces of the crystal in that the gas phase geometries match the crystal to within 0.01 It is interesting to see such contrast in the strength of intermolecular forces of organic crystals of the same molecule. It is known that C 2 axes are present in the crystalline form 43 The 1,5 diaxial 3,7 diequatorial chair chair form features a C 2 axis of symmetry, but the N N bond lengths do not match well, hence we defer to our previous conclusion that the crown C 2v form best matches the form. All of these conformers would be of interest in a solid state calculation. The potential exists that these conformers, as starting geometries of a lattice vector optimization, may yield new crystalline forms. Even if not, it would be of interest to observe how the intermolecular forces mold the gas phase conformers strongly in some cases, weakly in othe rs, to the solid state structure. Indeed, given how many RDX crystalline forms exist from a simple cyclohexane, it is surprising that only four HMX crystal forms have not yet been found. Previous work exists in the literature for the opposite process: st arting from the known crystal structures, new gas phase conformers were discovered at the time 51 It is a reasonable hope that the process might be fruitful in reverse.


48 CHAPTER 4 CL 20 CONFORMERS AND CRYSTAL STRUCTURES R elevant Literature Three B3LYP computational studies have been performed on gas phase conformers 159 161 The first such paper identified four conformers 159 As just discussed, i n chemically similar compounds, the harmonic vibrati onal frequencies are quite small 72 It would be u nwise to use a numerical grid for gradients or hessians, as is required by KS DFT. Numerical grids play havoc with the sign of harmonic vibrational frequencies, and the correct sign is needed to determine energetic minima vs. maxima. As such, the conform ers interpreted to be minima in that paper 159 cannot be known rigorously to be minima 132, 133 The other two papers chose one particular CL 20 conformer at random, without thought to the other conformers, and with a computational methodology insufficient to be trustworthy a priori One of tho se papers attempted to study the mechanism of CL 20 from an arbitrary conformer 160 This seems unwise, given the sensitivity of kinetic studies to small energy differences, as well as reliance on B3LYP/6 31G(d,p) ene rgies to be accurate to within 1 kcal/mol, which is known to be untrue 117 121, 136 As such, we enter the fray to answer questions of the conformers of CL 20 using high quality methodologies from the bottom to the top as with RDX and HMX. Thus we use par ameter less, ab initio methods in the final step of our analysis to confirm our studies, and experimental checks to confirm the validity of specific properties calculated (like crystal polymorphs, lattice constants, etc.) when using highly parameterized me thods. The value of using force fields in early steps, because of their computational cheapness, is not to be underestimated, but must be checked stringently in subsequent


49 steps. From such methodologies, we may know the potential energy surface of all co nformers of CL 20, and thereby predict new crystal structures as well as predict the chemical mechanism of detonation. M ethodology We began by assembling a structure of the CL 20 molecule from a rudimentary guess of common bond lengths and angles. This fo rm was then optimized using second order many body perturbation theory (MBPT(2))/6 311G(d,p)) 122,154,157, 158 Keep in mind that Many body perturbation theory is Mller Plesset perturbation theory (MP) if the reference is chosen to be Hartree Fock. Howeve r there are many instances in which the use of non Hartree Fock references is advantageous to circum navigate issues of degeneracy. A s such, we keep an open mind to the general formulation of applying perturbations atop an arbitrary reference state. We chose a triple zeta basis set, as MBPT(2) is known to be accurate to within +/ 0.01 for a single reference molecules 126 We emphasiz e the value of choosing perturbation theory. Because it is an ab initio method ( ab initio here meaning a parameter less methodology which is systematically correct able), we can expect it to work equally well on any single reference molecule. This is in contrast with KS DFT 139 142 for which the results can be highly dependent on the choice of functional. Again, w e use d the ACES III software for our calculations, as it is optimally designed to be a massively parallelized software for calculations using M BPT(2) 102 From the optimized molecular geometry, we began a Monte Carlo force field based search for other CL 20 conformers. Specifically, with respect to the initial CL 20 conformer, we varied bond lengths, angles, and dihedrals in a Metropolis Monte C arlo fashion using the PCModel92 software and the MMX force field. The standard MMX


50 force field does not have parameters for the nitramine functi onal group, as it is somewhat exotic for standard organic chemistry. We substituted an amine grou p in place o f the nitro group. F or the sake of finding a rough estimate of the conformers, the amino group is similar in size and bond length, such that the sterics of what is possible is preserved. This is a trivial assumption for the conformational search All as pects of the search were left at default parameters except for the following. The duration of the searches is indefinite until either 10,000 minima structures are found or a minimum structure is found 500 times. The search was executed 3 times from vario us random seeds. We allow ed for energy ranges of 3.5, 20, 50, 100, and 400 kcal/mol energy increases for possible generated structures to be minimized relative to the conformational mini mum. The search never generated more than 14 possible conformers, wh ich generally only differ in the axial/equatorial position of the amines. This seems reasonable; the caged structure of CL 20 (see below) is such that few relaxations seem possible. These 14 conformers were subject ed to a MBPT(2)/6 311G(d,p) optimization of geometry using a Newton Raphson algorithm 107 Of the 14, only 6 converged toward an ener getic extremum. The other 8 were discarded from further analysis; their geometric optimizations did not converge ( the RMS force remains quite high ) The 6 which d id converge did so relatively quickly for a molecule of this size. The harmo nic vibrational frequencies then were calculated using the same many body methodology. Five are energetic minima, the sixth a maximum. The energy levels then were considered fro m these MBPT(2)/6 311G(d,p) calculations.


51 T he Six Conformers The structure of CL 20 (shown in Figure 4 1 ) is that of two cyclopentanes bridged together by one carbon in each of the pentagons, in addition to which th e bases of the pentagons constitute a cyclohexane to form this bridged bicyclic compound. Figure 4 1 Side View of Equatorial/Axial CL 20 Figure 4 2 Side View of Diaxial CL 20.


52 There are two major differences that distinguish the CL 20 conformers from one another: the axial vs. equatorial nature of the cyclohexane nitramine groups, and the angle of the cyclopentane nitramine groups with respect to one another. Figure 4 2 shows t he differences in axial vs. equatorial positions in the bridging cyclohexane base. The below Figure 4 3 illustrates the differences in how the cyclopentane nitramines orient with respect to one another. Figure 4 3 Figur e 4 4 Conformer; axial, equatorial. See text for notation


53 Figure 4 5 Conformer; di equatorial. Figure 4 6 Figure 4 7 equatorial


54 Figure 4 8 Transition State Conformer; diaxial The conformers are so named based on CL 20 crystal structures The actual crystal structure data have varying public availability 99 The basic lattice constants and other bulk properties are available for four of the five polymorphs ( but not ) 99 162 166 Refined data, such as molecular level data, however, are not available for the anhydrous crystals or crystals as a matter of restricted information 99 162 165 The hydrated crystal data is incomplete 99 The and crystal structures, including molecular level properties, are public knowledge and well known 166 The form exists only at high pressures, and no data exists beyond its basic zeta crystalline assignment 163 From the myriad of sources referred to above, we have endeavored to establis h which of our conformers belong to which crystalline polymorph, if any. The first conformer has the cyclopentane nitramine orientation and axial/equatorial combination consistent with the crystal, and is thus named the conformer 99 The and crys tals are known to have a cyclopentane nitramine orientation consistent with both of the next two conformers 99, 166 The crystal has a form of CL 20 which is diequatorial, indicating that the third conformer ought to be designated the conformer 99 166 The crystal structure is not so resolved that we can claim that the


55 second conformer is the basis for the crystal structure; nonetheless, it is known to be different from the but still have the same cyclopentane nitramine orientation 99 We thus p ropose that the second conformer based on our evidence is likely the form. To emphasize that this is a reasonable interpretation of data but not experimentally confirmed, we call the second conformer It would be of great interest for our experiment al colleagues to investigate this conjecture. The conformers differ with respect to one another only in that has an axial and an equatorial pair of nitramines on the cyclohexane, whereas the conformer has two equatorial nitramines. The next p air of conformers, and again only differ in cyclohexane nitramines; the conformer has one axial, one equatorial, conformer being diequatorial. The conformer matches the crystal structure in nitramine cyclopentane orientation and diequatoria l nature 18,49 The conformer is of C s symmetry. That the conformer has a relationship to the crystal is very speculative; we offer this purely based on the fact that there is one remaining crystal structure not accounted for. It is more than possible that the crystal basis is one of the other conformers again or a completely unstable conformation in the gas phase alone (for example, the RDX crystal is actually the transition state between the diaxial, equatorial to triple axial conformatio ns, completely unstable in the gas phase on its own) 6 However, no nomenclature has existed before th is point in previous CL 20 work. W e offer this first attempt at taxonomy to help organize previous data in terms of our own more complete survey of the c onformational landscap e of CL 20, as well as to relate our data better to the solid state community. The transition state conformer is so named because it is the only one that is an energetic maximum via harmonic frequency analysis. It is also of C s symm etry.


56 We emphasize the success of our methodology for finding the molecular bases of the crystal line polymorphs Our Monte Carlo force field search combined with MBPT(2) discriminating against the false positives has produced 5 minima, in accordance with 5 crystal structures. It is suggestive (albeit not deterministic) that our methodology is efficient if we wish to apply it to similar chemicals to find the relevant structures. We also emphasize the extent to which these geometries are a true benchmark f or the current best estimate of the gas phase geometries of these conformers; the only method expected to produce a better geometry would be CCSD(T), or higher coupled cluster excitations 108, 123,124, 128,141, 152 154 such as presented in the preceding cha pters. Energetics and Geometric Speculation Therein The energetics of the compounds are displayed in Table 4 1. Table 4 1 Energetics of the CL 20 Conformers kJ/mol eV Conformer 0 0 Conformer 16.4 0.17 Conformer 17.1 0.18 Conformer 10.6 0.11 Conformer 23.1 0.24 Transition State 23.3 0.24 MBPT(2)/6 311G(d,p). All energies are reported to one decimal in kJ/mol, two in eV for the sake of qualitative comparison. We emphatically do not claim that this method is accurate to this many decimals in energy. The conformer is of lowest energy relative to the methodology used, but we adaman tly make no statement as to which conformer is the minimum. The only methodology capable of being accurate to


57 within 0.1eV consistently, a priori is pr obably CCSD(T) with an extrapolation of single particle basis sets or use of an explicit R 12 formalism 141 We had originally intended to perform CCSD(T)/cc pVTZ calculations on these systems, but changed course in the face of the foregoing evidence. Shor t of performing a basis extrapolation, no real insight would be gained about the reaction or conf ormer populations from just those calculation alone. We may say that the energy levels a re near degenerate, but it would be wildly speculative to read even th e qualitative ordering as being correct when the energy levels are this close. Rather, we can see from the data that all of the minima may be reasonably populated at STP and the reaction mechanism might proceed from any of these starting points. In order to consider the kinetics accurately, one must consider the activation barriers from any of these starting points. Note that this point is in sharp disagreement with previous kinetic studies 161 Rather, given that the energies are not quantitative enough to compare against one another, we attempted to find a geometric, intuitive means to judge the conformers bond lengths to 0.001; a clear trend mig ht be seen easily if manifest at 0.01. In Table 4 2, we list the bond lengths of various parts of the CL 20 molecule. We consider the C C, C N, and N N lengths of various parts of the hexagons and pentagons in the bicyclic structure. If a given conformer did not have a cer tain type of nitramine, we cyclopentanes which joins together the rings. The table lists the conformers in order of increasing energy from left to right, for ease of examining a tre nd within the lengths.


58 However, no clear trend emerges There are nonetheless clear trends in the comparative bond lengths of CL 20 vs. other nitramines. Table 4 2 Bond Lengths of CL 20 Conformers Conformer Transition N N Axial 1.43 1.43 1.43 N/A N/A 1.42 N N Equ. 1.44 1.44 1.45 1.44, 1.43 1.43 N/A N N Penta 1.43 1.43 1.44 1.45, 1.40 1.43 1.41 C C Hexa 1.57 1.57 1.57 1.57 1.59 1.57 C N Penta 1.48 1.47 1.45 1.47 1.47 1.46 C C Bridge 1.57 1.57 1.57 1.57 1.58 1.57 The two major qualities that distinguish the CL 20 conformers are the axial/equatorial combination and the cyclopentane nitramin e orientation. We considered whether equatorial is inherently more favorable in this caged structure than is axial, and vice versa. This com parison also shows no clear trend in stability, as there are mostly axial/equatorial mixtures, rather than diequatorial or diaxial. The orientation of the nitramines does not, to our analysis, lend itself to an obvious trend. We thus infer that no simple geometric measure can give guidance on the energy minimum. G iven the small energy differences in the conformers (0.2eV), this is not surprising. Implications Toward Reactivity The harmonic vibrational frequencies of the CL 20 conformers lend credence to the idea that homolytic dissociation of the N N bond in the nitramine is a major mechanism of decomposition. The harmonic vibrational frequencies are given in Table 4 3. We do not propose that the accuracy of the harmonic frequencies is within a wavenu mber; the absolute values of MP2 frequencies are not very good compared to an IR spectrum 108


59 Table 4 3 Lowest 4 Harmonic Vibrational Frequencies Calculated using MBPT(2)/6 311G(d,p); Conformer Transition 1 48 48 48 49 46 41 i 2 55 57 54 54 52 44 3 59 59 56 57 59 53 4 68 61 61 64 61 56 However, with an appropriate scaling factor, they do perform well 108 For our purposes, though, a qualitative insight into the mechanism is reliable from this methodology. We can see that the harmonic frequencies for all of the conformers are comparable, and all quite small. Visualization of these normal mo des show that all of them are various types of nitramine wags. Thus, the weakest bond is indeed the nitr amine N N bond. In particular, it is a nitramine wag on the cyclohexane nitramine. This implies that homolytic dissociation does not proceed from the cyclopentane nitramines, and considerably reduces the number of possible reaction paths to consider. Th e zero point energies are all 141 kcal/mol; the zero point energies only vary among con formers by up to 0.6 kcal/mol. We also compare the CL 20 N N length of RDX for the equatorial position is 1.41 72 ; for CL 20, it is 1.44. The CL 20 N length of RDX for the axial position is 1.44 72 ; for CL 20, it is 1.42. The CL 20 axial bond is shorter, and therefore stronger than in RDX. We know that CL 20's rate constant in solution phase decomposition is 3.5 greater than RDX 34 Assuming the reactivity trend is qualitatively true in the gas phase( not unreasonable), something must account for this difference in reactivity. The putative mechanism considered in


60 nitramine degredation is nitro group homolysis 1 10 If the CL 20 nitramine dissociation starts from the cyclohexane nitramines, we would infer that it is probably the equatorial nitramine bond that is breaking, given that the equatorial bond is weaker i n CL 20 than the axial.


61 CHAPTER 5 DEA CCSD CONQUERING PATHOLOGICAL MOLECULES The Multi reference Problem In this chapter, the singlet triplet splitting s of proto type diradicals methylene (CH 2 ), trimethylene methane (TMM), or th a, meta, and para benzynes, are c omputed with DEA methods. Such methods are commonly used to assess the performance of multi reference (MR) methods, and yet require more careful study unto themselves as a benchmark. All these diradical share the characte ristics of 2 hole 2 particle MR problem. Assuming that the underlying n 2 electron state is correctly described by single reference coupled cluster theory, the DIP and DEA methods by construction include all four determinants of a 2 hole and 2 particle pro blem of n electron system in an un biased manner. That is a prerequisite for the correct description of the multi determinantal n particle system (details below) Furthermore, instead of relying on widely used published geometries that were obtained at var ious levels of theories, a consistent set of geometries at the CCSD(T) level are obtained and used for all the single point energy calculations. Additionally, the SR MBPT(2), SR CC SD, SR also were computed for completeness. Multi refe rence Coupled Cluster Theory in Context The treatment of classes of multi reference coupled cluster (MR CC) problems can be achieved by exploiting single determinant reference states composed of a different number of particles from that for the MR target s tate, with subsequent addition or removal of excess electrons to define the N particle states of interest. The idea is not new, going back to the work by Nooijen and Bartlett on using the double ionization potential similarity transformed equation of moti on (DIP STEOM CC) and its


62 complementary double electron attached variant (DEA STEOM CC), to treat the vibrational frequencies of ozone and certain states of NO dimer. But such an equation of motion coupled cluster method for adding (DEA) or removing (DIP) two electrons from an underlying single reference CC method for an problem, with extension to three (TIP/TEA) or four (QIP/QEA), has several attractive features that make the approach worth pursuing. The ansatz has a global extensive part, and a local correlation, intensive part. The wavefunction for the n parti cle state system is written as where is the CI like right hand eigenvector in EOM CC. This has the advantage that instead of asking a fully e xtensive MR CC method to account for all important effects, the intensive part allows one to target the usually local multi reference behavior. When based upon a closed shell reference for |N2>, the target states are automatically spin eigenfunctions. Th e guarantee of a spin eigenfunction as is used in CI is not always achievable in most MR CC methods; this, however, is an exception. DIP/DEA EOM CC wavefunction is operationally single reference making it as easy to apply as single reference CC, with no de cisions for the user but basis set, level of correlation, and a choice of the one spatial active orbital to doubly occupy in the N+2 vacuum, or to un occupy the N 2 vacuum. The DIP/DEA EOM CC is invariant to active orbital rotations, by virtue of the activ e orbitals being in either the occupied or the unoccupied space. Multiple states can be obtained from the EOM matrix diagonalization, providing excited states as well as the ground state, so in the event that other occupied orbitals in DIP or virtuals in DEA


63 interact strongly with the chosen active orbital, then that solution occurs as well, and will appear as one of the eigenvectors. This helps to confirm the particular orbitals that manifest MR character are correct, and if there are more than two, it m ight suggest a subsequent TIP/TEA or QIP/QEA calculation. The price paid for these attractive features is that each stage of the calculation has to be converged: first the SR CC solution for the closed shell system, then the DEA/DIP EOM solution itse lf. Because of orbital dependence in these calculations, the DIP solution in particular can sometimes be difficult to converge. This contribution addresses the use of DIP/DEA CCSD single triplet splitting in diradicals, which have been considered a mult i reference test problem, with many different types of methods applied. These include generalized SR CC, like flip, and hybrids like Reduced Multi Reference. There are also more complete state universal (SU) forms: the state specific Brillouin Wigner (BW) form, Mulherjee, or internally contracted form. The singlet states involved tend to manifest most of the multi reference character, since they can vary from a two determinant open shell singlet, for which there are other met hods available, to a GVB when a bond starts to form, or to situations where all four determinants, as shown below, could have a large role in the wavefunction. The triplet, on the other hand, can frequently be described well by standard ROHF or UHF SR CC using the high spin determinant as the reference.


64 Theory The wavefunction ansatz for the MR N particle state is (5 1) The solution is a single reference CC one, but for electrons while the (5 2) and (5 3) restores the N particle solutions, by either removing or adding two electrons to the underlying SR CC solution, (5 4) for the Then inserting this ansatz into the Schrodinger equation using th e usual EOM CC strategy, we obtain (5 5) which provides the N particle k states of interest via the EOM CC diagonalization. The quantity is defined from the CC wavefunction (5 6) (5 7) In two steps, DIP/DEA EOM CC provides solutions for different k states via matrix diagonalization with two orbital, two electron MR character. To involve more orbitals and electrons in a MR description requires the TIP/TEA and QIP/QEA extensions. Because is non Hermitian, there is also a left hand eigenvector,


65 = which, unlike the right hand eigenvector, is not connected to Both are required in the treatment of properties. But the treatment of properties in EOM CC is well known, and the same methods can be used for the DIP/DEA EOM CC solutions. The above procedure is MR for the target states. For example, from the |N+2> vacuum, where the HOMO and HOMO 1 spatial orbitals are assumed to be quasi degenerate, the operator generates all four of the determinants that a MR space would demand. As each of these determinants will be weighted by the coefficients in the CI like operator, their contribution to the final wavefunction is determined by matrix diagonal ization allowing any coefficient to be as large as required in the MR description. Figure 5 1 Reference Determinant and Achieved Static Correlation in the DEA CCSD Scheme. Analogously, a N 2 vacuum state would consist of orbital J to I and A in this fig ure, but the same four determinants would be generated by adding two electrons in all possible ways into the two quasi degenerate orbitals. In these methods the MR space is naturally introduced after doing a SRCC calculation for the problem. The size extensive SRCC result handles dynamic correlation and the EOM diagonalization properly introduces the intensive or local part of the correlated wavefunction to account for non dynamic or static correlation. The former pertains to the left right corr elation for


66 bond breaking subject to an incorrectly separating RHF reference function, while the latter means the kind of correlation required to get the multiplets of complicated open shell systems like transition metal atoms right, of which being a spin eigenfunction might be as important as correlation. Thus the DEA CC has a right hand operator eigenfunction property of this method, R, that is CI like. If we started from a more traditional MR viewpoint, we would add to the four determinants shown above typically single and double excitation to build a MR CI wavefunction. The role of the 3h1p term in Equation 5 3 does exactly that for single excitations of the two hole MR space. So we view the r ij as the coefficients of the four determinant MR space, c, w ith the (r) corresponding to additional single excitations (s) among these four determinants, such that R=C(1+S). In principle, from the CI part of the ansatz, we could add double and higher excitations, leading to the full CI limit. The expone ntial part of the ansatz is SR CC theory. The reason the distinction between the CI part and the CC part is pertinent is that, in DIP EOM CC, contributions to the wavefunction involve three hole operators like that most naturally would arise f rom having T 3 in the ground state. A response theory like EOM CC focuses on the ground state, with all other states being derived from them. But these three hole terms cannot appear for a CCSD ground state, while they would from the MR CI viewpoint. It turns out that these terms are extremely important numerically, much more so than are the reference state triples. This importantly allows us to bypass triple excitations in the CC part of the ansatz. It also suggests the alternative designation MR EOM D IP CC, with the MR space created by the DIP operator, but then single excitations are taken among those functions as they


67 would be in MR CI while the CC part of the problem can be restricted to just CCSD. So from the MR CI viewpoint, we bypass a need for double excitations and from the CC viewpoint, using triple and higher excitations. The DIP method has a fundamental weakness: it is completely inapplicable to systems for which the 2 anion does not exist. If there is no bound state for the doubly negativ ely charged species, obviously the method is flawed. For example, the species does not exist, and thus DIP may not be used on If one attempts the calculation with small basis sets, one may converge the SCF and then use the EOM framework to return to the methylene species and attain accurate answers; however, this only works for a small basis set. Given enough basis functions, the attached electrons in the SCF are not bound and enter the continuum. This could hypothetically be overcome wit h some type of artificial environment that allows for a bound state to exist in the SCF, and whose effects are somehow overcome or perturbed away in the correlated calculation. There is no doubt that the large effects of single excitations are mostly due t o the need for orbital relaxation following the introduction of the vacuum. Any orbitals can be used to build this vacuum state, so there is no unusual dependence on e specific orbital dependence in practice. Methodology The geometries are optimized at the UHF CCSD(T)/6 311++G(2d,2p) level of theory. The harmonic vibrational frequencies are computed to verify that the optimized structures are minima and to obtain the zero point energy using MBPT(2)/6


68 311++G(2d,2p). In the case of para benzyne, the symmetry is lowered from the D 2h to C 2v subgroup to obtain a symmetry broken reference state that has partial diradical character (unpaired electrons localized on carbon ato ms at para positions). Additionally, the para benzyne harmonic vibrational calculation was performed using CCSD(T)/6 311++G(2d,2p) due to the extreme multi reference challenge of para benzyne. The single point energy calculations employ cc pVDZ, cc pVTZ a nd cc pVQZ basis sets. results are obtained with the extrapolation scheme proposed by Helgaker 126 The optimized geometries, the vibrational frequencies, CCSD, CCSD(T) energies are obtained with the ACES II program while the DEA results are obtained with ACES III. Spherical basis functions are used. Where core polarization calculations were undertaken, the protocol was to find the relative energy difference between singlet/triplet using CCSD(T)/cc pCVQZ with and without dropping the core. The difference in the singlet/triplet gap using dropped core orbitals vs. not is what we define to be core polarization. M ethylene The singlet/triplet gap of prototype dira dical methylene has been the su bject of numerous studies 167 176 and a historically significant narrative of the interplay between theory and the experiment The methylene radical has a triplet 3 B 1 ground state, which is single reference in nature. The mol ecule features a C 2v point group. The electronic configuration of M s =0 component of 3 B 1 state is The corresponding singlet configuration giving 1 B 1 state is higher in energy than the 1 A 1 state (The state labeling correspond to the C 2v irreducible


69 representations assignment to the x, y, z directions as shown in Fig 5 1). The singlet 1 A 1 state is multi reference in nature, ideally represented by two ref erence determinants and is also C 2v The experimentally measured S T gap is between the 1 A 1 and 3 B 1 states. The most acc urate experiments 172, 177 measure the singlet triplet splitting of methylene to be 9.0 (9.050.06 and 8. 9980.014 kcal/mol respectively ). The computed zero point energy, Born Oppenheimer, relativistic, and core polariza tion corrections are 0.49 172 0.1 173 0.04 178 0.4 (done in this work as well as confirmed in previous estimates 167 ) kcal/mol respectivel y. These corrections are subtracted from the experimental analog that can be compared with the computed single point energies. While it would have been more appropriate to correct the computed values, however, since we are also interested in comparing with the other published MR results; it became necessary to follow what is established in the literature. nergy splitting is is different from the often quoted experimental value, 9.37 kcal/mol which neglect the co re polarization effects 171 Unlike prior singlet/triplet gap calculations, which relied on CISD/DZP geometries, we have chosen to obtain geometries for each state using CCSD(T)/6 311 ++G(2d,2p) 174 It has been shown that CCSD(T) geometries on single reference cases are accurate to within 0.01 126 Obviously this estimated error does not directly translate into the estimated errors for systems that have MR character and should be used with caution in that context. The 6 311++G(2d,2p) basis is approximately of the


70 same quality of cc pVTZ in practice, but has some diffuse functions. There is indeed some variation in the geometries of the singlet A comprehensive collection of previously published S T gaps are presented in Table 5 1 and our DEA results ob tained at CCSD(T)/6 311++G(2d,2p) geometry, in various basis sets, are shown in Table 5 2. Table 5 1 Methylene Singlet Triplet Splittings, kcal/mol 3 B 1 1 A 1 S T kcal/mol Basis= Other cc pVDZ cc pVTZ cc pVQZ cc pV5Z CCSD 1 CCSD 12.77 11.20 10.58 10.39 CCSD 1 RMR 12.04 10.42 9.84 9.74 SU CCSD 1 RMR 12.20 10.71 10.17 10.08 SU CCSD 1 SU CCSD 12.04 10.52 9.92 9.76 CCSD(T) f 1,2 CCSD(T) 12.11 10.43 9.78 9.58 CCSD(T) se 1,3 CCSD(T) 12.05 10.35 9.70 9.51 CCSD(T) f 1 RMR(T) 11.89 10.19 9.55 9.40 CCSD(T) se 1 RMR(T) 11.81 10.11 9.48 9.33 SU CCSD(T) f 1 RMR(T) 12.20 10.76 10.22 10.10 SU CCSD(T) se 1 RMR(T) 12.05 10.55 10.00 9.88 CR CC(2,3) 1 CR CC(2,3) 11.91 10.29 9.68 9.53 EOM EE CCSD 4 12.13 EOM EE CC(2,3) 4 10.47 EOM SF CCSD 4 11.76 EOM SF CC(2,3) 4 11 CCSDt 5 CCSDt 9.87 CCSD(T) h 5 CCSD(T) h 9.14 CC(t;3) 5 CC(t;3) 9.24 MR BWCCSD MR BWCCSD 10.52 9.90 9.72 MR BWCCSDT MR BWCCSDT 10.18 9.51 TD CCSD TD CCSD 9.20 3 B 1 1 A 1 S T kcal/mol CCSDT 5 CCSDT 9.52 FCI 6 FCI 11.14 Experiment 7 8.25


71 We emphasize that our methodology, unlike many in Table 5 2, is the same for both singlet and triplet electronic state. By contrast, certain multi reference CC methods use a composite approach, and thereby treat the singlet with a different method than the triplet, as is the case in RMR style methodologies. Numbers are listed to two decimals for computational comparison sake, and are not intend ed to imply accuracy to so many decimals. If only one method is listed for the two states, the same method is used twice. 1. Reference 171 ; geometries come from reference 17 4, a CISD/DZP estimate. 2. diagona l Fock matrix elements are included 3. diagonal Fock matrix elements 4. All of these calculations are done us ing the TZ2P basis ; the calculatio ns were performed by reference 169 with one frozen core orbital and one frozen virtual orbital. Geometries come from FCI/TZ2P with one frozen core orbital and one frozen virtual orbital from reference 16 8. Spin flip calculations use a UHF reference. 5. These calculations we re performed by reference 17 9. All values are done with a basis set extrapolation based on aug cc pVTZ and aug cc pVQZ, except for the CCSDT calculation which used only aug cc pVTZ. All calculations used spherical d, f, and g functions, as well as a froz en core. Geometries come from FCI/TZ2P, using reference 8, one frozen core and one frozen virtual. The active space used was (2,2), corresponding to the HOMO LUMO orbitals. 6. This FCI comes from reference 16 8 below, using the TZ2P basis set, dropping one c ore and freezing one virtual. 7. This is based on the experimental result of reference 1 72 with computational estimates of the non adiabaticity and relativistic effects subt racted out based on references 178 and 17 3, respectively. The zero point energy diff erences are accounted for in the experimental protocol. For the discussion purpose, we can classify the S T gaps shown in Table 5 1 into two categories. In the first category, the S T gap is obtained as a difference of the two


72 separate individual energy calculations for the singlet and triplet each. Those individual singlet and triplet calculations do not necessarily have to be by the same method, and in fact, the majority of the results in Table 5 1 use composite methods in whi ch the singlet and triplet treated by two different methods (primarily the singlet as a MR state while the triplet as SR). In the other category, the S T gap is directly computed from a common reference state. The EOM, SF EOM and our DEA results belong to this category. Let us first focus on the results in the first category. The CCSD treatment for both singlet and triplet is inferior to the use of RMR to describe the multi references singlet. Nonetheless, energy remains off of the experimental answer by 1.5 kcal/mol, despite basis saturation. The various forms of CCSD(T) describe the S T gap more accurately. The use of RMR(T) on the singlet with a form of CCSD(T) on the triplet incrementally approaches the correct answer better than RMR, SU CCSD, or CCS D(T) alone. Curiously, forms of SU CCSD(T) for the triplet with RMR(T) on the singlet seem to make matters worse. CR CC(2,3) describes the gap correctly without requiring an active space. All of these calculations nonetheless remain off by 1 kcal/mol even at basis saturation, ZPE effects, polarization effects, and relativistic effects. All errors noted in the use of SU CCSD(T) might be attributed to issues of the geometry used. In the second category of results, the EOM EE and EOM SF methods may have the error they do purely as a function of the small basis set used; in all cases, some of the error is also undoubtedly the result of the CISD/DZP geometry used. The CCSDT/aug cc pVTZ calculation of Piecuch is close to the experiment, and would probably march closer to the 8.25 kcal/mol mark with a larger basis. A well chosen active space using CCSDt unsurprisingly fairs as well as the MR forms of coupled cluster in Table 5 1. The MR


73 BW styles perform as well as the other MR CC styles, hovering in the high 9 kcal/mol range. The two determinant CCSD performs better than almost all results in Table 5 1. The FCI presumably does not approach the experimental answer well purely as a function of the small basis used. It would be of interest to attempt a FCI in a standard Dunning basis and not a random basis like TZ2P. This would then allow for a more fair comparison of FCI results with standard coupled cluster techniques. It is of particular evidence that the lambda based methodologies are inferior; they do not perform any better than CR CC techniques, despite both techniques attempts to be a multi reference single reference coupled cluster methodology. One might say that there is indeed no free lunch without opting for a CI like framework. Table 5 2 Methylene Singlet/Triplet gap at CCSD(T)/6 311++G(2d,2p) Geometry. cc pVDZ cc pVTZ cc pVQZ 3 4 Extrap RHF / UHF 28.7 28.5 30.1 28.2 MBPT(2) 17.7 15.4 14.3 13.5 CCSD 12.8 11.2 10.6 10.2 CCSD(T) 12.1 10.3 9.7 9.2 12.1 10.4 9.8 9.3 DEA CCSD 11.6 10.4 9.6 8.9 Experiment 8.25 Methylene Singlet/Triplet Splittings, kcal/mol. Numbers are listed to one decimal purely for computational comparison sake, and are not intended to imply accuracy to so many decimals. Performed by the authors using CCSD(T)/6 311++G(2d,2p) for geometries and frequencies. Core orbitals were dropped for energy calculations. The experimental value has the same reference as in Table 5 1. As stated earlier, the DEA is direct method in the sense that it does not require t wo separate calculations. It is easy to use and there is no need to choose an active space. The valence space basis set extrapolation of our DEA CCSD value is 8.9


74 kcal/mol. Despite of the fact that there is not explicit inclusion of triples in the DEA CCS D, the results seems to be better agreement with experiment than the methods that include triples explicitly. The fact that DEA CCSD performs so well speaks to the accuracy of the methods The effect of geometry, in this case, is small; our CCSD value diff value). Our geometry using CCSD(T) has a trivial difference in energy (0.1 kcal/mol). D(T). Trimethylenemthane The singlet/triplet gap of TMM has been the su bject of numerous studies and a historically significant narrative of the interplay between theory and the experiment 169 171, 175, 179 182 The relevant singlet triplet gap for TMM is between the 1 A 1 singlet state and the 3 A 2 triplet state. This singlet state is C 2v and planar. This 3 A 2 may also be referred to as 3 B 2 ; the 3 B 2 state is specification of state in terms of its most symmetric Abelian group C 2v ; in reality, it is a D 3h molecule of electronic state 3 A 2 We specify what our Cartesian axes are; there has been confusion/contradiction/ambiguity in the literature in spectroscopic designation without this information. In Figures 5 1 and 5 2, we show our axes for the 1 A 1 stru cture and 3 A 2 (D 3h ) / 3 B 2 (C 2v ) structures, as these are the relevant states for the experimentally observed singlet/triplet gap. In Tables 5 3 and 5 4, we compare our CCSD(T)/6 311++G(2d,2p) geometries to other works. All lengths and angles are specified in terms of the atom numbering scheme in Figure 5 2. All dihedrals are 90 degrees by the symmetry of the system in question; all angles are 120 degrees by symmetry (we note that some authors have allowed a small relaxation from the symmetry governed angl e, which never relaxed by more than 1 degree). References are same as given in the singlet geometry


75 Figure 5 2 1 A 1 C 2v structure (a), general orbital diagram (b), and relationship among the TMM states (c). the occupation is represented by either | + > or | ->. We emphasize that we use the canonical choice of axes consistent with spectroscopy texts. The primary axis is always the z axis, perpendicular to the plane of the molecule. The x axis, in C 2v is the B 1 irreducible representation. We emphasize this choice for clarity, because previous pionee ring research on these molecules have chosen non standard (but legitimate) axes, and thus refer to the electronic states differently. A B


76 Figure 5 2 Continued It is not clear, a priori which geometric estimate is best. Clearly, the (10,10) active s pace of the CAS calculation encapsulates the static correlation in question, but with a double zeta basis set. It also begins to describe the dynamic correlation, but it is not clear how much. Our CCSD(T) methodology is clearly superior in dynamic correl ation and our basis set is superior, but it is unclear how well it fares in terms of static correlation. There have also been prior studies using CCSD/cc pVDZ geometries, but these would be inferior to the current estimates 171 Table 5 3 Geometry of 1 A 1 C 2v structure. Lengths CCSD(T)/6 311++G(2d,2p) MC SCF(10,10)/cc pVDZ SF DFT/6 31G(d) 12 1.081 1.080 1.074 24 1.491 1.496 1.453 48 1.342 1.370 1.338 89 1.082 1.080 1.077 Angles 248 121.9 121.1 120.9 476 120.3 120.2 120.4 489 121.2 121.2 121.5 475 120.7 120.9 121.2 C


77 The CCSD(T) result comes from using 6 311++G(2d,2p) basis; the MC SCF(10,10) comes from using cc pVDZ 180 The spin flip approach does not specify functional 170 All lengths and angles are specified in terms of the atom numbering scheme in Figure 5 1. All dihedrals are zero by the symmetry of the system in question. Figure 5 3 3 B 1 C 2v / 3 A 2 D 3h structure. Table 5 4 Geometry of 1/3 B 2 C 2v / 3 A 2 D 3h structure Lengths CCSD(T)/6 311++G(2d,2p) MC SCF(10,10)/cc pVDZ SF DFT/6 31G(d) 24 1.418 1.438 1.402 89 1.081 1.081 1.076 The SF DFT /6 31G(d) 170 studies should describe the static correlation properly as well; however, of the three, SF DFT/6 31G(d) differs most widely with the other two methodologies. The differences can be quite significant; the lengths differ by 0.04 in SF DFT compared to MC S CF(10,10) and CCSD(T). We note that there is no consistent agreement between even any pair of methods in describing the carbon carbon bond lengths; significant variation exists. It would be of great interest to see a carefully done MR CC geometry optimiz ation to settle this question. As long as such variation exists, a shadow of doubt may be cast on comparing the electronic energies with multi reference wavefunctions.


78 The experimental value for the singlet/triplet energy difference is 16.1 +/ 0.1 kcal/m ol 181 point energies (ZPE) on the electronic energy of the system. Some authors choose to subtract off a calculated ZPE from the above experimental value In doing so, one assumes that the calculated geometry and harmonic frequencies are the same as those in the experiment. Given the aforementioned differences in geometry, this may be un settling. DFT/6 31G(d) e stimate of 2.03 kcal/mol, with the singlet ZPE b eing lower in energy (r The CCSD/6 311++G(2d,2p) harmonic ZPE calculated on the CCSD(T)/6 311++G(2d,2p) geometry predicts a value of 2.05 kcal/mol. Specific values are listed in the supplementary informati on. The agreement may be fortuitous or not; given the geometric differences, it is hard to be optimistic. Additionally, one might object to geometry optimizations performed using one many body method which then rely on the ZPE from a completely different many body method. For better or for worse, our current kcal/mol. TMM has five relevant spectroscopic states for its quasi degenerate states: 1 A 1 1/3 B 2 and 1/3 B 1 The 1 A 1 singlet state is a two determinant There are some conflicting reports on the nature of this state. Using CCSD(T), this is a second order saddle MC SCF calculations 180 Ho wever, Krylov and coworkers claim that it is a minimum using SF DFT Experimentally, this state is observed to be an intermediate to the stable cyclic form (reference). We would write the electronic configuration of non Abelian 3 A 2 in the


79 Abelian C 2v sub group as for the 3 B 2 state. In D 3h symmetry, 1 A 1 and 1 B 2 are degenerate. The 1/3 B 1 is a non planar C 2v molecule of electronic configuration Table 5 5 lists various estimates of the electronic energy difference using a plethor a of many body methods. Table 5 5 TMM Singlet Triplet Splittings. 3 A 2 1 A 1 S T (kcal/mol) CISD 1 2RCISD 38.3 SU CCSD 1 SU CCSD 19.9 CCSD 1 CCSD 48.1 CCSD 1 RMR 30.6 CCSD(T)f 1 CCSD(T) 28.3 CCSD(T)f 1 RMR(T) 25.8 CCSD(T)se 1 CCSD(T) 26.3 CCSD(T)se 1 RMR(T) 23.8 CR CC(2,3) 1 CR CC(2,3) 31.9 MCSCF(10,10) 2 MCSCF(10,10) 18.9 SF CIS(D) 2 20.6 CASPT(2) CASPT(2) 19.1 EOM SF CCSD 3 21.5 EOM SF CC(2,3) 3 18.2 CCSDt 4 CCSDt 22.6 CCSD(T) h 4 CCSD(T) h 19.3 CC(t;3) 4 CC(t;3) 21.7 4R BWCCSD it 5 4R BWCCSD it 17.8 4R BWCCSDT 5 4R BWCCSDT 14.9 4R BWCCSDT 5 4R BWCCSDT 15.6 CCSDT 4 CCSDT 21.7 Experiment ZPE 6 18.1 Notes : TMM Singlet Triplet Splittings. 1.) From reference 1 71 ; all calculations done using cc pVTZ; triplet geometry from CCSD(T)/cc pVTZ, singlet from 2R RMR CCSD(T)/cc pVTZ. Note, however, that this only optimized carbon carbon bond lengths; carbon hydrogen bond len gths were taken from reference 180 (unspecified basis set and many body method). 2.) From reference 170 ; all calculations done usi ng cc pVTZ on carbons, cc pVTZ on hydrogens with a geometry derived from SF DFT/6 31G(d) 3.) From reference 169 ; all calculations done using cc pVTZ geometries derived from SF DFT/6 31G(d) 4.) From reference 179 ; all core orbitals were dropped; spherical d and f functions were used in the TZ2P basis set. The selected active space was (4,4). These use geometries derived fro m SF DFT/6 31G(d) of reference 170 5.) From reference 17 5; all calculations done using cc pVTZ, with geometries from CAS SCF(4,4)/cc pVDZ 6.) The ex perimental value is given in reference 181 ; the zero point energy (ZPE) comes from refer ence 170 (SF DFT/6 31G(d)) harmonic vibrational frequencies


80 Given that basis set extrapolations have not been done except in our DEA CCSD calculation, all reported ener gies might vary within ~2 kcal/mol (except DEA) due to valence basis completeness. This says nothing of the error intrinsic of imperfect geometries used. The CISD estimate unsurprisingly fails. We note that CCSD fails miserably far worse than in the ca se of the methylene. The static correlation is clearly insurmountable. It is of note that all combinations of CCSD, CCSD(T), RMR, RMR(T), and CR CC(2,3) are still very far from experimental answer. The MC SCF(10,10) and CASPT2(10,10) calculations are im pressively accurate, as is SU CCSD. The various spin flip methods are much closer to the mark; EOM SF CC(2,3) is the only method that finds the correct answer. It is surprising that an active space method such as CCSDt comes only as close as it does at 2 3 kcal/mol. The Brillouin Wigner estimates curiously get worse with the addition of triples, disconcertingly. The CCSDT result is close, relative to the smaller basis used. body methods are listed in Table 6 at our geometries, a CCSD. Table 5 6 TMM Singlet Triplet Splittings at CCSD(T)/6 311++G(2d,2p) Geometries cc pVDZ cc pVTZ cc pVQZ 3 4 Extrap RHF / UHF 102 99 98 97 MBPT(2) 40 37 35 34 CCSD 46 48 48 48 CCSD(T) 22.7 25.5 25.6 25.6 28.5 30.6 31.0 31.3 DEA CCSD 24.1 24.4 24.2 24.1 DIP CCSD 21.5 21.1 20.8 20.6 CCSDT 23.1 Experiment ZPE 18.1 TMM Singlet/Triplet Splittings, kcal/mol. Numbers are listed to one decimal purely for computational comparison sake, and are not intended to imply accuracy to so many


81 decimals. Performed by the authors using CCSD(T)/6 311++G(2d,2p) for geometries and frequencies. Core orbitals were dropped for energy calculations The DEA CCSD calculation provides a number consistent with the other multi reference coupled cluster methods. The curiously very far from the other coupled cluster methodologies. There is no real difference in CCSD a nd CCSD(T) methods at our geometries vs. previous geometries. It is unsurprising that SCFs and MBPT(2) fail miserably. There is a slight difference in our CCSDT/cc pVTZ calculation and the aforementioned CCSDT/aug cc pVTZ calculation at CCSD(T)/cc pVDZ g eometry by about a kcal/mol. The DIP results mirror the accuracy of the DEA calculation. Benzynes We consider the more complicated case of para benzyne first. Para benzyne is indeed infamous for its comp utational challenge 171, 179, 183 185 A complete d escription of its difficulties has been given by Crawford et al. 183 Experimental research has come to contra dictory conclusions 186, 187 The molecule is shown in Figure 5 3 along with its associated geometric parameters in Table 5 7. The multi referen ce nature of the wavefunction comes from the singly occupied orbitals on carbon atoms 1 and 4. This degeneracy is the root of the complexity. Figure 5 4 shows the relevant orbital combinations that lead to degeneracy.


82 Figure 5 4 Para Benzyne and our associated axis choice. Table 5 7 Geometric Parameters of Para Benzyne, Singlet followed by Triplet Lengths TCSCF MkCCSD/cc pVTZ CCSD(T)/cc pVDZ CCSD(T)/6 311++G(2d,2p) 12 1.363 1.390 1.377 23 1.430 1.427 1.418 CH 1.080 1.097 1.082 Angles 612 124.8 125.1 125.2 123 117.6 117.4 117.4 H23 118.4 120.2 119.7 Lengths TCSCF MkCCSD/cc pVTZ CCSD(T)/cc pVDZ CCSD(T)/6 311++G(2d,2p) 12 1.376 1.395 1.386 23 1.405 1.417 1.412 CH 1.081 1.098 1.083 Angles 612 127.0 126.3 126.4 123 116.5 116.8 116.8 H23 121.1 121.1 121.0 All lengths and angles are specified in terms of the atom numbering scheme in Figure 5 3. The letter H is specified for the hydrogens, given that they are all equivalent in this structure. The CCSD(T)/cc pVDZ 171 geometry differs significantly from our own. The MR CC geometry derives from private communication 188 In Figure 5 particle wavefunctions have a g symmetry, but in reality, the total wavefunction of the system will also have one determinant with two electrons of a g symmetry, and one determinant two electrons with


83 b 3u symmetry. This is because they are near degenerate. Linear combinations of these two orbitals is symmetry allowed with (a g ) and (b 3u ) orbitals corresponding to carbon carbon interactions. The result of these linear combinations is to de stabilize the carbon carbon interactions and greatly stabilize the b 3u C H bonds. The multi reference quality thus affects the global geometry of the molecule. It has been shown that a restricted Hartree Fock description tends to localize each electron, such that they do not interact 183 An unrestricted reference wavefunc tion tends to de localize the electronic density, better allowing for the interaction between the electrons to be described. Figure 5 5 Orbital diagram of Para Benzyne. A saving grace of para benzyne is that it does not have multiple low lying multiplet varied symmetry states to complicate analysis, as is the case in methylene and especially TMM. The singlet state is 1 A g in D 2h symmetry; its electron configuration


84 follows the above complexity. The triplet is 3 B 3u in D 2h symmetry, with configuration This designation follows our choice in axis, using the specification in Figure 5 3. We note this in particular because other good researchers made alternative axis choices, and call this state 3 B 1u via a non standard choice in primary ax is for point group symmetry. There are two reported values for the singlet/triplet g ap of para benzyne. Lineberger et al 187 estimates the gap to be 3 .8 +/ 0.3 kcal/mol; Leopold et al 186 estimates the gap to be 2.1 +/ 0.4 kcal/mol. The highest level calculations to date on the ZPE difference would be our own. We have performed optimization and harmonic vibrational frequency using CCSD(T)/6 311++G(2d,2p), producing a difference of 0.4 k cal/mol between the two states, with the triplet being higher in energy. Specific values are listed in the supplementary information. The TC SCF Mk CCSD/cc pVTZ geometry optimization of Allen et al 188 is likely to be the most trustworthy geometry yet av ailable on this heavily multi reference system; the broken symmetry CCSD(T) geometry produces a geometry in relatively close agreement with it, thus giving confidence to our estimate of the ZPE. Table 5 5 illustrates this concordance. We briefly discuss m eta benzyne and ortha benzyne; their multi reference quality is not so drastic. Figure 5 5 shows our numbering system for meta benzyne, and Table 5 8 gives the geometries calculated in our work using CCSD(T)/6 311++G(2d,2p) as well as TC SCF Mk CCSD/cc pVTZ. The agreement between the two methodologies is extremely high, excepting two singlet angles around the diradical centers. This gives considerable confidence to the geometries used. Meta benzyne singlet is a C 2v structure whose wavefunction follows 1 A 1 symmetry; its triplet is also C 2v


85 corresponding to 3 B 1 The electron configuration of the triplet is The ortha structure is not list ed below, but available in the Supplementary I nformation; the rationale is forthcoming. The ZPE dif ference corresponding to a MBPT(2)/6 311++G(2d,2p) harmonic vibrational frequency on the CCSD(T)/6 311++G(2d,2p) optimized meta benzyne structure yields a difference of 0.60 kcal/mol with the singlet ZPE being lower. Specific values are listed in the supp lementary information. Figure 5 6 Meta Benzyne. The same axis is used here as in all previous examples. Note that the use of the However one could choose to draw a Lewis structure alternatively describing the two lone particles forming a triple bond alone one line segment of the hexagon. The anti tumor properties of these compounds are reflected more properly in the diradical nature rather than an alkyne representation. Table 5 9 lists the values for the electronic energies calculated using a variety of methodologies. There is, again, a diversity of different geometry sources; we believe our geometries aid future work considerably, due to their higher qua lity. It is hard to know if errors from the commonly used CCSD(T)/cc pVDZ geometry influence the quality of these high level single point energy calculations. We can see the drastic failures of traditional CCSD in both cases. RMR acting on the para benz yne singlet is


86 also clearly insufficient, as it predicts the wrong ordering of states; CR CC(2,3) also fails badly for the same reason. RMR combined with SU CCSD underestimates the energy gap significantly. All the remaining computations of para benzyne including perturbative triples tend to favor the Lineberger experiment, being closer to 3.4 than 1.7. However, none of these are basis set extrapolated, and values might reasonably shift closer to the 1.7 estimate. Table 5 8 Geometries of Meta Benzyne Sin glet Lengths TC SCF MkCCSD/cc pVTZ CCSD(T)/6 311++G(2d,2p) 12 1.363 1.378 2H 1.076 1.077 34 1.373 1.381 4H 1.080 1.082 45 1.397 1.406 5H 1.084 1.086 Angles 123 96.0 99.4 234 138.5 135.7 345 116.6 117.4 456 113.7 114.4 H61 120.8 120.5 Triplet Lengths TCSCF MkCCSD/cc pVTZ CCSD(T)/6 311++G(2d,2p) 12 1.379 1.385 2H 1.082 1.083 34 1.377 1.383 4H 1.080 1.081 45 1.397 1.404 5H 1.082 1.084 Angles 123 115.0 115.2 234 124.9 124.7 345 116.8 116.9 456 121.4 121.5 H61 122.4 122.3 All units are in and degrees. All lengths are described by numbered line segments corresponding to Figure 5, as are the angles. The letter H indicates the hydrogen atom as a line segment point / angle vertex. All dihedrals are zero by sym metry. Reference 188 gives the MkCCSD geometry.


87 The DEA CCSD of Table 5 10 below for para benzyne is conspicuously on the mark for agreement w ith the 3.4 kcal/mol Given that our value has basis set extrapolation on a more refined geometry, our calculati on strongly supports, albeit not decisively, the L eopold experimental value of 3.4 kcal/mol. The CCSD calculations on our refined geometries are identical to the previous geometry. This might suggest that the old CCSD(T)/cc pVDZ is not so different; howe ver, the CCSD(T) single point energy calculation does march toward the 3.4 kcal/mol estimate on the improved geometry. CCSD features odd behavior at the cc pVQZ level; it either converge s to an excited state SCF or the larger basis set proves problematic for the 2 anion (because a bound state may not exist). Further study is needed. Regardless, the behavior is clearly not as good as DEA CCSD. For meta benzyne, given that methods in Tab le 5 9 have not performed a basis set extrapolation, one may reasonably believe that all of the proposed methods may reach 2 kcal/mol higher toward the experimental answer when perturbative triples are present. Without including perturbative triples, howe ver, the estimates are significantly off the mark. The DEA CCSD value is basis set extrapolated in the valence space and is comparable to other coupled cluster methods. It would seem triples are needed to account for the remaining energy difference betwe en the DEA result of 18.5 and experimental 20.4 kcal/mol. The same logic applies for DIP CCSD. This speaks to the ease of use of DEA CCSD: it achieves the same values in this case at a cheaper computational cost without using an active space.


88 Tab le 5 9 Meta and Para Benzyne Splittings Singlet Triplet E T S Meta (kcal/mol) E T S Para (kcal/mol) CCSD 1 CCSD 7.2 19.4 RMR 1 CCSD 13.8 2.1 RMR 1 SU CCSD 16.2 0.7 CCSD(T) 1 CCSD(T) f 17.1 3.1 CCSD(T) 1 CCSD(T) se 18.0 4.0 CR CC(2,3) 1 CR CC(2,3) 16.4 1.8 RMR(T) 1 CCSD(T) f 17.7 3.2 RMR(T) 1 CCSD(T) se 18.7 4.1 Delta ZPE 2 0.7 0.3 Delta ZPE 3 0.6 0.4 Expt. 1 4 21.0 +/ 0.3 3.8 +/ 0.3 Expt. 2 5 2.1 +/ 0.4 Expt. 1 ZPE 3 20.4 +/ 0.3 3.4 +/ 0.3 Expt. 2 ZPE 3 1.7 +/ 0.4 1.) From reference 1 71 ; all calculations use a modified cc pVTZ in which the f functions on carbons and the d functions on hydrogens are omitted. All geometries com e from reference 189 which used unrestricted CCSD(T)/cc pVDZ geometries and harmonic vibrational frequencies for para benzyne, and restricted CCSD(T)/cc geometries and harmonic vibrational frequencies. 2.) From reference 189 ; methodol ogy is the same as in Note 1 above 3.) From this work. Geometries come from CCSD(T)/6 311++G(2d,2p), harmonic freq uencies CCSD(T)/6 311++G(2d,2p) for para benzyne; geometries come from CCSD(T)/6 311++G(2d,2p), harmonic frequencies from MBPT(2)/6 311++G(2d,2p) for meta benzyne. 4.) From reference 186 5.) From reference 187 Our results on our improved geometries are given in Ta ble 5 10. Table 5 10 Meta and Para Benzyne Singlet/Triplet Splittings, kcal/mol. Meta cc pVDZ cc pVTZ cc pVQZ 3 4 Extrap RHF / UHF 50 49 49 49 MBPT(2) 55 60 62 63 CCSD 9.9 10.5 10.7 10.8 CCSD(T) 20.4 22.1 22.4 22.7 20.9 23.1 20.9 19.3 DEA CCSD 17.1 18.3 18.4 18.5 DIP CCSD 18.0 19.7 18.3 17.3 Experiment 20.4 +/ 0.3 Para cc pVDZ cc pVTZ cc pVQZ 3 4 Extrap RHF / UHF 90 90 89 89 MBPT(2) 25 27 28 29 CCSD 17.8 19.2 19.5 19.7 CCSD(T) 4.0 3.6 3.5 3.4 0.4 0.8 0.9 0.9 DEA CCSD 3.0 3.4 3.4 3.4 DIP CCSD 3.9 4.4 14.8 NA Experiment 3.4 +/ 0.3 1.7 +/ 0.4


89 Numbers are listed to one decimal purely for computational comparison sake, and are not intended to imply accuracy to so many decimals. Performed by the authors using geometries from CCSD(T)/6 311++G(2d,2p), harmonic frequencies CCSD(T)/6 311++G(2d,2p) for para benzyne; geometries come from CCSD(T)/6 311++G(2d,2p), harmonic frequencies from MBPT(2)/6 311++G(2d,2p) for meta benzyne Core orbitals w ere dropped for energy calculations We do not list the values of DEA CCSD and other many body methodologies for ortha benzyne. This is because DEA CCSD fails strongly in this case, illustrating a weakness of the method. The +2 charge ortho benzyne SCF is actually more multi reference than neutral ortha benzyne. The predicted singlet/triplet splitting is ~60 kcal/mol, nowhere near the correct value of ~37 kcal/mol. The same disaster occurs for the 2 charged ortho benzyne in the DIP CCSD calculation (pre dicted gap of ~80 kcal/mol). To examine this further, we examined the shifting of orbital energies in response to changing charg e in the benzynes. In Table 5 11 we list the cc pVDZ orbital energies (eV). A small HOMO LUMO gap, unto itself, is suggestiv e of multi reference character. Note that the biggest shifts in the HOMO LUMO gap belong to both charged species of ortho benzynes (about 6eV), and these are two systems for which the DIP/DEA methods suffer the most. We emphasize that this is a potential ly surmountable obstacle; in so much as we can converge to the ground electronic SCF wavefunction, we may still be able to use such methods. The issue of the proper metric to gauge the multi reference quality of a wavefunction (other than CI coefficents) remains difficult.


90 Table 5 11 HOMO/LUMO Gaps of Benzynes Para Para +2 Para 2 HOMO 7.3 HOMO 21.6 HOMO 9.1 LUMO 1.2 LUMO 13.3 LUMO 11.2 Gap 6.1 Gap 8.3 Gap 2 Meta Meta +2 Meta 2 HOMO 8.3 HOMO 22.1 HOMO 6.5 LUMO 1.4 LUMO 13.8 LUMO 14 Gap 9.7 Gap 8.2 Gap 7.6 Ortho Ortha +2 Ortha 2 HOMO 9.5 HOMO 22.1 HOMO 7.1 LUMO 1.8 LUMO 15.3 LUMO 13.8 Gap 11.3 Gap 6.7 Gap 6.7 Orbital energies, eV, for the singlet form of each molecule, cc pVDZ. Delta cation gap refers to the difference between the HOMO LUMO gap of the neutral and cation. The greater change in the delta gap is suggestive of a higher multi reference quality. Further Thoughts in Conclusion The DEA framework blossoms in utility in making complicated multi reference wavefunctions into straightforward single reference calculations. In choosing a single reference charged SCF, and adding/subtracting the appropriate particles to return to the multi reference system, we circumvent the challenges of static correlation. Choosing such pathological molecules as methylene, TMM, and benzynes, and still seeing such great success, illustrates the power of the methodology.


91 CHAPTER 6 MECHANISM OF RDX DECOMPOSITION The Stage Has Been Set, Let Us Dance A plethora of mechanisms have been previously proposed in litera ture 10 21 As stated previously, the conflict in interpretation of experiments combined with the suggestive but not deterministic nature of calculations begs resolution. Computational efforts to describe the condensed phase in an ab initio manner are currently im possible. However, we must walk before we can run. I n describing the gas phase with high accuracy, we can validate which experimental methods are correct (and rely on those styles of methods for the condensed phase). Additionally, the first order descri ption of a condensed phase system is the gas phase. Moreover, the system is generally not considered to have ionic mechanisms, decreasing the importance of environmental effects. Our efforts using coupled cluster techniques allow for a decisive answer to the gas phase mechanism in time. The largest effects, and most likely the true answer, is available from the employed methodology. Methodology Our protocol has general trends, but is ultimately adaptive to the various specific molecules in question. In general, we used MBPT(2)/cc pVTZ 102, 104 geometries for all energetic minima and their associated harmonic vibrational frequency calculations to confirm that they are minima. We used the ACES III 107 software for its efficiency in massively parallel ab init io calculations. We found transition states by first performing a PES scan with frozen coordinates except along the expected reaction path. This was done using B3LYP/6 31G(d,p) 148 using GAMESS. We found it necessary to calculate force constants at every step of the optimization in order to stay on trajectory to an


92 energetic maximum rather than a minimum; GAMESS offers automated calculation of force constants at every step. We chose to use KS DFT rather than MBPT(2) given the difficulty in transition sta s code for KS DFT as opposed to its inefficiency in wave function methodologies. If a suspected intermediate/transition state had small harmonic vibrational frequencies (below 30 cm 1 in magnitude), a CCSD /cc pVTZ Hessian was calculated to confirm the status as an energetic minimum/maximum. Without exception, single point energies have been perfor med using dropped core CCSD(T)/ cc pVTZ energies, choosing an unrestricted Hartree Fock reference where necessary. We chose an unrestricted reference because it offers spin polarization as well as being a variational method with no extra imposed constraints. We again used the ACES II I software. I n cases where both the reactant and tra nsition state were singlets, DEA CCSD/cc pVDZ calculations were employed as a test of the quality of CCSD(T)/cc pVDZ calculations on the multi reference transition states. This was only done on singlets because a code for doublets has not been written. Proposed Mechanisms Approximately 2 /5 of the mechanism given in this document are n ewly devised; the remainder has been described in some form in the literature. In almost all cases, this work represents the first time the mechanisms have been represented in the formal arrow pushing diagra ms, listing all the details. There are three major starting points often considered: N N homolysis, HONO 1, 6 2, and 6 3, respectively. There is experimental evidence for each, respect i vely We highlight the strongest evidence of each case, briefly


93 Figu re 6 1 N N Homolysis Mechanism Figure 6 2 HONO Elimination Mechanism Figure 6 3 Triple Whammy Mechanism The N N homolysis mechanism has been supported by the earliest studies of RDX. A gas phase kinetic studies demonstrated an alleged rate dependence on nitrogen dioxide 10 12 thus supporting the view. However, this is in contradiction to experiments which demonstrate a primary hydrogen isotope effect o n reaction


94 kinetics 14 ; this effect inherently supports the HONO reaction mechanism of Figure 6 2. In the HONO elimination, an oxygen of the nitramine group abstracts a hydrogen atom, causing the formation of a cyclic nitroalkene. Some mass spectrometry e xperiments 81 support the existence cyclic nitroalkene intermediates. However, this would not explain the nitrogen dioxide dependence previously observed 10 12 Lastly, one lone experimental group 15 found support for the triple whammy mechanism in Figure 6 3 No other experiments or calculations attemptin g to validate this mechanism have ever been found 10 21 The triple whammy mechanism corresponds to simultaneous bond break ing of 3 C N covalent bonds hetero lytically to form hydrogen cyanide. At a firs t glance, this mechanism seems silly; it would presumably be much harder to break 3 bonds at once, as opposed to the other mechanisms. The subsequent reactions for further degradation are important. Obviously, if a low barrier exists in N N homolysis, but subsequent reactions require higher barriers, N N homolysis is irrelevant. We emphasize these three first step possibilities simply as an organization tool to classify the many possible twists and turns of the mechanism. We emphasize that we explicitly d iscuss the reactions whose results were promising. We investigated many other reactions, available in Supporting Information, which proved to be completely un tenable (the intermediates were highly unstable). N N Homolysis The reactants and products of Figure 6 1 have been studied previously, as discussed in Chapter 2. We have refined our previous work by using aforementioned MBPT(2)/cc pVTZ geometries and frequencies for RDX (a basis set improvement), as well as optimizing RDR structures. We did not analyze the RDR structures for harmonic vib rational frequencies, earlier. W e took it for granted that the famous observed RDR


95 structures, such a simple devia tion from the RDX structure, were correct. Although the refined geom etries are not significantly different from the results of previous efforts, the harmonic vibrational frequencies indicate that the axial axial (AA) form is unstable; it collapses into the AE (axial equatorial) form readily. Table 6 1 gives the energies o f the RDX conformers as well as RDR. Table 6 1 RDX and RDR Energies kcal/mol AAA 1.3 AA 49.0 EEE 5.7 EE 50.6 AAE 0.0 AE 49.1 AEE 1.0 Boat 49.8 Twist 0.0 Twist 47.8 Boat 0.0 All structures are MBPT(2)/cc pVTZ optimized; energies calculated using CCSD(T)/cc pVTZ Numbers are listed to one decimal without implying accuracy to within tenths of kcal/mol; they are written to these decimals purely for QIC. We note t hat upon refining the structures of the conformers the energy levels are shifted even closer together. It is interesting to note that, over a decade of refining these energetic s of RDX, the conformers have consistently come closer together in energy 3, 72 The transition state for the reaction in Figure 6 1 does not involve any direct bond formation, in the chemistry sense, so much as the RDX geometry morphs to RDR until the transition state is met. The RDR conformers are again effectively degenerate. Regardless of the conformer pair studied between RDX and RDR, t he barrier to this proces s is ~ 59 kcal/mol (with one pairing being about ten kcal/mol higher, and this path obviously disfavored) From RDR, two meaningful possibilities exist: a CN bond fission to break up the ring, or yet another N N homolysis. We consider first subsequent N N homolysis first. Subsequent N N homolysis is demonstrated in Figure 6 4. Subsequent N N homolysis is not supported by isotopic scrambling experiments 17 which show only one


96 nitro group breaks apart homolytically. However, this does generate cyclic nitroa lkenes, which is supported by mass spectrometry 81 Figure 6 4 Secondary N N Homolysis Concurrent with the proposed N N homolysis is the radical abstraction of a proton in order to form the cyclic nitroalkene. The transition state involves a hydrogen migration between a carbon atom and a nitrogen atom, at a barrier of 51 kcal/mol. It is therefore a lesser barrier than the primary N N homolysis, and thus i s accessible kinetically. This total process releases 33 kcal/mol, as the product is very stable compared to RDR (as well as producing a gaseous nitrogen dioxide). From here, ring scission must begin. A form of triple whammy is unlikely; the ring is not symmetric such that it is implausible to think three different CN bonds would equally break at th e same time. The weakest CN bond cleavage is shown in Figure 6 5. Figure 6 5 CN Ring Scission from the Secondary N N Homolysis


97 This path has a barrier of 49 kcal/mol, and the products are 42 kcal/mol higher in energy than the reactant. This is not surprising given the diradical nature of the products. The barrier is again lower than the initial N N homolysis, and this sub reaction is possible if primary N N homolysis occurs. From here, any number of possible formation of gaseous products can occur and subsequent barriers are not a major concern. Once the ring is broken and products formed from which one can easily imagine making common gas products, the task is done. We now return to consider the possibility of CN bond breaking following primary N N homolysis (not the secondary N N homolysis route). This is shown in Figure 6 6. Figure 6 6 Primary CN Bond Cleavage of RDR This reaction has a surprisingly small barrier of only 29 kcal/mol. We normally associate formation of cyclohexanes as a particularly energetically favorable process, but in this case, the heterocyclic nature of the ring with heavy electron withdrawing groups is such that the CN bond cleavage is not kinetically difficult. The product is 29 kcal/mol uphill of the RDR startin g point The product is still a long chained radical, whose dissociation is not immediately obvious. We discuss two possible paths from here: isomerization vs immediate CN cleavage again.


98 Immediate CN bond fission (this is now the second CN bond fission !) is depicted in Figure 6 7. Trial and error demonstrated that this is the weakest CN bond. Figure 6 7 Second CN Cleavage of RDR The CN bond here is again weak; the transition state barrier is 13.5 kcal/mol. The product is only 16.5 kcal/mol above the reactants, again mak ing this kinetically accessible relative to the original N N homolysis. The nitrogen radical containing product is not a common species whose degradation is obvious; a plausible subsequent carbon cleavage is shown in Figure 6 8. Figure 6 8 CN Cleavage of Twice CN Cleaved RDR This particular path represents the cleavage of the weakest CN bond, and results in the formation of products which easily can form common gaseous products of the oxidation of RDX. The barrier is 15 kcal/mol and the products lie 11.6 kcal/mol uphill of the reactants.


99 We now return to the possibility of isomerization of the product of Figure 6 6. This is shown in Figure 6 9 Figure 6 9 Isomerization of CN Cleaved RDR This reaction is, in effect, a hydrogen migration. The net effect is to trade a nitrogen radical for a terminal carbon radical. The transition state barrier is a mere 11 kcal/mol above the reactant, and the reaction is 8 kcal/mol downhill from the reactant. Thus, this isomerization is the do minant p ath followed by N N homolysis, given the small barrier and how energetically favorable it is. Subsequent CN cleavage to make hydrogen cyanide is depicted in Figure 6 10. Figure 6 10 CN Cleavage of the Most Stable N N Homolysis Path The reaction is 19 kcal/mol uphill; this is unsurprising given that we are trading one terminal carbon radical for a smaller chained terminal carbon radical (less donation possibility to account for electron deficiency). The transition state barrier is a mere 8 kcal/m ol, making this reaction fast.

PAGE 100

100 If N N homolysis is the correct initial step, it is thus followed by CN bond fission rather than a secondary N N homolysis (confirming the isotope scrambling experiment 17 ); this product then undergoes isomerization shown in F igure 6 9. All other paths are significantly more expensive energetically. HONO Elimination The HONO elimination pathway of Figure 6 2 has a barrier of 47 kcal/mol to form its cyclopentane transition state. This barrier is thus smaller than that of N N h omolysis; if subsequent reaction steps feature barriers no larger, the putative N N homolysis mechanism of great literature popularity is untrue. It should be noted that because it is easier to do experiments/calculations on such a simple path, it has received undue popularity. Moreover, finding the fused cyclic ring transition state was particularly difficult. The product of Figure 6 2 is more stable than RDX by 6 kcal/mol. It also is a nitro cycloalkene, consistent with mass spectrometry experiment s. Performing yet another HONO elimination is shown in Figure 6 11. Figure 6 11 Second HONO Elimination HONO and formation of an increasingly aromatic cyclic nitroalkene is downhill by 8 kcal/mol. The transition state for the second HONO elimination is nearly identical to th e

PAGE 101

101 the transition state is 44 kcal/mol, close to the first elimination activation energy The third HONO elimination is shown in Figure 6 12, feat uring the formation of the very stable aromatic product TAZ. Figure 6 12 Third HONO Elimination The activation energy is again similar: 41 kcal/mol. The product (TAZ) is 27 kcal/mol more stable than the starting reactants, again favoring this path (mor e and more stable products, smaller and smaller activation energies). However, herein lies the challenge of this degradation path: how do you breakup such a stable aromatic structure? whammy rated in Figure 6 13.

PAGE 102

102 Figure 6 13 Simultaneous Heterolytic Bond Cleavage This degradation path is impossible, as the activation barrier is 88 kcal/mol. Breaking up the aromaticity all together in one step is not possible. Desperate efforts to perhaps break up the ring before total aromaticity forms by considering the product of Figure 6 11 is shown in Figure 6 14. Figure 6 14 Simultaneous Heterolytic Bond Cleavage on Pre Aromatic Product This path is easier, with a barrier of 66 kcal/mol, but still very much higher in energy than the original HONO step. We cannot consider starting from product of Figure 6 11, since it is implausible to simultaneously break 3 unequal bonds. Herein lies the challenge to this mechanism: how does one ever break up the very stable products one makes along the way? The HONO path seems to trade easier initial steps for very hard steps as one makes increasingly aromatic rings. Further work may try to

PAGE 103

103 imagine steps more kinetically accessible after the first HONO elimination; however, none are present in the literature nor imagined by the author. Triple Whammy The triple whammy mechanism of Figure 6 3 has not been supported by any experiment or calculation aside i ts original proposal 15 In all candor, the author h imself repeatedly scoffed at the possibility that all three bonds would break simultaneously, given the literature history of this mechanism. Its transition state is exceptionally difficult to isolate; the energetic peak has a very small rate of change. However, the activation energy of this process is a mere 30 kcal/mol! This is significantly easier than either HONO elimination or N N homolys is. On the basis of electronic energy, it is clear that the triple whammy mechanism dominates RDX decomposition. In particular, we note the very different estimate given previously using B3LYP/6 31G(d). The original work estimating this transition state barrier predicted a barrier of 59 kcal/mol (despite the paper claiming that B3LYP/6 31G(d) 37 is accurate to with in a kcal/mol or two). We are open minded that a CCSD(T) estimate of such a degenerate transition state may have some error, but this sort of difference is irreconcilable. Validation of Methodology on Transition States It has been discussed that transitio n states are inherently multi reference, and it should not be taken for granted that CCSD(T) will reliably describe the wavefunction. The DEA CCSD method is reliable to describe the static correlation of transition states which describe one bond breaking, one bond forming. Although this does not describe all of the aforementioned bond breaking processes, it would at the very least have more static correlation represented than a SR CCSD(T) calculation. It is not true that DEA CCSD is categorically better than CCSD(T); the presence of perturbative triples is an

PAGE 104

104 advantage of CCSD(T). Rather, general agreement between the methods boosts confidence that the numbers are not random. To examine this, we looked at the cc pVDZ energies of both methods on transiti on states where we could perform DEA CCSD in Table 6 2. Table 6 2 DEA CCSD and CCSD(T) cc pVDZ Energies on Transition States Reaction DEA CCSD CCSD(T) Fig 6 1 3 8 3 3 Fig 6 11 4 4 47 Fig 6 12 39 44 Fig 6 13 4 5 40 Fig 6 14 89 88 All units are kcal/mol; transition states can be found in the associated figures. We observe a general agreement on the energies. We do not see very drastic differences in the energies; differences on the order of a few kcal/mol are understandable given the differences in pertur bative triples. This gives confidence that we have some basic description of the wavefunction of the transition states. This is not meant as a rigorous test; there is, as of yet, no real way to measure the multi reference effectiveness of a coupled clust er method. Further Work We emphasize that the existing work represents the most complete study done to date on the mechanism of nitramines. The largest contributions toward the correct answer have been accounted for, and it is likely that the internal ene rgy described at the CCSD(T)/cc pVTZ level alone yields the correct answer. Nonetheless, there remain further corrections for which we hope to account in future work. The two largest effects that we plan to account for next are valence basis set extrapola tion and calculation of CCSD partition functions using cc pVTZ. The accuracy

PAGE 105

105 of the electronic energy can be off by a few kcal/mol if one does not extrapolate. This is realistically attainable for all the studied systems, computationally. As importantly it is desirable to calculate the Gibbs energies, rather than just internal energies. It is also realistically attainable to get the CCSD/cc pVTZ partition functions, computationally, albeit expensive. We have estimates of the Gibbs energies of the reac tion energies (not transition states) using MBPT(2)/cc pVTZ partition functions ; we did not consider them in this analysis only because they would not be consistent with our transition states. Due to the sensitivity of transition state optimizations, we h ave not had sufficient time to do CCSD optimizations; consequently, we have KS DFT partition functions. These may be adequate, but it does not make sense to compare them to MBPT(2) partition functions. The entropic effect is especially interesting in the triple whammy mechanism. Accomplishing these two goals would most certainly remove any ambiguity in the results. The next leading order effect would be the accuracy of the low frequency modes. This is a current area of intense research; leading order ef fects may be captured using 3 rd order polynomials for anharmonic effects, although many options are possible. The practicality of performing these calculations is unclear. The effects, numerically, would only be important if the result of the aforementio ned basis extrapolation and Gibbs energy calculations resulted in bringing transition state barriers to within a few kcal/mol of one another. The necessity of this will be clear in time. The most interesting future work would be to apply the results of th e RDX mechanistic study to HMX and CL 20. The HONO and N N homolysis mechanisms are likely to be near identical for HMX; we may likely borrow our RDX structures as starting

PAGE 106

106 points for the HMX efforts. The CL 20 mechanisms for HONO and N N homolysis are likely to be significantly different due to the caged structure. The triple whammy entropically. One might expect it to dominate HMX, as well; however, it obviously cannot occur for CL 20. Concluding Remarks The various paths of RDX decomposition have been explored from the three primary proposed initial steps: N N homolysis, HONO elimination, and triple whammy. We have found that the N N homolysis path proceeds via CN cleavage into an isomerization into another CN cleavage. The largest barrier in this path is the initial N N homolysis. The HONO path is more kinetically accessible than N N homolysis in the first few steps, but suffers from having successively more kinetically diffic ult steps due to the formation of aromatic products. The triple whammy pathway, the least supported pathway in the literature, as most dominant in the gas phase, by 19 kcal/mol compared to other paths. With a total activation barrier of 30 kcal/mol, it i s the most kinetically accessible, and the dominant mechanism in terms of electronic energies. Modulating the shock sensitivity would be best done by stabilizing or de stabilizing the transition states. It is not immediately obvious whether electron with drawing or donating groups would do the trick; trial and error is required.

PAGE 107

107 LIST OF REFERENCES (1) Choi, C. S.; Prince, E. Acta Crystallogr. 1972, B28, 2857 2862. (2) Karpowicz, R. J.; Sergio, S. T.; Brill, T. B. Ind. Eng. Chem. Prod. Res. Dev. 1983, 22, 363 365. (3) Rice, B. M.; Chabalowski, C. F. J. Phys. Chem. A 1997, 101, 8720 8726. (4) Harris, N. J.; Lammertsma, K. J. Am. Chem. Soc. 1997, 119, 6583 6589. (5) Byrd, E. F. C.; Scuseria, G. E.; Chabalowski, C. F. J. Phys. Chem. B 2004, 108, 13100 13106. (6) Forel, M. T.; Filhol, A.; Clement, C.; Paviot, J.; Rey Lafon, M.; Richoux, G.; Trinquecoste, C.; Cherville, J. J. Phys. Chem. 1971, 75, 2056 2060. (7) McCrone, W. Anal. Chem. 1950, 22, 954 955. (8) Shimojo, F.; Wu, Z.; Nakano, A.; Kalia, R. K.; Vashishta, P. J. Chem. Phys. 2010, 132, 094106. (9) Podeszwa, R.; Rice, B. M.; Szalewicz, K. Phys. Rev. Lett. 2008, 101, 115503. (10) Robertson, A. J. B. Trans. Faraday Soc. 1949, 45, 85 93. (11) Cosgrove, J. D.; Owen, A. J. Chem. Commun. 1968, 6, 286. (12) Rauch, F. C.; Fanelli, A. J. J. J. Phys. Chem. 1969, 73, 1604 1608. (13) Owens, F. J.; Sharma, J. J. Appl. Phys. 1980, 51, 1494 1497. (14) Bulusu, S.; Weinstein, D. I.; Autera, J. R.; Velicky, R. W. J. Phys. Chem. 1986, 90, 4121 4126. (15) Zhao, X.; Hintsa, E. J.; Lee, Y. T. J. Chem. Phys. 1988, 88, 801 810. (16) Behrens, R.; Bulusu, S. J. Phys. Chem. 1992, 96, 8877 8891. (17) Botcher, T. R.; Wight, C. A. J. Phys. Chem. 1994, 98, 5441 5444. (1 8 ) Robertson, A. J. B. The Th ermal Decomposition of Explosives. Trans. Faraday Soc. 1948, 44 (19 ) Smith, L. C.; Rogers, R. N. Estimation of Preexponential Factor f rom Thermal Decomposition Curve of an Unweighed Sample. Anal. (20 ) Daub, G. W.; Rogers R. N. Scanning Calorimetric Determination of Vapor Phase

PAGE 108

108 (21 ) Walker, F. E.; Shaw, R. Estimated Kinetics and Thermochemistry of Some Initial Unimolecular Reactions in the Thermal Decomposition of 1,3,5,7 Tetr anitro 1,3,5,7 tetraazacyclooctane in the Gas (22 ) Okamoto, Y.; Croce, M. Cationic Micellar Catalysis of the Aqueous Alkaline Hydrolyses of 1,3,5 Triaza 1,3,5 trinitrocyclohexane and 1,3,5,7 Tetraaza 1,3,5,7 tetra nitrocyclooctane. J. Org. Chem. 1979, (23 Octahydro 1,3,5,7 tetranitro 1,3,5,7 tetrazocine and Their Temperature Dependence. J. Phys. Chem. (24 ) Brill, T. B.; Landers, A. G. Pressure Temperature Dependence of Polymorph Interconversion in Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. (25 Thermal Analysis and Relationship to Propellants. J. Phys. Chem. 1982, (26 ) Ebinger, M. H.; Janney, J. L.; et al. Deuterium Isotope Effects in Condensed Phase Thermochemical Decomposition Reactions of Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. Chem. 1985, (27 ) Velicky, R. W.; Autera, J. R.; Weinstein, D. I.; Bulusu, S. Deuterium Kinetic Isotope Effect in the Thermal Decomposition of 1,3,5 Trinitro 1,3,5 triazacyclohexane a nd 1,3,5 7 tetranitro 1,3,5,7 tetracyclooctane: Its Use as an Experimental Probe for Their Shock (28 ) Haas, Y.; Greenblatt, G. D.; Zuckermann, H. OH Formation in the Infrared Multiphoton Decompositio n of Jet Cooled Cyclic Nitramines. J. Phys. Chem. 1987, 91, (29 ) Miller, P. J.; Block, S.; Piermarini, G. J. Effects of Pressure and Temperature on the Thermal Decomposition Rate and Reaction Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. (30 ) Behrens, R. Thermal Decomposition of Energetic Materials: Temporal Behaviors of the Rates of Formation of the Gaseous Pyrolysis Products from the Condensed Phase Decomposition of Octahydro 1,3,5 7 tetra nitro 1,3,5,7 tetrazocine. J. Phys. Chem. 1990, (31 ) Bulusu, S.; Behrens, R. Thermal Decomposition of Energetic Materials. 2. Deuterium Isotope Effects and Isotopic Scrambling on Condensed Phase Decomposition of Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. Chem. (32 ) Brower, K. R.; Naud, R. L. Pressure Effects on the Thermal Decomposition of Nitramines, Nitrosamines, and Nitrate Esters. J.

PAGE 109

109 (33 ) Williams, G. K.; Gongwer, P. E. ; Brill, T. B. Thermal Decomposition of Energetic Materials. 66. Kinetic Compensation Effects in HMX, RDX, and NTO. J. Phys. Chem. 12247. (34 ) Zheng, W.; Szekeres, R.; Kooh, A. B.; Oxley, J. C. Mechanisms of Nitramine Thermolysis. J. Phys. (35 ) Bharadwaj, R. K.; Smith, G. D. Quantum Chemistry Based Force Field for Simulations of HMX. J. Phys. Chem. B 1999, 103, (36 ) Voth, G. A.; Opdorp, K. V.; Glaesemann, K. R.; Lewis, J. P. Ab Initio Calculations of Octahydro 1,3,5 7 tetranitro 1,3,5,7 HMX). J. Phys. Chem. A 2000, 104, (37 ) Goddard, W. A.; Dasgupta, S.; Muller, R. P.; Chakraborty, D. The Mechanism for Unimolecular Decomposition of RDX (1,3,5 Trinitro 1,3,5 triazine), an ab Initio Study. J. Phys. Chem. A 2001, 105, (38 ) Lewis, J. P. Energetics of intermolecular HONO formation in condensed phase Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine (HMX). Chem. Phys. Lett. 2003, 371, 93. (39 ) Bernstein, E. R.; Greenfield, M.; Guo, Y. Q. Decomposition of nitramine energetic materials in excited electronic states: RDX and HMX. J. Chem. Phys. 2005, 122, (40 ) Bernstein, E. R.; Bhattacharya, A. B.; Greenfield, M.; Guo, Y. Q. On the excited electronic state dissociation of nitramine energetic materials and model systems. J. 154310. (41 ) Burnham, A. K.; Zaug, J. M.; Glascoe, E. A. Pressure Dependent Decomposition Kinetics of the Energetic Material HMX up to 3.6 GPa. J. Phys. Chem. A 2009, 113, (42 ) McCrone, W. C. Cyclotetramethylene Tetranitramine (HMX). Anal. Chem. 1950, 22, (43 ) Cromer, D. T.; Cady, H. H.; Larson, Larson, A. C. The Crystal HMX and a Refin HMX. (44 ) Boutin, H. P.; Choi, C. S. A Study of the Crystal structure of the Cyclotetramethylene Tetranitramine by Neutron Diffraction. Acta Crystallogr. 1970, B26, (45 ) Small, R W. H.; Cobbledick, R. E. The Crystal Structure of the Form of 1,3,5,7 Tetranitro 1,3,5,7 HMX). Acta Crystallogr. 1974, B30,

PAGE 110

110 (46 ) Small, R. W. H.; Cobbledick, R. E.; Main, P. Structure of the Fourth Form of 1,3,5,7 Tetr anitro 1,3,5,7 HMX), 2C4H8N8O8 5H20. Acta (47 ) Thompson, D. L.; Rice, B. M.; Sorescu, D. C. Isothermal Isobaric Molecular Dynamics Simulations of 1,3,5,7 Tetranitro 1,3,5,7 tetraazacyclooctane (H MX) Crystals J. Phys. Chem. B. 1998, 102, (48 (49 ) Thompson, D. L.; Rice, B. M.; Sorescu, D. C. Theoretical Studies of the Hydrostatic Compression of RDX, HMX, HNIW, and PETN Crystals. J. Phys. C hem. B 1999, 103, (50 ) Voth, G. A.; Evans, R. B.; Sewell, T. D.; Lewis, J. P. Electronic Structure Calculation of the Structures and Energies of the Three Pure Polymorphic Forms of Crystalline HMX. J. Phys. Chem. B 2000, 104, (51 ) Lip pert, T. K.; Pulay, P.; et al. Theoretical and Experimental Study of the Vibrational 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. Chem. B 2002, (52 ) Dickson, P. M.; Asay, B. W.; in the energetic nitramine octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine: (53 ) Chabalowsku, C. F.; Scuseria, G. E.; Byrd, E. F. C. An ab Initio S tudy of Solid Nitromethane, HMX, RDX, and CL20: Successes and Failures of DFT. J. Phys. Chem. (54 ) Peiris, S. M.; Gump, J. C. Isothermal equations of state of beta octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine at high temperature s. J. Appl. Phys. 2005, 97, (55 ) Xiao, H.; Zhao, F.; et al. First Principles Study of the Four Polymorphs of Crystalline Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. Chem. B 2007, (56 ) Rice, B. M.; Byrd, E. F. C. Ab Initio Study of Compressed 1,3,5,7 Tetranitro 1,3,5,7 tetraazacyclooctane (HMX), Cyclotrimethylenetrinitramine (RDX), 2,4,6,8,10,12 Hexanitrohexaazaisowurzitane (CL 20), 2,4,6 Trinitro 1,3,5 benzenetriamine (TATB), and Pentaerythritol Tetranitrate ( (57 ) Xiao, H.; Wei, T.; Zhang, X.; Zhu, W. DFT studies of pressure effects on structural and vibrational properties of crystalline octahydro 1,3,5,7 tetranitro 1,3,5,7 tetrazocine. Theor. Chem. Acc. 2009, 124,

PAGE 111

111 (58 ) Liu, Z.; Weck, P.; et al. A far and mid infrared study of HMX (octahydro 1,3,5,7 tetranitro 1,3,5,7 tetrazocine) under high pressure. (59 ) Jacobs, S. J.; Kamlet, M. J. Chemistry of Detonations. I. A Simpl e Method for Calculating Detonation Properties of C H N O Explosives. J. Chem. Phys. 1968, 48, (60 ) Ruggiero, A. J.; Fried, L. E. Energy Transfer Rates in Primary, Secondary, and 9791. (61 ) Kosh i, M.; Ye, S. Theoretical Studies of Energy Transfer Rates of Secondary (62 ) Sheffield, S. A.; Hooks, D. E.; et al. Isentropic loading experiments of a plastic bonded explosive and constituents. J. Appl. (63 ) Dickinson, C.; Rosen, J. M. Vapor Pressures and Heats of Sublimation of Some High Melting Organic Explosives. J. Chem. Eng. (64 cy clotetramethylene tetranitramine and some of its isotopic isomers. J. (65 ) Ferraro, J. R.; Brill, T. B.; Goetz, F. Pressure Dependence of the Raman and Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. (66 ) Marino, R. A.; Brill, T. B.; Landers, A. G. Electronic Effects and Molecular Motion in Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine Based on 14N Nuclear Quadrupole Resonance Sp ectroscopy. J. (67 ) Pace, M. D. EPR Spectra of Photochemical NO2 Formation in Monocyclic Nitramines and Hexanitrohexaazalsowurtzitane. J. Phys. (68 ) Nichols, A. L.; Chidester, S. K.; Tarver, C. M. Critical Conditions for Impact and Shock Induced Hot Spots in Solid (69 ) Ray, S. N.; Das, T. P.; Sahoo, R. P. Electronic Structure Investigation and Nuclear HMX. J. Phys. Chem. (70 ) Eckhardt, C. J.; Stevens, L. L. The elastic constants and related HMX determined by Brillouin scattering. J. Chem. Phys. (71 ) Hall, C. A.; Forbes, J. W.; et al. Isentropic compression of cyclotetramethylene tetranitramine (HMX) single crystals to 50 GPa. J. Appl. Phys. 2006, 99,

PAGE 112

112 (72 ) Molt, R. W.; Watson, T.; Lotrich, V. F.; Bartlett, R. J. RDX Geometries, Excited States, and Revised Energy Ordering of Conformers via MP2 an d CCSD(T) Methodologies: Insights into Decomposition Mechanism. J. Phys. Chem. A 2011, 115, (73 A.; Watson, T.; Bartlett, R. J. Conformers of CL 20 Explosive and ab Initio Refinement Using Perturbation Theory: Implications to Detonation Mechanisms. J. Phys. Chem. A (74 ) Rice, B. M.; Chabalowski, C. F. Ab Initio and Nonlocal Density Functional Study of 1,3,5 Trinitro s triazine (RDX) Conformers. J. (75 ) Goto, N.; Fuj ihisa, H.; Yamawaki, H.; Wakabayashi, K.; Nakayama, Y.; Yoshida, M.; Koshi, M. Crystal Structure of the High Pressure Phase of Hexahydro 1,3,5 trinitro 1,3,5 (76 ) Patterson, J. E.; Dreger, Z. A.; Hupta, Y. M. Shock Wave Induced Phase Transition in RDX Single Crystals. J. Phys. Chem. B (77 ) Davidson, A. J.; Oswald, I. D.; Francis, D. J.; Lennie, A. R.; Marshall, W. G.; Millar, D. I.; Pulham, C. R.; Warren, J. E.; Cumming, A. S. Explosives under pressure the RDX as determined by high pressure X ray and neutron diffraction. (78 ) Oswald, I. D.; Francis, D. J.; Marshall, W. G.; Millar, D. I.; Pulham, C. R.; Warren, J. E.; Cum ming, A. S. The crystal structure of RDX an elusive form of an explosive revealed. Chem. Commun. (79 ) Millar, D. I.; Oswald, I. D.; Barry, C.; Francis, D. J.; Marshall, W. G.; Pulham, C. R.; Cumming, A. S. Pressure cooking of explosives t RDX as determined by X ray and neutron (80 ) Dreger, Z. A.; Gupta, Y. M. Raman Spectroscopy of High Pressure High Temperature Polymorph of Hexahydro 1,3,5 trinitro 1,3,5 RDX). J. Phys. (81 ) Behrens, R. Thermal Decomposition of Energetic Materials: Temporal Behaviors of the Rates of Formation of the Gaseous Products from Condensed Phase Decomposition of Oxtayhydro 1,3,5,7 tetranitro 1,3,5,7 t etrazocine J. Phys. Chem. 1990, 94, 6709 6718. (82) Miller, R. S. In Decomposition, Combustion and Detonation Chemistry of Energetic Materials; Brill, T. B., Russell, T. P., Tao, W. C., Wardle, R. B., Eds.; Materials Research Society Symposium Proceeding s 418; Materials Research Society: Pittsburgh, PA, 1995; p 3. (83) Nielsen, A. T.; Chafin, A. P.; Christian, S.; Moore, D. W.; Nadler, M. P.; Nissan, R. 11812.

PAGE 113

113 (84) Rice, B. M.; Chabalowski, C. F. J. Phys. Ch em. A 1997, 101, (85) Goto, N.; Fujihisa, H.; Yamawaki, H.; Wakabayashi, K.; Nakayama, Y.; Yoshida, M.; Koshi, M. J. Phys. Chem. B. 2006, 110, (86) Patterson, J. E.; Dreger, Z. A.; Hupta, Y. M. J. Phys. Chem. B 2007, 111, 4. (87) Davidson, A. J.; Oswald, I. D.; Francis, D. J.; Lennie, A. R.; Marshall, W. G.; Millar, D. I.; Pulham, C. R.; Warren, J. E.; Cumming, (88) Oswald, I. D.; Francis, D. J.; Marshall, W. G.; Millar, D. I.; Pulham, C. R.; Warren, J. E.; Cumming, A. S. Chem. Commun. 2009, (89) Millar, D. I.; Oswald, I. D.; Barry, C.; Francis, D. J.; Marshall, W. G.; Pulham, C. R.; (90) Dreger, Z. A.; Gupta, Y. M. J. Phys. 7047. (91) Molt, R. W.; Watson, T.; Lotrich, V. F.; Bartlett, R. J. J. Phys. Chem. A 2011, 115, (92) Ashcroft, N. A. Mermin, N. D. Solid State Physics; Thomson Learning: Toronto, Canada, 1976; Chapter 24. (93) Owens, F. 5444. Fanelli, A. J.; (9 (97) Oxley, J. C.; Kooh, A. B.; Zheng, W. J. Phys. Chem. 1994, 98, 169. (99) Sorescu, D. R.; Rice, B. M.; Thompson, D. L. J. Ph ys. Chem. B 1998, 102, (100) Sorescu, D. R.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. B 1998, 102, (101) Byrd, E. F. C.; Scuseria, G. E.; Chabalowski, C. F. J. Phys. Chem. B 2004, 108, (102)Lotrich, V.; Flocke, N.; Ponto n, M.; Yau, A.; Perera, A.; Deumens, E.; Bartlett, R. J. Parallel implementation of electronic structure energy, gradient, and Hessian calculations. J. Chem. Phys. 2008 128 194104 194104 15.

PAGE 114

114 (103)Truhlar, D.G.; Leverentz, H.R., et al. Perspectives on Basis Sets Beautiful: Seasonal Plantings of Diffuse Basis Functions. J. Chem. Theory. Comput. 2011 7, 3027 3034. (104)Dunning, T.H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Che m. Phys. 1989 90 1007 1023. (105)Krishnan, R.; Binkley, J.; Seeger, R.; Pople, J. A. Self consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980 72 650 654. (106)T., C.; Chandrasekhar, J.; Schleyer, P. V. R. Efficient diffuse function augmented basis sets for anion calculations. III. The 3 21+G basis sets for first row elements, Li F. J. Comp. Chem. 1983 4 294 301. (107)Burden, R. and Faires, J. Numerical Analysis 5th Edition, 1993, p.56. (108)Adamo, C.; Barone, V. Exchange functionals with improved long range behavior and adiabatic connection methods without adjustable parameters: The mPW and mPW1PW models. J. Chem. Phys. 1998 108 664 675. (109)Bytheway, I.; Wong, M.W. The prediction of vibrational frequencies of inorganic molecules using density functional theory. Chem. Phys. Lett. 1998 282 219 226. (110)Arfken, G.B.; Weber, H.J. Mathematical Methods for Physicists 6th ed.; Elsevier: New York, 2005. (111)Hassani, S. Mathematical Physics ; Springer: New York, 1999. (112)Vladimiro ff T.; Rice, B.M. Reinvestigation of the Gas Phase Structure of RDX Using Density Functional Theory Predictions of Electron Scattering Intensities. J. Phys. Chem. A 2002 106 10437 10443. (113) Stoltz, C.; Mason, B. P.; Hooper, J. J. Chem. Phys. 2010, 107,103527. (114) Torres, P.; Mercado, L.; Cotte, I.; Hernandez, S. P.; Mina, N.;Santana, A.; Chamberlain, R. T.; Lareau, R.; Castro, M. E. J. Phys. Chem. B 2004, 108, 8799 8805. (115) Wu, C. J.; Fried, L E. J. Phys. Chem. A 1997, 101, 8675 8679. (116) Wade, L. G. Organic Chemistry, 5th ed.; Prentice Hall: NewYork, 2003. (117) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 1997, 106, 1063. (11 8 ) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 2000, 112, 7374.

PAGE 115

115 (119 ) Wang, N. X.; Wilson, A. K. J. Chem. Phys. 2004, 121, 7632. (120) Bartlett, R. J.; Lotrich, V. F.; Schweigert, I. V. J. Chem. Phys. 2005, 123, 062205. (121) Bartlett, R J.; Schweigert, I. V.; Lotrich, V. F. THEOCHEM 2006, 771, 1 8. (122) Bartlett, R. J.; Musial, M. Rev. Mod. Phys. 2007, 79, 291 352. (123) Coester, F.; Kummel, H. Nucl. Phys. 1960, 17 477. (124) Bartlett, R. J.; Stanton, J. Rev. Comput. Chem. 1994, 5, 65. (125) Urban, M.; Noga, J.; Cole, S.; Bartlett, R. J. J. Chem. Phys. 1985, 83, 4041. (126) Helgaker, T.; Gauss, J.; Jrgensen, P.; Olsen, J. J. Chem. Phys. 1997, 106, 6430 6440. (127) Krishnan, R.; Binkley, J.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (128) Bak, K. L.; Jrgensen, P.; Olsen, J.; Helgaker, T.; Klopper, W.J. Chem. Phys. 2000, 112, 9229. (129) Stanton, J.; Bartlett, R. J. Chem. Phys. 1993, 98, 7029. (130) Wiberg, K. B.; Hadad, C. M.; Foresman, J. B.; Chupka, W. A. J. Phys. Chem. 1992, 96, 10756 10768. (131) Wiberg, K.; Oliveira, A.; Trucks, G. J. Phys. Chem. A 2002, 106, 4192 4199. (132) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson G. A.; Montgomery, J. A.; Raghavachari, K.; Al Laham, M. A.; Zakrzewski, V. G.; Ortiz, J.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S .; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, revision b.1; Gaussian, Inc.: Pittsburgh, PA, 1995.

PAGE 116

116 (133) The Gaussian 09 manual has ch anged the default setting for numeric integration, in light of previous failings: Frisch, M. J.; Trucks,G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.;Caricato M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cr oss, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S .; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09, revision a.1; Gaussian, Inc.: Wallingford, CT, 2009. (134) Karpowicz, R.; Brill, T. J. Phys. Chem. 1984, 88, 348. (135) Im, H. S.; Bernstein, E. R. J. Ch em. Phys. 2000, 113, 7911 7918. (13 6 ) Cheron, N.; Denis Jacquemin, D.; Fleurat Lessard, P. Aqualitative failure of (137) Goddard, W. A.; Dasgupta, S.; et al. Mechanism for Unimolecular Decomposition of HMX (1,3,5,7 Tetranitro 1,3,5,7 tetraazocine), an ab Initio Study. J. Phys. Chem. A (138) Lippert, T. K.; Pulay, P.; et al. Theoretical and Experimental Study of the Octahydro 1,3,5 7 tetranitro 1,3,5,7 tetrazocine. J. Phys. Chem. B 2002, 106, (139) Kohn, W.; Sham, L. J. Self Consistent Equations Including Exchange and (140) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, (141) Cramer, C. Essentials of Computational Chemistry, 2nd ed.; John Wiley & Sons Ltd: West Sussex, U.K., 2008; Chapter 8. (142) Martin, R. M. Electronic Structure; Cambridge Uni versity Press:Cambridge, U.K., 2004.

PAGE 117

117 (143) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation 6. (144) MartinParr, R.G.; Tao, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, U.K., 1994; Chapter 6. (145) Vosko, S. H.; Wilk, L.; Nusair, M. Accurate spin dependent electron liquid correlation energies for local s pin density calculations: a critical analysis. Can. J. Phys. (146) Vela, A.; Trickey, S. B.; et al. Non empirical improvement of PBE and its hybrid PBE0 for general description of molecular properties. J. Chem. Phys. 2012, 136, 104108 (147) Zhao, Y.; Truhlar., D. G. The M06 suite of density functional for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06 class (148) Becke, A. D. Density functional thermochemistry. III. The role of exact exchange. (149) Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. (150) Yanai, T.; Tew, D.; Handy., N. A new hybrid exchange corr elation functional using the Coulomb (151) Grimme., S. Semiempirical GGA type density functional constructed with a long (1 5 2) Pople, J. A.; Head Gordon, M.; Raghavachari, K. Quadratic configuration interaction. A general technique for determining electron correlation energies. J. Chem. (1 53 ) Bartlett, R. J.; Musial, M. Coupled cluster theory in quantum chemistry. Rev. Mod. (1 54 ) Bartlett, R. J.; Shavitt, I. Many Body Methods in Chemistry and Physics;

PAGE 118

118 (1 55 ) Bartlett, R. J.; Shavitt, I. Many Body Methods in Chemistry and Physics; Cambridge University Press: Cambridge, U.K., 2009; Chapters (1 56 ) Truhlar, D. G.; Leverentz, H. R.; et al. Perspectives on Basis Sets Beautiful: Seasonal Plantings of Diffuse Basis Functions. J. Chem. Theory. Comput. 2011, 7, (1 57 ) Krishnan, R.; Binkley, J.; Seeger, R.; Pople, J. A. Self consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72, (1 58 ) Clark., T.; Chandrasekhar, J.; Spitznagel, G. W.; Schle yer, P. v.R. Efficient diffuse function augmented basis sets for anion calculations. III. The 3 21+G basis sets for first row elements, Li (1 59 ) Zhou, G.; Wang, J.; He, W.; Wong, N.; Tian, A.; Li, W. THEOCHEM 2002, 27 (1 60 ) Zhang, J.; Du, H.; Wang, F.; Gong, Z.; Huang, Y. J. Phys. Chem. A 2011, 115, (1 61 ) Okovytyy, S.; Kholod, Y.; Qasim, M.; Frederickson, H.; Leszczynski, J. J. Phys. (1 62 ) For a compilation of informati on about HNIW, see: Filliben, J. D., Ed. CPIA/M3 Solid Propellant Ingredients Manual; Publication Unit 89; Chemical Propulsion Information Agency: Columbia, MD. Distribution authorized to U.S. Government agencies only; test and evaluation, Sep 1994. Other requests for this document shall be referred to the Office of Naval Research, Code ONR 35, Arlington, VA 22217 5000. (1 63 ) Chan, M. L.; Carpenter, P.; Hollins, R.; Nadler, M.; Nielsen, A. T.; Nissan, R.; Vanderah, D. J.; Yee, R.; Gilardi, R. D. CPIA PUB 6 25, Apr 95, 17P, CPIA Abstract No. X95 07119, AD D606 761. Availability: Distribution authorized to the Department of Defenseand DoD contractors only; Critical Technology, Mar 17, 1995. Other requests shall be referred to the Naval Air Warfare Center Weapo ns Division (474000D), China Lake, CA 93555 6001. This document contains export controlled technical data. (1 64 ) Russell, T. P.; Miller, P. J.; Piermarini, G. J.; Block, S.; Gilardi, R.; George, C. AD C048 931 (92 0134), p 155, Apr 91, CPIA Abstract No. 9 2, 0149, AD D604 542, C D; Chemical Propulsion Information Agency: Columbia, MD. (1 65 ) Foltz, M. F.; Coon, C. L.; Garcia, F.; Nichols, A. L., III. ADC049 633L (93 0001), p 9, Apr 92, Contract W 7405 ENG 48, CPIA Abstract No. 93 0003, AD D605 199, U D; Che mical Propulsion Information Agency: Columbia, MD.

PAGE 119

119 (1 66 ) Bolotina, N. B.; Hardie, M. J.; Speer, R. L.; Pinkerton, A. A. J. Appl. Crystallogr. (1 67 ) Hirata, S., Yanai, T., Jong, W.A., et al J. Third order Douglas Kroll relativistic coupled cluster theory through connected single, double, triple, and quadruple substitutions: Applications to diatomic and triatomic hydrides. J. Chem. Phys. 120 (2004), 3297. (1 68 ) Sherill, C.D., Leininger, M.L, Van Huis, T. J., and Schaefer, H.F. Struct ures and vibrational frequencies in the full configuration interaction limit: Predictions for four electronic states of methylene using a triple zeta plus double polarization (TZ2P) basis. J. Chem. Phys. 108 (1998), 1040. (1 69 ) Slipchenko, L.V. and Krylov A.I. Spin conserving and spin flipping equation of motion coupled cluster method with triple excitations. J. Chem. Phys. 123 (2005), 084107. (1 70 ) Slipchenko, L.V. and Krylov, A.I. Singlet triplet gaps in diradicals by the spin flip approach: A benchmar k study. J. Chem. Phys. 117 (2002), 4694. (1 71 ) Li, X. and Paldus, J. Electronic structure of organic diradicals: Evaluation of the performance of coupled cluster methods. J. Chem. Phys. 129 (2008), 174101. (1 72 ) McKellar, A.R.W., Bunker, P.R., Sears, T. J., et al. Far infrared laser magnetic resonance of singlet methylene:Singlet triplet perturbations, singlet triplet transitions, and the singlet triplet splittinga). J. Chem. Phys. 79 (1983), 5251. (1 73 ) Davidson, E.R., Feller, D., and Phillips, P. Re lat ivistic Corrections for Methy lene. Chem. Phys. Lett. 76 (1980), 416. (1 74 ) Bauschlicher, C.W. And Shavitt, I. Accurate ab Initio Calculations on the Singlet Triplet Separations in Methylene. J. Am. Chem. Soc. 100 (1978), 739. (1 75 ) Demel, O. and Pittner, J. Multireference Brillouin Wigner coupled cluster method with singles, doubles, and triples: Efficient implementations and comparison with approximate approaches. J. Chem. Phys. 128 (2008), 104108. (176) Balkova, A. and Bartlett, R.J. On the singlet tri plet separation in methylene: A critical comparison of single versus twodeterminant (generalized valence bond) coupled cluster theory. J. Chem. Phys. 102 (1995), 7116. (177) P. Jensen and P. R. Bunker, J. Chem. Phys. 89 1327 1988. (178 ) N. C. Handy, Y. Yamaguchi, and H. F. Schaefer, III, J. Chem. Phys. 84 ,4481 1986.

PAGE 120

120 (179 ) Shen, J and Piecuch, P. Merging Active Space and Renormalized Coupled Cluster Methods via the CC(P;Q) Formalism, with Benchmark Calculations for Singlet Triplet Gaps in Biradical Systems. J. Chem. Theory. Comp. 8 (2012), 4968. (180 ) C. J. Cramer and B. A. Smith, J. Phys. Chem. 100 9664 (1996). (181 ) P. G. Wenthold, J. Hu, R. R. Squires, and W. C. Lineberger, J. Am. Chem. Soc. 118 475 (1996). (18 2) Slipchenko, L.V. and Kryl ov, A. Electronic structure of the trimethylenemethane diradical in its ground and electronically excited states: Bonding, equilibrium geometries, and vibrational frequencies. J. Chem. Phys. 118 (2003), 6874. (183 ) Crawford, T.D., Kraka, E., Stanton, J. F., and Cremer, D. Problematic p benzyne: Orbital instabilities, biradical character, and broken symmetry. J. Chem. Phys. 114 (2001), 10638. (184 ) Wang, E.B., Parish, C.A., and Lischka, H. An extended multireference study of the electronic states of p ara benzyne. J. Chem. Phys. 129 (2008) 044306. (185 ) Marquardt, R., Balster, A., Sander, W., et al. p Benzyne. Angew. Chem. Int. Ed. 37 (1998). 955. (186 ) Leopold, D.G., Miller, A., Lineberger, W.C. Determination of the Singlet Triplet Splitting and Electron Affinity of o Benzyne by Negative Ion Photoelectron Spectroscopy. J. Am. Chem. Soc. 86 (1986) 1379. (187 ) P. G. Wenthold, R. R. Squires, and W. C. Lineberger, J. Am. Chem. Soc. 120 5279 1998. (188) Private communication with Dr. Wes Allen. (189) F. A. Evangelista, W. D. Allen, and H. F. Schaefer III, J. Chem. Phys. 127 024102 ( 2007).

PAGE 121

121 BIOGRAPHICAL SKETCH Robert William Molt Junior grew up in the yahoo town of Northbridge, Massachusetts. He earned a Bachelor of Arts in a liberal arts education at The College of the Holy Cross despite the fact that he is not a Christian. He had the pleasure of studying chemistry and physics, prima rily. In his sophomore year, he learned about the many body problem, and was endlessly fascinated. His research work focused on experimental issues of solid state physics (magnetic hysteresis and eddy currents in steels), computational electromagnetic s ( electric field homogeneity in the presence of surface defects of parallel plate capacitors), and computational atomic physics (calculation of oscillator strengths and energies of Rydberg atoms using various model potentials) He gratefully thanks Drs. Ded ra Demaree and Paul Oxley for being his mentors. Robert then began graduate school at the University of Minnesota, working under Drs. Truhlar and Cramer. He enjoyed his coursework in the physics department, and would like to especially thank Drs. Vuk Man dic and Joseph Kapusta for excellent education in mathematical methods in physics and electromagnetic s is forever haunted by John David Jackson as a co nsequence. Robert left with a Master of Science to pursue the study of many body methodol ogies without the use of Dr. Bartlett foolishly took on Robert as his graduate student at the University of Florida. In between research, he managed the Sanibel conference for three years, defeated can cer, and rescued cats. He received his Ph.D. from the University of Florida in December 2013.