Exploring Potential Energy Surfaces by Ab Initio Molecular Dynamics

1
Exploring Potential Energy Surfaces by Ab Initio
Molecular Dynamics
Molecular Quantum MechanicsAnalytical Gradients
and Beyond A Conference in Honor of Peter
Pulay May 29 June 3, 2007, Margitsziget,
Budapest
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
Early encounters with Peter Pulay and his work

• During my graduate studies, 1972-75
• HF gradients for s,p basis sets
• Pulay was the external examiner on PhD thesis
• While a postdoc with John Pople, 1977-78
• HF second derivatives and post-SCF gradients with
  Krishnan Raghavachari

3
Reactions and Dynamics

• Modeling reactive potential energy surfaces with
  empirical valence bond methods using distributed
  Gaussians
• energy derivatives, redundant internal
  coordinates, DIIS
• Dynamics of the electron density of molecules
  interacting with intense laser pulses
• TD-HF and TD-CIS simulations of the response
  prior to ionization

4
Empirical Valence Bond Models for Reactive
Potential Energy Surfaces Using Distributed
Gaussians

• H. B. Schlegel, J. L. Sonnenberg,
• J. Chem. Theor. Comp. 2006, 2, 905-911.
• J. L. Sonnenberg, H. B. Schlegel,
• Mol. Phys. (submitted).
• Supported by a grant from ONR to
• Voth, Miller, Case, Cheatham and Schlegel

5
Empirical Valence Bond Models for Potential
Energy Surfaces for Reactions

• For condensed phase and enyzmatic systems where
  extensive dynamics is required, QM/MM
  calculations may still be too costly.
• EVB provides a simple, systematic way to
  construct an empirical PES for reactions,
  calibrated to QM calculations
• PES represented by 2 (or more) valence bond
  configurations and empirical interaction matrix
  elements
• Initially employed by Warshel
• Improved by Chang and Miller
• MCMM method by Truhlar

6
Empirical Valence Bond Models for Potential
Energy Surfaces for Reactions

• V11 and V22 can be treated adequately by
  molecular mechanics
• Need to find a suitable functional form for V12
• Fit V12 to ab initio calculation around the
  transition state

7
Empirical Valence Bond Models for Potential
Energy Surfaces for Reactions

• Chang Miller approach
• Fit V12 to energy, gradient and Hessian at or
  near transition state

8
Empirical Valence Bond Models for Potential
Energy Surfaces for Reactions

• Our approach
• Fit V12 to energy, gradient and Hessian at one or
  more points around the transition state

9
Empirical Valence Bond Models for Potential
Energy Surfaces for Reactions
HCN?HNC isomerization PES effect of
coordinate system Cartesian
Harmonic internal Anharmonic internal
Actual PES
10
Empirical Valence Bond Model using Distributed
Gaussians
11
Potential Energy Surface for 2-Pyridone
2-Hydroxypyridine Tautomerization
MP2/6-311G(d,p) optimization, 159 kJ/mol
barrier EVB surface with distributed Gaussians
fitted energies, gradients and Hessians at one to
nine points along the reaction path Redundant
internal coordinates used 56 (all), 32 (only
bonds and angles), 25 (changes
greater than 0.01), 10 (changes greater than 0.1
bohr or rad)
12
V122 along the reaction path for pyridone
tautomerization
V122 (au)
V122 difference (au)
Long dashes 3 Gaussian fit (R, P, TS) Short
dashes 7 Gaussian fit (R, P, TS, ½, ¼ ) Red
ab initio V122 along the path
13
Energy along the reaction path for pyridone
tautomerization
Total energy (au)
Energy difference (au)
1 Gaussian fit at TS max error 7.3 kJ/mol Long
dashes 3 Gaussian fit (R, P, TS) max error
4.1 kJ/mol Short dashes 7 Gaussian fit (R, P,
TS, ½, ¼ ) max error 0.4 kJ/mol Red ab
initio energy along the path
14
Gradient along the reaction path for pyridone
tautomerization
Long dashes 3 Gaussian fit (R, P, TS) Short
dashes 7 Gaussian fit (R, P, TS, ½, ¼ ) Red
ab initio gradient norm along the path
15
Potential energy surface for pyridone
tautomerization
16
Potential Energy Surface for Pyridone Water
Tautomerization
B3LYP/6-31G(d,p) opt., 50 kJ/mol barrier 1
Gaussian fit V122 using all coord., max
error 4.0 kJ/mol V122 using only bonds,
max error 2.9 kJ/mol 3 Gaussian fit (R, TS,
P) V122 using all coord., max error 1.8
kJ/mol V122 using only bonds, max error
3.0 kJ/mol
17
Potential Energy Surface for theClaisen
Rearrangment
MP2/6-311G(d,p) optimization, 99 kJ/mol
barrier EVB surface with 7 Gaussians has a
maximum error of 3.3 kJ/mol along the reaction
path
18
Electronic Response of Molecules to Short,
Intense Laser Pulses

