Theoretical Chemical Dynamics Studies for Elementary Combustion Reactions

1. Theoretical Chemical Dynamics Studies for Elementary Combustion Reactions Donald Thompson, Gia Maisuradze, Akio Kawano, Yin Guo, Oklahoma State University Advanced Software for the Calculation of Thermochemistry, Kinetics, and Dynamics Stephen Gray, Ron Shepard, Al Wagner, Mike Minkoff, Argonne National Laboratory • • Interpolating Moving Least-Squares Method (IMLS): • Potential Energy Surface (PES) Fitting Project • Outline • • Motivation • • Method • • Applications

2. Motivation • Potential Energy Surface (PES) • • electronic energy of a molecular system as a function of molecular coordinates • • hypersurface of the internal degrees of freedom • fit to values calculated at discrete geometries • by expensive electronic structure calculations • We seek to develop an automatic PES generator for both: • structured applications • e.g., define PES everywhere below 50 kcal/mol above an equilibrium position • dynamic applications • e.g., on-the-fly calculation of PES for trajectory studies • Automatic PES generation: • Given seed points on the PES, computer alone determines • whether new points are necessary to refine the PES to the input accuracy • what the geometries are where the new points will be calculated Fitting Technique • not high performance computing • directs high performance computational electronic structure calculations

3. Motivation Current most popular PES fitting method: Modified Shepard • unModified Shepard = 0th degree Interpolative Moving Least Squares (IMLS) fit => fits values of PES only => poor derivative properties HN2 example 33 equally spaced points Solid line = exact derivative Spiked line = shepard derivative • Modified Shepard fits value, gradient, and hessian => accurate derivatives • Higher degree IMLS fits are more accurate in value,derivative => Systematic exploration of IMLS for higher degrees

4. Method: • We want V(x) when we know • {V(xi) | i = 1,…,N} calc. ab. initio points • b(x) basis: e.g., b(x) = x = cos(x) = e-ax • Fit by Taylor Series of degree m? • V(x) = Sj=0m aj b(x)j where aj from least squares fit to {V(xi)} • Fit by IMLS of degree m? • V(x) = Sj=0m aj(x) b(x)j where aj from weighted least squares fit to {V(xi)} weights = wi(x,xi), e.g, = e-a∆x2/[∆xn + e] where ∆x = x-xi

5. N m+1 weights basis unknown aj {V(xi)} Method: • a obeys B(xi)T W(x,xi) B(xi) a(x) = B(xi)T W(x,xi) V(xi) • final IMLS fit: V(x) = Sj=0m aj(x) b(x)j => non-linear fit • SVD solution method is best: - more stable - allows reduction in parameters if justified by data • Shepard fit on V (not ∂V/∂x or ∂2V/∂x2) = IMLS fit for m=0

6. Method: • ∂V(x)/∂x = Sj=0m [aj(x) jb(x)j-1∂b(x)/∂x + ∂aj(x)/∂x b(x)] • ∂a/∂x obeys: B(xi)T W(x,xi) B(xi) ∂a(x)/∂x = B(xi)T ∂W(x,xi)/∂x [V(xi)- B(xi)a(x)] unique right hand side same left hand side as equation for a(x) => reuse decomposition of left hand side => direct derivatives (no finite differences) • Shepard has poor derivative properties because 0th IMLS => b(x)0 or derivative of basis does not contribute => only ∂aj(x)/∂x contribute ---sensitive to weights

7. Method: Automatic PES generation • Given some seed ab initio points, can fit method determine: - where to pick next ab initio points - when current fit is converged to a input accuracy • IMLS strategy - reasonable weights mean IMLS fits of all degrees are very close to PES at all ab initio points - away from ab initio points, different degree IMLS differ => let max. difference locate next ab initio point => let minimization of max. difference end generation

8. 100 kcal/mol range Results: 1D Applications Morse Oscillator (MO) 1D slice of HN2 spline PES by Koizumi et al.

9. MO example Equally spaced points IMLS degree 0 1 2 3 4 5 6 7 9 8 cubic spline HN2 example Equally spaced points IMLS degree 0 1 9 8 2 7 3 4 5 6 Results: 1D Applications RMS error in fitting values • compact fit capable of very high accuracy • increasing degree generally increases accuracy • oscillatory behavior at high degree degrades fit • non-linear fit => third degree better than cubic spline

