ARMY RESEARCH OFFICE Military University Research Initiative Oct 16, 2003

1. ARMY RESEARCH OFFICEMilitary University Research Initiative Oct 16, 2003 Rodney J. Bartlett Co-Workers Mr. Andrew Taube Mr. Josh McClellan Mr. Tom Hughes Mr.Luis Galiano Dr. Stefan Fau Dr. DeCarlos Taylor (ARL) Dr. Ariana Beste (ORNL) Quantum Theory Project Departments of Chemistry and Physics University of Florida Gainesville, Florida USA

2. QTP OBJECTIVES Identify and characterize the initial steps in nitramine and other detonation in the condensed phase. Progress requries NEW ab initio quantum mechancial techniques that have the accuracy and appicability to provide reliable results for unimolecular and bimolecular reaction paths. Study the series of molecules, nitramine (gas phase), methyl nitramine(liquid), dimethylnitramine(solid) which have (1) different reaction paths (2) different condensed phase effects Investigate their unimolecular, secondary, and bimolecular reaction mechanisms. Obtain definitive results for the comparative activation barriers for different unimolecular paths, particularly for RDX. Study the nitromethane molecule and its various isomers as a prototype for nitroalkanes. Generate ‘transfer Hamiltonians’ to enable direct dynamics simulations as a QM compliment to classical potentials, and to be able to reliably describe many units of a condensed phase explosive. Provide high-level QM results to facilitate the development of classical PES for large scale simulations. University of Florida: Quantum Theory Project

3. Quantum Mechanics I (Isolated gas phase molecules, 0K) Potential Energy Surface E(R) Different Unimolecular Decomposition Paths Activation Barriers Spectroscopic signatures for intermediates and products Quantum Mechanics II Bi(tri...)molecular reactions Long range (condensed phase, pressure) effects Activation Barriers, Spectroscopy Large Molecule QM-- Simplified Representation of H(R) Transfer Hamiltonian Electronic State Specific Classical Mechanics-- Representation E(R)

4. Reactions of one Me2N-NO2 Me2N-NO2 Me2N. + NO2.  MeN.-NO2 + Me.  H2C.-N(Me)-NO2 + H. [ H2C=N+(Me)-N(O-)OH]  H2C=NMe + HONO

5. Reactions of two Me2N-NO2 Same reactions as before, causing a slight change in the interaction energy with the second Me2N-NO2. Additionally: Me2N-NO2 + Me2N-NO2 Me2N-ONO + Me2N. + NO2.  Me2N-NMe2 + NO2. + NO2.  Me2N-Me + Me-N.-NO2 + NO2.  Me2N-N(Me)-NO2 + Me. + NO2.  Me2N-H + H2C.-N(Me)-NO2 + NO2.  Me2N-CH2-N(Me)-NO2 + H. + NO2. products from CH3., H.(HONO, H-Me, …)

6. In reposne to Bob Shaw’s comment about error bars in theoretical applications…. It is clear that if we want to know the right answer for activation barriers competitive decompisiton paths.. We need a high-level of theory like CCSD(T), with a large basis set.

7. MP2 CCSD +(T) D 6.0 -8.3 -1.0 Dstd 7.5 4.5 0.5 Coupled Cluster Calculation of De’s r(De) De (kcal/mol) From K. L. Bak et al., J. Chem. Phys. 112, 9229-9242 (2000)

8. Last year… • Reported the detailed theory for the compressed (SVD) CC, which ‘contracts’ the CC amplitudes in an optimum way to make it possible to perform much higher level CC calculations for large molecules.

9. 0 -50 -100 (E+100)*1000 (a.u.) -150 MRCI CCSD -200 CCSD(T) CCSDT-1 COMP.SDT-1 -250 -300 0.5 1 1.5 2 2.5 3 3.5 4 r(HF)/r(eq) Potenrial Energy Curve (1) (HF, aug-cc-pVDZ, HF-bond stretcing) r(eq)=1.733 bohr a=0.25

10. 100 0 (E+76)*1000 (a.u.) -100 MRCI CCSD CCSD(T) -200 COMP.SDT-1 CCSDT-1 -300 0.5 1 1.5 2 2.5 3 3.5 r(OH)/r(eq) Potential Energy Curve (2) (H2O, aug-cc-pVDZ, OH-bonds stretching) r(eq)=1.809 bohr a=0.25

11. Detailed a new extrapolation procedure for energies and forces that • has a mean error of nearly zero, and a maximum error of • 0.75 kcal/mol for nitramne and its components. • Reported on a series of CC studies of nitramine to assess its • decomposition paths. • This year….