• Phys. Rev. A. 2003, 68, 011402(R),
• Phys. Rev. A. 2004, 69, 013401,
• Phys. Chem. Chem. Phys. 2005, 7, 233-239,
• J. Phys. Chem. A 2005 109 5176-5185,
• J. Phys. Chem. A 2005, 109, 10527-10534,
• J. Phys. Chem. A 2007 (accepted),
• J. Chem. Phys. 2007 (accepted).
• Supported by a grant from NSF

19
TD-CI and TD-HF simulation of molecules in 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 h 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 h dP(t)/dt F(t), P(t)
• For constant F, can use a unitary transformation
  to integrate analytically
• P(ti1) 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)

20
H2 in an intense laser fieldTD-HF/6-311G(d,p)
Emax 0.10 au (3.5 ? 1014 W/cm2) ?
0.06 au (760 nm)
Test Case
21
H2 in an intense laser fieldTD-HF/6-311G(d,p)
Emax 0.12 au (5.0 ? 1014 W/cm2) ?
0.06 au (760 nm)
Laser pulse
Test Case
(a)
Instantaneous dipole response
(b)
(c)
Fourier transform of the residual dipole response
22
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)
23
Butadiene in an intense laser field(8.75 x 1013
W/cm2 760 nm)
HF/6-31G(d,p) Dt 0.0012 fs
24
The Charge Response of Neutral Butadiene
25
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
26
Hexatriene in an intense laser fieldTD-CIS/6-31G(
d,p), 200 singly excited states? 0.06 au (760
nm)
Fourier transform of the residual dipole
Excited state weights in the final wavefunction
27
Polyacenes in Intense Laser Pulse (Levis, R. J.
et al. Phys. Rev. A 2004, 69, 013401)
1 ? 1014 Wcm-2
Ion Signal, normalized
Time-of-flight, ms
28
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

29
Anthracene Dipole Response
I 1.58 x 1013 W/cm2
? 1.55eV, 760 nm
30
Tetracene Dipole Response
I 3.38 x 1012 W/cm2
? 1.55eV, 760 nm
31
Naphthalene in an intense laser
fieldTD-CIS/6-31G(d,p), 200 singly excited
states? 0.06 au (760 nm)
Fourier transform of the residual dipole
Excited state weights in the final wavefunction
32
Anthracene1Dependence on the Field Frequency
Emax 0.0183 au
33
Anthracene1Dependence 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
34
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.
  J. Phys. Chem. A (submitted)

35
Recent Group Members
36
Current Group Members
37
Acknowledgements

• Collaborators
• Dr. T. Vreven, Gaussian Inc.
• Dr. M. J. Frisch, Gaussian Inc.
• Prof. John SantaLucia, Jr., WSU
• Raviprasad Aduri (SantaLucia group)
• Prof. G. Voth, U. of Utah
• Prof. David Case, Scripps
• Prof. Bill Miller, UC Berkeley
• Prof. Thom Cheatham, U. of Utah
• Prof. S.O. Mobashery, Notre Dame U.
• Prof. R.J. Levis, Temple U.
• Prof. C.H. Winter, WSU
• Prof. C. Verani, WSU
• Prof. E. M. Goldfield, WSU
• Prof. D. B. Rorabacher, WSU
• Prof. J. F. Endicott, WSU
• Prof. J. W. Montgomery, U. of Michigan
• Prof. Sason Shaik, Hebrew University
• Prof. P.G. Wang, Ohio State U.

