Exploring Potential Energy Surfaces for Chemical Reactions

1. UC - Davis, March 14, 2007 Exploring Potential Energy Surfaces for Chemical Reactions Prof. H. Bernhard Schlegel Department of Chemistry Wayne State University Current Research Group Dr. Jason Sonnenberg Dr. Peng Tao Barbara Munk Jia Zhou Michael Cato Jason Sonk Brian Psciuk Recent Group Members Dr. Xiaosong Li Dr. Hrant Hratchian Dr. Stan Smith Dr. Jie (Jessy) Li Dr. Smriti Anand Dr. John Knox

2. Research overview • With molecular orbital calculations, it is possible to investigate details of chemical reactions and molecular properties that are often difficult to study experimentally • Our group is involved in both the development and the application of new methods in ab initio molecular orbital (MO) methods.

3. Features of Potential Energy Surfaces

4. Development of new algorithms • energy derivatives for geometry optimization • searching for transition states • following reaction paths • computing classical trajectories for molecular dynamics directly from the MO calculations. • spin projection methods to obtain more accurate energetics for open shell systems (radicals) • simultaneous optimization of the wavefunction and the geometry

5. Applications Organic systems Inorganic systems Biochemistry Study of materials Dynamics • Nickel catalyzed three component couplings • Reactions involving nitric oxides • Organo-metallic complexes • Interactions in the active site of enzymes • Guanine oxidation • Amber parameters for modified RNA and DNA bases • CVD studies on TiN and ZnO • Organic LED materials • Molecules in a nanotube • Blue shifted hydrogen bonds • One transition state serving two mechanisms • Molecules in intense laser fields • Two and three body photo dissociation reactions

6. Ni Catalyzed Three Component Coupling The mechanism for a family of nickel catalyzed three component coupling reactions has been studied experimentally by Prof. Montgomery (WSU). MO studies provide additional insight into the mechanism. Montgomery, Acc. Chem. Res.2000, 33, 467-473. Hratchian, Chowdhury, Gutierrez-Garcia, Amarasinghe, Heeg, Schlegel, Montgomery, Organomet.2004, 23, 4636-4646, 5652.

7. Reaction Profile for L=H2NCH2CH2NH2

8. Ligand Exchange with ZnMe2

9. OXA-10 β-lactamase • X-ray structure shows carboxylated Lys70 • Modified Lys70 has mechanistic role • Removes proton from Ser67 • Leads to acylation of Ser67 by substrate • Enzyme shows biphasic kinetics during substrate turnover J. Li; J. B. Cross; T. Vreven; S. O. Meroueh; S. Mobashery; H. B. Schlegel; Proteins 2005, 61, 246-257

10. ONIOM QM/MM Method • The active site region is treated using high-level molecular orbital theory, while the most distant parts of the enzyme are treated using low-cost molecular mechanics.

11. Carboxylation in Gas Phase and in Solution (B3LYP/6-31G(d,p)) Rx1: CH3NH2+CO2 Rx2: CH3NH2+HCO3- Rx3: CH3NH2+CO2+H2O Rx4: CH3NH2+HCO3-+H2O

12. Lys-70 Lys-70 Trp-154 Trp-154 Ser-67 Ser-67 QM/MM Calculations of the Transition State for Lys-70 Carboxylation Stereoview of the TS showing a molecule of water catalyzing the addition of carbon dioxide to the side chain of Lys-70

13. Carboxylation in the QM/MM model of the active site of OXA-10

14. QM/MM Calculations of the Reactants, TS and Products for Lys-70 Carboxylation QM/MM values in normal text X-ray values in italics

15. OXA-10 β-lactamase - Discussion • A water molecule in the active site can catalyze carboxylation of Lys70 with CO2 • X-ray structure is most likely the deprotonated carboxylation product • Carboxylation is accompanied by deprotonation • Re-protonation of carbamate nitrogen results in barrierless loss of CO2, accounting for biphasic kinetics of enzyme

16. Oxidative Damage to DNA Transformations of 8-hydroxy guanine radicalB3LYP/6-31+G(d) gas phase optimizationIEF-PCMB3LYP/aug-cc-pVTZsolution phase energies B. M. Munk, C. J. Burrows, H. B. Schlegel, Chem. Res. Toxicol. (accepted)

17. Transformations of 8-hydroxy guanine radicalPath 1: reduction followed by tautomerization and ring opening

18. Transformations of 8-hydroxy guanine radicalPath 2: tautomerization followed by ring opening and reduction

19. Transformations of 8-hydroxy guanine radicalPath 3: ring opening followed by reduction and tautomerization

20. Transformations of 8-hydroxy guanine radicalPath 4: ring opening followed by tautomerization and reduction

