Large amplitude quantum mechanics in polyatomic hydrides and multistate electronic
potential energy surfaces of highly electronegative F + HX reactive systems
byMichael Deskevich
B.S., Juniata College, 2000

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
Of the requirement for the degree ofDoctor of Philosophy
Chemical Physics
2007

This thesis entitled:
Large amplitude quantum mechanics in polyatomic hydrides and multistate electronic
potential energy surfaces of highly electronegative reactive F + HX systems
written by Michael Deskevich
has been approved for the Department of Chemistry and Biochemistry________________________________________________
Dr. David J. Nesbitt________________________________________________
Dr. Robert P. Parson

Date_____________

The final copy of this thesis has been examined by the signatories, and we find thatboth the content and the form meet acceptable presentation standards of scholarly
work in the above mentioned discipline.

iii
Deskevich, Michael P. (Ph.D., Chemical Physics)
Large amplitude quantum mechanics in polyatomic hydrides and multistate electronic
potential energy surfaces of highly electronegative reactive F + HX systems
Thesis directed by Professor David J. Nesbitt

In this work two different areas of study are discussed. First, a novel model ispresented for the quantum mechanics of large amplitude motion in floppy hydrides.
A framework is developed for converged quantum mechanical calculations on large
amplitude dynamics in polyatomic hydrides (XHn) based on a relatively simple, but
computationally tractable, “particles-on-a-sphere” (POS) model for the intramolecularmotion of the light atoms. The model assumes independent 2D angular motion of H
atoms imbedded on the surface of a sphere with an arbitrary interatomic angularpotential, which permits systematic evolution from “free rotor” to “tunneling” to“quasi-rigid” polyatomic molecule behavior for small but finite values of total
angular momentum J. Simple tri- and tetra- atom hydrides act as a test suite for the
POS model. After successfully modeling these systems, the model is used on systems
with 4 and 5 hydrogens surrounding a central heavy atom, with the final focus on thetheoretically and experimentally challenging CH5+ molecule.
Next, novel methods are introduced for computing multistate potential energy
surfaces for the abstraction of hydrogen by fluorine in two different experimentally
important systems: F(2P) + HCl Æ HF + Cl(2P) and F(2P) + H2O Æ HF + OH(2Σ). A
novel method of dynamically adjusted weighting factors in state-averaged
multiconfigurational self-consistent field calculations (SA-MCSCF) is developed andiv
tested on the F(2P) + H2O → HF + OH (2Π) reaction. Using the DW-MCSCF
approach a new multistate electronic potential energy surface for the F(2P) + HCl Æ
HF + Cl(2P) reaction is calculated in full dimensionality. The thesis concludes withnonadiabatic quantum nuclear dynamics calculations for the F(2P) + HCl Æ HF +Cl(2P) reaction.

v

Acknowledgments
I would like to acknowledge several people who provided valuable guidanceand support for the work in this dissertation. First and foremost, I would like to thankmy advisor David Nesbitt, whose rigorous standards for excellence in science haveshaped me into a critical observer of all that I encounter. I also wish to thank mycollaborators, including Robert Parson, Anne McCoy, Hans-Joachim Werner, andMillard Alexander, and all of the members of the Nesbitt lab. I am grateful to have
had the opportunity to work with such intelligent people.

vi

Contents
LIST OF TABLES ................................................................................................................................X
CHAPTER 2 .......................................................................................................................................... X
CHAPTER 3 .......................................................................................................................................... X
CHAPTER 4 ......................................................................................................................................... XI
CHAPTER 5 ........................................................................................................................................ XII
LIST OF FIGURES: ........................................................................................................................ XIII
CHAPTER 2 ....................................................................................................................................... XIII
CHAPTER 3 ........................................................................................................................................ XV
CHAPTER 4 ..................................................................................................................................... XVIII
CHAPTER 5 ........................................................................................................................................ XX
CHAPTER 6 ...................................................................................................................................... XXIICHAPTER 1: INTRODUCTION ......................................................................................................25
CHAPTER 2: LARGE AMPLITUDE QUANTUM MECHANICS IN POLYATOMICHYDRIDES: I. A PARTICLES-ON-A-SPHERE MODEL FOR XHN ...........................................30
I. INTRODUCTION ...............................................................................................................................30
II. THEORETICAL BACKGROUND ........................................................................................................35
A. Primitive and coupled basis sets ..............................................................................................36
B. Basis set symmetrization and permutation inversion ...............................................................38
C. Matrix element evaluation .......................................................................................................39
III. TWO PARTICLES-ON-A-SPHERE (XH2) ........................................................................................44
IV. THREE PARTICLES-ON-A-SPHERE(XH3) ......................................................................................47
V. DISCUSSION ...................................................................................................................................51
VI. SUMMARY AND CONCLUSION ......................................................................................................53
TABLES ..............................................................................................................................................55
FIGURES .............................................................................................................................................59

vii
REFERENCES ......................................................................................................................................67
CHAPTER 3: LARGE AMPLITUDE QUANTUM MECHANICS IN POLYATOMICHYDRIDES: II. A PARTICLE-ON-A-SPHERE MODEL FOR XHN (N = 4,5) ............................72
I. INTRODUCTION ...............................................................................................................................72
II. THEORETICAL BACKGROUND ........................................................................................................77
A. Primitive and coupled basis sets .............................................................................................79
B. Basis set symmetrization and permutation inversion ..............................................................80
C. Matrix element evaluation ......................................................................................................81
D. Rigid-body DMC.....................................................................................................................84
III. FOUR PARTICLES-ON-A-SPHERE(XH4) ........................................................................................86
IV. FIVE PARTICLES-ON-A-SPHERE (XH5) ........................................................................................89
V. DISCUSSION ...................................................................................................................................93
VI. SUMMARY AND CONCLUSION ......................................................................................................97
TABLES ..............................................................................................................................................99
FIGURES ...........................................................................................................................................106
REFERENCES ....................................................................................................................................119
CHAPTER 4: DYNAMICALLY WEIGHTED MCSCF: MULTISTATE CALCULATIONSFOR F + H2O → HF + OH REACTION PATHS............................................................................124
I. INTRODUCTION .............................................................................................................................124
II. AB INITIO CALCULATIONS ...........................................................................................................129
III. DYNAMICALLY WEIGHTED MULTICONFIGURATION SELF CONSISTENT FIELDS METHOD (DWMCSCF) ..........................................................................................................................................133
IV. RESULTS ....................................................................................................................................138
V. SUMMARY AND CONCLUSION ......................................................................................................142
TABLES ............................................................................................................................................144
FIGURES ...........................................................................................................................................147
REFERENCES ....................................................................................................................................154

viii
CHAPTER 5: MULTIREFERENCE CONFIGURATION INTERACTION CALCULATIONSFOR THE F(2P) + HCL → HF + CL(2P) REACTION: A CORRELATION SCALED GROUNDSTATE (12A’) POTENTIAL ENERGY SURFACE .......................................................................160
I. INTRODUCTION .............................................................................................................................160
II. DETAILS OF THE AB INITIO CALCULATIONS.................................................................................164
III. BENCHMARK CALCULATIONS.....................................................................................................166
A. Exothermicity .........................................................................................................................166
B. HF, HCl, and ClF asymptotic potentials................................................................................170
C. Transition state ......................................................................................................................170
D. Reaction path .........................................................................................................................172
E. Entrance and exit channel van der Waals wells.....................................................................174
IV. POTENTIAL ENERGY SURFACES .................................................................................................174
A. Two body terms ......................................................................................................................175
B. Three body terms....................................................................................................................176
V. DISCUSSION .................................................................................................................................177
VI. SUMMARY AND CONCLUSION ....................................................................................................183
TABLES ............................................................................................................................................186
FIGURES ...........................................................................................................................................189
REFERENCES ....................................................................................................................................199
CHAPTER 6: MULTIREFERENCE CONFIGURATION INTERACTION CALCULATIONSFOR F(2P) + HCL → HF + CL(2P) REACTION: CORRELATION-SCALED DIABATIC
POTENTIAL ENERGY SURFACES. .............................................................................................205
I. INTRODUCTION .............................................................................................................................205
II. AB INITIO BACKGROUND .............................................................................................................209
III. ADIABATIC SURFACE CALCULATIONS ........................................................................................213
IV. DIABATIC POTENTIAL ENERGY SURFACES .................................................................................217
A. Diabatic Surface and Non-Adiabatic Coupling Calculations ................................................217

ix
B. Spin-orbit Calculations ..........................................................................................................220
C. Analytical fits of the diabatic PES .........................................................................................222
D. Diabatic coupling and spin-orbit fits ....................................................................................224
V. DISCUSSION .................................................................................................................................226
VI. SUMMARY AND CONCLUSION ....................................................................................................234
FIGURES ...........................................................................................................................................236
REFERENCES ....................................................................................................................................247

x

List of tables
Chapter 2
Table 2.1. Basis set size and matrix calculation time scaling as a function of n, J, andjmax.

Table 2.2. Calculated and experimental H2O energy levels (cm-1) for the bending
mode and rotational constants (cm-1).

Table 2.3. Calculated and experimental NH3 energies and splittings.

Table 2.4. Calculated and experimental H3O+ energies and splittings.

Chapter 3
Table 3.1. The Legendre coefficients for the α = 1 CH4 and CH5+ potentials

Table 3.2. Time required to create an orthonormal basis set of one irreducible
representation in the G240 molecular symmetry group.

Table 3.3. First five J = 0 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.

Table 3.4. First five J = 1 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.

xi

Table 3.5. First five J = 2 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.

Table 3.6. The expectation value of the energy of a 1D rotor in the CH5+ α = 1
pairwise potential.

Table 3.7. The error in ground state energies of the POS model and the 1D Legendre
expansion of the pair correlation function for the converged n = 3 NH3 calculation.

Chapter 4
Table 4.1. Product, reactant energies, and reaction exothermicities as a function oflevel of theory and basis set size.

Table 4.2. Product, reactant energies, and reaction exothermicities as function of
optimizing the AVTZ basis set

Table 4.3. The geometry of the transition state at the SA-MCSCF and MRCI+Q level.
The lengths and angles are defined in Figure 4.1. Bond lengths are in angstroms,
angles in degrees. The basis set used is AVTZ[6s4p3d]/VTZ[3s2p] in a full valence
active space.

xii

Chapter 5
Table 5.1: ΔErxn in kcal/mol of F(2P) + HCl Æ HF + Cl(2P) calculated by various ab
initio methods. Calculated energies are corrected for zero-point and spin-orbit which
are not included in the single surface ab initio calculations. The experimental ΔE is-33.060 +/- 0.001 kcal/mol.

Table 5.2: Diatomic constants from ab initio PES calculated at MRCI+Q/AVnZ
compared with the correlation scaled values and experiment.

Table 5.3: The geometries and energies of the transition state for a range of methods
and basis sets. In particular, the transition state geometry is strongly bent at all levelsof theory.

xiii

List of figures:
Chapter 2
Figure 2.1. Schematic picture of the particles-on-a-sphere (POS) method, which isbased on a modeling of polyatomic hydrides (XHn) as one heavy atom (X) that
constrains the motion of n H atoms to the surface of a sphere.

Figure 2.2. EJ = 0 energy level correlations for XH2, plotted for the first few lowest
energies vs. potential stiffness, χ. At small χ, a simple rigid rotor pattern withincreasing spacings (4b, 8b, 12b, etc.) is clearly exhibited (the rotational levels are
reproduced with greater detail in the inset), but it evolves smoothly to a harmonic
pattern of constant energy spacing (hω0) in the limit of large χ.

Figure 2.3. EJ = 0,1,2 energy level correlations for XH2 in the ground vibrational state,plotted for the first few lowest energies versus potential stiffness, χ. At small χ, theenergy separations reflect uncorrelated free rotation of the hydrogen atoms. At large
χ, however, this smoothly transitions into the characteristic JKaKc patterns of a rigid
asymmetric top.

Figure 2.4. 1D wave function probability distribution for the XH2 ground vibrational
state plotted as a function of increasing χ. As the potential stiffness is increased, thewidth of the HXH angular distribution narrows toward γ = γ0, i.e., the equilibriumxiv
geometry. As shown in the inset, the first hydrogen is defined to be at (θ = 0, φ = 0),
with the |amplitude|2 for the second hydrogen plotted for (θ = γ, φ ).

Figure 2.5. EJ = 0 energy level correlation as a function of intermolecular stiffness for
XH3, plotted for the first few energies of the umbrella mode. At small χ, the levels
correspond to three independent rigid HX rotors (displayed with greater detail in theinset) but eventually correlate with inversion doubled tunneling splittings in the limitof large χ.

Figure 2.6. Ev = 0 energy level correlation for XH3 J = 0, 1, 2 as function ofintermolecular stiffness, χ. For small χ, the pattern reflects three freely rotatinghydrogens, which for large χ smoothly transitions into the rotational energy patterncharacteristic of an oblate symmetric top.

Figure 2.7. 2D cuts (in θ, φ space) through full 6D XH3 wave functions for (i)symmetric ground state (left hand side) and (ii) the antisymmetric first excited state(right hand side) of the umbrella mode for increasing values of χ. The contour plot
(center bottom) displays the corresponding potential. At χ = 6, 60 and 600, thecontour spacings are 1.5 cm-1, 15 cm-1, and 150 cm-1, respectively. As shown in theinset, the first and second hydrogens are fixed at (θ = 0, φ = 0) and (θ = γ0, φ = 0),
with amplitude for the third hydrogen plotted as a function of (θ , φ ).

xv
Figure 2.8. The log of the fractional error (base 10) as a function of jmax for several
values of stiffness χ. Larger χ converges more slowly as the basis set is optimized for
floppy intermolecular potentials. For a given χ, it is worth noting that the n = 2 and n= 3 calculations converge almost at the same rate, which bodes favorably for scaling
to a larger number of H atoms.

Chapter 3
Figure 3.1. Particle-on-a-sphere rotation-bending energies for a) n = 2 (XH2) and b) n
= 3 (XH3) model systems as a function of basis set size (jmax). a) In XH2, theasymmetric top rotational patterns for J = 0,1 corresponding to JKaKc = 000, 101, 111, 110
are clearly evident, with the J = 0 energy rapidly converging to the expected DMC
limit (dotted line). b) Converged J = 0 energy levels for XH3 reveal inversion
tunneling splittings in good agreement with NH3 (after which the potential wasmodeled) and increasing rapidly with symmetric bend excitation. Note the only
exception to this pattern (near 5000 cm-1), arising from an interloper vibrational state
due to degenerate asymmetric bend excitation, which exhibits a tunneling splitting
very similar to the ground state (0.798 cm-1).

Figure 3.2. Angular pair correlation functions for a) XH2 and b) XH3 as a function ofpotential stiffness χ = k/b, where k is the angular well depth and b is the rotational
constant for a single X - H rotor. For comparison, the shaded function is calculated
from DMC, whereas the dots reflect the fully converged J = 0 POS calculation.

xvi
Figure 3.3. Convergence of J = 0, 1, 2, 3, 4 energies for XH4 as a function of jmax.
The POS ground state (J = 0) energy is compared against DMC predictions (dottedline).

Figure 3.4. a) XH4 rotational energies of J = 0 - 4 as a function of potential stiffness
scale factor α.

b) At higher resolution, centrifugal induced splittings of thenominally degenerate symmetric top levels for XH4 are shown for α = 0.033, which
are considerably larger but reveal the same qualitative patterns as experimentallyobserved.

Figure 3.5. Vibrational energies of XH4 as a function of the potential stiffness scale
factor α.

Figure 3.6. Minimum energy paths for CH5+ through the Cs and C2v transition states,
separating the Cs minimum for both the 10D relaxed potential and the pairwise fit.

Figure 3.7. Geometry comparison of CH5+ critical points. Note the surprisingly
quantitative agreement between the full 10D relaxed potential and the pairwise
additive fit used in the POS model.

Figure 3.8. Pair correlation functions for J = 0, XH5 from DMC calculations for aseries of potential models: a) 15D full CH5+ potential, b) the 10D relaxed potential,xvii
and c) the pairwise additive least squares fit. Note the excellent agreement at thewave function level.

Figure 3.9. Ground-state energies of CH5+ as a function of the potential stiffness
parameter (α), calculated by DMC and POS with both the 10D (DMC, open squares)
and pairwise (POS, dashed line) potentials. The inset shows the difference between
the POS and DMC (pairwise) calculations (i.e., the level of convergence) as afunction of α.

Figure 3.10. Nonzero statistical symmetries for J = 0, 1, and 2 of XH5 (with a = .033)
as a function of basis set size, jmax and Ncoupled (the number of coupled basis functionsin the basis set). Note the slow but steady convergence of J = 0 energies to the
RBDMC limit with the size of Ncoupled.

Figure 3.11. The lowest J = 0 internal rotor energies (with respect to the ground stateenergy, G2- Symmetry) of the POS model compared with the RTRF-like model of
Bunker and coworkers. Interestingly, the splittings for the two model treatments areroughly the same order of magnitude, despite the stiffness of the potentials varying bya large factor.

Figure 3.12.

DMC pair correlation function of CH5+ as a function of potential

stiffness scale (α).

xviiiFigure 3.13. Successive fits of the α = 1 pair correlation function of CH5+ to a sum ofLegendre functions. Such behavior suggests a much improved level of convergence
for the POS model of CH5+ by jmax = 7 - 9.

Chapter 4
Figure 4.1: The internal coordinates used for specifying F + H2O geometries. The
geometry is given for the 5-state MCSCF transition state.

Figure 4.2: A 2D projection (RFH1 and ROH1 bond distances) of the full 6D reaction
path calculated at the 5-state MCSCF level. This 6D reaction path is used throughout
the paper.

Figure 4.3: Calculated ground and excited state energies along the F + H2O reaction
coordinate for various levels of theory, a) 3-state MCSCF, b) 5-state MCSCF, c)
Dynamically weighted MCSCF.

The inserts at reactant and product asymptotes

indicate a magnification of the energy scale to highlight deviations from perfect 3fold (reactant) and 2-fold (product) degeneracies. Note that the dynamic weighting
method produces the smoothest surface, and also reproduces the correct degeneracy atboth asymptotes of the reactants.

Figure 4.4: Dependence of ground state reaction path energy on weighting function in
the transition state region. a) Exponential, Gaussian or sech2 functional form. Notethat the choice of the functional form has less of an effect on the surface than thenumber of states included in a traditional SA-MCSCF. b) Choice of scale parameter

xix
β.

If β-1 is too small, small changes in energies in excited states can cause

discontinuities in the ground state (e.g. point P), if β-1 is too large, the weights do notdecay to 0 in the asymptotes for the excited states. Moderate β-1 parameters (i.e.,2 3eV) produce smooth curves and converge to ground-state-only asymptotes.

Figure 4.5: The lowest 4 states of the smooth DW-MCSCF reaction path. Theenergies along the reaction path not only vary smoothly, but the degeneracies arecorrectly reproduced in the asymptotes. Also note that there are a number of statesinteracting in the transition state region.

Figure 4.6: Converged weights (normalized to sum to 1) for the DW-MCSCF with β-1= 3. The energies along the reaction path are reproduced at the top to help guide the

1 1 1
1 1
eye. The weights smoothly vary from ( : : ) to approximately ( : :0). During
3 3 3
2 2
the cross-over through the transition state a fourth and fifth state are briefly
introduced. DW-MCSCF automatically included the relevant states in the different
regions of the PES.

Figure 4.7: MRCI+Q with spin-orbit interactions calculated at the CASSCF reaction
path geometries. The smoothly varying DW-MCSCF wave functions are used as
reference configurations for MRCI+Q calculations.

Note that the lowest barrier(when zero-point energy is included) is experimentally accessible in crossed
supersonic jets, but the other barriers are far out of reach.

xx

Chapter 5
Figure 5.1: AVnZ (n = 2 - 5) extrapolation for F(2P) + HCl and HF + Cl(2P) energies
to the CBS limit, indicating clear convergence (see inset) in the exothermicity (ΔE) as
a function of n. The residual error in ΔE (dashed line) is due to active space size in
the MRCI+Q calculations, as included in correlation scaling.

Figure 5.2: Correlation energy (EMRCI - EMCSCF) recovered for each basis set, AVnZ.
The nearly constant ratio of recovered correlation energy relative to the CBS
extrapolation (see inset) illustrates how correlation scaling effectively mimics a larger
active space.

Figure 5.3:

F-H-Cl transition state geometry at various levels of theory, allconfirmingasharplybenttransition

state:(a)

MRCI+Q/AVTZ;(b)

MRCI+Q/CBS/Scaled; UCCSD(T)/AVQZ.

Figure 5.4: The DW-MCSCF energies (a) and weights (b) as a function of reaction
coordinate, with up to 6 states dynamically weighted. Note the strongly avoided
crossings of F- + HCl(2Π)+ and HF(2Π)+ + Cl- charge transfer (CT) states withsurfaces correlating to the ground state (2P) asymptotes. Note also the smooth weightadjustment from a 3-state calculation at either asymptote to > 5 states near thetransition state geometry.

xxi
Figure 5.5: 12A’ reaction path at the MCSCF, MRCI+Q, and correlation scaled levels,
indicating relatively minor effects due to correlation scaling. Such scaling makes up
for 2.18 kcal/mol in ΔE and lowers the barrier by 0.43 kcal mol.

Figure 5.6: 2D slice (θ = 123.5o) of the fitted 12A’ surface with global RMS of 0.79
kcal/mol.

Figure 5.7:

Ground reference state coefficient in the MRCI+Q wavefunction,

indicating strong admixture of excited state character near the transition stategeometries closest to the conical intersection (marked by asterisks).

Figure 5.8: 12Α’ PES RHF vs. RHCl contours for (a) transition state (θ = 123.5o) and(b) collinear (θ = 180o) bend angles. Conical intersection seams occur to each side of
the transition state; the seam location closest to the reaction path is marked by theasterisk. Contour spacing is 5 kcal/mol with respect to F + HCl(Re).

Figure 5.9: F-H-Cl bending potential corresponding to transition state (dashed) andconical intersection seam (solid) RHF and RHCl bond lengths, with eigenvalues
obtained from a rigid bender analysis. Note that the conical intersection is lower than
the transition state with only a ∼1 kcal/mol barrier to linearity from the reaction path.

Figure 5.10: a) 2D effective transition state barrier height (solid) and HCl bond length(dashed) as a function of F-H-Cl bending angle (dotted line is Re for HCl). As

xxiiθ decreases below ∼120o, both the HCl bond length and barrier height increase

dramatically from transition state values, characteristic of a shift from an “early” to“late” transition state. b) Trajectory analysis for F-HCl angular deflection (Δθ) as afunction of initial bend angle, θ. For J = 0, only a very narrow range (dark grey) of
incident scattering angles are successfully “prerotated” (light grey) into theappropriately bent transition state angle to react.

Chapter 6
Figure 6.1: Three lowest adiabats along the reaction path for the F(2P) + HCl Æ HF
+ Cl(2P) system, obtained by high-level DW-MCSCF, MRCI+Q and CBS methods.
The lowest two correlate with ground state F and Cl, while the highest one correlates
asymptotically with spin-orbit excited F* and Cl* species. In the adiabatic limit, only
the lowest spin-orbit ground state can react over a barrier accessible at experimental
collision energies.

Figure 6.2: 2D slices through two full 3D A’ adiabats. a) The region near the conicalintersection seam (solid black line) at θ=180o, which by symmetry becomes b) an
avoided crossing for any non-collinear geometry (θ=170o).

Figure 6.3: A 1D angular cut through full 3D adiabatic and diabatic surfaces in theexit channel, with RHF=Req and RHCl=6.6 a0.

The adiabats and diabats are

asymptotically equal at collinear geometries, but strongly avoid as the symmetry is
lowered from C∞v.

xxiiiFigure 6.4: Adiabats (dashed lines) and diabats (solid lines) calculated along thereaction path a) in full 3D and b) with θ constrained at 180o. In full 3D the adiabats
avoid where the diabats cross, whereas for a collinear geometry the adiabats anddiabats overlap perfectly at each point along the conical intersection seam due to zero
coupling.

Figure 6.5: A 2D RHF, RHCl contour plot for the ground state adiabatic surface at botha) transition state (θ=123o) and b) collinear (θ=180o) geometries. In b) the conicalintersection seams are shown with dashed lines. Contour spacing is 3 kcal/mol withrespect to zero at reactant entrance channel.

Figure 6.6: Grid of distributed Gaussian locations used to fit the β coupling surface.
Angular dependence is fitted at each point (RHF, RHCl) by a linear combination ofsin(nθ) functions.

Figure 6.7: 2D contours (at θ=123o) for the a) VΣ, b) VΠx, and c) VΠy diabats with 3
kcal/mol spacing. A sample histogram of residuals for fitting the VΣ function is also
shown in a). The values are not distributed normally, which is why the RMS (0.34
kcal/mol) is considerably larger than the FWHM (0.049 kcal/mol).

xxiv
Figure 6.8: Sample 2D slices through the full 3D β coupling surface between VΣ andVΠx diabats at a) θ=180o, b) θ=150o, and c) θ=120o. Note that for collinear geometry,
β vanishes identically but grows rapidly as the molecule becomes bent.

Figure 6.9: VΣ and VΠx diabats along the reaction path (scale to left), and the βcoupling between the diabats (scale to right). Note the strong peaking of diabatic
coupling immediately in the post transition state region, due to rapid change inelectronic character for the newly formed bond.

Figure 6.10: Sample results starting at the transition state (p=0 au, Δp=2 au) for
multistate 1D wavepacket propagation along the F + HCl → HF + Cl reaction
coordinate. The top panel displays adiabatic energies, the middle panels reveal
snapshots of |Ψ|2 out to 100 fs, and the bottom panel accounts for the total flux thatpassed by each point on each surface.

Figure 6.11: Reaction dynamics for ground F and spin-orbit excited F* as a functionof energy. a) Reaction probability of F and F* along the reaction path (solid lines),
with dotted lines representing probability if non-adiabatic coupling is set to 0. b)
Reaction probability of F and F* along a 2D reaction path with θ constrained to be
180o. The dramatically increased reactivity of F* vs F results from the presence of
conical intersection seams sampled at collinear geometries.

25

Chapter 1: Introduction
In the field of theoretical chemical physics, the main focus of research
generally falls along one of two paths: 1) the exact study of the very small or 2) the
ensemble study of the very large. In the former, the systems are small enough (fewer
than 3 particles) that the mathematics allow one to reach exact results. In the large
systems (greater than 50 - 100 particles), there are generally enough degrees of
freedom that measurements of the bulk can accurately be described by statistical
averages where the exact details of each particle are not relevant. These two branchesof study leave behind a middle ground of systems that are too big to study exactly, yet
are too small to be treated statically. The key to studying these middle-of-the-roadsystems is to find the correct approximations that capture most of the chemistry andphysics, while still being computationally tractable and providing the correct simple
physical picture. In this work two different areas of study are presented, both ofwhich fall into this middle ground.

First, I present a model for the quantum

mechanics of a floppy 6-atom system (CH5+) which aims to capture the low-energy
dynamics of a system with a large number of large-amplitude degrees of freedom.
Then I will introduce novel methods for computing multistate potential energy
surfaces for the abstraction of hydrogen by fluorine in two different experimentally
important systems: F(2P)+HCl Æ HF + Cl(2P) and F(2P)+H2O Æ HF + OH(2Σ).
For the CH5+ problem, I start by a developing framework for converged
quantum mechanical calculations on large amplitude dynamics in polyatomic
hydrides (XHn) based on a relatively simple, but computationally tractable, “particleson-a-sphere” (POS) model for the intramolecular motion of the light atoms. The

26
model assumes independent 2D angular motion of H atoms imbedded on the surface
of a sphere with an arbitrary interatomic angular potential, which permits systematic
evolution from “free rotor” to “tunneling” to “quasi-rigid” polyatomic moleculebehavior for small but finite values of total angular momentum J. Chapter 2 focuseson simple triatom (n = 2) and tetratom (n = 3) systems as a function of interatomic
potential stiffness, with explicit consideration of H2O, NH3, and H3O+ as limiting test
cases. The POS model also establishes the necessary mathematical groundwork for
calculations on dynamically much more challenging XHn species with n > 3 (e.g.,models of CH5+), where such a reduced dimensionality approach offers prospects forbeing quantum mechanically tractable at low J values (i.e., J = 0, 1, 2) characteristicof supersonic jet expansion conditions.
In Chapter 3, the POS framework focuses on systems with many degrees of
freedom (i.e., n = 4 and 5). The goal is to probe and elucidate the large amplitude
dynamics in low J states (i.e., J = 0, 1, and 2) that will dominate high resolution
spectra of such hydrides (e.g., CH4, CH5+) obtained under supersonic jet-cooled
conditions. Comparison with experimental data as well as Diffusion Quantum Monte
Carlo methods show how this simple reduced-dimensionality model provides
potential insight into multidimensional quantum rovibrational dynamics.
Switching gears from quantum dynamics of nuclei to the creation of potential
energy surfaces (on which quantum dynamics of nuclei can be modeled), Chapter 4
discusses a novel method of dynamically adjusted weighting factors in state-averaged
multiconfigurational self-consistent field calculations (SA-MCSCF) which isapplicable to systems of arbitrary dimensionality. The proposed dynamically27weighted approach automatically weights the relevant electronic states in each regionof the potential energy surface, smoothly adjusting between these regions with anenergy-dependent functional. This method is tested on the F(2P) + H2O → HF + OH(2Π) reaction, which otherwise proves challenging to describe with traditional SAMCSCF methods due to i) different asymptotic degeneracies of reactant (3-fold) andproduct (2-fold) channels, and ii) presence of low lying charge transfer configurations
near the transition state region. The smoothly varying wave functions obtained bydynamically weighted multiconfigurational self-consistent field (DW-MCSCF)
methods represent excellent reference states for high-level multireference
configuration interaction (MRCI) calculations and offer an ideal starting point for
construction of multiple state potential energy surfaces.
The F + H2O system is still too large to study in full dimensionality, however,the novel DW-MCSCF approach allows one to accurately study smaller reactive
systems.

Using this approach, Chapter 5 presents a new ground state (12A’)electronic potential energy surface for the F(2P) + HCl Æ HF + Cl(2P) reaction. The
ab inito calculations are done at the MRCI+Q/CBS level of theory by extrapolation of

the MRCI+Q/aug-cc-pVnZ (n = 2,3,4) energies. Due to low-lying charge transfer
states in the transition state region, the molecular orbitals are obtained by 6-state DWMCSCF methods. Additional perturbative refinement of the energies is achieved by
implementing simple one-parameter correlation energy scaling to reproduce theexperimental exothermicity (ΔE = -33.06 kcal/mol) for the reaction. Ab initio points
are fit to an analytical function based on sum of 2- and 3-body contributions, yielding
an RMS deviation of < 0.3 kcal/mol for all geometries below 10 kcal/mol above the28barrier. Of particular relevance to non-adiabatic dynamics, the MCSCF-MRCI
calculations show significant multireference character in the transition state region,
which is located 3.8 kcal/mol with respect to F + HCl reactants and features astrongly bent F-H-Cl transition state geometry (123.5o). Finally, the surface also
exhibits two conical intersection seams that are energetically accessible at low
collision energies. These seams arise naturally from allowed crossings in the C∞v
linear configuration that become avoided in Cs bent configurations of both thereactant and product and should be characteristic of all X-H-Y atom transfer reaction
dynamics between (2P) halogen atoms.
Finally, in Chapter 6, the DW-MCSCF approach is used in creating a multielectronic state PES for the F(2P) + HCl Æ HF + Cl(2P) reaction. This chapter
presents benchmark ab initio combined calculations and analytical fits for the lowest
three adiabatic (12A’,22A’,12A”), diabatic (12Σ, 12Πx, 12Πy) and spin-orbit corrected
potential energy surfaces associated with the fundamental F(2P) + HCl Æ HF +Cl(2P) reaction. These surfaces are based on high-level MOLPRO multireference
configuration interaction (MRCI) calculations performed at the complete basis set(CBS) level by extrapolation of MRCI+Q/aug-cc-pVnZ (n = 2,3,4) energies as afunction of RHF, RHCl, and θFHCl . Molecular orbitals for the MRCI calculations are
obtained by 6-state dynamically weighted MCSCF (DW-MCSCF) methods,
necessitated by the presence of low-lying (F- + HCl+ and HF+ + Cl-) charge transfer
states in the transition state region. Additional perturbative refinement of the adiabatic
potentials is achieved by implementing correlation energy scaling methods toreproduce the experimental exothermicity. Adiabatic surfaces are transformed into a29smooth diabatic representation by explicit maximization of orbital overlap withrespect to limiting θFHCl . = 0o and 180o configurations, which reveals seams ofconical intersections in both entrance and exit channels and which is energeticallyaccessible at typical reactive collision energies. The resulting diabatic and nonadiabatic coupling surfaces are analytically least-squares fitted with 2- and 3-bodyterms, with further inclusion of the Breit-Pauli matrix operator yielding spin-orbit
splittings in good agreement with experiment.

By way of test application,

wavepacket calculations along the reaction path for the set of 3D multiple surfaces
are used to predict the role of non-adiabatic surface hopping probabilities for reactionof F(2P 3/2, 1/2) + HCl to form Cl(2P 3/2, 1/2) + HF.

30

Chapter 2: Large amplitude quantum mechanics inpolyatomic hydrides: I. A particles-on-a-sphere modelfor XHn
I. Introduction
The quantum mechanics of large amplitude motion have become an increasing
focus of experimental and theoretical attention, fueled in part by impressive advancesin high-resolution rotational/vibrational spectroscopy of tunneling dynamics inmolecules, ions, and hydrogen bonded clusters. Classic paradigms in the field of large
amplitude quantum dynamics have been systems such as tunneling inversion inammonia1-8 (NH3) and hydronium ion9-16 (H3O+), as well as numerous demonstrationsof large amplitude tunneling phenomena in van der Waals clusters.17-20 In favorably
small n-atom systems, the dynamics can be treated at the level of exact quantum
calculations for nuclear motion in full (i.e., 3n-6) dimensionality. This has been
elegantly demonstrated, for example, in studies of HF and HCl dimers,21-24 whichhave permitted converged variational calculations in all inter- (4D) andintramolecular (2D) coordinates, as well reduced dimensionality extensions to H2O
dimer25-28 (6D intermolecular coordinates with semirigid monomers) for which
multiple H atom tunneling pathways of the rigid H2O monomers are clearly evident.
The level of quantitative success in such benchmark hydrogen-bonded clusters isextremely encouraging and bodes well for understanding tunneling dynamics in evenmore complicated polyatomic systems.

31
However, as these large amplitude molecular systems increase in size andcomplexity, they quickly become prohibitive to treat by standard quantum calculation
methods. This is unfortunate, as some of the most interesting dynamical processes
begin to be possible in these larger systems.

Examples include the proposed

concerted “race track” tunneling motion of the 3 H atoms around protonatedacetylene29 (C2H3+), and, as an even more longstanding spectroscopic challenge,
facile intramolecular H atom exchange in protonated methane.30, 31 Furthermore, thislevel of theoretical difficulty increases rapidly with J, which makes the interpretationof rotational patterns especially challenging under all but the lowest temperature
conditions. In such cases, a more realistic first spectroscopic goal is to learn how to
predict and interpret the experimentally observed energy level patterns and thereby
infer new insights about the actual intermolecular dynamics.
As background motivation for the present theoretical study, protonatedmethane (CH5+) represents an extreme case of large amplitude quantum dynamics
with a rich scientific history. Adding a single proton to the well-studied methane,
CH4, creates a highly fluxional molecule where simple Lewis octet-bonding motifs do
not apply, resulting in 5 H atoms connected by 4 electron pairs and requiring thepresence of “three-center-two-electron” bonds (3c-2e).32 Such 3c-2e bonding motifs
correspond to a special class of hypercoordinated carbocations, extremely important
reactive intermediates in acid-catalyzed electrophilic reactions,33, 34 for which CH5+
represents the simplest prototypic “superacid.” CH5+ is also believed to be involved in
the synthesis of polyatomic species in cold interstellar clouds,35, 36 fostering
astrophysical interest in this simple but spectroscopically elusive molecular ion.

32
From a theoretical perspective, CH5+ is interesting because its small size
makes quantum mechanical calculations of the potential computationally feasible.
Early calculations of the potential surface37-39 predicted the lowest energy equilibrium
configuration to exhibit CS symmetry. Later studies40 predicted an equilibrium
structure reminiscent of CH3+ “solvated” by a more distant H2 moiety, but as the levelof theory has improved, the equilibrium separation between the CH3+ and H2 hassteadily decreased. Most recent high-level calculations suggest that there are, in fact,
three low-lying energy structures within about 1 kcal/mol of each other.30, 41-43
Furthermore, if zero point energy is taken into account, these three configurations are
all extensively sampled by the ground-state wave function. Stated semiclassically, thebarrier to rearrangement between these low-lying minima is lower than the zero point
vibrational energy, thus promoting facile intramolecular scrambling of the hydrogenatoms.41
From a spectroscopic perspective, the highly fluxional, nonclassical nature of
CH5+ begins to account for long-standing difficulties in obtaining and interpreting its
high resolution spectrum. Indeed, despite its first observation as a highly abundant ion
in mass spectrometers44 in the early 1950s, optical detection and characterization of
CH5+ eluded spectroscopists for another 50 years. The breakthrough came in 1999
from Oka and co-workers, who obtained a spectrum45 in the CH stretch region
exploiting velocity modulation methods and made convincing arguments that theextensive, albeit unassigned, spectrum belonged predominantly to CH5+. Desire tohelp in this assignment process has naturally led to considerable emphasis on
calculating the near IR spectrum from theoretical first principles. However, despite

33
intense theoretical efforts 29, 30, 37-41, 46-53 directed toward CH5+ and the recent
availability of a high-level potential surface,42, 43 even a qualitatively correct, highresolution spectrum based on the rovibrational energy levels has proven challenging
to calculate.
The reasons for this are at least twofold. On the experimental side, large zero
point energies and lack of any substantial barriers between H atom interchange leadto extensive delocalization in the wave function. This delocalization, in turn, canresult in large tunneling splittings and rovibrational energy patterns profoundly
perturbed away from rigid rotor expectations, especially for the J states thermally
populated under discharge conditions. From a theoretical perspective, however, thenumber of degrees of freedom for CH5+ is already too large (i.e., 3n-3 = 15) toachieve fully converged quantum calculation of the rovibrational energy levels in full
dimensionality.
In this chapter and the next, I explore a conceptually simple and physicallymotivated particles-on-a-sphere (POS) model for studying large amplitude quantum
dynamics in systems as large as CH5+. This model is stimulated by ab initio
calculations by both Bunker50, 51 and Bowman,42, 43 which reveal that each of the C-Hbond lengths in CH5+ remain nearly constant throughout the full manifold of H atom
exchange pathways. Specifically, the distribution of RCH distances on the full
dimensional potential surface differ only by ~0.015Å between equilibrium and saddle
point configurations, respectively, which is already comparable to zero point
displacements (~0.018 Å) in the corresponding radial CH stretch coordinate.54 This
indicates a high degree of decoupling between radial stretch and bend/rotation