Current Research Group Dr. Jason Sonnenberg Dr.
Peng Tao Barbara Munk Michael Cato Jia
Zhou Jason Sonk Brian Psciuk Recent Group
Members Prof. Xiaosong Li, U of Washington Prof.
Smriti Anand, Christopher-Newport U. Dr. Hrant
Hratchian, Indiana U. (Raghavachari grp) Dr. Jie
Li, U. California, Davis (Duan group) Dr. Stan
Smith, Temple U. (Levis group) Dr. John Knox
(Novartis) Funding and Resources National
Science Foundation Office of Naval
Research NIH Gaussian, Inc. Wayne State U.
38
Molecular geometriesand orientation of the field
39
Effect of Charge and Geometry on the Dipole
Moment Response Butadiene
I 8.75 x 1013 W/cm2
? 1.55eV, 760 nm
40
Butadiene1 Fourier Analysis of Residual
Oscillations
Main Transition (TDHF Coefficient) Energy (eV) Transition Dipole (au) Oscil. Stren.
Neutral Geometry
HOMO ? SOMO (1.00) 2.53 1.70 0.12
SOMO ? LUMO (0.92) 4.87 1.75 0.38
Ion Geometry
HOMO ? SOMO (0.95) 4.03 1.94 0.39
HOMO ? LUMO (0.83) 5.69 0.23 0.01
4.10 eV
2.57 eV
4.90 eV
Transition Amplitude
Neutral Geometry
Ion Geometry
41
Polyene Cations Summary

• The monocations have lower energy excited states
  and show greater non-adiabatic behavior than the
  dications
• Relaxing the geometry increases the energy of
  the lowest excited states and decreases the
  non-adiabatic behavior
• Fourier transform of the residual oscillations
  in the dipole moment shows that the non-adiabatic
  excitation involves the lowest excited 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 2. Linear Polyene Cations. J.
  Phys. Chem. A 2005, 109, 10527-10534.

42
Ionization Probability using NME
Molecule Excited State Energy (?) Transition Dipole Moment (au) Ionization Probability
Benzene
Neutral 8.00 1.8700 0.0022
1 Neutral Geometry 5.59 1.0294 0.0052
1 Ion Geometry 5.60 1.1165 0.0054
2 Neutral Geometry 7.53 1.3739 0.00034
2 Ion Geometry 7.41 1.2511 0.00023
Naphthalene
Neutral 6.94 3.0104 0.0001
1 Neutral Geometry 6.99 1.7097 0.011
1 Ion Geometry 7.26 2.3353 0.00018
2 Neutral Geometry 6.63 2.3504 0.00025
2 Ion Geometry 6.35 2.3733 0.00026
Anthracene
Neutral 6.21 3.9917 0.00015
1 Neutral Geometry 6.37 2.7456 0.00047
2 Neutral Geometry 5.94 3.1693 0.00056
Tetracene
Neutral 5.70 4.8840 0.00006
1 Neutral Geometry 6.28 3.3935 0.0011
2 Neutral Geometry 5.44 3.8964 0.00095
43
Pulse Shaping and Sequencing
Test case HF/6-31G(d,p) calculations of
excitation energies, transition dipoles and time
dependent response of H2
Transition dipole for 1?g? to 1?u? is 1.00 au,
and for 1?u? to 2?g? is 1.57 au, but no
transition dipole between 1?g? to 2?g?
Ladder-type STIRAP experiment (stimulated Raman
adiabatic passage)
A single-photon excitation from 1?g? ? 1?u? and
a second single-photon excitation from 1?u? ?
2?g? should populate the 2?g? state with
little or no population in the 1?u? state if the
1?u? ? 2?g? is pumped first
44
Electric Field Profiles for H2 with
Counter-Intuitive Pulse Timing
Pump Pulse
Stokes Pulse
FWHMP 85 fs
FWHM 120 fs
150 fs delay between pulses
45
STIRAP Results for H2
46
Effects of Detuning in Ladder Type STIRAP for H2
47
Effects of Detuning for H2
48
Ehrenfest Dynamics

• Time-dependent HF or DFT propagation of the
  electron density
• Classical propagation of the nuclear degrees of
  freedom
• Novel integration method using three different
  time scales
• Li, X. Tully, J. C. Schlegel, H. B. Frisch,
  M. J. Ab Initio Ehrenfest Dynamics. J. Chem.
  Phys. 2005, 123, 084106

49
Potential energy curves for H2CNH2 torsion
50
Torsional dynamicsfor H2CNH2