Srinivasan S. Iyengar Department of Chemistry and Department of Physics, Indiana University

1. Atom-centered Density Matrix Propagation (ADMP): Theory and Applications Srinivasan S. Iyengar Department of Chemistry and Department of Physics, Indiana University

2. Outline • Brief discussion of ab initio molecular dynamics • Atom-centered Density Matrix Propagation (ADMP) • Nut-n-bolts issues • Some Results: • Novel findings for protonated water clusters • QM/MM generalizations: ion channels • Gas phase reaction dynamics

3. Molecular dynamics on a single potential surface • Parameterized force fields (e.g. AMBER, CHARMM) • Energy, forces: parameters obtained from experiment • Molecular motion: Newton’s laws • Works for large systems • But hard to parameterize bond-breaking/formation (chemical reactions) • Issues with polarization/charge transfer/dynamical effects • Born-Oppenheimer (BO) Dynamics • Solve electronic Schrödinger eqn (DFT/HF/post-HF) for each nuclear structure • Nuclei propagated using gradients of energy (forces) • Works for bond-breaking but computationally expensive • Large reactive, polarizable systems: Something like BO, but preferably less expensive.

4. References… CP: R. Car, M. Parrinello, Phys. Rev. Lett. 55 (22), 2471 (1985). ADMP: Schlegel, et al. JCP, 114, 9758 (2001). Iyengar, et al. JCP, 115,10291 (2001). Iyengar et al. Israel J. Chem. 7, 191, (2002). Schlegel et al. JCP 114, 8694 (2002). Iyengar and Frisch JCP 121, 5061 (2004). Extended Lagrangian dynamics • Circumvent Computational Bottleneck of BO • Avoid repeated SCF: electronic structure, not converged,but propagated • “Simultaneous” propagation of electronic structure and nuclei: adjustment of time-scales • Car-Parrinello (CP) method • Orbitals expanded in plane waves • Occupied orbital coefficients propagated • O(N3) computational scaling (traditionally) • O(N) with more recent Wannier representations (?) • Atom-centered Density Matrix Propagation (ADMP) • Atom-centered Gaussian basis functions • Electronic Density Matrix propagated • Asymptotic linear-scaling with system size • Allows the use of accurate hybrid density functionals • suitable for clusters

5. Lagrangian Constraint for N-representability of P: Idempotency and Particle number Nuclear KE “Fictitious” KE of P Energy functional Atom-centered Density Matrix Propagation (ADMP) • Construct a classical phase-space {{R,V,M},{P,W,m}} • TheLagrangian (= Kinetic minus Potential energy) • P: representedusing atom-centered gaussian basis sets

6. acceleration of density matrix,P “Fictitious” mass ofP Force onP Euler-Lagrange equations of motion for ADMP • Equations of motion for density matrix and nuclei • Classical dynamics in {{R,V,M},{P,W,m}}phase space • Next few slide: Forces, propagation equations, formal error analysis

7. Hellman-Feynman contributions Pulay’s moving basis terms S=UTU, Cholesky or Löwdin Contributions due to [F,P]  0. Part of non-Hellman-Feynman Nuclear Forces: What Really makes it work

8. Density Matrix Forces: • Use McWeeny Purified DM (3P2-2P3) in energy expression to obtain

9. Direction of Increasing Frequency • Consequence of m: P changes slower with time: characteristic frequency adjusted • Consequence of m: P changes slower with time: characteristic frequency adjusted • Consequence of m: P changes slower with time: characteristic frequency adjusted • But Careful - too largem:non-physical • But Careful - too largem:non-physical • Appropriate m: approximate BO dynamics meffects an adjustment of time-scales: • Bounds for m: From a Hamiltonian formalism • m: also related to deviations from the BO surface

10. Commutator of the electronic Hamiltonian and density matrix: boundedbymagnitude ofm Iyengar et al. Israel J. Chem. 7, 191, (2002). Reference… “Physical” interpretation of m: Bounds • Magnitude of m: represents deviation from BO surface • macts as an “adiabatic control parameter”

11. The Lagrangian • The Conjugate Hamiltonian (Legendre Transform) Iyengar et al. JCP. 115,10291 (2001). Reference… Bounds on the magnitude of m : Controlling m:Deviations from BO surface and adiabaticity

12. Born-Oppenheimer dynamics: Converged electronic states. Approx. 8-12 SCF cycles / nuclear config. dE/dR not same in both methods ADMP: Electronic statepropagatedclassically : no convergence reqd. 1 SCF cycle : for Fock matrix -> dE/dP Current: 3-4 times faster. Comparison with BO dynamics References… Iyengar et al. Israel J. Chem. 7, 191, (2002). Schlegel et al. JCP 114, 8694 (2002). Iyengar and Frisch JCP 121, 5061 (2004).

13. Propagation of W Propagation of P: time-reversible propagation • Velocity Verlet propagation of P • Classical dynamics in {{R,V},{P,W}}phase space • Li and Li+1 obtained iteratively: • Conditions: Pi+12 = Pi+1 and WiPi + PiWi = Wi (next two slides)

14. Idempotency (N-Representibility of DM): • Given Pi2 = Pi, need Li to find idempotent Pi+1 • Solve iteratively: Pi+12 = Pi+1 • Given Pi, Pi+1, Wi, Wi+1/2, need Li+1 to find Wi+1 • Solve iteratively: Wi+1Pi+1 + Pi+1 Wi+1=Wi+1

15. Idempotency: To obtain Pi+1 • Given Pi2 = Pi, need to find indempotent Pi+1 • Guess: • Or guess: • Iterate Pi+1 to satisfy Pi+12 = Pi+1 • Rational for choice PiTPi + QiTQi above:

