Quantum wavepacket ab initio molecular dynamics A computational

Quantum wavepacket ab initio molecular dynamics: A computational approach for quantum dynamics in large systems Srinivasan S. Iyengar Department of Chemistry and Department of Physics, Indiana University Funding: Group members contributing to this work: Jacek Jakowski (post-doc), Isaiah Sumner (Ph. D student), Xiaohu Li (Ph. D student), Virginia Teige (BS, first year student)

Iyengar Group, Indiana University Predictive computations: a few (grand) challenges Lipoxygenase: enzyme Ø Bio enzyme: Lipoxygenase: Fatty acid oxidation • • • Rate determining step: hydrogen abstraction from fatty acid KIE (k. H/k. D)=81 – Deuterium only twice as heavy as Hydrogen – generally expect k. H/k. D = 3 -8 ! weak Temp. dependence of rate Ø Nuclear quantum effects are critical Ø Conduction across molecular wires • Is the wire moving? Ion (proton) channels Reactive over multiple sites Ø Polarization due to electronic factor Ø Polymer-electrolyte fuel cells Ø Dynamics & temperature effects Ø

Iyengar Group, Indiana University Chemical Dynamics of electron-nuclear systems Ø Our efforts: approach for simultaneous dynamics of electrons and nuclei in large systems: • • • accurate quantum dynamical treatment of a few nuclei, bulk of nuclei: treated classically to allow study of large (enzymes, for example) systems. Electronic structure simultaneously described: evolves with nuclei Spectroscopic study of small ionic clusters: including nuclear quantum effects Ø Proton tunneling in biological enzymes: ongoing effort Ø

Iyengar Group, Indiana University Hydrogen tunneling in Soybean Lipoxygenase 1: Introduce Quantum Wavepacket Ab Initio Molecular Dynamics Catalyzes oxidation of unsaturated fat Expt Observations “Quantum” nuclei Ø Rate determining step: hydrogen abstraction from fatty acid Ø Weak temperature dependence of k Ø k. H/k. D = 81 • • • Deuterium only twice as heavy as Hydrogen, generally expect k. H/k. D = 3 -8. Remarkable deviation The electrons and the “other” classical nuclei

Iyengar Group, Indiana University Quantum Wavepacket Ab Initio Molecular Dynamics The “Quantum” nuclei [Distributed Approximating Functional (DAF) approximation to free propagator] Ab Initio Molecular Dynamics (AIMD) using: The electrons and the “other” classical nuclei Atom-centered Density Matrix Propagation (ADMP) OR References… Born-Oppenheimer Molecular Dynamics (BOMD) S. S. Iyengar and J. Jakowski, J. Chem. Phys. 122 , 114105 (2005). Iyengar, TCA, In Press. J. Jakowski, I. Sumner and S. S. Iyengar, JCTC, In Press (Preprints: author’s website. )

Iyengar Group, Indiana University 1. DAF quantum dynamical propagation Quantum Dynamics subsystem: • Quantum Evolution: Linear combination of Hermite functions: The “Distributed Approximating Functional” is a banded, Toeplitz matrix Time-evolution: vibrationally non-adiabatic!! (Dynamics is not stuck to the ground vibrational state of the quantum particle. ) Linear computational scaling with grid basis

2. Iyengar Group, Indiana University Quantum dynamically averaged ab Initio Molecular Dynamics • Averaged BOMD: Kohn Sham DFT for electrons, classical nucl. Propagation • Approximate TISE for electrons • Computationally expensive. • Quantum averaged ADMP: • Classical dynamics of {RC, P}, through an adjustment of time-scales acceleration of density matrix, P “Fictitious” mass tensor of P Force on P • V(RC, P, RQM; t) : the potential that quantum wavepacket experiences Ref. . Schlegel et al. JCP, 114, 9758 (2001). Iyengar, et. al. JCP, 115, 10291 (2001).

