# Presentation: 125085

## 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**

## 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**

## Current most popular PES fitting method: Modified Shepard

**Current most popular PES fitting method: Modified Shepard**

**• ****Higher degree IMLS fits are more accurate in value,derivative**
**  ****=> Systematic exploration of IMLS for higher degrees**

**• *****un*****Modified Shepard **
**  ****= 0****th**** degree Interpolative Moving Least Squares (IMLS) fit**
**  ****=> fits values of PES only**
**  ****=> poor derivative properties**

**• *****Modified***** Shepard fits value, **
**  ****gradient, and hessian **
**  ****=> accurate derivatives**

**Motivation**

**HN****2**** example**
**33 equally spaced points**

**Solid line = exact derivative**
**Spiked line = shepard derivative**

## Method:

**• ****We want V(x) when we know **
**	****•  ****{V(x****i****) | i = 1,...,N}****	****	****calc. ab. initio points**
**	****•  ****b(x)****	****	****	****	****basis:     e.g., b(x)****	****= x**
**	****	****	****	****	****	****	****	****= cos(x)**
**	****	****	****	****	****	****	****	****= e****-********x**** ****	**

**• ****Fit by IMLS of degree m? **
**	****• ****V(x) = ********j=0****m ****a****j****(x) b(x)****j**
**  ****	****	****  ****where a****j**** from *****weighted***** least squares fit to {V(x****i****)}**

**• ****Fit by Taylor Series of degree m?**
**	**
**	****• ****V(x) = ********j=0****m ****a****j ****b(x)****j**
**  ****	****	****  ****where a****j**** from least squares fit to {V(x****i****)}**

**	****	****          ****weights = w****i****(x,x****i****), e.g, = e****-********∆x****2****/[∆x****n ****+ ********]**
**	****	****	****	****	****	****        ****where ∆x = x-x****i**

## Method:

**• ****a obeys****	****B(x****i****)****T**** W(x,x****i****) B(x****i****) a(x) = B(x****i****)****T**** W(x,x****i****) V(x****i****)**
**	**
**	****	**
**• ****final IMLS fit:****	****V(x) ****	****= ********j=0****m ****a****j****(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 ∂****2****V/∂x****2****) = IMLS fit for m=0**

## Method:

**• ∂****V(x)/∂x = ********j=0****m ****[a****j****(x) jb(x)****j-1****∂b(x)/∂x + ∂a****j****(x)/∂x b(x)]**
**• ∂****a/∂x obeys:**
**	**
** ****B(x****i****)****T**** W(x,x****i****) B(x****i****) ∂a(x)/∂x = B(x****i****)****T**** ∂W(x,x****i****)/∂x [V(x****i****)- B(x****i****)a(x)] **

**same left hand side**
**as equation for a(x)**

**unique right hand side**

**  ****=> reuse decomposition of left hand side**

** ****=> direct derivatives (no finite differences)**

**• ****Shepard has poor derivative properties because **
**  ****0****th**** IMLS => b(x)****0**** or derivative of basis does not contribute**
**	****      ****=> only ∂a****j****(x)/∂x contribute ---sensitive to weights**

## Method:

**Automatic PES generation****	**

**• ****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**

**• ****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**

## Morse Oscillator (MO)

**Morse Oscillator (MO)**

**1D slice of HN****2**** spline PES **
**by Koizumi et al. **

**Results: 1D Applications**

## Results: 1D Applications

**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**

## Results: 1D Applications

**Results: 1D Applications**
**RMS error in fitting derivatives**
**• ****0th degree (i.e., Shepard)**
**  ****improves poorly with more**
**  ****points**
**• ****higher degrees have **
**  ****qualitatively improved **
**  ****accuracy**

## 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**

## 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:  R****i**** = f****ni**** R****o****i**** for i = 1,6  where f>1 determines increment**
**• ****RMS error by MC or Grid:**
**  ****Sampling method matters**
**	**** **

**Random**

**EMS**

**COMB**

**GRID**

##     + 100 random ab initio pts.	 

**EMS-EMS**

**EMS-EMS**

**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****.****	**** **

**Results: HOOH 6-D Applications**

**Fitting directly to PES:**
**	**** **

## Results: 6-D Applications

**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**

## Results: 6-D Applications

**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**

## 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**

## Model variable dimensional PES

**Model variable dimensional PES**
**V(x****1****,x****2****,x****3****,...) = V****Eckart****(x****1****) + ********i****NDOF****{V****MO****(x****1****,x****i****)}**
**	**
**     ****where V****Eckart****(x****1****) =>**
**     ****where V****MO****(x****1****,x****i****) =>**

**Parameter values:**
**• ****10 kcal/mol barrier for thermoneutral reaction**
**• ****100 kcal/mol global dissociation energy**
**• ****MO width chosen randomly within a range**

**Fit constraints:**
**• ****fit V < 40 kcal/mol**
**• ****know turning points at 40 kcal/mol for all x****i**

## Results: up to 15D Applications

**• ****point selection:**
**  ****- on single diagonal**
**    ****(...,x****i****,...)  i = 1,N**
**  ****- on double diagonal**
**    ****(...,x****i****,...,y****i****,...) 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)**

## 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**
**  **

## 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)**