12. QTP OUTLINE I.A NEW APPROACHE TO HIGH-LEVEL COUPLED-CLUSTER THEORY FOR LARGER MOLECULES. A. A natural orbital coupled-cluster method. B. Comparisons toCompressed coupled-cluster theory. II. NUMERICAL ILLUSTRATIONS FOR DMNA, DMNA DIMERS AND RDX. Iii. CONCEPT OF A TRANSFER HAMILTONIAN AS A MEANS TO DESCRIBE COMPLEX SYSTEMS WITH QUANTUM MECHANCIAL FORCES FOR MD APPLICATIONS. A.Illustration for nitromethane and its isomers. University of Florida: Quantum Theory Project

13. How can we retain the accuracy of CCSD(T), ie an ~n7 method, but make it applicable to large molecules? Dimer (bimolecular sytem) is 28 times as difficult as the monomer, without modification. I. Natural orbital coupled-cluster theory II. Compressed coupled-cluster using Singular Value Decomposition.

14. Frozen Natural Orbital Coupled Cluster Theory • Dependence on size of virtual sector of basis set limits high-level CC to small molecules: • CCSD ~ V4; CCSD(T) ~ V4; CCSDT ~ V5 • Natural Orbitals (NOs) are known to be the best possible set of orbitals to truncate • Too costly to get exact – use approximate MBPT(2) NOs • To maintain advantages of a HF reference, only perform truncation in virtual space – leaving occupied space alone – Frozen Natural Orbitals (FNOs)

15. FNO Procedure • SCF gives Molecular Orbitals U • Construct MBPT(2) Density Matrix D in MO Basis • Solve DV=Vn • Truncate V to V’ throwing out less occupied virtuals, as measured by their occupation numbers. • Construct new Fock Matrix in FNO Virtual Space • Diagonalize for new orbital energies • Perform higher level (CC) calculation in truncated virtual space • For an estimate of truncation error, define: ∆MBPT(2) = MBPT(2) (Full) – MBPT(2) (Truncated)

16. Note: On the following slides, the percentage listed indicated the number of FNOs retained. For example, 20% DZP means that there are only 20% of the original number of DZP virtual orbitals left. Computational Details: Calculations were performed on an IBM RS/6000 375 MHz POWER3 processor with 3 GB RAM and 18 GB disk using the ACES II electronic structure program.

17. Three Different Bases. Same Number of Orbitals. FNO RHF CCSD(T) PES for Hydrogen Fluoride Reference

18. Average Timings for Determination of CCSD(T) PES for Six Small Molecules* with FNO Truncation *F2, HF, CO, N2, NH3, H2O – PES of symmetric dissociation *Relative to 20% DZP Basis Calculation

19. Truncated Large Bases are Better than Similar-sized Small Bases Points are in increments of 20% of basis set * Total Correlation Energy Determined by 100% cc-pVQZ calculation ** Times Relative to 20% DZP Calculation

20. Dimethylnitramine • Model compound for RDX and HMX. • Two possible conformers: • C2v – X-ray crystallography / gas-phase e- diffraction1; • Cs – predicted by theory at SCF, MBPT(2), DFT, QCISD levels with moderate basis sets2. It has been argued that theoretical predictions are more accurate than experimental values. • Dimer interactions are important to model dominant interaction in the solid phase. Given X-ray structure, best to use C2v monomer. Calculations have been done at SAPT level on fixed monomers3. • DMNA can undergo decomposition via NO2 and HONO elimination. HONO Elimination estimated to be exothermic by ~1-3 kcal/mol at standard conditions4, therefore need high-level calculations, zero-point energy corrections, etc.These have been investigated theoretically locating transition states at the QCISD, MBPT(2) and DFT levels2. 1) Stolevik, Rademacher, Acta Chem Scand196923 672 2) Harris, Lammertsma, JPCA1997101 1370; Smith, et al., JPCB1999103 705; Johnson, Truong, JPCA1999103 8840 3) Bukowski, Szalewicz, Chabalowski, JPCA1999103 7322 4) Shaw, Walker, JPC197781 2572

21. CCSD(T) RHF Drop Core C2v Cs C2v – Cs (kcal/mol) DZP – 100% -338.8157 -338.8302 9.09 DZP – 100% Time 0.39 hr 1.22 hr DZP – 60% -338.7549 -338.7671 7.66 DZP – 60% Time 0.06 hr 0.19 hr DZP – 60% + ∆MBPT(2) -338.8078 -338.8220 8.91 cc-pVTZ – 100% -339.1252 -339.1399 9.24 cc-pVTZ – 100% Time 13.1 hr 39.5 hr cc-pVTZ – 60% -339.1032 -339.1180 9.29 cc-pVTZ – 60% Time 1.80 hr 7.00 hr cc-pVTZ – 60% + ∆MBPT(2) -339.1283 -339.1431 9.29 Cs** DMNA Equilibrium Structure C2v* Total Energies in Hartree *Experimental geometry with HCH angles optimized at MBPT(2) Bukowski et.al. JPCA1999103 7322 **Theoretical prediction QCISD cc-pVDZ basis Johnson & Truong JPCA1999103 8840