21. Transformations of 8-hydroxy guanine radical(a) Pathways 2 and 4 are preferred(b) Barriers for ring opening and tautomerization are lower for the radical than for the closed shell molecule

22. AMBER Force Field Parameters for the Naturally Occurring Modified Nucleosides in RNAR. Aduri, B. T. Psciuk, P. Saro, H. B. Schlegel, J. SantaLucia Jr. J. Chem. Theor. Comp. (submitted) N2-methylguanosine 1-methylpseudouridine 5-carboxymethylamino methyluridine 4-demethylwyosine galactosyl-queuosine queuosine

23. Protocol for Determining Atom-centered Partial Charges

24. Modular approach to fitting RESP charges •  The C3’ endo sugar charge was obtained by multi equivalencing the four natural nucleosides •  Two stage RESP was used to fit the ESP of the modified bases and sugars •  Atom types and parameters available in GAFF were sufficient for almost all 103 modifications • The “prepin” and “frcmod” files generated for all 103 modifications • Parameters can be downloaded from http://ozone3.chem.wayne.edu:8080/Modifieds/index.jsp

25. 5MU PSU 2MG 5MC DHU 7MG 2MG M2G 5MC PSU WBG MRC MRG Test Application of the Parameters for Modified RNA Bases tRNAPhe with and without modified bases (1EHZ) H. Shi and P. B. Moore RNA (2000) 6 1091-2000

26. Ab Initio Molecular Dynamics (AIMD) • AIMD – electronic structure calculations combined with classical trajectory calculations • Every time the forces on the atoms in a molecule are needed, do an electronic structure calculation • Born – Oppenheimer (BO) method: converge the wavefunction at each step in the trajectory • Extended Lagrangian methods: propagate the wavefunction along with the geometry • Car-Parrinello – plane-wave basis, propagate MO’s • ADMP – atom centered basis, propagate density matrix

27. Ab Initio Classical Trajectory on theBorn-Oppenheimer Surface Using Hessians Calculate the energy, gradient and Hessian Solve the classical equations of motion on a local 5th order polynomial surface Millam, J. M.; Bakken, V.; Chen, W.; Hase, W. L.; Schlegel, H. B.; J. Chem. Phys. 1999, 111, 3800-5.

28. A Reaction with Branching after the Transition State • Previous work with S. Shaik (JACS1997, 119, 9237 and JACS2001, 123, 130): • Common TS for inner sphere ET and SUB(C) reactions. • Long C-C bond in TS (ca. > 2.45 Å ) favors ET; shorter favors SUB( C ). • A less electronegative halide switches the mechanism from SUB( C ) to ET. • Poorer electron donors of radical anions favor SUB( C ). • More bulkier alkyl halide or more strained TS favor ET.

29. SUB(C) and ET Reaction Paths for CH2O.- + CH3Cl ET CH2O + CH3 + Cl- Sub(C) OCH2CH3 + Cl- D(C-Cl) (bohr) TS D(C-C) (bohr)

30. Energetics at the UHF/6-31G(d) level of theory

32. Temperature dependence of the branching ratio Li, J.; Li, X.; Shaik, S.; Schlegel, H. B. J. Phys. Chem. A 2004, 108, 8526-8532.

33. Energetics at the G3 level of theory

34. Improved Potential Energy Surfaces using Bond Additivity Corrections (BAC) • The most important correction needed for this reaction are C-C and C-Cl bond stretching potentials. • BAC (bond additivity correction) • add simple corrections to get better energetics for the reaction • E = E′+ ∆E • ∆E = AC-CExp{-αC-C RC-C} + AC-ClExp{- αC-Cl RC-Cl} • add the corresponding corrections to gradient and hessian • G = G′+ ∂(∆E)/∂x • H = H′+ ∂2(∆E)/∂x2 • A and α are parameters obtained by fitting to G3 energies

35. BAC-UHF Dynamics Results Table 2. Branching ratios at different levels of theory. Li, J.; Shaik, S.; Schlegel, H. B.; J. Phys. Chem. A2006, 110, 2801-2806..

36. Electronic Response of Molecules Short, Intense Laser Pulses • For intensities of 1014 W/cm2, the electric field of the laser pulse is comparable to Coulombic attraction felt by the valence electrons – strong field chemistry • Need to simulate the response of the electrons to short, intense pulses • Time dependent Schrodinger equations in terms of ground and excited states  =  Ci(t) i i ħ dCi(t)/dt =  Hij(t) Ci(t) • Requires the energies of the field free states and the transition dipoles between them • Need to limit the expansion to a subset of the excitations – TD-CIS, TD-CISD • Time dependent Hartree-Fock equations in terms of the density matrix i ħ dP(t)/dt = [F(t), P(t)] • For constant F, can use a unitary transformation to integrate analytically P(ti+1) = V P(ti)  V† V = exp{ i t F } • Fock matrix is time dependent because of the applied field and because of the time dependence of the density (requires small integration step size – 0.05 au)