34
coordinates, which I exploit by focusing explicitly on the internal CH bend/rotational
degrees of freedom. Simply stated, the model considers the large amplitude quantum
dynamics of an effective XHn molecule (where mX >> mH) where n identical H atoms
are radially constrained to the surface of sphere, with angular motion dictated by agiven H-H interaction potential, thus motivating the POS description (see Figure 2.1).
This approach reduces the dimensionality of the problem to ∼2n (for a sufficiently
massive central atom) and therefore makes numerically converged bending/internal
rotor/tunneling calculations relatively straightforward for triatomic (XH2) andtetratomic (XH3) hydrides, as discussed herein. Most importantly, the method makes
calculations for J = 0,1,2 computationally feasible even up to systems as large as XH4
and XH5 (as discussed in the accompanying paper II).
Though this approximation may at first seem either simplistic or Draconian, it
turns out to be neither. As this first paper attempts to illustrate, such a model
accurately captures the essential large amplitude bend/rotational dynamics of small
XHn systems specifically corresponding to n = 2,3. Indeed, by systematically varying
the interatomic H-H potential “stiffness” and internal rotor constant for a single H
atom, the effective molecular system can be tuned from a nearly “free internal rotor,”
where large amplitude dynamics correspond to fully delocalized wave functions, to amore “rigid molecule” limit with energy-level-splitting patterns characteristic ofrovibrational tunneling. It is worth stressing at the outset that spectroscopicallyaccurate predictions for rotational/tunneling level splittings are not the goal of thiswork, nor are they even realistic with any known potentials or current computational
methods. Rather my modest aim is threefold: (i) to develop a general computational

35
formalism for n-particle systems in reduced dimensionality, (ii) to establish numericalconfidence in these methods for simple (n = 2,3) systems by comparison with
spectroscopically well known triatomic and tetratomic test case molecules, and (iii) to
exploit these methods to develop semi-quantitative intuition for internal
rotor/tunneling dynamics in larger systems (i.e., n = 4,5), specifically for the lowest
few J levels (J = 0,1,2). These low J predictions should complement and provide
particularly valuable guidance in the acquisition and analysis of jet-cooled spectra,
currently obtained under supersonic expansion conditions in a high-resolution IR
laser slit discharge spectrometer.55
The organization of this chapter is as follows: Section II describes the general
theoretical background necessary for efficiently solving the POS problem, focusing
on methods used in this work for the two- and three-particle system. Sections III and
IV discuss application of the POS model to two test systems, n = 2 (XH2) and n = 3(XH3), respectively. In Section V, I compare these n = 2,3 POS predictions withexperimental results for large amplitude systems such as H2O, NH3, and H3O+. For n= 2,3, I compare these results with predictions from full dimensional theoretical
calculations, which confirm the remarkable quantitative accuracy of the POS model.
Concluding comments are briefly summarized in Section VI.

II. Theoretical BackgroundAs mentioned above, the model restricts the total dimensionality by
constraining particles radially, while permitting them to execute large amplitude
angular motion on the surface of a sphere. The resulting Hamiltonian for n particles
can be simply expressed as

36n

H POS = ∑ bi ˆji2 + Vˆ (Ω ) ,

(2.1)

i =1

where ĵi is the angular momentum of the ith hydrogen with respect to the central
stationary atom and bi is the rotational constant for motion on the sphere. Vˆ (Ω ) is thepotential describing the H atom interactions, where Ω is a 2n dimensional vector of
all H atom angular coordinates. Simply stated, this Hamiltonian is expanded in asuitable basis set, and the matrix elements of HPOS are calculated, with the resulting
matrix diagonalized to generate the desired eigenvalues and eigenfunctions. Details
that make such a model practicable to implement, particularly in larger atom systems,
are described below.

A. Primitive and coupled basis sets
The choice of a primitive basis set is dictated by the large amplitude nature of
the problems I wish to address. Freezing the n X-H stretch coordinates yields a much
simpler 2n dimensional Hamiltonian, but one which nevertheless permits arbitrarilystrong coupling between overall rotation and H-X-H bending vibration. In thecompletely “floppy” limit of vanishing interaction between H atoms (V(Ω) = 0), thequantum mechanics is simply that of n free particles constrained to a sphere. In thecorresponding “stiff” extreme of strong interatomic interactions, one anticipates thequantum behavior to be described by localized harmonic oscillator wave functions ineach of the H-X-H angular coordinates. Since the goal of this work is to provide
insight into low-frequency rovibrational modes of molecules in the floppy limit (suchas CH5+), I choose a primitive basis diagonal in the kinetic energy operator and treat

37
the potential perturbatively. Such a basis is the direct product of rigid rotor basis
functions:
n

∏ jm
i

i

,

(2.2)

i =1

where ji is the angular momentum of the ith hydrogen and mi is the projection of ji on
the z axis.
The n-particle Hamiltonian is diagonal in total angular momentum, J, so it isadvantageous to transform this primitive basis into the coupled representation,56 Thispermits matrix construction and diagonalization for each quantum number J, withcorresponding reduction in computational expense. The coupled basis for the XH2
system is denoted as( j1 j 2 )JM(2.3)

and for the XH3 system as(( j1 j 2 ) j12 j3 )JM

,

(2.4)

where ji is the angular momentum of the ith particle, jij is the vector sum of ji and jj,
and the parentheses denote which angular momenta are coupled. J is the total angularmomentum, and M is the projection of J on the z axis, with standard sums over 3-J
symbols56 to transform between coupled and uncoupled representations. Whencreating a basis set for a given J, the primitive free particle angular momenta are
restricted to ji = [0, jmax], where jmax is chosen to achieve maximum convergence
while minimizing computational time. In the absence of external fields, all M states
for a given J are degenerate; thus M can be taken to be 0 for an additional 2J + 1
reduction of the Hamiltonian matrix.

38

B. Basis set symmetrization and permutation inversion
For a crucial further reduction in computational effort, the Hamiltonian isblock diagonalized using group theoretical methods.57, 58 Since these systems
exemplify large amplitude motion, I use permutation inversion (PI) theory58, 59 tocreate symmetry-adapted linear combinations (SALCs) of the coupled basis functions
that transform according to each PI symmetry group. The PI group for XH2 is G4(isomorphic to C2v), which yields four irreducible representations (A1, A2, B1, andB2), all nondegenerate. For XH3, on the other hand, the permutation inversion groupis G12 (isomorphic to D3h), characterized by six irreducible representations (A1’, A2’,
E’, A1’’, A2’’, and E’’), two of which (E’, E”) are intrinsically twofold degenerate. As
clearly outlined by Bunker and Jensen,59 standard quantum projection operator
techniques can be used to systematically generate SALCs of the coupled basis
functions represented by Equations (2.3) and (2.4) that transform according to asingle irreducible representation. This permits block diagonalization of the
Hamiltonian and a four- and sixfold reduction in the basis set size, respectively, forcalculation of matrix elements. Furthermore, for systems with three or more identicalparticles and therefore intrinsic degeneracies, one need only calculate these
eigenvalues once. Thus, by symmetrizing the basis set such that only one basis
function for each degenerate level is included,60, 61 the matrix size for a given
symmetry can be further reduced. Though not a limiting issue for n = 2 and 3
systems, this size reduction becomes important for larger n. For example, in CH5+ thestate degeneracy can be as high as six, providing an additional 36-fold reduction in
matrix element evaluations for a given state symmetry.

39
The use of permutation inversion-symmetrized basis sets also automaticallypermits nuclear spin statistics for the identical particles to be implemented.
Specifically for XHn, one constructs a reducible representation of the n identicalhydrogen atoms (spin ½ fermions), which is reduced according to standard grouptheoretical procedures59 to obtain nuclear spin weights. For n = 2 and n = 3, theseweights are readily found to be 3:1:3:1 (for Γ = A1:A2:B1:B2) and 2:2:1:2:2:1 (for Γ =A1’:A2’:E’:A1”:A2’’: E”) respectively. As the energy eigenvalues are automaticallysorted by symmetry, this permits immediate assignment of the nuclear spin states andweights.

C. Matrix element evaluation
The use of symmetrized basis functions drastically decreases the matrix size
for large n, thereby greatly enhancing the calculational efficiency. However, there is
still significant computational cost in evaluation of the resulting matrix elements,
which scales rapidly with the complexity of the system. Indeed, to more clearly
illustrate the convergence issues and rapid scaling with molecular size, Table 2.1
displays typical results for a given n-atom system, total angular momentum (J) andjmax internal rotor value in the uncoupled representation. Also shown in Table 2.1 fora series of n, J and jmax values are (i) Nsym, the number of irreducible representations
for a given n, (ii) Ncoupled, the number of coupled free rotor basis states consistent witha given total J and all ji ≤ jmax, (iii) NSALC, the average number of linearly independent
SALCs per symmetry species, (iv) MSALC, the average number of coupled basis statesin a given SALC, and (v) Nintegral and Ntot, the total number of integral evaluations
required to generate the matrix per atom pair or for all pairs, respectively.

40
The resulting matrix size (NSALC) is decreased by as much as 50-fold from thefull-coupled basis set (Ncoupled); thus diagonalization and storage are not the limiting
issues. Nevertheless, the number of integral evaluations to fill this matrix increases as
Nintegral ≈ 0.5*(NSALC )2 x (MSALC)2 per pair of atoms, which yields Ntot ≈ 0.25*(NSALC)2 x (MSALC)2 x n2 for the full sum over all atom pairs. As shown in Table 2.1, thenumber of integral evaluations for the two- and three-particle systems considered inthis work is quite manageable, even up to the large jmax values required for numericalconvergence with nearly rigid intermolecular anisotropies (e.g., models for H2O andNH3). However, the scaling behavior is so steep that an n = 5 system (such as CH5+)
will require as many as 6 × 1010 integral evaluations for only J = 2 jmax = 3. To makelarger n calculations feasible with this model, however, these evaluations need to be
exceedingly efficient.
Since both the primitive and coupled basis sets are diagonal in ji2, the kinetic
energy operator matrix elements satisfy this requirement automatically. For XH2,
these matrix elements simply are:

( j1 j 2 )J Tˆ ( j1' j 2' )J ' = δ JJ 'δ j j δ j j b[ j1 ( j1 + 1) + j 2 ( j 2 + 1)] ,
'
1 1

'
2 2

(2.5)

with obvious extension to XH3 (and higher particle) systems:(( j1 j2 ) j12 j3 )J Tˆ (( j1' j2' ) j12' j3' )J ' =δ JJ 'δ j j δ j j δ j j δ j j b[ j1 ( j1 + 1) + j2 ( j2 + 1) + j3 ( j3 + 1)].
'
1 1

'
2 2

'
12 12

(2.6)

'
3 3

The potential matrix elements, however, require substantially greater effort. Matrix
element evaluation by standard quadrature methods increases times by an additional
factor of Q2n, where Q is the number of points on a potential grid, i.e., unacceptable
for n > 3. By considering the potential as a multibody expansion, i.e.,

41
n

n

n

i< j

i< j <ki < j < k <lV (Ω) = ∑ V ( 2) (Ω i , Ω j ) + ∑ V ( 3) (Ω i , Ω j , Ω k ) + ∑ V ( 4 ) (Ω i , Ω j , Ω k , Ω l ) + ...

,

(2.7)
where each term represents a sum over all two-body, three-body, and four-bodycontributions, the work required to evaluate the potential elements can be greatly
reduced. Given the rapidly growing number of evaluations required in Table 2.1,
however, I restrict this expansion to only the first term, i.e., I consider arbitrary linear
combinations of pairwise additive intermolecular potentials between each particle.
Such pairwise additive potential depends only on the angle between pairs of particles:
V (Ω) = ∑ V ( 2 ) (γ ij )
n

,

(2.8)

i< j

where γij is the H-X-H angle between the ith and jth hydrogen.
Although not essential for the n = 2,3 test systems described in the current
study, such a pairwise additive approximation provides a crucial advance for treatinglarger systems, making theoretical studies of XH4 and XH5 tractable. Specifically, if I
expand any pairwise additive potential in a sum of Legendre polynomials, the matrix
element for each term can be calculated analytically via Clebsh-Gordan angularmomentum algebra.62 For the XH2 problem, this yields:( j1 j 2 )J Vˆ ( j1' j 2' )J ' = δ JJ ' ∑ f λ ( j1 j 2 ; j1' j 2' ; J ) ,

(2.9)

λ

where fλ are the Percival-Seaton coefficients:(f λ j1 j 2 ; j1' j 2' ; J) = ( j j )J P (cos γ ) ( j j )J[(2 j + 1)(2 j + 1)(2 j + 1)(2 j + 1)]= (− 1)
λ

1 2

' '
1 2

12

'

j2 + j2' + J⎛j× ⎜⎜ 1⎝0

λ0

1j ⎞⎛ j 2⎟⎜
0 ⎟⎠⎜⎝ 0
'
1

λ

0

2

j ⎞⎧ j 2⎟⎨
0 ⎟⎠⎩ j1'
'
2

'
1

λ
J

'
2

j ⎫
⎬
j1 ⎭
'
2

1
2

(2.10)

42⎛ . . .⎞ ⎧. . .⎫56⎟⎟ , ⎨and ⎜⎜⎬ are 3-J, 6-J symbols, respectively. For the n = 2 particle⎝ . . .⎠ ⎩. . .⎭system, calculation of these potential matrix elements is very fast, exploiting efficient
3-J and 6-J evaluation routines.63
The extension to a larger number of particles is straightforward. The n = 3
potential is expressed as a sum of pairwise interactions, i.e.,
V (Ω ) = V ( 2 ) (γ 12 ) + V ( 2 ) (γ 13 ) + V ( 2) (γ 23 ) ,(2.11)

where the matrix element of the first term is analogous to the two particle matrix
element and easily calculated as:(( j1 j 2 ) j12 j3 )J Vˆ ( 2) (γ 12 ) (( j1' j 2' ) j12' j3' )J ' =(δ JJ 'δ j j ' δ j j ' ∑ f λ j1 j 2 ; j1' j 2' ; j12
3 3

12 12

λ

)

,

(2.12)The other two terms simply require transformation to a different coupling scheme.56
For example, the V(2)(γ13) term is calculated by recoupling the j1 and j3 angularmomenta first to produce j13, i.e.,(( j1 j2 ) j12 j3 )J Vˆ ( 2) (γ 13 ) (( j1' j2' ) j12' j3' )J ' =∑ C j13 C j13' (( j1 j3 ) j13 j2 )J Vˆ ( 2) (γ 13 ) (( j1' j3' ) j13' j2' )J ' ,(2.13)'
j13 , j13

where the Ci coefficients can be expressed as 6-J symbols and each term evaluated asin Eq. 2.12.
It may seem restrictive to limit the potential to a pairwise additive form;
however, it is reasonable to believe that the approximation will work well for the
floppy CH5+. The large zero-point motion of the bending/rotation degrees of freedom
and the low barriers to interconversion of the hydrogens in CH5+ are likely the43
predominant factors influencing the spectroscopic patterns.

The high-quality

potential43 of Bowman and co-workers contains a Cs minimum and two low-lying
transitions states: A Cs transition state at 42 cm-1 and a C2v transition state at 190 cm1. Preliminary calculations on the Bowman surface with a restricted RCH coordinate

result in a C2v minimum and one Cs transition state at 115cm-1. A fit of a pairwise
additive potential to the Bowman surface yields a similar C2v minimum and a similar
Cs transition state 300cm-1 higher.

All of the transition states are an order ofmagnitude smaller than the zero-point energy of about 11,000cm-1 (or 4,100cm-1
when freezing the C-H stretch coordinates). The effects of the pairwise potential isinvestigated in more detail in Chapter 3; however, the large-amplitude coupledbending/rotation patterns should not change significantly. As will be demonstrated
later in this chapter, a pairwise additive potential reproduces the full dimensional andbending levels of both NH3 and H3O+ remarkably well.
The pairwise additive potential only depends on the angle between each pairof particles, which can be expanded in Legendre polynomials with no loss ofgenerality. In the interest of simplicity, I take this form of the potential for this paperto be
V ( 2 ) (γ ij ) =[]k
1 − cos(γ ij − γ 0 ) ,2(2.14)

which exhibits a minimum at γ = γ0. The three free parameters in the full Hamiltonian
are: (i) b, the rotational constant for a single H atom, (ii) k, the depth of the potential
well, and (iii) γ0, the equilibrium H-X-H angle. I further introduce the parameter χ =k/b, the ratio of the potential well to the rotational constant, which provides a44
convenient dimensionless metric for characterizing the “floppiness” vs. “stiffness” ofthe pairwise interaction.
I conclude this section with an explicit timing example. By exploiting efficient
3-J/6-J evaluations, Equations (2.11) and (2.12), and logical extensions thereof, eachpotential energy matrix element for a given coupled basis state can be calculated on asingle 2.4 GHz Intel processor in 6.5 μs. This rate translates into ∼2 hours per 109
integral evaluations, or, by inspection of Table 2.1, ∼100 hours for a system with n =5, J = 2, jmax = 3. This is encouraging and suggests that such POS calculations may be
tractable even past n = 5, at least for systems with intermolecular potentials
sufficiently “floppy” to be adequately represented by small values of jmax. Symmetry
and nuclear-spin-sorted results will be presented in Chapter 3 for n = 4,5 atom
systems, benchmarked against CH4, and will offer an approximate but insightfulmodel of quantum level structure for the highly fluxional CH5+ species. To firstestablish the feasibility of these methods, however, I focus herein on explicit
examples of the POS model for the computationally straightforward n = 2,3 systems.

III. Two Particles-On-A-Sphere (XH2)
The first test system is that of two POS, i.e., XH2, with Hamiltonianparameters (bi = 18.00 cm-1, γ0 = 104o) chosen to resemble the equilibrium structureof H2O in the limit of large anisotropy in the intermolecular pair potential. The
strategy is to investigate the large amplitude eigenenergies as a function of thedimensionless stiffness parameter, χ = ki / bi, for fixed bi. With only two particles, one
can readily anticipate the limiting behaviors in each of two regimes. Specifically the χ<< 0 (“floppy”) limit should yield rotation-like energy level patterns characteristic of45
two free internal rotors, with characteristic spacings as a function of j1, j2. Conversely,
the energy levels in the χ >> 1 (“stiff”) limit should begin to approach the asymmetric
top rotor levels for a harmonic oscillator in the H-X-H bend coordinate.
Figure 2.2 shows a sample energy eigenvalue correlation diagram for two
POS (J = 0, A1 symmetry), as the stiffness of the potential is increased. When χ issmall, the progression scales like ji(ji + 1) for each of the two independent rigid
rotors. For the specific case of J = 0, however, the angular momenta of the twohydrogen atoms are perfectly correlated in the coupled basis, i.e., j1 = j2. Thus, in thelimit of low potential anisotropy, the states can be approximately labeled by |j1j2>sym
for a given symmetry and one thus observes energies in Figure 2.2 for χ ≈1 increasing
uniformly as Erot = 2bj(j + 1), for j = 0,1,2… In the stiff limit, the intermolecularpotential is large enough that the two atomic displacements are now correlated andbest described as a small amplitude vibration. The appropriate notation in this limit istherefore v sym , where the ket denotes the number of quanta and symmetry in thebending vibrational mode. As expected, the correlation diagram in Figure 2.2 for χ≈1000 (J = 0, A1 symmetry) demonstrates harmonic behavior, i.e., approximately

equal vibrational spacings in the first few energy levels.
If one examines multiple J values for the lowest vibrational state (v = 0) as afunction of intermolecular stiffness, an analogous but quite different trend isanticipated. By way of example, Figure 2.3 shows the energies (relative to EJ = 0) of J
= 0, 1, and 2 for XH2 as χ is increased from 1 to 10000. At low χ, the energy spacing
is 2b, reflecting uncorrelated rotation of the two hydrogen atoms. As χ is increased,
however, evolution to a more localized rotational structure of XH2 is clearly evident.

46
Specifically, for χ ≈ 10000, the energy levels for XH2 converge on those of a rigid
molecule, with J = 0, 1 and 2 states taking on the characteristic Ka, Kc patterns of anasymmetric top. This is most readily evident in the set of JKaKc = 101, 111, and 110
levels, which systematically converge (for high χ) on the characteristic asymmetric
top values of Erot = B + C, A + C, A + B, respectively.
Additional confirmation for the two-particle dynamics can be obtained fromthe corresponding wave functions, representative 2D slices for the J = 0 XH2 groundstate which are displayed in Figure 2.4 as a function of χ. To visualize the anisotropyinduced correlation between the hydrogen atoms, I rotate the coordinate system sothat one hydrogen is always pointing upwards at θ = 0 and φ = 0. For an n = 2 system,
the probability of finding the second particle is then plotted as a function of γ, theangle between the two hydrogen atoms. Figure 2.4 clearly illustrates the expectedbehavior for the XH2 system; the wave functions become more strongly localized as
the potential is stiffened, slowly converging on the equilibrium angular separation (γ0)
between the two hydrogen atoms as χ increases.
As a final test of the n = 2 model, the XH2 system can be mapped onto thewell-studied H2O molecule, for which the full 3D intermolecular potential has been
determined to exceptionally high accuracy by Tennyson and co-workers.64 In thespirit of the POS model, I first fix the two OH bond lengths at equilibrium values(ROH = 0.958 Å) and then fit the reduced dimensionality-bending H2O potential to a
15-term Legendre expansion. The resulting 1D potential is then inserted into the
Hamiltonian with b = 18.4 cm-1 (obtained from equilibrium bond lengths) anddiagonalized to yield the converged rovibrational eigenstates and eigenfunctions for47
H2O in the POS limit. To facilitate direct quantitative comparison, the equivalent
asymmetric top rotational constants are extracted from the three lowest J = 1 levels,
as indicated in Figure 2.3. Similarly, the HOH vibrational bending levels are
estimated via differences in the lowest three J = 0 states of A1 symmetry, as illustrated
in Figure 2.2. The resulting rovibrational predictions are compared with known
experimental values for H2O in Table 2.2, which indicate nearly quantitative
agreement (< 0.1% for vibrations, < 1 - 5 % for rotational constants). Given thecalculational simplicity of this POS model for n = 2, this level of agreement for H2O
appears promising, especially considering the complete absence in such a model ofany stretch-bend coupling.

IV. Three Particles-On-A-Sphere(XH3)
With these encouraging results for the two-particle Hamiltonian system, I add
one more atom to the model to form XH3. In anticipation of numerical comparison
with experimental data, the intermolecular structural parameters for this system havebeen chosen to best approximate NH3, with b = 16.63 cm-1 and γ0 = 1080. Since thismodel potential now predicts an equilibrium nonplanar C3v structure with a double
minimum along the inversion coordinate, one expects to see evidence of (i)symmetric top rotational behavior and (ii) quantum mechanical tunneling through theinversion barrier.
Analogous to the two-particle system, I consider appropriate state labels as afunction of low and high intermolecular anisotropy. In the floppy limit (χ ≈ 1), thestates can be well designated by

j1 j 2 j3 sym where the symmetry label uniquely

48
identifies the proper linear combination of the three-particle permutations. In thesemirigid molecular limit for XH3 , there are three vibrational modes: one is theumbrella mode of A1’/A2” symmetry (νumb) with the remaining pair forming adegenerate bend (νbend) of E’/E” symmetry. The most convenient notation in the stiff
limit is therefore vumb vbend sym , with the corresponding number of quanta in either the
A umbrella or the E bend vibration explicitly denoted.
As in the case of XH2, I first examine the correlation diagram for energy
eigenvalues between floppy and stiff potential limits. Figure 2.5 displays
eigenenergies of the first few A1’/A2” symmetry states of XH3 with J = 0. In the
floppy limit (χ ≈ 1), a simple progression of rigid rotor energy levels for XH3 can beeasily recognized, although more complicated than for XH2. This is because there are
multiple ways to add the three H atom angular momenta (j1,j2,j3) and still get total J =0. Nevertheless, the lowest several states of A1’/A2” symmetry (with respect to EJ = 0)
are found to be consistent with the free rotor limit of Erot = 4b, 6b, 10b, 12b, etc. The
corresponding lowest energy level behavior in the stiff limit (χ >>1) clearly correlateswith a simple progression in the umbrella mode, split by tunneling through the planar
barrier. (The XH3 bending modes are present in this energy range as well but only
appear as doubly degenerate vibrations of E’/E” symmetries.) Note that themagnitude of these tunneling splittings increases rapidly with vumb, corresponding
semiclassically to vibrational promotion of barrier penetration as one energeticallyaccesses the planar transition state.
By way of confirmation, one can similarly examine low J levels in the groundvibrational v = 0 state as a function of χ . Figure 2.6 displays the correlation diagram

49
of rotational energies (relative to EJ = 0) for J = 0, 1, and 2 states of XH3 as a functionof increasing χ. At χ ≈ 1, the free rotor pattern emerges, with characteristic 2b energy
spacings. In the limit of χ >> 1, however, the correlation diagram correctly predicts aconvergence toward the familiar energy levels of a rigid symmetric top, specifically,EJK = B J (J + 1) + (C - B) K2. (Note that this requires inclusion of E’ symmetry torepresent the intrinsically degenerate K ≠ 0 states.) Since the equilibrium structuralparameters of the potential are modeled after NH3, the energy levels are those of anoblate symmetric top, i.e., energy decreasing with K for a fixed J value.
Further examination of the model behavior can be obtained from 2D slices
through the corresponding eigenfunctions. As the total wave function, Ψtot, is notpossible to plot in 6D, I examine a simpler 2D function of one hydrogen coordinate(θ, φ), with the other two hydrogens held fixed at (θ = 0, φ = 0) and (θ = γ0 φ = 0),
respectively (see insert in Figure 2.7 for relevant details). A series of slices throughΨtot is shown in Figure 2.7 as a function of χ, with the lowest A1’ and first excited

A2” states for J = 0 on the left and right, with the equivalent slice through thepotential contour shown at the bottom. As expected for an increase in potential
stiffness (i.e., increasing χ), the wave function probability for the last H atom
progressively localizes on the two equilibrium positions (γ = γ0 = 1080, φ = 1200,
2400). Furthermore, the wave function amplitudes clearly exhibit the expected
alternation between symmetric and antisymmetric wave function characteristic of adouble well potential.
As a logical benchmark test of these n = 3 POS calculations, I look to thelimiting case of ammonia (NH3), for which the intermolecular potential and energy

50
term values are well known from previous spectroscopic studies. Currently, the model
is only capable of representing a pairwise additive potential, whereas the mostaccurately determined NH3 potential65 requires small three-body terms to fully
describe the interaction of the three hydrogens. As I will be unable to include suchnonpairwise additive terms for n > 3 (e.g., CH5+), I choose instead to fit an effective
pairwise potential that most accurately reproduces the true inversion potential forNH3. Specifically, the k and γ0 parameters in Eq. 2.6 are fit to reproduce (i) theexperimentally known barrier height at planarity and (ii) the equilibrium HNH bend
angle for NH3. Note that this does not represent a fit to the complete tunneling path,as was done in the previous section for H2O. However, this approach does offer
additional insight into the importance of such three-body terms, for which sufficientcomputational resources will not be available for larger particle systems.
The results of these test calculations for NH3 are summarized in Table 2.3 for
J = 0 vibrational levels up to 2000 cm-1, along with the experimental values. Both thevibrational energies and inversion tunneling splittings for the umbrella mode match
experiment remarkably well (< 3–5%), particularly for such a simple two-parameter
representation of the ammonia potential. Though still respectable, the agreement is
not quite as good (indeed, systematically underestimated by ∼9%) in the degeneratebending vibrational predictions. However, this is not entirely surprising, since only
barrier height information has been included in formulating the approximate NH3
surface.

51

V. Discussion
I first begin with a brief discussion of the convergence properties for mycalculations. In Figure 2.8, the fractional error for the ground-state energy is plotted
as a function of jmax for several values of χ. As expected, the error (relative to fully
converged results at jmax = 50) decreases with increasing jmax and basis set size.
Furthermore, the size of jmax necessary to converge to a given level increases withstiffness in the potential. This is simply the result of a basis set designed to bediagonal at V(Ω) = 0; i.e., as the interatomic potential is increased, the free rotor
picture provides a progressively worse zeroth-order description. Also shown in Figure
2.8 for comparison are corresponding convergence tests for the XH3 system.
Particularly noteworthy is that the convergence rate for a given value of χ isremarkably similar for both XH2 and XH3. This is extremely encouraging forextension of these studies to n = 4,5 particle systems, for which the basis set size (andhence the Hamiltonian) grows so quickly that calculations for asymptotic values ofjmax may not be feasible to establish a rigorous level of convergence for a given test
potential. For the n = 2,3 systems discussed herein, however, all the calculations havebeen performed with jmax > 20, i.e., sufficiently high to ensure negligible fractional
errors (< 0.0001 %).
Now that I have demonstrated the capability of the POS model for reasonably
accurate predictions in “stiff” XH2 and XH3 systems such as H2O and NH3, I next
turn my attention to the floppier molecular ion, H3O+. The hydronium ion has been
extensively studied over the past several decades with both high level theoretical andexperimental methods. As a result, the intermolecular potential for this ion is known

52
to nearly spectroscopic accuracy. The barrier height in H3O+ (∼ 650 cm-1) is
considerably lower than for NH3 (∼ 1792 cm-1), which leads to tunneling splittings in
the ground state of nearly 60 cm-1. A test study of H3O+ via POS methods is therefore
helpful in developing further confidence in calculations for other large amplitude
tunneling systems, specifically other closed-shell-protonated molecular ions such as
CH5+.
I start by fitting the theoretical H3O+ 1D symmetric inversion potential66 to a
15-term Legendre expansion, again without including any three-body terms. Fully
converged rovibrational energy levels for the resulting potential are then obtained
from the POS Hamiltonian treatment, with b = 17.6 cm-1. The results are listed inTable 2.4. Similar to what was observed for the NH3 test problem, the predicted H3O+
energies agree reasonably well with both experiment and exact full 6D calculations
by Bowman et al. Indeed, the typical error between (i) the POS and (ii) full
dimensional calculations is on the order of 3–5%, which is certainly adequate for sucha simple model based on (i) an infinitely heavy central atom and (ii) fixed moment ofinertia for each HX rotor. Also similar to what was observed in the NH3 studies, there
are noticeably larger errors (8–12%) for predictions in the degenerate-bending
vibrational manifold. This may be simply attributable to the fact that the POS
potential is only fit to (and thus might be expected to better reproduce) the H3O+
potential along the inversion coordinate. Nevertheless, the results are reasonably
accurate for such a fit and may indeed be further improved by a better modeling of
the potential shape. Most importantly, however, the overall energy level patterns forthe umbrella inversion, tunneling, and degenerate bending modes are accurately and53semiquantitatively reproduced. Since one aim of the POS approach has been todevelop qualitative intuitions in highly fluxional tunneling systems, this may alreadyprove of substantial value in more challenging molecular ions such as CH5+.

VI. Summary and Conclusion
In this chapter, a conceptually simple method has been introduced forapproximate but reasonably quantitative treatment of large amplitude quantum
motion in small floppy hydrides such as XHn. This POS method takes advantage ofthe approximate independence of XH bond length on angular coordinates to motivate
a reduced dimensionality Hamiltonian treatment of hydrogenic motion constrained tothe surface of a sphere. Despite the simple nature of this approximation, the method isfound to work surprisingly well for test n = 2 and n = 3 particle systems. Furthermore,
the method yields good quantitative agreement in tests against experimental results
for stiff systems with known intermolecular potentials (e.g., H2O and NH3) as well as
more “floppy,” large amplitude systems (e.g., hydronium ion, H3O+). The agreementwith full dimensional theoretical calculations is extremely promising, particularly inthe quantitative reproduction of all patterns relevant to spectroscopic assignment andanalysis of highly quantum mechanical systems.
Most importantly, however, the numerical methods developed with thisapproach lead to favorable scaling properties that offer encouragement towardapplications with even larger n atom systems. The test cases shown hereindemonstrate the feasibility of extension of these methods up to n = 5, at least forangular momentum states (J = 0, 1 and 2) relevant to spectroscopic studies in a lowtemperature slit supersonic expansion. As all n = 2,3 test molecules studied with POS

54
methods in this paper reflect relatively stiff intermolecular potentials, I feel confident
that significant progress can be made toward larger systems. In later chapters, I apply
these POS methods to n = 4,5 particle systems along with any extra optimizations
needed to converge these exponentially much harder problems. The aim is to provide
reduced dimensionality tools to develop the necessary insight and intuition for large
amplitude dynamical systems, for which calculations in full dimensionality aresimply not feasible. The long range goal is to incorporate these reduced
dimensionality approaches with high-level ab initio potential surfaces to provide
valuable information for decoding high-resolution spectra of highly fluxional species
such as CH5+.

55

TablesTable 2.1. Basis set size and matrix calculation time scaling as a function of n, J, and jmax.
n
2

J
0
1
2

3

0
1
2

4

0
1
2

5

0
1
2

jmax
50

15

6

3

Nsyma

Ncoupled
51
150
246

NSALC

MSALC
1
2
2

Nintegral
1,326
5,050
7,750

Ntot

51
50
62

1,326
5,050
7,750

Time
0.01 s0.03 s0.05 s

6

2,056
6,120
10,050230
680
1.1 x103510
16

0.7 x106
23 x106
159 x106

2 x10669 x106
478 x1060.22 m7.54 m51.97m10

5,719
16,68826,320244
690
1.1 x103

28
75
113

23 x106
1,339 x106
7,753 x106

140 x1068,034 x106
46,520 x106

0.30 h
14.6 h84.3 h14

4,269
11,925
17,225

67
186
268

101
276
400

0.02 x109
1.3 x109
5.7 x109

0.2 x109
13 x10957 x1090.40 h
23.90 h104.1 h4aNsym refers to the number of irreducible symmetry representations; Ncoupled is the number of coupled basis set states for agiven J, jmax; NSALC is the average number of symmetry-adapted linear combinations (SALC) per symmetry species; MSALC
is the average number of coupled basis states in a given SALC; Nintegral is the number of integral evaluations per atom pair;
Ntot is the total number of integral evaluations over all atom pairs.

56
Table 2.2. Calculated and experimental H2O energy levels (cm-1) for the bending
mode and rotational constants (cm-1).

POS

a,b

3-D QMa

Experimentb

1 A1

1588

1594.71594.72 A13149

3151.53151.6

A
B
C

25.67
14.89
9.19

Values from Polansky et al.64

24.27
14.81
9.20

57
Table 2.3. Calculated and experimental NH3 energies and splittings.

State

POS
EnergySplitting

00 A '

0

10 A ''

0.798

20 A '

895.6

30 A ''

930.9

01 E '

1482

11 E ''

1483

40 A '

1536

50 A ''

1813

60 A '

1

2

1

2

Full QMa
Energy Splitting0} 0.7980.790
934.06} 35.3}35.51

00.793
932.43968.12

1627.91

1626.28

1600.79

1597.47}0.793}35.691885.06}284.27 1882.18 }284.71

2299

2387.82

2384.15

21 E '

2413

2543.43

2540.53

70 A ''

2795

2899.622895.51

02 A '

2938

3220.133216.102 E '

2953

3244.26

3240.4480 A '

3345

3467.24

3462

1

2

1

2

1

1

a,b

} 277

969.57}0.790Experimentb
Energy Splitting

Values from Rajamaki et al.67

58
Table 2.4. Calculated and experimental H3O+ energies and splittings.

State

POS
Energy Splitting

Full QMa
Energy Splitting

Experimentb
Energy Splitting

00 A '

0

10 A ''

54.33

20 A '

529.9

30 A ''

876.7

40 A '

1377

1445

1475.84

01 E '

1421

1624

1625.95

11 E ''

1498

1681

1693.87

50 A ''

1918

2005

1882.18

1

2

1

2

1

2

a,b

}54.33}347

Values from Huang et al.66

0
46
580
934

}46
}354

0
55.35
581.17954.4}55.35}373.23

59

FiguresFigure 2.1. Schematic picture of the particles-on-a-sphere (POS) method, which isbased on a modeling of polyatomic hydrides (XHn) as one heavy atom (X) that
constrains the motion of n H atoms to the surface of a sphere.

60

Fig 2.2. EJ = 0 energy level correlations for XH2, plotted for the first few lowest
energies vs. potential stiffness, χ. At small χ, a simple rigid rotor pattern withincreasing spacings (4b, 8b, 12b, etc.) is clearly exhibited (the rotational levels are
reproduced with greater detail in the inset), but which evolves smoothly to a harmonic
pattern of constant energy spacing (hω0) in the limit of large χ.

61

Fig 2.3. EJ = 0,1,2 energy level correlations for XH2 in the ground vibrational state,plotted for the first few lowest energies versus potential stiffness, χ. At small χ, theenergy separations reflect uncorrelated free rotation of the hydrogen atoms. At large
χ, however, this smoothly transitions into the characteristic JKaKc patterns of a rigid
asymmetric top.

62

Figure 2.4. 1D wave function probability distribution for the XH2 ground vibrational
state plotted as a function of increasing χ. As the potential stiffness is increased, thewidth of the HXH angular distribution narrows toward γ = γ0, i.e., the equilibriumgeometry. As shown in the inset, the first hydrogen is defined to be at (θ = 0, φ = 0),
with the |amplitude|2 for the second hydrogen plotted for (θ = γ, φ ).

63

Figure 2.5. EJ = 0 energy level correlation as a function of intermolecular stiffness for
XH3, plotted for the first few energies of the umbrella mode. At small χ, the levels
correspond to three independent rigid HX rotors (displayed with greater detail in theinset) but eventually correlate with inversion doubled tunneling splittings in the limitof large χ.

64Figure 2.6. Ev = 0 energy level correlation for XH3 J = 0, 1, 2 as function ofintermolecular stiffness, χ. For small χ, the pattern reflects three freely rotatinghydrogens, which for large χ smoothly transitions into the rotational energy patterncharacteristic of an oblate symmetric top.

65

Figure 2.7. 2D cuts (in θ, φ space) through full 6D XH3 wave functions for (i)symmetric ground state (left hand side) and (ii) the antisymmetric first excited state(right hand side) of the umbrella mode for increasing values of χ. The contour plot
(center bottom) displays the corresponding potential. At χ = 6, 60 and 600, thecontour spacings are 1.5 cm-1, 15 cm-1, and 150 cm-1, respectively. As shown in theinset, the first and second hydrogens are fixed at (θ = 0, φ = 0) and (θ = γ0, φ = 0),
with amplitude for the third hydrogen plotted as a function of (θ , φ ).

66

Figure 2.8. The log of the fractional error (base 10) as a function of jmax for several
values of stiffness χ. Larger χ converges more slowly as the basis set is optimized for
floppy intermolecular potentials. For a given χ, it is worth noting that the n = 2 and n= 3 calculations converge almost at the same rate, which bodes favorably for scaling
to a larger number of H atoms.

67

References
1

V. Spirko, J. Mol. Spectrosc. 101, 30 (1983).

