Three-dimensional potentials ˜ state of the OH–He complex

Three-dimensional potentials
of the X̃ 2Π state of the OH–He
complex
Master Thesis
Dick Tanis
May, 2012
Institute for Molecules and Materials
Summary
Hydroxyl (OH) molecules play an important role in atmospheric chemistry and
astrophysics. To understand this role it is important to investigate the collision of
OH radicals with other species. Especially the collisions of OH with rare gases
have been extensively investigated. However, the quantum mechanical effects of
vibrationally excited OH radicals, colliding with rare gas atoms, are yet virtually
unexplored. Therefore, this research will make a start to study the effect of vibrational excitations of OH when colliding with helium (He) atoms. The aim of this
work is to develop three-dimensional potential energy surfaces (PESs) of the X̃ 2 Π
state of the OH–He complex which could be used to study the effect of vibrational
excitations on the collision dynamics between OH and He.
For the development of the potential, ab initio points have been calculated
using the open-shell spin-restricted coupled cluster [RCCSD(T)] method. The
resulting diabatic surfaces have successfully been fitted by an analytical function
and have recently been used for quantum mechanical inelastic scattering calculations of OH (X 2 Π3/2 , J = 3/2, f ) with He. The results of these calculations show
that our potential is realistic and that the calculated relative cross sections from
our potential have improved the agreement with experimental data. The developed
potential and the results of these scattering calculations have been published [1].
iii
Dankwoord
Als eerste wil ik mijn begeleider Gerrit Groenenboom bedanken voor alle deskundige hulp en uitleg tijdens mijn Master Stage. Verder dank ik hem voor het
vertrouwen dat hij in me gesteld heeft.
Er zijn diverse mensen die direct of indirect aan de totstandkoming van dit verslag hebben meegeholpen. Bij deze bedank ik Koos en Liesbeth voor de begeleiding bij het schrijven van dit verslag. Ook wil ik Koos in het bijzonder bedanken
voor de samenwerking die er toe geleid heeft dat de ontwikkelde potentiaal een
onderdeel is geworden van een publicatie. Verder wil ik naast de hier al eerder genoemde personen, Ad, Anna, Eugenie, Mark, Paul, Rob en Vivike bedanken voor
de gezellige en leerzame periode bij de afdeling Theoretische chemie.
Mijn vrienden Arjan, Boudewijn, Gerbe en Tom wil ik speciaal bedanken voor
de steun en de hulp die ze me hebben geboden de afgelopen jaren.
Aangezien na de oplevering van dit verslag mijn studie aan de Radboud Universiteit Nijmegen (RUN) bijna heb voltooid, wil ik bij deze ook mijn mede RUNners Audry, Bart, Fieke, Harm, Irmgard, Jordy, Kors, Marcel, Martien, Michiel,
Mirjam, Nearchos, Pieter, Sanne, Stefanie, Teun en Willem bedanken voor de gezellige studie tijd.
Tot slot bedank ik mijn ouders voor wat zij hebben gedaan om mijn studie
mogelijk te maken.
v
Contents
Summary
iii
Dankwoord
v
1 Introduction
1.1 Potential energy surface of OH–He . . . . . . . . . . . . . . . . .
1.2 Earlier work of OH–He PESs . . . . . . . . . . . . . . . . . . . .
1.3 About this work . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
2
3
3
2 Ab initio calculations
2.1 The OH–He system . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Electronic structure calculations . . . . . . . . . . . . . . . . . .
2.3 RHF breakdown . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
5
6
7
3 The interaction potential
3.1 The grid setup . . . . . . . . . .
3.2 The analytic form . . . . . . . .
3.3 The fit procedure . . . . . . . .
3.3.1 The sum potential . . . .
3.3.2 The difference potential
3.4 Results . . . . . . . . . . . . . .
.
.
.
.
.
.
9
9
10
11
11
14
15
4 Inelastic scattering of OH by He
4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 QM calculations and results . . . . . . . . . . . . . . . . . . . . .
21
21
22
5 Conclusions
25
Bibliography
27
A MOLPRO script RCCSD(T) calculations
33
B FORTRAN code OH–He potential
37
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
vii
1
Introduction
Hydroxyl (OH) molecules are very important intermediates in atmospheric chemistry. Without them, the composition of the atmosphere would be totally different
and hazardous to life on Earth. The OH molecule is a radical and because of the
high reactivity, it transforms attacked molecules to increased oxidation and solubility [2, 3]. For this reason the OH radical is sometimes referred to as the atmosphere
cleaning agent. In the upper stratosphere and mesosphere, highly vibrational (up
to ν = 9) OH radicals are produced via reactive depletion of ozone [4–8]. The excited radicals are responsible for the near-infrared night-time air-glow [9]. These
emissions are called the Meinel bands and have also been observed from artificial aurora at higher altitudes, offering new possibilities to study interactions in
the ionosphere [10]. Another field where OH plays an important role is in astrophysics. For example, radio astronomers measure the stimulated emission of radio waves (MASER) by OH molecules near supernovae, to investigate supernova
shock waves [11]. Other applications of OH MASERs include the investigation
of pulsars [12], star-forming regions [13] and envelopes of late-type stars [14]. In
our own solar system comets also produce highly excited radicals. When they approach the sun, water vaporizes, which dissociates under influence of sunlight to
produce OH radicals. Some of these radicals will radiate and their emissions are
used as a tracer for water [15].
In atmospheric and interstellar chemistry, the OH radical has many collision
partners which are interesting to investigate. For example, in the mesosphere
the highly vibrationally excited radical is mainly deactivated by N2 , O2 , and O,
[4, 8, 16–18] and in interstellar OH MASERs, the required population inversion is
induced by particles like He and H2 [19]. There are many experimental studies on
the deactivation rates of OH by scattering with its collision partners but the more
detailed experimental studies on scattering on the quantum level, have been done
for OH with simple atoms. Especially the collisions of OH with rare gases like He,
Ar, and Xe have been extensively investigated [19–27] because these systems have
emerged as paradigms of scattering of an open-shell molecule by an atom. These
experimental and theoretical studies of scattering of OH with rare gas atoms mainly
focus on energy transfer between rotational levels. The quantum mechanical (QM)
effects of vibrationally excited OH radicals colliding with rare gas atoms are yet
virtually unexplored. Therefore, this research will make a start to study the effect
1
CHAPTER 1. INTRODUCTION
of vibrational excitations on the collision dynamics. In this thesis we will use the
OH–He system as our model system where the OH molecule is in the electronic
ground state.
1.1 Potential energy surface of OH–He
In order to perform QM scattering calculations on the OH–He system, potential
energy surfaces are required because these govern the dynamics. The potential energy surface (PES) is the interaction energy of a system as a function of all the
coordinates describing the relative position of molecules or atoms in the system.
For various purposes it is easy to think the PES in three dimensions: a vertical
dimension for the energy, and the two horizontal dimensions which are representative of the coordinates [28]. In our system, we have 3 coordinates which we will
discuss later in chapter 2. Furthermore, the PES can be seen as a landscape which
may contain minima, maxima, and saddle points. On this surface there is always a
‘global minimum’ that is lower than any of the other minima. Sometimes, due to
symmetry the surface contains several global minima of equal energy. The other
minima are called ‘local minima’ [28]. The global minimum can be seen as the
absolute ground state and the other minima are the metastable states.
The OH(2 Π)–He complex can be seen as an OH(2 Π) radical which is perturbed
by a helium atom. In our system the OH radical is in the X 2 Π ground state. Therefore, the state of our complex is denote by X̃ 2 Π. The ground state of the OH
radical has 2 Π3/2 electronic symmetry and the electron configuration is nominally
1σ 2 2σ 2 3σ 2 1π 3 [29]. Due to the presence of the perturbing He atom the cylindrical
symmetry C∞,v of OH is lowered to Cs symmetry and this results in two states
of A′ and A′′ symmetry. In this system the OH molecule and the He atom lie in
the Cs mirror plane. In the A′ state, the singly occupied π orbital lies in the Cs
mirror plane and the doubly occupied π orbital is perpendicular to the plane (see
Fig. 1.1a). In Fig. 1.1b we see the A′′ state where the doubly occupied π orbital lies
FIG. 1.1: The electronic configuration of the two states of the OH(2 Π)–He complex. Panel (a)
displays the A′ symmetry state and panel (b) the A′′ symmetry state. The OH molecule and the He
atom lie both in the plane. The gray dashed line is the interatomic axis of the OH molecule. The
black arrows are the electrons of the π orbitals and have a different configuration in the two states.
2
1.2. EARLIER WORK OF OH–HE PESS
in the plane and the singly occupied π orbital is perpendicular to the plane. Hence,
two PESs are needed for a complete description of the collision dynamics between
OH and He.
1.2 Earlier work of OH–He PESs
Earlier studies have already reported the development of potentials for the OH–He
complex. In particular, Lee et al. [30] developed two-dimensional (2D) PESs for
the X̃ 2 Π state and the low-lying excited Ã2 Σ+ state of the OH–He complex, using
the restricted open-shell coupled cluster [RCCSD(T)] method. In their calculation
they set the OH bond length to a fixed value which is equal to the experimentally
determined equilibrium length of 1.8324 a0 . These 2D potentials have been used
in various theoretical studies. For example, various cross sections for rotational
inelastic scattering calculated by Kłos et al. [19] and Scharfenberg et al. [27] are
in good agreement with experimental observations. Dagdigian et al. [25] computed tensor cross sections for collisions of OH(X 2 Π) with helium and predicted
rate constants for collisional depolarization of specific rotational fine-structure levels. These constants were in good agreement with spectroscopy experiments of
Paterson et al. [31]. Also bound states involved in the A2 Σ+ -X 2 Π transition of
OH in the OH–He complex were calculated. The level of agreement between the
predicted resonances and observed resonances also give a good indication that the
ground-state potential of Lee et al. [30] is realistic [32].
In this study we want to describe vibrationally excited states accurately, therefore the bond length of OH needs to be varied. This will produce three-dimensional
(3D) potentials which depend on the 3 coordinates of the OH–He system. Degli
Esposti et al. [33] also developed a 3D version of the X̃ 2 Π surface, using the
coupled electron pair approximation (CEPA). They used only three OH distances
around the equilibrium length to develop potentials which were vibrationally averaged for further calculations. These 3D potentials are also unsuitable for describing
vibrationally excited states because they only describe the ground state of OH near
equilibrium and are less accurate than the RCCSD(T) potentials. As mentioned
earlier, the 2D potentials of Lee et al. [30] are of high quality. So in this study we
will use the same method to create 3D OH(2 Π)–He potentials as it will give a good
indication that our potentials will also produce realistic results.
1.3 About this work
The aim of this work is to develop 3D PESs of the X̃ 2 Π state of the OH–He complex which could be used to study the effect of vibrational excitations on the collision dynamics between OH and He. This thesis also presents the results of calculations on resonances in rotationally inelastic scattering of OH (X 2 Π3/2 , J = 3/2, f )
with He. The cross sections, calculated from the newly developed 3D potentials
3
CHAPTER 1. INTRODUCTION
compare favorably with the recent crossed beam scattering experiments of Kirste
et al. [34].
This thesis is structured as follows. Chapter 2 describes the procedure and the
method used to the calculate the ab initio points of the A′ and A′′ surface. Chapter 3
describes the fitting procedure of the 3D potential for the diabatic states and reports
the results of the fitted A′ and A′′ PESs. Chapter 4 presents the first results of inelastic scattering calculations which are based on our developed potential. The results
of these QM calculations are also compared with experimental results. Finally,
chapter 5 reports the conclusions of this thesis.
4
2
Ab initio calculations
2.1 The OH–He system
The coordinates of the OH–He complex are given in Fig. 2.1. The system is defined
by the vector r with length r (OH bond length), vector R (center of mass of OH to
He) with length R, and the angle θ between the latter vectors. The expression of
these Jacobi coordinates are
r = rH − rO
mH rH + mO rO
R = rHe −
mO + mH
r·R
,
θ = arccos
rR
(2.1)
(2.2)
(2.3)
when rH , rO , rHe are the positions and mH and mO are the masses of the individual
atoms. Values of θ = 0◦ and θ = 180◦ correspond to the linear OH–He and He–OH
configurations of the complex, respectively.
FIG. 2.1: The coordinates for the OH–He complex. The geometry of the system is defined by 3
coordinates R, θ, and r. The extra coordinates Ra , θa , Rb , θb are used to describe the short range
behavior of the interaction potential near the H and O atom. The point “X” marks the location of the
bond orbitals used in the ab initio calculations.
5
CHAPTER 2. AB INITIO CALCULATIONS
We also define the vector Ri with length Ri and θi (i = a, b) as
Ra = rHe − rH
r · Ra
θa = arccos
rRa
Rb = rHe − rO
r · Rb
.
θb = π − arccos
rRb
(2.4)
(2.5)
(2.6)
(2.7)
These coordinates are convenient to describe the short range behavior near the H
and O atom separately and are used for the expression of the analytic expansion of
the PESs (see chapter 3).
2.2 Electronic structure calculations
The OH–He complex is a van der Waals complex because the interaction of the OH
radical with helium is only a few tenths of cm−1 . This weak binding arises mainly
from dispersion forces. To account adequately for these, electron correlation must
be described as accurately as possible. Therefore, we use the spin-restricted openshell coupled cluster approach by Knowles et al. [35, 36]. This method takes the
molecular orbitals from the restricted Hartree–Fock (RHF) calculation to construct
multi-electron wavefunctions using an excitation cluster operator. In our calculations the cluster operator is restricted to single and double excitations and perturbation theory is used to estimate the effect of triple excitations [RCCSD(T)].
For our calculations we use the triple-zeta correlation-consistent polarized basis set (aug-cc-pVTZ) [37–39], augmented with an additional (3s3p2d2f 1g) set of
bond functions to recover the dispersion energy of the van der Waals interaction.
These bond orbitals are placed on the vector R at a distance of min(Ra , Rb )/2
from helium (see Fig. 2.1) while the exponents are geometry dependent using the
procedure as described in Ref. [40]. The center exponent for the s and p functions
is defined as αb = 0.34(5.3/min(Ra , Rb ))2 while the other s and p exponents
are defined as αa = 2.7647αb and αc = 0.3530αb . We set the exponents of the
d and f functions to βa = 1.882αb and βb = 0.6762αb . For the g function we
choose γ = 1.0293αb . Note that the exponent-ratio of all bond functions is equal
to the one described in Ref. [41]. This procedure prevents the bond functions from
overlapping to much with the atom orbitals for small values of R.
All calculations are performed using the MOLPRO 2002 package [42]. The
interaction energies of the system are determined via the supermolecule method
as V (R, θ, r) = EOH−He − EOH − EHe , where the energies of the separating
fragments are computed in the same one-electron basis as the complex to avoid
introducing a basis set superposition error (BSSE) [43]. Appendix A reports the
general MOLPRO script used for the various energy calculations.
6
2.3. RHF BREAKDOWN
2.3 RHF breakdown
As already mentioned in section 2.2, for the ab initio calculation, we initially generate RHF molecular orbitals which are subsequently used for the coupled cluster
calculation. The coupled cluster method is a single-reference method because it
involves only one reference wavefunction. When r becomes large and the OH
bond breaks, the RHF wavefunction which consists of a single Slater determinant
is not able to describe the separate O and H atoms anymore. This is easily seen
in Fig. 2.2 where a single-reference and a multi-reference interaction potential of
OH(X 2 Π) are shown. At large r these interaction potentials should go to zero
which is the case for the multi-reference CASSCF+MRI method but at r > 4.6 a0
the single-reference RHF+CCSD(T) potential becomes unphysical. Figure 2.2 also
reports the T1 diagnostic of the RHF+CCSD(T) calculations using the open-shell
T1 formalism of Jayatilaka and Lee [44]. This T1 diagnostic value is a qualitative
measure of multireference character: the larger the T1 value, the less reliable the
results of the single-reference coupled cluster wave function. This T1 open-shell
diagnostic is consistent with the formalism of Lee and Taylor [45] for closed-shell
systems. Lee and Taylor suggest that values greater than 0.02 indicates the need for
a multi-reference wavefunction. However, Jayatilaka and Lee suggest that open6
0.08
0.05
0
4
0.07
−0.05
0.06
−0.1
−0.15
0.05
4
5
6
7
0
0.04
X2Π
−2
0.03
RHF+CCSD(T)
CASSCF+MRI
T [RHF+CCSD(T)]
−4
1
−0.2
T diagnogstics
energy (eV)
2
0.02
0.01
1
−6
1
1.5
2
2.5
3
r (a )
3.5
4
4.5
0
5
0
FIG. 2.2: Interaction potentials of the OH X 2 Π state with different calculation methods. The color
code is as follows: dot dashed blue: RHF+CCSD(T), AVTZ basis, dashed red: CASSCF+CI AV6Z
basis, solid black: T1 diagnostic of RHF+CCSD(T) calculations. Note that for r > 5.39 a0 the
RHF+CCSD(T) data is not plotted because this method does not convergence in that region.
7
CHAPTER 2. AB INITIO CALCULATIONS
shell T1 values may be larger than those of closed-shell systems. Kiracofe et al.
[46] reported that QM calculations with T1 values up to 0.044 gave reliable results
for their system. In our case the T1 diagnostic is 0.02 at r = 3.1 a0 and becomes
greater than 0.044 for r > 4.3 a0 . Hence, our RCCSD(T) calculations are definitively reliable up to r = 3.1 a0 and at r = 4.3 a0 the RCCSD(T) method breaks
down.
8
3
The interaction potential
3.1 The grid setup
For the development of the 3D interaction potential we have to define a grid which
depends on the coordinates R, θ, and r (see Fig. 2.1). A point on this grid represents the interaction energy of the OH–He system at Rℓ , θj , rk . For the short and
the intermediate range, where the interaction potential can vary substantially, we
require a dense grid to describe the potential energy surface accurately. Whereas
for the long range asymptotic behavior, a less dense grid will be sufficient.
For the OH interatomic distance r we used a grid from 0.75 to 4.5 (all distances are given in a0 ). Note that in the previous chapter, we concluded that the
RCCSD(T) method to calculate the interaction energy breaks down for r > 4.3.
We keep this in mind when performing the fit. We come back to the range of
validity for our obtained fit later on.
From Fig. 3.1 we see that the probability density for the ν = 11 vibration
wavefunction for r > 4.3 is less than 10−5 . Hence, we expect our potential to be
0
10
log probability
−2
−4
−6
−8
−10
1
2
3
r (a0)
4
5
ν=0
ν=1
ν=2
ν=3
ν=4
ν=5
ν=6
ν=7
ν=8
ν=9
ν=10
ν=11
ν=12
FIG. 3.1: The probability of the vibrational states, calculated with the monomer potential of OH in
the X 2 Π state. This potential is computed with RCCSD(T) and aug-cc-pVTZ as basis.
9
CHAPTER 3. THE INTERACTION POTENTIAL
suitable for collisional studies of OH–He up to at least ν = 11.
The Jacobi angle θ only needs to be varied between 0◦ and 180◦ . We distinguish two regions of R:
1. The short and intermediate range which has an interval of 3 ≤ R ≤ 12.5.
We used an equally spaced grid with ∆R = 0.25, ∆r = 0.25, and ∆θ = 12◦ ;
2. The long range which has an interval of 12.5 < R ≤ 20. We used 4 logarithmically spaced points for R. For r an equally spaced grid with ∆r = 0.50 is
used which is cut off at r = 4.5. The grid of the Jacobi angle θ was equally
spaced with 9 points.
In total for 6382 geometries the interaction energy was calculated. Points with total
energy of 2.6896 Eh above the energy of OH at equilibrium distance and points
where the interaction energy was larger than 0.1 Eh , were not used in the fit.
3.2 The analytic form
For scattering calculations it is convenient to use diabatic potentials which can be
expanded in spherical harmonics. The diabatic potentials Vsum and Vdiff are the
sum and difference potentials [47] of the two adiabatic potentials VA′ and VA′′ and
are expressed as
1
Vsum (R, θ, r) = [VA′′ (R, θ, r) + VA′ (R, θ, r)]
2
(3.1)
1
Vdiff (R, θ, r) = [VA′′ (R, θ, r) − VA′ (R, θ, r)],
2
(3.2)
which both are written as
X
V (R, θ, r) =
Vsr (Ri , θi , r) + Vlr (R, θ, r).
(3.3)
i=a,b
This representation is convenient, because the coordinates Ri and θi , i = a, b (see
Fig. 2.1) are ideally suited to describe the short range behavior near the H and
O atom, respectively, while the coordinates of the second term are convenient to
describe the long range behavior. The short range of the potential is expanded as
lmax1
Vsr (Ri , θi , r) =
X
(i)
e−βi Ri Pl (cos θi )sl
l=0
nmaxi kmaxi lmaxi
+
X X X
Rin e−βi Ri Pl (cos θi )
n=0
×
10
k=0 l=0
k −αi r 3 (i)
r e
snlk .
(3.4)
3.3. THE FIT PROCEDURE
TABLE 3.1: The limits of the summation terms used in Eq. (3.4) for the sum potential and the
difference potential.
Vsum
1
3
7
8
3
5
8
lmax1
nmaxa
lmaxa
kmaxa
nmaxb
lmaxb
kmaxb
Vdiff
2
5
6
5
5
6
4
For the long range the expansion is
13 n−4
X
X fn (βR)
Pl (cos θ)cInl (r)
Vlr (R, θ, r) =
Rn
n=10 l=0
9 n−4
X
X
+
n=6 l=0
fn (βR)
Pl (cos θ)cII
nl (r).
Rn
(3.5)
Only terms with l + n is even are included and
cinl (r) = cinl +
3
X
i 3
r k e−α r cinlk ,
for i = I, II.
(3.6)
k=0
The functions fn are Tang-Toenies damping functions [48] defined by
−x
fn (x) = 1 − e
n
X
xk
k=0
k!
,
(3.7)
which suppress the singularity of the long range function as R → 0. Since the A′
and A′′ surfaces are degenerate in linear geometries the difference potential is zero
at θ = 0◦ , 180◦ . We note that for Vdiff similar expressions hold. The only difference with Vdiff is that the Legendre polynomials Pl must be replaced by associated
(2)
Legendre functions Pl [49, 50]. Note that both types of Legendre functions are
Racah normalized in our expansion. The upper summation limits in Eq. (3.4) of
both diabatic potentials are described in section 3.3.1 and 3.3.2. These limits are
tabulated in table 3.1.
3.3 The fit procedure
3.3.1 The sum potential
Initially, we followed the global fit procedure as described in Ref. [51]. However,
fitting all points at once in the short and intermediate regions gave problems, so we
11
CHAPTER 3. THE INTERACTION POTENTIAL
adapted the procedure. This was done as follows. We start with the long range and
fit the ab initio points with R ≥ 16 a0 using only the n = 6, 7, 8, and 9 terms of
the long range function Vlr [see Eq. (3.5)]. The damping function fn is set to one
II
and the coefficients cII
nl and cnlk are determined in a weighted linear least squares
(WLLS) fit using the weight function w∞ (R) = R6 . In the next step the fitted
long range is substracted from all the ab initio points. For the resulting surface we
switch to a sequential fitting method which is described below in this section.
As we want to fit ab initio points that vary from the order of 1 µEh around
R = 15 a0 to about 0.1 Eh in the repulsive region we construct a weight function
w(R, θ, r) such that w(R, θ, r)V (R, θ, r) is in the order of 1 everywhere. In the
short and long range the function w = |V |−1 would work well, but in the intermediate range the interaction potential may go through zero. So, we factorize the
weight function w = wsr × wlr and choose factors such that