37. H2 in an intense laser fieldTD-HF/6-311++G(d,p)Emax = 0.10 au (3.5  1014 W/cm2) = 0.06 au (760 nm) Test Case

38. Laser pulse H2 in an intense laser fieldTD-HF/6-311++G(d,p)Emax = 0.12 au (5.0  1014 W/cm2) = 0.06 au (760 nm) Test Case (a) Instantaneous dipole response (b) (c) Fourier transform of the residual dipole response

39. Hydrogen Molecule aug-pVTZ basis plus 3 sets of diffuse sp shells Emax = 0.07 au (1.7  1014 W/cm2),  = 0.06 au (760 nm) (b) (a) (c) TD-CIS TD-CISD TD-HF (b) (d) (c) (e) (f)

40. Butadiene in an intense laser field(8.75 x 1013 W/cm2 760 nm) HF/6-31G(d,p) Dt = 0.0012 fs

41. The Charge Response of Neutral Butadiene

42. Butadiene in an intense laser fieldTD-CIS/6-31G(d,p), 160 singly excited states = 0.06 au (760 nm) Fourier transform of the residual dipole Excited state weights in the final wavefunction

43. 2  1014 W·cm-2 • 6  1013 W·cm-2 • 5.4  1013 W·cm-2 • 2.7  1013 W·cm-2 • 2.4  1013 W·cm-2 • 5.0  1012 W·cm-2 • 4.5  1012 W·cm-2 0 10 20 30 40 0 10 20 30 40 Polyacenes in Intense Laser Pulse (Levis et al. Phys. Rev. A 69, 013401 (2004)) • 1  1014 W·cm-2 Ion Signal, normalized Time-of-flight, ms

44. TDHF Simulations for Polyacenes • Polyacenes ionize and fragment at much lower intensities than polyenes • Polyacene experimental data shows the formation of molecular +1 cations prior to fragmentation with 60 fs FWHM pulses • Time-dependent Hartree-Fock simulations with 6-31G(d,p) basis, Dt = 0.0012 fs, ω=1.55 eV and 5 fs FWHM pulse • Intensities chosen to be ca 75% of the experimental single ionization intensities • Intensities of 8.75 x 1013, 3.08 x 1013, 2.1 x 1013 and 4.5 x1012 for benzene, naphthalene, anthracene and tetracene • Nonadiabatic multi-electron excitation model was used to check that these intensities are non-ionizing

45. Tetracene: Dipole Response I = 3.38 x 1012 W/cm2 ω = 1.55eV, 760 nm

46. ω = 1.55eV, 760 nm Naphthalene+1:Dependence on the Field Strength

47. 1.1 eV 1.1 eV E = 0.029 au E = 0.034 au E = 0.0155 au 1.1 eV 7.1 eV 4.5 eV Transition Amplitude 4.5 eV 7.1 eV 7.1 eV 8.95 eV 3.1 eV (2x1.55 eV) Energy (eV) Energy (eV) Energy (eV) ω = 1.55eV, 760 nm Naphthalene+1:Dependence on the Field Strength

48. Anthracene+1:Dependence on the Field Frequency Emax = 0.0183 au

49. Anthracene+1:Dependence on the Field Frequency Emax = 0.0183 au 3.63 eV ω = 1.00 eV ω = 3.00 eV ω = 2.00 eV 25 1.95 eV 1.95 eV 50 40 2.79 eV 20 40 4.61 eV 3.63 eV 3.63 eV 30 Transition Amplitude 15 5.58 eV 30 4.95 eV 7.97 eV 20 10 6.32 eV 6.32 eV 20 7.79 eV 10 5 10 10.23 eV 7.97 eV 9.57 eV Energy Energy Energy

50. Polyacenes: Summary • Non-adiabatic behavior increases with length • Non-adiabatic behavior is greater for monocation • Increasing the field strength increases the non-resonant excitation of the states with the largest transition dipoles • Increasing the field frequency increases the non-resonant excitation of higher states Smith, S. M.; Li, X.; Alexei N. Markevitch, A. N.; Romanov, D. A.; Robert J. Levis, R. J.; Schlegel, H. B.; Numerical Simulation of Nonadiabatic Electron Excitation in the Strong Field Regime: 3. Polyacene Neutrals and Cations. (JPCA submitted)