2

S. Urban, R. Dcunha, K. N. Rao, and D. Papousek, Can. J. Phys. 62, 1775
(1984).

3

S. L. Coy and K. K. Lehmann, Spectroc. Acta Pt. A-Molec. Biomolec. Spectr.
45, 47 (1989).

4

V. Spirko and W. P. Kraemer, J. Mol. Spectrosc. 133, 331 (1989).

5

S. Urban, N. Tu, K. N. Rao, and G. Guelachvili, J. Mol. Spectrosc. 133, 312(1989).

6

I. Kleiner, G. Tarrago, and L. R. Brown, J. Mol. Spectrosc. 173, 120 (1995).

7

I. Kleiner, L. R. Brown, G. Tarrago, Q. L. Kou, N. Picque, G. Guelachvili, V.
Dana, and J. Y. Mandin, J. Mol. Spectrosc. 193, 46 (1999).

8

C. Cottaz, I. Kleiner, G. Tarrago, L. R. Brown, J. S. Margolis, R. L. Poynter,
H. M. Pickett, T. Fouchet, P. Drossart, and E. Lellouch, J. Mol. Spectrosc.
203, 285 (2000).

9J. Tang and T. Oka, J. Mol. Spectrosc. 196, 120 (1999).

10

D. J. Liu, T. Oka, and T. J. Sears, J. Chem. Phys. 84, 1312 (1986).

11

P. B. Davies, S. A. Johnson, P. A. Hamilton, and T. J. Sears, Chem. Phys.
108, 335 (1986).

12

M. H. Begemann and R. J. Saykally, J. Chem. Phys. 82, 3570 (1985).

13

M. Gruebele, M. Polak, and R. J. Saykally, J. Chem. Phys. 87, 3347 (1987).

14

H. Petek, D. J. Nesbitt, J. C. Owrutsky, C. S. Gudeman, X. Yang, D. O.
Harris, C. B. Moore, and R. J. Saykally, J. Chem. Phys. 92, 3257 (1990).

68
15

M. Araki, H. Ozeki, and S. Saito, J. Chem. Phys. 109, 5707 (1998).

16

M. Araki, H. Ozeki, and S. Saito, Mol. Phys. 97, 177 (1999).

17

M. Farnik, S. Davis, and D. J. Nesbitt, Faraday Discuss. 118, 63 (2001).

18

M. J. Weida and D. J. Nesbitt, J. Chem. Phys. 110, 156 (1999).

19

D. T. Anderson, M. Schuder, and D. J. Nesbitt, Chem. Phys. 239, 253 (1998).

20

D. T. Anderson, S. Davis, and D. J. Nesbitt, J. Chem. Phys. 104, 6225 (1996).

21

D. H. Zhang, Q. Wu, J. Z. H. Zhang, M. Vondirke, and Z. Bacic, J. Chem.
Phys. 102, 2315 (1995).

22

G. W. M. Vissers, G. C. Groenenboom, and A. van der Avoird, J. Chem.
Phys. 119, 277 (2003).

23

Y. H. Qiu and Z. Bacic, J. Chem. Phys. 106, 2158 (1997).

24

Y. H. Qiu, J. Z. H. Zhang, and Z. Bacic, J. Chem. Phys. 108, 4804 (1998).

25

P. Eggert, A. Viel, and C. Leforestier, Comput. Phys. Commun. 128, 315(2000).

26

N. Goldman, R. S. Fellers, C. Leforestier, and R. J. Saykally, J. Phys. Chem.
A 105, 515 (2001).

27

F. N. Keutsch, N. Goldman, H. A. Harker, C. Leforestier, and R. J. Saykally,
Mol. Phys. 101, 3477 (2003).

28

N. Goldman, C. Leforestier, and R. J. Saykally, J. Phys. Chem. A 108, 787(2004).

29

J. S. Tse, D. D. Klug, and K. Laasonen, Phys. Rev. Lett. 74, 876 (1995).

30

P. R. Schreiner, S. J. Kim, H. F. Schaefer, and P. V. Schleyer, J. Chem. Phys.
99, 3716 (1993).

69
31

H. Muller and W. Kutzelnigg, Phys. Chem. Chem. Phys. 2, 2061 (2000).

32

D. Marx and M. Parrinello, Nature 375, 216 (1995).

33

G. A. Olah, N. Hartz, G. Rasul, and G. K. S. Prakash, J. Am. Chem. Soc. 117,
1336 (1995).

34

G. A. Olah and G. Rasul, Accounts Chem. Res. 30, 245 (1997).

35

E. Herbst, S. Green, P. Thaddeus, and W. Klemperer, Astrophys. J., 503(1977).

36

D. Talbi and R. P. Saxon, Astron. Astrophys. 261, 671 (1992).

37

A. Gamba, G. Morosi, and M. Simonetta, Chem. Phys. Lett. 3, 20 (1969).

38

G. A. Olah, G. Klopman, and R. H. Schlosberg, J. Am. Chem. Soc. 91, 3261
(1969).

39

W. T. A. M. Van der Lugt and P. Ros, Chem. Phys. Lett. 4, 389 (1969).

40

W. A. Lathan, W. J. Here, and J. A. People, J. Am. Chem. Soc. 93, 808(1971).

41

H. Muller, W. Kutzelnigg, J. Noga, and W. Klopper, J. Chem. Phys. 106,1863 (1997).

42

A. Brown, B. J. Braams, K. Christoffel, Z. Jin, and J. M. Bowman, J. Chem.
Phys. 119, 8790 (2003).

43

A. Brown, A. B. McCoy, B. J. Braams, Z. Jin, and J. M. Bowman, J. Chem.
Phys. 121, 4105 (2004).

44

V. L. Tal'rose and A. K. Lyubimova, Dokl. Akad. Nauk SSSR 86, 909 (1952).

45

E. T. White, J. Tang, and T. Oka, Science 284, 135 (1999).

46

P. R. Schreiner, Angew. Chem.-Int. Edit. 39, 3239 (2000).

7047

A. Komornicki and D. A. Dixon, J. Chem. Phys. 86, 5625 (1987).

48

W. Klopper and W. Kutzelnigg, J. Phys. Chem. 94, 5625 (1990).

49

S. J. Collins and P. J. Omalley, Chem. Phys. Lett. 228, 246 (1994).

50

M. Kolbuszewski and P. R. Bunker, J. Chem. Phys. 105, 3649 (1996).

51

A. L. L. East, M. Kolbuszewski, and P. R. Bunker, J. Phys. Chem. A 101,
6746 (1997).

52

P. R. Bunker, B. Ostojic, and S. Yurchenko, J. Mol. Struct. 695, 253 (2004).

53

G. A. Olah, Angew. Chem. Int. Ed. Engl. 34, 1393 (1995).

54

A. B. McCoy, B. J. Braams, A. Brown, X. C. Huang, Z. Jin, and J. M.
Bowman, J. Phys. Chem. A 108, 4991 (2004).

55

C. Savage, F. Dong, and D. J. Nesbitt, work in progress.

56

R. N. Zare, Angular Momentum: Understanding Spatial Aspects in Chemistry
and Physics. (John Wiley & Sons, New York, 1988).

57

F. A. Cotton, Chemical Applications of Group Theory, 3rd ed. (Wiley, New
York, 1990).

58

P. R. Bunker, Molecular Symmetry and Spectrosopy, 2nd ed. (NRC Research
Press, Ottawa, 1998).

59

P. R. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy, 2nd ed.
(NRC Research Press, Ottawa, 1998).

60

J. F. Cornwell, Group Theory in Physics: An Introduction. (Academic Press,
San Diego, 1987).

61

I. G. Kaplan, Symmetry of Many-electron Systems. (Academic Press, New
York, 1975).

71
62

J. M. Hutson, in Advances in Molecular Vibrations and Collision Dynamics,edited by J. M. Bowman (JAI Press, 1990), Vol. 1A.

63

K. Schulten and R. G. Gordon, J. Math. Phys. 16, 1961 (1975).

64

O. L. Polyansky, P. Jensen, and J. Tennyson, J. Chem. Phys. 105, 6490
(1996).

65

T. Rajamaki, A. Miani, and L. Halonen, J. Chem. Phys. 118, 10929 (2003).

66

X. C. Huang, S. Carter, and J. Bowman, J. Chem. Phys. 118, 5431 (2003).

67

T. Rajamaki, A. Miani, and L. Halonen, J. Chem. Phys. 118, 6358 (2003).

72

Chapter 3: Large amplitude quantum mechanics inpolyatomic hydrides: II. A particle-on-a-sphere modelfor XHn (n = 4,5)
I. Introduction
The quantum mechanics of large amplitude motion has been the focus ofexperimental and theoretical work for many years. Advances in high-resolution
rotational/vibrational spectroscopy have given much experimental evidence of large
amplitude dynamics in molecules and ions. The tunneling inversion in ammonia1-8(NH3) and hydronium ion9-16 (H3O+) are two such examples of well-studied systems
that exhibit large-amplitude motion. In systems with few enough degrees of freedom(N ≤ 6), the exact quantum nuclear dynamics can be solved in full dimensionality. For
example, this has been achieved in landmark studies of HF and HCl dimers17-20, as
well as reduced-dimensionality extensions to H2O dimer21-24 using six-dimensional(6D) intermolecular coordinates with semirigid H2O molecules.
However, as molecular systems increase in size and complexity, standard
quantum-calculation methods rapidly become too difficult to permit treatment in full
dimensionality. This is particularly unfortunate from a chemical physics perspective,
as some of the most interesting dynamical processes, such as the facile exchange ofidentical H atoms predicted to occur in protonated methane25,26, begin to emerge inthese larger systems. Furthermore, the level of theoretical difficulty increases rapidly
with total angular momentum, J, which makes the high-resolution interpretation of

73
rovibrational patterns especially challenging under all but the lowest temperature
conditions. In such cases, a more realistic first spectroscopic goal is to learn how to
predict and interpret the experimentally observed energy level patterns and thereby
infer new insights about the actual intermolecular dynamics.
As a specific target of interest, protonated methane (CH5+) represents an
extreme case of large-amplitude quantum dynamics with a rich scientific history. The
addition of a single proton to relatively “rigid” methane, CH4, creates the highly
fluxional CH5+ molecule where simple Lewis octet bonding motifs do not apply. Five
H atoms are connected by 4 electron pairs and require the presence of “three-centertwo-electron bonds” (“3c-2e”).27 Such 3c-2e bonding motifs correspond to a specialclass of hypercoordinated carbocations, which are extremely important reactive
intermediates in acid-catalyzed electrophilic reactions,28,29 for which CH5+ represents
the simplest prototypic “superacid.” CH5+ is also believed to be a key intermediate in
the synthesis of polyatomic organic species in cold interstellar clouds,30,31 further
motivating astrophysical interest in this simple, but spectroscopically elusive,
molecular ion.
From a theoretical perspective, CH5+ is interesting because its relatively small
size and high symmetry make quantum mechanical calculations of the potential
energy surface computationally tractable, though still quite challenging. The mostrecent high-level calculations suggest that there are, in fact, three low-lying energy
structures within about 1 kcal/mol of each other.25,32-34 Furthermore, if zero point
energy is taken into account, these three configurations are all extensively sampled by
the ground-state wave function. Stated simply, the barrier to rearrangement between

74
these low-lying minima is considerably lower than the quantum zero-point vibrationalenergy, thus permitting facile intramolecular scrambling of the hydrogen atoms32 totake place.
From a spectroscopic perspective, the highly fluxional, nonclassical nature of
CH5+ begins to account for long-standing difficulties in interpreting a high-resolution
spectrum. Despite its initial observation as a highly abundant ion in mass
spectrometers in the early 1950s,35 optical detection and characterization of CH5+
eluded spectroscopists for another 50 years. The breakthrough came in 1999 fromOka and coworkers, who obtained a spectrum36 in the CH-stretch region by exploiting
velocity modulation methods. They made convincing arguments that the extensive,
albeit rovibrationally unassigned, spectrum belonged predominantly to CH5+.
Recently, exciting progress has been made in obtaining broadband IR action spectraof CH5+ by Asvany and coworkers.37 These spectra yield unresolved vibrational
structure from 3500 down to < 1000 cm-1.

Our group has also obtained high-resolution jet-cooled IR spectra38 in the 2820 - 3050 cm-1 region.
Desire to help in this assignment process has naturally led to considerable
emphasis on direct calculation of near IR spectra from first principles. However,despite intense theoretical efforts directed toward CH5+25,32,39-51 and the recent
availability of a high-level potential surface,33,34 high resolution spectra based on fully
converged exact rovibrational energy levels have proven extremely challenging tocalculate. The reasons behind this challenge are at least twofold. On the experimentalside, large zero-point energies and lack of any substantial barriers between H atominterchange lead to extensive delocalization in the wave function. This in turn can

75
result in large tunneling splittings and rovibrational energy patterns profoundly
perturbed away from rigid rotor expectations, especially for the J states thermally
populated under discharge conditions. From a theoretical perspective, however, thenumber of degrees of freedom for CH5+ is already too large (15D) to achieve fully
converged quantum calculation of the rovibrational energy levels in full
dimensionality.
In Chapter 2,52 the particle-on-a-sphere (POS) model was introduced as analternative method for calculating the energies and wave functions for molecules withlarge amplitude motion in a large number of degrees of freedom. The model wastested on relatively “floppy” molecules (i.e., exhibiting large amplitude bending
dynamics) with considerable success for di- and tri-hydride species. In the XH2 case,the convergence of the energy levels was readily achieved for J ≤ 2, with therovibrational levels converged to near spectroscopic (< 10-4 cm-1) accuracy. By way
of example, Figure 3.1a displays rovibrational energies of the lowest J = 0 and J = 1
states for XH2 as a function of jmax. Here jmax is the maximum angular momentum fora single H basis state and therefore a measure of the variational basis set size, where
specific parameters for the POS potential and equilibrium bond lengths/bond angles
for XH2 are closely modeled after H2O. As expected, the variational energies decrease
with increasing basis set size, converging on the complete basis set limit to theground, fundamental, and 1st overtone excited bending states. Agreement at higher
resolution is also shown (see Figure 3.1a inset) with the predicted rotational energies
for a rigid asymmetric top with H2O equilibrium bond lengths/angles. As a moresignificant challenge to the model, energy levels could be similarly converged for

76
floppy XH3 molecules modeled after H3O+. Indeed, the expected progression ofumbrella inversion levels was reproduced even for potentials modeled after the
relatively high-tunneling barrier case of NH3. This pattern is summarized in Figure
3.1b, where the observed tunneling splittings are given for each pair of converged J =0 energy levels and exhibit a dramatic increase with energy above the barrier. It is
worth noting that apparent deviations from this pattern simply arise from a doubly
degenerate asymmetric bend state occurring near 5000 cm-1, which, as expected,
exhibits a “normal” tunneling splitting comparable to the ground state.
Most important, the POS model also provides access to rovibrational quantum
eigenfunctions as well as eigenenergies. By way of example, Figure 3.2a displays a

1D projection of the full 4D wavefunction for XH2 onto the γ coordinate, where γ isthe H-X-H bending angle between the two hydrogens. For comparison, I also show acorresponding slice through the J = 0 wave function (shaded in grey) obtained from
Diffusion Monte Carlo (DMC) methods, where the stiffness (i.e., anisotropy of thepotential) χ of the intermolecular bend potential is systematically incremented
between each plot (χ ≈ 3700 corresponds to the bending potential in H2O).
Agreement between POS and DMC predictions for XH2 is excellent, limited only by
the intrinsic statistical nature of the DMC wave functions. Similar levels of agreement
are evidenced between 1D representations (solid points) of the full 6D XH3 wave
functions and statistical DMC calculations (shaded in grey), which are plotted in
Figure 3.2b as a function of increasing stiffness of the intermolecular potential. Hereγ refers to the angle between any pair of H’s with respect to the central atom, andthus both the POS and DMC Ψ(γ) should be more appropriately thought of as

77
pairwise angular-distribution functions. As expected, each of the hydrogen pairs
become more tightly correlated with increasing potential stiffness, with the angle
strongly localizing around the equilibrium value (104o for XH2 and 108o for XH3).
However, the potential promise in such a POS approach lies in permitting one
to theoretically move past n = 2 and n = 3 to more challenging systems of higher
dimensionality, which represent the focus of the present work. The organization ofthis chapter is as follows. Section II reviews the general theoretical backgroundnecessary for efficiently solving the POS problem and focuses explicitly on methods
used in this work for the n = 4 and n = 5 particle system. Sections III and IV discuss
applications of the particle-on-a-sphere model to two test systems, n = 4 (XH4) and n= 5 (XH5), respectively. I specifically model the rotation-bending potentials on
suitably scaled versions of i) relatively “stiff” species such as methane (CH4) and ii)
highly fluxional species such as protonated methane (CH5+). In Section V, I compare
the n = 5 POS predictions with other large amplitude theoretical results for CH5+.
Concluding comments are summarized in Section VI.

II. Theoretical BackgroundAs discussed more fully in Chapter 2, the POS model assumes a massive
central atom (mX >> mH) and effectively averages over the high-frequency radialstretching of the X-H bond lengths, so that the resulting 2D angular coordinates canbe simply represented as motion constrained to the surface of a sphere. The reductionof dimensionality can be simply achieved by freezing the radial coordinates in the full
dimensional potential at some average value, or more reasonably, by adiabaticallyadjusting the radial coordinates to a minimum energy as a function of the remaining

78
2n angular coordinates. The justification for the latter approach is clearly superior,and in the case of CH5+, takes advantage of both the greater than twofold difference in
CH stretch vs. torsional / bending frequencies as well as the relatively small spread (<0.015 Å) in CH equilibrium bond lengths observed for H atoms at typical zero-point
energies. The resulting Hamiltonian for the large-amplitude angular motion of nparticles can then be simply expressed asn

H POS = ∑ bi ˆji2 + Vˆ (Ω ) ,

(3.1)

i =1

where ĵi is the angular momentum of the ith hydrogen with respect to the central
stationery atom, and bi is the rotational constant for motion on the sphere. Vˆ (Ω ) is
the potential describing the H atom interactions, where Ω is a 2n dimensional vectorof all n H atom angular coordinates.
One key advantage of this approach is that there is no distinction between
angular bending and end-over-end tumbling degrees of freedom. Thus the solutions
treat vibration and rotation on a completely equivalent footing (i.e., there are no
perturbative assumptions about the magnitude of Coriolis interactions). The methodshould work best in the limit of extreme large-amplitude quantum motion (e.g., CH5+)
and relatively weak interparticle interactions, as opposed to the standard decouplingof Ψtot into a product of rotational and vibrational wave functions describing small
amplitude displacements with respect to a zeroth-order body-fixed frame. For
example, such a method would be particularly well suited for describing the quantum
motion of n He atoms solvating a heavy central atom, as in van der Waals clusters of
XeHen. However, the corresponding cost of equality with respect to rotation/vibration

79
internal coordinates is a much more rapid scaling of the variational problem withnumber of particles. Thus, although n = 2 (4D) and n = 3 (6D) systems aretheoretically tractable from a variety of approaches, extension of the POS methods ton = 4 and 5 invoke special considerations, which are described below.

A. Primitive and coupled basis sets
As discussed earlier,52 the goal of the POS model is to describe floppy
molecules such as CH5+, where the Coriolis coupling between overall rotation and HX-H angular vibrations is large. In the “floppy” limit, where the potential V(Ω) canbe viewed as a perturbation of freely rotating molecules, a natural basis in which toexpand the Hamiltonian is simply a direct product of rigid rotor functions:n∏ ji mi ,

(3.2)

i =1

The full n-particle Hamiltonian is diagonal in total J2 and Jz, so it is advantageous totransform this primitive basis into an equivalent coupled representation53 that permits
matrix construction and diagonalization for each quantum number J. For increased
speed in the potential matrix element evaluation, my choice of a coupled basis for the
XH4 system is(( j1 j2 ) j12 ( j3 j4 ) j34 )JM ,(3.3)

and for the XH5 system,(( j1 ( j 2 j3 ) j23 ) j123 ( j4 j5 ) j 45 )JM ,(3.4)

In these expressions, ji represents the angular momentum for the ith particle, jij is thevector sum of ji and jj, parentheses denote the coupling order, and J and M represent
the total angular-momentum quantum number and projection on the z axis,

80
respectively.

In the absence of external fields, all M states for a given J are

degenerate; thus M can be taken to be 0 for a further factor of 2J + 1 reduction of the
Hamiltonian matrix size. To create a basis set for a given J, the primitive free-particle
angular momenta are restricted to j = [0, jmax].

B. Basis set symmetrization and permutation inversion
To achieve a crucial reduction in computational effort, the Hamiltonian isblock diagonalized using group theoretical methods.54,55 This involves the use of
permutation inversion theory55 to create symmetry-adapted linear combinations(SALCs) of the coupled basis functions that transform according to each permutation
inversion symmetry group. The permutation inversion group for XH4 is G48(isomorphic to the direct product Td symmetry group with the inversion operation),
which yields 10 irreducible representations: A1+, A2+, E+, F1+, F2+, A1-, A2-, E-, F1-,and F2- with nuclear spin weights of 5:1:3 for A1:E:F2, respectively, for spin ½ Hs.
For XH5, the permutation inversion group is G240, characterized by 14 irreducible
representations: A1+, A2+, G1+, G2+, H1+, H2+, I+, A1-, A2-, G1-, G2-, H1-, H2-, and I- with
nuclear spin weights of 6:4:2 for A2:G2:H2 symmetries,50 respectively.

Using

standard projection operator techniques,55 the coupled basis functions are transformed
into orthonormal SALCs of each irreducible representation. For XH5 in particular,
these irreducible representations reflect high levels of intrinsic degeneracy,
specifically E = 2, F = 3, G = 4, H = 5, and I = 6. By exploiting this degeneracy andincluding only one of the resultant SALCs for each irreducible representation,56,57 I
therefore dramatically reduce the basis set size.

81
It is important to note that while such a symmetrization need only be
performed once, the creation of these SALCs currently represents the limiting step in
POS extension to even larger basis sets. The nested loops required to symmetrize andorthogonalize the SALCs quickly grow extremely rapidly with increasing n and jmax.
Specifically, the time required to create an orthonormal SALC basis set of a single
irreducible representation increases as (jmax)7. By way of example, in order to create
an orthonormal basis set of SALCs for one symmetry with n = 5, J = 0, jmax = 5(Ncoupled ≈ 75,000) requires about 1 month to symmetrize and orthogonalize on a 2.4
GHz Intel CPU. This translates into a current upper limit for the XH5 problem of jmax
= 5, 4 and 3 for J = 0, 1, and 2, respectively. However, this process may besignificantly accelerated by parallelization of the code, which represents one ofseveral ongoing directions of further exploration.

C. Matrix element evaluation
The use of symmetrized basis functions significantly decreases matrix size forlarge n and jmax, reducing both memory requirements to store the matrix as well as
computational effort for matrix diagonalization. However, calculating the
Hamiltonian matrix in a symmetrized basis still requires significant CPU time. Matrix
element evaluation in the SALC basis requires calculation of the integral i Hˆ i' , thatis expanded into a sum of integrals over the direct product of each term in the SALCs.
Naively, when the basis function i is a SALC created from a sum of coupled basis
functions and the Hamiltonian matrix is too large to store in memory, integral

82
evaluation might be expected to occur via a process that scales quadratically with
SALC size:
M SALC M SALCi Hˆ i' = ∑k =1∑ C k*Cl k Hˆ l ,

(3.5)

l =1

where k and l are indices over each of the MSALC coupled basis functions in each of
the NSALC SALC basis functions.

Consequently, creation of the symmetrized

Hamiltonian matrix would increase roughly as the 4th power, i.e., O(N4), of basis set
size. More specifically, the number of evaluation of integrals required is ½ x NSALC2
x MSALC2, where NSALC is the number of SALC basis function for a given symmetry
and MSALC is the number of coupled basis functions in the linear combination thatcomprises the SALC. Even for small values of jmax, this number becomes quite large,
requiring for J = 0, jmax = 5, and n = 5 on the order of 109 integral evaluations to create
the Hamiltonian matrix; this calculation requires 1 to 2 weeks on a 2.4 GHz Intel
CPU, even with analytic matrix elements.
However, by changing the order of evaluation, the matrix element evaluation
can be expressed as an O(N3) operation, which is significantly faster when there are alarge number of symmetries into which the SALCs can be distributed. Rather thancalculating the full NSALC x NSALC Hamiltonian matrix ( H Γ ) from direct products of
SALCs, the matrix can be expressed as H Γ = S TΓ ⋅ H ' ⋅ S Γ , where H ' is the nonsparse

Ncoupled x Ncoupled unsymmetrized Hamiltonian, and S Γ is a NSALC x Ncoupled sparse
matrix that transforms the coupled (unsymmetrized) basis functions into the SALC
basis set. The advantage is that S Γ is quite sparse (e.g., for XH5, ∼95% of the matrix
elements are zero); thus computing terms in H ' (which itself is too large to store for83
typical Ncoupled ≈ 75,000), multiplied by nonzero terms in S Γ is easily accomplished.
Ultimately, this rearrangement of operations results in Ncoupled x NSALC x MSALC +NSALC2 x MSALC matrix element evaluations, i.e., reducing to an O(N3) algorithm. Thisis approximately (NSALC x MSALC)/(Ncoupled + NSALC) times faster than the naive direct
product approach, resulting in a 50- to 100-fold increase in speed for large basis setswith many symmetries.
By careful choice of the potential, the evaluation of both the kinetic energy
and potential energy matrix elements is analytic. Since the basis set is diagonal in ji2,
the kinetic energy operator matrix elements satisfy this requirement automatically.
For XHn, these matrix elements simply are:ni Tˆ i ' = b ∑ δ ii ' j k ( j k + 1) ,

(3.6)

k =1

where b is the internal X-H rotational constant and jk is the angular momentum of the
kth X-H rotor. Calculation of the corresponding potential matrix elements involves
substantially greater effort.

Evaluation of the potential by standard quadrature

methods requires Q2n calculations per matrix element, where Q is the number ofpoints on a potential grid and n is the number of hydrogens on the sphere, whichscales too rapidly when n > 3. To proceed, therefore, I consider the potential as amultibody expansion, i.e.,
n

n

i< j

i< j <kV (Ω ) = ∑ V ( 2 ) ( Ω i , Ω j ) + ∑ V ( 3) ( Ω i , Ω j , Ω k ) +n∑ Vi < j < k <l

( 4)

,

(3.7)

(Ω i , Ω j , Ω k , Ω l ) + ...

where each successive term represents a sum over all 2-body, 3-body, 4-body, etc.,
contributions. As a major simplification, I restrict this expansion to the first term, i.e.,

84
only considering linear combinations of pairwise additive intermolecular potentials
between each particle that depends only on the angle between pairs of particles:n( )V (Ω) = ∑ V ( 2) γ ij ,
i< j

(3.8)

where γij is the H-X-H angle between the ith and jth hydrogen and V(2) is a sum ofLegendre polynomials in γij. By truncating after the first term in the expansion, thework required to calculate a potential matrix element only grows as n2 (the number of
pairwise interaction for n particles), as opposed to Q2n with numerical quadrature (Qbeing the number of quadrature points per dimension). As a further bonus, if I expand
any pairwise additive potential in a sum of Legendre polynomials, the matrix element
for each term can be calculated analytically via Clebsh-Gordan angular momentum
algebra.52,58 The matrix element between two coupled particles and a Legendre
function are the Percival-Seaton coefficients58, which are simply sums over 3- and 6-J
symbols. Thus one has gained the advantage of both an n2 vs. Q2n scaling of integral
evaluations as well as a much more efficient analytic evaluation of each matrix
element. Although such a pairwise additive approximation to the potential was notessential for calculations in n = 2 and n = 3 systems, it proves absolutely crucial toextension to the larger systems of interest such as XH4 and XH5.

D. Rigid-body DMC
In Chapter 2, direct comparison between the POS model and exact full
dimensional quantum mechanical calculations proved feasible for low J states in XH2
and XH3, due to an acceptably small number of degrees of freedom. For the much
more computationally demanding XH4 and XH5 systems, exact full dimensional

85
quantum calculations for J > 0 are simply not available for benchmarking the POS
calculations. However, calculations for the nodeless ground rovibrational state (J = 0)
are feasible for n = 4 and 5 via DMC methods, which scale much more favorably thanvariational methods for systems of high dimensionality (roughly linearly in 2n).
Thus, I can test the convergence of the J = 0 POS methods against DMC results,
specifically exploiting the “rigid body” formulation (RBDMC) to constrain the radialstretching coordinate for the diffusing 2n dimensional “walkers.”59-62
By way of example, Figures 3.1a and 3.1b show the energies of XH2 and XH3
as a function of jmax (i.e., basis set size). As jmax is increased, the energy of each state
decreases variationally, with the ground state J = 0 energy converging on the DMC
result (dotted line). The insets show a more detailed picture of agreement between
variational POS and DMC POS methods. For both XH2 and XH3 calculations, theconverged variational calculation equals the DMC result within 1σ statisticaluncertainty, where the DMC error estimate reflects 10 ensembles of 2000 walker
trajectories starting with different initial conditions.
The distribution of walker angular coordinates also provides a statisticalsampling of the wave function. Figure 3.2a displays the DMC wave function forwalkers at an angle γ apart from each other in XH2, binned in 1o increments andnormalized to unity; the black dots represent the wave function results calculated by
the variational POS method. Note that when the angular stiffness (χ) is low, the two
H atoms rotate independently of each other and the distribution reduces to sin(γ), i.e.,a uniform sampling of interangular coordinates on a sphere. However, as the stiffness
χ is increased, the wave function becomes localized around the equilibrium angle

86
corresponding to H2O (γo = 104o). A similar pair of POS and DMC wavefunction
plots for XH3 is shown in Figure 3.2b, where γij represents the angle between any twoof the particles. Here the equilibrium angle corresponds to that of NH3 (γo = 108o),and thus the projection of the wave function peaks near 108o for a stiff (i.e., large χ)
potential and becomes less localized as the potential is softened. Agreement between
the DMC and variationally calculated POS wave functions for each value of χ isexact within statistical uncertainty limits.

III. Four Particles-On-A-Sphere(XH4)
Building on the above successes for XH2 and XH3, the next test system toinvestigate with the POS model is XH4. The rotor and bending potential parameters(bi = 16.63 cm-1 and γ0 = 109.47o) are specifically chosen to recapitulate theequilibrium structure of methane, CH4, with a V(2)(γ) stiffness in the pairwise
potential chosen to yield the correct CH-bending frequencies. As I will not in general
be able to converge POS calculations for such a strongly localized intermolecularpotential, it will be useful to introduce a scaled function, Vscaled(γ) = αV(2)(γ), toreflect an angular potential “softer” than methane by a multiplicative factor of α. In
the stiff limit (α = 1), my POS model of methane should yield a spherical top rotor
with two ν2 and ν4 bending states corresponding to the E (~1526 cm-1) and F2 (~1306
cm-1) symmetry vibrations, respectively. At higher resolution, the 2J + 1 K level
degeneracy of a rigid spherical top will be lifted by centrifugal interactions inrotationally excited states, which will lead to additional rotational fine structure for J

87> 0. See Table 3.1 for a list of the coefficients used in the Legendre expansion for theCH4 potential.
I begin the POS investigation of XH4 by exploring the rate of convergence asa function of basis set size. In Figure 3.3, the vibrationless state energies for J = 0 - 4
are shown as a function of jmax, based on a scaled methane potential with a stiffness of
α = 0.033. By way of comparison, the dashed line in Figure 3.3 represents the J = 0
DMC ground-state energy for the same potential; the inset illustrates agreement of
eigenenergies within the 1σ error of the DMC calculation. The size of the basis set
increases extremely rapidly with jmax. For example, for J = 4, Ncoupled increases from
16,429 to 39,046 to 80,866 for jmax = 5, 6, and 7, respectively. As a result, jmax = 7 isthe largest calculation that can be performed for the 4 particle POS model with J = 4(for J = 0, jmax can be as large as 9). This choice of α reflects a bending potential thatis ~30-fold softer than actual methane, but nevertheless permits systematicconvergence of POS calculations for a series of rotationally excited states within
computational limits of time and resources. The plots in Figure 3.3 demonstrate thatas jmax is increased up to the maximum value, the variational POS energies decrease
monotonically to an asymptotically converged value. If I empirically characterize
convergence by the incremental decrease in POS eigenenergies for an incremental
jmax- 1 to jmax increase in basis set size, then jmax = 7 corresponds to a 0.001 cm-1
convergence for the lowest rotational state J = 0, with only slightly reduced (0.01 cm1) levels of convergence for excited rotational states J = 1 - 4. At an even finer levelof detail, the jmax = 7 energies clearly demonstrate converged spherical top rotational

88fine structure, with centrifugal distortion lifting the nominally degenerate 2J + 1 K
energy levels for each J.
The progression of the v = 0 energies as a function of α is shown in Figure
3.4a, as the CH4-like potential is increased from α = 0.0001 (i.e., nearly a free rotor
limit) up to the stiffest calculations (α = 0.033) that can be currently converged. In the
α = 0.0001 limit, the 4 individual C-H rotors rotate more or less independently,
resulting in energy-level spacings characteristic of a large single C-H rotor constant.
As the potential is stiffened, these states converge onto what approximates a ~J(J + 1)
spherical top progression of energies, with the fine structure splittings in Figure 3.4a(shown in detail in Figure 3.4b) due to centrifugal bend-rotation coupling.
Noteworthy is that while the effective B for α = 0.033 is only 20% above the true
value, the fine structure splittings are ∼500 times larger, albeit in the correct
symmetry order with even qualitatively correct splitting ratios. This implies that evenfor a 30-fold scaled-down internal anisotropy, the end-over-end tumbling structure of
the 4 particle system is already well defined, with the fine structure effects requiringmuch larger basis sets (as well as the inclusion of radial degrees of freedom) toconverge.
In addition to comparison of POS vs. experimental rotational energies in the
ground vibrational state, I can also inspect the low frequency vibrational states ofmethane (ν2,E and ν4,F2) in the rotationless J = 0 level (see Figure 3.5). Since it waspossible to use a larger jmax for J = 0 calculations, the energies have been convergedwith a 10-fold stiffer α = 0.33 than for J > 0, but still somewhat short of the fullmethane potential. As the potential approaches the full methane stiffness, the energy

89of the ν2,E and ν 4,F2 vibrations increase monotonically towards experimental
values,63 with the order of the ν2,E and ν 4,F2 vibrations maintained even for
substantially floppier potentials. Similar qualitative agreement is noted with twoquanta of bending energy, as shown in Figure 3.5 for 2 ν 2,A1 + E, ν 2 + ν 4,F1 + F2,
and 2 ν 4,A1 + E + F2. Analogous to the calculations presented earlier, the ordering of
the dyads is correctly predicted, with the vibrational fine structure splittings
decreasing as the potential stiffens . Although 4 particle POS calculations cannot be
fully converged up to α = 1, it is important to note that many of the qualitative
patterns are already quite visible and correctly predicted at potentials substantially
softer than the full potential.

IV. Five Particles-On-A-Sphere (XH5)
The overall target of the POS model is to expand to 5 particles on a sphere,
which I hope will provide initial insight into the challenging fluxional CH5+ problem.
I model the individual C-H rotor with b = 12.81cm-1, i.e., consistent with expectation
values of the C-H bond length from full 15D DMC calculations.34,64 As discussed
earlier in this paper, my choice of i) a pairwise additive potential with ii) each term
expressed as a sum of Legendre functions makes matrix element calculation of thepotential operator analytic. To take advantage of this, the full 15D MP2/cc-pVTZ
Bowman CH5+ potential34 needs to be expressed in this sum of pairwise potentials.
To accomplish this, the radial coordinates of the 15D potential are first adiabaticallyrelaxed in the “fast” CH stretch coordinate to create a 10D potential that depends only
on the relative angular coordinates. The critical point geometries are shown in Figure

903.6. The minimum energy structure of CH5+ on the Bowman potential is of Cs
symmetry, consisting of what looks like an H2 perched on top of a CH3+ structure.
However, at the wave function level, there is rapid quantum exchange between “H2”or “CH3” moieties, resulting in complete delocalization of each of the H atoms
between the 120 equivalent geometries separated by low barriers. For example, theH2 (or CH3+) moiety can “rotate” through a Cs transition state (E = 43 cm-1) or “flip”
H’s between H2 and CH3+ through a slightly higher C2v transition state (E = 192 cm1). Most importantly, this process ensures that angular geometries and energies of all

critical points for the full 15D and reduced 10D potential are identical.
The second step is to represent the 10D potential by a sum of pairwise
potentials, which is achieved by least-squares fitting 10,000 randomly selected points
with energies below 1,000 cm-1. The energies and angular geometries of these points
are then fit to a sum of pairwise terms, with each term described by a sum ofLegendre functions (See Table 3.1 for a list of the coefficients used in the Legendre
expansion for the CH5+ potential). Of crucial relevance is the reproduction quality forcritical points and the “flip” and “torsion” reaction paths between C2v and Cs minima.
The critical points of the angular coordinates are shown in Figure 3.6 for i) the 10D
radially relaxed version of the full 15D Bowman potential and for ii) the pairwise fit
potential. Given the enormously simplifying nature of the pairwise additive
approximation, the agreement between 10D and pairwise critical points is quite
remarkable. Similar accuracy is also demonstrated in Figure 3.7 for both the 10-D
potential and the pairwise fit along the Cs (internal rotation) and C2v (“H2 flip”)

91
reaction paths. Again, the differences between the 10D potential and the pairwise fit
are surprisingly small.
The availability of DMC methods provides one more opportunity to test thevalidity of the pairwise approximation for both the 10D radially relaxed as well as thefull 15D potentials. Specifically, DMC is able to calculate the ground state wave
function at all three levels of dimensionality (POS, 10D and 15D), each of which canthen be projected on to the γ coordinate (where γ is the angle between any twoparticles) for more quantitative comparison. The results from such a comparison aresummarized in Figure 3.8, which reveals quite good agreement at the wave functionlevel. (Note that the 10D and 15D histograms are offset from 0 for better visualclarity with respect to the pairwise potential histogram.) The main peak in each of the
wave functions near 110o corresponds to unresolved features due to i) the H-C-H
angles in the CH3+ motif, as well as ii) the angle between one of the H2 hydrogens anda CH3+ hydrogen.

The smaller wavefunction shoulder near 60o corresponds tocontributions from iii) the smaller H-C-H central angle between the H2 hydrogens.
Treatment at the POS, 10D and 15D levels yield reduced interatomic distributions
that are nearly identical, and indeed indistinguishable within statistical uncertainty forthe pairwise fit and full 10D potential versions of the DMC wave functions. Though
clearly not exact, this suggests that the use of a pairwise potential still captures the HC-H angular correlations at a high level and, at the very least, is not the significant
source of error in the reported calculations.
In analogy to the strategy for the n = 4 particle problem, I first look at the
ground state (J = 0) zero-point energy for n = 5 particle POS calculations as the

