- Email: [email protected]
On the sampling of microcanonical distribution for rotating triatomic molecules
Chemical
Physics ELSEVIER
Chemical Physics 213 (1996) 243-252
On the sampling of microcanonical distribution for rotating triatomic molecules I. Rosenblum, E.I. Dashevskaya, E.E. Nikitin *, I. Oref Department of Chemistry, Technion - Israel Institute of Technology, Haifu 32000, Israel
Received 6 May 1996
Abstract
For a rotating harmonic molecule, we suggest a successive sampling method from a microcanonical ensemble which allows one to specify the energies of any number of normal modes, and the energies of active and adiabatic rotations. The method is checked by comparing the calculated histogram distributions of various dynamical variables with those predicted by the rigid rotor-harmonic oscillator model. We also compared the distributions over rotational energies for a model of free SO 2 molecule, and the mean energies transferred in SO2-Ar collisions as calculated with the sampling of initial conditions by our method and by the Markov chain sampling procedure. The successive sampling method can be recommended for specifying initial conditions in studies of the collisional energy transfer under condition of local excitation.
1. I n t r o d u c t i o n
A sampling of initial conditions is an essential prerequisite for the trajectory calculations of the energy transfer in molecular collisions. It is usually assumed that the phase space of the internal motion of a molecule (vibrations and overall rotation) before a collision is populated uniformly under certain additional constraints such as a given total internal energy E and the total angular momentum j. Because of these constraints, the selection of initial conditions becomes a non-trivial task. A review of different methods of sampling and the description of a new method for microcanonical sampling with preselected values of E and j is given in Ref. [1]. These methods of selection of the initial conditions suffice for studying the problems of collisional energy trans-
* Corresponding author.
fer for which one is confident that the intramolecular energy redistribution has brought the molecule to the statistical distribution by the time when a next collision is going to occur. However, for some studies one would need the sampling which properly accounts for the energy fluctuations within the isolated molecule. For instance, if one or several modes of a molecule are weakly coupled to the rest of the modes, one should take into account the possibility of incomplete equilibrium between the modes. Then, the initial conditions for initiating the trajectory calculation do not correspond to a statistical distribution. A similar problem in the sampling of initial conditions arises when one wants to study the kinetics of intramolecular energy exchange between a specific mode and the rest of the molecule. One example of a mode with weak coupling to the vibrational reservoir of a molecule is the so-called active rotation of the molecule. This mode is of special
0301-0104/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PII S0301-0104(96)00268-6
244
L Rosenblum et al. / Chemical Physics 213 (1996) 243-252
interest for triatomic molecules, since it was found recently that it represents a kind of doorway step in transferring the kinetic energy of the colliding pair into the vibrational energy of the molecule [2]. Finally, even if the distribution over molecular states is indeed microcanonical, one may be interested in studying special kind of collisions in which molecules with locally excited modes participate [3]. The problem of sampling with the specification of the energies of different normal modes of a harmonic molecule has been addressed by Hase and Buckowski [4]. In the present paper we generalize this approach by way of example of a non-linear triatomic molecule to incorporate the adiabatic and active rotations into the sampling. The method is tested directly for a free molecule by comparison with theoretical predictions for a rigid rotorharmonic oscillator (RRHO) model, and indirectly by comparing the mean energy transferred in collisions calculated by the use of this method and that by the method of Schranz et al. [5,6] under the conditions when the specified quantities are the total energy and the total angular momentum. We realize that the assumption about harmonic vibrations of a molecule is a restrictive one, but in this respect we join the others who use a harmonic model in the sampling of the initial conditions when studying the collision-induced energy transfer (see, e.g., the recent papers [7,8]).
2. Distribution over adiabatic rotational energies for a microcanonical ensemble of nearly symmetric top molecules The rotational and vibrational motion of a triatomic molecule are not separable, and therefore, the definition of the rotational energy is not unambiguous. However, as was shown in our previous work [2], the rotational energy, E R, can be reasonably defined as the averaged time-dependent instantaneous rotational energy ER(t). The latter is taken to be identical with the energy of a rotating rigid body of instantaneous configuration possessing angular momentum j and tensor of inertia I [9,10]. The averaging is performed over many vibrational periods and over the representatives of the subensemble with a fixed value of j. For a symmetric top molecule,
which is an acceptable approximation for SO 2 molecule [11,12], E R c a n be represented in the form which is formally identical with the rotational energy of a symmetric top: g R = Beffj 2 -k ( Aeff - Beff) k 2 ,
(l)
where Aeff and Beff are the effective rotational constants, and k is the mean projection of j onto the line connecting the two oxygen atoms. The use of the mean projection k in Eq. (1) rather than the instantaneous projection Jz of j onto this line is related to the fact that for a slightly asymmetric top the projection Jz oscillates with a small amplitude about its mean value. These oscillations are related to the exchange of the energy between rotations about three principal axes of inertia. One can ignore this phenomenon provided that deviations of Jz from k are small. On a much longer time-scale that corresponds to thousands of vibrational periods, the dynamical quantity k suffers a change. This change is related to the intramolecular vibrational-rotational energy redistribution (IVRR) which is due to the coupling between vibrational and rotational motion. The existence of three nonoverlapping time scales, related to the characteristic period of vibration rv~b, period of rotation frot and the time of IVVR, rWR R implies the inequality rvi b << 7"rot << 7"IVRR .
(2)
The first inequality allows one to define effective rotational constants, the second to speak about the rotational energy as a conserved quantity during a time period that is longer than the rotational period but substantially shorter than the time of the vibrational-rotational energy redistribution. For a longer time, a part of the rotational energy which is proportional to k 2 changes while that proportional to j2 remains constant. For this reason, the rotational energy E R is split into the adiabatic part Ej = Befff and the active part E~ = ( A e f f - Beff)k 2. The latter is called the active rotation, and it is usually assumed, in the framework of statistical theory, to constitute a part of the vibrational reservoir of the molecule. Note that the quantity TIVRR introduced above has no direct relation to another important quantity, the time of the intramolecular vibrational redistribution (IVR), rWR" r~v R is expected to be much shorter that TIVRR, implying that IVRR occurs under condition of
245
1. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
statistical distribution of the vibrational energy over the vibrational degrees of freedom. A subensemble of molecules with definite values of E and j can be equally well defined by two other quantities, the total energy E and the energy of adiabatic rotation, Es. Therefore, if the (E, j) subensembles are selected by binning the microcanonical ensemble with respect to the values of the total angular momenta j, one gets a distribution function over the energy of adiabatic rotation, Fe(Ej)dE j. The general shape of this function can be found by studying the RRHO model. In this model, the total internal energy is represented by the sum of the rotational, Ej + E k, and vibrational, Eo, energies, and the statistical weights of these contributions are identified with the statistical weight of the linear (two-dimensional) rigid rotor p~2o)t(Ej),one-dimensional rigid rotor Prot(Ek) (*) and a system of s harmonic oscillators ~'~b -('>(E ~). The distribution function F~RU°(Ei)dEj is obtained by integrating the microcanonical distribution over all the energies except E /
k 2 cannot exceed j2. These two conditions put together read:
F U " ° ( es)dE
where
E k < min{E - E i,
REj},
R = ( aef f -- Beff)//Bef f .
(6) For a slightly non-linear triatomic molecule whose rotation can be approximated by the rotation of a symmetric top, R exceeds unity, and at small Ej the condition (6) differs from the condition E k < E - Ej. Calculation of the integral in Eq. (3) under the condition (6) leads to a lengthy expression. For the particular case s = 3 we get:
P~R"°(Ej)dEj ,
=
E \2,
E \'/2
+ 2(RE, c l a E , "5\ -E']
l-~
4
l
f ° r E J - < I+R --E
16( 5 'J2 ej eor5 "~1,
N
-
E]
is
NE
the
EJt(REj~3/2
E >-
(7)
I+R
normalization
constant,
N
= -~os~R/(I +R). ~( e - ej - E~-E")P~°t(Ej)dEj (2>
= ffE (1)
(s)
× Prot (E~) d E k P~ib(Eo)
dE~,
(3)
where pm,(Ei) (2) is independent of Ej, Prot (l) (Ek) is proportional to 1 / ~ - k , and Pvib (,) (Eo) is proportional to ( E , ) ' - 1 (see, e.g., Ref. [13]). If one integrates the r.h.s, of Eq. (3) over all possible values of E k and Eo compatible with a given value of E - E j one arrives at the following formula:
r[""° (ej) dej ( e - ej) s- ,/2 aE:.
A qualitative difference between distribution functions given by Eq. (5) and Eq. (7) is that the former decreases monotonically with increase of E s, while the latter shows a maximum at Ej = E(12 + 3R - {24R - 12 ) / ( 3 R 2 + 16R + 28). Since the existence of this maximum has quite general origin related to the conservation of the angular momentum, we believe that this maximum will be present also for more realistic anharmonic models of the molecule.
(4) 3. S u c c e s s i v e s a m p l i n g m e t h o d
For a nonlinear triatomic molecule s = 3, and the normalized distribution function (4) becomes
7(
1-
e
(5)
However, the integration over E k in Eq. (3) is limited not simply by the condition of the energy requirement, E - E i, but also by the condition that
In this section we present the 'successive sampiing' approach which is suitable for specifying initially the energies of both active and adiabatic rotations, E k and E:. This method is based on the known analytical probabilities of the energy redistribution over different modes and therefore is exact for separable systems, when the notion of a mode is well-defined. With nonzero coupling between modes, this
246
1. Rosenblum et a l . / Chemical Physics 213 (1996) 2 4 3 - 2 5 2
method is approximate. However, if the coupling is weak in the sense that it does not affect drastically the pattern of the spectrum of a separable system, the successive sampling method can still be used for trajectory calculations. The following presents a description of the successive sampling method. (i) We begin with the distribution function P(E~; E) over the vibrational energy E~ for a rotating system of oscillators (three vibrations and three rotations) for a given total energy E, calculated considering the limitation connected with the conservation of the total angular momentum. This distribution function is calculated by integrating the product of densities of energy levels for oscillators, active rotation and adiabatic rotation:
Pe( Eo) dEo - E2 dE--''~vrj (e-C,)(R/(R N Ek= 0 ×
- E,,- Ek
rE %=0
~
(E-Eo-Ek-Ei)~dEj
105 E ~ / - 16 ~ V 1
E~ dEo E E
+
1))
d Ek ¢E~
(8)
(ii) In sampling the vibrational energies E~,, E~ and Eo3 of the three normal modes we follow Hase and Buchowski [4]. Since in their work rotations were not considered, we used the total vibrational energy Eo instead of the full molecular energy E. In short, the first two energies were sampled successively from the corresponding distributions
The simplest way to do this, as well as to obtain the conjugate velocities qi simultaneously, is [4] to choose them as
qi
=
2 - v,lki
qi :
cos(27rr), sin(2rrr),
(12)
where /z i is the effective mass of the normal mode and r is a random number uniformly distributed in the interval (0,1). (iv) The next step is the transformation of normal modes coordinates {q} and of conjugated velocities {q} to Cartesian Jacobi (i.e. mass-center-separated ' m c s ' ) coordinates {x mcs} and to the corresponding vibrational velocities {v rib mc~}.The procedure for the determination of the transformation matrix 12 is similar to that described elsewhere (see, for example, Ref. [14]); still, it contains some distinctions connected with the introduction of instantaneous, i.e. non-equilibrium, Eckart conditions. For this reason we would like to dwell on it. L e t cifl A be the coefficients of the matrix 12, which, transforms the normal modes coordinates and momenta into 3 n - 3 Cartesian ' m c s ' coordinates and velocities:
X i ~cs =
3n-6 X 0 ; cs "~ E h=l
Cil3hqA,
(13)
3n-6
Pe,(Eo,)dE~ = 2 ( E ~ - E v , ) d E ~ J E ~ ,
(9)
Pe,.e,.,(Eo~)dEo = d E o J ( E ~ - E v , ) .
(10)
When E~, and Eo2 are fixed, Ev3 is determined as
Eo~ = Eo - Eo, - Eo . (iii) Now, the normal coordinates qi should be sampled for individual oscillators, possessing force constants k~, so as to satisfy their theoretical distributions for harmonic oscillators:
2dqi PE,,(qi) dqi = 7r¢(2Ev+/ki)_q 2
(11)
viTS=
E
ci~xdlA,
(14)
A=l
where xg~ s are the equilibrium ' m c s ' coordinates, q~ and 0h are the chosen coordinates of the normal modes and their derivatives, the subscript i = 1...n - 1, corresponds to the number of Jacobian ' m c s ' particles (see Ref. [14]), and the subscript /3 = x, y, z to the number of Cartesian axes. The summation is accomplished over all 3n - 6 normal modes A. These 3 ( n - 1 ) ( 3 n - 6) coefficients ci~ A should be found from three types of conditions. First, since q;~ are the normal modes, the inverse transformation from {x mcs} to {q} must put the potential energy to
247
1. Rosenblum et al./ Chemical Physics 213 (1996) 243-252
n - 61 % 2q ,2. This requirement the form U = ( / z 0 ,/ 2) x E 3,= gives (3n - 6)(3n - 7 ) / 2 equations: E
iak[3
"~ mcs-"'~. . . . OXia
Xk~ ] q=O
CiaACkfl~t = O,
A,/x= 1..... 3n-6,
A4:/.t. (15)
Second, the condition that the vibrational kinetic energy T v m u s t a s s u m e the f o r m Tv = ( 3n-6 /z0/2)Y"A= 1 (I~A)2 leads to the requirement of orthogonalization of these coefficients:
El.ZiCi[3ACifl~ = ]£Ot~AV,, ifl
/~, ].L = 1 . . . . . 3n - 6,
(16) w h e r e / z 0 is the normalization constant of the dimension of mass, and the summation is accomplished over all i = 1. . . . . n - 1 and /3 = x, y, z. The number of equations is equal to the number of all combinations of normal modes, i.e. to (3n - 6X3n - 5 ) / 2 . The third type of conditions arises from the requirement to separate the coupled rotations and vibrations as much as possible, i.e. to minimize the
25
I
I
t
L
t
I
vibrational angular m o m e n t u m L = E i r mcS × /ziv mcS. Since the vibrational angular momentum depends on normal modes, it is impossible make them vanish for all possible configurations; however, it is possible to do so for a certain configuration. Usually it is done for equilibrium configuration [14]; the corresponding conditions we call the equilibrium Eckart conditions. However, in our case it is more convenient to make the vibrational angular momentum equal to zero for the current state for which the normal modes {q} were sampled. This leads to the so-called instantaneous Eckart conditions providing the general 'collinearity' of the velocity and position vectors: n~l
[ o mcs omcs ]'Zi Xifl CiyA -- Xiy CiflA
i
+
3n-6 ] E (Ci~l~CiyA -- Ci~/xCi~X)q~ = O, #=1
/3,y=x,
y, z; A = 1 . . . . . 3 n -
6.
(17)
These 3(3n - 6) conditions complete the set of 3(n 1)(3n - 6) equations needed for the determination -
I
I
I
20
15 IX. l0
0
-o.1 0.0 O.l -o.1 0.o 0.1 -O.l 0.0 o.1 Symmetric stretch, ~, Bending, ~2, Asymmetric stretch, ~3, reduced units reduced units reduced units Fig. 1. Distribution of normal mode coordinates for the bending mode (the middle column) and two stretching modes (two side columns). Full lines corresponds to theoretical RRHO model, histograms to the simulation by the successive sampling method.
248
I. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
0.20 0.3
> 0.2 0.10
0.1
0.0 0.0
2.0
4.0 Ev),
6.0
Ev 2 or Ev3,
8.0
000
10.0
0.0
2.0
4.0
kcal/mol
6.0
Ev,
8.0
10.0
kcal/mol
Fig. 2. Distribution of normal mode energies calculated by the successive sampling method (histogram) and by the analytical formula for RRHO model.
Fig. 3. Distribution of the vibrational energy E,., calculated the successive sampling method. Full lines corresponds to theoretical RRHO model.
of the coefficients ci~ ,. Thus, Eqs. (15)-(17) produce all the necessary coefficients ci~ a, and, therefore, provide the transformation of normal modes and velocities to the Jacobi coordinates and velocities by the use of Eqs.
(13) and (14). The further transformation to the Cartesian coordinates {x n} and vibrational velocities {Gvib } of individual atoms is trivial, if the condition of the mass center being located in the origin is imposed. The arbitrariness which is still left is elimi-
I
I
I
I
I
I
I
I
I
I
I
I
h
I
1
I
E
q'--
0.4
b 0.3
0.2
0.1
0.0
0
1
2
3
4
5
6
Ej, kcal/mol
7
8
9
10
1
2
3
4
5
6
7
8
9
10
Ej, kcal/mol
Fig. 4. Distribution of adiabatic rotational energy for a harmonic model of SO 2 molecule for the total energy of 10 kcal/mol. (a) Markov chain sampling with the weighting function from Eq. (26); (b) result of the successive sampling method. Full curve corresponds to Eq. (7).
I. Rosenblum et al./ ChemicalPhysics213 (1996) 243-252 nated by three consecutive rotations of the whole molecule through three random Euler angles. (v) Afterwards, two rotational energies Ej, and Ej~ are sampled successively according to the probability densities
active parts of rotational energy are equal, correspondingly, to
e, = 8of, J2= 8eff(J 1
R + 1Ej~'
=E,, +Ejc+
dEj, Pe'e"(g")dEi"= 2 t E j . ( E - E È )
'
(20)
=eR - ej= (eA + G +gc)
(18)
-
dG Pe'e,'eJ~(Ejc)dEjc= 7r E~jcI E - E
249
- Ej - E j ~
E1 + E ~ c + R + I E J A
=R+IEJ.
(vi) The above set of rotational energies determine the components of the angular momentum about the principal axes of inertia:
(19)
J,~= 21/~I,,, When Ej. and Ejc are fixed, the third component of rotational energy, E j , is determined from the equation Ei~ = E - E v - Ej, - E j . The adiabatic and the
I
I
I
I
I
I
I
I
(22)
where m = A,B,C, and I m are the instantaneous values of the principal momenta of inertia.
I
]
I
I
I
E
I
~
I
I
I
0.04
0.03
0.02
0.01
0.00
010
20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 7 0 8 0 Ej, kcal/mol E), kcal/mol Fig. 5. The same as in Fig. 5,
(21)
but for a total energy of 100 kcal/mol.
90100
250
1. Rosenblum et a L / Chemical Physics 213 (1996) 2 4 3 - 2 5 2
(vii) With the above values of angular momentum components, the projections of the angular velocities onto the principal axes of inertia are calculated: (23)
to,~ = J , , / I m .
The directions of these axes in the space-fixed frame are found by the diagonalization of the instantaneous tensor of inertia. Thus, the angular velocity vector to was completely determined. The vector of angular velocity obtained in this way is used to calculate the rotational contribution to the linear velocity of an atom n: rot :
Pn
(24)
¢0 X r n ,
where G is the radius-vector of the nth atom from the origin (such choice renders the zero velocity of the mass center). The resulting set of rotational velocities {vr°t} makes the total rotational energy E r equal to E - E v. Besides, the instantaneous Eckart conditions ensure that the so-called vibrational velocities {Vnrib} are collinear to the radius-vectors r, and, therefore, orthogonal to the rotational velocities {v~°t} determined by Eq. (24). Consequently, the total velocity of each atom (defined as v i = v ivib + vfOt) yields the required value of the kinetic energy T = E - U, U being the molecular potential energy: m i v 2 : ]~ m i ( v / r i b ) 2
T=
- 2
2
= E~ + r v = E - r:.
mi(t)r°t) 2
+E
2 (25)
A comparison of distributions obtained by the successive sampling method with predictions of RRHO model are shown in Figs. 1-3. Since the agreement is very good we conclude that this method of sampling can be recommended for selection of initial conditions from the microcanonical distribution.
4. C o m p a r i s o n o f the successive s a m p l i n g m e t h o d with the M a r k o v chain s i m u l a t i o n for equiprobable distribution o f 'vibrational' energy over vibrations and active rotations
In this section we present a comparison of results obtained with the help of successive sampling method and the Markov chain simulation. We emphasize that
the latter, as described in Refs. [5,6], generates the microcanonical sampling for prescribed values of the total energy E and the total angular momentum j. Therefore, the comparison is carried on under these conditions. On the other hand, we compared both 'collisionless' and 'collision-induced' quantities: the histogram distributions for rotational adiabatic energies for a free SO 2 molecule and the mean energy transferred in the collisions of SO 2 molecule with Ar atoms.
4.1. Collisionless events
In Figs. 4 and 5 we compare the distributions of the adiabatic energy E~ for a model SO 2 molecule for total energies E of 10 and I00 kcal/mol. The 'effective' vibrational energy which is defined as the difference between the total energy and the energy of adiabatic rotation [2,11] actually includes the energy of vibrational modes and the energy of the active rotation. The distributions shown refer to the equiprobable distribution of the 'effective' vibrational energy, E,.ef f = E - E i, over the above modes. In these calculations, we use the value of B e f f = 0 . 0 0 9 4 ( g / t o o l - Az) - I . With this value of Beef, the maximum in the theoretical RRHO curve corresponds to R = 5.336. It should be noted here that the aforementioned value of Beff is slightly different from that used by Troe et al. [11,12] which was found from the condition Beff=E/J . 2. . . . B e f f = 0 . 0 1 ( m o l / g . A2). Our way of determination of Beff was as follows. We defined the instantaneous Beff as (j2/l B +JZ/1c)/(J2 + J ~ ) , where JB and arc are the projections of the total angular momentum on the instantaneous axes of inertia and I~ and I c are the corresponding principal moments of inertia. This value oscillates slightly due to the non-rigidity of the molecule. We used its mean value, obtained by the averaging over approximately 100 vibrations; it is close to the equilibrium value Beq = ( l ~ q + I~q ) / 4 I ~ q 1~q. The panels a in Figs. 4 and 5 correspond to the Markov chain simulation [1], where the sampling of coordinates was accomplished with the use of the weighting function W ( q ) ot [ E -
U( q)] ~/2+W2.
(26)
1. Rosenblum et al./ Chemical Physics 213 (1996) 243-252
251
Table 1 Comparison of the transfer of different components of energy in collisions of SO 2 molecule with A r for two sampling methods. Total molecular energies 10, 50 and 100 k c a l / m o l , heat bath temperature 300 K, harmonic intramolecular potential. All values are in k c a l / m o l Energy transferred
E = 10 k c a l / m o l
E = 50 k c a l / m o l
E = 100 k c a l / m o l
SS method
SNF method
SS method
SNF method
SS method
SNF method
- 0.382
- 0.390
- 1.432
- 1.507
- 2.633
- 2.484
(AEj)
-0.139
-0.142
-0.417
-0.467
-0.603
-0.673
(AE,,eff)
-0.243
-0.248
- 1.015
- 1.039
-2.031
- 1.811
( A E R)
- 0.176
- 0.177
- 0.570
- 0.571
- 0.843
- 0.834
( A E, ) (AEk)
- 0.206 -0.038
- 0.213 -0.035
- 0.862 -0.153
- 0.936 -0.104
- 1.789 -0.240
- 1.650 -0.210
This weighting function is appropriate for rotating nonlinear molecule with s vibrational degrees of freedom; in our case U(q) represents a harmonic potential over the set of q coordinates which describe three vibrational modes, s = 3. It is seen that the histogram obtained is well reproduced by the theoretical curve 1. The panels b in Figs. 4 and 5 correspond to the successive sampling method. Once again, we see a good agreement of the histogram distribution with theoretical predictions for RRHO model. One also sees that both methods yield very close results.
4.2. Collision-induced events We have calculated, by the standard trajectory method as implemented in our previous work [2], the mean energy transferred in collisions of SO 2 molecule with heat-bath Ar atoms using two different procedures for the sampling of initial conditions. The model used in these calculations, is described in detail in [2]. Briefly, the potential of the free SO 2 molecule was constructed as a superposition of atom-atom pair harmonic potentials. The interaction between SO 2 and Ar was taken as a sum pairwise Lennard-Jones potentials between Ar and the atoms of the molecule (the so-called SLJ model [ 11 ]). The
We note in passing that the histogram distribution over Ej presented in Fig. 2 of Ref. [1 l] does not show a well-developed m a x i m u m and can be approximately described by the distribution (5) over the whole range of energies. As indicated by Schranz [15], this m a y be due to poor statistics.
relative velocity of collider was sampled from the three-dimensional Maxwell flux distribution [ 16]. The trajectories started from a mass center separation of 5 A, and were calculated during approximately 25 ps. The average transfer of the total energy, ( A E ) , of the adiabatic rotational energy, (AEj), and of the 'effective' vibrational energy, ( A E v . e f f ) = ( A E ) ( A E j ) , in these collisions was calculated. Also, using the definition of the rotational energy [9,10], we determined the transfer of this total rotational energy, (AER), and of its 'active' part, (AEk). These values were compared to those obtained with the Markov chain sampling procedure [5,6]. The results for molecular energies 10, 50 and 100 kcal/mol and translational temperature T = 300 K are presented in Table 1. The consistency of the quantities, obtained by the new, successive sampling, method ('SS'), with the results of the traditional method of Schranz et al. ('SNF') [1] confirms the applicability of our method.
5. Conclusions
In this paper we have presented a method of sampling of initial conditions for rotating triatomic molecules, which allows one to specify the energies of normal modes, and also the adiabatic and active parts of the rotational energy (the so-called successive sampling method). Within this method, the distribution of over vibrational phases and the Euler
252
1. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
angles that describe the rotation is equiprobable. We checked this sampling against exact analytical results for the rigid rotor-harmonic oscillator model. Also, we have compared some results of this sampling with the appropriate counterparts derived by using the Markov chain sampling procedure, when the only specified quantities are the total energy and the total angular momentum. The good agreement, found in both kinds of comparison, indicates that the successive sampling procedure can be used for a more detailed study of intra- and intermolecular dynamics, e.g. collisionless intramolecular energy redistribution and collision-induced energy transfer under condition of non-uniform distribution of the vibrational energy over molecular modes.
Acknowledgements We acknowledge very fruitful discussions with Professor J. Troe and Dr. H.W. Schranz. This research was supported by the Fund for the Promotion o f R e s e a r c h at t h e T e c h n i o n .
References [1] G. Nyman, S. Nordholm and W. Schranz, J. Chem. Phys. 93 (1990) 6767. [2] I. Koifman, E.I. Dashevskaya, E.E. Nikitin and J. Troe, J. Phys. Chem. 99 (1995) 15348. [3] I. Oref, Chem. Phys. 187 (1994) 163. [4] L. Hase and D.G. Buckowski, Chem. Phys. Lett. 74 (1980) 284. [5] H.W. Schranz, S. Nordholm and F. Freaser, Chem. Phys. 108 (1986) 69. [6] H.W. Schranz, S. Nordholm and G. Nyman, J. Chem. Phys. 94 (1991) 1487. [7] T. Lenzer, K. Luther, J. Troe, R.G. Gilbert and K.F. Lim, J. Chem. Phys. 103 (1995) 626. [8] V. Bernshtein and I. Oref, J. Chem. Phys. 104 (1996) 1958. [9] J. Jellinek and D.H. Li, Phys. Rev. Lett. 62 (1989) 241. [10] J. Jellinek and D.H. Li, Chem. Phys. Lett. 169 (1990) 380. [11] H. Hippler, H.W. Schranz and J. Troe, J. Phys. Chem. 90 (1986) 6158. [12] H.W. Schranz and J, Troe, J. Phys. Chem. 90 (1986) 6168. [13] E.E. Nikitin, Theory of elementary atomic and molecular processes in gases (Clarendon Press, Oxford, 1974). [14] E.B. Wilson, J.C. Decius and P.C. Cross, Molecular vibrations (McGraw-Hill, New York, 1955). [15] H.W. Schranz, private communication. [16] K.F. Lim and R.G. Gilbert, J. Phys. Chem. 94 (1990) 72.
Physics ELSEVIER
Chemical Physics 213 (1996) 243-252
On the sampling of microcanonical distribution for rotating triatomic molecules I. Rosenblum, E.I. Dashevskaya, E.E. Nikitin *, I. Oref Department of Chemistry, Technion - Israel Institute of Technology, Haifu 32000, Israel
Received 6 May 1996
Abstract
For a rotating harmonic molecule, we suggest a successive sampling method from a microcanonical ensemble which allows one to specify the energies of any number of normal modes, and the energies of active and adiabatic rotations. The method is checked by comparing the calculated histogram distributions of various dynamical variables with those predicted by the rigid rotor-harmonic oscillator model. We also compared the distributions over rotational energies for a model of free SO 2 molecule, and the mean energies transferred in SO2-Ar collisions as calculated with the sampling of initial conditions by our method and by the Markov chain sampling procedure. The successive sampling method can be recommended for specifying initial conditions in studies of the collisional energy transfer under condition of local excitation.
1. I n t r o d u c t i o n
A sampling of initial conditions is an essential prerequisite for the trajectory calculations of the energy transfer in molecular collisions. It is usually assumed that the phase space of the internal motion of a molecule (vibrations and overall rotation) before a collision is populated uniformly under certain additional constraints such as a given total internal energy E and the total angular momentum j. Because of these constraints, the selection of initial conditions becomes a non-trivial task. A review of different methods of sampling and the description of a new method for microcanonical sampling with preselected values of E and j is given in Ref. [1]. These methods of selection of the initial conditions suffice for studying the problems of collisional energy trans-
* Corresponding author.
fer for which one is confident that the intramolecular energy redistribution has brought the molecule to the statistical distribution by the time when a next collision is going to occur. However, for some studies one would need the sampling which properly accounts for the energy fluctuations within the isolated molecule. For instance, if one or several modes of a molecule are weakly coupled to the rest of the modes, one should take into account the possibility of incomplete equilibrium between the modes. Then, the initial conditions for initiating the trajectory calculation do not correspond to a statistical distribution. A similar problem in the sampling of initial conditions arises when one wants to study the kinetics of intramolecular energy exchange between a specific mode and the rest of the molecule. One example of a mode with weak coupling to the vibrational reservoir of a molecule is the so-called active rotation of the molecule. This mode is of special
0301-0104/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PII S0301-0104(96)00268-6
244
L Rosenblum et al. / Chemical Physics 213 (1996) 243-252
interest for triatomic molecules, since it was found recently that it represents a kind of doorway step in transferring the kinetic energy of the colliding pair into the vibrational energy of the molecule [2]. Finally, even if the distribution over molecular states is indeed microcanonical, one may be interested in studying special kind of collisions in which molecules with locally excited modes participate [3]. The problem of sampling with the specification of the energies of different normal modes of a harmonic molecule has been addressed by Hase and Buckowski [4]. In the present paper we generalize this approach by way of example of a non-linear triatomic molecule to incorporate the adiabatic and active rotations into the sampling. The method is tested directly for a free molecule by comparison with theoretical predictions for a rigid rotorharmonic oscillator (RRHO) model, and indirectly by comparing the mean energy transferred in collisions calculated by the use of this method and that by the method of Schranz et al. [5,6] under the conditions when the specified quantities are the total energy and the total angular momentum. We realize that the assumption about harmonic vibrations of a molecule is a restrictive one, but in this respect we join the others who use a harmonic model in the sampling of the initial conditions when studying the collision-induced energy transfer (see, e.g., the recent papers [7,8]).
2. Distribution over adiabatic rotational energies for a microcanonical ensemble of nearly symmetric top molecules The rotational and vibrational motion of a triatomic molecule are not separable, and therefore, the definition of the rotational energy is not unambiguous. However, as was shown in our previous work [2], the rotational energy, E R, can be reasonably defined as the averaged time-dependent instantaneous rotational energy ER(t). The latter is taken to be identical with the energy of a rotating rigid body of instantaneous configuration possessing angular momentum j and tensor of inertia I [9,10]. The averaging is performed over many vibrational periods and over the representatives of the subensemble with a fixed value of j. For a symmetric top molecule,
which is an acceptable approximation for SO 2 molecule [11,12], E R c a n be represented in the form which is formally identical with the rotational energy of a symmetric top: g R = Beffj 2 -k ( Aeff - Beff) k 2 ,
(l)
where Aeff and Beff are the effective rotational constants, and k is the mean projection of j onto the line connecting the two oxygen atoms. The use of the mean projection k in Eq. (1) rather than the instantaneous projection Jz of j onto this line is related to the fact that for a slightly asymmetric top the projection Jz oscillates with a small amplitude about its mean value. These oscillations are related to the exchange of the energy between rotations about three principal axes of inertia. One can ignore this phenomenon provided that deviations of Jz from k are small. On a much longer time-scale that corresponds to thousands of vibrational periods, the dynamical quantity k suffers a change. This change is related to the intramolecular vibrational-rotational energy redistribution (IVRR) which is due to the coupling between vibrational and rotational motion. The existence of three nonoverlapping time scales, related to the characteristic period of vibration rv~b, period of rotation frot and the time of IVVR, rWR R implies the inequality rvi b << 7"rot << 7"IVRR .
(2)
The first inequality allows one to define effective rotational constants, the second to speak about the rotational energy as a conserved quantity during a time period that is longer than the rotational period but substantially shorter than the time of the vibrational-rotational energy redistribution. For a longer time, a part of the rotational energy which is proportional to k 2 changes while that proportional to j2 remains constant. For this reason, the rotational energy E R is split into the adiabatic part Ej = Befff and the active part E~ = ( A e f f - Beff)k 2. The latter is called the active rotation, and it is usually assumed, in the framework of statistical theory, to constitute a part of the vibrational reservoir of the molecule. Note that the quantity TIVRR introduced above has no direct relation to another important quantity, the time of the intramolecular vibrational redistribution (IVR), rWR" r~v R is expected to be much shorter that TIVRR, implying that IVRR occurs under condition of
245
1. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
statistical distribution of the vibrational energy over the vibrational degrees of freedom. A subensemble of molecules with definite values of E and j can be equally well defined by two other quantities, the total energy E and the energy of adiabatic rotation, Es. Therefore, if the (E, j) subensembles are selected by binning the microcanonical ensemble with respect to the values of the total angular momenta j, one gets a distribution function over the energy of adiabatic rotation, Fe(Ej)dE j. The general shape of this function can be found by studying the RRHO model. In this model, the total internal energy is represented by the sum of the rotational, Ej + E k, and vibrational, Eo, energies, and the statistical weights of these contributions are identified with the statistical weight of the linear (two-dimensional) rigid rotor p~2o)t(Ej),one-dimensional rigid rotor Prot(Ek) (*) and a system of s harmonic oscillators ~'~b -('>(E ~). The distribution function F~RU°(Ei)dEj is obtained by integrating the microcanonical distribution over all the energies except E /
k 2 cannot exceed j2. These two conditions put together read:
F U " ° ( es)dE
where
E k < min{E - E i,
REj},
R = ( aef f -- Beff)//Bef f .
(6) For a slightly non-linear triatomic molecule whose rotation can be approximated by the rotation of a symmetric top, R exceeds unity, and at small Ej the condition (6) differs from the condition E k < E - Ej. Calculation of the integral in Eq. (3) under the condition (6) leads to a lengthy expression. For the particular case s = 3 we get:
P~R"°(Ej)dEj ,
=
E \2,
E \'/2
+ 2(RE, c l a E , "5\ -E']
l-~
4
l
f ° r E J - < I+R --E
16( 5 'J2 ej eor5 "~1,
N
-
E]
is
NE
the
EJt(REj~3/2
E >-
(7)
I+R
normalization
constant,
N
= -~os~R/(I +R). ~( e - ej - E~-E")P~°t(Ej)dEj (2>
= ffE (1)
(s)
× Prot (E~) d E k P~ib(Eo)
dE~,
(3)
where pm,(Ei) (2) is independent of Ej, Prot (l) (Ek) is proportional to 1 / ~ - k , and Pvib (,) (Eo) is proportional to ( E , ) ' - 1 (see, e.g., Ref. [13]). If one integrates the r.h.s, of Eq. (3) over all possible values of E k and Eo compatible with a given value of E - E j one arrives at the following formula:
r[""° (ej) dej ( e - ej) s- ,/2 aE:.
A qualitative difference between distribution functions given by Eq. (5) and Eq. (7) is that the former decreases monotonically with increase of E s, while the latter shows a maximum at Ej = E(12 + 3R - {24R - 12 ) / ( 3 R 2 + 16R + 28). Since the existence of this maximum has quite general origin related to the conservation of the angular momentum, we believe that this maximum will be present also for more realistic anharmonic models of the molecule.
(4) 3. S u c c e s s i v e s a m p l i n g m e t h o d
For a nonlinear triatomic molecule s = 3, and the normalized distribution function (4) becomes
7(
1-
e
(5)
However, the integration over E k in Eq. (3) is limited not simply by the condition of the energy requirement, E - E i, but also by the condition that
In this section we present the 'successive sampiing' approach which is suitable for specifying initially the energies of both active and adiabatic rotations, E k and E:. This method is based on the known analytical probabilities of the energy redistribution over different modes and therefore is exact for separable systems, when the notion of a mode is well-defined. With nonzero coupling between modes, this
246
1. Rosenblum et a l . / Chemical Physics 213 (1996) 2 4 3 - 2 5 2
method is approximate. However, if the coupling is weak in the sense that it does not affect drastically the pattern of the spectrum of a separable system, the successive sampling method can still be used for trajectory calculations. The following presents a description of the successive sampling method. (i) We begin with the distribution function P(E~; E) over the vibrational energy E~ for a rotating system of oscillators (three vibrations and three rotations) for a given total energy E, calculated considering the limitation connected with the conservation of the total angular momentum. This distribution function is calculated by integrating the product of densities of energy levels for oscillators, active rotation and adiabatic rotation:
Pe( Eo) dEo - E2 dE--''~vrj (e-C,)(R/(R N Ek= 0 ×
- E,,- Ek
rE %=0
~
(E-Eo-Ek-Ei)~dEj
105 E ~ / - 16 ~ V 1
E~ dEo E E
+
1))
d Ek ¢E~
(8)
(ii) In sampling the vibrational energies E~,, E~ and Eo3 of the three normal modes we follow Hase and Buchowski [4]. Since in their work rotations were not considered, we used the total vibrational energy Eo instead of the full molecular energy E. In short, the first two energies were sampled successively from the corresponding distributions
The simplest way to do this, as well as to obtain the conjugate velocities qi simultaneously, is [4] to choose them as
qi
=
2 - v,lki
qi :
cos(27rr), sin(2rrr),
(12)
where /z i is the effective mass of the normal mode and r is a random number uniformly distributed in the interval (0,1). (iv) The next step is the transformation of normal modes coordinates {q} and of conjugated velocities {q} to Cartesian Jacobi (i.e. mass-center-separated ' m c s ' ) coordinates {x mcs} and to the corresponding vibrational velocities {v rib mc~}.The procedure for the determination of the transformation matrix 12 is similar to that described elsewhere (see, for example, Ref. [14]); still, it contains some distinctions connected with the introduction of instantaneous, i.e. non-equilibrium, Eckart conditions. For this reason we would like to dwell on it. L e t cifl A be the coefficients of the matrix 12, which, transforms the normal modes coordinates and momenta into 3 n - 3 Cartesian ' m c s ' coordinates and velocities:
X i ~cs =
3n-6 X 0 ; cs "~ E h=l
Cil3hqA,
(13)
3n-6
Pe,(Eo,)dE~ = 2 ( E ~ - E v , ) d E ~ J E ~ ,
(9)
Pe,.e,.,(Eo~)dEo = d E o J ( E ~ - E v , ) .
(10)
When E~, and Eo2 are fixed, Ev3 is determined as
Eo~ = Eo - Eo, - Eo . (iii) Now, the normal coordinates qi should be sampled for individual oscillators, possessing force constants k~, so as to satisfy their theoretical distributions for harmonic oscillators:
2dqi PE,,(qi) dqi = 7r¢(2Ev+/ki)_q 2
(11)
viTS=
E
ci~xdlA,
(14)
A=l
where xg~ s are the equilibrium ' m c s ' coordinates, q~ and 0h are the chosen coordinates of the normal modes and their derivatives, the subscript i = 1...n - 1, corresponds to the number of Jacobian ' m c s ' particles (see Ref. [14]), and the subscript /3 = x, y, z to the number of Cartesian axes. The summation is accomplished over all 3n - 6 normal modes A. These 3 ( n - 1 ) ( 3 n - 6) coefficients ci~ A should be found from three types of conditions. First, since q;~ are the normal modes, the inverse transformation from {x mcs} to {q} must put the potential energy to
247
1. Rosenblum et al./ Chemical Physics 213 (1996) 243-252
n - 61 % 2q ,2. This requirement the form U = ( / z 0 ,/ 2) x E 3,= gives (3n - 6)(3n - 7 ) / 2 equations: E
iak[3
"~ mcs-"'~. . . . OXia
Xk~ ] q=O
CiaACkfl~t = O,
A,/x= 1..... 3n-6,
A4:/.t. (15)
Second, the condition that the vibrational kinetic energy T v m u s t a s s u m e the f o r m Tv = ( 3n-6 /z0/2)Y"A= 1 (I~A)2 leads to the requirement of orthogonalization of these coefficients:
El.ZiCi[3ACifl~ = ]£Ot~AV,, ifl
/~, ].L = 1 . . . . . 3n - 6,
(16) w h e r e / z 0 is the normalization constant of the dimension of mass, and the summation is accomplished over all i = 1. . . . . n - 1 and /3 = x, y, z. The number of equations is equal to the number of all combinations of normal modes, i.e. to (3n - 6X3n - 5 ) / 2 . The third type of conditions arises from the requirement to separate the coupled rotations and vibrations as much as possible, i.e. to minimize the
25
I
I
t
L
t
I
vibrational angular m o m e n t u m L = E i r mcS × /ziv mcS. Since the vibrational angular momentum depends on normal modes, it is impossible make them vanish for all possible configurations; however, it is possible to do so for a certain configuration. Usually it is done for equilibrium configuration [14]; the corresponding conditions we call the equilibrium Eckart conditions. However, in our case it is more convenient to make the vibrational angular momentum equal to zero for the current state for which the normal modes {q} were sampled. This leads to the so-called instantaneous Eckart conditions providing the general 'collinearity' of the velocity and position vectors: n~l
[ o mcs omcs ]'Zi Xifl CiyA -- Xiy CiflA
i
+
3n-6 ] E (Ci~l~CiyA -- Ci~/xCi~X)q~ = O, #=1
/3,y=x,
y, z; A = 1 . . . . . 3 n -
6.
(17)
These 3(3n - 6) conditions complete the set of 3(n 1)(3n - 6) equations needed for the determination -
I
I
I
20
15 IX. l0
0
-o.1 0.0 O.l -o.1 0.o 0.1 -O.l 0.0 o.1 Symmetric stretch, ~, Bending, ~2, Asymmetric stretch, ~3, reduced units reduced units reduced units Fig. 1. Distribution of normal mode coordinates for the bending mode (the middle column) and two stretching modes (two side columns). Full lines corresponds to theoretical RRHO model, histograms to the simulation by the successive sampling method.
248
I. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
0.20 0.3
> 0.2 0.10
0.1
0.0 0.0
2.0
4.0 Ev),
6.0
Ev 2 or Ev3,
8.0
000
10.0
0.0
2.0
4.0
kcal/mol
6.0
Ev,
8.0
10.0
kcal/mol
Fig. 2. Distribution of normal mode energies calculated by the successive sampling method (histogram) and by the analytical formula for RRHO model.
Fig. 3. Distribution of the vibrational energy E,., calculated the successive sampling method. Full lines corresponds to theoretical RRHO model.
of the coefficients ci~ ,. Thus, Eqs. (15)-(17) produce all the necessary coefficients ci~ a, and, therefore, provide the transformation of normal modes and velocities to the Jacobi coordinates and velocities by the use of Eqs.
(13) and (14). The further transformation to the Cartesian coordinates {x n} and vibrational velocities {Gvib } of individual atoms is trivial, if the condition of the mass center being located in the origin is imposed. The arbitrariness which is still left is elimi-
I
I
I
I
I
I
I
I
I
I
I
I
h
I
1
I
E
q'--
0.4
b 0.3
0.2
0.1
0.0
0
1
2
3
4
5
6
Ej, kcal/mol
7
8
9
10
1
2
3
4
5
6
7
8
9
10
Ej, kcal/mol
Fig. 4. Distribution of adiabatic rotational energy for a harmonic model of SO 2 molecule for the total energy of 10 kcal/mol. (a) Markov chain sampling with the weighting function from Eq. (26); (b) result of the successive sampling method. Full curve corresponds to Eq. (7).
I. Rosenblum et al./ ChemicalPhysics213 (1996) 243-252 nated by three consecutive rotations of the whole molecule through three random Euler angles. (v) Afterwards, two rotational energies Ej, and Ej~ are sampled successively according to the probability densities
active parts of rotational energy are equal, correspondingly, to
e, = 8of, J2= 8eff(J 1
R + 1Ej~'
=E,, +Ejc+
dEj, Pe'e"(g")dEi"= 2 t E j . ( E - E È )
'
(20)
=eR - ej= (eA + G +gc)
(18)
-
dG Pe'e,'eJ~(Ejc)dEjc= 7r E~jcI E - E
249
- Ej - E j ~
E1 + E ~ c + R + I E J A
=R+IEJ.
(vi) The above set of rotational energies determine the components of the angular momentum about the principal axes of inertia:
(19)
J,~= 21/~I,,, When Ej. and Ejc are fixed, the third component of rotational energy, E j , is determined from the equation Ei~ = E - E v - Ej, - E j . The adiabatic and the
I
I
I
I
I
I
I
I
(22)
where m = A,B,C, and I m are the instantaneous values of the principal momenta of inertia.
I
]
I
I
I
E
I
~
I
I
I
0.04
0.03
0.02
0.01
0.00
010
20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 7 0 8 0 Ej, kcal/mol E), kcal/mol Fig. 5. The same as in Fig. 5,
(21)
but for a total energy of 100 kcal/mol.
90100
250
1. Rosenblum et a L / Chemical Physics 213 (1996) 2 4 3 - 2 5 2
(vii) With the above values of angular momentum components, the projections of the angular velocities onto the principal axes of inertia are calculated: (23)
to,~ = J , , / I m .
The directions of these axes in the space-fixed frame are found by the diagonalization of the instantaneous tensor of inertia. Thus, the angular velocity vector to was completely determined. The vector of angular velocity obtained in this way is used to calculate the rotational contribution to the linear velocity of an atom n: rot :
Pn
(24)
¢0 X r n ,
where G is the radius-vector of the nth atom from the origin (such choice renders the zero velocity of the mass center). The resulting set of rotational velocities {vr°t} makes the total rotational energy E r equal to E - E v. Besides, the instantaneous Eckart conditions ensure that the so-called vibrational velocities {Vnrib} are collinear to the radius-vectors r, and, therefore, orthogonal to the rotational velocities {v~°t} determined by Eq. (24). Consequently, the total velocity of each atom (defined as v i = v ivib + vfOt) yields the required value of the kinetic energy T = E - U, U being the molecular potential energy: m i v 2 : ]~ m i ( v / r i b ) 2
T=
- 2
2
= E~ + r v = E - r:.
mi(t)r°t) 2
+E
2 (25)
A comparison of distributions obtained by the successive sampling method with predictions of RRHO model are shown in Figs. 1-3. Since the agreement is very good we conclude that this method of sampling can be recommended for selection of initial conditions from the microcanonical distribution.
4. C o m p a r i s o n o f the successive s a m p l i n g m e t h o d with the M a r k o v chain s i m u l a t i o n for equiprobable distribution o f 'vibrational' energy over vibrations and active rotations
In this section we present a comparison of results obtained with the help of successive sampling method and the Markov chain simulation. We emphasize that
the latter, as described in Refs. [5,6], generates the microcanonical sampling for prescribed values of the total energy E and the total angular momentum j. Therefore, the comparison is carried on under these conditions. On the other hand, we compared both 'collisionless' and 'collision-induced' quantities: the histogram distributions for rotational adiabatic energies for a free SO 2 molecule and the mean energy transferred in the collisions of SO 2 molecule with Ar atoms.
4.1. Collisionless events
In Figs. 4 and 5 we compare the distributions of the adiabatic energy E~ for a model SO 2 molecule for total energies E of 10 and I00 kcal/mol. The 'effective' vibrational energy which is defined as the difference between the total energy and the energy of adiabatic rotation [2,11] actually includes the energy of vibrational modes and the energy of the active rotation. The distributions shown refer to the equiprobable distribution of the 'effective' vibrational energy, E,.ef f = E - E i, over the above modes. In these calculations, we use the value of B e f f = 0 . 0 0 9 4 ( g / t o o l - Az) - I . With this value of Beef, the maximum in the theoretical RRHO curve corresponds to R = 5.336. It should be noted here that the aforementioned value of Beff is slightly different from that used by Troe et al. [11,12] which was found from the condition Beff=E/J . 2. . . . B e f f = 0 . 0 1 ( m o l / g . A2). Our way of determination of Beff was as follows. We defined the instantaneous Beff as (j2/l B +JZ/1c)/(J2 + J ~ ) , where JB and arc are the projections of the total angular momentum on the instantaneous axes of inertia and I~ and I c are the corresponding principal moments of inertia. This value oscillates slightly due to the non-rigidity of the molecule. We used its mean value, obtained by the averaging over approximately 100 vibrations; it is close to the equilibrium value Beq = ( l ~ q + I~q ) / 4 I ~ q 1~q. The panels a in Figs. 4 and 5 correspond to the Markov chain simulation [1], where the sampling of coordinates was accomplished with the use of the weighting function W ( q ) ot [ E -
U( q)] ~/2+W2.
(26)
1. Rosenblum et al./ Chemical Physics 213 (1996) 243-252
251
Table 1 Comparison of the transfer of different components of energy in collisions of SO 2 molecule with A r for two sampling methods. Total molecular energies 10, 50 and 100 k c a l / m o l , heat bath temperature 300 K, harmonic intramolecular potential. All values are in k c a l / m o l Energy transferred
E = 10 k c a l / m o l
E = 50 k c a l / m o l
E = 100 k c a l / m o l
SS method
SNF method
SS method
SNF method
SS method
SNF method
- 0.382
- 0.390
- 1.432
- 1.507
- 2.633
- 2.484
(AEj)
-0.139
-0.142
-0.417
-0.467
-0.603
-0.673
(AE,,eff)
-0.243
-0.248
- 1.015
- 1.039
-2.031
- 1.811
( A E R)
- 0.176
- 0.177
- 0.570
- 0.571
- 0.843
- 0.834
( A E, ) (AEk)
- 0.206 -0.038
- 0.213 -0.035
- 0.862 -0.153
- 0.936 -0.104
- 1.789 -0.240
- 1.650 -0.210
This weighting function is appropriate for rotating nonlinear molecule with s vibrational degrees of freedom; in our case U(q) represents a harmonic potential over the set of q coordinates which describe three vibrational modes, s = 3. It is seen that the histogram obtained is well reproduced by the theoretical curve 1. The panels b in Figs. 4 and 5 correspond to the successive sampling method. Once again, we see a good agreement of the histogram distribution with theoretical predictions for RRHO model. One also sees that both methods yield very close results.
4.2. Collision-induced events We have calculated, by the standard trajectory method as implemented in our previous work [2], the mean energy transferred in collisions of SO 2 molecule with heat-bath Ar atoms using two different procedures for the sampling of initial conditions. The model used in these calculations, is described in detail in [2]. Briefly, the potential of the free SO 2 molecule was constructed as a superposition of atom-atom pair harmonic potentials. The interaction between SO 2 and Ar was taken as a sum pairwise Lennard-Jones potentials between Ar and the atoms of the molecule (the so-called SLJ model [ 11 ]). The
We note in passing that the histogram distribution over Ej presented in Fig. 2 of Ref. [1 l] does not show a well-developed m a x i m u m and can be approximately described by the distribution (5) over the whole range of energies. As indicated by Schranz [15], this m a y be due to poor statistics.
relative velocity of collider was sampled from the three-dimensional Maxwell flux distribution [ 16]. The trajectories started from a mass center separation of 5 A, and were calculated during approximately 25 ps. The average transfer of the total energy, ( A E ) , of the adiabatic rotational energy, (AEj), and of the 'effective' vibrational energy, ( A E v . e f f ) = ( A E ) ( A E j ) , in these collisions was calculated. Also, using the definition of the rotational energy [9,10], we determined the transfer of this total rotational energy, (AER), and of its 'active' part, (AEk). These values were compared to those obtained with the Markov chain sampling procedure [5,6]. The results for molecular energies 10, 50 and 100 kcal/mol and translational temperature T = 300 K are presented in Table 1. The consistency of the quantities, obtained by the new, successive sampling, method ('SS'), with the results of the traditional method of Schranz et al. ('SNF') [1] confirms the applicability of our method.
5. Conclusions
In this paper we have presented a method of sampling of initial conditions for rotating triatomic molecules, which allows one to specify the energies of normal modes, and also the adiabatic and active parts of the rotational energy (the so-called successive sampling method). Within this method, the distribution of over vibrational phases and the Euler
252
1. Rosenblum et al. / Chemical Physics 213 (1996) 243-252
angles that describe the rotation is equiprobable. We checked this sampling against exact analytical results for the rigid rotor-harmonic oscillator model. Also, we have compared some results of this sampling with the appropriate counterparts derived by using the Markov chain sampling procedure, when the only specified quantities are the total energy and the total angular momentum. The good agreement, found in both kinds of comparison, indicates that the successive sampling procedure can be used for a more detailed study of intra- and intermolecular dynamics, e.g. collisionless intramolecular energy redistribution and collision-induced energy transfer under condition of non-uniform distribution of the vibrational energy over molecular modes.
Acknowledgements We acknowledge very fruitful discussions with Professor J. Troe and Dr. H.W. Schranz. This research was supported by the Fund for the Promotion o f R e s e a r c h at t h e T e c h n i o n .
References [1] G. Nyman, S. Nordholm and W. Schranz, J. Chem. Phys. 93 (1990) 6767. [2] I. Koifman, E.I. Dashevskaya, E.E. Nikitin and J. Troe, J. Phys. Chem. 99 (1995) 15348. [3] I. Oref, Chem. Phys. 187 (1994) 163. [4] L. Hase and D.G. Buckowski, Chem. Phys. Lett. 74 (1980) 284. [5] H.W. Schranz, S. Nordholm and F. Freaser, Chem. Phys. 108 (1986) 69. [6] H.W. Schranz, S. Nordholm and G. Nyman, J. Chem. Phys. 94 (1991) 1487. [7] T. Lenzer, K. Luther, J. Troe, R.G. Gilbert and K.F. Lim, J. Chem. Phys. 103 (1995) 626. [8] V. Bernshtein and I. Oref, J. Chem. Phys. 104 (1996) 1958. [9] J. Jellinek and D.H. Li, Phys. Rev. Lett. 62 (1989) 241. [10] J. Jellinek and D.H. Li, Chem. Phys. Lett. 169 (1990) 380. [11] H. Hippler, H.W. Schranz and J. Troe, J. Phys. Chem. 90 (1986) 6158. [12] H.W. Schranz and J, Troe, J. Phys. Chem. 90 (1986) 6168. [13] E.E. Nikitin, Theory of elementary atomic and molecular processes in gases (Clarendon Press, Oxford, 1974). [14] E.B. Wilson, J.C. Decius and P.C. Cross, Molecular vibrations (McGraw-Hill, New York, 1955). [15] H.W. Schranz, private communication. [16] K.F. Lim and R.G. Gilbert, J. Phys. Chem. 94 (1990) 72.