Tutorial: Quantum Dynamics of Four-Atom Reactions within the Real Wave Packet Framework

1. Tutorial: Quantum Dynamics of Four-Atom Reactions within the Real Wave Packet Framework Stephen Gray Chemistry Division Argonne National Laboratory Argonne, Illinois 60439

2. Outline • Overview • Representation and Evaluation of H q • Initial Conditions, Propagation, Analysis • Computational Experiments

3. Overview • We will learn the “nuts and bolts” of six-dimensional AB + CD  ABC + D reactive scattering • Within the RWP framework but many issues apply to other propagation schemes • System: H2(v1, j1) + OH(v2, j2)  H2O + H using the old, but much studied “WDSE” potential

4. Representation

5. R, r1 : H2 -- OH distance and H2 internuclear distance : Large, evenly spaced grids r2 : OH distance -- remains bound so a PODVR is convenient since it is a (small) set of grid points consistent with a set of vibrational states Rotational basis: j1 = 0, 2, 4, .. ; j2 = 0, 1,.. ; k1 = determined o be consistent. E.g., for K = 0, even parity, k1= 0, 1, .., min(j1, j2)

6. H q = (Td+Trot) q + V q • Use Dispersion Fitted Finite Differences (DFFD’s) to evaluate R and r1 terms in Td If C(iR, i1, i2, j1, j2, k1) denotes the wave packet, the R kinetic energy would involve do over k1, j2, j1, i2, i1 do iR = 1, NR do s = -n, n C’’(iR, ..) = C’’(iR, ..) + d(s) C(iR+s, ..)

7. PODVR for r2 : A small kinetic energy matrix in the PODVR points acts on the r2 part of C Trot: Not strictly diagonal with our rotational basis -- tridiagonal in k1 -- however this is irrelevant in terms of actual numerical effort which is dominated by Td and V Numerical effort for Tdq ? Nrot {N1 N2 NR (2n+1) + N2 NR N1 (2n+1) + NR N1 N22 } = Ntot { 4n+ 2 + N2 }, near linear scaling

8. V q • Diagonal in the three radial distances, so loop over them and (i) transform from rotational basis to angular grid :

9. Multiply by a diagonal V and convert back to basis : Effort : about 2 NR N1 N2 Nj1 Nj2 Nk1 (Nj1+Nj2+Nk1) -- again near linear

10. Propagation and Analysis

11. Analysis Use H2 distance to separate reactants from products : S(r1) = r1 – r1* = 0 Write out q(R, r1*, r2, j1,j2,k1) and q(R, r1, r2, j1,j2,k1)/ r1 | r1=r1* every L time (or iteration) steps

13. Computational Experiments • Incoming Gaussian wave packet centered at R = 10 ao, with e = 0.25 eV and a width, s = 0.3 ao • Reactants in ground vibration-rotation states • r1* = 3 ao for analysis line • Absorption in last 3 ao of R and r1 grids

14. Run 1 : Radial representation R/ ao : 1 – 14, NR = 130 points r1/ao : 0.5 – 6 , N1 = 35 points . 2 PODVR points in r2, based, on diagonalization of a primitive grid Hamiltonian with r2/ao : 1 – 4, 32 points (r2e = 1.85 ao; the PODVR pts are 1.79 and 2.07 ao) Rotational Representation Nj1 = 5 H2 rotational states, j1 = 0, 2, .., 4 Nj2= 10 OH rotational states, j2 = 0, 1, .., 9 A total of 180 rotational states. (The number of angular grid points should be about 10 for each angle.)

15. Reaction Probability

16. Run 1 requires 2.6 hrs on a 1 GHz Linux workstation (and 65 MB) The reaction probabilitiy P(E) is essentially converged over the E = 0.5 to 1.1 eV total energy range (energy relative to separated H2 + OH) General aim of the experiments: to experiment with the representation details to see how they affect P(E) and if smaller grids or bases can be used for some purposes

17. Run 2 : Try just one PODVR point, N2 = 1 all else the same. Result: good, but still rather long for our purposes. (PODVR pt 1.90 ao, CPU, memory both halved.)

18. Run 3: NR = 80 and N1 = 25, N2 = 1 CPU time has been further reduced by 0.5 to 0.6 hrs or 36 minutes. Result -- still reasonably good

19. Run 4 : LikeRun 3 butNj1 = 4 (j2 = 0, 2, 4 and 6) and Nj1 = 9 (j1 = 0,1,..,8) the total number of allowed rotational states decreases from 180 to 94. • Result: CPU time of just 19 minutes. Noticeable high error in P(E) :

20. Run 5 :Like Run 4 but now we have reduced the number of R points from NR = 80 to NR = 60. This leads to about a 33% speedup and the calculation requires about 14 min.

21. Questions and Further Runs • Consider a setup like either Run 4 or Run 5. Experiment further with reductions in the number of grid points in R and r1. • Investigate, with Run 4 (or 5), how the quality of the calculated reaction probability varies with DFFD approximation. • Experiment, with Run 4 (or 5), the role of rotational excitation in the reactants

22. Appendix I: Making and Running the Codes --See also “README” files • To make the propagation program, abcd.x : make -f makefile.abcd • To run it : abcd.x <abcd.run5.in > abcd.run5.out& • Making, running the (flux-based) probability program : g77 -O3 prodflux.f -o prodflux.x prodflux.x <prodflux.run5.in >prodflux.run5.out

23. Appendix II : DFFD files • Subdirectory dffd has various (2n+1) DFFD’s of various overall accuracies e E. g., fd11.e-3 is a (2n+1) = 11 DFFD with accuracy 10-3. See Gray-Goldfield paper for more details [JCP 115, 8331 (2001).]

24. Appendix III: PODVR’s (Echave and Clary; Wei and Carrington) • Define a finite grid representation of some 1D potential problem in x: Ho = To+ Vo • Diagonalize, obtain numerical eigenstates {<xi|v>, Ev} , xi = grid, v = 0, 1,.. • Now represent “x” in a finite vibrational basis, xv,v’ = <v| x | v’>, v,v’ = 0,..,Npo-1 • The eigenvalues of the xv,v’ matrix are the PODVR grid pts

25. Can use, e.g., T0 = H0 - V0 where H0 = Npo x Npo rep. of H0 in PODVR eigenstates and V0= diagonal potential in PODVR to approximate KE operator for the x degree of freedom within the PODVR.