92potential (asymptotically the pairwise 10D fit of the full 15D Bowman CH5+
potential) is stiffened from α = 0.0001 up to α = 1.0. Figure 3.9 illustrates thevariationally calculated zero-point energy (solid line) from the POS model, withconverged DMC energies and uncertainties for the 10D pairwise (solid circles) andnonpairwise 10D potentials (open squares). (The 15D zero-point energy can also bereadily calculated from DMC, but is irrelevant for the present comparison due to an
additional ∼7000 cm-1 in the radial stretch coordinates). The agreement between POS
eigenvalue predictions for the pairwise and exact 10D results is excellent at low α,
confirming the utility of the pairwise approximation to the CH5+ potential. Of morequantitative relevance, however, this comparison illustrates the degree of
convergence (or lack thereof) of the POS calculations at a spectroscopic level as afunction of potential stiffness. Up to α = 0.01, the ground-state energies appear to bequite converged to < 1 cm-1, with this quality of convergence degrading substantiallyas α approaches unity.
More specifically, I can monitor the convergence for a particular value of α
(0.01) with respect to basis set size (i.e., jmax). In Figure 3.10, the ground stateenergies for each of the 6 nonzero nuclear spin weight symmetries50 for CH5+ are
plotted as a function of jmax. Note that the totally symmetric ground state (A1+),
which is the only energy that can be rigorously compared with DMC results (J = 0), is
not represented, has zero statistical weight in CH5+, and is not shown. Also included
along the ordinate axis is the size of the basis set for each jmax for J = 0, 1, and 2.
Note the rapidly growing size of Nprim as jmax increases, which is why jmax = 5 is thelargest basis set that can be implemented at the current time for a J = 0 calculation,

93
with jmax = 4 and jmax = 3 for J = 1 and 2, respectively. As mentioned previously, thisconstraint is due to the maximum size coupled basis set that can currently besymmetrized on a single processor with non-parallelized code.

V. Discussion
The difficult theoretical issues posed by CH5+ have been addressed bynumerous ab initio potential surface investigations over the years; however, efforts topredict detailed rovibrational energy levels for a given potential surface have proveneven more challenging. Landmark work by Bunker and co-workers,50 using arotation-contortion Hamiltonian,49 has been used to predict the J = 0 Å 1 microwave
spectrum of CH5+, exploiting symmetry correlations between the two extremes
associated with a rigid Cs and C2v potential minima. In this work, considerable effort
was made to calculate energies for states of nonzero statistical weight at a number ofdifferent flip and torsion rearrangement barriers, based on a Hamiltonian which
adiabatically separates one large amplitude internal contortional mode from theremaining 11 “fast” vibrational modes. The Hamiltonian is a sum Hrot (end-over-end
tumbling of the molecule), Hτ (motion along the contortion degree of freedom), andHrot,τ (coupling between tumbling and contortion), with all coupling neglected
between the “fast” vibrations and overall rotation. As an alternative formulation, thePOS model neglects the contributions from only the 5 “fastest” CH-stretchingvibrations but does include the 7 low-energy bending vibrations and allows full
coupling of these vibrations to end-over-end tumbling of the CH5+ molecule. The
corresponding disadvantage of the POS model, however, is the significant additional
challenge to achieve convergence for realistically stiff versions of the potentials. This

94
motivates a brief discussion of the POS predictions with the rotation-contortion
results, with an eye toward establishing additional insights into the dynamical trends
anticipated for CH5+.
The 53.9 cm-1 torsional and 213.9 cm-1 flip barriers from least squares
pairwise fits to the full 15D Bowman surface (i.e., α = 1) most closely approximate
the reduced dimensional rotation-contortion model barriers of 50 cm-1 and 200 cm-1,
respectively. I am currently unable to converge CH5+ calculations for α = 1 andtherefore compare results as a function of increasing potential stiffness. In this spirit,
converged POS eigenenergies for the lowest 10 J = 0 states of nonzero nuclear spin
statistical weight are presented in Figure 3.11 for the progression of α = 0.0001, α =0.001, and α = 0.01. The corresponding J = 0 energy levels obtained from the fullrotation-contortion model for the same series of nuclear-spin symmetry states are also
shown in Figure 3.11. In each column, the lowest energy state is defined to be zero;thus the relative energies are presented as a function of model and potential stiffness.
At the simplest level, the essentially unhindered free-rotor energy patterns
observed at low α evolve to more complex internal rotor splittings with increasing
potential stiffness. A general monotonically increasing trend in rotor energy spacings
with stiffness is clear, with an energy ordering pattern that is largely (but not entirely)
maintained with respect to α. The final extrapolation from α = .01 to the rotationcontortion model represents a large change in both the potential and level ofdynamical approximation, though much of the energy ordering is maintained as theenergy separations increase. However, arguably the most striking aspect is that theconverged energy level differences for stiffest (α = 0.01) POS calculations, though

95
still far from the full potential (α = 1), have already achieved magnitudes
qualitatively comparable (i.e., within 2- to 3-fold) with respect to the rotationcontortion model predictions. We find this result particularly surprising. However,such behavior is likely to be more consistent with more free internal rotor statebehavior rather than localized vibrational or librational motion, for which thespacings would presumably continue to scale with the square root of the potential
stiffness, α1/2.
Since the convergence of the CH5+ and CH4 calculations is limited by the sizeof the basis set into which the POS Hamiltonian is expanded, the size of the basis set
is limited by the amount of time needed to create the orthonormal SALCs from thecoupled basis sets. Creating the SALC from standard methods54,55 requires creating
linear combinations of all of the permutations of identical particles in the basis set.
For the 4 and 5 particle problems, the length of the linear combinations of these
SALCs, MSALC, can be on the order of 1,000 coupled basis functions. For a given
symmetry, each of these SALCs must be orthogonalized with respect to each other.
As an example, in Table 3.2 the actual time required to create a J = 0 SALC for XH5
is given as a function of jmax. The jmax 6 and 7 values are estimated from the timerequired for jmax = 1 - 5. This reinforces the previous assertion that the time required
to create the symmetrized basis sets, rather than creation or diagonalization of theresulting Hamiltonian matrix, is currently the rate-limiting factor for achieving
converged calculations at higher stiffness parameters. We are currently pursuingseveral methods to increase the efficiency of this process with parallelization of the
code.

96
As a final note, however, we can roughly estimate the size of jmax which will
be necessary to achieve convergence for CH5+ in the POS model. In Figure 3.12, the
RBDMC wave function is projected onto the γ coordinate for α = 0.001, 0.01, 0.1,
and 1. As would be expected, the α = 0.001 wave function looks almost identical to
the sin(γ) Jacobian that would be expected for independent particles.

As α isincreased to 1, more structure in the wave function appears, which requires sufficientangular flexibility in the multidimensional wave function. Therefore, a crude measureof the jmax value required to converge the α = 1 POS calculation is to expand thecorresponding 1D DMC pair correlation function in the H-C-H angle γ as a sum ofLegendre functions Pl(γ) and to probe the degree of convergence with respect to l.
The physical notion is that angular flexibility in a sum of Legendre functions with l≤ jmax mirrors the angular flexibility for a sum of SALCs obtained for the same jmax.
The assumption is that the SALCs are simply sums of spherical harmonics which are
the 2D analog to 1D Legendre polynomials; thus with an expansion of the same size,
the flexibility of the two functions should be similar. Figure 3.13 displays the 1D
DMC wavefunction along with several expansions in Legendre polynomials,
indicating a rapid visual convergence with increasing jmax. More quantitatively, if the
expectation value of the energy for a 1D rotor in the α = 1 pairwise potential iscalculated using the wave functions expanded in Legendre polynomials of increasingorder (jmax), the convergence of this energy as a function of jmax should mirror thenumber of Legendre functions required to accurately represent the 1D wave function(see Table 3.6). To converge the energy to about 1% requires jmax ≈ 9, which withcurrent basis set creation algorithms would take about 10 CPU-years to create.

97
However, once a basis is created, the POS method should do exceptionally well
reproducing the correct dynamics. The validity of this simple model can be tested
with a converged POS wave function of lower dimensionality. For example, the POS
model for NH3 is converged to better than 0.01 cm-1 for a jmax = 15; by expanding thisconverged n = 3 1D pair correlation function in a sum of Legendre polynomials up tol = 15 and calculating the expectation value of the energy, we can compare the rate ofconvergence of the energy with respect to jmax. Table 3.7 shows the error in theenergy of converged POS calculation as a function of jmax compared to the error
expectation value of the energy of the 1D expansion of the pair correlation function;
both energies converge at about the same rate, implying that the 1D expansion gives agood first order estimate of the size of basis set needed to converge a full POS
calculation.

VI. Summary and Conclusion
In this paper, a simple, yet surprisingly powerful method for the approximate
quantum treatment of large-amplitude motion in systems with many degrees of
freedom has been tested with large molecules. The POS method takes advantage ofthe approximate independence of the XH bond length on the angular coordinates tomotivate a reduced dimensionality Hamiltonian treatment of the motion of nhydrogens constrained to the surface of a sphere. For n = 2 and 3, the problem isquite tractable and solutions of very stiff systems are easily converged, although not
as efficiently as a method designed for rigid molecules. For the harder n = 4 and 5,
the problem space grows so quickly that converged results are obtainable only forvery floppy systems. However, the important energy level patterns discovered in the

98
floppy system have been shown to correlate well to the stiff systems and give insightto the dynamics of these complicated systems. The aim of the POS method is toprovide a new tool to gain insight into the quantum dynamics of systems that display
large amplitude motion. This reduced dimension method is necessary for systems
such as the highly fluxional CH5+ where the full dimension quantum calculations are
not possible.

99

TablesTable 3.1. The Legendre coefficients for the α = 1 CH4 and CH5+ potentials
l
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

-1

-1

CH4 al (cm )CH5+ al (cm )3892.709
4999.699
6942.05701561.9630705.05270403.3715
0261.5983
0183.51920135.90250104.7094083.15934067.64614

4190.442
3060.609
14790.49
1821.198
9008.753
2796.952
5626.363
1216.296
2174.996
238.9801
381.5193

100
Table 3.2. Time required to create an orthonormal basis set of one irreducible
representation in the G240 molecular symmetry group.
jmax
1
2
3
4
5
6*
7*

*

Time
50 s
100 m
20 hr
3d
30 d
110 d
1 yr

Predicted time based on scaling observed in previous calculations.

101
Table 3.3. First five J = 0 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.
α

0.001

0.01

0.1

246.86 (0.02)
322.08 (0.01)
375.57 (0.11)
406.43 (0.20)
425.34 (0.00)
298.44 (0.01)
325.38 (0.21)
406.30 (0.16)
425.88 (0.11)
458.43 (0.37)
250.83 (0.03)
255.60 (0.02)
306.40 (0.18)
322.84 (0.02)
329.75 (0.02)
108.95 (0.04)
192.39 (0.05)
216.28 (0.01)
277.00 (0.04)
280.09 (0.04)
144.94 (0.08)
173.88 (0.04)
217.60 (0.06)
245.34 (0.03)
253.30 (0.04)
191.93 (0.05)
198.73 (0.05)
214.81 (0.02)
276.26 (0.05)
279.29 (0.04)

479.0 (1.7)
545.0 (1.7)
609.6 (7.2)
622.0 (7.7)
630.1 (1.9)
522.5 (3.0)
551.7 (12.0)
632.5 (8.5)
645.5 (9.8)
698.0 (15.1)
472.5 (1.9)
498.5 (3.3)
546.8 (7.6)
568.1 (3.8)
583.1 (10.2)
304.2 (3.1)
399.0 (2.9)
425.4 (3.1)
506.8 (2.6)
527.7 (4.8)
379.7 (5.7)
410.1 (3.1)
464.4 (5.2)
474.5 (2.6)
496.5 (4.0)
384.1 (1.8)
422.3 (3.4)
458.1 (5.2)
507.9 (1.4)
509.5 (5.4)

1713 (195)1836 (160)1906 (307)2011 (300)2113 (341)
1673 (201)1836 (328)2039 (330)2070 (355)2207 (348)1493 (198)
1669 (312)1777 (257)1793 (266)1837 (270)1189 (200)1355 (194)1605 (247)
1631 (256)1783 (270)1297 (154)1479 (180)1579 (332)1682 (252)1749 (228)1273 (133)1408 (265)1623 (283)1712 (254)1730 (307)

1

Γ

A2+

A2-

G2+

G2-

H2+

H2-

10752 (3074)
12125 (3959)
13120 (3636)
13789 (3665)
14343 (4237)
10602 (2997)
12916 (4125)
13676 (4974)
14415 (4500)
14963 (5199)
9475 (2988)
10203 (4937)
11363 (3947)
11957 (3638)
12049 (3963)
7695 (2827)
8547 (2962)
10392 (3839)
10699 (3950)
11417 (3899)
8420 (2015)
8686 (3034)
9909 (4362)
11050 (3639)
11391 (3880)
8234 (2177)
8564 (3655)
9937 (4450)
11039 (3623)
11471 (3856)

102
Table 3.4. First five J = 1 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.
α

0.001

0.01

0.1

1

Γ

A2+

A2-

G2+

G2-

H2+

H2-

240.19 (0.29)
327.76 (0.60)
349.42 (0.06)
373.88 (0.24)
375.02 (0.33)
269.56 (0.25)
292.05 (0.28)
323.56 (0.18)
327.48 (0.13)
374.95 (0.34)
135.60 (0.03)
163.77 (0.25)
220.12 (0.22)
240.29 (0.21)
241.92 (0.31)
189.27 (0.17)
197.61 (0.25)
221.84 (0.41)
269.87 (0.29)
271.25 (0.28)
138.65 (0.05)
163.11 (0.26)
171.42 (0.31)
218.70 (0.18)
221.80 (0.22)
116.78 (0.05)
162.66 (0.01)
190.46 (0.20)
196.12 (0.29)
200.71 (0.27)

430.4 (18.1)
553.6 (22.4)
572.6 (10.8)
590.6 (15.3)
600.9 (27.6)
472.4 (15.4)
486.0 (16.0)
548.9 (17.2)
569.6 (9.9)
597.6 (26.1)
331.7 (4.2)
373.3 (15.7)
427.5 (14.4)
442.6 (15.4)
448.6 (14.7)
392.2 (10.7)
432.8 (17.3)
440.2 (16.2)
479.9 (15.3)
482.4 (15.1)
347.4 (9.4)
374.1 (11.5)
403.0 (16.7)
442.6 (10.9)
447.4 (11.5)
335.6 (5.2)
382.5 (3.4)
402.8 (11.2)
428.8 (22.4)
442.3 (14.5)

1579 (511)1931 (445)2138 (447)2170 (487)2284 (491)1679 (371)1966 (299)2004 (375)2077 (569)2161 (509)1417 (289)1456 (388)1584 (508)1691 (417)1746 (485)
1441 (358)
1588 (433)1716 (372)1730 (484)1913 (369)1419 (346)1482 (357)1581 (456)1641 (451)1683 (503)1453 (226)1493 (321)1589 (456)
1666 (423)1714 (491)11703 (6311)
14787 (5301)
16197 (5755)
16667 (6033)
17719 (6483)
12379 (4635)
15134 (3903)
15332 (4738)
15715 (6675)
16574 (6083)
10280 (4230)
11225 (4249)
11699 (6016)
12515 (5395)
12873 (5780)
10441 (4471)
11709 (5004)
12319 (5300)
12723 (5926)
14012 (5019)
10574 (3859)
10824 (4942)
11644 (5437)
11902 (5899)
12318 (6064)
10587 (3398)
11154 (3898)
11684 (4898)
12037 (5743)
12817 (6154)

103
Table 3.5. First five J = 2 energy levels of each nonzero statistical weighted symmetry
for n = 5. Energies are given in cm-1, with the estimated error given in parentheses.
α

0.001

0.01

238.13 (0.24)
244.33 (0.33)
322.43 (0.39)
327.85 (0.48)
340.42 (0.13)
185.17 (0.46)
289.83 (0.14)
296.19 (0.20)
320.42 (58.92)
322.52167.36 (0.45)
214.69 (0.59)
239.70 (0.38)
244.57 (0.58)
245.48 (0.65)
187.94 (0.62)
190.88 (0.65)
192.94 (0.98)
196.67 (0.93)
216.75 (0.30)
139.78 (1.07)
166.38 (0.60)
169.84 (0.62)
175.52 (0.81)
217.08 (0.89)
113.59 (1.01)
190.62 (0.66)
191.94 (0.73)
196.95 (1.03)
198.97 (1.11)

435.4 (19.2)
484.7 (26.8)
517.3 (6.7)
547.7 (45.2)
560.8 (78.7)
378.8 (32.5)
481.4 (10.9)
521.2 (37.7)
540.4 (131.0)
561.8
388.6 (21.1)
432.2 (35.0)
451.4 (27.8)
474.7 (40.0)
489.1 (40.1)
390.0 (35.1)
411.6 (36.5)
423.7 (37.8)
441.1 (29.5)
452.6 (26.1)
355.1 (43.9)
397.9 (39.5)
406.6 (37.7)
441.9 (33.5)
443.1 (50.1)
334.3 (46.3)
404.2 (26.6)
425.1 (35.7)
435.7 (36.2)
441.8 (45.9)

0.1

1

Γ

A2+

A2-

G2+

G2-

H2+

H2-

1900 (362)2203 (442)2295 (729)2389 (1100)2452 (1327)
1869 (537)2105 (637)2312 (871)2480 (1147)
2593
1835 (576)1912 (882)2039 (833)2078 (896)2159 (835)1840 (448)1855 (695)1929 (747)2121 (699)2169 (813)1760 (667)1851 (742)1877 (822)2093 (793)2149 (792)1816 (512)1838 (755)1923 (782)2055 (923)2121 (877)

15832 (3983)
18616 (5716)
19598 (8367)
20110 (12502)
20806 (14431)
16106 (5628)
17819 (8012)
19908 (9529)
21549 (11662)
22287
15057 (7170)
16113 (9158)
16907 (9736)
17175 (10362)
18068 (9713)
14873 (5662)
15635 (7528)
16458 (8076)
17916 (8268)
18322 (9620)
14741 (7174)
15580 (8266)
15685 (9396)
17325 (9406)
18050 (9285)
14786 (6297)
15940 (8209)
16289 (8588)
17238 (10168)
17855 (10423)

104
Table 3.6. The expectation value of the energy of a 1D rotor in the CH5+ α = 1
pairwise potential.
jmax
1
2
3
4
5
6
7
8
9
10
RBDMC

<E> (cm-1)
13104.011278.2
4638.33644.04543.53322.03019.03097.72820.52783.72777.7

Error
371.8%
306.0%
67.0%
31.2%
63.6%
19.6%
8.7%
11.5%
1.5%
0.2%

105

Table 3.7. The error in ground state energies of the POS model and the 1D Legendre
expansion of the pair correlation function for the converged n = 3 NH3 calculation.
jmax

POS Model

1D expansion

2
4
6
8
10
12
14

330%
110%
41.2%
13.7%
3.53%
0.589%
0.00148%

367%
133%
28.0%
3.17%
0.469%
0.0842%0.00348%

106

FiguresFigure 3.1. Particle-on-a-sphere rotation-bending energies for a) n = 2 (XH2) and b) n
= 3 (XH3) model systems as a function of basis set size (jmax). a) In XH2, theasymmetric top rotational patterns for J = 0,1 corresponding to JKaKc = 000, 101, 111, 110
are clearly evident, with the J = 0 energy rapidly converging to the expected DMC
limit (dotted line). b) Converged J = 0 energy levels for XH3 reveal inversion
tunneling splittings in good agreement with NH3 (after which the potential wasmodeled) and increasing rapidly with symmetric bend excitation. Note the only
exception to this pattern (near 5000 cm-1) arising from an interloper vibrational state
due to degenerate asymmetric bend excitation, which exhibits a tunneling splitting
very similar to the ground state (0.798 cm-1).

107

Figure 3.2. Angular pair correlation functions for a) XH2 and b) XH3 as a function ofpotential stiffness χ = k/b, where k is the angular well depth and b the rotational
constant for a single X-H rotor. For comparison, the shaded function is calculated
from DMC, whereas the dots reflect the fully converged J = 0 POS calculation.

108

Figure 3.3. Convergence of J = 0, 1, 2, 3, 4 energies for XH4 as a function of jmax.
The POS ground state (J = 0) energy is compared against DMC predictions (dottedline).

109

Figure 3.4. a) XH4 rotational energies of J = 0-4 as a function of potential stiffness
scale factor α.

b) At higher resolution, centrifugal induced splittings of thenominally degenerate symmetric top levels for XH4 are shown for α = 0.033, which
are considerably larger but reveal the same qualitative patterns as experimentallyobserved.

110

Figure 3.5. Vibrational energies of XH4 as a function of the potential stiffness scale
factor α.

111

Figure 3.6. Minimum energy paths for CH5+ through the Cs and C2v transition states,
separating the Cs minimum for both the 10D relaxed potential and the pairwise fit.

112

Figure 3.7. Geometry comparison of CH5+ critical points. Note the surprisingly
quantitative agreement between the full 10D relaxed potential and the pairwise
additive fit used in the POS model.

113

Figure 3.8. Pair correlation functions for J = 0, XH5 from DMC calculations for aseries of potential models: a) 15D full CH5+ potential, b) the 10D relaxed potential,
and c) the pairwise additive least squares fit. Note the excellent agreement at thewave function level.

114

Figure 3.9. Ground-state energies of CH5+ as a function of the potential stiffness
parameter (α), calculated by DMC and POS with both the 10D (DMC, open squares)
and pairwise (POS, dashed line) potentials. The inset shows the difference between
the POS and DMC (pairwise) calculations (i.e., the level of convergence) as afunction of α.

115

Figure 3.10. Nonzero statistical symmetries for J = 0, 1, and 2 of XH5 (with a = .033)
as a function of basis set size, jmax and Ncoupled (the number of coupled basis functionsin the basis set). Note the slow but steady convergence of J = 0 energies to the
RBDMC limit with the size of Ncoupled.

116

Figure 3.11. The lowest J = 0 internal rotor energies (with respect to the ground stateenergy, G2- Symmetry) of the POS model compared with the RTRF-like model of
Bunker and coworkers. Interestingly, the splittings for the two model treatments areroughly the same order of magnitude, despite the stiffness of the potentials varying bya large factor.

117

Figure 3.12.

DMC pair correlation function of CH5+ as a function of potential

stiffness scale (α).

118

Figure 3.13. Successive fits of the α = 1 pair correlation function of CH5+ to a sum ofLegendre functions. Such behavior suggests a much improved level of convergence
for the POS model of CH5+ by jmax = 7-9.

119

References
1

V. Spirko, J. Mol. Spectrosc. 101, 30 (1983).

2

S. Urban, R. Dcunha, K. N. Rao, and D. Papousek, Can. J. Phys. 62, 1775
(1984).

3

S. L. Coy and K. K. Lehmann, Spectroc. Acta Pt. A-Molec. Biomolec. Spectr.

45, 47 (1989).
4

V. Spirko and W. P. Kraemer, J. Mol. Spectrosc. 133, 331 (1989).

5

S. Urban, N. Tu, K. N. Rao, and G. Guelachvili, J. Mol. Spectrosc. 133, 312(1989).

6

I. Kleiner, G. Tarrago, and L. R. Brown, J. Mol. Spectrosc. 173, 120 (1995).

7

I. Kleiner, L. R. Brown, G. Tarrago, Q. L. Kou, N. Picque, G. Guelachvili, V.
Dana, and J. Y. Mandin, J. Mol. Spectrosc. 193, 46 (1999).

8

C. Cottaz, I. Kleiner, G. Tarrago, L. R. Brown, J. S. Margolis, R. L. Poynter,
H. M. Pickett, T. Fouchet, P. Drossart, and E. Lellouch, J. Mol. Spectrosc.

203, 285 (2000).
9J. Tang and T. Oka, J. Mol. Spectrosc. 196, 120 (1999).

10

D. J. Liu, T. Oka, and T. J. Sears, J. Chem. Phys. 84, 1312 (1986).

11

P. B. Davies, S. A. Johnson, P. A. Hamilton, and T. J. Sears, Chem. Phys.

108, 335 (1986).
12

M. H. Begemann and R. J. Saykally, J. Chem. Phys. 82, 3570 (1985).

13

M. Gruebele, M. Polak, and R. J. Saykally, J. Chem. Phys. 87, 3347 (1987).

120
14

H. Petek, D. J. Nesbitt, J. C. Owrutsky, C. S. Gudeman, X. Yang, D. O.
Harris, C. B. Moore, and R. J. Saykally, J. Chem. Phys. 92, 3257 (1990).

15

M. Araki, H. Ozeki, and S. Saito, Mol. Phys. 97, 177 (1999).

16

M. Araki, H. Ozeki, and S. Saito, J. Chem. Phys. 109, 5707 (1998).

17

D. H. Zhang, Q. Wu, J. Z. H. Zhang, M. Vondirke, and Z. Bacic, J. Chem.
Phys. 102, 2315 (1995).

18

G. W. M. Vissers, G. C. Groenenboom, and A. van der Avoird, J. Chem.
Phys. 119, 277 (2003).

19

Y. H. Qiu, J. Z. H. Zhang, and Z. Bacic, J. Chem. Phys. 108, 4804 (1998).

20

Y. H. Qiu and Z. Bacic, J. Chem. Phys. 106, 2158 (1997).

21

P. Eggert, A. Viel, and C. Leforestier, Comput. Phys. Commun. 128, 315(2000).

22

N. Goldman, R. S. Fellers, C. Leforestier, and R. J. Saykally, J. Phys. Chem.
A 105, 515 (2001).

23

F. N. Keutsch, N. Goldman, H. A. Harker, C. Leforestier, and R. J. Saykally,
Mol. Phys. 101, 3477 (2003).

24

N. Goldman, C. Leforestier, and R. J. Saykally, J. Phys. Chem. A 108, 787(2004).

25

P. R. Schreiner, S. J. Kim, H. F. Schaefer, and P. V. Schleyer, J. Chem. Phys.

99, 3716 (1993).
26

H. Muller and W. Kutzelnigg, Phys. Chem. Chem. Phys. 2, 2061 (2000).

27

D. Marx and M. Parrinello, Nature 375, 216 (1995).

121
28

G. A. Olah, N. Hartz, G. Rasul, and G. K. S. Prakash, J. Am. Chem. Soc. 117,
1336 (1995).

29

G. A. Olah and G. Rasul, Accounts Chem. Res. 30, 245 (1997).

30

E. Herbst, S. Green, P. Thaddeus, and W. Klemperer, Astrophys. J. 503,
(1977).

31

D. Talbi and R. P. Saxon, Astron. Astrophys. 261, 671 (1992).

32

H. Muller, W. Kutzelnigg, J. Noga, and W. Klopper, J. Chem. Phys. 106,1863 (1997).

33

A. Brown, B. J. Braams, K. Christoffel, Z. Jin, and J. M. Bowman, J. Chem.
Phys. 119, 8790 (2003).

34

A. Brown, A. B. McCoy, B. J. Braams, Z. Jin, and J. M. Bowman, J. Chem.
Phys. 121, 4105 (2004).

35

V. L. Tal'rose and A. K. Lyubimova, Dokl. Akad. Nauk SSSR 86, 909 (1952).

36

E. T. White, J. Tang, and T. Oka, Science 284, 135 (1999).

37

O. Asvany, P. Kumar, B. Redlich, I. Hegemann, S. Schlemmer, and D. Marx,
Science 309, 1219 (2005).

38

X. C. Huang, A. B. McCoy, J. M. Bowman, L. M. Johnson, C. Savage, F.
Dong, and D. J. Nesbitt, Science 311, 60 (2006).

39

J. S. Tse, D. D. Klug, and K. Laasonen, Phys. Rev. Lett. 74, 876 (1995).

40

A. Gamba, G. Morosi, and M. Simonetta, Chem. Phys. Lett. 3, 20 (1969).

41

G. A. Olah, G. Klopman, and R. H. Schlosberg, J. Am. Chem. Soc. 91, 3261
(1969).

42

W. T. A. M. Van der Lugt and P. Ros, Chem. Phys. Lett. 4, 389 (1969).

122
43

W. A. Lathan, W. J. Here, and J. A. People, J. Am. Chem. Soc. 93, 808(1971).

44

P. R. Schreiner, Angew. Chem.-Int. Edit. 39, 3239 (2000).

45

A. Komornicki and D. A. Dixon, J. Chem. Phys. 86, 5625 (1987).

46

W. Klopper and W. Kutzelnigg, J. Phys. Chem. 94, 5625 (1990).

47

S. J. Collins and P. J. Omalley, Chem. Phys. Lett. 228, 246 (1994).

48

M. Kolbuszewski and P. R. Bunker, J. Chem. Phys. 105, 3649 (1996).

49

A. L. L. East, M. Kolbuszewski, and P. R. Bunker, J. Phys. Chem. A 101,
6746 (1997).

50

P. R. Bunker, B. Ostojic, and S. Yurchenko, J. Mol. Struct. 695, 253 (2004).

51

G. A. Olah, Angew. Chem. Int. Ed. Engl. 34, 1393 (1995).

52

M. P. Deskevich and D. J. Nesbitt, J. Chem. Phys. 123, (2005).

53

R. N. Zare, Angular momentum: Understanding spatial aspects in chemistry
and physics. (John Wiley & Sons, New York, 1988).

54

F. A. Cotton, Chemical applications of group theory, 3rd ed. (Wiley, New
York, 1990).

55

P. R. Bunker and P. Jensen, Molecular symmetry and spectrosopy, 2nd ed.
(NRC Research Press, Ottawa, 1998).

56

I. G. Kaplan, Symmetry of many-electron systems. (Academic Press, New
York, 1975).

57

J. F. Cornwell, Group theory in physics: An introduction. (Academic Press,
San Diego, 1987).

58

J. M. Hutson, edited by J. M. Bowman (JAI Press, 1990), Vol. 1A, pp. 1-45.

123
59

V. Buch, J. Chem. Phys. 97, 726 (1992).

60

J. B. Anderson, J. Chem. Phys. 63, 1499 (1975).

61

M. A. Suhm and R. O. Watts, Phys. Rep.-Rev. Sec. Phys. Lett. 204, 293(1991).

62

A. Vegiri, M. H. Alexander, S. Gregurick, A. B. McCoy, and R. B. Gerber, J.
Chem. Phys. 101, 2577 (1994).

63

G. Herzberg, Infrared and raman spectra of polyatomic molecules. (D. Van
Nostrand Company, Inc, New York, 1945).

64

A. B. McCoy, B. J. Braams, A. Brown, X. C. Huang, Z. Jin, and J. M.
Bowman, J. Phys. Chem. A 108, 4991 (2004).

124

Chapter 4: Dynamically Weighted MCSCF: MultistateCalculations for F + H2O → HF + OH Reaction Paths
I. Introduction
Chemical reactions represent a complex network of elementary reaction steps,
each of which typically proceeds via highly reactive, open shell radical species.
Accurate characterization of these elementary chemical reactions from first principles
has been a long standing focus for both quantum chemists and dynamicists. This
represents a significant challenge on several fronts. As a necessary first step in thisprocess, one must calculate accurate adiabatic (i.e., Born-Oppenheimer) potential
energy surfaces as a function of fixed nuclear configurations. There are now many
examples of high-level ab initio potentials for simple atom + diatom chemical
reaction systems1, though this rapidly becomes a daunting task in higher (i.e.,3N - 6)
dimensionality. Once an adiabatic potential has been obtained, the next challenge is toperform exact quantum state-to-state dynamics calculations on such a surface. This is
an area in which there have been significant advances in the last decade, and it is now
more or less straightforward to obtain these calculations for triatomic systems
reacting on a single electronic surface. However, extension of these exact quantum
scattering methods to polyatomic reaction systems (N > 3) represents the current
state-of-the art for chemical dynamics and has thus far only been heroically achieved
for special 4 atom systems such as OH + H2 system with only one “heavy” (i.e., nonhydrogenic) atom2 3.

125
Unfortunately, state-to-state chemical reaction dynamics may not be so simply
described by quantum wave packet propagation on a single electronic potential
energy surface. As a result, quantum reactive scattering even for “simple” A + BC
triatomic systems can become substantially more complicated by the presence of nonadiabatic interactions between various Born-Oppenheimer surfaces. Indeed, there is a
growing body of experimental and theoretical evidence (with particular focus on
prototypical F + H24,5 and Cl + H26 atom abstraction systems) that non-adiabatic
effects can play a finite role in the reactive scattering dynamics. The presence of nonadiabatic interactions even in these relatively simple atom + diatom systems raises
obvious questions about the role of such nonadiabatic effects in atom + polyatom
systems. This highlights the crucial importance of high-level ab initio methods forcalculating reaction paths and potential energy surfaces for multiple electronically
excited states, which is of particular relevance to the present work.
A recent thrust in our group has been the use of high resolution UV lasers incrossed supersonic jets to study quantum state-resolved reaction dynamics forbenchmark prototypic 4 atom reactions. The focus of particular interest has been theatom + triatom system F(2P) + H2O → HF + OH (2Π). There are many aspects of thisreaction that make it particularly appropriate for combined experimental andtheoretical study. First of all, it is strongly exothermic (ΔE ≈17.6 kcal/mol7) andproceeds rapidly at thermal collision energies over a relatively low barrier (∼5kcal/mol), which makes it experimentally convenient for study under crossed
molecular beam conditions. Second, the products HF and OH can each be detectedwith full quantum state resolution, exploiting either LIF 8, or high resolution IR5,9,

126
techniques available that are in the group. Thirdly, from a theoretical perspective, thenumber of electrons (19) is small enough for obtaining accurate reaction paths,
potential energy surfaces, and nonadiabatic coupling matrix elements via high-level,
multireference ab initio methods. Most importantly, such a “heavy + light-heavylight” system offers the next level of challenge to quantum reactive scattering
methods, as well as eventually providing the experimental data with which tofacilitate detailed comparison with theory.
In the context of non-adiabatic reaction dynamics, the F(2P) + H2O systemrepresents an additional challenge for several reasons. First of all, in the asymptotic
regions of the potential, there are nearly degenerate surfaces corresponding to F(2P) +H2O reagent (3-fold degenerate) and HF + OH(2Π) product (2-fold degenerate) states.
In principle, such close proximity over long reaction path distances could promote
significant multiple state mixing in the entrance/exit channels, even for relatively
modest non-Born-Oppenheimer coupling matrix elements. Alternatively, strong nonadiabatic coupling could arise in the transition state region due to avoided crossings
between electronic surfaces. Indeed, the chemical bond rearrangement occurring near
the transition state reflects a rapid evolution of the wave function along the adiabatic
reaction path, which as elegantly shown by Butler and coworkers10, can often
dominate the shape and location of the transition state barrier. As also shown by
Allison et. al. 11, close avoided crossings lead to strong non-adiabatic behavior. The
complex nature of these curve crossings near transition states become especially richin systems with high electron affinity atoms, where the presence of charge transfer
states can be sufficiently low in energy to be important near the transition state127
region. This proves to be particularly true for the F + H2O system, where thetransition state region is dominated by avoided crossings between states correlating
asymptotically with ionic F- + H2O+ and HF+ + OH- surfaces. As this work will
illustrate, such strong avoidance of ground and excited state potential surfaces
generates special demands on calculating reaction paths via high-level, multireference
ab initio methods.
In order to investigate the role of non-adiabatic effects on reaction dynamics,
one must evaluate matrix elements of the nuclear kinetic energy operator, which
requires derivative calculations for ground and multiply excited electronic state wave
functions12. Such non-adiabatic matrix elements are calculated from the first andsecond derivatives of the wave function Ψi (Q1, Q2,…) for the ith electronic state,
where Ql represents motion along a specific intramolecular coordinate (l = 1, 3N - 6).
The full nuclear kinetic energy matrix elements are then obtained from thesederivative elements by summing with G matrix analysis13 and then added to the
traditional Born-Oppenheimer (BO) matrix to yield the full non-adiabatic
Hamiltonian12. The key point is that the reliable calculation of these derivative
coupling matrix elements therefore requires smoothly varying wave functions andadiabatic energies for each of these ground/excited states, providing crucially
important constraints on the fidelity of any high-level multireference calculations.
The conventional way to calculate high-level multireference SCF energies andwave functions for a series of ground and excited states is by a state-averaged
variational calculation14. In essence, this approach minimizes the weighted average of
MCSCF energies for each of the states in question,

128
⎛ n
⎞
E SA− MCSCF = ⎜ ∑ wi Ei ⎟⎝ i =1⎠(4.1)

where wi is a constant weight factor for state i. The advantage of such an approach isthat each of the electronic surfaces is treated equivalently in a variational sense,
which is obviously relevant to achieving a balanced wave function description in thelimit of asymptotically degenerate reactant (F + H2O) or product (HF + OH) states.
This can lead to small but significant shifts in the ground state energies. Whenvariationally optimizing the ground state(s), any included upper state will lead to aless correct description of the ground state(s).
The fundamental complication of such “statically weighted” ab initio methods
arises when there is a different electronic degeneracy for the product and productasymptotes, as in F + H2O → HF + OH, where the reactants (products) are triply(doubly) degenerate in the absence of spin-orbit coupling. In this case, a 3-state
MCSCF calculation will yield the correct exact 3-fold asymptotic degeneracy of theF(2P) reactants, which is broken by inclusion of higher excited states unless alldegenerate higher states are included with the same weight. Similarly, a 2-state
MCSCF calculation can correctly capture the asymptotic 2-fold degeneracy of theOH(2Π) products, which is again lifted by variationally averaging in higher state
contributions. A 3-state calculation will not describe the 2-fold degenerate OH(2Π)
ground state as well as a 2-state calculation. This situation is further exacerbated in
the F + H2O transition state region, where rapid avoided crossings of determinants in
the incomplete MCSCF wave function can result in abrupt discontinuities in both
energies and wave functions along the reaction path.

129
Simply stated, the challenge of such a state-averaged MCSCF calculation is toinclude a sufficient number of states in the chemically interesting transition region
and yet correctly capture the smooth evolution of the wave functions to asymptotic 3fold (reactant) and 2-fold (product) degeneracy, without introducing discontinuities inpotential surface. Toward this end, I propose a simple yet surprisingly robust method,“dynamically weighted” MCSCF (DW-MCSCF), whereby the state weighting is
explicitly varied as a continuous, damped function of the relative state energies.
⎛ n⎞
E DW − MCSCF = ⎜ ∑ wi ( Ei − E1 ) Ei ⎟⎝ i =1⎠(4.2)

For a broad range of damping functions, this yields excellent results andautomatically achieves the desired goals of i) smoothly varying energies and wave
functions along any locally smooth path (e.g. the reaction path) and ii) correct
degeneracy behavior out in the asymptotic reagent/product region by including only
the relevant states in each region.
The organization of this chapter is as follows. Section II describes details of
how the ab initio transition state and reaction path calculations are performed,
specifically focusing on the F + H2O system as a test case. Further details of the
dynamic weighting algorithm are described in Sec III, with Section IV presenting
results using dynamic weighting MCSCF along the F + H2O reaction path.
Concluding comments are summarized in Section V.