Iyengar Group, Indiana University Quantum Wavepacket Ab Initio Molecular Dynamics: The pieces of the puzzle The “Quantum” nuclei [Distributed Approximating Functional (DAF) approximation to free propagator] Simultaneous dynamics Ab Initio Molecular Dynamics (AIMD) using: The electrons and the “other” classical nuclei Atom-centered Density Matrix Propagation (ADMP) OR Born-Oppenheimer Molecular Dynamics (BOMD) S. S. Iyengar and J. Jakowski, J. Chem. Phys. 122 , 114105 (2005) J. Jakowski, I. Sumner, S. S. Iyengar, J. Chem. Theory and Comp. In Press

Iyengar Group, Indiana University So, How does it all work? • A simple illustrative example: dynamics of Cl. HCl • • • Chloride ions: AIMD Shared proton: DAF wavepacket propagation Electrons: B 3 LYP/6 -311+G** • As Cl- ions move, the potential experienced by the “quantum” proton changes dramatically. • The proton wavepacket splits and simply goes crazy!

Iyengar Group, Indiana University Spectroscopic Properties Ø The time-correlation function formalism plays a central role in non-equilibrium statistical mechanics. When A and B are equivalent expressions, eq. (18) is an autocorrelation function. Ø The Fourier Transform of the velocity autocorrelation function represents the vibrational density of states. Ø

Vibrational spectra including quantum dynamical effects Iyengar Group, Indiana University Cl. HCl- system: large quantum effects from the proton Ø Simple classical treatment of the proton: Ø • • • Ø Geometry optimization and frequency calculations: Large errors Dimensionality of the proton is also important: – 1 D, 2 D and 3 D treatment of the quntum proton provides different results. Mc. Coy, Gerber, Ratner, Kawaguchi, Neumark … In our case: Use the wavepacket flux and classical nuclear velocities to obtain the vibrational spectra directly: • Includes quantum dynamical effects, temperature effects (through motion of classical nuclei) and electronic effects (DFT). In good agreement with Kawaguchi’s IR spectra References… J. Jakowski, I. Sumner and S. S. Iyengar, JCTC, In Press (Preprints: Iyengar Group website. )

Iyengar Group, Indiana University The Main Bottleneck: The work around: Time-dependent Deterministic Sampling (TDDS) • Consider the phenol amine system Need the quantum mechanical Energy at all these grid points!! • • However, some regions are more important than others? Addressed through TDDS, “on-the-fly”

Iyengar Group, Indiana University The Main Bottleneck: The quantum interaction potential 1. Quantum Dynamics subsystem: 2. AIMD subsystem (ADMP for example) The Interaction Potential: A major computational bottleneck • • The potential for wavepacket propagation is required at every grid point!! • Expensive from an electronic structure perspective And the gradients are also required at these grid points!!

Iyengar Group, Indiana University Time-dependent deterministic sampling 1) Importance of each grid point (RQM) based on: - large wavepacket density - r - potential is low - V - gradient of potential is high - V 2) So, the sampling function is: Ir , IV’ --- adjust importance of each component

TDDS - Haar wavelet decomposition Iyengar Group, Indiana University

Iyengar Group, Indiana University Generalization to multidimensions - Haar wavelet decomposition

Iyengar Group, Indiana University TDDS/Haar: How well does it work? The error, when the potential is evaluated only on a fraction of the points is really negligble!!! 1 m. Eh = 0. 0006 kcal/mol = 2. 7 * 10 -5 e. V Hence, PADDIS reproduces the energy: Computational gain three orders of magnitude!!

Iyengar Group, Indiana University TDDS/Haar: Reproduces vibrational properties? The error in the vibrational spectrum: negligible These spectra include quantum dynamical effects of proton along with electronic effects!

Iyengar Group, Indiana University Hydrogen tunneling in biological enzymes: The case for Soybean Lipoxygenase 1 Lipoxygenase: enzyme Enzyme active site shown Ø Catalyzes the oxidation of unsaturated fat! Ø Rate determining step: hydrogen abstraction Ø Weak temperature dependence of k Ø Hydrogen to deuterium KIE is 81 Ø • • Deuterium is only twice as larger as Hydrogen, Generally expect k. H/k. D = 3 -8.

