THE JOURNAL OF CHEMICAL PHYSICS 139, 194309 (2013)
Exact quantum scattering calculations of transport properties for the H2 O–H system Paul J. Dagdigian1,a) and Millard H. Alexander2,b) 1
Department of Chemistry, The Johns Hopkins University, Baltimore, Maryland 21218-2685, USA Department of Chemistry and Biochemistry and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742-2021, USA 2
(Received 28 September 2013; accepted 28 October 2013; published online 19 November 2013) Transport properties for collisions of water with hydrogen atoms are computed by means of exact quantum scattering calculations. For this purpose, a potential energy surface (PES) was computed for the interaction of rigid H2 O, frozen at its equilibrium geometry, with a hydrogen atom, using a coupled-cluster method that includes all singles and doubles excitations, as well as perturbative contributions of connected triple excitations. To investigate the importance of the anisotropy of the PES on transport properties, calculations were performed with the full potential and with the spherical average of the PES. We also explored the determination of the spherical average of the PES from radial cuts in six directions parallel and perpendicular to the C2 axis of the molecule. Finally, the computed transport properties were compared with those computed with a Lennard-Jones 12-6 potential. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829681] I. INTRODUCTION
The chemical kinetic modeling of combustion systems, such as flames and combustion engines, requires not only temperature-, and for some reactions, pressure-dependent reaction rate constants, but also transport properties because the species concentrations will have a spatial dependence in such systems. While considerable effort has been made to compute transport properties, as discussed in a recent review,1 there is a need for more accurate values of these parameters. The overwhelming majority of transport property calculations have employed a classical treatment with parameterized, analytic potentials, such as the Lennard-Jones (LJ) 12-6 potential. There have been a number of quantum scattering calculations of molecular transport properties employing accurate ab initio potentials.2–11 Other quantum scattering transport property calculations have employed the spherical average of such potentials.12–15 Improvements in computational methods and resources permit the calculation of state-of-the-art potential energy surfaces (PESs) for systems involving both stable species and reactive intermediates.16, 17 In recent work,10, 11 we have explored the calculation of transport properties for two atom-molecule systems involving reactive intermediates [OH(X2 )–He and CH2 (X˜ 3 B1 , a˜ 1 A1 )– He] through exact quantum scattering calculations with stateof-the-art anisotropic PES’s. For the former system, only slight differences (3%–5%) in the computed transport properties were found between those computed with the full PES and truncated potentials (spherical average of the sum potential Vsum or the difference potential Vdif set to zero).10 More significant differences were found for this system between the calculations with the full PES and those computed with a LJ 12-6 potential. a) Electronic mail: [email protected] b) Electronic mail: [email protected]
0021-9606/2013/139(19)/194309/8/$30.00
The methylene CH2 –He system offered an opportunity to investigate in more detail the effect of the PES on transport properties since the anisotropies of the CH2 (X˜ 3 B1 )–He and CH2 (a˜ 1 A1 )–He systems are dramatically different. The ˜ is very anisotropic because of the strong PES for CH2 (a)–He interaction of the electrons on the helium atom with the unoccupied CH2 orbital perpendicular to the molecular plane.18 By ˜ is much contrast, the anisotropy of the PES for CH2 (X)–He smaller because this orbital is singly occupied in this electronic state.19 These differences were found to yield transport properties of different magnitudes for the two CH2 electronic states, ranging from 3% to 15% higher for the a˜ state at the lowest and highest temperatures considered (200 to 1500 K, respectively).11 More significant differences were observed for the transport properties computed with the full PES’s and those computed with LJ potentials whose parameters were derived by combination rules.20 Here, we carry out exact quantum scattering calculations for collisions of the water molecule with the hydrogen atom. Both species have significant concentrations in combustion media. Moreover, the transport of hydrogen is expected to be rapid because of its small mass. The PES for the H2 O–H interaction was computed more than 20 years ago by Zhang et al.21 using several different methods [the unrestricted open-shell Hartree-Fock method, second-order Møller-Plesset perturbation theory (MP2), and fourth-order Møller-Plesset perturbation theory (MP4)]. These authors determined a spherically averaged, isotropic potential from radial cuts in six directions parallel and perpendicular to the C2 axis of the molecule. In our approach, we have computed the interaction energy for many more angular orientations of the H atom and have determined a full PES, using a coupled-cluster theory with a reasonably large basis set. Note that our study focusses on the repulsive ground electronic state of the H3 O system, not the bound excited state in which one of the electrons is promoted into a Rydberg orbital on the O atom.22
139, 194309-1
© 2013 AIP Publishing LLC
194309-2
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
At sufficiently high collision energy (or, equivalently, temperature) or with high vibrational excitation of the molecule, H + H2 O can react to form H2 + OH products. The activation energy for this reaction is high, and the rate constant is only approximately 3 × 10−15 cm3 molecule−1 s−1 at 1000 K.23 Two different types of experimental studies of the reaction dynamics have been carried out on this system, namely reaction of H2 O in its ground vibrational level with photolytically generated hyperthermal H/D atoms24–26 and reaction of thermal H/D atoms with vibrationally excited water molecules.27–29 The theoretical work has involved calculation of a potential energy surface30, 31 for the reaction and quasiclassical trajectory and time-dependent quantum mechanical wave packet calculations of cross sections.26, 32, 33
II. CALCULATION OF TRANSPORT PROPERTIES
In this section, a brief review of the calculation of transport properties of molecule immersed dilutely in a bath of a collision partner is presented. Transport properties can be computed from collision integrals (n, s) (T), which can be obtained by taking an appropriate integral over the collision energy and carrying out a Boltzmann state average.4, 20, 34, 35 Reversing the usual order of carrying out this integration and state averaging, we can express the collision integral as a Boltzmann state average over state-dependent collision integrals (n,s) ji (T ): 1 (2ji + 1) exp(−εi /kB T )(n,s) (n,s) (T ) = ji (T ). (1) qR j i
In Eq. (1), εi is the energy of the ith rotational level, qR is the rotational partition function, and kB is the Boltzmann constant. The state-dependent collision integrals are the integrals of state-dependent effective cross sections Q(n) ji (E) over the collision energy E: (n,s) ji
kB T 1/2 1 2π μ (kB T )s+2 ∞ × E s+1 exp(−E/kB T )Q(n) ji (E)dE.
1 = 2
(2)
0
Here, μ is the atom-molecule reduced mass. The state-dependent effective cross section in Eq. (2) is a sum over final levels jf of state-to-state effective cross sections: Q(n) ji (E) =
Q(n) ji →jf (E).
(3)
where Rˆ = (θ, φ) is the orientation of the Jacobi vector R. The weighting factors in Eq. (4) for n = 1 and 2 are4, 11, 34 1 (E) = 1 − (E /E)1/2 P1 (cos θ ), 2 (E) =
5 6
− 16 (E /E)2 − 23 (E /E)P2 (cos θ ),
(5) (6)
where E is the relative translational energy after scattering into level jf . It should be noted that for integral cross sections the weighting factor in Eq. (4) is 0 (E) = 1. Details about the evaluation of the Legendre moments in the effective cross sections in Eq. (4) can be found in the literature.36, 37 The effective cross sections in Eq. (4) can be written as the sum of several low-order Legendre moments of the differential cross sections: Q(1) ji −jf =
Q(2) ji −jf =
π (−1)ji −jf A0 − 13 (E /E)1/2 A1 , (2ji + 1)ki2
π (−1)ji −jf 5 1 − 6 (E /E)2 A0 − 6 2 (2ji + 1)ki
(7)
2 (E /E)A2 15
.
(8) In Eqs. (7) and (8), ki is the initial wave vector and Aλ equals Aλ = Z(l1 J l2 J ; ji λ)Z(l1 J l2 J ; jf λ) J J l1 l2 l1 l2
∗ × TjJi l1 ,jf l TjJi l2 ,jf l . 1
2
(9)
In Eq. (9), T designates a T-matrix element, expressed in the space-fixed frame, and the Z coefficients36, 38 are given by 1 1/2
Z(lJ l J ; j λ) = (−1) 2 (λ−l+l )+J +J [λ][l][l ][J ][J ]
×
l l λ 000
l J j , J l λ
(10)
where [x] = 2x + 1, ·· ·· ·· is a 3j symbol, and ·· ·· ·· is a 6j symbol.39 We employed Eqs. (3) and (7)–(10) to compute the statedependent effective cross sections in the scattering calculations described in Sec. IV. These cross sections were used in Eq. (2) to derive state-dependent collision integrals. A Boltzmann state average of these integrals was finally carried out [Eq. (1)] in order to compute collision integrals. We can relate the collision integrals to two transport properties, namely the binary diffusion coefficient Dab and the quantity ηab that determines the curvature of the mixture viscosity as a function of concentration. The relations are20, 34, 35
jf
The state-to-state effective cross sections on the righthand side of Eq. (3) are weighted angle averages of the ji → jf differential cross section: Q(n) ji →jf (E)
=
dσ d
ˆ n (E)d R, ji →jf
(4)
Dab =
3kB T , 16N μ(1,1)
(11)
5kB T , 8(2,2)
(12)
ηab =
where N is the total number density of the gas.
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
III. POTENTIAL ENERGY SURFACE
−4
6
−5 0 −53 −5
120
−474 −4
−56
−41
−5 −60 6 −47 −50 −5 3 −4 4
−41
−4 4
−50
−44
−53
−47
60 −53
−44
−47
φ / degrees
−53 −5 −47 0
−56
−56
180
−60
−47 3 −50 −5
−41
4
−
4
−41
240 56
−4
−53
−44
−50
−47
−47
−50
300
−53
The coordinates used to describe the H2 O–H PES are the same as those used previously in our study of the interaction of the He atom with both H2 O and CH2 .18, 19, 40 As depicted in Fig. 1 of Ref. 18, the body-frame z axis is chosen to lie along the a inertial axis of the molecule, and the molecule lies in the xz plane. The water molecule was taken to be rigid, with bond length and bond angle (r = 1.8361 bohr and θ = 104.69◦ , respectively) corresponding to the vibrationally averaged geometry reported by Mas and Szalewicz.41 The PES for the H2 O–H system was computed with the MOLPRO 2010.1 suite of computer codes,42 by means of restricted coupled cluster calculations with inclusion of single and double excitations, as well as perturbative contributions of connected triple excitations [RCCSD(T)].43, 44 An atomcentered avqz atomic-orbital basis was used,45 with the addition of a set of bond functions located midway between the atom and the molecule.46, 47 A counterpoise correction was used at all geometries to correct for basis set superposition error (BSSE).48 No hint of access to the reactive channel leading to OH + H2 products was found in the calculations; this is due to the fact that the H2 O geometry was held fixed. Since we are interested in collisions at thermal energies and only slightly above, it is reasonable to neglect the reactive channel. The PES was determined on a grid of 20 values of the atom-molecule separation R ranging from 3.0 to 10 bohr in steps of 0.5, with additional points at 11, 12, 13, 15, and 20 bohr. By exploiting symmetry we can limit the calculations to the angular range 0◦ ≤ θ ≤ 90◦ and 0◦ ≤ φ ≤ 180◦ , both varied in steps of 10◦ , for a total of 190 orientations. The total number of nuclear geometries at which the interaction energy was computed was 3800. The computed potential as a function of R and the orientation was fit to the following expansion in spherical harmonics: Vλμ (R)(1 + δμ0 )−1 V (R, θ, φ) =
360
−50
194309-3
0 0
30
60
90 120 θ / degrees
150
180
z H
R y
θ=30o
x H O
φ=180o H
FIG. 1. (Top panel) Dependence of the potential energy (in cm−1 ) on the orientation of the hydrogen atom collision partner with respect to the H2 O molecule for an atom-molecule separation R = 6.5 bohr. (Bottom panel) Arrangement of the atoms at the H2 O–H equilibrium geometry [R = 6.50 bohr, θ = 30◦ , φ = 180◦ ].
λ,μ≥0
× [Yλμ (θ, φ) + (−1)μ Yλ,−μ (θ, φ)]. (13) Symmetry considerations restrict the sums in Eq. (13) to terms with λ + μ even.18 To get an acceptable fit of the PES, we needed to include all terms with λ ≤ 10 and μ ≤ 8, for a total of 33 anisotropic terms. Figure 1 presents a contour plot of the dependence of the potential energy upon the orientation of the H atom collision partner at an atom-molecule separation R = 6.50 bohr, the separation at the global minimum. The anisotropy of the PES at this value of R is seen to be small. This is true for most values of R, except at small atom-molecule separations for which there is significant repulsion between the H atom collision partner and the H atoms on the water molecule. The global minimum of the PES has an energy of 61.0 cm−1 below the H2 O + H asymptote, at a geometry of R = 6.50 bohr, θ = 30◦ , φ = 180◦ . The hydrogen atom collision partner thus lies within the molecular plane at the global minimum and close to one of the H atoms of the water molecule. Because of the small anisotropy of the PES, the angular orientation for the minimum interaction energy varies somewhat
as a function of R. The binding energy for H2 O–H is somewhat larger than for the related H2 O–He system, for which the binding energy is computed to be 35 cm−1 , while the atom-molecule separation R for the latter system is slightly smaller (5.92 bohr).49 In their calculations on H2 O–H, Zhang et al.21 considered that their best estimate for the well depth of this system (53 ± 6 cm−1 ) was intermediate between the values computed in MP4 calculations without and with the counterpoise correction (69.9 and 36.8 cm−1 , respectively). The presently computed well depth of 61 cm−1 for H2 O–H is slightly larger than the estimate of Zhang et al.21 Figure 2 presents a plot of the larger expansion coefficients Vλμ as a function of the atom-molecule separation R in the region of the van der Waals well. The largest coefficient in the region of the well is the isotropic V00 term. We see that for values of R smaller than that of the global minimum some of the anisotropic coefficients become compatible to or larger than the V00 term. In their investigation of the interaction of H2 O with H, Zhang et al.21 derived a spherically averaged, isotropic potential by taking the average of radial cuts along six orientations
194309-4
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
with a small nonlinear polyatomic molecule by averaging cuts along six orthogonal approach directions of the perturber to the molecule, at least for the H2 O–H system. A similar approach has been employed by Stallcop et al.12, 14 for several atom-molecule systems. Jasper and Miller50 have considered methods to spherically average potentials involving larger, floppy molecules. Also plotted in Fig. 3 is the spherical average obtained by Zhang et al.21 by taking the average of MP4 calculations with a BSSE correction along the ±x, ±y, and ±z directions. This potential energy curve has a slightly shallower well and longer equilibrium atom-molecule separation than the present calculation.
200 150
energy / cm
−1
100
V20
V31
50
V11
V22
0 −50 −100 V00
−150 −200 4
5
6
7 R / bohr
8
9
10
FIG. 2. Dependence of the larger expansion coefficients Vλμ [defined in Eq. (13)] on the atom-molecule separation R.
parallel and perpendicular to the C2v axis of the molecule. With the availability of our full PES for this system, we can check the accuracy of this method to generate a spherically averaged H2 O–H potential. Figure 3 presents a plot of the spherical average of our full PES [which equals (4π )−1/2 V00 ]. A plot of the average of the radial cuts of our PES along the ±x, ±y, and ±z directions is indistinguishable from the spherical average on the scale of this figure. At the minimum of the potential, the two curves differ by 0.5 cm−1 , while the well depth of the former equals 48.6 cm−1 . At 4 bohr, the potentials differ by 31 cm−1 , as compared to an interaction energy of 1240 cm−1 . As an indication of the small anisotropy for this system, we note that the well depth for the spherically averaged potential is only slightly smaller than the well depth for the full PES. This comparison of the two ways to obtain a spherical average of the potential suggests that it may be reasonable to compute an isotropic potential for the interaction of an atom
100
energy / cm−1
present work Zhang et al.
50
0
−50 4
5
6
7 R / bohr
8
9
10
FIG. 3. Spherically averaged H2 O–H potentials determined by averaging radial cuts along the ±x, ±y, and ±z directions. Red: present work; blue: MP4 calculation with counterpoise correction by Zhang et al. (Ref. 21).
IV. SCATTERING CALCULATIONS
Close-coupling calculations were carried out with the HIBRIDON suite of programs,51 which has recently been extended to allow the determination of the Q1 and Q2 effective cross sections [Eqs. (7) and (8)]. The cross sections were checked for convergence with respect to the inclusion of a sufficient number of partial waves and energetically closed channels. At the highest energies for which calculations were carried out (4000 cm−1 ), the rotational basis included all levels whose energy was less than 4600 cm−1 , and the scattering calculations included all total angular momenta J ≤ 133 ¯. In previous work,11 we checked that the computed effective cross sections agreed with values obtained by numerical integration of the differential cross sections weighted by the 1 (E) and 2 (E) factors [Eqs. (5) and (6)]. The high-energy tails of the integrals over collision energy [Eq. (2))] were evaluated by exponential extrapolation of the integrand. The sum over rotational levels ji in Eq. (1) included all levels whose energy was less than 3000 cm−1 . For levels whose energy is between 1200 and 3000 cm−1 , the total effective cross section [Eq. (4)] was taken to be equal to that of the highest rotational level of the same ka projection number with energy less than 1200 cm−1 .
V. RESULTS
The rotational levels of an asymmetric top like water are labeled jka kc , where j is the rotational angular momentum, ka is its projection along the principal axis in the prolate limit, and kc is its projection along the principal axis in the oblate limit.52 Each rotational level exists for only one of the nuclear spin modifications (ortho and para). The two nuclear spin modifications cannot be easily interconverted in nonreactive collisions. Figure 4 displays plots of the collision energy dependence of the state-dependent effective cross sections Q(1) ji and [Eq. (3)] for collisions of para-H O in its lowest and a Q(2) 2 ji higher rotational level [000 and 551 , respectively]. The rotational energy of the latter is 758 cm−1 . The elastic contribution to these effective cross sections is also plotted. As we ˜ a)–He ˜ collisions,10, 11 the observed in OH(X)–He and CH2 (X, (1) Qji effective cross sections decrease monotonically with increasing collision energy.
P. J. Dagdigian and M. H. Alexander
total
V
00
term only
16
i
3
/ cm 10
18 17
Ω 11
(d) 551 V00
(b) 1500 K
s
–11 –1
12
V00 term only
3
13
i
elastic
/ cm 10
total
–11 –1
19
14
(1,1) j
V00
20
(a) 300 K
s
(2)
(c) 551
i
elastic
Q j effective cross section / Å2
total
i
(1)
Q j effective cross section / Å2
V00
V00
15
(b) 000
(1,1) j
(a) 000
J. Chem. Phys. 139, 194309 (2013)
Ω
194309-5
15 14
0 2 4 6 8 10 rotational angular momentum j
0 2 4 6 8 10 rotational angular momentum j
total elastic
elastic
collision energy / cm–1
FIG. 4. Dependence on collision energy of the state-dependent effective (1) (2) cross sections Qji [panels (a) and (c)] and Qji [panels (b) and (d)] for collisions of H2 O in the 000 and 551 initial levels with hydrogen atoms. The elastic contribution to the effective cross section is indicated in each panel. In addition, the effective cross section computed with the spherical average of the potential is displayed for comparison.
At low collision energies, elastic collisions make by far the dominant contribution to the Q(1) ji effective cross section, while at higher energies inelastic scattering plays a much larger role. As we have noted previously,10, 11 this contrasts with the dominant role of elastic scattering in total integral cross sections. We showed in our study of CH2 –He transport properties11 that the differing importance of elastic scattering in the total integral and transport effective cross section results from the angular weighting of the differential cross section in the latter [see Eq. (4)]. Because of the weighting factor n (E) [Eqs. (5) and (6)], transport cross sections, particularly for elastic transitions, are much smaller than the corresponding integral cross sections, as discussed previously.11 Also included in the plots in Fig. 4 are the effective cross sections computed with the spherical average of the PES [(4π )−1/2 V00 ]. Of course, only elastic collisions can take place with this potential. We see that the Q(1) ji effective cross sections for this isotropic potential and the full PES are very similar in magnitude for all collision energies. We consider below the effect of the anisotropy of the PES through examination of state-dependent collision integrals [Eq. (1)]. It can be seen in Fig. 4 that the Q(2) ji effective cross section for the 000 level displays a similar collision energy dependence as for Q(1) ji for this level, namely a monotonically decreasing magnitude vs. collision energy. By contrast, the Q(2) ji effective cross section for the 551 level becomes very large and negative at low collision energies. As noted previously for Q(2) ji effective cross sections for higher rotational levels,10, 11 this is a consequence of large values for the factor E /E in the 2 (E) weighting factor in Eq. (6) for superelastic collisions (E > E) at small collision energies. As we have discussed previously,10, 11 this singularity does not cause problems in the evaluation of the Q(2) ji effective cross sections for the higher rotational levels, and hence also the determination
(1,1)
FIG. 5. State-dependent collision integrals ji as a function of the rotational angular momentum j for H2 O–H at temperatures of (a) 300 K and (b) 1500 K. Red and blue symbols denote ortho and para levels, respectively. Levels with different values of the body-frame projection quantum number are denoted by the following symbols: ka = 0, circles; ka = 1, plus signs; ka = 2, diamonds; ka = 3, upward-pointing triangles; ka = 4, downwardpointing triangles; ka = 5, left-pointing triangles; and ka = 6, right-pointing triangles. In each panel, the value of the collision integral computed with the spherical average of the potential is indicated with a dashed line.
of the (2,2) collision integral. As with the Q(1) ji effective cross sections, elastic collisions make a decreasing contribution to Q(2) ji effective cross sections with increasing collision energy. Finally, we note that the values of Q(2) ji computed with the spherical average of the potential are very similar to the values computed with the full PES, except at the lowest collision energies for the higher rotational levels. From the state-dependent effective cross sections we can use Eq. (1) to compute the (1,1) and (2,2) collision integrals. We plot in Fig. 5 the variation with rotational level of the statecollision integrals at two temperatures. The dependent (1,1) ji (2,2) ji collision integrals display a similar dependence upon the rotational quantum numbers and hence are not plotted is plotted are those here. The rotational levels for which (1,1) ji for which collision-energy dependent effective cross sections have been explicitly computed (see Sec. IV). The most probable value of j is ∼3 at 300 K and ∼7 at 1500 K. We see that on the rotational quantum at 300 K the dependence of (1,1) ji numbers is weak. This dependence is stronger at 1500 K. In values are seen in Fig. 5(b) to show a departicular, the (1,1) ji crease with increasing j. Presumably, this reflects the greater anisotropy of the PES at higher energies, because of the repulsion of the hydrogen atom collision partner by the hydrogen atoms on the water molecule. We also plot in Fig. 5 the effective cross sections under the assumption that the interaction is purely spherical (only the V00 term); in this case there is no dependence on rotational level. As shown in Fig. 5, it is easy to use the spherically averaged potential to predict collision integrals. Figure 6 displays the temperature dependence of these quantities. We see that the collision integrals computed with the spherically averaged potential are slightly larger than those computed with the full PES. Specifically, the collision integrals computed with the former are ∼1.4% larger at 300 K and 6% larger at 1500 K. These differences are smaller than what we found ˜ a)–He ˜ systems, likely a reflection of the larger for the CH2 (X, anisotropies for the PESs for these systems.
194309-6
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
2.5
full PES
−1 4
full PES
Dab / 10 Torr cm s
V00 only
2 −1
2
V00 only 1.5
1
LJ 12–6 0.5
0 200
full PES
FIG. 6. The collision integrals (1,1) and (2,2) as a function of temperature for the H2 O–H system. These quantities have been computed for the full PES and for the spherically averaged potential.
VI. DISCUSSION
We have employed quantum scattering calculations to compute collision integrals for the H2 O–H system, from which transport properties can readily be predicted. As an example we will consider the H2 O–H diffusion coefficient. We can also compare these values with a conventional calculation using a parameterized LJ 12-6 potential. The parameters specifying this potential were estimated through conventional combining rules: To estimate the well depth ε and the length parameter σ for H2 O–H, we took the geometric and arithmetic means, respectively, of the corresponding parameters for the like systems [ε = 572.4 K and σ = 2.605 for H2 O– H2 O and ε = 145 K and σ = 2.050 for H–H].53 Since H2 O possesses a dipole moment, there is an induction contribution to the long-range attractive interaction. Using the dipole moment of water (1.844 D)53 and the polarizability of the H atom (0.663 3 ),54 we follow the procedure of Hirschfelder et al.20 to correct the LJ parameters; this correction is quite small (
Exact quantum scattering calculations of transport properties for the H2 O–H system Paul J. Dagdigian1,a) and Millard H. Alexander2,b) 1
Department of Chemistry, The Johns Hopkins University, Baltimore, Maryland 21218-2685, USA Department of Chemistry and Biochemistry and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742-2021, USA 2
(Received 28 September 2013; accepted 28 October 2013; published online 19 November 2013) Transport properties for collisions of water with hydrogen atoms are computed by means of exact quantum scattering calculations. For this purpose, a potential energy surface (PES) was computed for the interaction of rigid H2 O, frozen at its equilibrium geometry, with a hydrogen atom, using a coupled-cluster method that includes all singles and doubles excitations, as well as perturbative contributions of connected triple excitations. To investigate the importance of the anisotropy of the PES on transport properties, calculations were performed with the full potential and with the spherical average of the PES. We also explored the determination of the spherical average of the PES from radial cuts in six directions parallel and perpendicular to the C2 axis of the molecule. Finally, the computed transport properties were compared with those computed with a Lennard-Jones 12-6 potential. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4829681] I. INTRODUCTION
The chemical kinetic modeling of combustion systems, such as flames and combustion engines, requires not only temperature-, and for some reactions, pressure-dependent reaction rate constants, but also transport properties because the species concentrations will have a spatial dependence in such systems. While considerable effort has been made to compute transport properties, as discussed in a recent review,1 there is a need for more accurate values of these parameters. The overwhelming majority of transport property calculations have employed a classical treatment with parameterized, analytic potentials, such as the Lennard-Jones (LJ) 12-6 potential. There have been a number of quantum scattering calculations of molecular transport properties employing accurate ab initio potentials.2–11 Other quantum scattering transport property calculations have employed the spherical average of such potentials.12–15 Improvements in computational methods and resources permit the calculation of state-of-the-art potential energy surfaces (PESs) for systems involving both stable species and reactive intermediates.16, 17 In recent work,10, 11 we have explored the calculation of transport properties for two atom-molecule systems involving reactive intermediates [OH(X2 )–He and CH2 (X˜ 3 B1 , a˜ 1 A1 )– He] through exact quantum scattering calculations with stateof-the-art anisotropic PES’s. For the former system, only slight differences (3%–5%) in the computed transport properties were found between those computed with the full PES and truncated potentials (spherical average of the sum potential Vsum or the difference potential Vdif set to zero).10 More significant differences were found for this system between the calculations with the full PES and those computed with a LJ 12-6 potential. a) Electronic mail: [email protected] b) Electronic mail: [email protected]
0021-9606/2013/139(19)/194309/8/$30.00
The methylene CH2 –He system offered an opportunity to investigate in more detail the effect of the PES on transport properties since the anisotropies of the CH2 (X˜ 3 B1 )–He and CH2 (a˜ 1 A1 )–He systems are dramatically different. The ˜ is very anisotropic because of the strong PES for CH2 (a)–He interaction of the electrons on the helium atom with the unoccupied CH2 orbital perpendicular to the molecular plane.18 By ˜ is much contrast, the anisotropy of the PES for CH2 (X)–He smaller because this orbital is singly occupied in this electronic state.19 These differences were found to yield transport properties of different magnitudes for the two CH2 electronic states, ranging from 3% to 15% higher for the a˜ state at the lowest and highest temperatures considered (200 to 1500 K, respectively).11 More significant differences were observed for the transport properties computed with the full PES’s and those computed with LJ potentials whose parameters were derived by combination rules.20 Here, we carry out exact quantum scattering calculations for collisions of the water molecule with the hydrogen atom. Both species have significant concentrations in combustion media. Moreover, the transport of hydrogen is expected to be rapid because of its small mass. The PES for the H2 O–H interaction was computed more than 20 years ago by Zhang et al.21 using several different methods [the unrestricted open-shell Hartree-Fock method, second-order Møller-Plesset perturbation theory (MP2), and fourth-order Møller-Plesset perturbation theory (MP4)]. These authors determined a spherically averaged, isotropic potential from radial cuts in six directions parallel and perpendicular to the C2 axis of the molecule. In our approach, we have computed the interaction energy for many more angular orientations of the H atom and have determined a full PES, using a coupled-cluster theory with a reasonably large basis set. Note that our study focusses on the repulsive ground electronic state of the H3 O system, not the bound excited state in which one of the electrons is promoted into a Rydberg orbital on the O atom.22
139, 194309-1
© 2013 AIP Publishing LLC
194309-2
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
At sufficiently high collision energy (or, equivalently, temperature) or with high vibrational excitation of the molecule, H + H2 O can react to form H2 + OH products. The activation energy for this reaction is high, and the rate constant is only approximately 3 × 10−15 cm3 molecule−1 s−1 at 1000 K.23 Two different types of experimental studies of the reaction dynamics have been carried out on this system, namely reaction of H2 O in its ground vibrational level with photolytically generated hyperthermal H/D atoms24–26 and reaction of thermal H/D atoms with vibrationally excited water molecules.27–29 The theoretical work has involved calculation of a potential energy surface30, 31 for the reaction and quasiclassical trajectory and time-dependent quantum mechanical wave packet calculations of cross sections.26, 32, 33
II. CALCULATION OF TRANSPORT PROPERTIES
In this section, a brief review of the calculation of transport properties of molecule immersed dilutely in a bath of a collision partner is presented. Transport properties can be computed from collision integrals (n, s) (T), which can be obtained by taking an appropriate integral over the collision energy and carrying out a Boltzmann state average.4, 20, 34, 35 Reversing the usual order of carrying out this integration and state averaging, we can express the collision integral as a Boltzmann state average over state-dependent collision integrals (n,s) ji (T ): 1 (2ji + 1) exp(−εi /kB T )(n,s) (n,s) (T ) = ji (T ). (1) qR j i
In Eq. (1), εi is the energy of the ith rotational level, qR is the rotational partition function, and kB is the Boltzmann constant. The state-dependent collision integrals are the integrals of state-dependent effective cross sections Q(n) ji (E) over the collision energy E: (n,s) ji
kB T 1/2 1 2π μ (kB T )s+2 ∞ × E s+1 exp(−E/kB T )Q(n) ji (E)dE.
1 = 2
(2)
0
Here, μ is the atom-molecule reduced mass. The state-dependent effective cross section in Eq. (2) is a sum over final levels jf of state-to-state effective cross sections: Q(n) ji (E) =
Q(n) ji →jf (E).
(3)
where Rˆ = (θ, φ) is the orientation of the Jacobi vector R. The weighting factors in Eq. (4) for n = 1 and 2 are4, 11, 34 1 (E) = 1 − (E /E)1/2 P1 (cos θ ), 2 (E) =
5 6
− 16 (E /E)2 − 23 (E /E)P2 (cos θ ),
(5) (6)
where E is the relative translational energy after scattering into level jf . It should be noted that for integral cross sections the weighting factor in Eq. (4) is 0 (E) = 1. Details about the evaluation of the Legendre moments in the effective cross sections in Eq. (4) can be found in the literature.36, 37 The effective cross sections in Eq. (4) can be written as the sum of several low-order Legendre moments of the differential cross sections: Q(1) ji −jf =
Q(2) ji −jf =
π (−1)ji −jf A0 − 13 (E /E)1/2 A1 , (2ji + 1)ki2
π (−1)ji −jf 5 1 − 6 (E /E)2 A0 − 6 2 (2ji + 1)ki
(7)
2 (E /E)A2 15
.
(8) In Eqs. (7) and (8), ki is the initial wave vector and Aλ equals Aλ = Z(l1 J l2 J ; ji λ)Z(l1 J l2 J ; jf λ) J J l1 l2 l1 l2
∗ × TjJi l1 ,jf l TjJi l2 ,jf l . 1
2
(9)
In Eq. (9), T designates a T-matrix element, expressed in the space-fixed frame, and the Z coefficients36, 38 are given by 1 1/2
Z(lJ l J ; j λ) = (−1) 2 (λ−l+l )+J +J [λ][l][l ][J ][J ]
×
l l λ 000
l J j , J l λ
(10)
where [x] = 2x + 1, ·· ·· ·· is a 3j symbol, and ·· ·· ·· is a 6j symbol.39 We employed Eqs. (3) and (7)–(10) to compute the statedependent effective cross sections in the scattering calculations described in Sec. IV. These cross sections were used in Eq. (2) to derive state-dependent collision integrals. A Boltzmann state average of these integrals was finally carried out [Eq. (1)] in order to compute collision integrals. We can relate the collision integrals to two transport properties, namely the binary diffusion coefficient Dab and the quantity ηab that determines the curvature of the mixture viscosity as a function of concentration. The relations are20, 34, 35
jf
The state-to-state effective cross sections on the righthand side of Eq. (3) are weighted angle averages of the ji → jf differential cross section: Q(n) ji →jf (E)
=
dσ d
ˆ n (E)d R, ji →jf
(4)
Dab =
3kB T , 16N μ(1,1)
(11)
5kB T , 8(2,2)
(12)
ηab =
where N is the total number density of the gas.
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
III. POTENTIAL ENERGY SURFACE
−4
6
−5 0 −53 −5
120
−474 −4
−56
−41
−5 −60 6 −47 −50 −5 3 −4 4
−41
−4 4
−50
−44
−53
−47
60 −53
−44
−47
φ / degrees
−53 −5 −47 0
−56
−56
180
−60
−47 3 −50 −5
−41
4
−
4
−41
240 56
−4
−53
−44
−50
−47
−47
−50
300
−53
The coordinates used to describe the H2 O–H PES are the same as those used previously in our study of the interaction of the He atom with both H2 O and CH2 .18, 19, 40 As depicted in Fig. 1 of Ref. 18, the body-frame z axis is chosen to lie along the a inertial axis of the molecule, and the molecule lies in the xz plane. The water molecule was taken to be rigid, with bond length and bond angle (r = 1.8361 bohr and θ = 104.69◦ , respectively) corresponding to the vibrationally averaged geometry reported by Mas and Szalewicz.41 The PES for the H2 O–H system was computed with the MOLPRO 2010.1 suite of computer codes,42 by means of restricted coupled cluster calculations with inclusion of single and double excitations, as well as perturbative contributions of connected triple excitations [RCCSD(T)].43, 44 An atomcentered avqz atomic-orbital basis was used,45 with the addition of a set of bond functions located midway between the atom and the molecule.46, 47 A counterpoise correction was used at all geometries to correct for basis set superposition error (BSSE).48 No hint of access to the reactive channel leading to OH + H2 products was found in the calculations; this is due to the fact that the H2 O geometry was held fixed. Since we are interested in collisions at thermal energies and only slightly above, it is reasonable to neglect the reactive channel. The PES was determined on a grid of 20 values of the atom-molecule separation R ranging from 3.0 to 10 bohr in steps of 0.5, with additional points at 11, 12, 13, 15, and 20 bohr. By exploiting symmetry we can limit the calculations to the angular range 0◦ ≤ θ ≤ 90◦ and 0◦ ≤ φ ≤ 180◦ , both varied in steps of 10◦ , for a total of 190 orientations. The total number of nuclear geometries at which the interaction energy was computed was 3800. The computed potential as a function of R and the orientation was fit to the following expansion in spherical harmonics: Vλμ (R)(1 + δμ0 )−1 V (R, θ, φ) =
360
−50
194309-3
0 0
30
60
90 120 θ / degrees
150
180
z H
R y
θ=30o
x H O
φ=180o H
FIG. 1. (Top panel) Dependence of the potential energy (in cm−1 ) on the orientation of the hydrogen atom collision partner with respect to the H2 O molecule for an atom-molecule separation R = 6.5 bohr. (Bottom panel) Arrangement of the atoms at the H2 O–H equilibrium geometry [R = 6.50 bohr, θ = 30◦ , φ = 180◦ ].
λ,μ≥0
× [Yλμ (θ, φ) + (−1)μ Yλ,−μ (θ, φ)]. (13) Symmetry considerations restrict the sums in Eq. (13) to terms with λ + μ even.18 To get an acceptable fit of the PES, we needed to include all terms with λ ≤ 10 and μ ≤ 8, for a total of 33 anisotropic terms. Figure 1 presents a contour plot of the dependence of the potential energy upon the orientation of the H atom collision partner at an atom-molecule separation R = 6.50 bohr, the separation at the global minimum. The anisotropy of the PES at this value of R is seen to be small. This is true for most values of R, except at small atom-molecule separations for which there is significant repulsion between the H atom collision partner and the H atoms on the water molecule. The global minimum of the PES has an energy of 61.0 cm−1 below the H2 O + H asymptote, at a geometry of R = 6.50 bohr, θ = 30◦ , φ = 180◦ . The hydrogen atom collision partner thus lies within the molecular plane at the global minimum and close to one of the H atoms of the water molecule. Because of the small anisotropy of the PES, the angular orientation for the minimum interaction energy varies somewhat
as a function of R. The binding energy for H2 O–H is somewhat larger than for the related H2 O–He system, for which the binding energy is computed to be 35 cm−1 , while the atom-molecule separation R for the latter system is slightly smaller (5.92 bohr).49 In their calculations on H2 O–H, Zhang et al.21 considered that their best estimate for the well depth of this system (53 ± 6 cm−1 ) was intermediate between the values computed in MP4 calculations without and with the counterpoise correction (69.9 and 36.8 cm−1 , respectively). The presently computed well depth of 61 cm−1 for H2 O–H is slightly larger than the estimate of Zhang et al.21 Figure 2 presents a plot of the larger expansion coefficients Vλμ as a function of the atom-molecule separation R in the region of the van der Waals well. The largest coefficient in the region of the well is the isotropic V00 term. We see that for values of R smaller than that of the global minimum some of the anisotropic coefficients become compatible to or larger than the V00 term. In their investigation of the interaction of H2 O with H, Zhang et al.21 derived a spherically averaged, isotropic potential by taking the average of radial cuts along six orientations
194309-4
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
with a small nonlinear polyatomic molecule by averaging cuts along six orthogonal approach directions of the perturber to the molecule, at least for the H2 O–H system. A similar approach has been employed by Stallcop et al.12, 14 for several atom-molecule systems. Jasper and Miller50 have considered methods to spherically average potentials involving larger, floppy molecules. Also plotted in Fig. 3 is the spherical average obtained by Zhang et al.21 by taking the average of MP4 calculations with a BSSE correction along the ±x, ±y, and ±z directions. This potential energy curve has a slightly shallower well and longer equilibrium atom-molecule separation than the present calculation.
200 150
energy / cm
−1
100
V20
V31
50
V11
V22
0 −50 −100 V00
−150 −200 4
5
6
7 R / bohr
8
9
10
FIG. 2. Dependence of the larger expansion coefficients Vλμ [defined in Eq. (13)] on the atom-molecule separation R.
parallel and perpendicular to the C2v axis of the molecule. With the availability of our full PES for this system, we can check the accuracy of this method to generate a spherically averaged H2 O–H potential. Figure 3 presents a plot of the spherical average of our full PES [which equals (4π )−1/2 V00 ]. A plot of the average of the radial cuts of our PES along the ±x, ±y, and ±z directions is indistinguishable from the spherical average on the scale of this figure. At the minimum of the potential, the two curves differ by 0.5 cm−1 , while the well depth of the former equals 48.6 cm−1 . At 4 bohr, the potentials differ by 31 cm−1 , as compared to an interaction energy of 1240 cm−1 . As an indication of the small anisotropy for this system, we note that the well depth for the spherically averaged potential is only slightly smaller than the well depth for the full PES. This comparison of the two ways to obtain a spherical average of the potential suggests that it may be reasonable to compute an isotropic potential for the interaction of an atom
100
energy / cm−1
present work Zhang et al.
50
0
−50 4
5
6
7 R / bohr
8
9
10
FIG. 3. Spherically averaged H2 O–H potentials determined by averaging radial cuts along the ±x, ±y, and ±z directions. Red: present work; blue: MP4 calculation with counterpoise correction by Zhang et al. (Ref. 21).
IV. SCATTERING CALCULATIONS
Close-coupling calculations were carried out with the HIBRIDON suite of programs,51 which has recently been extended to allow the determination of the Q1 and Q2 effective cross sections [Eqs. (7) and (8)]. The cross sections were checked for convergence with respect to the inclusion of a sufficient number of partial waves and energetically closed channels. At the highest energies for which calculations were carried out (4000 cm−1 ), the rotational basis included all levels whose energy was less than 4600 cm−1 , and the scattering calculations included all total angular momenta J ≤ 133 ¯. In previous work,11 we checked that the computed effective cross sections agreed with values obtained by numerical integration of the differential cross sections weighted by the 1 (E) and 2 (E) factors [Eqs. (5) and (6)]. The high-energy tails of the integrals over collision energy [Eq. (2))] were evaluated by exponential extrapolation of the integrand. The sum over rotational levels ji in Eq. (1) included all levels whose energy was less than 3000 cm−1 . For levels whose energy is between 1200 and 3000 cm−1 , the total effective cross section [Eq. (4)] was taken to be equal to that of the highest rotational level of the same ka projection number with energy less than 1200 cm−1 .
V. RESULTS
The rotational levels of an asymmetric top like water are labeled jka kc , where j is the rotational angular momentum, ka is its projection along the principal axis in the prolate limit, and kc is its projection along the principal axis in the oblate limit.52 Each rotational level exists for only one of the nuclear spin modifications (ortho and para). The two nuclear spin modifications cannot be easily interconverted in nonreactive collisions. Figure 4 displays plots of the collision energy dependence of the state-dependent effective cross sections Q(1) ji and [Eq. (3)] for collisions of para-H O in its lowest and a Q(2) 2 ji higher rotational level [000 and 551 , respectively]. The rotational energy of the latter is 758 cm−1 . The elastic contribution to these effective cross sections is also plotted. As we ˜ a)–He ˜ collisions,10, 11 the observed in OH(X)–He and CH2 (X, (1) Qji effective cross sections decrease monotonically with increasing collision energy.
P. J. Dagdigian and M. H. Alexander
total
V
00
term only
16
i
3
/ cm 10
18 17
Ω 11
(d) 551 V00
(b) 1500 K
s
–11 –1
12
V00 term only
3
13
i
elastic
/ cm 10
total
–11 –1
19
14
(1,1) j
V00
20
(a) 300 K
s
(2)
(c) 551
i
elastic
Q j effective cross section / Å2
total
i
(1)
Q j effective cross section / Å2
V00
V00
15
(b) 000
(1,1) j
(a) 000
J. Chem. Phys. 139, 194309 (2013)
Ω
194309-5
15 14
0 2 4 6 8 10 rotational angular momentum j
0 2 4 6 8 10 rotational angular momentum j
total elastic
elastic
collision energy / cm–1
FIG. 4. Dependence on collision energy of the state-dependent effective (1) (2) cross sections Qji [panels (a) and (c)] and Qji [panels (b) and (d)] for collisions of H2 O in the 000 and 551 initial levels with hydrogen atoms. The elastic contribution to the effective cross section is indicated in each panel. In addition, the effective cross section computed with the spherical average of the potential is displayed for comparison.
At low collision energies, elastic collisions make by far the dominant contribution to the Q(1) ji effective cross section, while at higher energies inelastic scattering plays a much larger role. As we have noted previously,10, 11 this contrasts with the dominant role of elastic scattering in total integral cross sections. We showed in our study of CH2 –He transport properties11 that the differing importance of elastic scattering in the total integral and transport effective cross section results from the angular weighting of the differential cross section in the latter [see Eq. (4)]. Because of the weighting factor n (E) [Eqs. (5) and (6)], transport cross sections, particularly for elastic transitions, are much smaller than the corresponding integral cross sections, as discussed previously.11 Also included in the plots in Fig. 4 are the effective cross sections computed with the spherical average of the PES [(4π )−1/2 V00 ]. Of course, only elastic collisions can take place with this potential. We see that the Q(1) ji effective cross sections for this isotropic potential and the full PES are very similar in magnitude for all collision energies. We consider below the effect of the anisotropy of the PES through examination of state-dependent collision integrals [Eq. (1)]. It can be seen in Fig. 4 that the Q(2) ji effective cross section for the 000 level displays a similar collision energy dependence as for Q(1) ji for this level, namely a monotonically decreasing magnitude vs. collision energy. By contrast, the Q(2) ji effective cross section for the 551 level becomes very large and negative at low collision energies. As noted previously for Q(2) ji effective cross sections for higher rotational levels,10, 11 this is a consequence of large values for the factor E /E in the 2 (E) weighting factor in Eq. (6) for superelastic collisions (E > E) at small collision energies. As we have discussed previously,10, 11 this singularity does not cause problems in the evaluation of the Q(2) ji effective cross sections for the higher rotational levels, and hence also the determination
(1,1)
FIG. 5. State-dependent collision integrals ji as a function of the rotational angular momentum j for H2 O–H at temperatures of (a) 300 K and (b) 1500 K. Red and blue symbols denote ortho and para levels, respectively. Levels with different values of the body-frame projection quantum number are denoted by the following symbols: ka = 0, circles; ka = 1, plus signs; ka = 2, diamonds; ka = 3, upward-pointing triangles; ka = 4, downwardpointing triangles; ka = 5, left-pointing triangles; and ka = 6, right-pointing triangles. In each panel, the value of the collision integral computed with the spherical average of the potential is indicated with a dashed line.
of the (2,2) collision integral. As with the Q(1) ji effective cross sections, elastic collisions make a decreasing contribution to Q(2) ji effective cross sections with increasing collision energy. Finally, we note that the values of Q(2) ji computed with the spherical average of the potential are very similar to the values computed with the full PES, except at the lowest collision energies for the higher rotational levels. From the state-dependent effective cross sections we can use Eq. (1) to compute the (1,1) and (2,2) collision integrals. We plot in Fig. 5 the variation with rotational level of the statecollision integrals at two temperatures. The dependent (1,1) ji (2,2) ji collision integrals display a similar dependence upon the rotational quantum numbers and hence are not plotted is plotted are those here. The rotational levels for which (1,1) ji for which collision-energy dependent effective cross sections have been explicitly computed (see Sec. IV). The most probable value of j is ∼3 at 300 K and ∼7 at 1500 K. We see that on the rotational quantum at 300 K the dependence of (1,1) ji numbers is weak. This dependence is stronger at 1500 K. In values are seen in Fig. 5(b) to show a departicular, the (1,1) ji crease with increasing j. Presumably, this reflects the greater anisotropy of the PES at higher energies, because of the repulsion of the hydrogen atom collision partner by the hydrogen atoms on the water molecule. We also plot in Fig. 5 the effective cross sections under the assumption that the interaction is purely spherical (only the V00 term); in this case there is no dependence on rotational level. As shown in Fig. 5, it is easy to use the spherically averaged potential to predict collision integrals. Figure 6 displays the temperature dependence of these quantities. We see that the collision integrals computed with the spherically averaged potential are slightly larger than those computed with the full PES. Specifically, the collision integrals computed with the former are ∼1.4% larger at 300 K and 6% larger at 1500 K. These differences are smaller than what we found ˜ a)–He ˜ systems, likely a reflection of the larger for the CH2 (X, anisotropies for the PESs for these systems.
194309-6
P. J. Dagdigian and M. H. Alexander
J. Chem. Phys. 139, 194309 (2013)
2.5
full PES
−1 4
full PES
Dab / 10 Torr cm s
V00 only
2 −1
2
V00 only 1.5
1
LJ 12–6 0.5
0 200
full PES
FIG. 6. The collision integrals (1,1) and (2,2) as a function of temperature for the H2 O–H system. These quantities have been computed for the full PES and for the spherically averaged potential.
VI. DISCUSSION
We have employed quantum scattering calculations to compute collision integrals for the H2 O–H system, from which transport properties can readily be predicted. As an example we will consider the H2 O–H diffusion coefficient. We can also compare these values with a conventional calculation using a parameterized LJ 12-6 potential. The parameters specifying this potential were estimated through conventional combining rules: To estimate the well depth ε and the length parameter σ for H2 O–H, we took the geometric and arithmetic means, respectively, of the corresponding parameters for the like systems [ε = 572.4 K and σ = 2.605 for H2 O– H2 O and ε = 145 K and σ = 2.050 for H–H].53 Since H2 O possesses a dipole moment, there is an induction contribution to the long-range attractive interaction. Using the dipole moment of water (1.844 D)53 and the polarizability of the H atom (0.663 3 ),54 we follow the procedure of Hirschfelder et al.20 to correct the LJ parameters; this correction is quite small (