II. Ab Initio Calculations
All calculations are performed using the MOLPRO suite of ab initio
programs15. As my focus is on non-adiabatic dynamics, this requires a very high130level of theory capable of accurately describing multiple electronic states:
specifically, I perform complete active space self-consistent field calculations(CASSCF)14 using a full valence active space, followed by inclusion ofmultireference configuration interaction (MRCI)16 and Davidson correction(MRCI+Q)17. At this level of theory, the choice of basis set is somewhat constrainedby computational expense, particularly toward eventual calculation of a 4 atom
potential surface. For reference, a single point CASSCF+MRCI+Q energy calculation
with Dunning’s aug-cc-pVDZ18 basis takes approximately 3 CPU-hours on an AMD
Athlon 2.4 GHz. Timing/convergence tests for absolute energies of reagents andproducts, as well as harmonic zero point corrected reaction exothermicities for augcc-pVDZ, aug-cc-pVTZ, aug-cc-pVQZ, aug-cc-pV5Z basis sets are listed in Table
4.1. This series also permits extrapolation to the complete basis set (CBS) limit19,
based on a conventional 1/n3 plot with respect to the zeta basis set index. The timing
studies indicate a 2.6-fold increase in computational expense between aug-cc-pVDZ
and aug-cc-pVTZ, with steeper incremental penalties (4.2- and 5.2- fold) for each
additional zeta. However, acceptably good convergence (< 0.1 kcal/mole with respect
to CBS) is already achieved in predicted zero point corrected exothermicities at theaug-cc-pVTZ level, which dictates my choice.
Since my spin-orbit program cannot handle general contractions, a segmented
basis set was generated from the correlation-consistent triple-zeta basis set of
Dunning. For O and F, the exponents were taken from the aug-cc-pVTZ basis18, andthe first six s and three p-functions were contracted using the 1s and 2p contraction
coefficients, respectively. For H, the exponents of the cc-pVTZ basis18 were used and131the first three s-functions were contracted. Furthermore, the f-functions on O,F andthe d-functions on H were neglected. This resulted in a basis [6s4p3d] for O,F and[3s2p] for H (110 contractions). A number of basis contractions were tested for
calculational speed and accuracy before the final choice was made. In Table 4.2 arelisted a number of different contractions with convergence and timing results. These
data indicate the AVTZ[6s3p3d]/VTZ[3s2p] basis demonstrates similar levels of
convergence (0.34 kcal/mol above CBS) compared with a standard aug-cc-pVTZ
basis (0.12 kcal/mol above CBS), but with an additional 2.1-fold reduction in time. In
effect, the contracted aug-cc-pVTZ basis performs comparably for F + H2O to thestandard aug-cc-pVTZ basis, but at a computational speed only 20% slower than augcc-pVDZ. For perspective, also shown in Table 4.2 are results including spin-orbit
and zero-point energy corrections, which indicate residual discrepancies of 0.7
kcal/mol (4%) between the CBS limit and experimental determined exothermicities.
This discrepancy reflects correlation still not taken into account, and is consistent
with the residual 0.43 kcal/mol errors in exothermicity for the exhaustively studied F+ H2 potential surface of Stark and Werner20.
As a longer range goal, I am interested in high-level multistate ab initio
calculations for the full F + H2O PES, based on sampling approaches and Shephardinterpolation methods outlined by Collins et al21. The thrust of the current work,
however, is to illustrate the advantages of dynamically weighted MCSCF, which canbe adequately elucidated with calculations restricted along the reaction path. I start,
therefore, by obtaining the F + H2O transition state on the lowest electronic surface
by CASSCF geometry optimization programs in MOLPRO22. The crux of non-

132
adiabatic interactions for the F + H2O system arise from rapidly avoided crossings
between excited charge transfer states in the transition state region. To take thesestates into account in zeroth order, I search for this lowest saddle point in a traditional(i.e., non dynamically weighted) 5-state SA-MCSCF calculation; Figure 4.1 shows
the geometry of the transition state. Note that the transition state is significantly nonplanar (by nearly 60 degrees), and therefore no symmetry element is conserved
throughout

thereaction

coordinate.

This

most

likelyreflects

additional

hyperconjugative interaction between the radical p hole on F and the sp3 lone pairs on
the H2O subunit, which is further supported by the strongly non-collinear (140degree) F-H-O bond angle. For comparison, I have also converged the F + H2O
transition state geometry at MRCI+Q level, which yields similar parameters listed inTable 4.3.
From this transition state, a conventional reaction path is mapped in bothproduct and reagent directions by following the path of steepest decent for the lowest
state at the 5-state SA-MCSCF level23. When the gradient in the asymptotic region
becomes too small to follow, the reaction path is defined by simply extending the F –H2O or HF – OH bonds toward the reactants or products. Figure 4.2 shows theevolution of the nascently formed and broken bonds along this minimum energy path,in bond breaking ( ROH1 ) and making ( RFH1 ) coordinates; this curve will be defined as
the reaction coordinate for the rest of the paper. It is worth noting that, due to lack ofanalytic derivatives in the MOLPRO based MRCI code, the reaction coordinate has
been obtained only at the CASSCF level, which I next tackle with higher (MRCI+Q)
levels of theory. This will lead to additional errors (overestimation by about 0.9133kcal/mol) in the predicted barrier height. However, for the present purposes, the
CASSCF reaction path provides an excellent first approximation, which can later becompared with results from a more extensive MRCI+Q mapping of the transition
state region.
As a final step, spin-orbit interactions partially lift the degeneracy of theasymptotic F(2P) and OH(2P) states. To capture nonadiabatic effects in the F + H2O
reaction, therefore, spin-orbit interactions are calculated by diagonalizing the 6x6
Breit-Pauli Hamiltonian24 using the zero order non-spin-orbit wave functions as abasis. Although it is possible to use either SA-MCSCF or MRCI wave functions at
this point, I find that the SA-MCSCF wave functions reproduce the experimental
spin-orbit splitting as well as the MRCI wave functions, albeit with a factor of 100
reduction in calculation time. As a result, I use the SA-MCSCF wave functions to
calculate the spin-orbit interaction along the reaction path. Since the SA-MCSCF
wave functions are used as both reference configurations for MRCI calculations andas basis functions for calculating the spin-orbit interaction, a method is clearly needed
to calculate SA-MCSCF wave functions and eigenvalues that vary smoothly
throughout the entire PES. The next section outlines such a method for obtaining thisin a non-biased manner and is applicable to PES calculation in any dimensionality.

III. Dynamically Weighted Multiconfiguration Self ConsistentFields Method (DW-MCSCF)
The central issue with the conventional state averaged MCSCF approach to
the F + H2O → HF + OH reaction is easily identified. The fluorine atom has a triply
degenerate p-hole, thus requiring a 3-state SA-MCSCF calculation to correctly

134
describe the nature of the reactant channel. However, the lowest OH(2Π) productchannel is only doubly degenerate, which, if I attempt to describe with a 3-state
MCSCF calculation, will partially mix in character from the first excited state (5 eV),giving a less accurate description of the OH(2Π) ground state. By way of example,
Figures 4.3a and 4.3b show how the number of states in a SA-MCSCF calculation
affects the energies of levels for both reagents and products. In Figure 4.3a, the 3state calculation achieves the correct degeneracy for both the reactants and products,
but closer inspection indicates that the reaction path energies also exhibit strong
discontinuities in the transition state region. These are predominantly due to neglectof asymptotically higher energy charge transfer states that decrease rapidly with
distance and exhibit multiple avoided crossings with lower lying covalent states in thetransition state region. This can be improved by including more state averaging (forexample, 5-state results are shown in Fig 4.3b). This produces a much smoother set of
surfaces, but with a serious cost; equal variational weighting of these charge transfer
states now leads to significantly nondegeneracies for both reactant ground state (110cm-1).
To correctly maintain asymptotic degeneracies, one requires 3- (or 2-) statecalculations for the reactants (or products).

However, in the conventional SA-

MCSCF approach, this would necessitate an ad hoc “seam” region between the two
calculations, creating non-physical discontinuities in both the PES and wave
functions, and would thereby prohibit extraction of non-adiabatic coupling matrix
elements. An improved method by design would ideally yield a smooth transition
from a 3-state to a 2-state calculation so that the weights are correct asymptotically,

135
and yet would also include sufficient states in the transition state region to adequately
describe the avoided crossing dynamics. One idea is to set the weights with respect toa smoothly varying function of the nuclear coordinates. This is relatively easy todefine for a one-dimensional reaction path
(4.3)

wi = f i (s ) ,where wi is the weight of state i and if is a smoothly varying function of s, thereaction coordinate25.However, this becomes much harder for any higher dimensional
potential energy surface, and in addition, requires a priori knowledge of the numberof states important to the calculation in each region of configuration space.
As I care most about the lowest lying surface(s) correlating with reactants andproducts, a more physically motivated approach is to adjust weights based on energy
separation with respect to the ground state at a given geometry. Toward this end, I
define an energy dependent weighting function

wi
=
f (0)
=
lim f (ΔE ) →ΔE → ∞f (ΔE )
1 ,
0

(4.4)

where ΔE is the energy with respect to the ground state. By definition, the weights forasymptotically degenerate states will be identical (as required to achieve a proper
description of these states), and yet will appropriately decay to zero for sufficiently
highly excited states.
The algorithm is implemented iteratively, with initial estimates of energy
dependent weights chosen from a closely neighboring geometry (or a priori
knowledge of the electronic states in the region); the process continues until energies
are converged to a specific target value. For the present example, I chose 10-5 Hartree136(~0.006 kcal/mol or 2 cm-1), as it is about 1% of the asymptotic spin-orbit splittings.
Although the SA-MCSCF program is called more than once per point, the wholeprocess is actually quite efficient; calculations converge within ∼2-4 calls of the SAMCSCF program. This is because the weights change very little along the reaction
path, and the starting orbitals are already well optimized from the last SA-MCSCF
calculation. It is worth noting that even for bad starting guesses for the weights, e.g.
all weight set equal, only 1 or 2 more SA-MCSCF calls are needed for convergence tothe same threshold. In this application, the SA-MCSCF calculation is a 6-statecalculation with the weights of the first 6 states determined by the weighting function.
Throughout the entire reaction path the weight of the 6th state never exceeds
approximately 10-5, so increasing to a larger calculation will be of no benefit.
We next investigate the behavior of three possible weighting functions thatmatch the criteria in Eq. 4.3. An obvious first choice for a damping functional form
is f(ΔE) ∝ exp(-β (E - E0)), but this leads to slower iterative convergence due to a
non-zero slope in the exponential at E = E0. An improved functional form is a
Gaussian [f(ΔE) ∝ exp(-(β (E - E0))2)], which has a zero derivative at E = E0; thisensures that all nearly degenerate ground states receive equal weighting and therefore
the process converges efficiently (∼ 2 - 3 iterations). However, for Gaussian damping,
the relative weights decay too rapidly as near degeneracies begin to lift; this proves
undesirable for describing multiple levels near the transition state region accurately.
For example, this can result in state weighting coefficients changing nearly
discontinuously over relatively short reaction path distances, and thereby introducing
non-physical curvature in the calculated potential curves. A good compromise

137
between these two damping functional forms is empirically found to be f(ΔE)∝ sech2(-(β (E - E0)), which has both a zero derivative near E = E0 and yet exhibits amore gradual exponential vs. Gaussian drop off for large ΔE.
The dependence on choice of weighting function is illustrated in Fig 4.4a,
explicitly focusing on the ground state surface. Figure 4.4a displays DW-MCSCF
results for exponential, Gaussian, and sech2 weighting functions for a fixed value of
β-1 = 2eV. Noteworthy is the relatively weak sensitivity to choice of weighting
function, even near the transition state region. Also included on the plot are theenergies of the ground state for the 3- and 5-state constant weight methods. The
sensitivity (in shape and absolute energy of the ground state) to weighting function ismuch smaller than the sensitivity on number of states included in the calculation. On
closer inspection, one sees a small but systematic increase in variational energy forexponential weighting. This is due to the slower drop-off in excited state weighting
coefficients with increasing βΔE values, which therefore tends to include more statesat the expense of ground state accuracy. Conversely, considerably fewer differences
are evident between Gaussian and sech2 functional weighting, due to more equivalentweighting of states with βΔE ≈ 0.
An additional degree of freedom is the choice of β, which characterizes therange of energies weighted appreciably. If β is too large, then the calculations
effectively collapse to ∼1 (i.e., ground) state MCSCF, which is therefore unable toaccount for rapidly avoiding adiabatic curve crossings near the transition state.
Conversely, for small β, the weights do not decay fast enough with energy in theasymptotic reactant and product regions, which will not satisfy the requirement of a

138
3-state calculation in the reactant channel and 2-state in the product channel. The
sensitivity to choice of energy scale parameter β-1 is illustrated in Figure 4.4b for a
fixed sech2(βΔE) weighting function, with β-1 varying from 1eV to 4eV. Note again
the relatively weak dependence of the reaction path for small scale parameters (β-1 ≈ 2- 3 eV) , gradually degrading with increasing β-1 as more excited states are included.
However, also worth noting is the impact of too small a scale parameter, which canlead to distortions in the reaction path near the transition state (e.g. point P). Sucheffects occur because the first excited state energy is rapidly increasing near thetransition state, which for small β-1 implies an abrupt decrease in the number ofappreciably weighted excited states. More efficient iterative convergence properties
may be achieved by better tailored weighting functions for a given reaction system. In
practice, however, this simple first-order treatment already makes it straightforward
to obtain reaction paths relatively insensitive to scale parameter and functional
weighting, as described in the following section.

IV. Results
We next demonstrate how the DW-MCSCF procedure works when applied to
the F + H2O Æ HF + OH system. DW-MCSCF calculations with a sech2 weighting
function (β-1 = 3 eV) have been performed along the CASSCF reaction path. The
results for the 4 lowest states are summarized in Figure 4.5 and indicate several
features that are worth noting.
First of all, the potential curves highlight the presence of multiple barriers,
avoided crossings and a bound excited state well in the transition state region. Wave
function analysis indicates the dominant orbital configurations of the bound state well

139
to be F- + H2O+ and HF+ + OH- charge transfer states, on the reactant and productsides of the transition state, respectively. Though only adiabatic curves can be
rigorously defined, there are also clearly indications of a diabatic crossing between
the first excited F + H2O reactant state and the HF+ + OH- surface, as well as between
the F- + H2O+ charge transfer channel and the first excited HF + OH(2Π) productchannel. There are also indications of a second isolated curve crossing between thethird state and fifth state (not shown). These multiple crossings illustrate why suchhigh state averaging proved necessary (at least local to the transition state) and further
confirm expectations of strongly non-adiabatic dynamics. Second, despite thepresence of multiple crossings, the potential curves (and wave functions) can beobtained as a smooth function of the reaction path. Third, reactant and productchannels now achieve the correct degeneracy in the asymptotic limits, i.e., 3-fold and2-fold degenerate for F + H2O and HF + OH, respectively.
A quantitative display of dynamic weighting as a function of distance along
the reaction path is shown in Figure 4.6, with the adiabatic potentials reproduced in
the top panel. As desired, the figure clearly indicates a smooth evolution from i) a 3state calculation for triply degenerate F + H2O reactants, ii) through a dynamic
mixture near the transition state, into iii) a ∼2-state calculation for doubly degenerate
HF + OH products. Upon closer inspection of the product asymptotic region, there isa small but finite weighting of the third state (HF + OH(A2Σ)). This could be further
minimized by choice of β-1 scale parameter, but the effect is so small to not incur asignificant change in the OH(2Π) energies, and using β-1 produces the smoothest
reaction path.

140
To correctly account for spin-orbit effects in the F + H2O reaction, the BreitPauli Hamiltonian24 is diagonalized using the zero order non-spin-orbit wave
functions as a basis. As mentioned previously, the MCSCF wave functions reproduce
the experimental spin-orbit splitting as well as the higher level MRCI+Q wave
functions, albeit with a 100-fold reduction in computational expense. Using the DWMCSCF results as basis functions for the spin-orbit Hamiltonian, the 6x6 spin-orbit
matrix is therefore evaluated and diagonalized as a function of reaction pathgeometry. The asymptotic spin-orbit splittings in the reactant channel due to F(2P3/2,-11/2) are calculated to be 390 cm , i.e., in good agreement withthe experimental

values of 404.10 cm-1. Similarly, the spin-orbit splittings for product OH(2Π 1/2, 3/2)
are found to be 134 cm-1, in similarly good agreement with the 139.27 cm-1
experimental value.
As a final stage, results for the lowest three states from a full MRCI+Q
reaction path calculation with spin-orbit interactions are presented in Figure 4.7 andindicate several final points worth noting. First of all, the spin-orbit splitting is now
visible on the scale of the plot, and is in good agreement with experiment. Second, thereaction barrier for the ground state has now dropped to about 7 kcal/mol. Whencorrected for zero point effects at the reactant and transition state, this translates into abarrier around 5 kcal/mol, which is in reasonable agreement with experiment.
As a final note, these adiabatic correlation diagrams offer special relevance forfurther studies of non-adiabatic reaction dynamics. Specifically, the reactants startingout in the ground spin-orbit F(2P3/2) and crossing over the ground state barriercorrelate exclusively with OH (2Π3/2) products in the ground spin-orbit state. On the141other hand, the first excited spin-orbit state F(2P1/2) correlates adiabatically over a
very high barrier (∼60 kcal/mol) with the electronically excited OH (A2Σ) product,
and therefore should be non-reactive in the pure Born-Oppenheimer limit. Thisbehavior is quite similar to Born-Oppenheimer predictions for F/F* spin-orbit
reactivity for the well studied F + H2, 4,5and F + Cl2, 6systems. However, what is novel
about the set of F + H2O surfaces is the presence of a relatively high barrier (∼25
kcal/mol) for reactions adiabatically forming the first excited spin-orbit stateOH(2Π1/2). This suggests that at typical collision energies in a crossed molecularbeam, reaction dynamics in the purely adiabatic limit should produce exclusively
ground state OH (2Π3/2) and that the presence of any excited OH(2Π1/2) products must
therefore arise from non-adiabatic “surface hopping” dynamics. We are currently
investigating this prediction experimentally, taking advantage of high sensitivity laser
induced fluorescence (LIF) to probe the spin-orbit and rotational distributions of thenascent product OH. From a theoretical perspective, I will also use the smooth
adiabats calculated by the DW-MCSCF method as reference functions for high-level
MRCI+Q calculations and thereby compute the non-adiabatic coupling matrix
elements relevant for the F + H2O reaction system. Of special interest will be
determining the region of the potential where these non-adiabatic couplings dominate
and whether any surface hopping dynamics between OH (2Π3/2) and OH (2Π1/2)
surfaces are restricted to i) long range entrance or exit channel interactions26, in fact
are ii) more strongly localized in the transition state vicinity27, or are iii) a kinetic
effect where angular momentum of molecular rotation and spin-orbit states couple,
mixing the two states28.

142

V. Summary and Conclusion
DW-MCSCF is an addition to the standard SA-MCSCF calculation that
adjusts the weights of included states automatically to account for changes in theelectronic spectrum in different regions of a global PES. DW-MCSCF is useful inreactions where i) the degeneracy in the asymptotes is different and ii) there are anumber of low-lying states in the transition state region, because it will smoothly
adjust the weights among the regions of interest.
The DW-MCSCF algorithm is an approach where a starting guess for theweights is improved upon through iterations of a SA-MCSCF calculation.

The

process is quite efficient, and converges in only a small number of iterations. Since
the MRCI+Q calculations are orders of magnitude longer in computational time, it isworth the extra MCSCF iterations to ensure a smooth set of starting reference wave
functions and energies before taking the next steps.
DW-MCSCF is general and can be applied to calculate other PESs where
there are a number of important states varies throughout configuration space.
Moreover, DW-MCSCF can be easily applied to a PES of any number of dimensions
since the weighting function is only a function of the energy of the states (as opposed
to nuclear coordinates). The smoothly varying wave functions obtained by DWMCSCF can then be used as reference functions for MRCI+Q calculations, and nonadiabatic calculations which require derivatives of the wave function.
The example reaction F(2P) + H2O → HF + OH(2Π) is prototypical of thereactions that are good candidates for DW-MCSCF:i) the reactants are 3-folddegenerate and products are 2-fold, and ii) there are avoided crossings involving low-

143
lying charge transfer states in the transition state region. Non-adiabatic effects arelikely very important in this reaction, so it is imperative that smoothly varying multistate wave functions are calculated, as the non-adiabatic corrections are functions ofderivatives of the wave function. The full 6-D PES of this reaction will be calculatedusing DW-MCSCF followed by MRCI+Q along with spin-orbit and non-adiabatic
corrections to the lowest 3 surfaces.

144

TablesTable 4.1. Product, reactant energies and reaction exothermicities as a function of level of theory and basis set size.

Basis
AVDZ
AVTZ
AVQZ
AV5Z
CBS

F + H2O (Hartree)
CASSCFMRCI
MRCI+Q0a -0.325174 -0.349268-0.044082 -0.459217 -0.490918-0.057301 -0.502245 -0.536145-0.060628 -0.516598 -0.551124-0.064583 -0.525946 -0.561136

HF + OH (Hartree)
CASSCFMRCI
MRCI+Q-0.009777 -0.346794 -0.372539-0.054397 -0.481240 -0.514744-0.067061 -0.524344 -0.560164-0.070337 -0.538542 -0.574991-0.074483 -0.548033 -0.585155

CASSCF
-6.13
-6.47
-6.12
-6.09
-6.21

a) Absolute energies stated relative to AVDZ reactant energy: -175.471563 Hartree.

ΔE (kcal/mol)
MRCI
MRCI+Q-13.57-14.60-13.82-14.95-13.87-15.07-13.77-14.98-13.86-15.07

rel. time
1.0
2.6
11.0
57.3
-

145
Table 4.2. Product, reactant energiesa and reaction exothermicities as function of optimizing the AVTZ basis set

Basis
AVTZ[5s4p3d2f]/AVTZ[4s3p2d]
AVTZ[6s4p3d2f]/AVTZ[3s3p2d]
AVTZ[6s4p3d]/VTZ[3s2p1d]
AVTZ[6s4p3d2f]/AVTZ[3s3p]
AVTZ[6s4p3d]/VTZ[3s2p]+Spin-Orbit+ZPE
ExperimentF + H2O MRCI+Q (Hartree)a-0.490918-0.475371-0.474516-0.424399-0.423947-0.424561-0.386702-

HF + OH MRCI+Q (Hartree)a-0.514744-0.499164-0.498290-0.447828-0.447425-0.447633-0.429889-

a) Absolute energies relative to AVDZ reactant energy: -175.471563 Hartree.

ΔE (kcal/mol)-14.95-14.93-14.92-14.70-14.73-14.48-16.59-17.61

rel. time
2.6
2.6
2.0
1.4
1.2
-

146
Table 4.3. The geometry of the transition state at the SA-MCSCF and MRCI+Q level. The lengths and angles are definedin Figure 4.1. Bond lengths are in angstroms, angles in degrees. The basis set used is AVTZ[6s4p3d]/VTZ[3s2p] in a fullvalence active space.

SA-MCSCF
MRCI+Q

RHF1
1.218
1.350

ROH1
1.093
1.031

ROH2
0.974
0.971

Α
101.378102.551β139.771118.407Γ
58.15869.552

147

FiguresFigure 4.1: The internal coordinates used for specifying F + H2O geometries. The
geometry is given for the 5-state MCSCF transition state.

148

3.0

2.5

1

Å

2.0

OH
1.5
‡

1.0
1.0

1.5

2.0

2.5

3.0

R FH /Å
1

Figure 4.2: A 2D projection (RFH1 and ROH1 bond distances) of the full 6D reaction
path calculated at the 5-state MCSCF level. This 6D reaction path is used throughout
the paper.

149

4

3-State

2

0.02
0.00
-0.02

-0.26
-0.28
-0.30

0
4

5-State

2

0.74
0.72
0.70

0.42
0.40
0.38

0
4

Dynamic Weight

2

0.02
0.00
-0.02

-0.26
-0.28
-0.30

0

Reaction Coordinate
Figure 4.3: Calculated ground and excited state energies along the F + H2O reaction
coordinate for various levels of theory, a) 3-state MCSCF, b) 5-state MCSCF, c)
Dynamically weighted MCSCF.

The inserts at reactant and product asymptotes

indicate a blow up of the energy scale to highlight deviations from perfect 3-fold(reactant) and 2-fold (product) degeneracies. Note that the dynamic weighting method
produces the smoothest surface, and also reproduces the correct degeneracy at both
asymptotes of the reactants.

150

1.2

3-State5-State
0.9
0.6

exp
Gaussian
2
sech

0.3
0.0

(a)

Reaction Coordinate/

1.2

β -1=2eVP

0.9

β =3eV-1β =4eV-1

0.6

β =1eV
-1

0.3
0.0-0.3 (b)Reaction Coordinate
Figure 4.4: Dependence of ground state reaction path energy on weighting function in
the transition state region. a) Exponential, Gaussian or sech2 functional form. Notethat the choice of the functional form has less of an effect on the surface than thenumber of states included in a traditional SA-MCSCF. b) Choice of scale parameterβ.

If β-1 is too small, small changes in energies in excited states can cause

discontinuities in the ground state (e.g. point P), if β-1 is too large, the weights do notdecay to 0 in the asymptotes for the excited states. Moderate β-1 parameters (i.e.,23eV) produce smooth curves and converge to ground-state-only asymptotes.

151

6

-

F +H 2O++

HF +OH-

5

HF+OH(2Σ)

4
3
2
1
0

F(2P)+H 2O

HF+OH(2Π)

Reaction Coordinate
Figure 4.5: The lowest 4 states of the smooth DW-MCSCF reaction path. Theenergies along the reaction path, not only vary smoothly, but the degeneracies arecorrectly reproduced in the asymptotes. Also note that there are a number of statesinteracting in the transition state region.

152

1.0

0.8

State 3

State 5
State 4

0.6

State 20.4

0.2

State 1

0.0

Reaction Coordinate
Figure 4.6: Converged weights (normalized to sum to 1) for the DW-MCSCF with β-1= 3. The energies along the reaction path are reproduced at the top to help guide the

1 1 1
1 1
eye. The weights smoothly vary from ( : : ) to approximately ( : :0). During
3 3 3
2 2
the cross-over through the transition state a fourth and fifth state are briefly
introduced. DW-MCSCF automatically included the relevant states in the different
regions of the PES.

153

4
HF+OH(2Σ)

3

2

1
F(2P1/2)+H2O
0

HF+OH(2Π 1/2)

F(2P3/2)+H2O

HF+OH(2Π 3/2)

-1

Reaction Coordinate
Figure 4.7: MRCI+Q with spin-orbit interactions calculated at the CASSCF reaction
path geometries. The smoothly varying DW-MCSCF wave functions are used as
reference configurations for MRCI+Q calculations.

Note that the lowest barrier(when zero-point energy is included) is experimentally accessible in crossed
supersonic jets, but the other barriers are far out of reach.

154

References

1

J. M. Bowman and G. C. Schatz, Annual Review of Physical Chemistry 46,
169 (1995); S. C. Althorpe and D. C. Clary, Annual Review of Physical
Chemistry 54, 493 (2003).

2

D. H. Zhang and J. Z. H. Zhang, Journal of Chemical Physics 99 (7), 5615
(1993); D. H. Zhang and J. Z. H. Zhang, Journal of Chemical Physics 101 (2),1146 (1994); D. H. Zhang and J. Z. H. Zhang, Journal of Chemical Physics

100 (4), 2697 (1994); Y. C. Zhang, D. S. Zhang, W. Li, Q. G. Zhang, D. Y.
Wang, D. H. Zhang, and J. Z. H. Zhang, Journal of Physical Chemistry 99(46), 16824 (1995); D. H. Zhang, J. Z. H. Zhang, Y. C. Zhang, D. Y. Wang,
and Q. G. Zhang, Journal of Chemical Physics 102 (19), 7400 (1995); D. H.
Zhang and J. Z. H. Zhang, Chemical Physics Letters 232 (4), 370 (1995); D.
H. Zhang and J. C. Light, Journal of Chemical Physics 105 (3), 1291 (1996);
D. H. Zhang and J. C. Light, Journal of Chemical Physics 104 (12), 4544
(1996); W. Zhu, J. Q. Dai, J. Z. H. Zhang, and D. H. Zhang, Journal ofChemical Physics 105 (11), 4881 (1996); Y. C. Zhang, Y. B. Zhang, D. S.
Zhang, Q. G. Zhang, D. H. Zhang, and J. Z. H. Zhang, Chinese Science
Bulletin 42 (2), 116 (1997); D. H. Zhang and J. C. Light, Journal of the
Chemical Society-Faraday Transactions 93 (5), 691 (1997); D. H. Zhang and
J. C. Light, Journal of Chemical Physics 106 (2), 551 (1997); J. C. Light andD. H. Zhang, Abstracts of Papers of the American Chemical Society 215,
U201 (1998); D. H. Zhang, J. C. Light, and S. Y. Lee, Journal of Chemical

155
Physics 109 (1), 79 (1998); R. P. A. Bettens, M. A. Collins, M. J. T. Jordan,
and D. H. Zhang, Journal of Chemical Physics 112 (23), 10162 (2000); D. H.
Zhang, M. A. Collins, and S. Y. Lee, Science 290 (5493), 961 (2000); Y. M.
Li, M. L. Wang, J. Z. H. Zhang, and D. H. Zhang, Journal of Chemical
Physics 114 (16), 7013 (2001); M. H. Yang, D. H. Zhang, M. A. Collins, andS. Y. Lee, Journal of Chemical Physics 115 (1), 174 (2001); M. H. Yang, D.
H. Zhang, M. A. Collins, and S. Y. Lee, Journal of Chemical Physics 114(11), 4759 (2001).
3

M. Brouard, I. Burak, S. Marinakis, D. Minayev, P. O'Keeffe, C. Vallance, F.
J. Aoiz, L. Banares, J. F. Castillo, D. H. Zhang, D. Xie, M. Yang, S. Y. Lee,and M. A. Collins, Physical Review Letters 90 (9) (2003); M. Brouard, I.
Burak, D. Minayev, P. O'Keeffe, C. Vallance, F. J. Aoiz, L. Banares, J. F.
Castillo, D. H. Zhang, and M. A. Collins, Journal of Chemical Physics 118(3), 1162 (2003).

4

M. H. Alexander, D. E. Manolopoulos, and H. J. Werner, Journal of Chemical
Physics 113 (24), 11084 (2000); M. H. Alexander, Abstracts of Papers of theAmerican Chemical Society 218, U306 (1999); M. H. Alexander, H. J.
Werner, and D. E. Manolopoulos, Journal of Chemical Physics 109 (14), 5710
(1998); S. A. Nizkorodov, W. W. Harper, and D. J. Nesbitt, Faraday
Discussions (113), 107 (1999); J. F. Castillo, D. E. Manolopoulos, K. Stark,and H. J. Werner, Journal of Chemical Physics 104 (17), 6531 (1996); F. J.
Aoiz, L. Banares, B. MartinezHaya, J. F. Castillo, D. E. Manolopoulos, K.

156
Stark, and H. J. Werner, Journal of Physical Chemistry A 101 (36), 6403
(1997).
5

S. A. Nizkorodov, W. W. Harper, W. B. Chapman, B. W. Blackmon, and D. J.
Nesbitt, Journal of Chemical Physics 111 (18), 8404 (1999).

6

M. H. Alexander, G. Capecchi, and H. J. Werner, Science 296 (5568), 715(2002); M. Alagia, N. Balucani, L. Cartechini, P. Casavecchia, G. G. Volpi, F.
J. Aoiz, L. Banares, T. C. Allison, S. L. Mielke, and D. G. Truhlar, Physical
Chemistry Chemical Physics 2 (4), 599 (2000).

7

S. A. Harich, D. W. H. Hwang, X. F. Yang, J. J. Lin, X. M. Yang, and R. N.
Dixon, Journal of Chemical Physics 113 (22), 10073 (2000); W. T. Zemke,
W. C. Stwalley, J. A. Coxon, and P. G. Hajigeorgiou, Chemical Physics
Letters 177 (4-5), 412 (1991).

8

O. Votava, D. F. Plusquellic, and D. J. Nesbitt, Journal of Chemical Physics

110 (17), 8564 (1999); J. R. Fair, O. Votava, and D. J. Nesbitt, Journal ofChemical Physics 108 (1), 72 (1998); D. F. Plusquellic, O. Votava, and D. J.
Nesbitt, Journal of Chemical Physics 109 (16), 6631 (1998); D. F. Plusquellic,
O. Votava, and D. J. Nesbitt, Journal of Chemical Physics 107 (16), 6123
(1997); O. Votava, J. R. Fair, D. F. Plusquellic, E. Riedle, and D. J. Nesbitt,
Journal of Chemical Physics 107 (21), 8854 (1997).
9

W. W. Harper, S. A. Nizkorodov, and D. J. Nesbitt, Journal of Chemical
Physics 116 (13), 5622 (2002); W. W. Harper, S. A. Nizkorodov, and D. J.
Nesbitt, Chemical Physics Letters 335 (5-6), 381 (2001); W. W. Harper, S. A.
Nizkorodov, and D. J. Nesbitt, Journal of Chemical Physics 113 (9), 3670

157
(2000); W. B. Chapman, B. W. Blackmon, S. Nizkorodov, and D. J. Nesbitt,
Journal of Chemical Physics 109 (21), 9306 (1998).
10

T. L. Meyers, D. C. Kitchen, B. Hu, and L. J. Butler, Journal of Chemical
Physics 105 (7), 2948 (1996); G. C. G. Waschewsky, P. W. Kash, T. L.
Myers, D. C. Kitchen, and L. J. Butler, Journal of the Chemical SocietyFaraday Transactions 90 (12), 1581 (1994); P. W. Kash, G. C. G.
Waschewsky, L. J. Butler, and M. M. Francl, Journal of Chemical Physics 99(6), 4479 (1993).

11

T. C. Allison, S. L. Mielke, D. W. Schwenke, and D. G. Truhlar, Journal of
the Chemical Society-Faraday Transactions 93 (5), 825 (1997).

12

J. C. Tully, in Dynamics of Molecular Collisions: Part B, edited by W. H.
Miller (Plenum Press, New York, 1976).

13

J. H. Frederick and C. Woywod, Journal of Chemical Physics 111 (16), 7255
(1999); B. R. Johnson and W. P. Reinhardt, Journal of Chemical Physics 85(8), 4538 (1986).

14

H. J. Werner and P. J. Knowles, Journal of Chemical Physics 82 (11), 5053
(1985); P. J. Knowles and H. J. Werner, Chemical Physics Letters 115 (3),259 (1985).

15

R. D. Amos, A. Bernhardsson, A. Berning, P. Celani, D. L. Cooper, M. J. O.
Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, P. J. Knowles, T.
Korona, R. Lindh, A. W. Lloyd, S. J. McNicholas, F. R. Manby, W. Meyer,
M. E. Mura, A. Nicklass, P. Palmieri, R. Pitzer, G. Rauhut, M. Schutz, U.
Schumann, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and H. J.

158
Werner, MOLPRO, a package of ab initio programs designed by

H. J.

Wernerand P. J. Knowles, version 2002.1 (Birmingham, UK, 2002).
16

H. J. Werner and P. J. Knowles, Journal of Chemical Physics 89 (9), 5803
(1988); P. J. Knowles and H. J. Werner, Chemical Physics Letters 145 (6),514 (1988); P. J. Knowles and H. J. Werner, Theoretica Chimica Acta 84 (12), 95 (1992).

17

S. R. Langhoff and E. R. Davidson, International Journal of Quantum
Chemistry 8 (1), 61 (1974).

18

T. H. Dunning, Journal of Chemical Physics 90 (2), 1007 (1989).

19

T. Helgaker, P. Jorgensen, and J. Olsen, Molecular Electronic-Structure
Theory. (John Wiley & Sons, LTD, New York, 2000).

20

K. Stark and H. J. Werner, Journal of Chemical Physics 104 (17), 6515
(1996).

21

M. A. Collins, Abstracts of Papers of the American Chemical Society 223,
U493 (2002); K. C. Thompson, M. J. T. Jordan, and M. A. Collins, Journal ofChemical Physics 108 (20), 8302 (1998); K. C. Thompson, M. J. T. Jordan,and M. A. Collins, Journal of Chemical Physics 108 (2), 564 (1998); K. C.
Thompson and M. A. Collins, Journal of the Chemical Society-Faraday
Transactions 93 (5), 871 (1997); M. J. T. Jordan, K. C. Thompson, and M. A.
Collins, Journal of Chemical Physics 103 (22), 9669 (1995); M. J. T. Jordan,
K. C. Thompson, and M. A. Collins, Journal of Chemical Physics 102 (14),159
5647 (1995); J. Ischtwan and M. A. Collins, Journal of Chemical Physics 100(11), 8080 (1994).
22

F. Eckert, P. Pulay, and H. J. Werner, Journal of Computational Chemistry 18(12), 1473 (1997).

23

F. Eckert and H. J. Werner, Theoretical Chemistry Accounts 100 (1-4), 21(1998).

24

A. Berning, M. Schweizer, H. J. Werner, P. J. Knowles, and P. Palmieri,
Molecular Physics 98 (21), 1823 (2000).

25

K. Fukui, Accounts of Chemical Research 14 (12), 363 (1981); K. Yamashita,
T. Yamabe, and K. Fukui, Chemical Physics Letters 84 (1), 123 (1981).

26

A. R. Offer and G. G. Balintkurti, Journal of Chemical Physics 101 (12),10416 (1994).

27

T. L. Myers, N. R. Forde, B. Hu, D. C. Kitchen, and L. J. Butler, Journal ofChemical Physics 107 (14), 5361 (1997).

28

R. N. Zare, in Angular Momentum (John Wiley & Sons, New York, 1988), pp.
302.

160

Chapter 5:

Multireference configuration interaction