V0

for small R(V → ∞)

V
wsr ∼
(3.8)


1
for large R(V → 0),
wlr ∼
We take
 1



 V0
for small R(V → ∞)
(3.9)

6


 R
C6
for large R(V → 0).
wsr = [ln(eV /V0 + e − 1)]−1
"
wlr = 1 +
Ri
R0
6 #
V0−1 ,
(3.10)
(3.11)
where Ri = min(Ra , Rb ), Ra and Rb are defined in Fig. 2.1 and V0 = C6 /R06 .
The parameter V0 determines where the short range factor of the weight function
“switches on”. We set it equal to V0 = 5|E0 |, where E0 = −203 µEh which is the
most attractive point on the grid. For C6 we take the rounded value of C6 = 7 (see
section 3.4) which gives R0 = 4.36 a0 .
To fit the remainder of the sum potential we follow a three-step procedure. In
the first step, for each θj and rk of our grid we use
V (R, θj , rk ) =
3
X
Rn e−β1 R cn (θj , rk )
n=0
+
3
X
m=0
12
Rm e−β2 R cm (θj , rk ),
(3.12)
3.3. THE FIT PROCEDURE
to perform a WLLS fit in R. The long range has another grid spacing for θ and r
than the short and intermediate range, so there are at certain θj and rk , grid points
missing in R which we need for this fit. Therefore, we use one-dimensional (1D)
interpolation in R to calculate additional points in order to create the same grid
spacing for both the short, intermediate and long range. With β1 = 0.805 and
β2 = 1.57 from the WLLS fit we find for all θj , rk that the largest absolute error
in the attractive region of the potential (V < 0) is 2.22 µEh . The largest relative
error for the long range (R > 12.5) and the repulsive part (V > V0 ) is 2.12%. All
errors in this section hold for the grid points that are used in the various WLLS fits.
In the second step, for each Rℓ of our grid at θj we use
V (Rℓ , θj , r) = c0 (Rℓ , θj ) +
7
X
3
r k e−αr ck (Rℓ , θj ).
(3.13)
k=0
We perform again a WLLS fit and find α = 0.00529. In this case, the grid points θj
and Rℓ rise to a relative maximum deviation of 2.01%. The largest absolute error
is 0.39 µEh for geometries where V < 0. In the third step, at each Rℓ , rk , we use
an expansion in Legendre polynomials
V (Rℓ , θ, rk ) =
12
X
Pl (cos θ)cl (Rℓ , rk ),
(3.14)
l=0
to carry out the final fit using WLLS. In the repulsive area near R = 3 and r = 4.5
the RCCSD(T) method did not convergence. Therefore we generated additional
points in this region by means of 1D extrapolation in r. These additional points
helped us to prevent the fit become unphysical (going to −∞ instead of ∞) in the
strongly repulsive area. This fit has a maximum absolute error of 6.20 µEh for
V < 0. For the repulsive part the maximum relative error is around 27% and for
the long range region the maximum relative error is 2.56%.
With this three-step procedure, we have been able to determine nonlinear fit
parameters β1 , β2 , and α. To evaluate the potential at an arbitrary R, θ, r, we
keep the nonlinear parameters fixed at their determined values. We follow the
three-step procedure again using Eq. (3.12) and the previously determined values
β1 , β2 , cn (θj , rk ), cm (θj , rk ) to find V (R, θj , rk ). In step two, we use Eq. (3.13)
to determine c0 (R, θj ) and ck (R, θj ) by performing the linear fit for given R and
all θj . Consequently, we evaluate V (R, θj , r). In the last step, using Eq. (3.14)
we perform the corresponding linear fit to determine cl (R, r). This allows us to
evaluate V (R, θ, r). The total fit for the Vsum surface using this method has an
average absolute error of 0.11 µEh in the attractive region. The average relative
error for all grid points with R > 12.5 (long range region) is 0.39% and for V > V0
(the repulsive part) the average relative error is 2.01%.
The performed sequential fit has some disadvantages: first, some Van der Waals
minima are not described accurately because the absolute errors can go up to 6 µEh
in that region. Second, above interaction energies of 0.01 µEh the relative errors
13
CHAPTER 3. THE INTERACTION POTENTIAL
can become very large up to values of 27%. Third, evaluating the potential at a
certain point takes two WLLS fits which is quite time consuming. However, from
the sequential fit we gained quite some knowledge. So we now know the minimum
and maximum number of summation terms needed to fit the surface of the sum
potential. With this knowledge we return to Eqs. (3.4)–(3.6) that define a global
fit. Note that the factors that depends on θ and r in the global fit and the sequential
fit are very similar. Due to the experience obtained with the sequential fit we were
finally also able to perform the global fit.
From the global fit using WLLS we find that βa = 2.57, βb = 1.48, β = 1.24,
αa = 0.0285, αb = 0.0171, αI = 0.0160, and αII = 0.0464. The determined linear parameters are available in Appendix B which reports the FORTRAN code for
generating the interaction potential V (R, θ, r). The fit has the following properties:
(a) For all points with V > V0 (short range), the largest relative error is 6.67%
and the average relative error is 0.32%;
(b) For all points with V < 0, the largest absolute error is 2.37 µEh and the
average absolute error is 0.13 µEh ;
(c) For all points with R > 12.5 (long range), the largest relative error is 3.23%
and the average relative error is 0.49%.
(d) The root means square (RMS) of the absolute errors for all weighted grid points
is 0.00398 [the weight w = wsr × wlr is defined by Eqs. (3.10) and (3.11)].
3.3.2 The difference potential
We started the fitting procedure by using the Eqs. (3.3)–(3.6) defining the global fit.
For the weight function we used the same formulas and coefficients as described in
section 3.3.1. The long range part with the n = 6, 7, 8, and 9 terms [see Eq. (3.6)]
II
was easily fitted. In the next step, the determined coefficients {cII
nl , cnlk , n ≤ 9}
and nonlinear parameters αII and β were kept fixed. To determine the number
of summations terms needed for the global fit, we used the same sequential fitting
method as for the sum potential. In the following step, we determined the remaining nonlinear and linear parameters in a WLLS fit and we obtained an accurate fit
up to r = 3 a0 .
The determined nonlinear parameters are βa = 0.846, βb = 2.43, β = 1.24,
αa = 0.0130, αb = 0.0058, αI = 0.0506, and αII = 0.0526. Appendix B
reports the FORTRAN code for generating the interaction potential V (R, θ, r) and
contains also the determined linear parameters. The fit of the difference surface,
up to r = 3, where we consider only the range 0.75 ≤ r ≤ 3 has the following
deviation errors for the used grid points:
(a) For all points with V > V0 , the largest relative error is 1.10% and the average
relative error is 0.04%;
14
3.4. RESULTS
(b) For all points with V < 0, the largest absolute error is 0.01 µEh and the
average absolute error is 0.0001 µEh ;
(c) For all points with R > 12.5 (long range), the largest relative error is 3.32%
and the average relative error is 0.90%.
(d) The RMS of the absolute errors for all weighted grid points is 1.25 · 10−4 [the
weight w = wsr × wlr is defined by Eqs. (3.10) and (3.11)].
3.4 Results
Figure 3.2 shows four cuts of the 3D OH–He interaction potential of the A′ surface
with the OH distance r fixed at r = 0.75 a0 , r = re , r = 3 a0 , and r = 4.5 a0 ,
where re is the experimentally determined equilibrium distance of the X 2 Π state
[52]. The horizontal and vertical axes give the Cartesian coordinates of the helium
atom in the xy plane. The OH molecule lies on the horizontal axis with the oxygen
left of the origin and the hydrogen right from the origin. Each geometry of the
complex (defined by xHe , yHe , r in a0 ) leads to an interaction energy. The resulting
interaction energies are plotted as contours in reciprocal centimeters. The four cuts
of the A′′ surface are displayed in Fig. 3.3.
We find that the global minimum of our interaction potential in the range of
0.75 ≤ r ≤ 4.5 a0 is at R = 5.142 a0 , θ = 0◦ , r = 0.75 a0 . The interaction energy
in this minimum is −45.45 cm−1 (−207.1 µEh ) which is the global minimum of
both surfaces. Note that our potential does not include the contribution of the
OH monomer interaction potential. Therefore, the mentioned global minimum is
not the absolute ground state of the OH–He complex. The A′′ surface has also a
local minimum of −17.40 cm−1 (−79.28 µEh ) at R = 6.083 a0 and θ = 180◦
when r = 0.75 a0 . When r is enlarged we see that at re on the A′ surface a long
trench appears at 0 ≤ θ ≤ 92◦ which has two minima. The shallowest one with
an interaction energy of −30.02 cm−1 (−136.8 µEh ) which has a bent geometry
(R = 5.685 a0 and θ = 69.2◦ ) and the linear geometry (R = 6.544 a0 and θ = 0◦ )
which has a dept of −27.19 cm−1 (−123.9 µEh ). At re the A′′ surface still has
two minima that have moved away from the origin but stay at their linear locations.
The well at R = 6.553 a0 , θ = 0◦ has become less deep with an interaction
energy of −27.19 cm−1 (−123.9 µEh ) while the other well near the oxygen side
is getting larger and deeper. Its energy is −21.68 cm−1 (−98.82 µEh ). At r = 3
a0 the trench of the A′ surface has moved from the H side to the O side and is less
shallow. The T-shaped minimum at R = 5.621 a0 , θ = 89.0◦ has a dept of −32.26
cm−1 (−147.0 µEh ). Note that the local minimum of the linear geometry near the
hydrogen atom is starting to vanish. For the A′′ surface a new local minima has
appeared at R = 6.838 a0 , θ = 65.5◦ which has a dept of −15.23 cm−1 (−69.40
µEh ). The well near the oxygen atom at R = 8.466 a0 is still getting larger and
deeper and has now an interaction energy of −24.63 cm−1 (−112.2 µEh ) while
the other well near the hydrogen atom is also starting to vanish. Although our
15
CHAPTER 3. THE INTERACTION POTENTIAL
(a) r=0.75 a
0
−16
−10
00
−5
−4 −8
0
200
−16
2
8
−2
4
−20
−16
1
10 00
10 00
00
0
6
Y
He
0
(a )
8
.2 5
−0 −0.5 −1
−2
10
−36
−
− 0.2
−1 0.5 5
−
2
−8 −4
12
0
5
10
(b) r=r =1.8324 a
e
0
20
00
0
7
2
−2
4
−24 −26 −28
10
100 0
0
1000
0
−2
2
6
Y
He
0
(a )
8
5
.2
.5
−0
−0 −1 −2
−4 −8
−16
10
−0
.2
−
5
−1 0.5
−
−4 2
−
−16 8
−20
12
0
−10
−5
0
5
10
(c) r=3 a
0
−0
.5
−0
−1
−
−16 8
−2
0
−2
−4
Y
4
300
0
00
2
100
1000
1000
0
4
4
−2
0
(a )
−28
−1 −2 −4 −8
6
−1
He
.5
8
−0
10
.25
12
−10
−5
0
5
10
(d) r=4.5 a
0
12
−1
−2
0
6
−1
2
−10
−5
−24
1
1 00
10 000
00
0
Y
4
0
6
−1
−8 −4
6
−2 −4
(a )
8
He
.5
−0
.2
5
−0
.5
−1
−0
10
0
XHe (a0)
5
10
FIG. 3.2: A′ surface cut through of the RCCSD(T) interaction potential for four values of the OH
distance (r). The energies of the contour levels are given in cm−1 . The origin corresponds to the
center of mass of OH and the position of the He atom is given by (xHe , yHe ) = (R cos θ, R sin θ).
potential is less accurate from r > 3.0 a0 , we still decided to make a cut of the
PES at r = 4.5 a0 . These cuts are displayed in panel (d) of Figs. 3.2 and 3.3. We
16
3.4. RESULTS
(a) r=0.75 a
0
10
6
Y
0
−8
−16
−10
500
00
−5
0
−4 − 8
2
28
−16 −
4
−8
10
10 0
10 00
00
0
He
0
(a )
8
.25
−0 −0.5 1 2
− −
−0
.2
−
− 0.5 5
−4 −2 1
12
5
10
(b) r=r =1.8324 a
e
0
12
−0
.25
−
−1 0.5
30
−20
−26
00
0
−16
Y
−8
−4 −8
−22
−10
6
2
100
100
0
1000
0
−1
4
0
−8
−4 −2
6
−2
He
0
(a )
8
5
.2
5
−0 −0. −1
10
−5
0
5
10
(c) r=3 a
0
−0
.5
−0
−1
0
−1
6
−1
−8
00
300
−20
2
−10
−5
10
0
00
10
10
00
−8
4
−14
−8
4
−4 −2
0
(a )
Y
6
−1 −2 −4
8
He
.5
−0
10
.25
12
0
0
5
10
(d) r=4.5 a
0
5
−0
−0
.
−1
−8
−12
6
Y
100
1000
10000
−16
4
2
0
−10
−4
He
0
−4
−2
(a )
8
.5
−2
−1
−0
10
.25
12
−5
X
0
(a )
He
5
10
0
FIG. 3.3: A′′ surface cut through of the RCCSD(T) interaction potential for four values of the OH
distance (r). The energies of the contour levels are given in cm−1 . The origin corresponds to the
center of mass of OH and the position of the He atom is given by (xHe , yHe ) = (R cos θ, R sin θ).
see that our potential still behaves well at this distance and that the surface does
not contain any unphysical abnormalities in the regions of interest away from the
17
CHAPTER 3. THE INTERACTION POTENTIAL
TABLE 3.2: This table lists the global minimum, local minima and various saddle points of the
A′ and A′′ PES and the corresponding geometries (R, θ, r). For comparison the results of the 2D
potential of Lee et al. [30] are given in parentheses. The OH(X 2 Π) equilibrium distance is given by
re = 1.8324 a0 . Note that these minima are calculated in the interval 0.75 ≤ r ≤ 3.00 a0 , where
our fit of the PES is most accurate.
A′ ,
A′′
Global
Saddle A′ , local A′′
Local A′ , A′′
Local A′
Saddle A′ , local A′′
Local A′ , A′′
Local A′′
Local A′
Saddle A′ , local A′′
R(a0 )
5.142
6.083
6.553
(6.544)
5.685
(5.688)
6.091
(6.091)
8.466
6.838
5.621
6.033
θ(◦ )
0.0
180.0
0.0
(0.0)
69.2
(68.6)
180.0
(180.0)
0.0
65.5
89.0
180.0
r(a0 )
0.750
0.750
re
(re )
re
(re )
re
(re )
3.000
3.000
3.000
3.000
V (µEh )
-207.1
-79.28
-123.9
(-123.3)
-136.8
(-136.8)
-98.82
(-99.28)
-64.10
-69.40
-147.0
-112.2
highly repulsive region at short distance. On the A′ surface the minimum of the
T-shaped geometry moves about 20◦ to the right and becomes about 3.2 cm−1 (15
µEh ) less deep and at θ = 180◦ a new local minimum appears. For the A′′ surface
there is only one minimum left which also lies at θ = 180◦ . Note that for both
surfaces the minima at θ = 0◦ have vanished completely. For all eight cuts we
also show where the potential becomes unphysical in the repulsive area. As can be
seen, the first contour that leads to unphysical results can vary from 5 · 104 cm−1 to
1 · 104 cm−1 . This depends on the distance between helium and the nearest atom.
We found out that the unphysical region is given by R < 2.5 a0 or by Ri < 1.5 a0
(i = a, b).
In table 3.2 we report the numerical values for the global minimum, local minima and various saddle points in the range of 0.75 ≤ r ≤ 3.00 a0 . The results
reported by Lee et al. [30] are given in parentheses and are in good agreement
(maximum deviation is less than 1%) with our results. Note that they only calculated their PES at a fixed value of r = re , while we extended the PES to evaluate
the interaction energy for an arbitrary r.
The long range interaction is for large R proportional to R−6 . The proportionally constant is called C6 , which in general is angular dependent. However, the
main contribution in our case comes from the isotropic cII
6,0 coefficient of the sum
II
potential, so that we will identify C6 = −c6,0 [see Eq. (3.5)]. We expect that when
r increases, the dipole moment of the OH molecule will also increase which leads
to an increase of the induced dipole moment of the helium atom. These enlarged
dipoles will also increase the value of C6 . Another interaction effect which also
contributes to C6 , is called dispersion. This effect is a correlation effect between
18
3.4. RESULTS
9
C6,0 (Eha60)
8
7
6
5
4
1
1.5
2
2.5
3
r (a0)
3.5
4
4.5
FIG. 3.4: The long range parameter C6,0 as function of r. Note that C6,0 = −cII
6,0 [see Eq. (3.5)].
electrons which results in attractive forces. When monomers of a system increase
in polarizability, these dispersive forces become stronger [28]. Therefore, we expect that in first instance the dispersion interaction also enlarges the value of C6 for
increasing r. At very large values of r the OH molecule will dissociate into hydrogen and oxygen and the dipole moment of the OH molecule is lost. This means that
we expect that the C6 coefficient will decrease after a certain value of r is reached.
We thus conclude that the parameter C6 depends on r, and that for increasing r
it will first rise to a maximum and then will decrease again. Our plot of the long
range parameter C6 in Fig. 3.4 is consistent with these expected properties. The
long range parameter C6,0 (r) has its maximum value of 8.8 Eh a60 at r = 3.1 a0
which is about 1.8 times the equilibrium distance of the OH ground state.
In table 3.3 we list the isotropic long range coefficient C6,0 and the anisotropic
coefficients C6,2 , C7,1 , and C7,3 . The coefficients are obtained from the fit of the
sum potential for various values of the OH interatomic distance r.
TABLE 3.3: Dependence of the long range coefficients Cn,l (r) in Eh an
0 for the OH–He interaction
at certain OH interatomic distances. Coefficients are obtained from the global fit of the sum potential.
Note that Cn,l (r) = −cII
n,l (r) [see Eq. (3.5)].
n, l
6,0
6,2
7,1
7,3
r = 0.75 a0
4.86
0.398
37.2
0.277
r = re
6.83
1.42
64.6
6.98
r = 3 a0
8.80
2.60
92.7
18.8
19
4
Inelastic scattering of OH by He
The three-dimensional OH(X 2 Π)–He potential presented in chapter 3 has recently
been used to perform scattering calculations. The potential and the results for the
scattering calculations have been published, see Ref. [1]. In this chapter we briefly
explain how the scattering calculations have been performed and we also briefly
discuss the results of these calculations. We find that when the OH vibrational
motion is included in the scattering, our new potential improves the agreement
with recent high-accuracy collision experiments.
4.1 Theory
The scattering of 2 Π-state molecules with 1 S-state atoms is well described in the
literature [49]. Therefore, we only summarize briefly the relevant theory used for
the quantum scattering treatment of rotationally inelastic collisions of OH molecules with the He atom. More details of the theory can be found in Refs. [53, 54].
The Hamiltonian of the OH(X 2 Π)–He collision complex is defined by
Ĥ =
X
−~2 ∂ 2
L̂2
R
+
+
|Λ′ iVΛ′ ,Λ (R, θ)hΛ| + ĤOH ,
2µR ∂R2
2µR2
′
(4.1)
Λ ,Λ
where R and θ are the Jacobi coordinates (see Fig. 2.1), µ is the reduced mass of
the complex, L̂ is the orbital angular momentum operator corresponding to endover-end rotation of the OH–He complex, and ĤOH is the Hamiltonian of the OH
molecule in the X 2 Π state. The 2 Π electronic state of OH has two degenerate
components Λ = ±1 which are the projections of the orbital electronic angular
momentum on the internuclear axis. The operators |Λ′ iVΛ′ ,Λ (R, θ)hΛ| represents
the electronic interaction of the complex and couple the different electronic states
Λ and Λ′ . Furthermore, the Hamiltonian of OH consists of the following parts: the
rotational kinetic energy, spin-orbit coupling and Λ-doubling [55, 56]. For the QM
calculations, the OH rotational constant B = 18.5487 cm−1 , the spin-orbit coupling constant A = −139.21 cm−1 , and Λ-doubling parameters p = 0.235 cm−1
and q = −0.0391 cm−1 have been used [57]. For more information about diatomic
molecules in a 2 Π state, see the appropriate Hamiltonian and corresponding energy
levels in Ref. [58].
21
CHAPTER 4. INELASTIC SCATTERING OF OH BY HE
As already explained in section 1.1, the degeneracy within the Π state of the
OH radical is lifted when the helium atom approaches. This results in two PESs
which can be expanded in Racah normalized spherical harmonics Cl,m
V1,1 = V−1,−1 =
V1,−1 = V−1,1 =
X
VA′ + VA′′
=
vl,0 (R)Cl,0 (θ, 0),
2
l
X
VA′′ − VA′
=
vl,2 (R)Cl,2 (θ, 0),
2
(4.2)
l
where A′ and A′′ refer to the reflection symmetry of the electronic states which are
described in section 1.1. The surfaces V1,1 and V1,−1 are also known as the sum
Vsum and difference Vdiff potential energy surface, respectively (see section 3.2).
The Racah normalized spherical harmonics of Eq. (4.2) are related to the Legendre
functions used in Eq. (3.4) and Eq. (3.5) as follows
Cl,0 (θ, 0) = Pl (cos θ),
1
(l − 2)! 2 (2)
Cl,2 (θ, 0) =
Pl (cos θ).
(l + 2)!
(4.3)
4.2 QM calculations and results
For the scattering calculations, the 3D potential has been applied using the following procedure. First, the OH vibrational motion in the potential which is generated by the monomer potential energy VOH (r) [59] and the potential energy of
the complex V (R, θ, r), is determined. Namely, for fixed R and θ, we have a
one-dimensional problem that can be solved easily by standard numerical methods, like the discrete variable representation based on sinc-functions (sinc-DVR)
[60]. Second, the resulting ground state energy for each R and θ is subtracted from
the vibrational energy of the OH monomer in the absence of the He atom. The
resulting adiabatic 2D PES is checked to be very similar to the potential that is
obtained when first the lowest vibrational level of the monomer potential VOH (r)
is calculated and then an averaging of the interaction V (R, θ, r) over this ground
state is performed. Because of the well separated higher vibrational levels of OH,
these levels should only have a small influence, therefore the two methods are expected to give similar results. Finally, we note that in Ref. [1] a correction to the
interaction potential on the level of the AV5Z basis set was calculated. However
this correction was shown to only have a slight influence on the scattering results.
For calculating the OH monomer eigenfunctions it is convenient to use a parity
adapted Hund’s case (a) basis set which is labeled by |Ω, J, MJ , pi where J is the
total angular momentum of the OH molecule, Ω and MJ are the projections on
the molecular and space-fixed quantization axes, and p is the wave function parity
under inversion of the coordinates. Note that |Ω| is a nearly good quantum number of the exact eigenfunctions of OH [61]. The total angular momentum of the
22
4.2. QM CALCULATIONS AND RESULTS
OH–He complex is given by the operator F̂ = Ĵ + L̂, where Ĵ is the angular momentum operator of the OH monomer and L̂ is the angular momentum operator of
the end-over-end rotation of the complex. The eigenfunctions of F̂ are obtained by
coupling the monomer basis with the special harmonics |L, ML i = YL,ML (ϑ, ϕ),
where ϑ and ϕ are the space-fixed spherical coordinates of vector R. When we
write the scattering wave functions as products of radial and angular functions they
can be expressed as
,MF ,P
ΨFβ,L
=
1 X F ,MF ,P
F ,P
χ ′ ′
(R)ψβF′,M
(R̂, r̂),
,L′
R ′ ′ β ,L ←β,L
(4.4)
β ,L
where R̂ and r̂ are unit vectors in the direction of R and r respectively. The Greek
letter β is a shorthand notation for the monomer quantum numbers (Fi , J) and
index i is used to distinguish between the F1 and F2 spin-orbit manifolds of the
OH eigenstates [61]. Note that the during the collision process the total angular
momentum F , its space-fixed projection MF and the parity of the complex P =
p(−1)L are conserved. The cross sections, which are the experimental observables,
are expressed in terms of the scattering matrix and can be obtained via standard
asymptotic matching procedures [62]. The obtained scattering matrix or S-matrix
is then related to the scattering amplitudes, which in turn determine the differential
cross sections [63].
For achieving convergence in the calculations of the cross sections, a basis set
was used which includes all OH rotational states up to an angular momentum of
J = 21/2. Also all partial wave contributions up to a total angular momentum
of F = 241/2 have been taken into account. The Numerov method was used for
the propagation of wavefunction, starting at 4 a0 and continuing up to 35 a0 . The
evaluation of the cross sections was done on a energy grid which had an interval
spacing of 1 cm−1 . This is well below the experimental energy resolution in all
cases.
The inelastic cross sections of collisions of OH (X 2 Π3/2 , J = 3/2, f ) radicals
with He atoms were calculated with the developed adiabatic potential. The results
of these calculations, including the experimental data of Ref. [27] and the scattering
results with the potential of Lee et al. [30] are shown in Fig. 4.1. The theoretical
data are convoluted with the experimental resolution. This is done by taking a
Gaussian energy distribution with a standard deviation that is a function of the
energy. The standard deviation has a value of 24 cm−1 at low collision energies
which changes to 59 cm−1 at the highest collision energies. Note that the results
show the relative cross sections, rather than the absolute cross sections, because
these relative cross sections are experimentally measured. More details about this
theoretical study and the experimental studies can be found in Refs. [1, 26, 27, 34].
When we compare the results of the two theoretical potentials with the results of
the experiment, it is noticed that the overall agreement with experimental data has
improved significantly with our adiabatic potential.
23
relative cross section / %
relative cross section / %
(a)
(b)
50
20
0
1
0
(c)
(d)
5
0.5
0
100
200
300
Collision energy / cm-1
0
400 100
200
300
Collision energy / cm-1
400
FIG. 4.1: Relative state-to-state inelastic cross sections for collisions of OH (X 2 Π3/2 , J = 3/2, f )
with He as a function of the collision energy. On the vertical axes of the plots, 100% refers to the
total inelastic cross section. The experimental data from Ref. [27] are plotted with dots, while the
theoretically calculated cross sections of Lee et al. [30] are shown as solid curves. The dashed lines
are the cross sections calculated from our developed adiabatic potential. Panel (a) shows the cross
sections of the populated J = 3/2, F1 , e (black), J = 5/2, F1 , e (purple), and J = 5/2, F1 , f
(orange) states; (b) the J = 1/2, F2 , e (brown) and J = 3/2, F2 , f (pink) states; (c) the J =
1/2, F2 , f (red) and J = 3/2, F2 , e (blue) states; (d) the J = 7/2, F1 , e (green) and J = 7/2, F1 , f
(yellow) states.
5
Conclusions
We have developed three-dimensional potential energy surfaces of the X̃ 2 Π state
of the OH–He complex which are suitable to study the effect of vibrational excitations on the OH + He collision dynamics. The developed OH(X 2 Π)–He potential
is most accurate for OH interatomic distances of 0.75 ≤ r ≤ 3.00 a0 and becomes
unreliable at r > 4.5 a0 . This implies that the potential can be used for collisional
studies of OH(X 2 Π)–He up to ν = 3 (see Fig. 3.1). For studies with highly excited
(ν ≥ 11) OH radicals, an OH(X 2 Π)–He potential needs to be developed which describes the bond cleavage of the OH radical accurately. For this end, multireference
coupled cluster (MRCC) theory will be required. Several MRCC methods have already been developed and recently, the state-specific MRCC approach suggested
by Mukherjee and co-workers (Mk-MRCC) [64] emerged as the best candidate
theory. However, this method still does not provide the high accuracy of single
reference coupled cluster methods [65].
From our potential a 2D adiabatic vibrationally averaged ground state potential
has been developed which has recently been used for inelastic scattering calculations of OH (X 2 Π3/2 , J = 3/2, f ) with He. The results show that the potential
is realistic and that the calculated relative cross sections from our potential have
improved the agreement with experimental data [34] significantly. The potential
and the results of these scattering calculations have been published in Ref. [1].
25
Bibliography
[1] Koos B. Gubbels, Qianli Ma, Millard H. Alexander, Paul J. Dagdigian,
Dick Tanis, Gerrit C. Groenenboom, Ad van der Avoird, and Sebastiaan
Y. T. van de Meerakker. Resonances in rotationally inelastic scattering of
OH(X 2 Π) with helium and neon. J. Chem. Phys., 136:144308–144320, 2012.
doi: 10.1063/1.3697816.
[2] L. Zhang and A.J.C. Varandas. Dynamics of the OH(ν = 1, 2, 4) +
O3 atmospheric reaction. Phys. Chem. Chem. Phys., 3:1439–1445, 2001.
doi: 10.1039/B010149O.
[3] T.E. Graedel and P.J. Crutzen. Atmospheric Change: An Earth System Perspective. W. H. Freeman & Co, New York, NY, 1993. ISBN 0716723344.
[4] J. W. Meriwether. A review of the photochemistry of selected nightglow
emissions from the mesopause. J. Geophys. Res., 94:14629–14646, 1989.
doi: 10.1029/JD094iD12p14629.
[5] P. O. Wennberg, R. C. Cohen, R. M. Stimpfle, J. P. Koplow, J. G.
Anderson, R. J. Salawitch, D. W. Fahey, E. L. Woodbridge, E. R. Keim,
R. S. Gao, C. R. Webster, R. D. May, D. W. Toohey, L. M. Avallone,
M. H. Proffitt, M. Loewenstein, J. R. Podolske, K. R. Chan, and S. C.
Wofsy. Removal of stratospheric O3 by radicals: In situ measurements
of OH, HO2 , NO, NO2 , ClO, and BrO. Science, 266:398–404, 1994.
doi: 10.1126/science.266.5184.398.
[6] M. E. Summers, R. R. Conway, D. E. Siskind, M. H. Stevens, D. Offermann,
M. Riese, P. Preusse, D. F. Strobel, and J. M. Russell III. Implications of
satellite OH observations for middle atmospheric H2 O and ozone. Science,
277:1967–1970, 1997. doi: 10.1126/science.277.5334.1967.
[7] P. O. Wennberg, T. F. Hanisco, L. Jaeglé, D. J. Jacob, E. J. Hintsa, E. J.
Lanzendorf, J. G. Anderson, R.-S. Gao, E. R. Keim, S. G. Donnelly, L. A. Del
Negro, D. W. Fahey, S. A. McKeen, R. J. Salawitch, C. R. Webster, R. D. May,
R. L. Herman, M. H. Proffitt, J. J. Margitan, E. L. Atlas, S. M. Schauffler,
F. Flocke, C. T. McElroy, and T. P. Bui. Hydrogen radicals, nitrogen radicals,
and the production of O3 in the upper troposphere. Science, 279:49–53, 1998.
doi: 10.1126/science.279.5347.49.
[8] James A. Dodd, Steven J. Lipson, John R. Lowell, Peter S. Armstrong,
William A. M. Blumberg, Richard M. Nadile, Steven M. Adler-Golden,
William J. Marinelli, Karl W. Holtzclaw, and B. David Green. Analysis of hydroxyl earthlimb airglow emissions: Kinetic model for stateto-state dynamic of OH (ν, N ). J. Geophys. Res., 99:3559–3586, 1994.
doi: 10.1029/93JD03338.
27
BIBLIOGRAPHY
[9] A.B. Meinel. OH emission bands in the spectrum of the night sky. Astrophys.
J., 111:555–564, 1950. doi: 10.1086/145296.
[10] L. M. Kagan, M. J. Nicolls, M. C. Kelley, H. C. Carlson, V. V.
Belikovich, N. V. Bakhmet’eva, G. P. Komrakov, T. S. Trondsen, and
E. Donovan. Observation of radio-wave-induced red hydroxyl emission
at low altitude in the ionosphere. Phys. Rev. Lett., 94:095004, 2005.
doi: 10.1103/PhysRevLett.94.095004.
[11] Mark Wardle and Farhad Yusef-Zadeh.
Supernova remnant OH
masers: Signposts of cosmic collision. Science, 296:2350–2354, 2002.
doi: 10.1126/science.1068168.
[12] Joel M. Weisberg, Simon Johnston, Bärbel Koribalski, and Snezana
Stanimirović. Discovery of pulsed OH maser emission stimulated by a pulsar.
Science, 309:106–110, 2005. doi: 10.1126/science.1112494.
[13] Malcolm D. Gray, Kevin N. Jones, and David Field. OH masers and the
nature of massive star-forming regions. J. Chem. Soc., Faraday Trans., 89:
2231–2237, 1993. doi: 10.1039/FT9938902231.
[14] P. F. Bowers, K. J. Johnston, and J. H. Spencer. Microwave OH maser emission in the circumstellar envelopes of late-type stars. Nature, 291:382–385,
1981. doi: 10.1038/291382a0.
[15] Boncho P. Bonev, Michael J. Mumma, Neil Dello Russo, Erika L. Gibb,
Michael A. DiSanti, and Karen Magee-Sauer. Infrared OH prompt emission
as a proxy of water production in comets: Quantitative analysis of the multiplet near 3046 cm−1 in comets c/1999 h1 (Lee) and c/2001 a2 (LINEAR).
Astrophys. J., 615:1048–1053, 2004. doi: 10.1086/424587.
[16] Steven Adler-Golden. Kinetic parameters for OH nightglow modeling consistent with recent laboratory measurements. J. Geophys. Res., 102:19969–
19976, 1997. doi: 10.1029/97JA01622.
[17] Mark J. Dyer, Karen Knutsen, and Richard A. Copeland. Energy transfer in
the ground state of OH: Measurements of OH(ν=8,10,11) removal. J. Chem.
Phys., 107:7809–7815, 1997. doi: 10.1063/1.475094.
[18] Raji Viswanathan, Michelle Dolgos, and Robert J. Hinde. Quasiclassical trajectory study of the vibrational quenching of hydroxyl radicals
through collision with O atoms. J. Phys. Chem. A, 111:783–792, 2007.
doi: 10.1021/jp0667947.
[19] J. Kłos, F. Lique, and M. Alexander. Temperature dependence of rotational
excitation rate coefficients of OH(X 2 Π) in collision with He. Chem. Phys.
Lett., 445:12–16, 2007. doi: 10.1016/j.cplett.2007.07.035.
28
BIBLIOGRAPHY
[20] M. C. van Beek, J. J. ter Meulen, and M. H. Alexander. Rotationally inelastic
collisions of OH(X 2 Π) + Ar I. state-to-state cross sections. J. Chem. Phys.,
113:628, 2000. doi: 10.1063/1.481839.
[21] Joop J. Gilijamse, Steven Hoekstra, Sebastiaan Y. T. van de Meerakker,
Gerrit C. Groenenboom, and Gerard Meijer. Near-threshold inelastic collisions using moleculear beams with a tunable velocity. Science, 313:1617–
1620, 2006. doi: 10.1126/science.1131867.
[22] Manuel Lara, John L. Bohn, Daniel Potter, Pavel Soldan, and Jeremy M.
Hutson. Ultracold Rb-OH collisions and prospects for sympathetic cooling.
Phys. Rev. Lett., 97:183201, 2006. doi: 10.1103/PhysRevLett.97.183201.
[23] L. González-Sánchez, E. Bodo, and F. A. Gianturco. Quantum scattering of OH(X 2 Π) with He(1 S): Propensity features in rotational
relaxation at ultralow energies.
Phys. Rev. A, 73:022703, 2006.
doi: 10.1103/PhysRevA.73.022703.
[24] Zoran Pavlovic, Timur V. Tscherbul, Hossein R. Sadeghpour, Gerrit C.
Groenenboom, and Alex Dalgarno. Cold collisions of OH(2 Π) molecules
with He atoms in external fields. J. Phys. Chem. A, 113:14670–14680, 2009.
doi: 10.1021/jp904512r.
[25] Paul J. Dagdigian and Millard H. Alexander. Tensor cross sections and collisional depolarization of OH(X 2 Π) in collisions with helium. J. Chem. Phys.,
130:164315, 2009. doi: 10.1063/1.3119978.
[26] Ludwig Scharfenberg, Jacek Kłos, Paul J. Dagdigian, Millard H. Alexander,
Gerard Meijer, and Sebastiaan Y. T. van de Meerakker. State-to-state inelastic scattering of Stark-decelerated OH radicals with Ar atoms. Phys. Chem.
Chem. Phys., 12:10660–10670, 2010. doi: 10.1039/C004422A.
[27] Ludwig Scharfenberg, Koos B. Gubbels, Moritz Kirste, Gerrit C.
Groenenboom, A. van der Avoird, Gerard Meijer, and Sebastiaan Y. T. van de
Meerakker. Scattering of Stark-decelerated OH radicals with rare-gas atoms.
Eur. Phys. J. D, 65:189–198, 2011. doi: 10.1140/epjd/e2011-20009-4.
[28] Anthony J. Stone. The Theory of Intermolecular Forces. Oxford University
Press, Oxford, 1996.
[29] David R. Yarkony. A theoretical treatment of the predissociation of the individual rovibronic levels of OH/OD(A2 Σ+ ). J. Chem. Phys., 97:1838–1850,
1992. doi: 10.1063/1.463172.
[30] Hee-Seung Lee, Anne B. McCoy, Rafał R. Toczyłowski, and Sławomir M.
Cybulski. Theoretical studies of the X̃ 2 Π and Ã2 Σ+ states of the HeOH and Ne-OH complexes.
J. Chem. Phys., 113:5736–5749, 2000.
doi: 10.1063/1.1290605.
29
BIBLIOGRAPHY
[31] Grant Paterson, Sarantos Marinakis, Matthew L. Costen, Kenneth G.
McKendrick, Jacek Kłos, and Robert Toboła. Orientation and alignment depolarization in OH(X 2 Π)+Ar/He collisions. J. Chem. Phys., 129:074304,
2008. doi: 10.1063/1.2967861.
[32] Jiande Han and Michael C. Heaven. Bound states and scattering resonances
of OH(A)-He. J. Chem. Phys., 123:064307, 2005. doi: 10.1063/1.1993587.
[33] Alessandra Degli Esposti, Andreas Berning, and Hans-Joachim Werner.
Quantum scattering studies of the Λ doublet resolved rotational energy transfer of OH(X 2 Π) in collisions with He and Ar. J. Chem. Phys., 103:2067–
2082, 1995. doi: 10.1063/1.469682.
[34] Moritz Kirste, Ludwig Scharfenberg, François Lique Jacek Kłos, Millard H.
Alexander, Gerard Meijer, and Sebastiaan Y. T. van de Meerakker. Lowenergy inelastic collisions of OH radicals with He atoms and D2 molecules.
Phys. Rev. A, 82(4):042717, 2010. doi: 10.1103/PhysRevA.82.042717.
[35] P. J. Knowles, C. Hampel, and H.-J. Werner. Coupled cluster theory for high
spin, open shell reference wave function. J. Chem. Phys., 99:5219–5227,
1993. doi: 10.1063/1.465990.
[36] P. J. Knowles, C. Hampel, and H.-J. Werner.
Erratum: “Coupled
cluster theory for high spin, open shell reference wave function” [J.
Chem. Phys. 99, 5219 (1993)]. J. Chem. Phys., 112:E3106–3107, 2000.
doi: 10.1063/1.480886.
[37] Thom H. Dunning. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys., 90:
1007, 1989. doi: 10.1063/1.456153.
[38] R. A. Kendall, T. H. Dunning, and R. J. Harrison. Electron affinities of the
first-row atoms revisited. Systematic basis sets and wave functions. J. Chem.
Phys., 96:6796, 1992. doi: 10.1063/1.462569.
[39] David E. Woon and Thom H. Dunning. Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response
properties. J. Chem. Phys., 100:2975, 1994. doi: 10.1063/1.466439.
[40] G. C. Groenenboom and N. Balakrishnan. The He-CaH (2 Σ+ ) interaction. I.
Three-dimensional ab initio potential energy surface. J. Chem. Phys., 118:
7380–7385, 2003. doi: 10.1063/1.1562946.
[41] Sławomir M. Cybulski and Rafał R. Toczyłowski. Ground state potential
energy curves for He2 , Ne2 , Ar2 , He–Ne, He–Ar, and Ne–Ar: A coupledcluster study. J. Chem. Phys., 111:10520, 1999. doi: 10.1063/1.480430.
30
BIBLIOGRAPHY
[42] MOLPRO, a package of ab initio programs designed by H.-J. Werner and P.
J. Knowles.
[43] S. F. Boys and F. Bernardi. The calculation of small molecular interactions
by the differences of separate total energies. Some procedures with reduced
errors. Mol. Phys., 19:553–566, 1970. doi: 10.1080/00268977000101561.
[44] Dylan Jayatilaka and Timothy J. Lee. Open-shell coupled-cluster theory. J.
Chem. Phys., 98:9734, 1993. doi: 10.1063/1.464352.
[45] Timothy J. Lee and Peter R. Taylor. A diagnostic for determining the quality
of single-reference electron correlation methods. Int. J. Quantum Chem., 36:
199–207, 1989. doi: 10.1002/qua.560360824.
[46] Jonathan C. Rienstra-Kiracofe, Wesley D. Allen, and Henry F. Schaefer. The
C2 H5 + O2 reaction mechanism: High-Level ab initio characterizations. J.
Phys. Chem. A, 104:9823–9840, 2000. doi: 10.1021/jp001041k.
[47] Millard H. Alexander and Gregory C. Corey. Collision induced transitions between 2 Π and 2 Σ states of diatomic molecules: Quantum theory and collisional propensity rules. J. Chem. Phys., 84:100, 1986.
doi: 10.1063/1.450831.
[48] K. T. Tang and J. Peter Toennies. An improved simple-model for the van
der Waals potential based on universal damping functions for the dispersion
coefficients. J. Chem. Phys., 80:3726, 1984. doi: 10.1063/1.447150.
[49] Millard H. Alexander. Rotationally inelastic collisions between a diatom
molecule in a 2 Π electronic state and a structureless target. J. Chem. Phys.,
76:5974, 1982. doi: 10.1063/1.442951.
[50] Millard H. Alexander. Quantum treatment of rotationally inelastic collisions
involving molecules in Π electronic states: New derivation of the coupling
potential. Chem. Phys., 92:337, 1985. doi: 10.1016/0301-0104(85)85029-1.
[51] G. C. Groenenboom and I. M. Struniewicz. Three-dimensional ab initio potential energy surface for He-O2 . J. Chem. Phys., 113:9562–9566, 2000.
doi: 10.1063/1.1321311.
[52] K. P. Huber and G. Herzberg. Molecular Spectra and Molecular Structure.
IV. Constants of Diatomic Molecules. Van Nostrand Reinhold, New York,
1979.
[53] H. Klar.
Theory of collision induced rotational energy transfer in
the Π-state of diatomic molecules.
J. Phys. B, 6:2139–2149, 1973.
doi: 10.1088/0022-3700/6/10/025.
31
BIBLIOGRAPHY
[54] Timur V. Tscherbul, Gerrit C. Groenenboom, Roman V. Krems, and
Alexander Dalgarno. Dynamics of OH(2 Π)-He collisions in combined
electric and magnetic fields.
Faraday Discuss., 142:127–141, 2009.
doi: 10.1039/b819198k.
[55] J.M. Brown, M. Kaise, C.M.L. Kerr, and D.J. Milton. A determination of
fundamental Zeeman parameters for the OH radical. Molecular Physics, 36:
553–582, 1978. doi: 10.1080/00268977800101761.
[56] J. M. Brown and A. J. Merer. Lambda-type doubling parameters for molecules in Π electronic states of triplet and higher multiplicity. J. Mol. Spectrosc., 74:488–494, 1979. doi: 10.1016/0022-2852(79)90172-3.
[57] J.P. Maillard, J. Chauville, and A.W. Mantz. High-resolution emission spectrum of OH in an oxyacetylene flame from 3.7 to 0.9 µm. J. Mol. Spec., 63:
120–141, 1976. doi: 10.1016/0022-2852(67)90139-7.
[58] Peter F. Bernath. Spectra of atoms and molecules. Oxford University Press,
New York, 2005. ISBN 0-19-517759-2.
[59] Mark P. J. van der Loo and Gerrit C. Groenenboom. Theoretical transition
probabilities for the OH Meinel system. J. Chem. Phys., 126:114314, 2007.
doi: 10.1063/1.2646859.
[60] Gerrit C. Groenenboom and Daniel T. Colbert. Combining the discrete variable representation with the S-matrix Kohn method for quantum reactive scattering. J. Chem. Phys., 99:9681–9696, 1993. doi: 10.1063/1.465450.
[61] S.Y.T. van de Meerakker. Deceleration and Electrostatic Trapping of OH
radicals. PhD thesis, University of Nijmegen, 2006.
[62] Bruce R. Johnson.
The multichannel log-derivative method
for scattering calculations.
J. Comp. Phys., 13:445–449, 1973.
doi: 10.1016/0021-9991(73)90049-1.
[63] M. S. Child. Molecular Collision Theory. Academic, New York, 1974.
[64] U. S. Mahapatra, B. Datta, and D. Mukherjee. A state-specific multi-reference
coupled cluster formalism with molecular applications. Mol. Phys., 94(1):
157–171, 1998. doi: 10.1080/002689798168448.
[65] Monika Musiał, Ajith Perera, and Rodney J. Bartlett. Multireference
coupled-cluster theory: The easy way. J. Chem. Phys., 134:114108, 2011.
doi: 10.1063/1.3567115.
32
A
MOLPRO script RCCSD(T)
calculations
Listing A.1: Main MOLPRO script, CCSD_X2PI.sh
# ! / b i n / t c s h −f x
s e t o u t d i r = / s c r a t c h / $USER / s r a n g e / # F o r m o l p r o o u t p u t s
m k d i r −p $ o u t d i r
s e t e n v RUNNR $1
s e t e n v ANGLE ‘ nawk ’NR== n r { p r i n t } ’ n r =$RUNNR . / a n g l e s . d a t ‘
ech o t h e t a =$ANGLE
t i m e / v o l / t h ch em / m o l p r o 2 0 0 2 . 6 / b i n / m o l p r o << EOF >> CCSD_theta$ {ANGLE} _R12 .5−3 _ r 2 . 2 5 . o u t
∗∗∗ He −−> (H−−>O) CCSD ( T )
memory , 1 0 0 ,m
g t h resh , e n e r g y = 1 . 0 d −12;
gprint , basis , o r b i t a l
! Intermolecular coordinates :
theta
= $ANGLE
lR_s
= [ 1 2 . 5 d0 1 2 . 2 5 d0 1 1 . 7 5 d0 1 1 . 5 d0 1 1 . 2 5 d0 1 0 . 7 5 d0 1 0 . 5 d0 1 0 . 2 5 d0 10 d0 9 . 7 5 d0
9 . 5 d0 9 . 2 5 d 9 d0 8 . 7 5 d0 0 8 . 5 d0 8 . 2 5 d0 8 d0 7 . 7 5 d0 7 . 5 d0 7 . 2 5 d0 7 d0 6 . 7 5 d0
6 . 5 d0 6 . 2 5 d0 6 d0 5 . 7 5 d0 5 . 5 d0 5 . 2 5 d0 5 d0 4 . 7 5 d0 4 . 5 d0 4 . 2 5 d0 4 d0 3 . 7 5 d0
3 . 5 0 d0 3 . 2 5 d0 3 d0 ] AU
sr_s
= [ 2 . 2 5 d0 ] AU
! Properties
ispin
=
nelctot =
nelcHO =
n el cH e =
co m p l ex
1 ! spin doublet
11 ! t o t a l number o f e l e c t r o n s
9 !
2 !
R
= lR_s ( 1 )
rOH
= sr_s (1)
i n c l u d e co o r d _ Be . i n p
dummy , Be
! s e t t i n g e x p o n e n t s o f bond f u n c t i o n s
a_sp2
a_sp1
a_sp3
a_df1
a_df2
a_gf1
=
=
=
=
=
=
0 . 3 4 ∗ ( 5 . 3 d0 / R_atomOH ) ^ 2
2 .7 6 4 7 ∗ a _ s p 2
a_sp2 / 2 . 8 3 3
a _ s p 2 ∗ 1 .8 8 2
a_df1 / 2 .7 8 3
a_sp1 / 2 . 6 8 6
basis
H= a v t z ,
sp , Be ,
d f , Be ,
g , Be ,
O= a v t z , He= a v t z
a_ sp 1 , a_ sp 2 , a _ s p 3 ;
a_ d f 1 , a _ d f 2 ;
a_gf1 ;
end
RHF; wf , n e l c t o t , 1 , i s p i n ; ! RHF # 0 , i n p u t a t d e n s , o u t p u t 2 1 0 2 . 2
start , atdens
! Overwrite the d e f a u l t i n p u t 2100.2
sa v e , 2 1 0 2 . 2
! Overwrite the d e f a u l t output 2100.2
RHF; wf , n e l c t o t , 2 , i s p i n ; ! RHF # 1 , i n p u t a t d e n s , o u t p u t 2 1 0 3 . 2
start , atdens
! Overwrite the d e f a u l t i n p u t 2101.2
sa v e , 2 1 0 3 . 2
! Overwrite the d e f a u l t output 2101.2
33
APPENDIX A. MOLPRO SCRIPT RCCSD(T) CALCULATIONS
ij = 0
! Loop o v e r l a r g e R
do i =1 ,# l R _ s
do j =1 ,# s r _ s ! l o o p o v e r rOH , i f r e q u i r e d
ij
= ij + 1
R
= lR_s ( i )
rOH
= sr_s ( j )
i n c l u d e co o r d _ Be . i n p
a_sp2
a_sp1
a_sp3
a_df1
a_df2
a_gf1
=
=
=
=
=
=
0 . 3 4 ∗ ( 5 . 3 d0 / R_atomOH ) ^ 2
2 .7 6 4 7 ∗ a _ s p 2
a_sp2 / 2 . 8 3 3
a _ s p 2 ∗ 1 .8 8 2
a_df1 / 2 .7 8 3
a_sp1 / 2 . 6 8 6
basis
H= a v t z , O= a v t z , He= a v t z
sp , Be , a_ sp 1 , a_ sp 2 , a _ s p 3 ;
d f , Be , a_ d f 1 , a _ d f 2 ;
g , Be , a _ g f 1 ;
end
t e x t , He−HO l o w e s t A’ s t a t e
dummy , Be
RHF; wf , n e l c t o t , 1 , i s p i n ;
! RHF # 2 , i n p u t 2 1 0 2 . 2 o u t p u t 2 1 0 2 . 2
maxit , 9 0 ;
RCCSD( t ) ; maxit , 1 9 0
! RCCSD # 0 , i n p u t p r e v i o u s RHF
ccsdt1 ( i j ) = energt (1)
t e x t , He−HO l o w e s t A’ ’ s t a t e
dummy , Be
RHF; wf , n e l c t o t , 2 , i s p i n ;
! RHF # 3 , i n p u t 2 1 0 3 . 2 o u t p u t 2 1 0 3 . 2
maxit , 9 0 ;
RCCSD( t ) ; maxit , 1 9 0
! RCCSD # 1 , i n p u t p r e v i o u s RHF
ccsdt2 ( i j ) = energt (1)
t e x t , HO g r o u n d A’ s t a t e , He−g h o s t
dummy , He , Be ;
RHF; wf , nelcHO , 1 , i s p i n ;
! RHF # 4 , i n p u t 2 1 0 4 . 2 o u t p u t 2 1 0 4 . 2
maxit , 9 0 ;
RCCSD( t ) ; maxit , 1 9 0
! RCCSD # 2 , i n p u t p r e v i o u s RHF
HO_1 ( i j ) = e n e r g t ( 1 )
t e x t , HO g r o u n d A’ ’ s t a t e , He−g h o s t
dummy , He , Be ;
RHF; wf , nelcHO , 2 , i s p i n ;
! RHF # 5 , i n p u t 2 1 0 5 . 2 o u t p u t 2 1 0 5 . 2
maxit , 9 0 ;
start ,2105.2
sa v e , 2 1 0 5 . 2
RCCSD( t ) ; maxit , 1 9 0
! CCSD # 3 , i n p u t p r e v i o u s RHF
HO_2 ( i j ) = e n e r g t ( 1 )
t e x t , He g r o u n d s t a t e , HO−g h o s t
dummy , Be , H, O ;
RHF; wf , n el cH e , 1 , 0 ; maxit , 9 0 ;
CCSD; maxit , 1 9 0
He_1 ( i j ) = e n e r g c
! RHF # 4 , i n p u t 2 1 0 6 . 2 o u t p u t 2 1 0 6 . 2
! CCSD # 4 , i n p u t p r e v i o u s RHF
! He , no t r i p l e s so u s e CCSD , r e s u l t s t o r e d i n e n e r g c
lRR ( i j ) = R
s r r ( i j ) = rOH
th ( i j ) = theta
D1 ( i j ) = c c s d t 1 ( i j )−HO_1 ( i j )−He_1 ( i j )
D2 ( i j ) = c c s d t 2 ( i j )−HO_2 ( i j )−He_1 ( i j )
t a b l e , lRR , s r r , t h , c c s d t 1 , c c s d t 2 , D1 , D2
format , ’ ( f 1 0 . 2 , f 1 0 . 2 , 1 f 1 3 . 6 , 4 f 2 0 . 1 2 ) ’
enddo
enddo
t a b l e , lRR , s r r , t h , c c s d t 1 , c c s d t 2 , D1 , D2
format , ’ ( f 1 0 . 2 , f 1 0 . 2 , 1 f 1 3 . 6 , 4 f 2 0 . 1 2 ) ’
t i t l e , R[AU] , r [AU] , t h e t a ( H−−>O ) , RCCSD( T ) ( A’ ) , RCCSD( T ) ( A " ) , D(A ’ ) , D(A " )
−−−
EOF
Listing A.2: Bondorbital coordinates, coord_Be.inp
! G eo m et ry
OH
mH
= 1.0078250321
34
mO
mtot
! rHO
rHO
rHC
rOC
b2a
= 15.9949146221
= mH+mO
= 0 . 9 6 9 6 6 ! o p t i m a l b o n d l e n g t h f o r X2PI s t a t e i n a n g s t r o m
= rOH
= mO∗rHO / m t o t
= mH∗rHO / m t o t
= 0.529177249
R_H
R_O
= s q r t ( R^2+rHC^2−2∗R∗rHC∗ c o s ( t h e t a ) )
= s q r t ( R^2+rOC^2−2∗R∗rOC∗ c o s (180− t h e t a ) )
! Cartesian c o o r d i n a t e s ( angstrom ) :
Ox
Oy
Oz
= −rOC∗b2a
= 0 d0
= 0 d0
H1x
H1y
H1z
=
=
=
rHC∗b2a
0 d0
0 d0 ! word Hz i s r e s e r v e d by MOLPRO
IF ( R_H . LT . R_O ) THEN
shorterbond_H = 1
R_atomOH = R_H
Bex = ( ( R−(R_H / 2 d0 ) ) ∗ c o s (
Bey = 0 d0
Bez = ( ( R−(R_H / 2 d0 ) ) ∗ s i n (
ELSE
shorterbond_O = 1
R_atomOH = R_O
Bex = ( ( R−(R_O / 2 d0 ) ) ∗ c o s (
Bey = 0 d0
Bez = ( ( R−(R_O / 2 d0 ) ) ∗ s i n (
! Bez = R / 2 d0∗b2a
ENDIF
Hex
Hey
Hez
=
=
=
t h e t a ) ) ∗ b2a
t h e t a ) ) ∗ b2a
t h e t a ) ) ∗ b2a
t h e t a ) ) ∗ b2a
( R∗ c o s ( t h e t a ) ) ∗ b2a
0 d0
( R∗ s i n ( t h e t a ) ) ∗ b2a
geomtyp = xyz
g e o m e t r y ={
4
(H−−>O) −−> He , H−−>O i n o r i g i n
He , Hex , Hey , Hez
H , H1x , H1y , H1z
O , Ox , Oy , Oz
Be , Bex , Bey , Bez
}
Listing A.3: Angles data, angles.dat
0.0001
12
24
36
48
60
72
84
96
108
120
132
144
156
168
180.0001
Listing A.4: Frontend, start.sh
# ! / b i n / t c s h −e f
@ i = 1
w h i l e ( $ i <= 16 )
ech o $ i
. / CCSD_X2PI . sh $ i &
@ i = $i + 1
end
35
B
FORTRAN code OH–He potential
Listing B.1: Frontend, main.f
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
c T h i s program c a l c u l a t e s t h e i n t e r a c t i o n e n e r g i e s o f OH−He PES
c From s t a n d a r d i n p u t i t r e a d s t h e OH−He c o o r d i n a t e s and p r i n t s
c o u t V_A ’ , V_A ’ ’ , ( V_A’+V_A ’ ’ ) / 2 and ( V_A’’−V_A ’ ) / 2 i n cm−1
c
c
mO= 1 5 .9 9 4 9 1 4 6 2 2 1 , mH= 1 .0 0 7 8 2 5 0 3 2 ( NIST . gov )
c
R : d i s t a n c e f ro m He t o c e n t e r o f mass o f OH−H ( b o h r ) .
c
z : e q u a l s c o s ( t h e t a ) , t h e t a =0 c o r r e s p o n d s t o t h e l i n e a r
c
g e o m e t r y OH−He
c
r1 : O−H bond d i s t a n c e ( b o h r ) .
c
c
The i n p u t f o r m a t i s :
c
c l i n e 1: npoint
( t h e number o f p o i n t s , max 1 0 0 0 0 0 )
c l i n e 2 : R_1 , z_1 , r1 _ 1
( coordinates of the f i r s t points )
c l i n e 3 : R_2 , z_2 , r2 _ 1
( i d em f o r t h e s e c o n d p o i n t )
c ....
c l i n e n p o i n t +1: R _ n p o i n t , z _ n p o i n t r 1 _ n p o i n t ( i d em f o r t h e l a s t p o i n t )
c
c V e r s i o n 1 . 1 1 (25−Mar−2012)
c A u t h o r s : D . T a n i s & G . Groenenboom
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
program OHHePES
i m p l i c i t none
c
i n t e g e r i , j , i p o t , mp , n p o i n t
parameter ( mp=1 0 0 0 0 0 )
c
r e a l ∗8 Rs ( mp ) , z s ( mp ) , r 1 s ( mp ) , p o t s ( mp , 4 ) , cm1
r e a l ∗8 R t s t , z t s t , r 1 t s t , v0sum , v 0 d i f
c
cm1 = 1 d0 / 2 1 9 4 7 4 . 6 3 1 3 7 0 9 8 d0
c
c
Check i n t e g r i t y o f t h e OH−He p o t e n t i a l
R t s t = 6 . 7 d0
z t s t = 0 . 5 d0
r 1 t s t = 2 d0
c
v0sum = −0.914094522474498d−04 ! 26 March 2 0 1 1 , g77 ( gcc − 3 . 4 . 6 ) mO a s REAL∗4
c
v 0 d i f = −0.682289339301687d−04 ! 26 March 2 0 1 1 , g77 ( gcc − 3 . 4 . 6 ) mO a s REAL∗4
v0sum = −0.914094522474498d−04 ! 25 March 2 0 1 2 , g f o r t r a n ( gcc −4.4.3−4 u b u n t u 5 . 1 )
v 0 d i f = −0.682289340828068d−04 ! 25 March 2 0 1 2 , g f o r t r a n ( gcc −4.4.3−4 u b u n t u 5 , 1 )
c a l l ohhe_pot (1 , Rtst , z t s t , r 1 t s t , 1 , pots ( 1 , 1 ) )
i f ( d a b s ( p o t s ( 1 , 1 ) / v0sum−1d0 ) . g t . 1 d −10) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n vOH−He : i n t e g r i t y c h e c k n o t p a s s e d ’
w r i t e ( ∗ , ∗ ) ’∗∗ v_Ap ( ’ , R t s t , ’ , ’ , z t s t , ’ , ’ , r 1 t s t , ’ ) = ’ , p o t s ( 1 , 1 )
w r i t e ( ∗ , ∗ ) ’∗∗ t h i s s h o u l d be ’ , v0sum
s t o p 230
endif
c a l l ohhe_pot (1 , Rtst , z t s t , r 1 t s t , 2 , pots ( 1 , 2 ) )
if
( d a b s ( p o t s ( 1 , 2 ) / v 0 d i f −1d0 ) . g t . 1 d −10) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n vOH−He : i n t e g r i t y c h e c k n o t p a s s e d ’
w r i t e ( ∗ , ∗ ) ’∗∗ v_App ( ’ , R t s t , ’ , ’ , z t s t , ’ , ’ , r 1 t s t , ’ ) = ’ , p o t s ( 1 , 2 )
w r i t e ( ∗ , ∗ ) ’∗∗ t h i s s h o u l d be ’ , v 0 d i f
s t o p 230
endif
c
c
I n t e g r i t y p a s s e d . L et ’ s do some c a l c u l a t i o n s
37
APPENDIX B. FORTRAN CODE OH–HE POTENTIAL
c
rea d ( ∗ , ∗ ) n p o i n t
do i =1 , n p o i n t
rea d ( ∗ , ∗ ) Rs ( i ) , z s ( i ) , r 1 s ( i )
enddo
do i p o t = 1 , 4
c a l l o h h e _ p o t ( n p o i n t , Rs , zs , r 1 s , i p o t , p o t s ( 1 , i p o t ) )
end do
w r i t e ( ∗ , ∗ ) ’R i n b o h r , z= c o s ( t h e t a ) , OH−He p o t e n t i a l s i n cm−1’
w r i t e ( ∗ , ’ ( a8 , a8 , a8 , 4 a18 ) ’ ) ’R , ’ , ’ z , ’ , ’ r 1 ’ ,
+ ’V_Ap ’ , ’V_App , ’ , ’ ( V_Ap+V_App ) / 2 , ’ , ’ ( V_App−V_Ap ) / 2 ’
do i =1 , n p o i n t
w r i t e ( ∗ , ’ ( f 8 . 4 , f 8 . 4 , f 8 . 4 , 4 e18 . 9 ) ’ ) Rs ( i ) , z s ( i ) , r 1 s ( i ) ,
+ ( p o t s ( i , j ) / cm1 , j = 1 , 4 )
enddo
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
Listing B.2: Main program, oh-he.f
c
D i ck T a n i s , v e r s i o n 1 . 1 1 (25−Mar−2012)
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
s u b r o u t i n e p o t e n t i a l 3 d ( n , Rs , zs , r1 s , i p o t , p o t s )
c
c
function :
c
t o e v a l u a t e OH−He p o t e n t i a l s f o r an a r r a y o f n p o i n t s
c
c
input :
c
n : number o f p o i n t s
c
Rs ( n ) : d i s t a n c e atom − c e n t e r o f mass d i a t o m ( a_0 )
c
z s ( n ) : i n p u t : c o s ( t h e t a ) , t h e t a =J a c o b i a n g l e ( t h e t a =0 OH−−He )
c
i p o t : p o t e n t i a l s e l e c t i o n number
c
i p o t = 1 : V_A ’
c
i p o t = 2 : V_A ’ ’
c
i p o t = 3 : ( V_A’+V_A ’ ’ ) / 2
c
i p o t = 4 : ( V_A’’−V_A ’ ) / 2
c
c
output :
c
p o t s : OH( X^2 Pi)−He p o t e n t i a l e n e r g y ( H a r t r e e ) o f p o t e n t i a l n
c
c
error :
c
s t o p 1 0 2 : no v a l i d p o t e n t i a l s e l e c t e d
c
s u b r o u t i n e o h h e _ p o t ( n , Rs , zs , r 1 s , i p o t , p o t s )
i m p l i c i t none
integer n , ipot , i
r e a l ∗8 Rs ( n ) , z s ( n ) , r 1 s ( n ) , p o t s ( n ) , vs , vd
c
i f ( i p o t . eq . 3 ) t h en
do i =1 , n
c a l l ohhe_vsum ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , p o t s ( i ) )
enddo
e l s e i f ( i p o t . eq . 4 ) t h en
do i =1 , n
c a l l o h h e _ v d i f ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , p o t s ( i ) )
enddo
e l s e i f ( i p o t . eq . 1 ) t h en
do i =1 , n
c a l l ohhe_vsum ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , v s )
c a l l o h h e _ v d i f ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , vd )
p o t s ( i ) = vs−vd
enddo
e l s e i f ( i p o t . eq . 2 ) t h en
do i =1 , n
c a l l ohhe_vsum ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , v s )
c a l l o h h e _ v d i f ( Rs ( i ) , z s ( i ) , r 1 s ( i ) , vd )
p o t s ( i ) = v s +vd
enddo
else
s t o p 102
endif
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c s u b r o u t i n e ohhe_vsum ( R , z , r1 , v )
c
c function :
c
e v a l u a t e V_s a t c o o r d i n a t e R , z , r1
c
c input :
38
c
R : d i s t a n c e f ro m He t o c e n t e r o f mass o f OH ( a_0 ) .
c
z : e q u a l s c o s ( t h e t a ) , t h e t a =0 c o r r e s p o n d s t o t h e l i n e a r
c
g e o m e t r y OH−He
c
r1 : OH i n t e r n a l bond d i s t a n c e ( a_0 ) .
c
c output :
c
v : e n e r g y ( H a r t r e e ) o f t h e d i a b a t i c p o t e n t i a l V_s
c
at supplied coordinates
c
c wa rn i n g :
c
p o t e n t i a l i s l e s s a c c u r a t e f o r 3 . 0 0 < r1 <= 4 . 5
c
c error :
c
s t o p 2 3 1 : c o s t h e t a i s n o t b e t w e e n −1 o r 1 )
c
s t o p 2 3 2 : R < 2 . 5 o r Ra / Rb < 1 . 5 a_0 , i n t h e s e r e p u l s i v e r e g i o n s ,
c
p o t e n t i a l can become u n p h y s i c a l l a r g e n e g a t i v e
c
s t o p 2 3 3 : f o r r1 < 0 . 7 5 o r r1 > 4 . 5 a_0 , p o t e n t i a l i s u n r e l i a b l e
c
s u b r o u t i n e ohhe_vsum ( R , z , r 1 , v )
i m p l i c i t none
r e a l ∗8 R , z , r 1 , v
c
l o g i c a l d eb u g
r e a l ∗8 b e t a _ a , b e t a _ b , b e t a , a l p h a _ a , a l p h a _ b , a l p h a _ I , a l p h a _ I I
r e a l ∗8 amH , amO
r e a l ∗8 Ra , za , Rb , zb , p l a ( 8 ) , p l b ( 6 ) , p l ( 1 0 ) ,
+ Ran ( 4 ) , Rbn ( 4 ) , r 1 n _ I ( 4 ) , r 1 n _ I I ( 4 ) , r 1 a n ( 9 ) , r 1 b n ( 9 ) , Rn ( 8 ) ,
+ f _ l o n g , sum
i n t e g e r i , j , k , i n d , kmax , j 0
c
i n c l u d e ’ l c_ v su m . f ’ ! l i n e a r c o e f f i c i e n t s i n a r r a y c _ v s
d a t a b e t a _ a / 2 . 5 6 6 3 3 1 6 9 4 8 5 7 4 7 d0 / , b e t a _ b / 1 . 4 7 6 4 0 9 2 3 6 8 5 2 9 0 d0 /
d a t a b e t a / 1 . 2 4 d0 /
d a t a a l p h a _ a / 0 . 0 2 8 5 1 5 9 1 1 3 8 0 0 9 d0 / , a l p h a _ b / 0 . 0 1 7 1 4 2 9 6 5 2 8 8 4 1 d0 /
d a t a a l p h a _ I / 0 . 0 1 6 0 3 1 0 7 3 9 8 8 8 9 d0 / , a l p h a _ I I / 0 . 0 4 6 4 d0 /
d a t a amH / 1 . 0 0 7 8 2 5 0 3 2 d0 / , amO / 1 5 . 9 9 4 9 1 4 6 2 2 1 d0 / ! ( NIST . gov 2 0 0 5 )
c
d a t a d eb u g / . f a l s e . /
c
v =0 d0
c
i f ( r 1 . l t . 0 . 7 5 d0 . o r . r 1 . g t . 4 . 5 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n ohhe_vsum : p o t e n t i a l i s u n r e l i a b l e f o r
+ r 1 < 0 . 7 5 a_0 o r r 1 > 4 . 5 a_0 . ’
w r i t e ( ∗ , ∗ ) ’∗∗ R = ’ , R , ’ z = ’ , z , ’ r 1 = ’ , r 1
s t o p 233
endif
c
i f ( r 1 . g t . 3 . 0 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ WARNING f r o m ohhe_vsum : p o t e n t i a l i s l e s s a c c u r a t
+e f o r 3 < r 1 <= 4 . 5 a_0 . ’
w r i t e ( ∗ , ∗ ) ’∗∗ R = ’ , R , ’ z = ’ , z , ’ r 1 = ’ , r 1
endif
c
i f ( d a b s ( z ) . g t . 1 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n ohhe_vsum : i l l e g a l c o s ( t h e t a ) = z = ’ , z
s t o p 231
endif
c
c a l l j a c 2 y u 2 ( amH , amO , R , z , r 1 , Ra , za , Rb , zb )
i f ( d eb u g ) w r i t e ( ∗ , ’ (A , f 1 2 . 9 , A , f 1 2 . 9 , A , f 1 2 . 9 , A , f 1 2 . 9 ) ’ )
+ ’ Ra= ’ , Ra , ’ Rb= ’ , Rb , ’ za = ’ , za , ’ zb = ’ , zb
i f (R . l t . 2 . 5 d0 . o r . Ra . l t . 1 . 5 d0 . o r . Rb . l t . 1 . 5 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n ohhe_vsum : p o t e n t i a l c o u l d p r o d u c e ’
w r i t e ( ∗ , ∗ ) ’∗∗ n e g a t i v e v a l u e i n r e p u l s i v e r e g i o n f o r g e o m e t r y : ’
w r i t e ( ∗ , ∗ ) ’∗∗ R = ’ , R , ’ z = ’ , z , ’ r 1 = ’ , r 1
s t o p 232
endif
c
c a l l p _ l e g ( p l a , 7 , za ) ! p l a ( 1 ) = P_0 ( za ) . . p l a ( 8 ) = P_7 ( za )
c a l l p _ l e g ( p l b , 5 , zb ) ! p l b ( 1 ) = P_0 ( zb ) . . p l b ( 6 ) = P_5 ( zb )
Ran ( 1 ) = d ex p ((− b e t a _ a )∗ Ra )
Rbn ( 1 ) = d ex p ((− b e t a _ b )∗Rb )
c
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ Ran ( 1 ) = ’ , Ran ( 1 ) , ’ Rbn ( 1 ) = ’ , Rbn ( 1 )
r 1 n _ I ( 1 ) = d ex p ((− a l p h a _ I )∗ r 1 ∗∗3)
r 1 n _ I I ( 1 ) = d ex p ((− a l p h a _ I I )∗ r 1 ∗∗3)
r 1 a n ( 1 ) = d ex p ((− a l p h a _ a )∗ r 1 ∗∗3)
r 1 b n ( 1 ) = d ex p ((− a l p h a _ b )∗ r 1 ∗∗3)
do i =2 ,4
Ran ( i ) = Ra∗Ran ( i −1)
Rbn ( i ) =Rb∗Rbn ( i −1)
r 1 n _ I ( i ) = r 1 ∗ r 1 n _ I ( i −1)
r 1 n _ I I ( i ) = r 1 ∗ r 1 n _ I I ( i −1)
end do
39
APPENDIX B. FORTRAN CODE OH–HE POTENTIAL
do i =2 ,9
r 1 a n ( i ) = r 1 ∗ r 1 a n ( i −1)
r 1 b n ( i ) = r 1 ∗r 1 b n ( i −1)
end do
c
c
c
c
c
c
c
c
c
c
t h e r1 i n d e p e n d e n t p a r t o f t h e s h o r t r a n g e
i n d =0
sum=0 d0
v=0 d0
do i =1 ,2 ! l = 0 . . 1
i n d = i n d +1
sum=sum+ p l a ( i )∗ c _ v s ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ p l a ( i ) = ’ , p l a ( 1 )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ p l a ( i )∗ Ran ( 1 ) = ’ , p l a ( i )∗ Ran ( 1 )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ sum = ’ , sum
end do
v=Ran ( 1 ) ∗ sum
sum=0 d0
do i =1 ,2 ! l =0 t o 1
i n d = i n d +1
sum=sum+ p l b ( i )∗ c _ v s ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ sum = ’ , sum
end do
v=v+Rbn ( 1 ) ∗ sum
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v= ’ , v
t h e r1 d e p e n d e n t p a r t o f t h e s h o r t r a n g e
do i =1 ,4
! n =0..3
do j =1 ,8
! l =0..7
sum=0 d0
do k =1 ,9 ! k = 0 . . 8
i n d = i n d +1
sum=sum+ c _ v s ( i n d )∗ r 1 a n ( k )
end do
v=v+Ran ( i )∗ p l a ( j )∗ sum
end do
end do
do i =1 ,4
! n =0..3
do j =1 ,6
! l =0..5
sum=0 d0
do k =1 ,9 ! k = 0 . . 8
i n d = i n d +1
sum=sum+ c _ v s ( i n d )∗ r 1 b n ( k )
end do
v=v+Rbn ( i )∗ p l b ( j )∗ sum
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
c
c
c
the long range
c a l l p _ l e g ( p l , 9 , z ) ! p l ( 1 ) = P_0 ( z ) . . p l ( 1 0 ) = P_9 ( z )
do i =1 ,8
Rn ( i ) = f _ l o n g ( R , b e t a , 5 + i ) ! Rn ( 1 ) = f n ( B e t a∗R ) / R ^ 6 . . Rn ( 8 ) = f n ( B e t a∗R ) / R^13
end do
c
c
c
t h e l o n g ra n g e , t e r m s n = 1 0 . . 1 3
kmax=4
! k = 0..3
j 0 =1
do i =1 ,4
! n = i +9 ( n = 1 0 . . 1 3 −> i = 1 . . 4 )
j 0=1− j 0
! l m i n =0 ( ev en n ) , 1 ( odd n )
do j =1+ j 0 , i +6 ,2
! lmax = n−4 = i +5
i n d = i n d +1
sum= c _ v s ( i n d )
do k =1 , kmax
i n d = i n d +1
sum=sum+ r 1 n _ I ( k )∗ c _ v s ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i = ’ , i , ’ j = ’ , j , ’ k = ’ , k
end do
v=v+Rn ( i +4)∗ p l ( j )∗ sum ! s t a r t a t Rn ( 5 ) = f n ( B e t a∗R ) / R^10
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
c
c
c
t h e l o n g ra n g e , t e r m s n = 6 . . 9
j 0 =1
do i =1 ,4
j 0=1− j 0
do j =1+ j 0 , i +2 ,2
i n d = i n d +1
40
! n = i +5 ( n = 6 . . 9 −> i = 1 . . 4 )
! l m i n =0 ( ev en n ) , 1 ( odd n )
! lmax = n−4 = i +1
sum= c _ v s ( i n d )
do k =1 , kmax
i n d = i n d +1
sum=sum+ r 1 n _ I I ( k )∗ c _ v s ( i n d )
end do
v=v+Rn ( i )∗ p l ( j )∗ sum ! s t a r t a t Rn ( 1 ) = f n ( B e t a∗R ) / R^6
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c s u b r o u t i n e o h h e _ v d i f ( R , z , r1 , v )
c
c function :
c
e v a l u a t e V_d a t c o o r d i n a t e R , z , r1
c
c input :
c
R : d i s t a n c e f ro m He t o c e n t e r o f mass o f OH ( a_0 ) .
c
z : e q u a l s c o s ( t h e t a ) , t h e t a =0 c o r r e s p o n d s t o t h e l i n e a r
c
g e o m e t r y OH−He
c
r1 : OH i n t e r n a l bond d i s t a n c e ( a_0 ) .
c
c output :
c
v : e n e r g y ( H a r t r e e ) o f t h e d i a b a t i c p o t e n t i a l V_d
c
at supplied coordinates
c
c wa rn i n g :
c
p o t e n t i a l i s l e s s a c c u r a t e f o r 3 . 0 0 < r1 <= 4 . 5 a_0
c
c error :
c
s t o p 2 4 1 : c o s t h e t a i s n o t b e t w e e n −1 o r 1 )
c
s t o p 2 4 2 : R < 2 . 5 a_0 o r Ra / Rb < 1 . 5 a_0 , i n t h e s e r e p u l s i v e
c
r e g i o n s , p o t e n t i a l can become u n p h y s i c a l l a r g e n e g a t i v e
c
s t o p 2 4 3 : f o r r1 < 0 . 7 5 a_0 o r r1 > 4 . 5 a_0 , p o t e n t i a l i s u n r e l i a b l e
c
subroutine o h h e_ v d i f (R, z , r1 , v )
i m p l i c i t none
r e a l ∗8 R , z , r 1 , v
c
l o g i c a l d eb u g
r e a l ∗8 b e t a _ a , b e t a _ b , b e t a , a l p h a _ a , a l p h a _ b , a l p h a _ I , a l p h a _ I I
r e a l ∗8 amH , amO
r e a l ∗8 Ra , za , Rb , zb , p l a ( 5 ) , p l b ( 5 ) , p l ( 8 ) ,
+ Ran ( 6 ) , Rbn ( 6 ) , r 1 n _ I ( 4 ) , r 1 n _ I I ( 4 ) , r 1 a n ( 6 ) , r 1 b n ( 5 ) , Rn ( 8 ) ,
+ f _ l o n g , sum
i n t e g e r i , j , k , i n d , kmax , j 0
c
i n c l u d e ’ l c _ v d i f . f ’ ! l i n e a r c o e f f i c i e n t s i n a r r a y c_ v d
data
data
data
data
data
b e t a _ a / 0 . 8 4 6 3 9 6 3 2 5 0 2 4 5 6 3 d0 / , b e t a _ b / 2 . 4 2 7 3 5 4 1 0 3 1 1 3 3 6 1 d0 /
b e t a / 1 . 2 4 d0 /
a l p h a _ a / 0 . 0 1 2 9 5 7 0 9 1 4 4 1 1 6 8 d0 / , a l p h a _ b / 0 . 0 0 5 7 9 8 8 3 0 0 0 8 4 6 3 d0 /
a l p h a _ I / − 0 .0 5 0 6 0 8 8 6 9 800 03 7d0 / , a l p h a _ I I / 0 . 0 5 2 6 3 6 7 1 8 7 5 0 0 0 0 d0 /
amH / 1 . 0 0 7 8 2 5 0 3 2 d0 / , amO / 1 5 . 9 9 4 9 1 4 6 2 2 1 d0 / ! ( NIST . gov 2 0 0 5 )
c
d a t a d eb u g / . f a l s e . /
c
v =0 d0
c
i f ( r 1 . l t . 0 . 7 5 d0
w r i t e ( ∗ , ∗ ) ’∗∗
+ r 1 < 0 . 7 5 a_0 o r
w r i t e ( ∗ , ∗ ) ’∗∗
s t o p 243
endif
. o r . r 1 . g t . 4 . 5 d0 ) t h en
ERROR i n o h h e _ v d i f : p o t e n t i a l i s u n r e l i a b l e f o r
r 1 > 4 . 5 a_0 . ’
R = ’ , R, ’ z = ’ , z , ’ r1 = ’ , r1
c
i f ( r 1 . g t . 3 . 0 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ WARNING f r o m o h h e _ v d i f : p o t e n t i a l i s l e s s a c c u r a t
+e f o r 3 < r 1 <= 4 . 5 a_0 . ’
w r i t e ( ∗ , ∗ ) ’∗∗ R = ’ , R , ’ z = ’ , z , ’ r 1 = ’ , r 1
endif
c
i f ( d a b s ( z ) . g t . 1 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n o h h e _ v d i f : i l l e g a l c o s ( t h e t a ) = z = ’ , z
s t o p 241
endif
c a l l j a c 2 y u 2 ( amH , amO , R , z , r 1 , Ra , za , Rb , zb )
i f ( d eb u g ) w r i t e ( ∗ , ’ (A , f 1 2 . 9 , A , f 1 2 . 9 , A , f 1 2 . 9 , A , f 1 2 . 9 ) ’ )
+ ’ Ra= ’ , Ra , ’ Rb= ’ , Rb , ’ za = ’ , za , ’ zb = ’ , zb
i f (R . l t . 2 . 5 d0 . o r . Ra . l t . 1 . 5 d0 . o r . Rb . l t . 1 . 5 d0 ) t h en
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n o h h e _ v d i f : p o t e n t i a l c o u l d p r o d u c e ’
w r i t e ( ∗ , ∗ ) ’∗∗ n e g a t i v e v a l u e i n r e p u l s i v e r e g i o n f o r g e o m e t r y : ’
41
APPENDIX B. FORTRAN CODE OH–HE POTENTIAL
w r i t e ( ∗ , ∗ ) ’∗∗ R = ’ , R , ’ z = ’ , z , ’ r 1 = ’ , r 1
s t o p 242
endif
c
c
c
c
c
c
c
c
c
c
c
c
c a l l p _ a s s l e g 2 ( p l a , 6 , za ) ! p l a ( 1 ) = P2_2 ( za ) . . p l a ( 5 ) = P2_6 ( za )
c a l l p _ a s s l e g 2 ( p l b , 6 , zb ) ! p l b ( 1 ) = P2_2 ( zb ) . . p l a ( 5 ) = P2_6 ( zb )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) p l a ( 1 ) , p l a ( 2 ) , p l a ( 3 ) , p l a ( 4 ) , p l a ( 5 )
Ran ( 1 ) = d ex p ((− b e t a _ a )∗Ra )
Rbn ( 1 ) = d ex p ((− b e t a _ b )∗Rb )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ Ran ( 1 ) = ’ , Ran ( 1 ) , ’ Rbn ( 1 ) = ’ , Rbn ( 1 )
r 1 n _ I ( 1 ) = d ex p((− a l p h a _ I )∗ r 1 ∗∗3)
r 1 n _ I I ( 1 ) = d ex p((− a l p h a _ I I )∗ r 1 ∗∗3)
r 1 a n ( 1 ) = d ex p((− a l p h a _ a )∗ r 1 ∗∗3)
r 1 b n ( 1 ) = d ex p((− a l p h a _ b )∗ r 1 ∗∗3)
do i =2 ,6
Ran ( i ) = Ra∗Ran ( i −1)
Rbn ( i ) = Rb∗Rbn ( i −1)
r 1 a n ( i ) = r 1 ∗ r 1 a n ( i −1)
end do
do i =2 ,5
r 1 b n ( i ) = r 1 ∗r 1 b n ( i −1)
end do
do i =2 ,4
r 1 n _ I ( i ) = r 1 ∗ r 1 n _ I ( i −1)
r 1 n _ I I ( i ) = r 1∗ r 1 n _ I I ( i −1)
end do
t h e r1 i n d e p e n d e n t p a r t o f t h e s h o r t r a n g e
i n d =0
sum=0 d0
v=0 d0
i =1 ! l = 2 . . 2
i n d = i n d +1
sum=sum+ p l a ( i )∗ c_ v d ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ Ran ( 1 ) = ’ , Ran ( 1 )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ p l a ( i ) = ’ , p l a ( i )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ p l a ( i )∗ Ran ( 1 ) = ’ , p l a ( i )∗ Ran ( 1 )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ sum= ’ , sum
v=Ran ( 1 ) ∗ sum
i =1 ! l = 2 . . 2
i n d = i n d +1
sum=sum+ p l b ( i )∗ c_ v d ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ p l b ( i )∗ Rbn ( 1 ) = ’ , p l b ( i )∗ Rbn ( 1 )
v=v+Rbn ( 1 ) ∗ sum
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v= ’ , v
t h e r1 d e p e n d e n t p a r t o f t h e s h o r t r a n g e
do i =1 ,6
! n =0..5
do j =1 ,5
! l =2..6
sum=0 d0
do k =1 ,6 ! k = 0 . . . 5
i n d = i n d +1
sum=sum+ c_ v d ( i n d )∗ r 1 a n ( k )
end do
v=v+Ran ( i )∗ p l a ( j )∗ sum
end do
end do
do i =1 ,6
! n =0..5
do j =1 ,5
! l =2..6
sum=0 d0
do k =1 ,5 ! k = 0 . . 4
i n d = i n d +1
sum=sum+ c_ v d ( i n d )∗ r 1 b n ( k )
end do
v=v+Rbn ( i )∗ p l b ( j )∗ sum
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
c
c
c
the long range
c a l l p _ a s s l e g 2 ( p l , 9 , z ) ! p l ( 1 ) = P2_2 ( z ) . . p l ( 8 ) = P2_9 ( z )
do i =1 ,8
Rn ( i ) = f _ l o n g ( R , b e t a , 5 + i ) ! Rn ( 1 ) = f n ( B e t a∗R ) / R ^ 6 . . Rn ( 8 ) = f n ( B e t a∗R ) / R^13
end do
c
c
c
t h e l o n g ra n g e , t e r m s n = 1 0 . . 1 3
kmax=4
j 0 =1
do i =1 ,4
j 0=1− j 0
42
! k = 0..3
! n = i +9 ( n = 1 0 . . 1 3 −> i = 1 . . 4 )
! l m i n =0 ( ev en n ) , 1 ( odd n )
c
c
c
c
c
do j =1+ j 0 , i +4 ,2
! lmax = n−4 = i +4
i n d = i n d +1
sum= c_ v d ( i n d )
do k =1 , kmax
i n d = i n d +1
sum=sum+ r 1 n _ I ( k )∗ c_ v d ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i = ’ , i , ’ j = ’ , j , ’ k = ’ , k
end do
v=v+Rn ( i +4)∗ p l ( j )∗ sum ! s t a r t a t Rn ( 5 ) = f n ( B e t a∗R ) / R^10
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
t h e l o n g ra n g e , t e r m s n = 6 . . 9
j 0 =1
do i =1 ,4
! n = i +5 ( n = 6 . . 9 −> i = 1 . . 4 )
j 0 =1− j 0
! l m i n =0 ( ev en n ) , 1 ( odd n )
do j =1+ j 0 , i , 2
! lmax = n−4 = i
i n d = i n d +1
sum= c_ v d ( i n d )
do k =1 , kmax
i n d = i n d +1
sum=sum+ r 1 n _ I I ( k )∗ c_ v d ( i n d )
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i = ’ , i , ’ j = ’ , j , ’ k = ’ , k
end do
v=v+Rn ( i )∗ p l ( j )∗ sum ! s t a r t a t Rn ( 1 ) = f n ( B e t a∗R ) / R^6
end do
end do
i f ( d eb u g ) w r i t e ( ∗ , ∗ ) ’ i n d = ’ , i n d , ’ v = ’ , v
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c s u b r o u t i n e j a c 2 y u 2 ( ama , amb , R , z , r1 , Ra , za , Rb , zb )
c
c function :
c
t h i s r o u t i n e r e t u r n s t h e a n i s o t r o p i c v a l u e s Ra , za , Rb , Rz
c
o f a 3 atom s y s t e m w i t h 3 v a r i a b l e c o o r d i n a t e s : R , r , z=c o s ( t h e t a )
c
c input :
c
ama : u n i t mass o f atom a
c
amb : u n i t mass o f atom b
c
R:
d i s t a n c e f ro m atom c t o c e n t e r o f mass o f m o l e c u l e ab ( a_0 ) .
c
z:
e q u a l s c o s ( t h e t a ) , t h e t a =0 c o r r e s p o n d s t o t h e l i n e a r
c
g e o m e t r y ab−c
c
r1 : bond d i s t a n c e o f m o l e c u l e ab ( a_0 ) .
c
c output :
c
Ra : d i s t a n c e f ro m atom a t o atom c ( a_0 )
c
Rb : d i s t a n c e f ro m atom b t o atom c ( a_0 )
c
za : c o s ( t h e t a _ a ) , t h e t a _ a i s t h e a n g l e b e t w e e n Ra and
c
t h e i n t e r m o l e c u l a r a x i s o f m o l e c u l e ab
c
t h e t a _ a =0 c o r r e s p o n d s t o l i n e a r g e o m e t r y ab−c
c
zb : c o s ( t h a t a _ b ) , t h e t a _ b i s t h e a n g l e b e t w e e n Rb and
c
t h e i n t e r m o l e c u l a r a x i s o f m o l e c u l e ab
c
t h e t a _ b =0 c o r r e s p o n d s t o l i n e a r g e o m e t r y c−ab
c
s u b r o u t i n e j a c 2 y u 2 ( ama , amb , R , z , r 1 , Ra , za , Rb , zb )
i m p l i c i t none
r e a l ∗8 ama , amb , R , r 1 , z , Ra , za , Rb , zb
r e a l ∗8 r 1 a , r 1 b , a1 , b1 , tmpa , tmpb , x , y2 , a m t o t
c
a m t o t =ama+amb
r 1 a = r 1 ∗amb / a m t o t
r 1 b = r 1 ∗ama / a m t o t
i f (R . eq . 0 d0 ) t h en
za=−1d0
zb=−1d0
else
i f ( z . eq . 1 d0 ) t h en
za =1 d0
zb=−1d0
e l s e i f ( z . eq .−1d0 ) t h en
za=−1d0
zb =1 d0
else
a1 = r 1 a / R
b1 = r 1 b / R
tmpa= d s q r t ( 1 d0−2d0∗a1∗z+ a1∗a1 )
tmpb = d s q r t ( 1 d0 +2 d0∗b1∗z+b1∗b1 )
i f ( tmpa . eq . 0 d0 ) t h en
za =1 d0
else
43
APPENDIX B. FORTRAN CODE OH–HE POTENTIAL
za =( z−a1 ) / tmpa
end i f
i f ( tmpb . eq . 0 d0 ) t h en
zb =1 d0
else
zb =(−( z+b1 ) ) / tmpb
end i f
end i f
x=R∗z
y2=R∗R∗(1 d0−z∗z )
Ra= d s q r t ( ( x−r 1 a )∗∗2 + y2 )
Rb= d s q r t ( ( x+ r 1 b )∗∗2 + y2 )
endif
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c s u b r o u t i n e p _ l e g ( p , lmax , z )
c
c function :
c
r o u t i n e r e t u r n s t h e L e g e n d r e p o l y n o m i a l s P_ l ( z ) f o r l = 0 . . l m a x
c
c input :
c
p : a r r a y f o r P_ l ( z ) , s i z e s h o u l d be l m a x+1
c
l m a x : u p p e r l i m i t o f l , g e n e r a t e P_0 ( z ) . . P_lmax ( z )
c
z : t h e d e p e n d e n t v a r i a b l e o f P_ l ( z )
c
c output :
c
p ( 1 ) = P_0 ( z ) , p ( 2 ) = P_1 ( z ) . . p ( l m a x+1)= P_lmax ( z )
c
c error :
c
when l m a x < 0
c
s u b r o u t i n e p _ l e g ( p , lmax , z )
i m p l i c i t none
i n t e g e r lmax , i
r e a l ∗8 p ( 0 : lmax ) , z , d i
c
i f ( lmax . l t . 0 ) s t o p ’ l l e s s t h a n z e r o , r o u t i n e p _ l e g ’
p ( 0 ) = 1 . d0
i f ( lmax . eq . 0 ) r e t u r n
p(1) = z
d i = 0 d0
do i =1 , lmax−1
d i = d i + 1 d0
p ( i +1 ) = ( ( 2 . d0∗ d i + 1 . d0 ) ∗ z ∗ p ( i )
+
− d i ∗ p ( i −1) ) / ( d i + 1 d0 )
end do
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c s u b r o u t i n e p _ a s s l e g 2 ( p , lmax , z )
c
c function :
c
r o u t i n e r e t u r n s Racah n o r m a l i z e d a s s o c i a t e d L e g e n d r e f u n c t i o n s
c
P2 _ l ( z ) f o r l = 2 . . l m a x
c
c input :
c
p : a r r a y f o r P2 _ l ( z ) , s i z e s h o u l d be lmax−1
c
l m a x : u p p e r l i m i t o f l , g e n e r a t e P2_2 ( z ) . . P2_lmax ( z )
c
z : t h e d e p e n d e n t v a r i a b l e o f P2 _ l ( z )
c
c output :
c
p ( 1 ) = P2_2 ( z ) , p ( 2 ) = P2_3 ( z ) . . p ( lmax −1) = P2_lmax ( z )
c
c error :
c
when l m a x < 2
c
s u b r o u t i n e p _ a s s l e g 2 ( p , lmax , z )
i m p l i c i t none
i n t e g e r lmax , i
r e a l ∗8 p ( 2 : lmax ) , z
i f ( lmax . l t . 2 ) s t o p ’ l s m a l l e r t h a n 2 , r o u t i n e p _ a s s l e g 2 ’
do i =2 , lmax
call sf_cl2 ( i , z , p ( i ))
end do
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
c
r e t u r n r a c a h n o r m a l i z e d s p h e r i c a l h a rm o n i c , C_LM ( Z , p h i )
c
f o r M=2 , p h i =0 , and Z=c o s ( t h e t a )
44
c
c
c
GCG, 11−Jun −2006
subroutine s f _ c l 2 ( l , z , cl2 )
i m p l i c i t none
integer l , m
r e a l ∗8 z , p , c l 2 , a l
d a t a m/ 2 /
c
c a l l a s s l e g p l ( l , m, z , p )
al = dble ( l )
c l 2 = p / d s q r t ( a l ∗( a l ∗ a l −1d0 ) ∗ ( a l +2 d0 ) )
return
end
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
s u b r o u t i n e a s s l e g p l ( L ,M, X , P )
C
Compute a s s o c i a t e d L e g e n d r e p o l y n o m i a l s .
C
From N u m e r i c a l R e c i p e s , p . 182
C
IMPLICIT REAL∗8(A−H, O−Z )
C
IF ( (M. LT . 0 ) . OR. (M. GT . L ) . OR . (DABS (X ) . GT . 1 . D0 ) )
.
s t o p ’ bad a r g u m e n t s , a s s l e g p l ’
PMM= 1 . D0
IF (M. GT . 0 ) THEN
SOMX2=SQRT ( ( 1 . D0−X ) ∗ ( 1 . D0+X ) )
F f = 1 . D0
DO 11 I =1 ,M
PMM=(−PMM)∗ F f∗SOMX2
F f = F f + 2 . D0
11
CONTINUE
ENDIF
IF ( L . EQ .M) THEN
P=PMM
ELSE
PMMP1=X∗DBLE(2∗M+1)∗PMM
IF ( L . EQ .M+1 ) THEN
P=PMMP1
ELSE
DO 12 LL=M+2 ,L
PLL =(X∗DBLE(2∗LL−1)∗PMMP1−DBLE( LL+M−1)∗PMM) / DBLE( LL−M)
PMM=PMMP1
PMMP1=PLL
12
CONTINUE
P=PLL
ENDIF
ENDIF
RETURN
END
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
r e a l ∗8 f u n c t i o n f _ l o n g ( r , b e t a , n )
i m p l i c i t none
r e a l ∗8 r , b e t a
integer n
c
c
f _ l o n g ( r , b e t a , n )= d_n ( b e t a ∗r , n ) / r ^ n
c
c
Tang−T o e n n i e s damping f u n c t i o n :
c
d_n ( x , n ) = 1 − exp(−x ) ∗ sum_k =0:n x ^ k / k !
c
c
For s m a l l x ( when d_n <1d−8) u s e :
c
d_n ( x , n ) = exp(−x )∗ sum_k=n +1: i n f t y x ^ k / k !
c
r e a l ∗8 sum , t er m , dex , r n , dn , x
integer i
c
i f ( r . eq . 0 d0 ) t h en
f _ l o n g =0 d0
return
endif
c
x= b e t a ∗ r
sum=1 d0
t e r m =1 d0
r n =1 d0
do i =1 , n
rn = rn∗r
t e r m = t e r m∗x / d f l o a t ( i )
sum=sum+ t e r m
end do
45
APPENDIX B. FORTRAN CODE OH–HE POTENTIAL
dex =d ex p(−x )
dn =1 d0−dex∗sum
i f ( d a b s ( dn ) . l t . 1 d−8) t h en
c
c
c
a l t e r n a t i v e formula f o r small x
sum=0 d0
do i =n +1 ,1 0 0 0
t e r m = t e r m∗x / d f l o a t ( i )
sum=sum+ t e r m
i f ( t e r m / sum . l t . 1d−8) go t o 10
enddo
w r i t e ( ∗ , ∗ ) ’∗∗ ERROR i n F_LONG : no c o n v e r g e n c e , x= ’ , x , ’ n= ’ , n
s t o p 232
continue
dn=sum∗dex
endif
f _ l o n g =dn / r n
end
10
c
∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
c
Listing B.3: Vsum lin. coef. data, lc_vsum.f (20 lines)
r e a l ∗8 c _ v s ( 6 4 8 )
data c_vs /
+
3 .7 4 1 4 6 3 5 8 7 6 1 4 9 8 1 0 6 5 6 0 d +0 2 ,
+
3 .2 7 0 7 5 8 6 3 8 1 5 9 8 3 1 8 3 2 5 8 d +0 1 ,
+
6 .0 3 2 6 4 8 5 3 7 9 2 3 2 6 6 1 2 1 6 5 d +0 4 ,
+
6 .0 0 1 9 4 0 5 1 8 0 5 8 2 5 2 0 5 5 2 0 d +0 5 ,
+
4 .5 0 1 2 5 0 9 6 9 9 1 9 2 1 9 4 9 6 6 6 d +0 5 ,
+
4 .3 5 6 1 7 2 8 0 2 5 0 7 9 9 5 2 2 0 3 7 d +0 4 ,
+
3 .1 2 9 7 6 0 5 7 4 5 4 9 0 1 4 7 9 1 0 5 d +0 2 ,
+
7 .6 5 1 2 2 8 8 9 3 8 7 8 0 0 1 2 5 4 0 5 d +0 5 ,
−4.34676007946692891437 d +0 2 ,
−8.92734935833879177380 d +0 1 ,
−2.95500960259121959098 d +0 5 ,
−6.74209337257708190009 d +0 5 ,
−1.81842675528502295492 d +0 5 ,
−5.69456226750027053640 d +0 3 ,
−1.53956676653732836712 d +0 5 ,
−1.57017876715056924149 d +0 6 ,
+
+
+
+
+
+
+
+
+
+
1 .5 3 5 4 5 4 7 5 7 6 5 0 4 3 1 6 4 5 2 9 d +0 2 ,
−1.42255938894606089207 d +0 2 ,
−9.81441775217014765076 d +0 3 ,
1 .8 0 3 0 6 0 3 6 0 7 4 1 4 0 5 0 2 0 2 4 d +0 4 ,
3 .9 1 1 4 3 5 4 6 8 1 8 8 9 5 1 7 4 5 1 1 d +0 3 ,
6 .2 9 8 5 9 2 1 7 8 0 0 6 6 2 0 9 9 7 6 2 d +0 3 ,
−1.17055436756619510561 d +0 4 ,
−1.65767060106257235930 d +0 3 ,
1 .1 0 0 8 6 6 7 8 7 9 7 3 9 4 3 4 7 3 8 4 d +0 3 ,
3 .2 8 5 7 6 0 7 3 9 5 7 0 5 7 8 1 5 2 7 8 d + 0 2 /
−2.70358512491298313307 d +0 2 ,
2 .2 9 9 8 2 4 3 0 7 0 2 2 0 0 9 5 7 7 3 9 d +0 2 ,
3 .8 5 7 9 0 0 6 8 7 4 9 8 3 3 3 5 2 5 1 3 d +0 1 ,
1 .0 0 9 2 0 6 4 9 5 1 9 1 5 4 1 0 3 4 7 3 d +0 4 ,
−1.16956379732742789201 d +0 4 ,
−1.25611957111754491052 d +0 4 ,
1 .5 7 0 6 9 5 8 9 1 6 2 3 6 5 2 8 2 3 9 1 d +0 4 ,
3 .4 6 9 3 0 3 5 9 8 2 4 1 2 4 8 8 9 2 0 1 d +0 3 ,
1 .4 2 3 4 2 8 2 2 4 1 9 4 0 6 3 6 5 9 5 2 d +0 3 ,
−9.90865270108422237172 d +0 2 ,
Listing B.4: Vdiff lin. coef. data lc_vdif.f (20 lines)
46
r e a l ∗8 c_ v d ( 4 3 2 )
d a t a c_ v d /
+
−1.58614844142719556286 d −02,
+
3 .0 3 7 5 2 5 7 1 0 4 3 6 9 8 3 4 9 4 8 4 d −01,
+
−2.60879984118389030012 d −01,
+
−2.01792537772616485670 d −02,
+
6 .7 8 4 6 8 5 8 9 9 1 0 6 7 0 0 2 2 4 2 5 d −02,
+
1 .7 6 7 5 0 2 1 7 3 5 8 0 4 4 9 4 2 0 5 1 d −01,
+
−9.64116191052472037981 d −03,
+
−3.07034882653377284856 d −02,
−1.66979092649802311144 d +0 2 ,
1 .2 6 0 1 6 0 2 4 9 6 1 5 8 7 6 7 4 8 5 7 d −01,
1 .1 4 3 0 6 7 7 5 0 5 2 3 0 1 1 2 1 1 2 1 d −01,
1 .0 3 2 1 5 3 4 4 7 7 8 8 4 1 9 2 5 5 1 6 d −03,
−3.23690623283119982556 d −01,
−2.06713783527393765715 d −02,
2 .4 5 4 4 2 4 5 5 5 9 2 5 8 8 7 1 1 0 8 9 d −03,
1 .9 7 2 2 6 3 0 4 8 2 1 0 5 7 2 4 8 0 8 5 d −02,
+
+
+
+
+
+
+
+
+
+
9 .7 9 5 5 1 4 2 2 7 0 0 0 7 4 8 5 1 9 4 6 d +0 1 ,
1 .0 0 4 3 5 9 2 4 7 6 3 6 0 0 6 2 4 4 8 7 d +0 2 ,
8 .4 6 6 1 8 5 5 7 6 3 5 8 3 0 6 5 9 7 0 2 d −02,
2 .4 4 8 9 5 4 5 6 4 9 4 2 8 8 0 5 5 1 8 3 d +0 1 ,
4 .7 8 0 4 1 2 0 9 9 2 9 8 6 5 8 2 9 0 5 4 d +0 0 ,
−2.75567466972600598751 d +0 1 ,
−3.35659188798050877267 d +0 2 ,
2 .8 2 8 2 0 7 8 4 2 2 5 4 8 1 9 9 3 7 0 1 d +0 1 ,
1 .0 6 8 5 8 6 3 4 9 3 1 9 5 1 8 4 0 4 8 5 d +0 2 ,
2 .0 0 4 3 1 5 3 4 9 8 9 0 7 5 4 6 6 0 5 8 d + 0 1 /
−1.12747648523908182483 d +0 2 ,
−1.65917477688537047698 d +0 2 ,
−2.46078771062091270494 d +0 1 ,
−1.02474839715459395961 d +0 1 ,
−1.82363851106500476362 d +0 1 ,
−1.92418330360481888874 d +0 2 ,
3 .6 9 4 3 4 2 5 6 4 3 0 2 4 8 0 6 9 5 9 1 d +0 2 ,
6 .4 6 5 1 1 0 7 3 8 0 3 8 9 1 9 1 8 0 3 3 d +0 1 ,
−6.80343405903647209243 d +0 1 ,
−8.50705387138675916958 d +0 1 ,
Three-dimensional potentials of the X̃ 2 Π state of the OH–He complex, May 8, 2012