10. IMLS degree 0 1 2 3 4 5 6 7 9 8 MO example Equally spaced points cubic spline IMLS degree 0 9 1 8 2 3 4 5 6 7 HN2 example Equally spaced points Results: 1D Applications RMS error in fitting derivatives • 0th degree (i.e., Shepard) improves poorly with more points • higher degrees have qualitatively improved accuracy

11. Results: 1D Applications Automatic PES generation: 1D Morse Oscillator Example • IMLS degrees for 17 points • max differences where there are no points • contrast of automatic PES generation: 5 seed points + a point at a time where degree difference is maximum to repeated halving of grid increment

12. Results: HOOH 6-D Applications • Tom Rizzo’s 6D HOOH PES • Coordinate representation in terms of 6 interatomic distances ROH, RO’H’, ROH’, RO’H, RHH’, ROO’ • Ab initio sampling - 89 points in the vicinity of HOOH minimum, HOOH hindered rotation barrier, HO--OH reaction path - augmented by Monte Carlo (MC) or Grid sampling up to 100 kcal * MC: (EMS or Random or Combination (EMS+Random) * Grid: Ri = fni Roi for i = 1,6 where f>1 determines increment • RMS error by MC or Grid: Sampling method matters EMS Random GRID COMB

13. Results: HOOH 6-D Applications Fitting directly to PES: Fitting to Differences: • Develop a qualitative fit Vo Apply IMLS to V-Vo • HOOH example - simple functional form - 89 predetermined ab initio + 100 random ab initio pts. EMS-EMS EMS-EMS

14. Results: 6-D Applications • sampling techniques make noticeable differences in rms error • higher degree usually implies higher accuracy • fit cross terms uncoupled to reaction coordinate have negligible effects • fit cross terms coupled to reaction coordiate have noticeable effects

15. Results: 6-D Applications Automatic PES Generation • substantially improves accuracy • works well with modest numbers of seed points • 1 kcal/mol accuracy for 1000 points => 3.2 points/dimension in a 6D grid

16. Results: 6-D Applications Rate constant convergence: • data point selection - 5 points on reaction path - 20 points near HOOH equilibrium - extra points randomly selected • fit: fourth degree with only 6 cross terms • trajectories: - 500 for each case - zero angular momentu • results: - rates from trajectories converge much faster than rms error on the surface

17. Results: up to 15D Applications Model variable dimensional PES V(x1,x2,x3,…) = VEckart(x1) + SiNDOF{VMO(x1,xi)} where VEckart(x1) => where VMO(x1,xi) => number of degrees of freedom rxn barrier width (i dependent) VMO reactants Local Diss. energy (x1 dependent) products x1 = rxn path 0 - Local Diss. Energy + VEckart = fixed global Diss. Energy xi = deviation off rxn path Fit constraints: • fit V < 40 kcal/mol • know turning points at 40 kcal/mol for all xi Parameter values: • 10 kcal/mol barrier for thermoneutral reaction • 100 kcal/mol global dissociation energy • MO width chosen randomly within a range

18. Results: up to 15D Applications • point selection: - on single diagonal (…,xi,…) i = 1,N - on double diagonal (…,xi,…,yi,…) i = 1, N - points accepted if V < Vmax • basis set: Third Degree IMLS without cross terms • Results: - reasonable accuracy - very few points (uniform grid would have very few points per dimension)

19. Results: up to 15D Applications • Effect of cutoff: weight = 0 if weight/max-weight < input limit • cutoff: • • reduces effective # of points • time/evaluation goes linearly • • at extremes, increases error effective # of points

20. Conclusions IMLS: is interesting • PROs - non-linear, flexible, easy extension of Shepard - gradients and hessian not necessary but can be used - efficient direct derivatives - compact, black box code for any dimension PES • user cleverness in basis selection - automatic point selection encouraging - sensitivity to weight selection seems minor • CONs - least squares evaluation every time - every ab initio point “touched” every evaluation unless weight-based screening of points • Future - perfect “black box” code - develop parallel IMLS drivers for electronic structure and trajectory automatic surface generation (collaboration with other SciDAC efforts)