calculations for the F(2P) + HCl → HF + Cl(2P) reaction:
A correlation scaled ground state (12A’) potential
energy surface
I. Introduction
A first principles theoretical understanding of elementary reaction dynamics at
the quantum state resolved level has remained a major challenge to the chemical
physics community over the last several decades. As a result of this intense interest,
there have been impressive advances in both ab initio development of accurate
potential energy surfaces, as well as methods for numerically exact quantum
dynamical calculations on these surfaces. One system that has represented awatershed for reaction dynamics has been the F(2P) + H2 Æ HF + H(2S) reaction, forwhich there now exists an exceptionally high level potential surface1,2 and numerous
theoretical3-8 and experimental9-20 papers on state-to-state reactive scattering. Of
particular dynamical interest has been the subtly but importantly different “heavy +light-heavy” isotopic variation on this reaction system, F(2P) + HD Æ HF + D(2S),for which long lived scattering resonances in the transition state region have beenboth predicted21-31 and experimentally32 observed. Although there are many
contributing factors, inspection of the transition state wave functions offers a simple
physical picture of such resonances, classically corresponding to highly excited H
atom motion in a quasi bound HF(v = 3)--D state. In essence, this state executes rapid161
H atom motion between the F and D atoms, which eventually predissociates into thecontinuum of states corresponding to HF(v = 2,J) + D. Such efforts have further
served to elucidate the considerable complexity and dynamical richness present evenin simple H atom transfer reactions and continue to offer challenges for firstprinciples theoretical calculations in larger systems.
The clear observation of transition state resonances in F + HD suggests that
there may be similar resonance dynamical effects in the corresponding F(2P) +hydrogen halide systems. Of particular experimental and theoretical interest is theF(2P) + HCl Æ HF + Cl(2P) system, the high-level ab initio potential surface which
represents the focus of this paper. This reaction has a nearly identical exothermicity(ΔEHF+Cl = -33.1 kcal/mol vs. ΔEHF+ H = -31.2 kcal/mol), a thermally accessible
barrier height, and yet is still within reach of modern high-level ab initio efforts formapping ground and excited potential surfaces out in full dimensionality. In thiswork, I present a new benchmark ground state potential energy surface (PES) for theF(2P) + HCl Æ HF + Cl(2P) reaction, to accompany experimental direct IR absorption
and ion imaging studies currently in progress. In addition to being a prototypicallyheavy-light-heavy system, the F(2P) + HCl Æ HF + Cl(2P) reaction is also importantin understanding hydrogen-halide chemical lasers33 and provides novel opportunities
to study non-adiabatic dynamics on multiple electronic surfaces.
Previous theoretical studies on this surface have included quasiclassicaltrajectory (QCT) calculations performed on an empirical LEPS surface34, QCT35andfull dimensional time-dependant wavepacket 36 calculations obtained on a newer ab
initio (PUMP2/6-311G(3d2f,3p2d)) surface37, and a LEPS surface created from162
recent experimental kinetics information38. All previous calculations have focused
specifically on the ground (12A’) electronic state, though study of other F atom
abstraction systems has revealed evidence for avoided crossings, charge transfer
states, seams of conical intersection, etc. For example, it has been shown that there issignificant non-adiabatic behavior for F(2P) + H2 HF + H(2S)

in the entrance

channel3,23, and there is also direct experimental evidence in 4 atom systems(i.e.,(2P)F + H2O Æ HF + OH(2Π) that exhibit strong non-adiabatic behavior in boththe entrance and exit channels39. It is likely that similar effects are present in F(2P) +HCl Æ HF + Cl(2P), the dynamics of which clearly necessitate a high-level
multireference electronic treatment. In this work, I present a new ground electronic
state (12A’) PES for the F(2P) + HCl reaction, based on high-level multireference
calculations and representing a first step toward multiple fitted surfaces with nonadiabatic coupling, spin-orbit interactions explicitly included.
The F(2P) + HCl Æ HF + Cl(2P) reaction presents many of the samechallenges seen earlier in the F(2P) + H2O Æ HF + OH(2Π) reaction40, namely thepresence of highly electronegative atoms significantly lower the charge transfer (CT)states (F- + HCl+(2Π) Æ HF+(2Π) + Cl-) in the transition state region. This leads tonon-adiabatic coupling and avoided crossing behavior between the 12A’ and 22A’,
which requires multireference calculations to describe accurately. In light of these
considerations, ab initio calculations for this PES are obtained via internallycontracted multireference configuration interaction (MRCI)41-43 using reference statestaken from a dynamically weighed multiconfigurational self-consistent field (DWMCSCF)40 calculation. In addition, extrapolation to the complete basis set (CBS)163
limit44 and correlation energy scaling45,46 are used to further refine the PES shape andenergetics accurately.
As a further source of dynamical richness, the F(2P) + HCl Æ HF + Cl(2P)
surface is also expected to exhibit conical intersections in both entrance and exit
channels. Indeed, this can be anticipated from rather elementary chemistry
considerations, as described by Hutson and coworkers47-49 for open shell radical XHX complexes. Specifically, well outside the region of chemical reaction, the surface
will be asymptotically dominated by dipole-quadrupole and dipole-induced dipole
interactions. This results in a state energy ordering of E(Σ) > E(Πx) = E(Πy) for the
linear approach geometry. In the vicinity of the transition state, however, the lowest
energy collinear surface arises from the unfilled p-orbital of the bare F (or Cl) atom inan Σ configuration with respect to the HCl reactant (or HF product), i.e., yielding areverse energy ordering of E(Σ) < E(Πx) = E(Πy). Since the symmetries of thecoplanar Πx and Σ surfaces are reduced to non-crossing 12A’ and 22A’ surfaces forany non-collinear geometry, this necessarily implies the presence of a 3N-5 = 1D
“seam” of conical intersections. Of special interest for F + HCl reaction dynamics,
these conical intersection seams are both energetically low enough (Ecrossing = 3-4
kcal/mol) with respect to collision energy (Ecom = 5 kcal/mol) as well as the F(2P) +HCl transition state (ETST ≈ 3.8 kcal/mol) that they are likely play an important role.
The organization of this chapter is as follows: Section II describes the detailsof the ab initio calculations specifically focusing on the 12A’ electronic state. In
Section III, I describe the benchmark calculations, with emphasis on the CBS
extrapolation and correlation energy scaling that is used to provide an accurate PES.

164
Section IV presents the details of the analytical fit to the ab initio points. In Section
V I discuss the interesting aspects of this PES followed by a summary and concluding
remarks in Section VI.

II. Details of the Ab Initio Calculations
All calculations are performed using the MOLPRO suite of ab initio
programs50. Since I explicitly need to capture non-adiabatic as well as single surface
effects, electronically excited states of the F(2P) + HCl Æ HF + Cl(2P) reaction are
required. Thus, the ab initio methods utilize multireference methods to obtainwavefunctions for non-adiabatic calculations and characterization of the conicalintersections. In addition, other 3- and 4-atom systems40 involving F have indicatedlarge multireference character in the ground state wavefunctions.. In light of theserestrictions and past success with the dynamically weighted multiconfigurational selfconsistent field (DW-MCSCF)40 method, the natural choice for this reaction is to use
DW-MCSCF along with internally contracted multireference CI (MRCI) 41-43 and themultireference Davidson correction (MRCI+Q)51-53.

Two separate MRCI+Q

calculations are done for each symmetry: one for the 12A’ and 22A’ states and one forthe 12A” state. In this work, I restrict the focus on the ground electronic (12A’)
surface; a more detailed description of all three adiabatic surfaces, as well as the
diabatic representation and non-adiabatic coupling between these surfaces will bepresented elsewhere.54
Residual ab initio error due to incomplete basis sets can be corrected byextrapolating to the complete basis set (CBS) limit44 using the standard aug-cc-pVnZ(n = 2,3,4) basis sets of Woon and Dunning55-57. Others46,58 have shown that an

165
additional tight (large zeta) d function to the Cl basis set improves the accuracy, butbenchmark tests (vide infra) show no significant improvement in the F(2P) + HCl Æ
HF + Cl(2P) system, as calibrated against the precise experimental reaction
exothermicity. I therefore elect to use the standard aug-cc-pVnZ (AVnZ, n = 2-4)
basis functions and extrapolate to the CBS limit. Based on the benchmark studies ofPeterson, this extrapolation is obtained via the three parameter fit44,46
2

E (n) = ECBS + Be − ( n −1) + Ce − ( n −1) .

(5.1)

This function yields an excellent extrapolation to the CBS limit (see Figure 5.1) as
explicitly confirmed by additional calculations for reagent and product asymptotic
geometries with basis set extensions up to n = 5.
The MRCI+Q calculations are performed using a full valence complete active
space (CAS)59 reference function with 15 active orbitals (12A’,3A”), 6 of which(5A’,1A”) are uncorrelated.

The orbitals are obtained by a state-averaged

multiconfiguration SCF (SA-MCSCF)60,61 calculation with 15 active orbitals (12A’,
3A”) with 7 (6A’, 1A”) closed to excitations. A total of 6 states are then included in
the SA-MCSCF calculation (4A’, 2A”), with the weights determined by the
dynamically weighted-MCSCF algorithm.

In order to account for the triply

degenerate p-hole in both the F and Cl atoms at the asymptote, the SA-MCSCF
calculations include all 3 states (2A’, 1A”); however, low-lying charge transfer (CT)states (e.g. F- + HCl+(2Π) Æ HF+(2Π) + Cl-) require the inclusion of 2 more A” and 1
more A’ state (a total of 6 states) to fully describe the transition state region. As
previously demonstrated40, root flipping as the CT surfaces cross other states along
the reaction path makes it extremely challenging to obtain a smooth accurate PES for166
a fixed set of SA-MCSCF weights. In DW-MCSCF, on the other hand, the weights
for each state are chosen as a function of its energy with respect to the ground stateenergy at each point in configuration space. The function for the weight based on theenergy difference from the ground state iswi ∝ sech 2 (− β (Ei − E 0 ))(5.2)

The decay coefficient (β−1) for DW-MCSCF is empirically optimized to be 3eV,though the results prove relatively insensitive to this choice. Most importantly, thismethod permits one to eliminate spurious discontinuities from root flipping andthereby obtain a smooth reaction path and surface while accurately describing theasymptotes and the transition state. The orbitals obtained from the DW-MCSCF
calculation are used as reference states for the internally contracted MRCI
calculation.

III. Benchmark Calculations
A. Exothermicity
By way of calibration, I have performed a number of single point calculationsat product and reagent asymptotes (F + HCl and HF + Cl) to benchmark how well
MRCI+Q predicts the reaction exothermicity, ΔErxn, with respect to previous lower
level methods. Table 5.1 contains a summary of reaction exothermicities calculatedby these methods, each for an AVnZ basis set and extrapolated to the CBS limit. To
facilitate most the direct comparison with experiment, the energies in Table 5.1 havebeen corrected for experimental zero-point (5.98 kcal/mol and 4.31 kcal/mol for HF
and HCl, respectively) and spin-orbit splittings (1.16 kcal/mol and 2.52 kcal/mol for

167
F and Cl), which are not included the ab initio calculations reported here. The resultsin Table 5.1 indicate rapid convergence with respect to basis set size (ΔE ≈ 0.06
kcal/mol between AVQZ and the CBS limit values). Indeed, Figure 5.1 clearly showsthat energies for both reactant and product asymptotes are independently converging
to the CBS limit. By way of confirmation, the dotted line in Figure 5.1 represents the
CBS extrapolation for AVnZ basis sets with n = 2 - 4, which in fact very closely
reproduces the single point calculated values for AV5Z.
It is important to note that although the calculated reaction exothermicity
converges nicely to ∼-31.0 kcal/mol, the true exothermicity (dashed line) is in fact ∼2kcal/mol lower (-33.06 kcal/mol62,63).

This is not a limitation in the CBS

extrapolation but rather the size of the “complete” active space. This can be readily
seen in the incremental shifts in ΔE as a function of basis set size (inset in Figure 5.1),
which indicate the CBS extrapolation converging upward to a value higher than theexperimental energy. This arises from incomplete recovery of the full correlationenergy at the MRCI+Q level and is reminiscent of the comparable (∼ 0.4 kcal/mol)
exothermicity discrepancies from benchmark calculations for the F + H2 surface.
Indeed, as demonstrated by Stark and Werner for F + H21, 3pF orbitals must be
included in the active space to get results accurate to < 0.1 kcal/mol. For the F + HCl
reaction, however, including 3pF orbitals would energetically require also including
all 5 of the 3dCl orbitals, increasing the number of occupied orbitals from 15 (12A’,
3A”) to 24 (18A’, 6A”), and making the calculation prohibitively expensive. By way
of example, single point AVTZ calculations at reagent and product geometries withthis enhanced active space decrease the exothermicity error down to ∼0.5 kcal/mol;

168
however, even this requires more than 2 days on a single 2.4GHz processor, i.e., toolong to map out an entire 3D PES. I address below how to systematically compensate
for these modest residual discrepancies in correlation energy, but for the moment,
simply note that the calculations are quite well converged with respect to basis set by
CBS extrapolation for n = [2 - 4].
Building on ideas by Peterson46 and Truhlar45, I employ a simple global
correlation energy scaling procedure calibrated to match the experimentally known
reaction exothermicity. Specifically, I uniformly scale the correlation energy by aconstant factor near unity that reproduces the exothermicity exactly while retaining
the shape of the PES and not introducing any discontinuities.

First, I definecorrelation energy as the difference between the MRCI+Q energy and the MCSCF
energy:E ( R) CORR = E ( R) MRCI + Q − E ( R) MCSCF ,(5.3)

where R is a vector of the intermolecular coordinates. The scaled energy is just thesum of the MCSCF energy plus the correlation energy scaled by a factor, γ:
E ( R) SCALED = E ( R) MCSCF + γE ( R) CORR .

(5.4)

In the limit of γ = 0, the scaled energy is just the MCSCF energy (i.e., no extra
correlation energy is recovered), when γ = 1 the scaled energy is the MRCI+Q
energy (i.e., all correlation energy is recovered). Allowing γ to be greater than unity
permits recovery of additional correlation energy, in effect, approximating an increasein active space.

The value for γ is empirically chosen to match the reaction

exothermicity, i.e.,E ( prod ) SCALED − E (reac) SCALED = ΔE EXPERIMENT .

(5.5)

169
The scaling factors (γ) for AVnZ (n = 2 - 4), as well as the CBS extrapolation,
monotonically decrease as the basis set becomes more complete. Indeed, at the CBS
limit it is close to unity (γ = 1.027) for the A’ symmetry, i.e., one only needs torecover an extra 2.7% of the correlation energy to reproduce the experiment. It is
worth noting that the scaling values for A’ (γ = 1.027) and A” (γ = 1.046)
symmetries differ; this is because in the final MRCI+Q calculation, the two lowest
2

A’ states and the one lowest 2A” state are calculated, resulting in a larger A’reference space.

Permitting γ to differ provides additional flexibility for eachsymmetry species to make up for finite active space size.
By way of example, Fig 5.2 shows the results for correlation energy as
defined in Eq. 5.3 along the F + HCl reaction coordinate (described later) for each of
the AVnZ basis sets and for the CBS limit. The shaded portion reflects the “extra”
amount of energy recovered from the CBS values by correlation scaling. Such apattern of correlation energy as a function of basis set and reaction coordinate makes
it clear that correlation energy scaling is approximately equivalent to increasing thebasis size. More quantitatively, the inset in Figure 5.2 shows the ratio of correlationenergy for a given basis set to that inferred from a CBS extrapolation. The constancyof this ratio for each basis set provides further support that additional correlationenergy is recovered by an almost constant scale factor along the full reaction path.
Furthermore, this indicates that the shape of the PES is already quite accurately
determined at the MRCI+Q level, with correlation scaling resulting in only minor
adjustment of the reaction energetics to match experiment.

170

B. HF, HCl, and ClF asymptotic potentials
As a further test of the ab initio accuracy, the three asymptotic contributions
for the F(2P) + HCl Æ HF + Cl(2P) surface corresponding to HF, HCl, and ClF havebeen benchmarked with respect to both basis size and correlation energy scaling.
Table 5.2 summarizes the dissociation energies (kcal/mol), bond lengths (A),harmonic frequencies (cm-1), anharmonicities (cm-1), and rotational constants (cm-1)
for HF, HCl, and ClF calculated by DW-MCSCF/MCRI+Q for the various AVnZ (n= 2 - 4). The bond dissociation energies for each of the reactant and product channels
are both in very satisfactory agreement (∼0.1 - 0.2 kcal/mol) with experiment, andare significantly improved by correlation scaling. It is noteworthy that such scaling
has relatively little influence on the vibration/rotation diatomic constants, providing
additional support that correlation scaling results in only very minor changes in theoverall PES topography. Also remarkable is that the single scaling factor chosen to
match the experimental exothermicity for F + HCl reproduces properties of the ClF +H channel quite well, despite being quite remote to the reaction path for HF + Cl
formation.

C. Transition state
We begin with gradient saddle point searching for the F-H-Cl transition stategeometry and barrier height (See Figure 5.3 for transition states at multiple levels oftheory and Table 5.3 for a summary of transition state geometries). Previous ab initio
PES37 calculations have reported a strongly bent F-H-Cl angle (137o) and a classicalbarrier height of 4 - 6 kcal/mol.

At the current improved level of theory(MRCI+Q/CBS/Scaled), the transition state geometry is bent even further (123.5o),171with a classical barrier lowered to 3.8 kcal/mol. I note for later discussion that this isqualitatively different than the slightly bent transition state observed for F + H2,
which has important consequences for the resulting reaction dynamics. Furthermore,
the transition state is relatively loose, resulting in zero-point energy along the reaction
path that lowers the classical barrier height by ΔE ≈ 0.9 kcal/mol. This effectively
yields an adiabatic MRCI barrier of ∼2.9 kcal/mol.
By way of additional comparison, a single critical point UCCSD(T)/AVQZ
calculation64 also yields a highly bent transition state (118o), so this novel degree ofbending is clearly a robust feature of the true potential surface (see Figure 5.3).
Interestingly, such UCCSD(T) calculations suggest an even smaller classical barrierheight of 2.2 kcal/mol, though this is obtained using a spin unrestricted basis set.
There have been previous temperature dependent studies65 of F + HCl kinetics by
Houston and coworkers, which could in principle help resolve the differences in the
MRCI+Q/CBS/Scaled vs. UCCSD(T)/AVQZ barrier values. Unfortunately, the F +HCl experimental results exhibit strongly curved non-Arrhenius behavior, which
complicates comparison of empirical activation energies with actual ab initio
transition state barriers. However, if I examine the CI vector from the MRCI
calculations, the MCSCF ground state contribution accounts for only ∼60% of thetotal wave function probability at the transition state, with single excited state
amplitudes as large as 0.53. This indicates strong multireference character in thetransition state wavefunction and considerable mixing of higher electronic states notreflected in the UCCSD(T) calculation. Furthermore, inclusion of ∼0.9 kcal/mol
zero point effects would presumably similarly decrease the UCCSD(T)/AVQZ

172
adiabatic barrier to ∼ 1.3 kcal/mol. Such an adiabatic barrier height prediction would
already be lower than for F + HD (∼ 1.5 kcal/mol1), though experimental studies66
indicate F + HCl to have at least 3 - 4 fold lower cross sections for similar collisionenergies. Therefore, although further high-level surface calculations at the CCSD(T)level could be interesting to pursue, one can be cautiously confident from the present
MRCI/Q results that a strongly bent transition state (∼123o) with a vibrationallyadiabatic barrier height of ∼2.9 kcal/mol (before spin-orbit effects) represents a very
good approximation to the true values.

D. Reaction path
The reaction path provides a convenient set of benchmark points with which
to further verify that my choices of method, basis, correlation scaling, etc are valid
throughout configuration space and, in particular, away from the transition state.67 For
simplicity, the reaction path geometries are obtained from a SA-MCSCF ((3 state)calculation, for which analytical gradients are available for efficient reaction pathfollowing (MOLPRO does not yet provide analytic gradients for full MRCI+Q) in
both product and reagent directions along the path of steepest decent for the lowest
state. When the potential gradients in the asymptotic region become too small tofollow reliably, the reaction path is simply defined by increasing the F - HCl or HF Cl center of mass separations toward reactants or products.
Using the coordinates mapped out along the SA-MCSCF reaction path, I next
perform a 6 state dynamically weighted MCSCF calculation. In Figure 5.4, theenergies (a) and normalized dynamical weights (b) along the reaction path are shown
as a function of the reaction coordinate68. Notice that in each asymptotic region, the173
DW-MCSCF method effectively reduces to a ∼three-state calculation, i.e., sufficient
for accurately describing the triply degenerate p-hole in the lone F or Cl halogen
atoms. As the transition state region is entered, however, the two additional F- +HCl+(2Π) and HF+(2Π) + Cl- charge transfer states rapidly decrease in energy andbecome significantly weighted in the DW-MCSCF calculation and exhibit strong
effects on the reaction path curvature (see Figure 5.4a). As evident in Figures 5.4a, b,
dynamical weighting of the 6 lowest states more than ensures that all 5 of the relevant
open shell and charge transfer interactions are accounted for throughout the reaction
path, and in the transition state region in particular.
An internally contracted MRCI+Q calculation is then performed using thereference states from the DW-MCSCF calculations for each AVnZ n = [2,4] basis set.
The energies of both the DW-MCSCF and MRCI+Q calculations are extrapolated asa function of basis set to the CBS limit, with final refinement by scaling of thecorrelation energies by a constant γ to match experimental exothermicity. As noted
earlier, this correlation scaling method only makes minor adjustments in theenergetics of the system without influencing the shape of the PES. By way of
example, Figure 5.5 presents the energy profile along the reaction path for the seriesof MCSCF/CBS, MRCI+Q/CBS, and MRCI+Q/CBS/Scaled procedures. The large
shift between the MCSCF and MRCI+Q curves is due to the large amount ofcorrelation energy recovered by the MRCI+Q method, in contrast with the very smallchange (see magnification in Figure 5.5) between MRCI+Q and the MRCI+Q/Scaled
methods due to minor additional correlation energy that scaling recovers. Most
importantly, this dynamically weighted approach results in smoothly changing

174
energies along the reaction path, which represents a small but obviously key subset ofgeometries for dynamics on the full PES.

E. Entrance and exit channel van der Waals wells
As noted in the previous surface calculations37, there are 12A’ van der Waals
complexes in both entrance and exit channels. The entrance channel exhibits arelatively shallow F--HCl well depth (relative to the F + HCl asymptote) of ∼0.43
kcal/mol for a linear geometry with RHF = 4.50 a0 and RHCl = 2.44 a0. The effective
well depth in a full collision should decrease with rotational averaging; this is in good
agreement with crossed molecular beam experiments69, which place the rotationallyaveraged F + HCl van der Waals well depth at 0.32 kcal/mol. As expected from thelarger dipole-quadrupole and dipole-induced dipole interactions70, the exit channel
exhibits a stronger van der Waals well of 1.52 kcal/mol, for a linear geometry with
RHF = 1.75 a0 and RHCl = 4.73 a0. Interestingly, the CI vectors for both entrance andexit channel wells show significant contributions (∼10% ) from excited state
amplitudes, further supporting the importance of high-level multireference ab initio
methods in constructing a full potential energy surface.

IV. Potential Energy Surfaces
In order to construct the full 3D potential energy surface, calculations havebeen performed along a grid of points with the ranges RHF = [1.3 - 15.0] a0 and RHCl =[1.8 - 15.0] a0 and θF-H-Cl = [180, 150, 120, 90, 60, 30, 0], sampling only energies
lower than 50 kcal/mol above the entrance channel energy. These 3230 points arechosen on a grid such that the density is maximized in important regions of the175
surface (transition state, van der Walls wells) while still sampling a large enoughregion to get an accurate description of the PES. Also included in the full set are 121
points sampling directly along the reaction path. Following the methods outlined by
Aguado and Paniagua,71 these ab initio points are fit to an analytical multibody
expansion:
V (RHF , R HCl , RClF ) = V (1) +∑ V (R ) + V ( R , R( 2)i∈{ HF , HCl ,ClF }

( 3)

i

HF

HCl

, RClF ) . (5.6)

where V(1) is the energy of the separated atoms which was subtracted from the abinitio energies, V(2) are the diatomic potentials, and V(3) reflects the remaining 3-bodyinteractions.

A. Two body terms
The two body term is based on a sum of modified Rydberg functions71:
V ( 2 ) (R ) =

(

)

i
(2)
C 0 − β1( 2 ) R Ne+ ∑ Ci R e − β2 R .
R
i =1

(5.7)

where R represents the pairwise interatomic distances between the 3 atoms. Ab initio
points for each of the diatomics are fit to a 9th order polynomial (N = 9). For the Cl +HF reagent channel, 19 ab initio points between RHF = 1.2 a0 and RHF = 3.0 a0 havebeen used, yielding fits with a global RMS deviation of 0.039 kcal/mol. Similarly for
HCl, 26 points between RHCl = 1.6 a0 and RHCl = 4.0 a0 have been fit to a global RMS
of 0.015 kcal/mol. Although the ClF region of the surface is not sampled directly inthe reactive scattering experiments, the ClF potential is sampled for 21 points
between RClF = 2.5 a0 and RClF = 4.5 a0, yielding a global diatomic fit with a verysatisfactory RMS deviation of ∼0.046 kcal/mol.

176

B. Three body terms
The three body term simply reflects the remaining difference between ab
initio points after subtracting the diatom V(2) and atom V(1) energies. Since this term
must vanish in any of the diatomic asymptotes, a commonly used functional form is
simply a linear combination of two-body product functions, i.e.,M(

V (3) (R1 , R2 , R3 ) = ∑ Dijk R1e − β1 R1( 3)) (R ei2− β 2( 3 ) R2

) (R ej3( 3)− β 31
R3

),
k

(5.8)

i , j ,kwhere the total order of the polynomial is restricted to i + j + k ≤ M. The explicit
numerical least squares weighting strategy for obtaining this fit is modeled after
efforts to develop the Stark-Werner F + H2 surface1. Specifically, the ab initio points
along the reaction path are given a weight of 100. Of the 3230 grid points, the 2200
geometries with energies less than 15 kcal/mol above the entrance channel areweighted by 10, with the remaining 1030 points with energies between 15 and 50
kcal/mol above the entrance channel receiving a weight of 1. The resulting global
RMS deviation for least squares fit with a 9th order polynomial (M = 9) is ∼0.72
kcal/mol, i.e., 3.5 times better than the previous surface37, while the RMS deviation
for points along the reaction path improves to ∼0.29 kcal/mol. A 2D cut through the
PES at θ = 123.5o (i.e., the MRCI+Q/CBS/ Scaled transition state bend angle ) isshown in Figure 5.6. As I am restricting my initial focus to the ground state, this does
not yet reflect spin-orbit coupling effects, though they will obviously be incorporated
in a full non-adiabatic treatment of the lowest three surfaces54.
Although this fit is already quite good, it is important to note that RMS
deviation for F + HCl surface may tend to overestimate how well the surface is in fact

177
represented. Specifically, the surface has two conical intersection seams in theentrance (F(2P) + HCl) and exit (HF + Cl(2P)) channels, indeed, even at energies
comparable to the transition state.

Although a proximity metric in multiple

dimensions is not uniquely defined, the points on the conical intersection seam thatare “closest” to the reaction path (as simply evaluated by sum of Euclidian
displacements for each atom) are RHF = 3.578 a0, RHCl = 2.459 a0, E = 3.55 kcal/mol in
the entrance channel and RHF = 2.166 a0, RHCl = 2.770 a0, E = -3.23 kcal/mol in the exit
channel.

These conical intersection seams introduce regions of derivative

discontinuities in the adiabatic PES, which is precisely where the least squares
multibody expansion deviates most from the ab initio values. Indeed, in somerespects, these are the most dynamically interesting regions of the multiple surfaces
and are clearly responsible for the role of non-adiabaticity in the reaction dynamics. A
more complete treatment of such non-adiabatic effects therefore requires
simultaneous diabitization and fitting of all three lowest surfaces, which is outside thescope of this paper and will be presented elsewhere.54

V. Discussion
As has been noted in other PESs for reaction systems with open highly
electronegative open shell atoms40, the electronic wavefunction can be highly
multireferential, i.e., significant amplitudes of excited reference state configurations
become mixed with the ground state.

The absence of such mixing (i.e., single

reference behavior) can be quantified by χ, the amplitude of the ground reference
configuration in the final MRCI+Q wavefunction. To further demonstrate the crucial
importance of multireference methods in the F + HCl PES, Figure 5.7 plots the

178
magnitude of the ground state reference configuration (i.e.,χ2) as a function of thereaction coordinate. The strong and highly structured decrease in χ2 clearly indicatesthat excited reference configurations constitute a significant portion of the wave
function in the broad vicinity around the transition state. Interestingly, there are two
structured dips in the χ2 plot where the excited state admixture is greatest; these
correspond quite nicely to points along the reaction coordinate closest to the conicalintersection seams (marked with a *).

Clearly, significant mixing of excitedconfigurations is taking place over the entire barrier region of the surface, further
underscoring that single reference ab initio methods cannot be expected to adequately
describe the electronic wavefunction.
The importance of these conical intersection seams in F + HCl dynamics will
depend on how intimately the full 3D wavefunction samples these regions on the way
toward reaction. The extent of this sampling is not immediately obvious, since theconical intersections occur only for a collinear configuration (θ = 0o, 180o), and yet
the transition state is quite strongly bent (θ = 123.5o). To provide more insight into
this question, 2D contour plots as a function of RHF and RHCl are shown in Fig, 5.8,
for both collinear and transition state bending angles. The dashed lines in Fig 5.8b
represent the coordinates of the conical intersection seams, with the asterisks
indicating the positions closest to the (bent) reaction path. These points clearly occur
high up on the reaction barrier and for RHF and RHCl bond lengths quite comparable tothose of the transition state.
Based on the above discussion, one relevant question is therefore the“stiffness” of the bending potential, evaluated near the transition state bond lengths.

179
Figure 5.9 displays 1D cuts through the full potential as a function of the bend angle
θ. The dashed line cut is for RHF and RHCl coordinates fixed at the optimized
transition state values. The strongly bent nature of the transition state is clearly
evident with a modest ∼3 kcal/mol higher barrier at the linear configuration.
However, this substantially underestimates the energetic proximity to the conicalintersection seam; specifically, if one also adjusts RHCl and RHF to the conicalintersection values represented by the two asterisks, the energies drop to 3.55
kcal/mol and -3.23 (referenced to the F + HCl asymptotic entrance channel). Thus
even the higher of these two points on the seam is lower than the transition state by∼0.25 kcal/mol. Simply stated, all reagents with sufficient energy to surmount theclassical reaction barrier have energetic access to non-adiabatic seam regions for both
reagent and product channel sides of the transition state.
A more relevant angular slice through the potential is therefore at the RHCl
and RHF geometries associated with the two conical intersection seam values closest
to the transition state. Such a slice is also displayed in Fig 5.9 for the (higher energy)
reagent channel seam region, which now indicates a much softer bending potentialand even lower barrier (∼ 1 kcal/mol) to the linear configuration than at the transition
state bond lengths. Note that the barrier at linearity on this second curve now
explicitly corresponds to a point directly on the conical intersection seam.
Furthermore, this curve intersects with the reaction path (point marked “P”) at anenergy only 0.9 kcal/mol below this linear conical intersection, underscoring theenergetic proximity of reactants to these regions of high non-adiabaticity. A morequantitative metric for how accessible the conical intersection region is to reagents

180
along the reaction path can be obtained by solving for the ground state 1D angularwave function at fixed RHCl and RHF. Since the F-H-Cl moment of inertia becomessingular at the linear configuration, this requires taking advantage of the rigid bender
methods developed by Hougen, Bunker and Johns. Using these methods, the lowest
several eigenenergies have been calculated and are displayed for both the transition
state and conical intersection curves in Figure 5.9. Note that linear conicalintersection region becomes classically accessible even for a few quanta of bend
excitation. A more detailed investigation of these conical intersections with full spinorbit effects and fitting of the resulting diabatic 22A’ and 12A” surfaces is explored in
the next chapter.54
The most striking aspect of the F(2P) + HCl → HF + Cl(2P) PES is the
sharply bent MRCI+Q/CBS/Scaled transition state geometry of 123.5o. This feature
is qualitatively reiterated at any level calculation, from previous single reference
surface calculations (θ ≈137o) to single point UCCSD(T)/AVQZ (118o) calculations.
Fig 5.8 exhibits the typical “early” barrier expected for a strongly exothermic
reaction.

However, both the barrier height and “early” vs. “late” nature of thetransition state changes quite dramatically with the F-HCl bend angle θ. As noted
previously by Skodje and coworkers,64 this strongly bent transition state can have asignificant influence on the reaction dynamics, specifically in reducing the reaction
cross section for JHCl » 0 reagents.
A simple picture for such behavior can be obtained as follows. 2D transition
state geometries and barrier heights have been obtained from the full 3D surface at a
fixed angle, varying only the two bond length coordinates RHCl and RHF. These

181
effective transition state barrier heights and HCl bond lengths are plotted as afunction of θ in Figure 5.10a. For more collinear geometries (θ ≈ 180), the 2D barrierheight and transition state HCl bond length prove to be relatively insensitive to θ.
Furthermore, rHCl remains quite close to the free equilibrium value, as characteristicof “early barrier” behavior even out to near the transition state bend angle (θ ≈ 123.5).
For more strongly bent geometries (θ < 120), however, the 2D transition state HClbond length grows dramatically with decreasing θ, implying a rapid shift from “early”
to “late” barrier dynamics, accompanied by a similarly precipitous increase in the 2D
barrier height. By way of example, the shaded portion in Fig 5.10a displays therelatively small range of F-H-Cl angles that correspond to 2D barriers below a typical
5 kcal/mol experimental collision energy. Alternatively stated, the strongly benttransition state nature of the F + HCl PES translates into a strongly reduced angular
acceptance window for reaction, resulting in anomalously small cross sections evenfor an “early” barrier process at energies well above the transition state threshold.
Indeed, the strongly bent nature of the F-H-Cl transition state generates
additional dynamical restrictions on the reaction probabilities. As a simple physicalpicture for this, consider a collision of F with HCl oriented initially in a bent
configuration. Particularly for such a “heavy + light-heavy” mass combination, thispresents a large moment arm for prerotating the HCl reagent in the entrance channel,
thereby decreasing θ below the optimum transition state angle. From Fig 5.10a, this“entrance channel rotation” translates into higher effective barrier heights, thereby
shutting off the probability of a successful reaction. To quantify how large these“entrance channel rotation” effects might be during a F + HCl collision, preliminary

182
classical trajectory simulations have been performed on the current PES. Initialconditions are chosen such that the F atom initial velocity is directly toward the Hatom at 5.0 kcal/mol center of mass collision energy (i.e., 1.2 kcal/mol above
threshold) for an initial F-HCl bend angle between 60o and 180o. The resulting F-HCl bending angle is then recorded at the location of closest approach of the F atom to
the HCl center of mass. Results of this sample trajectory analysis are shown in Figure
10b, where the change in angle Δθ at the closest approach distance has been plottedagainst the initial F-H-Cl angle. Noteworthy is that for all incident angles, Δθ is quite
large and negative (i.e., the F-H-Cl angle becomes more bent). In conjunction with
Figure 5.10a, this implies a large increase in the effective barrier height, which begins
to explain the exceptionally low probability of reactive trajectories observed even
above threshold.
Indeed, based on this elementary classical analysis, a successful reaction
would require collision at very nearly collinear geometries to “prerotate” the HCl into
the appropriate transition state bend angle. This range of angles is indicated by the
steeply rising edge feature in Fig 5.10b, which translates into a considerable
narrowing of reactive acceptance window at 5.0 kcal/mol, specifically from θ ≈ 33o toθ ≈ 2o. Furthermore, when one takes the correct sin(θ) area element into account, thismodel predicts a reactive event with only 0.21% of the hard sphere cross section. The
dashed line in Figure 5.10b shows similar predicted deflection plots for trajectories
with the H atom replaced with D. As expected, the deflection is somewhat decreased
due to the more massive D atom but still quite large and sampling a region of theconfiguration space where the barrier is too large to react. Note that such a dynamical

183
picture due to pre-rotation into the transition state is most appropriate for non-rotating
HCl reagent molecules. This is a likely explanation for observations obtained fromclassical and full 3-D quantum dynamics on this surface64 , which show very small
cross sections for reactions starting in JHCl = 0 but which increase dramatically with
JHCl>0. Whether such a simple picture for the role of transition state bending angle
effects will be dynamically relevant in other reaction systems represents aninteresting question for further exploration.

VI. Summary and Conclusion
This chapter presents results from high-level ab initio calculations on the 12A’electronic state for the F(2P) + HCl → HF + Cl(2P) reaction, along with an analytical
fit to this surface in full 3D. The ab initio data for the analytical fit are calculated at a
grid of 3D points using

dynamically weighted MCSCF methods with 6 states,

followed by high-level multireference (MRCI+Q) methods for a series of correlationconsistent basis sets (AVnZ, n = 2,3,4). These results are then extrapolated to thecomplete basis set limit (CBS), with a single global correlation scaling parameter tomatch theoretical and experimental exothermicities. Sequential fits of this surface toatom, pairwise, and full three body terms have been obtained and reported for the
ground adiabatic state. This represents a necessary first step toward obtaining all
multistate diabatic electronic surfaces for use in quantum scattering calculations andrigorous comparison with state-resolved experimental data.
The successful extraction of this full 3D surface demonstrates the robust
applicability of DW-MCSCF method to obtain multiple smooth potentials in systems
with strongly avoided crossings and charge transfer character at the transition state.

184
Furthermore, the nominally “ground” electronic state wave function proves to behighly multireference in character, particularly in the chemically important region
near the transition state. This underscores the importance of multireference methods
in open shell reaction dynamics, particularly for highly electronegative atoms withstrong charge transfer contributions and permits extraction of a much improved
surface for comparison with previous single reference state methods. A simple oneparameter correlation scaling procedure is introduced that recovers relatively minor
correlation energy effects due to frozen core orbitals, benchmarked to match theexperimentally measured exothermicity.
This surface exhibits several intriguing dynamical features for further
exploration. First of all, elementary chemical arguments are presented for open shell
halogen/H atom abstraction reactions that predict the appearance of conicalintersections close to the transition state region. The presence of two such seams ofconical intersections is confirmed in the ab initio calculations, occurring at RHCl and
RHF bond lengths close to and energies slightly below the transition state values. Most
importantly, this may suggest a major role for multiple electronic surfaces and nonadiabatic effects for F + HCl even under thermal reaction conditions. In fact, thetemperature dependent studies by Houston and coworkers indicate a curiously nonArrhenius kinetic behavior. One plausible speculation is that this could arise from
non-adiabatic hopping effects due to temperature dependent access to the entrance
and exit conical intersection seams on the surface. This would be extremelyinteresting to explore, both by full multisurface dynamics calculations as well asrovibrational and spin orbit state-resolved scattering experiments.

185
Second, the transition state is quite strongly bent, with a F-H-Cl bend angle ofθ ≈ 123.5o. This is in contrast with the more nearly linear transition states observedfor F + H2 surfaces and may have a significant impact on the reaction dynamics inmultiple ways. For example, 2D slices (RHF, RHCl) of the potential surface for a fixed
bend angle indicate a dramatic increase in the effective transition state barrier and
RHCl bond length with decreasing θ below 120. This is consistent with rapid shift
from early to late barrier dynamics as the F-H-Cl angle becomes more bent, which,particularly for a heavy + light-heavy mass combination, is likely to occur naturally inthe entrance channel. This is further verified by trajectory calculations, which reveal astrong “entrance channel rotation” of the HCl reagent, greatly narrowing the angular
window of acceptance for scaling the transition state and achieving successful
reaction. Such strongly bent transition states clearly offer a novel platform and testing
ground for predictions of reaction dynamics as a function of rotational and vibrational
quantum state, one of many directions which may prove interesting to pursue with theavailability of the present surface.

186

TablesTable 5.1: ΔErxn in kcal/mol of F(2P) + HCl Æ HF + Cl(2P) calculated by various ab
initio methods. Calculated energies are corrected for zero-point and spin-orbit which
are not included in the single surface ab initio calculations. The experimental ΔE is-33.060 +/- 0.001 kcal/mol.

Basis
AVDZ

AVTZ

AVQZ

CBS

Error
(CBS)

-19.3
-23.1
-30.41-31.33-18.83-22.87-29.98-31.04-18.67-22.75-29.91-30.94-18.57-22.68-29.87-30.88-14.49-10.38-3.19-2.18

Method
HF
MCSCF
MRCI
MRCI+Q

187
Table 5.2: Diatomic constants from ab initio PES calculated at MRCI+Q/AVnZ compared with the correlation scaled
values and experiment.

HF

HCl

ClF

Basis

De
(kcal/mol)Dokcal/mol)

ZPE(kcal/mol)

AVDZ
AVTZ
AVQZ
CBS
Scaled
expt.
AVDZ
AVTZ
AVQZ
CBS
scaled
expt.
AVDZ
AVTZ
AVQZ
CBS
scaledexpt.

133.84138.08139.61140.50141.58141.63100.29104.18105.75106.67107.35107.36
51.44
56.86
59.16
60.54
61.73
61.50

128.10132.24133.74134.61135.69135.6596.08
99.95
101.52102.44103.11103.05
50.34
55.74
58.02
59.39
60.57
60.30

5.74
5.84
5.87
5.89
5.89
5.98
4.21
4.23
4.23
4.23
4.24
4.31
1.09
1.12
1.14
1.15
1.16
1.13

re(A)

ωe
(cm-1)

ωexe(cm-1)

0.92600.92170.91850.91650.91640.91681.29381.27901.27751.27691.27661.27461.68471.6455
1.6362
1.63151.6304
1.62834038.14118.64134.64143.24143.04138.32966.92975.52976.22978.22980.22990.9
770.2
787.7
803.2
811.7
816.3
786.2

76.8
77.1
77.1
77.1
76.4
89.9
49.0
47.5
46.1
45.6
45.4
52.8
6.5
6.0
5.9
5.8
5.8
6.2

Be
(cm-1)

αe(cm-1)

20.6086
20.8440
20.9766
21.059921.0586
20.9557
10.3740
10.6048
10.6291
10.6388
10.6433
10.5934
0.48280.50600.51180.51480.51550.5165

0.707
0.695
0.700
0.705
0.701
0.798
0.299
0.301
0.297
0.295
0.294
0.307
0.004
0.005
0.004
0.004
0.004
0.004

188
Table 5.3: The geometries and energies of the transition state for a range of methods
and basis sets. In particular, the transition state geometry is strongly bent at all levelsof theory.

Method
37

UMP2
37
PUMP2
37
PUMP4
UCCSD(T)64
MRCI+QBasis
6-311G(3d2f,3p2d)
6-311G(3d2f,3p2d)
6-311G(3d2f,3p2d)
AVQZ
AVDZ
AVTZ
AVQZ
CBS
ScaledBend Angle
ΔE‡(degree) (kcal/mol)
137.4
137.4
137.4
118.0
126.2
126.4
125.9
125.7
123.5

6.2
4.7
4.0
2.2
4.2
4.2
4.2
4.2
3.8

189

FiguresFigure 5.1: AVnZ (n = 2-5) extrapolation for F(2P) + HCl and HF + Cl(2P) energies
to the CBS limit, indicating clear convergence (see inset) in the exothermicity (ΔE) as
a function of n. The residual error in ΔE (dashed line) is due to active space size in
the MRCI+Q calculations, as included in correlation scaling

190

Figure 5.2: Correlation energy (EMRCI-EMCSCF) recovered for each basis set, AVnZ. The nearly constant ratio of recovered correlation energy relative to the CBS
extrapolation (see inset) illustrates how correlation scaling effectively mimics a larger
active space.

191

Figure 5.3:

F-H-Cl transition state geometry at various levels of theory, allconfirmingasharplybenttransition

MRCI+Q/CBS/Scaled; UCCSD(T)/AVQZ.

state:(a)

MRCI+Q/AVTZ;

(b)

192

Figure 5.4: The DW-MCSCF energies (a) and weights (b) as a function of reaction
coordinate, with up to 6 states dynamically weighted. Note the strongly avoided
crossings of F- + HCl(2Π)+ and HF(2Π)+ + Cl- charge transfer (CT) states withsurfaces correlating to the ground state (2P) asymptotes. Note also the smooth weightadjustment from a 3-state calculation at either asymptote to > 5 states near thetransition state geometry.

193

Figure 5.5: 12A’ reaction path at the MCSCF, MRCI+Q, and correlation scaled levels,
indicating relatively minor effects due to correlation scaling. Such scaling makes up
for 2.18 kcal/mol in ΔE and lowers the barrier by 0.43 kcal mol.

194

Figure 5.6: 2D slice (θ = 123.5o) of the fitted 12A’ surface with global RMS of 0.79
kcal/mol.

195

Figure 5.7:

Ground reference state coefficient in the MRCI+Q wavefunction,

indicating strong admixture of excited state character near the transition stategeometries closest to the conical intersection (marked by asterisks).

196

Figure 5.8: 12Α’ PES RHF vs. RHCl contours for (a) transition state (θ = 123.5o) and(b) collinear (θ = 180o) bend angles. Conical intersection seams occur to each side of
the transition state; the seam location closest to the reaction path is marked by theasterisk. Contour spacing is 5 kcal/mol with respect to F + HCl(Re).

197

Figure 5.9: F-H-Cl bending potential corresponding to transition state (dashed) andconical intersection seam (solid) RHF and RHCl bond lengths, with eigenvalues
obtained from a rigid bender analysis. Note that the conical intersection is lower than
the transition state with only a ∼1 kcal/mol barrier to linearity from the reaction path.

198

Figure5.10: a) 2D effective transition state barrier height (solid) and HCl bond length(dashed) as a function of F-H-Cl bending angle (dotted line is Re for HCl). Asθ decreases below ∼120o, both the HCl bond length and barrier height increase
dramatically from transition state values, characteristic of a shift from an “early” to“late” transition state. b) Trajectory analysis for F-HCl angular deflection (Δθ) as afunction of initial bend angle, θ. For J = 0, only a very narrow range (dark grey) of
incident scattering angles are successfully “prerotated” (light grey) into theappropriately bent transition state angle to react.