22. DMNA Dimer Structures† - fixed monomer geometry M2 M1 M3 M4 † Monomer geometries are experimental methyl angles, optimized at MBPT(2) level. Dimer structures are minima of SAPT method – Bukowski, Szalewicz, Chabalowski JPCA1999103 7322

23. DMNA Dimer Interaction Energies † * * CCSD(T) FNO 60% DZP with core occupieds dropped and ∆MBPT(2) correction † Bukowski, Szalewicz, Chabalowski JPCA1999103 7322

24. Background information for RDX Decomposition Mechanisms: • NO2 Elimination: B3LYP/6-3311G** (1), DSC closed pan(liquid-like) (2) • HONO Elimination (3) • Triple Bond fission: IRMPD(4), DSC open pan(gas-like)(2) • Other mechanisms: NO elimination (MALDI)(5), internal ring formation, -OH loss Conformers of RDX: • Solid α-RDX: AAE (Cs) (6) • Solid β-RDX: AAA • Vapor phase(e- diffraction): AAA(C3v) • Gas/liquid dynamically averaged structure(7) The energy difference between minima is on the order of 1 Kcal/mol (B3LYP/6-311G**) with AAE being the most stable(1) 1)N. Harris, K. Lammertsma., J. Am. Chem. Soc. 119,6583 (1997) 2) G. Long, S. Vyazovkin, B. A. Brems, C.A. Wight, J. Phys. Chem. B., 104, 2570 (2000) 3) D. Chakraborty, R.P. Muller, S. Dasgupta, W. Goddard III, J. Phys. Chem. A., 104,2261 (2000) 4) X. Zhao, E. Hintsa, Y. Lee, J. Chem. Phys., 88, 2 ,801 (1988) 5) H.S. Im and E.R. Bernstein. J.Chem.Phys., 113, 18 ,7911 (2000) 6)B. Rice, C. Chabalowski, J. Phys. Chem. A., 101, 8720 (1997) 7) T. Vladimiroff, B. Rice J. Phys. Chem. A., 106, 10437 (2002)

25. CCSD(T) RHF Drop Core AAA Chair AAA Boat Chair-Boat(kcal/mol) DZP – 60% -895.1732 -895.1728 -0.25 DZP – 60% + ∆MBPT(2) -895.3311 -895.3289 -1.4 CC Time (hr) 30 145 Estimated Time for Full Basis (hr) 200 1000 RDX Minima* Total Energies in Hartree AAA Boat AAA Chair *Conformations determined by B3LYP 6-31G(d) Calculations – Chakraborty, et al JPCA2004104 2261

26. Triples at a Fraction of the Cost Compressed Coupled Cluster & FNO CC for DMNA CCSD(T) FNO Timings for AAA Chair Conformer of RDX CCSDT-1 RHF DZ Basis Frozen Core: 6.3 *FNO Speed-up ~ 8x faster with 50% truncation. For large numbers of occupied orbitals, FNO speed-up determined by o3v3 term in CCSD equations – Not by o2v4 term ** Estimated.

27. Transfer Hamiltonian for large clusters of molecules. Basic idea: Represent the CC Hamitonian in its one-particle form by a low-rank operator, that permits rapid generation of forces for MD, but can (hopefully) retain the accuracy of CC theory, in the process.

28. ITR TRANSFER HAMILTONIAN In CC theory we have the equations… exp(-T) Hexp(T) = Ĥ Ĥ|0 = E|0WhereEis the exact correlated energy m |Ĥ|0 =0Wherem|is a single, double, triple, etc excitation which provides the equations for the coefficients inT, ie tia, tijab, etc. (R)E(R) = F(R)Provides the exact forces (x)=0| exp(-T)(x-x’)exp(T) |0 gives the exact density andm| Ĥ|n  Ĥ andĤRk = kRkGives the excitation (ionization, electron attached) energieskand eigenvectorsRk University of Florida: Quantum Theory Project

29. COMPARATIVE APPLICABILITY OF METHODS 10-10 CP 10-8 TB 10-6 SE TH COST 10-4 DFT 10-2 CC 1 ACCURACY

30. TRANSITION FROM MANY-PARTICLE HAMILTONIAN TO EFFECTIVE ONE-PARTICLE HAMILTONIAN... Wavefunction Approach 0|{i†a}Ĥ|0=0= a| Ĝ |i=0 Ĝ|i=i|i i Parameterize Ĝ with a GA to satisfy E= 0|Ĥ|0, E=F(R), (r), (Fermi) = I Density Functional Approach Ĝ|i=i|i i where Ĝ =t+E/(x) and E[]=E, E=F(R), (r)= |i i|, (Fermi) = I Future? Remove orbital dependence and/or self-consistency?

