Fast algorithms for ab initio molecular dynamics John

Fast algorithms for ab initio molecular dynamics John M. Herbert Department of Chemistry University of California Berkeley, CA 94720

Sample timings for – (H 2 O)24 B 3 LYP / 6 -31(1+, 3+)G* (672 basis functions) BOMD timings: (13 SCF cycles) x (380 cpu sec / cycle) = 4940 cpu sec / time step 1 SCF gradient x (1160 cpu sec / grad) = 1160 cpu sec / time step 6100 cpu sec / time step d t ≈ 0. 36 fs ==> 200 cpu days / ps / classical trajectory

Born-Oppenheimer MD Anatomy of an ab initio MD calculation Guess MOs at current molecular geometry Construct Fock matrix, F, from MOs n SCF iterations No Diagonalize F to get new MOs Converged? Yes Calculate energy gradient w. r. t. nuclear coordinates Advance nuclei from t to t + dt

Born-Oppenheimer MD Extended-Lagrangian MD Anatomy of an ab initio MD calculation Guess MOs at current molecular geometry Construct Fock matrix, F, from MOs n SCF iterations No Diagonalize F to get new MOs Start with MOs at time t Construct Fock matrix, F, from MOs Diagonalize F to get new MOs Converged? Yes Calculate energy gradient w. r. t. nuclear coordinates Advance nuclei from t to t + dt Calculate energy gradient w. r. t. nuclei & MOs Advance nuclei from t to t + dt

Extended-Lagrangian MD Anatomy of an ab initio MD calculation Force / 10– 5 a. u. m = 45 a. u. Start with MOs at time t Construct Fock matrix, F, from MOs m = 360 a. u. Diagonalize F to get new MOs Calculate energy gradient w. r. t. nuclei & MOs time / fs Solid lines (color): ELMD force Broken lines (black): BOMD force Advance nuclei from t to t + dt

ELMD: Unsafe at any speed? Vibrational spectra for hydrogen fluoride (Note that 1 a. u. = 0. 024 fs) m = 360 a. u. dt = 14 a. u. m = 180 a. u. m = 90 a. u. dt = 10 a. u. dt = 7 a. u. m = 45 dt = 5 BOMD Vibrational frequency redshifts for diatomic molecules dw / cm – 1 dt / fs w / cm-1 Herbert & Head-Gordon, JCP 121, 11542 (2004)

Accelerating BOMD w/ Fock extrapolation The AO Fock matrix elements are oscillatory functions of time that we fit to a low-order polynomial, Given converged Fock matrices from the previous N time steps, we solve N linear equations for the M + 1 unknown extrapolation coefficients, for each m and n : Call this (N, M) extrapolation. N = number of saved Fock matrices M = degree of extrapolation polynomial Pulay & Fogarasi, CPL 386, 272 (2004)

d. E / millihartree Energy drift SAD = superposition of atomic densities SCF guess

Convergence error in the forces SAD guess SCF thresh = 10– 4 a. u. Force error / 10– 8 a. u. SCF thresh = 10– 6 a. u. (4, 2) extrapolation

Convergence error in the forces (4, 2) extrapolation error / 10– 10 a. u. Force error / 10– 8 a. u. SAD guess error / 10– 8 a. u. (12, 6) thresh = 4 (4, 2) thresh = 6

Extended basis set calculations HF/6 -311 G* simulations of C 2 F 4 at T = 500 K dt = 0. 48 fs thresh = 4 SCF guess NFB thresh = 6 d. E / m. Eh noise NFB drift / ps thresh = 8 d. E / m. Eh noise drift / ps NFB d. E / m. Eh noise drift / ps SAD 6. 0 11. 2 16. 9 9. 6 4. 1 2. 5 13. 3 3. 5 0. 3 Old MOs 2. 0 887. 7 4997. 0 6. 0 3. 7 16. 0 10. 2 3. 6 0. 4 (6, 3) 2. 0 24. 5 65. 8 2. 8 3. 6 8. 9 6. 0 3. 6 0. 2 (12, 6) 3. 1 25. 2 48. 3 3. 1 3. 6 0. 5 4. 1 3. 6 0. 2 (16, 8) 3. 5 26. 1 53. 9 3. 1 3. 6 0. 3 3. 9 3. 6 0. 2