Iyengar Group, Indiana University Soybean Lipoxygenase 1: Lipoxygenase: enzyme A slow time-scale process for AIMD Ø Improved computational treatment Ø through “forced” ADMP. • The idea is the donor atom is “pulled” slowly along the reaction coordinate Bottomline: Donor acceptor distance is not constant during the hydrogen transfer process. Ø The donor-acceptor motion reduces barrier height Ø

Iyengar Group, Indiana University Soybean Lipoxygenase 1: Proton nuclear “orbitals”: Look for the “p” and “d” type functions!! s-type p-type d-type These states are all within 10 kcal/mol Eigenstates obtained from Arnoldi iterative procedure

Reactant Iyengar Group, Indiana University Transition State Eigenstates obtained using: Instantaneous electronic structure (DFT: B 3 LYP) Ø finite difference approximation to the proton Hamiltonian. Ø Arnoldi iterative diagonalization of the resultant large (million by million) eigenvalue problem. Ø For Deuterium, the excited proton state contributions are about 10% For hydrogen the excited state contribution is about 3% Significant in an Marcus type setting.

Iyengar Group, Indiana University Transition state quantum H D classical

Iyengar Group, Indiana University Conclusions and Outlook Ø Quantum Wavepacket ab initio molecular dynamics: Seems Robust and Powerful • Quantum dynamics: efficient with DAF – Vibrational non-adiabaticity for free • AIMD efficient through ADMP or BOMD – Potential is determined on-the-fly! • Ø Importance sampling extends the power of the approach In Progress: • • • QM/MM generalizations: Enzymes generalizations to higher dimensions and more quantum particles: Condensed phase Extended systems (Quantum Dynamical PBC): Fuel cells

Iyengar Group, Indiana University Additional slides

Iyengar Group, Indiana University Optimization of w‘(RQM) with respect to a, b, g a=1 (IY) b=3 (IYP) g=1 (IChi) RMS error of intrepolation during a dynamics within mikrohartrees

Iyengar Group, Indiana University Computational advantages to DAF quantum propagation scheme • Free Propagator: is a banded, Toeplitz matrix: • Time-evolution: vibrationally non-adiabatic!! (Dynamics is not stuck to the ground vibrational state of the quantum particle. )

Group, Indiana University Quantum Wavepacket Ab Initio. Iyengar Molecular Dynamics: Working Equations Quantum Dynamics subsystem: Trotter Coordinate representation: • • The action of the free propagator on a Gaussian: exactly known Expand the wavepacket as a linear combination of Hermite Functions are derivatives of Gaussians Therefore, the action of free propagator on the Hermite can be obtained in closed form: • Coordinate representation for the free propagator. Known as the Distributed Approximating Functional (DAF) [Hoffman and Kouri, c. a. 1992] Wavepacket propagation on a grid •

Spreading transformation Iyengar Group, Indiana University We want to do potential evaluation for η fraction of grid -Density from ω(x) may be larger than current grid densityexceeding density is spread over low density grid area - for η 1 weighting ω(x) should tend to 1 Grid point for potential evaluation are deteminned by integrating [N*w(x)] Interpolation of potential Version of cubic spline interpolation - based on on potentials and gradients - easy to generalize in multidimensions - general flexible form

Group, Indiana University Another example: Proton transfer. Iyengar in the phenol amine system • • • Shared proton: DAF wavepacket propagation All other atoms: ADMP Electrons: B 3 LYP/6 -31+G** • C-C bond oscilates in phase with wavepacket Wavepacket amplitude near amine Scattering probability: References… S. S. Iyengar and J. Jakowski, J. Chem. Phys. 122 , 114105 (2005).

Iyengar Group, Indiana University Potential Adapted Dynamically Driven Importance Sampling (PADDIS) : basic ideas The following regions of the potential energy surface are important: -Regions with lower values of potential -That’s probably where the WP likes to be -Regions with large gradients of potential -Tunneling may be important here -Regions with large wavepacket density Consequently, the PADDIS function is: The parameters provide flexibility