199

References

1

K. Stark and H. J. Werner, J. Chem. Phys. 104, 6515 (1996).

2

M. Hayes, M. Gustafsson, A. M. Mebel, and R. T. Skodje, Chem. Phys. 308,
259 (2005).

3

Y. R. Tzeng and M. H. Alexander, J. Chem. Phys. 121, 5812 (2004).

4

Y. Zhang, T. X. Xie, K. L. Han, and J. Z. H. Zhang, J. Chem. Phys. 119,12921 (2003).

5

S. D. Chao and R. T. Skodje, J. Chem. Phys. 113, 3487 (2000).

6

M. H. Alexander, H. J. Werner, and D. E. Manolopoulos, J. Chem. Phys. 109,
5710 (1998).

7

F. J. Aoiz, L. Banares, B. MartinezHaya, J. F. Castillo, D. E. Manolopoulos,
K. Stark, and H. J. Werner, J. Phys. Chem. A 101, 6403 (1997).

8

J. F. Castillo, D. E. Manolopoulos, K. Stark, and H. J. Werner, J. Chem. Phys.

104, 6531 (1996).
9

S. A. Nizkorodov, W. W. Harper, W. B. Chapman, B. W. Blackmon, and D. J.
Nesbitt, J. Chem. Phys. 111, 8404 (1999).

10

S. A. Nizkorodov, W. W. Harper, and D. J. Nesbitt, Faraday Disc. 113, 107(1999).

11

W. B. Chapman, B. W. Blackmon, S. Nizkorodov, and D. J. Nesbitt, J. Chem.
Phys. 109, 9306 (1998).

200
12

W. B. Chapman, B. W. Blackmon, and D. J. Nesbitt, J. Chem. Phys. 107,
8193 (1997).

13

M. Baer, M. Faubel, B. Martinez-Haya, L. Rusin, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 110, 10231 (1999).

14

M. Baer, M. Faubel, B. Martinez-Haya, L. Y. Rusin, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 108, 9694 (1998).

15

L. Y. Rusin and J. P. Toennies, Phys. Chem. Chem. Phys. 2, 501 (2000).

16

G. Dharmasena, K. Copeland, J. H. Young, R. A. Lasell, T. R. Phillips, G. A.
Parker, and M. Keil, J. Phys. Chem. A 101, 6429 (1997).

17

G. Dharmasena, T. R. Phillips, K. N. Shokhirev, G. A. Parker, and M. Keil, J.
Chem. Phys. 106, 9950 (1997).

18

M. Faubel, L. Rusin, S. Schlemmer, F. Sondermann, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 101, 2106 (1994).

19

D. M. Neumark, A. M. Wodtke, G. N. Robinson, C. C. Hayden, and Y. T.
Lee, J. Chem. Phys. 82, 3045 (1985).

20

D. M. Neumark, A. M. Wodtke, G. N. Robinson, C. C. Hayden, K. Shobatake,
R. K. Sparks, T. P. Schafer, and Y. T. Lee, J. Chem. Phys. 82, 3067 (1985).

21

H. Kornweitz and A. Persky, J. Phys. Chem. A 108, 8599 (2004).

22

Y. R. Tzeng and M. H. Alexander, J. Chem. Phys. 121, 5183 (2004).

23

Y. R. Tzeng and M. Alexander, Phys. Chem. Chem. Phys. 6, 5018 (2004).

24

Y. Zhang, T. X. Xie, K. L. Han, and J. Z. H. Zhang, J. Chem. Phys. 120, 6000
(2004).

201
25

T. X. Xie, Y. Zhang, M. Y. Zhao, and K. L. Han, Phys. Chem. Chem. Phys. 5,2034 (2003).

26

S. C. Althorpe, J. Phys. Chem. A 107, 7152 (2003).

27

S. C. Althorpe, Chem. Phys. Lett. 370, 443 (2003).

28

S. H. Lee, F. Dong, and K. P. Liu, J. Chem. Phys. 116, 7839 (2002).

29

K. P. Liu, Nuc. Phys. A 684, 247C (2001).

30

R. T. Skodje, D. Skouteris, D. E. Manolopoulos, S. H. Lee, F. Dong, and K. P.
Liu, Phys. Rev. Lett. 85, 1206 (2000).

31

R. T. Skodje, D. Skouteris, D. E. Manolopoulos, S. H. Lee, F. Dong, and K.
Liu, J. Chem. Phys. 112, 4536 (2000).

32

W. W. Harper, S. A. Nizkorodov, and D. J. Nesbitt, J. Chem. Phys. 116, 5622
(2002).

33

O. D. Krogh and G. C. Pimentel, J. Chem. Phys. 67, 2993 (1977).

34

A. M. G. Ding, L. J. Kirsch, D. S. Perry, J. C. Polanyi, and J. L. Schreiber,
Faraday Disc. Chem. Soc. 55, 252 (1973).

35

R. Sayos, J. Hernando, R. Francia, and M. Gonzalez, Phys. Chem. Chem.
Phys. 2, 523 (2000).

36

B. Y. Tang, B. H. Yang, K. L. Han, R. Q. Zhang, and J. Z. H. Zhang, J. Chem.
Phys. 113, 10105 (2000).

37

R. Sayos, J. Hernando, J. Hijazo, and M. Gonzalez, Phys. Chem. Chem. Phys.

1, 947 (1999).
38

H. Kornweitz and A. Persky, J. Phys. Chem. A 108, 140 (2004).

202
39

D. J. Nesbitt, M. P. Deskevich, M. Wocjik, M. Ziemkiewicz, A. Zolot, and E.
Whitney, Abstracts Of Papers Of The American Chemical Society 227, U334
(2004).

40

M. P. Deskevich, D. J. Nesbitt, and H. J. Werner, J. Chem. Phys. 120, 7281
(2004).

41

H. J. Werner and P. J. Knowles, J. Chem. Phys. 89, 5803 (1988).

42

P. J. Knowles and H. J. Werner, Chem. Phys. Lett. 145, 514 (1988).

43

P. J. Knowles and H. J. Werner, Theoretica Chimica Acta 84, 95 (1992).

44

K. A. Peterson, D. E. Woon, and T. H. Dunning, J. Chem. Phys. 100, 7410
(1994).

45

F. B. Brown and D. G. Truhlar, Chem. Phys. Lett. 117, 307 (1985).

46

B. Ramachandran and K. A. Peterson, J. Chem. Phys. 119, 9590 (2003).

47

M. Meuwly and J. M. Hutson, J. Chem. Phys. 119, 8873 (2003).

48

M. L. Dubernet and J. M. Hutson, J. Chem. Phys. 101, 1939 (1994).

49

M. L. Dubernet and J. M. Hutson, J. Phys. Chem. 98, 5844 (1994).

50

H.-J. Werner, P. J. Knowles, R. Lindh, M. Schutz, P. Celani, T. Korona, F. R.
Manby, G. Rauhut, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper,
M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, A. W.
Lloyd, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R.
Pitzer, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni, and T. Thorsteinsson,
Molpro,version2002.6,apackageofabinitioprograms

http://www.Molpro.Net) (Birmingham, UK, 2003).
51

S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem. 8, 61 (1974).

(see

203
52

M. R. A. Blomberg and P. E. M. Siegbahn, J. Chem. Phys. 78, 5682 (1983).

53

J. Simons, J. Phys. Chem. 93, 626 (1989).

54

M. P. Deskevich and D. J. Nesbitt (work in progress).

55

D. E. Woon and T. H. Dunning, J. Chem. Phys. 98, 1358 (1993).

56

T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).

57

R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796
(1992).

58

T. H. Dunning, K. A. Peterson, and A. K. Wilson, J. Chem. Phys. 114, 9244
(2001).

59

L. M. Cheung, K. R. Sundberg, and K. Ruedenberg, Int. J. Quantum Chem.

16, 1103 (1979).
60

H. J. Werner and P. J. Knowles, J. Chem. Phys. 82, 5053 (1985).

61

P. J. Knowles and H. J. Werner, Chem. Phys. Lett. 115, 259 (1985).

62

J. D. D. Martin and J. W. Hepburn, J. Chem. Phys. 109, 8139 (1998).

63

W. T. Zemke, W. C. Stwalley, J. A. Coxon, and P. G. Hajigeorgiou, Chem.
Phys. Lett. 177, 412 (1991).

64

M. Y. Hayes, M. P. Deskevich, D. J. Nesbitt, K. Takahashi, and R. T. Skodje,
J. Chem. Phys.(in press) (2005).

65

E. Wurzberg and P. L. Houston, J. Chem. Phys. 72, 5915 (1980).

66

A. Zolot and D. J. Nesbitt (work in progress).

67

F. Eckert and H. J. Werner, Theor. Chem. Acc. 100, 21 (1998).

68

K. Fukui, Acc. Chem. Res. 14, 363 (1981).

204
69

V. Aquilanti, D. Cappelletti, F. Pirani, L. Y. Rusin, M. B. Sevryuk, and J. P.
Toennies, J. Phys. Chem. 95, 8248 (1991).

70

J. M. Hutson, J. Chem. Phys. 89, 4550 (1988).

71

A. Aguado, C. Tablero, and M. Paniagua, Comput. Phys. Comm. 108, 259(1998).

205

Chapter 6:

Multireference configuration interaction

calculations for F(2P) + HCl → HF + Cl(2P) reaction:
Correlation-scaled diabatic potential energy surfaces.
I. Introduction
Understanding reaction dynamics at the quantum state resolved level has beena major challenge to the theoretical and experimental chemical physics community
for many years. As a result of such interest, there has been much development
towards both accurate ab initio potentials and numerically exact quantum dynamical
calculations on these surfaces. One particularly well-studied system is the F(2P) + H2
Æ HF + H(2S) reaction, for which there now exists an exceptionally accurate

potential surface1,2 from the Werner group, as well as numerous theoretical3-8
calculations and experimental9-20 papers on state-to-state reactive scattering. Also ofinterest has been the similar but notably different “heavy + light-heavy” isotopic
variation on this reaction system, F(2P) + HD Æ HF + D(2S), for which long soughtafter quantum transition state resonances have been both predicted 21-31 andexperimentally32 observed. In particular, stabilization calculations at the resonanceenergy in the transition state region yield a resonance wave function that can be
described classically as an excited H atom motion in a quasi bound HF(v = 3)--D
state. As a simple physical picture, the light H atom executes rapid oscillations
between the relatively heavy F and D atoms, eventually predissociating into statescorresponding to HF(v = 2,J) + D. Ongoing theoretical and experimental work on this206
H atom transfer triatomic system has revealed the considerable dynamical
sophistication still to be found even in “simple” heavy + light-heavy reactions, which
continue to offer fundamental challenges for first principles theoretical calculations.
The clear observation of transition state resonance dynamics in F + HD
suggests that there may be similar resonances in other heavy + light-heavy reaction
systems such as F(2P) + HX. One such reaction of particular interest is the F(2P) +HCl Æ HF + Cl(2P) system, the high-level benchmark calculations of 3-state
adiabatic and diabatic potential energy surfaces which represent the primary focus of
this chapter. The reason for this choice is several-fold. First of all, F(2P) + HCl Æ HF
+ Cl(2P) has a very nearly identical exothermicity to the F(2P) + HD Æ HF + D(2S)
reaction (ΔEHF + Cl = -33.1 kcal/mol vs ΔEHF + H = -31.2 kcal.mol), as well as anexperimentally accessible barrier height of Ebarr = 3-5 kcal/mol. Secondly, the total
electron count (Nelec = 27) and number of nuclear degrees of freedom are challenging
but still manageable with modern high-level multireference configuration interaction(MRCI) ab initio efforts for mapping ground and excited potential surfaces out in full
dimensionality. Thirdly, simple electrostatic arguments predict the presence of nonadiabatic surface crossings and 1D seams of conical intersections that should beenergetically accessible at typical reactive collision energies. Finally, the open shell
halogen in the F(2P) + HCl Æ HF + Cl(2P) reaction has multiple spin-orbit
configurations of products and reactants in the entrance and exit channels, which
provides a potential experimental handle to distinguish surface hopping events as well
as novel opportunities to study non-adiabatic dynamics on multiple electronic
surfaces.

207
Previous theoretical studies on this surface have included quasiclassicaltrajectory (QCT) calculations performed on an empirical LEPS surface33, QCT34andfull dimensional time-dependant wavepacket 35 calculations obtained on a newer ab
initio (PUMP2/6-311G(3d2f,3p2d)) surface36, a LEPS surface created from recent
experimental kinetics information37, and a correlation-scaled MRCI surface from ourgroup38. All previous calculations33-35 have focused specifically on the ground (12A’)electronic state, though study of other F atom abstraction systems has revealed
evidence for avoided crossings, charge transfer states, seams of conical intersection,etc. For example, it has been shown that there is significant non-adiabatic behavior
for F(2P) + H2->HF + H(2S) in the entrance channel3,23, and there is also direct
experimental evidence in 4 atom systems (i.e. F (2P) + H2O Æ HF + OH(2Π) thatexhibit strong non-adiabatic coupling behavior in both the entrance and exit
channels39. It is likely that similar effects are present in F(2P) + HCl Æ HF + Cl(2P),
the dynamics of which clearly necessitate a high-level multireference electronic
treatment. In this work, I present a new 3-state (12Σ, 12Πx, 12Πy) diabatic PES for theF(2P) + HCl reaction, based on high-level multireference calculations including spinorbit contributions.
The F(2P) + HCl Æ HF + Cl(2P) reaction presents many of the samechallenges seen earlier in the F(2P) + H2O Æ HF + OH(2Π) reaction40, namely, thepresence of highly electronegative atoms that significantly lower the charge transfer(CT) state (i.e. F- + HCl + (2Π) and HF + (2Π) + Cl-) energies in the barrier region. In
addition, there is significant non-adiabatic coupling and avoided crossing behavior
between the 12A’ and 22A’ surfaces, which requires explicit multireference208
calculations to describe accurately.

In light of these considerations, ab initio

calculations for this PES are obtained via internally contracted multireference
configuration interaction (MRCI)41-43 methods using reference states taken from a
dynamically weighed multiconfigurational self-consistent field (DW-MCSCF)40
calculation. In addition, extrapolation to the complete basis set (CBS) limit44 andcorrelation energy scaling45,46 are used to further refine the PES shape as well as
accurately reproduce the experimentally determined energetics.
As a further source of dynamical richness, the F(2P) + HCl Æ HF + Cl(2P)
surfaces exhibit conical intersection seams in both entrance and exit channels. Thiscan be explained by elementary chemistry considerations, as described by Hutson andcoworkers47-49 for open shell radical X-HX complexes. This is most readily
understood by first considering the collinear limit, with z as the molecular axis andelectronic states classified in the C∞v symmetry group, followed by bending of thetriatomic away from linearity in the xz plane and consequent reduction to the Cs
symmetry group. Far outside the region of significant chemical reaction, the surfaces
will be asymptotically dominated by electrostatic dipole-quadrupole and dipoleinduced dipole interactions. For a collinear configuration of atoms, this results in anenergy ordering where the Σ state (with the unfilled p orbital parallel to the z-axis)
lies above the doubly degenerate Πx, Πy states (unfilled p-hole perpendicular to the zaxis). However, for the reaction to occur over the lowest barrier, the unfilled p-orbitalmust be approximately parallel to the reaction coordinate to participate in the H atomtransfer. Thus, in the vicinity of the transition state, the lowest energy barrier (again
in C∞v) must arise from a Σ configuration of the F (or Cl) atom with respect to the209
HCl (or HF), thereby reversing the asymptotic Σ and Π ordering of states with respect
to the reactant (or product) channels. Since the symmetries of the Πx and Σ surfaces
are reduced to non-crossing 12A’ and 22A’ surfaces for any non-collinear geometry,
this necessarily results in a 3N-5 = 1-D seam of conical intersections in both theentrance channel and exit channel. As will be shown in this work, these conicalintersection seams for the F + HCl Æ HF + Cl system are energetically quite low(Ecrossing = 3-4 kcal/mol) with respect to the collision energy (Ecom = 5 kcal/mol) as
well as the F(2P) + HCl transition state (ETST ≈ 3.8 kcal/mol). Thus, although sucharguments for the existence of conical intersections are quite generally true for 2P +closed shell systems, these intersections for F + HCl Æ Cl + HF are energeticallyaccessible in all reactive encounters and therefore likely to be of significant relevancein accurately describing the reaction dynamics.
The organization of this chapter is as follows: Section II describes the detailsof the ab initio calculations, including a discussion on the CBS extrapolation andcorrelation energy scaling that is used to provide an accurate PES. In Section III, I
discuss the adiabatic potential energy surfaces, followed by Section IV, detailing thecreation of and analytical fit to the diabatic energies with a discussion on the
inclusion of spin-orbit effects in the diabatic basis. In Section V, I discuss theinteresting aspects of this PES followed by a summary and concluding remarks in
Section VI.

II. Ab Initio BackgroundAll calculations are performed using the MOLPRO suite of ab initio
programs50. Since I explicitly need to account for non-adiabatic as well as single

210
surface effects, electronically excited states of the F(2P) + HCl Æ HF + Cl(2P)
reaction are required. Thus, the ab initio methods employed include multireference
methods to obtain wavefunctions for non-adiabatic calculations and characterizationof the conical intersections. In addition, other 3- and 4-atom systems involving F
have indicated large multireference character in the ground state wavefunctions38,40.
In light of these restrictions and past success with the dynamically weighted
multiconfigurational self-consistent field (DW-MCSCF) method, the natural choicefor this reaction is to use DW-MCSCF along with internally contracted
multireference CI (MRCI) and the multireference Davidson correction (MRCI+Q).
As the MRCI calculations only correlate electrons in orbitals of the same symmerty,
separate adiabatic MRCI+Q calculations are performed for each Cs symmetry: one forthe 12A’ and 22A’ states together and one for the 12A” state.
Residual ab initio error due to incomplete basis sets can be corrected byextrapolating to the complete basis set (CBS) limit44 using the standard aug-cc-pVnZ(n = 2,3,4) basis sets of Woon and Dunning51-53. Based on the benchmark surface
studies of O(3P) + HCl by Peterson44,46, this extrapolation is obtained via the threeparameter fit
2

E (n) = ECBS + Be − ( n −1) + Ce − ( n −1) .

(6.1)

This function yields an excellent extrapolation to the CBS limit as explicitly
confirmed by additional calculations for reagent and product asymptotic geometries
with basis set extensions up to n = 5. The absolute energy difference between the
CBS extrapolation for n = 2 - 4 and n = 2 - 5 was ~0.39 kcal/mol for both the entrance

211
and exit channel asymptotic calculations; a fuller description of the CBS
extrapolation procedure is given in previous work on the ground state PES.38
The adiabatic MRCI+Q calculations are performed using a full valencecomplete active space54 (CAS) reference function with 15 active orbitals (12A’,3A”),
6 of which (5A’,1A”) are uncorrelated. The orbitals are obtained by a state-averaged
multiconfiguration SCF (SA-MCSCF) calculation with 15 active orbitals (12A’, 3A”)
with 7 (6A’, 1A”) closed to excitations. A total of 6 states are then included in the
SA-MCSCF calculation (4A’, 2A”), with the weights determined by the DW-MCSCF
algorithm. In order to account for the triply degenerate p-hole in both the F and Cl
atoms at the asymptote, the SA-MCSCF calculations include all 3 states (2A’, 1A”);
however, low-lying charge transfer (CT) states (e.g. F- + HCl + (2Π) Æ HF + (2Π) + Cl) require the inclusion of 2 more A” and 1 more A’ state (a total of 6 states) to fully
describe the transition state region. As previously demonstrated40, root flipping as theCT surfaces cross other states along the reaction path makes it extremely challenging
to obtain a smooth accurate PES for a fixed set of SA-MCSCF weights. In DWMCSCF, on the other hand, the weights for each state are chosen as a function of its
energy with respect to the ground state energy at each point in configuration space.
The function for the weight based on the energy difference from the ground state iswi ∝ sech 2 (− β (Ei − E 0 ))(6.2)

The decay coefficient (β−1) for DW-MCSCF is empirically chosen to be 3eV. In the
SA-MCSCF calculations, at the transition state, the CT states are about 2eV above theground state, and from our group’s previous F + H2O calculations40 with similarenergetics, it is known that a 3eV DW-MCSCF decay works well. Most importantly,

212
this method permits one to eliminate spurious discontinuities from root flipping andthereby obtain a smooth reaction path and surface while accurately describing both
asymptotic and transition state regions. The orbitals obtained from the DW-MCSCF
calculation are then used as reference states for the internally contracted MRCI
calculation. The benefit of using the DW-MCSCF procure is that the lowest MRCI
states that arise from the DW-MCSCF orbitals are smoothly changing as a function ofreaction coordinate. However, it is impossible to completely eliminate discontinuities
in the upper MRCI states without an unreasonably large decay coefficient. As noted
in previous F + H2O work40, the relative smoothness of the DW-MCSCF energies iscontrolled by the β-1 parameter. A larger β-1, correlates to smoother curves, althoughaccuracy at the asymptotes suffers when β-1 is too large. In this work, the 3eV choiceof β-1 reflects the smoothest 12A’, 22A’, and 12A’’ curves that could be obtained
while still accurately describing the asymptotic energies. By the way of an example,
in the absence of spin-orbit contributions, the three lowest asymptotic F + HCl (andHF + Cl) states should be exactly degenerate. When β-1 is equal to 3eV, these threestates at both asymptotes are degenerate to within 2x10-6 kcal/mol; however, a slightly
larger β-1 of 5eV yields states that are only degenerate to 0.06 kcal/mol.
Building on ideas by Peterson and coworkers46 and Truhlar and coworkers45, I
employ a simple global correlation energy scaling procedure calibrated to match theexperimentally known reaction exothermicity. Specifically, the correlation energy isuniformly scaled by a constant factor near one that reproduces the exothermicity
exactly while retaining the shape of the PES and not introducing any discontinuities.

213
The scaled energy is simply the sum of the MCSCF energy plus the correlationenergy scaled by a factor, γ:
E ( R) SCALED = E ( R) MCSCF + γE ( R) CORR .

(6.3)

In the limit of γ = 0, the scaled energy is just the MCSCF energy (i.e., no extra
correlation energy is recovered), and when γ = 1, the scaled energy is the MRCI+Q
energy (i.e., all correlation energy is recovered). Allowing γ to be greater than one
permits recovery of additional correlation energy, in effect, approximating an increasein active space. The scaling factor (γ) at the CBS limit determined in previous work38
is close to unity (γ = 1.027) for the A’ symmetry; i.e., one only needs to recover an
extra 2.7% of the correlation energy to reproduce the experiment.

III. Adiabatic Surface Calculations
In the F(2P) + HCl Æ HF + Cl(2P) reaction, there are 3 low-lying adiabatic
potential energy surfaces that correlate to the spin-orbit split, F(2P3/2) + HCl andF*(2P1/2) + HCl in the entrance channel and HF + Cl(2P3/2) and HF + Cl*(2P1/2) in theexit channel. The overall calculational strategy is as follows. First of all, in the BornOppenheimer approximation, ab initio methods directly yield the adiabatic energies
for arbitrarily slowly moving nuclei, with non-adiabatic coupling between thesurfaces arising from derivatives of the wave function with respect to 3N-6
displacements. To facilitate later dynamical calculations, however, it proves
convenient to express these (3x3) adiabatic surfaces in an equivalent (3x3) diabatic
representation with off-diagonal coupling terms mixing the 1A’ and 2A’ states. This
choice of diabatic representation is in general non-unique, but can be calculated

214
simply in MOLPRO by constraining the diabats to converge onto the appropriate
adiabatic states for collinear and asymptotic geometries. Secondly, I need to include
spin-orbit interaction due to 2P orbital angular momentum in the F and Cl atoms.
Since this interaction is relatively weak, this is initially neglected this in thecalculation of both adiabatic and diabatic surfaces, but it is later included explicitly
via Breit-Pauli spin-orbit matrix elements55 to generate a 6x6 matrix representation of
the 3 doubly degenerate spin-orbit states.
In order to construct the full 3D potential energy surface, adiabatic
calculations have been performed along a grid of points with the ranges RHF = [1.3 15.0] a0 and RHCl = [1.8 - 15.0] a0 and θF-H-Cl = [180, 150, 120, 90, 60, 30, 0] degrees,
sampling only energies lower than 50 kcal/mol above the entrance channel. The
resulting 3230 points are chosen such that the density is maximized in the mostimportant regions of the surface (i.e., the transition state and van der Walls well
regions), while still sampling a large enough region to get an accurate global
description of the PES. Also included in this full set are 121 points sampling directly
along the reaction path.
Figure 6.1 shows the energies of these three adiabatic surfaces (with the spinorbit interactions included) as a function of the 12A’ reaction coordinate. The lowest
adiabat has a barrier of ~5 kcal/mol, which is lowered to 3.8 kcal/mol when
vibrationally adiabatic zero-point energies in orthogonal coordinates are taken into
account.38 By way of contrast, the 22A’ and 12A’’ surfaces both have much higher
barriers of about 25 kcal/mol. The two high barriers are the result of an avoided
crossing from the relatively low energy CT states that correspond to F- + HCl + (2Π)

215
Æ HF + (2Π) + Cl-, as clearly revealed in the DW-MCSCF calculations. Indeed, it isprecisely the interaction among these states that necessitated the inclusion of up to 6
states in the DW-MCSCF calculation described earlier. A key result of this analysis isthat at experimental collision energies of ~5 kcal/mol, only reactions over the lowest
barrier are energetically permitted. Such adiabatic surface correlations makeunambiguous predictions about the reaction dynamics. For example, in absence of
non-adiabatic surface hopping, spin-orbit excited F would be completely non-reactive
at 5 kcal/mol and any product Cl would be formed exclusively in the ground spinorbit state. As I shall show later, this is largely true but not entirely accurate due tonon-adiabatic coupling between the 12A’ and 22A’ surfaces.
For the special case of 3 atoms in a linear configuration, the three lowest
adiabatic states, 12A’, 22A’, and 12A’’, transform according to three different
symmetries, 12Σ, 12Πx, and 12Πy, respectively. By the Wigner non-crossing rule, the12A’ and 22A’ states couple and avoid for any non-collinear geometry. However,these surfaces can and indeed do cross as 12Σ and 12Πx states when the molecule iscollinear, thereby creating a 1D manifold of conical intersections. As noted in the
introduction, the F(2P) + HCl Æ HF + Cl(2P) reaction exhibits two conicalintersection seams (one in the entrance channel, and one in the exit channel). Figure
6.2 illustrates the nature of the conical intersection region in the exit channel. Whenthe θF-H-Cl is 180o (see Figure 6.2a), there is a 1-D “seam” of geometries (solid blackline) at which the 12Σ and 12Πx states cross (or, alternately, the 12A’ and 22A’ states
touch) and are exactly degenerate. As θF-H-Cl bends away from linearity (Figure 6.2b),
however, the symmetries are reduced from 12Σ and 12Πx to 12A’ and 22A’, thereby

216
avoiding one another in a classic conical intersection.

Of particular dynamicalinterest is that this seam occurs at energies lower than the ground transition statebarrier, and though clearly for a collinear as opposed to the bent transition stategeometry, not dramatically far away from the reaction path. Thus essentially all
reactive collisions (with the exception of tunneling events) have sufficient energy tosample the conical intersection region.
This energetic proximity of adiabatic electronic surfaces implies significant
derivatives in the wavefunctions with respect to motion along the reaction coordinate.
Specifically, when the molecule is bent, the 12A’ ( 22A’) state correlates adiabaticallywith the lower (higher) of the 12Σ and 12Πx states, respectively. As mentioned
previously, when the F atom is far away, the lower state is approximately “Π-like”,
but switches to a “Σ-like“ state as the F approaches HCl. (Σ and Π are not rigorous
quantum labels for anything except the collinear geometry, but nevertheless offer
useful qualitative descriptors.) The region in which the electrons rearrange from Πlike to Σ-like states is the region of the avoided crossing. The time required forelectron rearrangement compared to the time over which the nuclei approach governs
the role of non-adiabatic surface hopping dynamics, as first described in 1-D systems
by Landau and Zener. Simply stated, at sufficiently low energies, the electrons
quickly rearrange to remain on the same adiabatic surface. However, at higher
energies where nuclei move on time scales competitive with electron rearrangement,
non-adiabatic transitions to a different adiabatic surface can result.

Explicit

calculation of these transitions requires a full multiple electronic surface treatment todetermine the role of non-adiabatic contributions to the reaction dynamics, as will be217demonstrated later in Section V. However, these calculations prove to be more
tractable in the corresponding diabatic basis, the extraction of which is discussed in
the next section.

IV. Diabatic Potential Energy Surfaces
A. Diabatic Surface and Non-Adiabatic Coupling Calculations
In contrast to the adiabatic representation, the diabatic state reflects a specificrepresentation for which the wavefunctions do not change rapidly in the vicinity of acrossing. The key advantage of this approach is that off-diagonal coupling between
the diabatic states can be treated as independent of nuclear velocities. Though the
choice of such a representation is not unique, there are several methods for obtaining
well-behaved diabatic electronic and associated coupling surfaces from adiabatic ab
initio calculations. Shatz and coworkers56 have developed a method to approximate
the mixing angle, χ, from matrix elements of the electronic angular momentum
operator.

With this mixing angle, the transformation from adiabatic to diabatic

energies can be described by a simple 3x3 rotation⎡ Σ ⎤ ⎡cos χ − sin χ 0⎤ ⎡ 1A' ⎤
⎢
⎥ ⎢
⎥
⎥⎢
⎢ Π x ⎥ = ⎢ sin χ cos χ 0⎥ ⎢ 2 A' ⎥⎢Π ⎥ ⎢ 001⎥⎦ ⎢⎣ 1A' ' ⎥⎦⎣ y ⎦ ⎣

(6.4)

where 1A' , 2 A' , and 1A' ' are the adiabatic ab initio energies. The corresponding
potential in the diabatic basis is then reported as a 3x3 matrix57:

218
⎡
⎢
Σ
V = ⎢⎢
Π⎢ x⎢⎣ Π y

Σ
VΣ

Πx

β

VΠ x

0

0

β

Πy ⎤
⎥
0 ⎥
0 ⎥
⎥
VΠy ⎥⎦(6.5)

ΠyHere the diagonal terms represent diabatic energies ofΣ , Π x , andrespectively, where β is the coupling surface between Σand Π x states anddirectly related to the mixing angle by 2χ = tan-1[2β/(VΣ−VΠx)]. It is worth stressingthat the labels Σ , Π x , and Π yare rigorously valid only for collinear geometries,

which therefore have a unique axis around which the electronic angular momentum
can be defined. More precisely, these labels denote the corresponding 12A’, 22A’ and12A” diabatic surfaces with which these states correlate upon bending away fromcollinearity.
The mixing angle is determined through MOLPRO subroutines, whichcalculate χ by directly examining the electronic wavefunctions58-60.

Simply stated,

the adiabatic ab initio energies at each geometry are calculated for aug-cc-pVnZ
bases (n = 2,3,4), extrapolated to the CBS limit via Equation 6.1 and further refinedby correlation energy scaling via Equation 6.3. These adiabatic energies are thentransformed to the diabatic basis using the mixing angle obtained from the highest
level aug-cc-pVQZ ab initio calculation. However, sensitivity to basis set level isextremely small, with less than 1%.variation in mixing angles for all three AVnZ (n =2,3,4).
As the specific choice of diabatic representation is non-unique, I use thefollowing method to ensure that the resultant diabatic surfaces are internally self-

219
consistent. At linear geometries, the coupling among the three states must vanish
identically, since all three (12Σ, 12Πx, and 12Πy) are of different symmetries. For each
fixed RHF and RHCl geometry, therefore, I make a 1D angular sweep from θF-H-Cl =180o to 0 o. I start at θF-H-Cl = 180o with diabatic states, Σ , Π x , and Π y , defined
to be the adiabatic states, 1A' , 2 A' , and 1A' ' . In the spirit of the diabatic
representation, calculations at each subsequent non-collinear angle proceed first by
maximizing overlap with the previous orbital, determining the mixing angle via
MOLPRO, and iterating for a series of θF-H-Cl = [120o,90 o,60 o,30 o,0 o] angles. As atest of self-consistency of the diabitization process, the resultingΣandΠ x diabatic states can now in general cross, but must converge smoothly onto theadiabatic states, 1A' , 2 A' , at both θF-H-Cl = 180o and 0o collinear geometries. In
addition, at a few select geometries, θF-H-Cl was swept in 1o increments and comparedto the 30o steps at the same points. Variations in the mixing angle were less than 1%
between the 1o and 30o angular sweeps, indicating excellent convergence of χ.
Figure 6.3 shows sample adiabatic/diabatic results from calculations in theexit channel, where RHF is fixed at the HF equilibrium value (3.4 a0) for RHCl at 6.6 a0.
As expected, the coupling vanishes at both θF-H-Cl = 180o (and θF-H-Cl = 0o), resulting
in identical diabatic and adiabatic energies. At θF-H-Cl = 180o, the Π state is lower inenergy than the Σ state, due to dipole-quadrupole and dipole-induced dipole
interactions (i.e., the p-hole perpendicular to the partially positive H end of the HF).
47-49

As θF-H-Cl sequentially decreases from 180o, the Σ state (p-hole pointing

preferentially parallel to the reaction coordinate) drops in energy due to interaction

220with the partially negative F on the HF. Note that as θF-H-Cl is swept from 180o to 0o,the adiabatic 1A’ and 2A’ states exhibit an avoided crossing in the precise area wherethe Σ and Π states swap energy ordering.
Similarly, diabatic energy crossings (and adiabatic avoided crossings) of the2Σ and 2Π states also occur as a function of the reaction coordinate, as shown in

Figure 6.4a, Note that in the asymptotes, the low-energy states always correlate withthe Π (blue) surface, due to the same electrostatic arguments given above. As thereagents get closer toward the transition state, the Σ-like state must be lower in energy
for the H abstraction to occur.47-49 Notice again that in the region where the diabatic
surface are crossing, the adiabatic surfaces avoid strongly. It is in this region where,
adiabatically, the p-hole is moving from a Π-like state (pointing away from thediatomic) to a Σ-like state (pointing towards the diatomic). Figure 6.4b shows asimilar plot of the reaction path, but this time with the geometry artificiallyconstrained to be collinear (θF-H-Cl = 180o), for which the adiabatic and diabatic
energies are necessarily the same. The adiabats do not avoid in this case because they
are all of different symmetry (2Σ, 2Πx, 2Πy as opposed to the bent geometry case of2

A’ and 2A’’).

B. Spin-orbit Calculations
The spin-orbit terms are calculated identically as described by Alexander andcoworkers61 for the F + H2 Æ HF + H multistate surface. In this method, first theBreit-Pauli spin-orbit operator55 is used to calculate the adiabatic spin-orbit states62
using the unperturbed adiabatic states as a basis. The Breit-Pauli spin-orbit operator

221
generates an adiabatic spin-orbit matrix (HSO), which is then transformed into a
diabatic spin-orbit matrix using the known mixing angle. This procedure results into
two additional A and B surfaces as a function of coordinates (θ, RHF, RHCl), obtained
from the diabatic orbitals via:
A = i Π y H SO Π x(6.6)

B = Π x H SO Σ

(6.7)

and

where the absence (presence) of a bar over a state signifies an electron spin projectionof + 1/2 (-1/2) . After inclusion of spin-orbit terms, the final potential is represented asa 6x6 matrix:
⎡
⎢
⎢ Σ
⎢ Σ
⎢
V = ⎢ Πx
⎢
⎢ Πx
⎢Π
⎢ y
⎢ Πy
⎣

Σ

Σ

Πx

Πx

Πy

VΣ

0

β

B

0

0

VΠ x

B

β

0

β

B

VΠ y

0

A

B

β

0

VΣ

0

0

0

A

0

VΠ x x

0

0

0

−A

0

Πy ⎤
⎥
0 ⎥
0 ⎥
⎥
0 ⎥
⎥
−A⎥
0 ⎥
⎥
VΠ y ⎥⎦(6.8)

The 6 parameters obtained from this adiabatic to diabatic transformation ( VΣ , VΠ x ,VΠ y , β, A, and B) are then fit to analytical functions of θ, RHF, and RHCl as described

below. A FORTRAN subroutine and complete listing of least-squares coefficients forcalculating this matrix and the lowest three diabatic and adiabatic PES’s will be madeavailable though EPAPS63.

222

C. Analytical fits of the diabatic PES
Following the methods outlined by Aguado and Paniagua64, the diabatic
surface energies, VΣ, VΠx, and VΠy, are least-squares fit to a multibody expansion in
pairwise displacement coordinates:

V (RHF , R HCl , RClF ) = V (1) +∑ V (R ) + V ( R , R( 2)i∈{ HF , HCl ,ClF }

( 3)

i

HF

HCl

, RClF ) . (6.9)

where V(1) ( R HF , RHCl , RClF ) represents the energy zero for three separated atoms (i.e.,
subtracted from each ab initio energy), V(2) (Ri) are diatomic potentials for in theasymptotic pairwise limit, and V(3) ( R HF , RHCl , RClF ) reflects the remaining 3-bodyinteractions.
The two body interaction term is based on a sum of modified Rydberg
functions64:
V ( 2 ) (R ) =()

C 0 −α1R Nie+ ∑ C i R e −α 2 R . (6.10)
R
i =1

where R represents the pairwise interatomic distances between the 3 atoms. Ab initio
points for each of the diatomics are fit to a Nth order polynomial with N = 9. For the
Cl + HF reagent channel, 19 ab initio points between RHF = 1.2 a0 and RHF = 3.0 a0
have been used, yielding a least-squares fits for α1, α2 and Ci (i = 0,9) with a global
RMS deviation of 0.039 kcal/mol. Similarly for HCl, 26 points between RHCl = 1.6 a0
and RHCl = 4.0 a0 have been fit to a global RMS of 0.015 kcal/mol. Although the ClF
region is not energetically sampled in reactive scattering experiments at typicalcollision energies, the ClF potential is also obtained for 21 points between RClF = 2.5
a0 and RClF = 4.5 a0, yielding a global diatomic fit with a RMS deviation of ~0.046
kcal/mol.

223
The three body term simply reflects the remaining difference between ab
initio points after subtracting the diatom V(2) and atom V(1) energies. Since this term
must vanish in any of the diatomic asymptotes, a commonly used functional form is
simply a linear combination of two-body product functions, i.e.,M(

V (3) (R1 , R2 , R3 ) = ∑ Dijk R1e − β1 R1( 3)) (R ei2− β 2( 3 ) R2

) (R ej3( 3)− β 31R3),
k(6.11)

i , j ,kwith the total order of the polynomial restricted to i + j + k≤M.

The explicit

numerical least-squares weighting strategy for obtaining this fit is modeled after
efforts developing the Stark-Werner F + H2 surface. First of all, ab initio points along
the reaction path are the most important and given the highest weight of 100. Of theremaining 3230 grid points, the 2200 geometries with energies less than 15 kcal/mol
above the F + HCl entrance channel are weighted 10-fold, with the remaining 1030
points with energies between 15 and 50 kcal/mol receiving a weight of 1. For the VΣ
surface, the resulting global RMS deviation for a least-squares fit with a 12th orderpolynomial (M = 12) is ~0.31 kcal/mol, with RMS deviations for VΠx and VΠysurfaces of ~0.35 kcal/mol and ∼0.34 kcal/mol respectively. However, the errors insuch a fitting process prove to have non-Gaussian tails, for which a simple RMS
metric can substantially overestimate the importance of the deviations. By way of
example, the inset in Figure 6.7a displays the explicit error histogram for the VΣ
surface least-squares fit, which exhibits a substantially improved FWHM metric of~0.049 kcal/mol, i.e., nearly an order of magnitude better than the RMS value.
Furthermore, the surface fit quality for geometries along the reaction coordinate is~0.005 kcal/mol, reflecting improvement by yet another order of magnitude.

224
The above comparison speaks to the presence of residual errors in the leastsquares fitted surface, but in regions that are generally much higher in energy andaway from the reaction path sampled at typical collision conditions. There appear to
be two contributions. The first is due to the strongly exothermic nature of thepotential, for which the above weighting strategy necessarily includes product stateregions high up on the repulsive walls. Secondly, analysis of the residual errors
indicates that regions > 20 kcal/mol above the barrier are less well fit, due specificallyto VΠx and VΠy surfaces in the transition state region. These surfaces are much higher
in energy (~30 kcal/mol) than the VΣ surface (~5 kcal/mol) and are influenced bymultiple and rapid avoided crossings with excited charge transfer states. These effects
can be compensated for by including additional electronic states in the DW-MCSCF
prodecure and/or utilizing a faster decay parameter (~5-6 eV). However, both options
would degrade PES accuracy in the physically relevant regions corresponding totypical collision energies at or near the transition state.

D. Diabatic coupling and spin-orbit fits
Slices of the diabatic coupling, β, are shown in Figure 6.8 for three different
F-H-Cl bend angles, revealing the basic surface topography. Based on theseobservations, the surface is fit to a 3 body-like term given by the sum of i) a simple
polynomial that creates a steep repulsive wall as RHF and RHCl become small and ii) a
direct product of distributed Gaussians in RHF and RHCl and sin(nθ):β (RHF , RHCl ,θ ) =a0 sin(θ )( RHF + RHCl ) 4+ ∑ ain g i ( R HFi , RHCli , σ i ) sin(nθ ) . (6.12)
in