(H 2 O)6 at T = 300 K HF / 6 -31 G* thresh = 6 dt = 0. 24 fs SCF guess NFB B 3 LYP / 6 -31 G* thresh = 6 dt = 0. 24 fs d. E / m. Eh Rel. cpu time noise drift / ps SCF guess NFB d. E / m. Eh Rel. cpu time noise drift / ps SAD 9. 0 1. 00 8. 7 1. 8 SAD 9. 0 1. 00 6. 9 1. 9 (6, 3) 3. 1 0. 58 8. 8 92. 9 (6, 3) 3. 9 0. 60 6. 9 5. 3 (12, 6) 3. 2 0. 56 8. 7 0. 5 (12, 6) 3. 2 0. 55 6. 9 1. 7 (16, 8) 3. 2 0. 57 8. 7 0. 9 (16, 8) 3. 0 0. 54 6. 9 2. 0

[Fe(H 2 O)6]2+ B 3 LYP / 6 -31 G* thresh = 6 dt = 0. 24 fs SCF guess NFB SAD Relative cpu time SCF time/ grad time 13. 2 1. 00 (6, 3) 4. 3 (12, 6) (16, 8) d. E / m. Eh noise drift / ps 2. 71 8. 1 24. 4 0. 50 0. 86 7. 1 18. 4 5. 0 0. 54 0. 98 7. 0 1. 9 5. 3 0. 55 1. 03 7. 1 2. 3

(H 2 O)4– B 3 LYP / 6 -31(1+, 3+)G* thresh = 8; dt = 0. 36 fs 50– 100 trajectories per temperature T/K NFB NSAD d. E / m. Eh noise drift / ps 50 7. 3 ± 0. 9 12. 3 ± 0. 6 2. 7 ± 0. 4 3. 1 ± 4. 8 100 7. 3 ± 0. 8 12. 6 ± 0. 8 5. 4 ± 2. 9 ± 2. 5 150 8. 0 ± 0. 8 12. 9 ± 0. 9 7. 7 ± 3. 5 2. 7 ± 2. 9 200 8. 2 ± 0. 9 13. 1 ± 0. 9 10. 4 ± 5. 4 4. 0 ± 3. 9 250 8. 4 ± 1. 0 13. 4 ± 0. 9 16. 1 ± 8. 5 5. 9 ± 6. 7 300 8. 5 ± 1. 1 13. 5 ± 0. 9 16. 5 ± 8. 6 9. 5 ± 11. 5

Summary • Energy drift in BOMD is an artifact rather than a phenomenon – Use SCF convergence thresh = 10– 6 a. u. for energy conservation • Energy-conserving Fock matrix extrapolation reduces SCF iterations to 3– 5 per time step and reduces overall cost by ~45% in realistic applications • SCF convergence time is now comparable to SCF gradient time ELMD in Gaussian basis sets (2001– 2005)

Curvy-steps Extended Lagrangian MD Gradient d. E / d. Q is unconstrained and can be calculated without self-consistent iteration Herbert & Head-Gordon, JCP 121, 11542 (2004)

ELMD: Unsafe at any speed? Vibrational spectra for hydrogen fluoride (Note that 1 a. u. = 0. 024 fs) adiabatic gap = 5, 727 cm-1 m = 360 a. u. dt = 14 a. u. m = 180 a. u. m = 90 a. u. dt = 10 a. u. dt = 7 a. u. m = 45 dt = 5 BOMD w / cm-1 Herbert & Head-Gordon, JCP 121, 11542 (2004) vib. fundamental vib. overtone lowest “fictons” w / 103 cm-1

ELMD: Unsafe at any speed? Vibrational spectra for hydrogen fluoride (Note that 1 a. u. = 0. 024 fs) m = 360 a. u. dt = 14 a. u. m = 180 a. u. m = 90 a. u. dt = 10 a. u. dt = 7 a. u. m = 45 dt = 5 BOMD w / cm-1 Herbert & Head-Gordon, JCP 121, 11542 (2004) m = 45 a. u. m = 360 a. u.

ELMD: Unsafe at any speed? t(nuc) / t(elec) > 10 t(nuc) / t(elec) < 10 Herbert & Head-Gordon, JCP 121, 11542 (2004)