16. Idempotency: To obtain Wi+1 • Given WiPi + PiWi = Wi, find appropriate Wi+1 • Guess: • Iterate Wi+1 to satisfy Wi+1Pi+1 + Pi+1Wi+1 = Wi+1

17. How it all works … • Initial config.: R(0). Converged SCF: P(0) • Initial velocities V(0) and W(0) : flexible • P(Dt), W(Dt) : from analytical gradients and idempotency • Similarly for R(Dt) • And the loop continues…

18. References… ADMP treatment of protonated water clusters: Iyengar, et al. JCP, 123, 084309 (2005). Iyengar et al. Int. J. Mass Spec.241, 197 (2005). Iyengar JCP 123, 084310, (2005). Protonated Water Clusters • Important systems for: • Ion transport in biological and condensed systems • Enzyme kinetics • Acidic water clusters: Atmospheric interest • Electrochemistry • Experimental work: • Mass Spec.: Castleman • IR: M. A. Johnson, Mike Duncan, M. Okumura • Sum Frequency Generation (SFG) : Y. R. Shen, M. J. Schultz and coworkers • Lots of theory too: Jordan, McCoy, Bowman, Klein, Singer (not exhaustive by any means..) • Variety of medium-sized protonated clusters using ADMP

19. Protonated Water Clusters: Hopping via the Grotthuss mechanism True for 20, 30, 40, 50 and larger clusters…

20. (H2O)20H3O+: Magic number cluster • Hydronium goes to surface: 150K, 200K and 300K: B3LYP/6-31+G** and BPBE/6-31+G** • Castleman’s experimental results: • 10 “dangling” hydrogens in cluster • Found by absorption of trimethylamine (TMA) • 10 “dangling” hydrogens: consistent with our ADMP simulations • But: hydronium on the surface

21. (H2O)20H3O+: A recent spectroscopic quandry Theory Experiment J.-W. Shin, N. I. Hammer, E. G. Diken et al., Science 304, 1137 2004.

22. Spectroscopy: A recent quandry Water Clusters: Important in Atmospheric Chemistry Bottom-right spectrum From ADMP agrees well with expt: dynamical effects in IR spectroscopy Explains the experiments of M. A. Johnson

23. Spectroscopy: A recent quandry ADMP Spectrum!! Iyengar et al. JCP, 123 , 084309 (2005)

24. (H2O)20H3O+: Magic number cluster • Hydronium goes to surface: 150K, 200K and 300K: B3LYP/6-31+G** and BPBE/6-31+G** • Castleman’s experimental results: • 10 “dangling” hydrogens in cluster • Found by absorption of trimethylamine (TMA) • 10 “dangling” hydrogens: consistent with our ADMP simulations • But: hydronium on the surface

25. Larger Clusters and water/vacuum interfaces: Similar results

26. Predicting New Chemistry: Theoretically A Quanlitative explanation to the remarkable Sum Frequency Generation (SFG) of Y. R. Shen, M. J. Schultz and coworkers

27. Protonated Water Cluster: Conceptual Reasons for “hopping” to surface Hydrophobic and hydrophillic regions: Directional hydrophobicity (it is amphiphilic) H3O+ has reduced density around Reduction of entropy of surrounding waters Is Hydronium hydrophobic ? H2O coordination 4 H3O+coordination =3

28. References… P. B. Miranda and Y. R. Shen, J. Phys. Chem. B, 103, 3292-3307 (1999). M. J. Schultz, C. Schnitzer, D. Simonelli and S. Baldelli, Int. Rev. Phys. Chem. 19, 123-153 (2000) Experimental results suggest this as well • Y. R. Shen: Sum Frequency Generation (SFG) • IR for water/vapor interface shows dangling O-H bonds • intensity substantially diminishes as acid conc. is increased • Consistent with our results • Hydronium on surface: lone pair outwards, instead of dangling O-H • acid concentration is higher on the surface • Schultz and coworkers: acidic moieties alter the structure of water/vapor interfaces

29. QM/MM treatment: ONIOM ADMP Unified treatment of the full system within ADMP (I) (This talk will not overview the ONIOM scheme, but the interested reader should look at the reference below) N. Rega, S. S. Iyengar, G. A. Voth, H. B. Schlegel, T. Vreven and M. J. Frisch, J. Phys. Chem. B 108 4210 (2004).

30. Side-chain contribute to hop “Eigen” like configuration possible using protein backbone B3LYP and BLYP: qualitatively different results

31. HCHO photodissociation • Photolysis at 29500 cm-1 : To S1 state • Returns to ground state vibrationally hot • Product: rotationally cold, vibrationally excited H2 • And CO broad rotational distr: <J> = 42. Very little vib. Excitation • H2CO  H2 + CO: BO and ADMP at HF/3-21G, HF/6-31G**

32. Glyoxal 3-body Synchronous photo-fragmentation

33. Conclusions • ADMP: powerful approach to ab initio molecular dynamics • Linear scaling with system size • Hybrid (more accurate) density functionals • Smaller values for fictitious mass allow • treatment of systems with hydrogens is easy (no deuteriums required) • greater adiabatic control (closer to BO surface) • Examples bear out the accuracy of the method

34. Acknowledgment • The work has enormously benefited from my former advisors and collaborators: • Greg Voth • Berny Schlegel • Gus Scuseria • Mike Frisch • At IU, people contributing to this work are: • Jacek Jakowski (post-doc) • Isaiah Sumner (grad student) • Xiaohu Li (grad student) • Virginia E. Teige (Freshman)