225
As required, the functional dependence of β on sin(nθ) ensures identically zero
coupling between diabats for molecules in the linear configuration (θ = 0, 180). The
first term represents an empirical 4th order polynomial “wall” for non linear
molecules as RHF and RCl approach zero, with an a0 value least-squares fit to the data.
With the large dynamic range terms covered in Equation 6.12, the remaining surface
can then be efficiently fitted by normalized 2D Gaussians, g ( RHFi , RHCli , σ i ) ,centered on ( RHFi , RHCli ) with a width of σi.. The direct product grid of 134 Gaussians
with sin(nθ) (n = 1-3) functions results in a total of 402 basis functions for the leastsquares fit. The RHF, RHCl grid of Gaussians is shown explicitly in Figure 6.6, with σifor each grid point being 0.6411 of the grid spacing. This choice of σi arises from adesire to minimize the coupling (or overlap) of the Gaussian basis functions, whilestill having them wide enough to fill out the space between the grid points. The
weights for the 3230 ab initio points are chosen such that large weights correspond tolarge values of the diabatic coupling. Specifically, guided by simple perturbation
theory, the weight of each ab initio point is taken to be proportional to β/(VΣ-VΠξ)with a cutoff dynamic range of 0.1 and 1,000. The RMS residuals of such a leastsquares fit to the ab initio points are 0.064 kcal/mol on a range of 0 - 15 to kcal/mol.
Similar to the above analysis of the β surface, the two spin-orbit surfaces, A
and B, are fit with a product of distributed Gaussian and angular terms as in Equation
6.12. However since there is now spin-orbit interaction asymmetry between theasymptotic F + HCl and Cl + HF limits, an additional offset term S must be included:

A(orB)(RHF , RHCl ,θ ) = S ( R HCl − RHF ) + ∑ ail g ( R HF i , RHCli , σ i ) Pl (cos(θ i )) , (6.13)
il226
Here g is a normalized Gaussian, Pl is a Legendre function (with l = 0, 1) , and theoffset, S, is in the form of a switching function because the two asymptotes (F + HCl
and HF + Cl) have a constant non-zero value . Since these surfaces must converge on
F and Cl spin-orbit values at the reactant and product asymptotes, respectively, asimple switching function for A is used:
S (ΔR ) =( ACl − AF )( A + AF ),tanh(αΔR ) + Cl
2
2

(6.14)

with a correspondingly similar expression for B. In Equation 6.14, AX is theasymptotic value of i Π y H SO Π x for the X atom asymptote where AF = 124.14 cm1

and ACl = 272.29 cm-1. At the asymptotes, the B spin-orbit term, Π x H SO Σ , isconveniently identical to the A term61, so the corresponding switching functions for A
and B use the same parameters. Based on these correct spin-orbit asymptotic values
for the switching functions plus 202 distributed Gaussian functions in the interaction
region, the RMS residuals of the fit are 3.1 cm-1 for the A surface, and 2.8 cm-1 for the
B surface, both on a range of ~275 cm-1.

V. Discussion
The presence of low-lying conical intersections in the adiabatic surface38
prompts the calculation of diabatic surfaces to be used in full 3-D quantum dynamics
calculations. Figure 6.5b shows the contours (3 kcal/mol spacing) of the collinear
potential energy surface of the 12A’ state along with two dashed lines representing thelocation of the conical intersection, that is, the coordinates where the 12A’ and 22A’

227
states have equal energies at θ = 180o, or, equivalently, where the VΣ and VΠxsurfaces cross at θ = 180o.
Note that conical intersections in both the entrance and exit channels occur atenergies lower than that of the transition state. (Figure 6.5a shows the contours of the
12A’ surface at the transition state angle θ = 123.5o.) Specifically, the barrier toreaction on the 12A’ surface is 3.8 kcal/mol, with the conical intersection energies
dropping as low as 3.55 kcal/mol and -3.23 kcal/mol in the entrance and exit channels
respectively38. Stated simply, this means that any collision with enough energy tocross the barrier and react also has enough energy to access the conical intersection
and interact with the excited 22A’ surface. It is for precisely this reason that each of
the lowest three diabatic surfaces are calculated.
The diabatic surfaces are composed of transformations of the three lowest
adiabatic surfaces. VΣ (Σ-like) and VΠx (Π-like) are mixtures of the 12A’ and 22A’
surfaces while VΠy is identical to the 12A’’ surface, which, by symmetry, does nothave any interaction with the two A’ surfaces. The contours of VΣ, VΠx, and VΠy at
the transition state angle are shown in Figure 6.7. The coupling between VΣ and VΠxis represented in the β surface shown in Figure 6.8 at three angles. At θ = 180o (and

θ = 0o, not shown), the β surface is identically 0 because coupling between thesurfaces is prohibited by symmetry. As seen in Figure 6.8, θ = 150o and 120o, thevalue of the β surface is non-zero and significant in the region of the reaction path.
The magnitude of the coupling is more easily seen in Figure 6.9, where it is plotted asa function of the reaction coordinate. The two coupled diabatic curves are shown
along with the value of β (using the right y-axis) along the reaction path. The

228
coupling (β) is the largest in the steep descent out of the transition state towards thereactants. This is reminiscent of the large non-adiabatic coupling seen in the F + H2O
reaction40,where the OH spin-orbit excited state is observed in greater populations39
than would be predicted in an adiabatic model. In the F + HCl reaction, the two A’
states correlate to the Cl(2P3/2) and Cl*(2P1/2) spin-orbit states asymptotically; thus, asin the F + H2O reaction, non-adiabatic coupling between these two surfaces may
produce an experimentally detectable amount of Cl* even if the reaction only hasenough energy to cross the lowest adiabatic barrier.
To estimate the approximate coupling between the electronic surfaces at the~5 kcal/mol experimental collision energies, I start with the simplest approach anduse the well known Lendau-Zener model65 to predict the probability of surface
hopping along the 1-D F(2P) + HCl Æ Cl(2P) + HF reaction coordinate. In thiselementary calculation, it is assumed that the reaction has just enough energy to getover the lowest adiabatic barrier (similar to the experimental conditions). As thesystem travels along the reaction path, the probability of surface hopping (i.e., exiting
on a different surface than the molecules entered on) can be estimated based on i) the
semiclassical velocity, ii) slope difference between the crossing diabats and iii) thestrength of the coupling.65 The Lendau-Zener model assumes constant velocity andconstant coupling along the reaction path; this is approximated as the values of thevelocity and coupling at the points along the reaction path where the diabats cross.
Since the diabats cross in both the entrance and exit channels, the two crossing
probabilities are considered independently. In the entrance channel, this surface
hopping probability at 5 kcal/mol collision energy is ~0.1. Thus, ~10% of the229
reagents entering on the excited state (F(2P1/2*)) surface with sufficient energy doindeed cross to the lower state and react, with 90% remaining on the upper state
surface and being reflected from the higher adiabatic barrier. Similarly, the exit
channel non-adiabatic contributions can be interpreted by considering all reactions
that cross the lowest barrier. The Landau-Zener expressions predict products moving
over the transition state barrier to have a small but finite probability (5%) of emerging
on the excited state Cl*(2P1/2) surface. Despite the simplicity of the model, this rangeof numbers (5%-10%) is certainly qualitatively consistent with full 3D non-adiabatic
calculations of spin-orbit vs. ground state reactivity for the corresponding F/F* + H2
reaction.
The Landau-Zener model is exact only for i) constant semiclassical velocity
and ii) constant non-adiabatic coupling, which is not rigorously the case for any
realistic molecular system. To get a more quantitative feel for the non-adiabaticity of
the reaction, therefore, I present preliminary results from a more sophisticated model
based on a 1-D wavepacket constrained along the reaction coordinate. Simply stated,
I propagate a 1-D Gaussian wavepacket in time along the 1-D F(2P) + HCl Æ Cl(2P)+ HF reaction coordinate using a standard grid representation66. To propagate a wave
function in time, the relation

Ψ (t 0 + Δt ) = e (−iHΔt )Ψ (t 0 ) , (6.15)
ˆis expanded in a suitable Chebychev polynomial basis67,68 which operates on the wave
function at t = t0 and results in the wave function at t = t0 + Δt. As described by Rushand Wiberg69, the effective mass along the constrained 1D reaction coordinate can bereadily calculated at each point based on large amplitude motion (LAM) Hamiltonian

230
methods of Hougen, Bunker and Johns. An attenuating operator is used70,71 to preventreflection of the wavepacket at the end of the grid, with the probability flux calculated
for each of the electronic surfaces as the wavepacket evolves in time.
By using such wavepacket propagation methods, I am able to include thecoupling between not only between VΣ and VΠx surfaces, but also the additional
coupling with the VΠy surface that occurs upon inclusion of spin-orbit interactions.
Since the potential in Equation 6.8 is a non-Hermetian operator, I switch from theCartesian Σ, Πx, Πy basis to a the signed notation Σ, Π-1, and Π1 using the standard
relations72. In this basis, the potential is a Hermetian block diagonal 6x6 matrix61:
⎡
⎢
⎢ Σ
⎢ Π
⎢ 1
V = ⎢ Π −1
⎢
⎢ Σ
⎢ Π
⎢ 1
⎢Π
⎣ −1
, (6.16)

Σ

Π1

Π −1

− V1 − i 2 B
VΣV1− V1 + i 2 B
VΠ + A
V2
V1
V2
VΠ − A

Σ
0
0
0

0

0

0

VΣ

0

0

0

− V1 − i

0

0

0

V1

Π −1 ⎤
⎥
0
0 ⎥
0
0 ⎥⎥
0
0 ⎥
⎥
− V1 + i 2 B
V1 ⎥2B
VΠ + AV2 ⎥⎥V2VΠ − A⎥⎦Π1

where VΠ = (VΠx + VΠy) / 2, V1 = β / 21/2, and V2 = (VΠy - VΠx) / 2. However, only theupper (or lower) block of the matrix is necessary for wavepacket propagation, withthe remaining elements offering redundant dynamical information on the degenerate
α and β spin-orbit pair73.
By way of demonstration, I present sample results from two wavepacket
calculations on this 1-D grid. The first corresponds to a Gaussian wavepacket started
at the transition state with no momentum (p = 0 au, Δp = 2 au), which then evolves inboth forward and backward directions to reagents and products. Although, the231
calculation was performed in the signed diabatic basis, the results are shown in the
adiabatic basis. At each point the diabatic basis is diagonalized and the resulting
eigenvectors are used to transform the wave packed probability from diabatic to
adiabatic states. The results in Figure 6.10 start at t = 0 with all of the probability on
the lowest adiabatic surface at the transition state.

Subsequent panels reveal

snapshots of |Ψ|2 evolution up to 100 fs, with the bottom panel displaying the total
accumulated flux passing through each point for the full calculation. It is clear thatprobability is coupled primarily between the first and second lowest adiabats, both ofwhich correlate asymptotically with the ground Cl(2P3/2) state. There is transientminor coupling in the entrance well region of wavepacket flux into the third adiabat,but this asymptotically results in only ~1% probability for non-adiabatic transfer tothe spin-orbit excited Cl*(2P1/2) state (Cl:Cl* = 0.993:0.007). It is worth noting thatonce the wavepacket has propagated appreciably away from the transition state, theflux map indicates that probability for each of the adiabats is effectively frozen in.
The issue of non-adiabatic reactivity can be addressed more explicitly bystarting the wavepacket in one spin-orbit asymptote and measuring the flux thatemerges in each of the final spin-orbit states. For wavepackets starting on the F(2P3/2)
surface with enough momentum to cross the barrier (p = 22 au, Δp = 2 au), the total
normalized flux accumulated in the Cl(2P3/2) and Cl*(2P1/2) states state is 0.511 and
0.005, respectively. Thus half of the wavepacket reacts (the remaining probability
reflected back to reactants), again with about ~1% non-adiabatic surface hopping
probability. Conversely, I can start the wavepacket on the spin-orbit excited F*(2P1/2)
state and monitor reaction to form Cl and Cl* product channels. Note that the only

232
way for F*(2P1/2) to react is to undergo a non-adiabatic transition in the entrance
channel, since the barrier that adiabatically correlates with F*(2P1/2) + HCl is ~25
kcal/mol, or about 20 kcal/mol too high for the mean collision energy. As expected,
the non-adiabatic reaction probability is small, corresponding to 2.8% productformation in the Cl(2P3/2) (2.72%) and Cl*(2P1/2) (0.03%) states. In comparison withthe 52% reactivity for wavepackets starting on the ground spin-orbit state surface, thiswould indicate F*(2P1/2) to be ~5% as reactive as F(2P3/2) and in qualitative agreementwith the above Landau-Zener predictions.
This 1D reaction path model also permits mapping out the reactivity ofF(2P3/2) and F*(2P1/2) as a function of energy. Recall from Figure 6.10 that there are
two adiabatic states correlating asymptotically with ground spin order F(2P3/2) andCl(2P3/2), with barriers denoted in energetic order by E1 and E2. Conversely, the single
adiabatic state correlating with spin-orbit excited F*(2P1/2) and Cl*(2P1/2) exhibits abarrier height E3, where E3 ≈ E2 and both E2,E3 >> E1. For these calculations, I
therefore start the reactants in each of the three adiabatic reactant states and computethe total flux in any product (Cl(2P3/2) or Cl*(2P1/2)) exit channels. The results ofthese wavepacket calculations are shown in Figure 6.11a, where solid lines represent
the reaction probability (i.e. sum of the flux in the all of the three exit channels) andvertical dashed lines denote the positions of the three adiabatic energy barriers.
Careful inspection reveals the lowest barrier in this plot is to be E1 ≈ 5.1 kcal/mol, notthe 3.8 kcal/mol value stated earlier. This is because reaction path geometries arecalculated from analytic derivatives on the DW-MCSCF surface, and thus are not in
the exact same place as obtained with MRCI methods (for which MOLPRO does not

233
yet provide analytic gradients). However, this only shifts the reaction thresholds andhas no bearing on the qualitative behavior.
As expected, both reactants correlating with the ground spin-orbit state
F(2P3/2) start yielding significant reaction probability as the energy is increased above
the E1 barrier. The finite (~5-10%) F* reactivity at these low energies is due to nonadiabatic coupling that transfers reactants from the F* surface to (both) F surfaces
prior to the transition state. As the energy is increased further to the higher E2 and E3
barriers, all incoming states become chemically reactive, reaching 100% reactivity
asymptotically. For comparison, also included in Figure 6.11a are dotted lines that
represent purely adiabatic reactivity levels, for a calculation in which the nonadiabatic coupling is arbitrarily set to 0. The difference between dotted and solidlines therefore reflects the non-trivial changes in reactivity due to non-adiabatic
coupling.
It is worth noting that the above wavepacket calculations reflect a noncollinear reaction path and thus do not sample the energetically accessible conicalintersections.

I can therefore take these model calculations one step further toillustrate the possible importance of conical intersections on the reaction dynamics.
By way of example, Figure 6.11b shows the results of the same wavepacket
calculations as in 11a, however in these calculations, the 1-D path the wavepacket is
now constrained to a collinear reaction path. Note that a direct wavepacket samplingof the conical intersection seam results in a preferential reactivity of spin-orbit excited
F* vs F at low energies. This occurs because F* now adiabatically crosses down to
the F states and can therefore react over the lowest barrier. Conversely, a wavepacket

234
starting in one of the F ground spin-orbit states transfers adiabatically to F* at theconical intersection seam and therefore to react must cross one of the two higher
barriers.

Such dramatic differences between collinear versus the full reaction path

dynamics speak to the crucial role of conical intersection seams in the F + HCl
system on controlling the low energy reactivity of F*.

VI. Summary and Conclusion
This chapter presents results for high-level ab initio calculations on the threelowest diabatic electronic states for the F(2P) + HCl → HF + Cl(2P) reaction, alongwith an analytical fit to these surfaces in full 3D and inclusion of spin-orbit effects.
The ab initio data for the analytical fit are calculated at a grid of 3D points using
dynamically weighted MCSCF methods, followed by high-level multireference(MRCI+Q) methods for a series of correlation-consistent basis sets (AVnZ, n =2,3,4). These results are then extrapolated to the complete basis set limit (CBS), witha single global correlation scaling parameter to reproduce the experimentally wellknown exothermicity. The adiabatic surfaces reveal extensive seams of conicalintersections, as predicted from simple orbital considerations for X + HY halogen
atom abstraction processes. For use in multiple surface quantum reactive scattering
dynamical calculations, these surfaces have been systematically converted into a
diabatic representation. The uniqueness of this transformation is maintained byrequiring diabatic and adiabatic surfaces to match at both theta = 0 and 180
configurations, which is equivalent to forcing the coupling matrix elements between
the diabatic surfaces to vanish for collinear geometries. Sequential full 3D fits ofthese surfaces to atom, pairwise and full three body terms have been obtained and235
reported for the three lowest diabatic surfaces, as well as explicit inclusion of spinorbit effects via Breit-Pauli matrix operators that yield excellent agreement withexperiment. By way of preliminary tests, 1-D wavepacket calculations have been
performed along the F(2P) + HCl → HF + Cl(2P) reaction coordinate to estimate theextent of non-adiabatic coupling and reactivity for both ground F(2P3/2) and spin-orbit
excited F*(2P1/2) atomic species. The goal of this work has been a benchmark set of
multistate ab initio potential energy surfaces for the F(2P) + HCl → HF + Cl(2P)system, that should provide much stimulation for further theoretical and experimental
exploration.

236

FiguresFigure 6.1: Three lowest adiabats along the reaction path for the F(2P) + HCl Æ HF
+ Cl(2P) system, obtained by high-level DW-MCSCF, MRCI+Q and CBS methods.
The lowest two correlate with ground state F and Cl, while the highest one correlates
asymptotically with spin-orbit excited F* and Cl* species. In the adiabatic limit, only
the lowest spin-orbit ground state can react over a barrier accessible at experimental
collision energies.

237

Figure 6.2: 2D slices through two full 3D A’ adiabats. a) The region near the conicalintersection seam (solid black line) at θ=180o, which by symmetry becomes b) an
avoided crossing for any non-collinear geometry (θ=170o).

238

Figure 6.3: A 1D angular cut through full 3D adiabatic and diabatic surfaces in theexit channel, with RHF=Req and RHCl=6.6 a0.

The adiabats and diabats are

asymptotically equal at collinear geometries, but strongly avoid as the symmetry is
lowered from C∞v.

239

Figure 6.4: Adiabats (dashed lines) and diabats (solid lines) calculated along thereaction path a) in full 3D and b) with θ constrained at 180o. In full 3D the adiabats
avoid where the diabats cross, whereas for a collinear geometry the adiabats anddiabats overlap perfectly at each point along the conical intersection seam due to zero
coupling.

240

Figure 6.5: A 2D RHF, RHCl contour plot for the ground state adiabatic surface at botha) transition state (θ=123o) and b) collinear (θ=180o) geometries. In b) the conicalintersection seams are shown with dashed lines. Contour spacing is 3 kcal/mol withrespect to zero at reactant entrance channel.

241

Figure 6.6: Grid of distributed Gaussian locations used to fit the β coupling surface.
Angular dependence is fitted at each point (RHF, RHCl) by a linear combination ofsin(nθ) functions.

242

Figure 6.7: 2D contours (at θ=123o) for the a) VΣ, b) VΠx, and c) VΠy diabats with 3
kcal/mol spacing. A sample histogram of residuals for fitting the VΣ function is also
shown in a). The values are not distributed normally, which is why the RMS (0.34
kcal/mol) is considerably larger than the FWHM (0.049 kcal/mol).

243

Figure 6.8: Sample 2D slices through the full 3D β coupling surface between VΣ andVΠx diabats at a) θ=180o, b) θ=150o, and c) θ=120o. Note that for collinear geometry,
β vanishes identically but grows rapidly as the molecule becomes bent.

244

Figure 6.9: VΣ and VΠx diabats along the reaction path (scale to left), and the βcoupling between the diabats (scale to right). Note the strong peaking of diabatic
coupling immediately in the post transition state region, due to rapid change inelectronic character for the newly formed bond.

245

Figure 6.10: Sample results starting at the transition state (p=0 au, Δp=2 au) for
multistate 1D wavepacket propagation along the F + HCl → HF + Cl reaction
coordinate. The top panel displays adiabatic energies, the middle panels reveal
snapshots of |Ψ|2 out to 100 fs, and the bottom panel accounts for the total flux thatpassed by each point on each surface.

246

Figure 6.11: Reaction dynamics for ground F and spin-orbit excited F* as a functionof energy. a) Reaction probability of F and F* along the reaction path (solid lines),
with dotted lines representing probability if non-adiabatic coupling is set to 0. b)
Reaction probability of F and F* along a 2D reaction path with θ constrained to be
180o. The dramatically increased reactivity of F* vs F results from the presence of
conical intersection seams sampled at collinear geometries.

247

References
1

K. Stark and H. J. Werner, J. Chem. Phys. 104, 6515 (1996).

2

M. Hayes, M. Gustafsson, A. M. Mebel, and R. T. Skodje, Chem. Phys. 308,
259 (2005).

3

Y. R. Tzeng and M. H. Alexander, J. Chem. Phys. 121, 5812 (2004).

4

Y. Zhang, T. X. Xie, K. L. Han, and J. Z. H. Zhang, J. Chem. Phys. 119,12921 (2003).

5

S. D. Chao and R. T. Skodje, J. Chem. Phys. 113, 3487 (2000).

6

M. H. Alexander, H. J. Werner, and D. E. Manolopoulos, J. Chem. Phys. 109,
5710 (1998).

7

F. J. Aoiz, L. Banares, B. MartinezHaya, J. F. Castillo, D. E. Manolopoulos,
K. Stark, and H. J. Werner, J. Phys. Chem. A 101, 6403 (1997).

8

J. F. Castillo, D. E. Manolopoulos, K. Stark, and H. J. Werner, J. Chem. Phys.

104, 6531 (1996).
9

S. A. Nizkorodov, W. W. Harper, W. B. Chapman, B. W. Blackmon, and D. J.
Nesbitt, J. Chem. Phys. 111, 8404 (1999).

10

S. A. Nizkorodov, W. W. Harper, and D. J. Nesbitt, Faraday Disc. 113, 107(1999).

11

W. B. Chapman, B. W. Blackmon, S. Nizkorodov, and D. J. Nesbitt, J. Chem.
Phys. 109, 9306 (1998).

12

W. B. Chapman, B. W. Blackmon, and D. J. Nesbitt, J. Chem. Phys. 107,
8193 (1997).

248
13

M. Baer, M. Faubel, B. Martinez-Haya, L. Rusin, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 110, 10231 (1999).

14

M. Baer, M. Faubel, B. Martinez-Haya, L. Y. Rusin, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 108, 9694 (1998).

15

L. Y. Rusin and J. P. Toennies, Phys. Chem. Chem. Phys. 2, 501 (2000).

16

G. Dharmasena, K. Copeland, J. H. Young, R. A. Lasell, T. R. Phillips, G. A.
Parker, and M. Keil, J. Phys. Chem. A 101, 6429 (1997).

17

G. Dharmasena, T. R. Phillips, K. N. Shokhirev, G. A. Parker, and M. Keil, J.
Chem. Phys. 106, 9950 (1997).

18

M. Faubel, L. Rusin, S. Schlemmer, F. Sondermann, U. Tappe, and J. P.
Toennies, J. Chem. Phys. 101, 2106 (1994).

19

D. M. Neumark, A. M. Wodtke, G. N. Robinson, C. C. Hayden, and Y. T.
Lee, J. Chem. Phys. 82, 3045 (1985).

20

D. M. Neumark, A. M. Wodtke, G. N. Robinson, C. C. Hayden, K. Shobatake,
R. K. Sparks, T. P. Schafer, and Y. T. Lee, J. Chem. Phys. 82, 3067 (1985).

21

H. Kornweitz and A. Persky, J. Phys. Chem. A 108, 8599 (2004).

22

Y. R. Tzeng and M. H. Alexander, J. Chem. Phys. 121, 5183 (2004).

23

Y. R. Tzeng and M. Alexander, Phys. Chem. Chem. Phys. 6, 5018 (2004).

24

Y. Zhang, T. X. Xie, K. L. Han, and J. Z. H. Zhang, J. Chem. Phys. 120, 6000
(2004).

25

T. X. Xie, Y. Zhang, M. Y. Zhao, and K. L. Han, Phys. Chem. Chem. Phys. 5,2034 (2003).

26

S. C. Althorpe, J. Phys. Chem. A 107, 7152 (2003).

249
27

S. C. Althorpe, Chem. Phys. Lett. 370, 443 (2003).

28

S. H. Lee, F. Dong, and K. P. Liu, J. Chem. Phys. 116, 7839 (2002).

29

K. P. Liu, Nuc. Phys. A 684, 247C (2001).

30

R. T. Skodje, D. Skouteris, D. E. Manolopoulos, S. H. Lee, F. Dong, and K. P.
Liu, Phys. Rev. Lett. 85, 1206 (2000).

31

R. T. Skodje, D. Skouteris, D. E. Manolopoulos, S. H. Lee, F. Dong, and K.
Liu, J. Chem. Phys. 112, 4536 (2000).

32

W. W. Harper, S. A. Nizkorodov, and D. J. Nesbitt, J. Chem. Phys. 116, 5622
(2002).

33

A. M. G. Ding, L. J. Kirsch, D. S. Perry, J. C. Polanyi, and J. L. Schreiber,
Faraday Disc. Chem. Soc. 55, 252 (1973).

34

R. Sayos, J. Hernando, R. Francia, and M. Gonzalez, Phys. Chem. Chem.
Phys. 2, 523 (2000).

35

B. Y. Tang, B. H. Yang, K. L. Han, R. Q. Zhang, and J. Z. H. Zhang, J. Chem.
Phys. 113, 10105 (2000).

36

R. Sayos, J. Hernando, J. Hijazo, and M. Gonzalez, Phys. Chem. Chem. Phys.

1, 947 (1999).
37

H. Kornweitz and A. Persky, J. Phys. Chem. A 108, 140 (2004).

38

M. P. Deskevich, M. Y. Hayes, K. Takahashi, R. T. Skodje, and D. J. Nesbitt,
Journal Of Chemical Physics 124, (2006).

39

M. Ziemkiewicz, M. Wojcik, and D. J. Nesbitt, J. Chem. Phys. 123, 224307(2005).

250
40

M. P. Deskevich, D. J. Nesbitt, and H. J. Werner, J. Chem. Phys. 120, 7281
(2004).

41

H. J. Werner and P. J. Knowles, J. Chem. Phys. 89, 5803 (1988).

42

P. J. Knowles and H. J. Werner, Chem. Phys. Lett. 145, 514 (1988).

43

P. J. Knowles and H. J. Werner, Theoretica Chimica Acta 84, 95 (1992).

44

K. A. Peterson, D. E. Woon, and T. H. Dunning, J. Chem. Phys. 100, 7410
(1994).

45

F. B. Brown and D. G. Truhlar, Chem. Phys. Lett. 117, 307 (1985).

46

B. Ramachandran and K. A. Peterson, J. Chem. Phys. 119, 9590 (2003).

47

M. Meuwly and J. M. Hutson, J. Chem. Phys. 119, 8873 (2003).

48

M. L. Dubernet and J. M. Hutson, J. Chem. Phys. 101, 1939 (1994).

49

M. L. Dubernet and J. M. Hutson, J. Phys. Chem. 98, 5844 (1994).

50

H.-J. Werner, P. J. Knowles, R. Lindh, M. Schutz, P. Celani, T. Korona, F. R.
Manby, G. Rauhut, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper,
M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, A. W.
Lloyd, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R.
Pitzer, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni, and T. Thorsteinsson,
Molpro,version2002.6,apackageofabinitioprograms(see

http://www.Molpro.Net) (Birmingham, UK, 2003).
51

D. E. Woon and T. H. Dunning, J. Chem. Phys. 98, 1358 (1993).

52

T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).

53

R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. 96, 6796
(1992).

251
54

L. M. Cheung, K. R. Sundberg, and K. Ruedenberg, Int. J. Quantum Chem.

16, 1103 (1979).
55

A. Berning, M. Schweizer, H. J. Werner, P. J. Knowles, and P. Palmieri,
Molecular Physics 98, 1823 (2000).

56

A. J. Dobbyn, J. N. L. Connor, N. A. Besley, P. J. Knowles, and G. C. Schatz,
Physical Chemistry Chemical Physics957.

57

M. H. Alexander, J. Chem. Phys. 99, 6014 (1993).

58

H. J. Werner and W. Meyer, J. Chem. Phys. 74, 5802 (1981).

59

H. J. Werner, B. Follmeg, and M. H. Alexander, J. Chem. Phys. 89, 3139
(1988).

60

D. Simah, B. Hartke, and H. J. Werner, J. Chem. Phys. 111, 4523 (1999).

61

M. H. Alexander, D. E. Manolopoulos, and H. J. Werner, J. Chem. Phys. 113,11084 (2000).

62

A. Berning, M. Schweizer, H. J. Werner, P. J. Knowles, and P. Palmieri,
Molecular Physics 98, 1823 (2000).

63

See epaps document no xxx for a complete listing of the functional form andleast squares coefficients for the pes. This document can be reached via a
direct link in the online article's htms reference section or via the epaps
homepage (http://www.Aip.Org/pubservs/epaps.Html)

64

A. Aguado, C. Tablero, and M. Paniagua, Comput. Phys. Comm. 108, 259(1998).

65

M. S. Child, Molecular collision theory. (Academic Press, New York, 1974).

66

D. Kosloff and R. Kosloff, Journal Of Computational Physics 52, 35 (1983).

252
67

H. Talezer and R. Kosloff, Journal Of Chemical Physics 81, 3967 (1984).

68

J. R. Fair, D. Schaefer, R. Kosloff, and D. J. Nesbitt, J. Chem. Phys. 116,
1406 (2002).

69

D. J. Rush and K. B. Wiberg, Journal Of Physical Chemistry A 101, 3143
(1997).

70

V. A. Mandelshtam and H. S. Taylor, Journal Of Chemical Physics 103, 2903
(1995).

71

V. A. Mandelshtam and H. S. Taylor, Journal Of Chemical Physics 102, 7390
(1995).

72

R. N. Zare, Angular momentum. (John Wiley & Sons, New York, 1988).

73

M. Tinkham, Group theory and quantum mechanics. (McGraw-Hill, New
York, 1964).