J . Phys. Chem. 1993, 97, 2523-2526
2523
Rotational Relaxation and Transport Properties of Oxygen by Using the Importance Sampling Method Lichang Wang and Gert D. Billing' Department of Chemistry, H . C. 0rsted Institute, University of Copenhagen, DK-2100 Copenhagen, Denmark Received: September 9, 1992; I n Final Form: November 19, 1992
The coefficients of rotational relaxation and transport properties of oxygen are calculated by using the rigidrotor model with two analytical potentials in the temperature range 200-600 K and compared with experimental results. In the molecular dynamics simulations, the importance sampling method, which was used in the calculation of viscosity of nitrogen a t one temperature, is extended here to several properties in a range of temperatures.
1. Introduction
The 0 2 molecule has an open-shell structure, which makes it different from the closed-shell molecules such as N2, H2, etc. From the viewpoint of intermolecular interactions the potential of open-shell molecules depends not only on the intermolecular distance and theorientation of molecules but also on the orientation of the spins of the molecules. Also 0 2 is abundant in the atmosphere and plays an important role in many atmospheric and other chemical reactions. It is the above nature that makes the studies of its properties of great interest from both theoretical and practical aspects. Although the investigation of the properties of 02is interesting and important, the measurements of them are rare, especially of thermal conductivity. In 1973, Hanley and Elyl reviewed the situation and stated that the available data for oxygen viscosity are generally poor. The data for the thermal conductivity are even poorer than those for viscosity. There are few works on the measurement of the transport properties in the past 20 years, which has slowed down the improvement of the experimental data. On the other hand, the development of theoretical studies on the 0 2 in gas phase is also very slow. McCourt et al.2 reviewed the theoretical studies of transport properties recently but did not mention any calculations on 02. The transport properties and rotational relaxation of oxygen are calculated here in the temperature range 200-600 K. As in the previous calculations for hydrogen fluoride and hydrogen chloride gases3 we treat the oxygen molecule as a rigid-rotor and use the first-order classical limit of the Wang Chang-Uhlenbeck theory for which the classical mechanical treatment is satisfactory in the range of temperatures considered. Here we extend the importancesampling method, which weemployed in the molecular dynamics simulations of N2 to calculate vis~osity,~ further. In the present paper we have considered two analytical potential surfaces suggested by van der A ~ o i r d . ~They , ~ are described in section 2. The rigid-rotor model and importance sampling method are reviewed in section 3. In sections 4 and 5 we give the calculated results and the discussion, respectively. The conclusions are given in section 6.
with angular functions
The functions YL.,+,(Q)are spherical harmonics and the first factor in eq 3 is a 3 -jcoefficient. The angular function A and the potential depend on the polar angles of R and the molecular axes denoted by 52, uA,and uB,respectively. In their treatment, a special body-fixed frame with 52 = (@) = (O,O), U A = (OA,O), and W B = (OB,&,) was used. The expansion coefficients in eq 2 can be written as fL,.LB.L(R) = ~ L , + L , , L ~ L , . L ,
R-LA-LB-I
+E,,L,,L(R)
(4)
where the electrostatic coefficients are given by
with the multipole moments QL,,,QL,.
fLA,L,,L is given as
The exchange contributions
(R) = g o L ~ . L ~ . L(1 + y L ~ . L ~x) .L x e ~ p ( - ~ L A h
[email protected] ) (6) where x = ( R - Ro)/Ro. All above constants in eqs 2-6 come from ref 6. For convenience the potential represented by eq 1 is referred to as potential A and the other as potential B in the following text. The complete forms of both potentials A and B include the dispersion
f xL,.L,.L
term
2. Potentials For the interaction between two oxygen molecules we use two different analytical site-site potentials. A simplified form of the potential is suggested by van der AvoirdS from the a b initio calculations and it is formulated as where A = 346.72 X
between the atoms in different molecules. The quadruple moment Q is -0.2714 a ~ . ~ Another is also given by van der Avoird6 from the a b initio calculation but by virtue of a spherical expansion for the orientational dependence and exponential functions for the distance dependence of the expansion coefficients. The expression is as follows
V,, = A exp(-BR) (1) lo5K, B = 0.21 19ao-', and R is the distance 0022-3654 f 93f 2097-2523SO4.00 f 0
Vdisp= -C6f R6
(7)
where C,,= 61.57 au.' 3. Theory
a. Rigid Rotor Model. In the present paper, we calculate the transport properties and the rotational relaxation coefficients in 0 1993 American Chemical Society
2524
Wang and Billing
The Journal of Physical Chemistry, Vol. 97, No. 11, I993
the temperature range 200-600 K. Under these conditions, one can treat the 02 molecule as a rigid diatomic rotor which only has rotational motion. The classical limit of the Wang ChangUhlenbeck theory developed in ref 8 is as before also used here. For the process
where cfl,, = c I n t / kand
a = -3 %
CintT
a(((€, ;)k, - t , h 2 - (€, - 6,)YY’ cos 511)
(21)
b. Important Sampling Method. In order to calculate the transport properties, we have to compute the average collision integrals (9), which can be rewritten as
the average collision integral is written as
( { 1)
= JF(UJ-I exp(-Pv) d(PW
(22)
where 2Th2 F(U,T) = -
IIQ~Q~
where
gr
is the collision cross section gr
dO = 2~Pf)‘(b)bd b
(10)
and P,:’(b) is the transition probability as a function of impact parameter b, t, = E , / k T , y2 = ‘ / 4 m u 2 / k TQA , = &xp(-c,) and At = c k - t I- e,, m is the mass of 02,EJ is the internal energy, u is the relative velocity, and QA the rotational partition function. In order to calculate the properties in question, the following collision integrals are needed
+
(11)
where 6 is defined as the angle by which the direction of the rotational angular momentum vector of one of the colliding molecules is changed due to the collision, 5 is the scattering angle in the relative coordinate system, T~ the reduced velocity after collision, and the average internal energy for a molecule. The rotational and transport coefficients can be obtained by the first-order approximation of the Wang Chang-Uhlenbeck theory and the above collision integrals. (a) The relaxation collision number is given as
and { ]denotes a particular quantity corresponding to the transport properties and the rotational relaxation time. Qaand Qb are rotational partition functions for the diatomic molecules, p the reduced mass for the relative motion. jl,j2, and lare the rotational and the orbital angular momenta, respectively. Formula (9) or (22) is a multidimensional integral and is evaluated by Monte Carlo averaging technique or so-called straightforward sampling method. According to this method, the points at which the function is evaluated are chosen uniformly and no reference is made to the nature of the function. If the function has a large variation, the uncertainty from the Monte Carlo estimate will be large. Consequently, if the function is uniform, then the estimate will be most efficient. For practical purposes, the convergence of the straightforward sampling is too slow; for example, an increase of the sample size N by a factor of 100 results in a decrease of the uncertainty by a factor of only 10. In order to increase the speed of convergence, one can make the variance as small as possible, which is the idea of the importance sampling method. We used the importance sampling method in the calculation of viscosity of nitrogen at only two temperatures and the details were discussed in ref 4. In that article we only tested the method a t two temperatures, but they are calculated separately, i.e. we did not utilize the trajectories from one temperature to another. In this respect we can also continue to improve the approach. Here we use the method and extend it to several properties and temperatures. We choose the importance sampling function p(v> = AUY, as in ref 4, then eq 22 becomes ( ( 1) = s q I - ( v
(Pv)
where 7 is the viscosity and c,,~the internal heat capacity. (b) The cross section for depolarized Rayleigh spectral linewidths
where
in which x changes from 0 to 1 corresponding U E (0,a) and y( ) is the incomplete gamma function. So, let
(c) The shear viscosity 7’
+ 1) dx
0.0863 125 T
((r4sin’ F - ‘/2(At)2 sin2 5 + ‘/3(At)2))
(d) The coefficient of self-diffusion
(e) The thermal conductivity
(18) the integral (24) can be rewritten as
Instead of calculating the trajectories separately at each temperature, we select a part of the trajectories calculated and used at previous temperatures and only increase with a small number of trajectories. The main points for calculating a number of temperatures for a given property are
Transport Properties of Oxygen
The Journal of Physical Chemistry, Vol. 97, No. 11, 1993 2525
TABLE I: Collision Integrals in Units of 10-10 cm3 s-1 T/K eq 11 eq 12 eq 13 eq 14 eq 200 300 400 500 600
0.822 0.721 0.816 0.740 0.812 0.754 0.808 0.765 0.806 0.774
2.24 1.67 2.28 1.77 2.33 I .84 2.37 1.90 2.40 1.94
1.20 1.01 1.19 1.06 1.18 1.10 1.18 1.13 1.18 1.16
1.09 0.746 1.09 0.794 1.09 0.829 1.09 0.858 1.08 0.882
15
0.555" 0.470h 0.525 0.475 0.508 0.479 0.495 0.482 0.485 0.484
" The results from potential A. The results from potential B. (a) Assign a basic temperature To arbitrarily, here we choose 200 K. (b) Calculate at the basic temperature exactly as described in rei 4, and these trajectories compose a trajectory basis. (c) For a given temperature T , we need only to change the output of x in the trajectory basis obtained from step b to xd which is the value of x at the particular temperature. This can be done by changing @U = U/kTo in eq 25 to U / k T , where T , To are the particular and basic temperatures respectively, to get xd. It is not a CPU-costly step since we have Ufrom the trajectory basis. Then replace the xo as xd in our trajectory basis and it can be used in the following calculation. (d) Follow the steps described in ref 4 again. First we compare the x which was chosen randomly with xd. If they are nearly equal, then we use this trajectory in the calculations concerned, otherwise, we need to calculate a new trajectory. This principle is also applied toget different properties in the present calculation. The procedure only differs in step c where we change x because u is different for various properties at the temperature of interest. If we need to calculate a set of temperatures, then follow the steps c and d. Here it should be stressed that the trajectory basis can be enlarged with the need of extending the calculation to get the maximum utilization of trajectories by including the information from that for the temperatures of interest. It is sometimes useful to switch from one trajectory basis set to a new basis set in order to calculate different properties conveniently.
TABLE 11: Rotational Relaxation and Transport Properties as a Function of Temperature T IK
P
UDPRn
Xh
7'
#
200
2.77 2.36 2.85 2.44 2.92 2.48 2.99 2.52 3.03 2.55
86.3 73.1 66.6 60.3 55.8 52.7 48.7 47.4 43.5 43.4
9.22 11.8 13.7 16.7 18.1 21.5 22.4 26.1 26.7 30.6
7.70 10.3 11.4 14.6 14.8 18.8 18.2 22.7 21.6 26.7
0.0487 0.07 12 0.1 1 1 0.150 0.195 0.256 0.304 0.387 0.443 0.542
300 400 500 600
In units of A?. In units of I O ) W m In units of cm? s at 1 atm
'K
I.
In units of IO-h Pass.
d*
.
----V."
200
300
400
500
600
t e m p e r a t u r e (K)
Figure 1. Viscosity as a function of temperature. Dashed and solid lines represent the calculated results from potential A and B respectively. Experimental data from refs 10 and 1 1 are given by 0 and A,respectively. X
4. Results
In the present calculation, we use the formula shown in the above section to get the integral (22). Equation 27 can be evaluated further by , N,
In the first step we need to find the u, which is done by running trajectories at four energies, at each about 250 trajectories are used, and then choose x E [O,l] randomly to calculate the properties following the principle described in the above section. Tables I and I1 give the collision integrals and the properties considered respectively for the two potentials. For potential A , by using about 1800 trajectories at each temperature for a given property, the error bar of the calculation of collision integrals reduces from 14% (running 1000 trajectories) to less than 4.0%. On the contrary, by the straightforward sampling method, if we want to reduce the error bar to the same extend, we must run 16000 trajectories. For potential B we obtained results within 8% by running 1100 trajectories. 5. Discussion
As mentioned in subsection 3b, one of the main purposes of the present work is to find a way to extend the importance sampling method to a broad temperature range and to calculate several
'
5.0 200
I
1
1
1
300
400
500
600
temperature (K) Figure 2. Thermal conductivity as a function of temperature. Same meaning of dashed and solid lines as in Figure I . Experimental data from refs 12-1 5 are given by A, 0 , 0, and 0 , respectively.
properties by using as few trajectories as possible. The principle employed here solved the above puzzle successfully. We cannot add all the trajectories at T I to T2 which we are interested in because the contribution of energy U to collision integrals at different temperatures or the distribution of sampling points with the energy U is different for different temperatures. While calculating the integral at 300 K, we use 1790trajectories in total and here 965 trajectories come from the calculation at 200 K. That means more than 50% trajectories can be used and save CPU timeof about 50% compared with calculating separately at every temperature. The advantages can be exhibited especially when doing the calculation in the case of small temperature differences. We have compared the viscosity, thermal conductivity, and diffusion coefficient with experimental data which are shown in Figures 1, 2, and 3. From the figures we find the comparison of properties for potential A is in poor agreement with experiments.
The Journal of Physical Chemistry, Vol. 97, No. 11, 1993
2526
0.6 I h
20.5 0
v
4
3
0.4
.e
0.3 Q U
c 0.2 .-0
t 5 0.1 n
0.0
'
200
I
1
I
I
300
400
500
600
temperature (K) Figure 3. Diffusion coefficient a s a function of temperature. S a m e meaning of dashed and solid lines as in Figure I . Experimental d a t a
from refs 16 and 17 a r e given by 0 and A, respectively.
But for the calculations with potential B, the agreement is better than for potential A. Another fact is that the agreement with experimental data is better for diffusion coefficient and viscosity than for thermal conductivity. This may be explained by more inaccuracy in the experimental values of thermal conductivity since the recent measurements for thermal conductivity have shown that thevalues decrease for manysystems. Both potentials A and Bare derived from ab initio calculations which considered only the spin-independent situation. That may also be one of the reasons for the lack of agreement with experimental results. In fact, it is difficult to find other potentials in literature since the specificity of 0 2 increases the complexity of the ab initio or semiempirical determination of the potential. From the calculated results, we can also conclude that the rotational relaxation time for 0 2 is less than that for Nza9 6. Conclusion
From the present calculation, we can get the following conclusions: (a) The importance sampling function p(v) = AU' used for the calculation of transport properties can indeed reduce the variance. (b) Combining the principle given here with the general feature of importance sampling technique described previ~usly,~ we have completed the importance sampling method for calculating the transport properties. (c) We have found that the simple form of the potential given by eq 1 is less accurate for calculating transport properties than the more sophisticated choice eq 2. Although the agreement with experimental data was improved, it is still not good. With the same method and type of potentials we have previously been able to obtain in general good agreement with experimental data on N2. Thus for DPR cross section (not measured for 02) we obtained (see ref 8) UDPR = 40 f 2 and 50 f 3 A at 300 K using two different site-site potentials due to Berns and van der AvoirdIs and Hemert and Berns,19 respectively. Using a potential by Ling
Wang and Billing and Rigby,Zo we recently obtained 45 f 2 A2.21 Closer to the experimental value of 34.4 f 0.6 A2 (ref 22) came Turfa and K n a a ~ . ~They 3 obtained 29.4 f 0.5 A2 using 25 000 trajectories and a semiempirical site-site surface. Although the number of trajectories used in our previous works was much smaller, 1-3000, it allowed for sufficient accuracy for distinguishing among the various PES. For the rotational relaxation cross section Nyeland et ale8 obtained 10.3 f 1.0 and 7.9 & 0.8 A2 with the two surfaces mentioned above. Billing and Wang2' obtained 12.5 f 1.0 A2. The experimental data lie in the range 7.6-10.3 A2.24.25 Turfa et a1.26obtained with the same potential mentioned above 8.9 f 0.1 A2and 7.5 f 0.3 A2 including higher order Sonine correction terms. Thus for N2 the agreement with experimental data is much more satisfactory than for 0 2 . Since the surface used for 0 2 is of the same quality (from an ab initio point of view) one might state that either the experimental data for this system are lessaccurateand/or the fact that 0 2is an open shell type molecule makes it necessary to include this explicitely in the calculations.
Acknowledgment. The Danish Research Academy and the Danish Natural Science Research Council are acknowledged for support of this research. References and Notes ( I ) Hanley, H. J. M.; Ely, J. F. J . Phys. Chem. Ref, Dara 1973,2, 735. (2) McCourt, F. R. W.; Beenakker, J. J. M.; Kahler, W. E.; KuSPer, 1. Nonequilibrium phenomena in polyatomic gases; Clarendon Press: Oxford, 1990; Vol. 1 and 2. (3) (a) Nyeland, C.; Poulsen, L. L.; Billing, G . D. J . Phys. Chem. 1984, 88,5858. (b) Wang, L.; Billing, G . D. J . Chem. SOC.,Faraday Trans. 1992, 88, 163. (4) Billing, G . D.; Wang, L. Chem. Phys. Lerr. 1992. 188, 315. (5) Van der Avorid, A. Private communication. (6) Wormer, P. E. S.; Van der Avorid, A. J . Chem. Phys. 1984,81,1929. (7) Rijks, W.;Van Heeringen, M.; Wormer, P. E. S. J. Chem. Phys. 1989, 90, 6501. (8) Nyeland, C.; Poulsen, L. L.; Billing, G . D. J . Phys. Chem. 1984.88, 1216. (9) Nyeland, C.; Billing, G . D. J . Phys. Chem. 1988, 92, 1952. ( I O ) Touloukian, Y. S.,Ed. ThermophysicalproperriesofMatrer; Plenum Press: New York, 1970 and 1973; Vol. 3 and I I . (1 I ) Lavushchev, A. V.; Lyusternik. V. E. Zh. Fir. Khim. 1976,50,3009. (12) Vanicheva, N. A.; Zaitseva, L.S.;Yakush, L. V. Inrh. Fiz. Zh. 1985, 49, 94. (13) Franck, E. U.2.Elekrrochem. 1951, 55,636. (14) Todd, G. W. Proc. R . SOC.London 1909, A83, 19. (15) Geier, H.; Schafer, K. Allgem. Warmerck. 1961, 10, 70. (16) Winn, E. 8. Phys. Reu. 1950, 80, 1024. (17) Winter, E. R. S. Trans. Faraday Soc. 1951, 47, 342. (18) Berns, R. M.; van der Avoird, A. J . Chem. Phys. 1980, 72, 6107. (19) van Hemert, M. C.; Berns, R. M. J . Chem. Phys. 1982, 76, 354. (20) Ling, M. S. H.; Rigby, M. Mol. Phys. 1984, 51, 855. (21) Billing, G . D.; Wang, L. J . Phys. Chem. 1992, 96, 2572. (22) Keijser, R. A. J.; van der Hout, K. D.; de Grott, M.; Knaap, H. F. P. Physica 1974, 75, 5 15. (23) Turfa, A. F.; Knaap, H. F. P. Chem. Phys. 1981, 62, 57. (24) Prangsma, G . J.; Alberga, A. H.; Beenakker, J. J. M. Physica 1973, 64. 278. (25) Carnevale, E. H.; Corey, C.; Larsen, G . J . Chem. Phys. 1967, 47, 2829. (26) Turfa, A. F.; Knaap, H. P. F.; Thijsse, B. J.; Beenakker, J. J. M. Physica A 1982, 112, 18.