31. ITR RELATIONSHIP BETWEEN COUPLED-CLUSTER/DFT HAMILTONIAN AND SIMPLIFIED THEORY Second QuantizedĜ Ĝ =gpq{p†q} +ZAZB /RAB Transition from orbital based to atom based--  (hAA + AA)+hAB(R)+  AB(R) + Z’AZ’B /RAB{akAexp[-bkA(RAB-ckA)2] +akBexp[-bkB(RAB- ckB)2]} hAB(R)= (+ )KS(R) AB(R)= [( RAB)2 +0.25(1/AA+1/BB)2]-1/2 University of Florida: Quantum Theory Project

32. Nitromethane Background nitromethane (NMT) methylnitrite (MNT) aci-nitromethane • Proper description of NMT unimolecular rearrangement is required for • adequate description of combustion, detonation and pollution chemistry. • NMT is a model system for energetic materials, e.g. FOX-7, TNAZ. • NMT→ CH3∙ +NO2∙ most energetically favored, ~63 kcal/mol*. • NMT→MNT→ CH3O∙ +NO∙ second most energetically favored, ~69 kcal/mol*. • Molecular beam experiments** demonstrate NMT → MNT *CCSD(T)/cc-PVTZ, Nguyen et al., J. Phys. Chem. A 2003, 107, 4286 **Wodtke et al., J. Chem. Phys. 1986, 90, 3549

33. Unimolecular Decomposition Pathways of NMT and MNT CH3-NO2 CH3. + NO2.  TS1  CH3ONO Rearrangement CH3ONO  CH3O . + NO.  TS2  CH2O + HNO Rearrangement • G2MP2* • CCSD/TZP • CCSD(T)/cc-PVTZ** Energies relative to NMT * Hu et al., J. Phys. Chem. A 2002, 106, 7294 ** Nguyen et al., J. Phys. Chem. A 2003, 107, 4286

34. Nitromethane HT

35. Nitromethane PES for C-N rupture

36. Reference Data

37. Nitromethane Clusters • Our nitromethane dimer and trimer calculations used local minima found by Li, Zhao, and Jing in their application of BSSE corrected DFT/B3LYP with 6-31++G** basis*. This reference also concludes that a proper description of three body effects is needed for accurate determination of potential energy surfaces for bulk nitromethane. • The most stable dimer configuration was that in which two hydrogen bonds of length 2.427Å are formed while trimer involved a ring structure in which the three hydrogen bonds of lengths 2.329Å, 2.313Å, and 2.351Å. This configuration for the dimer minimum is also supported by CP corrected SDQ-MBPT/DZP in which the hydrogen bonding distance is found to be 2.25Å**. • With TH-CCSD we observe the formation of methoxy radical in the dimer and the formation of methylnitrite in the trimer, similar to predicted unimolecular mechanisms found at the G2MP2/B3LYP/6-311++G(2d,2p) level of theory ***. • * J. Li, F. Zhao, and F. Jing, JCC 24 (2003) 345. • ** S. J. Cole, K. Szalewicz, G. D. Purvis III, and R. J. Bartlett, JCP 12 (1986) 6833. • S. J. Cole, K. Szalewicz, and R. J. Bartlett, IJQC (1986) 695. • *** W.F. Hu, T.J. He, D.M. Chen, and F.C. Liu, JPCA 106 (2002) 7294.

38. Breaking of H-Bonds in Nitromethane Dimers with Frozen Monomers SAPT equilibrium OH distance*

39. Nitromethane Decomposition in Dimer H3CNO2* + H3CNO2 H3CNO + H3CO + NO2 UHF TH-CCSD predicts that rearrangement occurs when C-N bond is 2.42 Å UHF AM1 prediction * Indicates bond rupture

40. Nitromethane Rearrangement in Trimer H3CNO2* + 2 H3CNO2 3H3CONO UHF TH-CCSD predicts that rearrangement occurs when C-N bond is 1.94 Å UHF AM1 prediction * Indicates bond rupture

41. Equilibrium Geom of NMT

42. TNAZ C-N Bond Rupture *C-N bond breaking trans to N-NO2

43. SUMMARY • Derived and illustrated the NO-CC methodwith applications to • DMNA and RDX. Savings is ~7 out of theoretical 16. With • rapid processors, will make state-of-the-art CC resutls possible in • a week, instead of a year. • Awaiting analytical gradients (programmed) to enable geometry • and transition state searches. • Illustrated tranfer Hamiltonian approach to retain the accuracy of • CC theory, but for much more complicated representations of the • condensed phase. • Establsihed rigor of the theory. We are working on alternative, • and better realizations of the concept.

44. Nitramine: Dangers of DFT ⇒ Forces from DFT are qualitatively and quantitatively wrong at non-equilibrium geometries!