1 / 37

Comparison of Born-Oppenheimer and Atom-Centered Density Matrix Propagation Methods for Ab Initio Molecular Dynamics

Comparison of Born-Oppenheimer and Atom-Centered Density Matrix Propagation Methods for Ab Initio Molecular Dynamics. H. Bernhard Schlegel Dept. of Chemistry, Wayne State U. Detroit, Michigan 48202, USA. Approaches for Ab Initio Molecular Dynamics. Born-Oppenheimer Dynamics

jaunie
Download Presentation

Comparison of Born-Oppenheimer and Atom-Centered Density Matrix Propagation Methods for Ab Initio Molecular Dynamics

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Comparison of Born-Oppenheimer and Atom-Centered Density Matrix Propagation Methods for Ab Initio Molecular Dynamics H. Bernhard Schlegel Dept. of Chemistry, Wayne State U. Detroit, Michigan 48202, USA

  2. Approaches for Ab Initio Molecular Dynamics • Born-Oppenheimer Dynamics • to move on an accurate PES, converge the electronic structure calculation at each point • Extended Lagrangian Dynamics • propagate the electronic structure as well as nuclei using an extended Lagrangian (e.g. Car-Parrinello, ADMP, etc.)

  3. Born-Oppenheimer approach • to get an accurate PES at each step, converge wavefunction rather than propagate • gradient-based integrators • standard numerical methods • small step sizes • Hessian-based integrators • local quadratic surface - larger steps • Helgaker et al. CPL173, 145 (1990) • predictor – corrector based on a local 5th order surface • CPL228, 436 (1994), JCP111, 3800 and 8773 (1999)

  4. Ab Initio Classical Trajectory on theBorn-Oppenheimer Surface Using a Hessian-based Predictor-Corrector Method Calculate the energy, gradient and Hessian Solve the classical equations of motion on a local 5th order polynomial surface

  5. Logarithm of mass-weighted step size (amu bohr) 1/2 Energy Conservation as a Function of Step Size for Hessian Based Integrators 0 5th order fit (3.00) Rational fit (2.47) -2 2nd order fit (2.76) -4 Logarithm of the average error in the conservation of total energy (hartree) -6 -8 -10 -12 -2.5 -2 -1.5 -1 -0.5 0

  6. Updating methods for Hessian-based trajectory integration • Significant savings in CPU time if the Hessian can be updated for a few steps before being recalculated • BFGS and BSP updating formulas used in geometry optimization are not suitable here • Murtaugh-Sargent (Symmetric Rank 1) update is satisfactory

  7. Relative CPU Time as a Function of the number of Hessian Updates 1.20 C3N3H3 H2CO + CH3F 1.00 NCCHO + CH3Cl 0.80 Relative CPU usage 0.60 0.40 0.20 0.00 0 3 6 9 Number of Hessian updates

  8. Dynamics with Atom-centered Basis Functions and Density Matrix Propagation (ADMP) • Schlegel, Millam, Iyengar, Voth, Daniels, Scuseria, Frisch, JCP114, 9758 (2001), JCP115, 10291 (2001), JCP117, 8694 (2002) • Lagrangian in an orthonormal basis L = ½Tr(VTMV) + ½Tr((m1/4W m1/4)2) - E(R,P)- Tr[L(P2-P)] R, V, M – nuclear coordinates, velocities and masses P, W – density matrix and density matrix velocity in an orthonormal basis (Löwdin or Cholesky) m – fictitious electronic ‘mass’ matrix E(R,P) – energy (electronic + VNN) Tr[L(P2-P)]– constraint for idempotency and N-representability

  9. Euler-Lagrange Equations of Motion • equations of motion for the nuclei M d2R/dt2 = - ¶E/¶R|P • equations of motion for the density matrix d2P/dt2 = -m-½ [¶E/¶P|R+LP+PL-L] m-½ • integrate using velocity Verlet • iteratively solve for L to impose the idempotency constraints on P and W

  10. Behind the Scenes: Choices and Challenges • Basis functions for representing the electronic structure • Method of orthonormalization • Calculation of the energy derivatives • Integration of the equations of motion • Satisfying the idempotency constraints • Mass weighting

  11. Basis functions for representing the electronic structure: plane waves vs atom centered functions • Car-Parrinello method (CP)PRL55, 2471 (1985) • expand in plane waves (appropriate for condensed phase) • most integrals easy to calculate with fast Fourier transforms • only Hellmann-Feynman terms required for gradients • Atom-centered Density Matrix Propagation (ADMP) • expand in atom centered functions (e.g. gaussians) • use methods from standard molecular orbital calculations • far fewer basis functions needed than plane waves • fast routines for multi-centered integrals • gradients require Hellmann-Feynman and Pulay terms

  12. Representing the electronic structure • molecular orbitals vs density matrix • arbitrary rotations among occupied molecular orbitals do not change the energy or the density • density matrices become sparse for large systems and calculations scale linearly • atomic orbitals vs orthonormal orbitals • in an AO basis, changes in the density matrix reflect changes in bonding and changes in overlap • in an AO basis, effect of changes in the overlap must be handled by the propagation • in an orthonormal basis, changes in the overlap are handled separately and only changes in the bonding are handled by propagation • equations of motion simpler in an orthonormal basis

  13. Derivative of the energy wrt the density • energy calculated with the purified density in an orthonormal basis (like CG-DMS method) E = Tr[h P + ½G(P) P] + VNN P = 3 P2 – 2 P3 – McWeeny purified density h, G(P)– one and two electron integral matrices • derivative of the energy with respect to the density in an orthonormal basis ¶E/¶P|R= 3 F P + 3 P F – 2 F P2 – 2 P F P – 2 P2 F F = h+ G(P) – Fock matrix, Kohn-Sham matrix • derivative contains only occupied-virtual blocks (occupied-occupied and virtual-virtual blocks are zero) • equations of motion satisfy idempotency condition to first order

  14. Derivative of the energy wrt the nuclei ¶E/¶R|P = Tr[dh/dRP +1/2 ¶G(P)/¶R|PP] + dVNN/dR where h, G and P are matrices in the orthonormal basis • transformation to the orthonormal basis S’ = UTU, P = UP’UT, h = U-Th’U-1 , etc. where S’, h’, etc. are matrices in the atomic orbital basis ¶E/¶R|P = Tr[dh’/dR P’ +1/2 ¶G’(P’)/¶R|P’P’] + dVNN/dR - Tr[F’ U-1dU/dR P’+ P’ dUT/dR U-T F’] (obtained usingUU-1 = I,UdU-1/dR =dU/dR U-1 and P2=P) ¶E/¶R|P = Tr[dh’/dRP’+1/2 ¶G’(P’)/¶R|P’P’] + dVNN/dR - Tr[F’ P’dS’/dR P’] + Tr{[F, P] (Q dU/dRU-1-P U-TdUTdR)} • for a converged SCF calculation, [F, P]=0

  15. Derivatives of the transformation matrix • for the Löwdin orthonormalization, U= S’1/2 ¶U/¶R = Ssi (si1/2 + sj1/2)-1 (siT¶S’/¶R sj) sjT where s and s are the eigenvalues and eigenvectors of S’ • for the Cholesky basis, S’=UTU and U is upper triangular (¶U/¶RU-1)m,n= (U-T¶S’/¶RU-1)m,n for m < n = ½ (U-T¶S’/¶RU-1)m,n for m = n = 0 for m > n • for large, sparse systems, Cholesky is O(N) and the transformed F and P are sparse

  16. Velocity Verlet step for the density • symplectic integrators provide excellent energy conservation over long periods Pi+1=Pi+Wit - -½ [¶E(Ri,Pi)/¶P|R+LiPi+PiLi-Li] -½t2/2 Wi+½=Wi - -½ [¶E(Ri,Pi)/¶P|R+LiPi+PiLi-Li] - ½t/2 = (Pi+1-Pi)/Dt Wi+1=Wi+½ --½[¶E(Ri+1,Pi+1)/¶P|R+Li+1Pi+1+Pi+1Li+1-Li+1]-½ t/2 • Lagrangian multipliers chosen so that the idempotency condition and its time derivative are satisfied, but must not affect the conservation of energy, etc.

  17. Constraint for Idempotency: scalar mass m DP = P-P0 = W0t--1[¶E(R0,P0)/¶P|R+LP0+P0L-L] t2/2 • P0 is idempotent, choose Lso that P2 = P • W0 and ¶E(R0,P0)/¶P|Rcontain only occ-virt blocks but LP0+P0L-Lcontains only occ-occ and virt-virt blocks • closed form solution, P0DP P0 = -(I-(I-4AAT)1/2)/2, Q0DP Q0 = (I-(I-4ATA)1/2)/2, where A=P0{W0t - -1 [¶E(R0,P0)/¶P|R] t2/2}Q0 • iterative solution, B (A+AT)2 + B2, P0DP P0 = -P0B P0, Q0DP Q0 = Q0B Q0 • Wmust satisfy WP + PW = W W = W* - PW*P - QW*Q where W*=(P-P0)/Dt - -1 [¶E(R0,P0)/¶P|R] t/2

  18. Mass-weighting • in the course of a trajectory, core orbitals change more slowly than valence orbitals • second derivatives wrt density matrix element involving core orbitals are much larger than for valence orbitals • core functions can be assigned heavier fictitious masses so that their dynamics are similar to valence functions • in initial tests,  chosen to be a diagonal matrix  =  for valence orbitals ½ = ½ (2 (Fii + 2)1/2 +1) for core orbitals with Fii < -2 hartree • more general choice may be desirable, depending on the elements and the basis set

  19. Iterative solution for the Lagrangian multipliers for mass weighted case, m • Liis chosen so that Pi+12 = Pi+1 (1) Pi+1= Pi+Wit - -½¶E(Ri,Pi)/¶P|R-½t2/2 (2) P i+1= 3 Pi+12 – 2 Pi+13 (3) dLi = ½ (P i+1- Pi+1) ½ (4) Pi+1= Pi+1+ -½ (PidLi Pi+ QidLi Qi) -½ go to (2) if not converged • Li+1is chosen so that Wi+1 Pi+1+ Pi+1 Wi+1 = Wi+1 (1) Wi+1= (Pi+1-Pi)/t - -½¶E(Ri+1,Pi+1)/¶P|R-½t/2 (2) Wi+1= Wi+1-Pi+1Wi+1 Pi+1 -Qi+1Wi+1 Qi+1 (3) dLi+1= ½ (Wi+1- Wi+1) ½ (4) Wi+1= Wi+1+ -½ (Pi+1dLi+1 Pi+1+ Qi+1dLi+1 Qi+1) -½ go to (2) if not converged

  20. Comparison of Resource Requirements • relative timings can be estimated from the times for computing the Fock matrix, total energy, gradient and Hessian • BO with Hessian based predictor-corrector t = (t(energy) + t(Hessian)) / Dt • BO with Hessian based predictor-corrector with 5 updates t = (t(energy) + 1/6 t(Hessian) + 5/6 t(gradient)) / Dt • BO with gradient based integrators with one gradient eval per Dt t = (t(energy) + t(gradient)) / Dt • ADMP with velocity Verlet integrator t = (t(Fock) + t(gradient)) / Dt • time step, Dt, chosen so that energy is conserved to ca 10-5 Hartree (Dt = ca 0.7 fs for Hessian methods, ca. 0.1 for gradient methods)

  21. Relative timings for HF/6-31G(d) calculations on linear hydrocarbons

  22. Relative timings for B3LYP/6-31G(d) calculations on linear hydrocarbons

  23. Chloride – Water Cluster (0.10 amubohr2 (182 a.u.) fictitious ‘mass’ for the density)

  24. Energy Conservation for Cl- (H2O)25

  25. Comparison of the Fourier transform of the velocity-velocity autocorrelation function for Cl- (H2O)25 at PBE/3-21G*

  26. Fourier transform of the velocity-velocity autocorrelation function for Cl- (H2O)25 at B3LYP/6-31G* (~1 ps) (compare with 3620-3833 cm-1 for OH str in (H2O)2 and 3379 cm-1 in Cl--H2O)

  27. Effect of Fictitious Mass on Vibrational Frequencies • for CP calculations of ionic systems, effective mass of ion is the nuclear mass plus the fictitious mass of the electrons • vibrational frequencies depend on the fictitious mass, but can be rescaled by a factor depending on the electron-nuclear coupling (a) (b) Phonon density of states for (a) crystalline silicon and (b) crystalline MgO P. Tangney and S. Scandolo, JCP116, 14 (2002)

  28. Effect of Fictitious Mass on Vibrational Frequencies • for ADMP calculations, basis functions move with the nuclei • fictitious mass affects only the response of the density, not the dynamics of the nuclei • vibrational frequencies do not depend on the fictitious mass Vibrational motion of diatomic NaCl (fs)

  29. Ab initio classical trajectory study of H2CO ® H2 + CO dissociation. • Important test case, since studied intensively, both experimentally and theoretically. • Excitation to S1 followed by internal conversion to S0, with sufficient energy to dissociate. • Products are rotationally and vibrationally excited. • Born-Oppenheimer ab initio classical trajectory studies: W. Chen, W. L. Hase, H. B. Schlegel, CPL 228, 436 (1994) X. Li, J. M. Millam, H. B. Schlegel, JCP113, 10062 (2000)

  30. Effect of fictitious mass on energy conservation and adiabaticity mvalence = 0.40 amu bohr2 ..... kinetic energy of the density --- total energy __ real energy (total – kinetic) mvalence = 0.20 amu bohr2 mvalence = 0.10 amu bohr2

  31. Effect of the Fictitious Mass on the Vibrational Distributions for CO and H2

  32. Effect of the Fictitious Mass on the Rotational Distributions for CO and H2

  33. Three Body Photofragmentation of GlyoxalXiaosong Li, John M. Millam, H. B. SchlegelJCP 114, 8 (2001), JCP 114, 8897 (2001), JCP 115 6907 (2001) • Under collision free conditions, excitation to S1 followed by internal conversion to S0, with 63 kcal/mol excess energy. • High pressure limit for rate constant: Ea = 55 kcal/mol. • CO formed with Jmax = 42 and a broad distribution, but vibrationally cold • H2 has significant population in v=1, but low J

  34. Transition States for Glyoxal Unimolecular Dissociation Barrier enthalpies (CBS-APNO) glyoxal ® H2 + 2 CO 59.2 kcal/mol H2CO + CO 54.4 kcal/mol HCOH + CO 59.7 kcal/mol 2 HCO 70.7 kcal/mol

  35. Vibrational energy distribution ADMP HF/3-21G

  36. Acknowledgments • Group members • Xiaosong Li, Smriti Anand, Hrant Hratchian, Jie Li, Jason Cross, Stanley Smith, John Knox • Collaborators • S. S. Iyengar, G. E. Scuseria, G. A. Voth, J. M. Millam, M. J. Frisch, W. L. Hase, Ö. Farkas, C. H. Winter, M. A. Robb, S. Shaik, T. Helgaker, • Funding • National Science Foundation, DOE, Gaussian Inc., Wayne State University

  37. Xiaosong Li, John Knox, Hrant Hratchian, Stan Smith, Smriti Anand, HBS, Jie Li

More Related