Computer Simulations of a Rough Sphere Fluid

a Rough Sphere Fluid

J.W. Lyklema

Computer Simulations of

a Rough Sphere Fluid

I'"Hi; II III

iiiifiiif

IIIIIII ilUI nil IIII. III III III lillUtlli

o

o Ui C -J BIBLIOTHEEK TU Delft P 1789 4052 C 533180

a Rough Sphere Fluid

Proefschrlft ter verkrijging van

de graad van doctor in de

technische wetenschappen

aan de Technische Hogeschool Delft,

op gezag van de rector magnificus

prof. dr. ir. F.J. Kievits,

voor een commissie aangewezen

door het college van dekanen

te verdedigen op

woensdag 6 September 1978

te 16.00 uur door

Jans Wijbe Lyklema,

chemisch doctorandus

geboren te Groningen

<^^>? ^OS^

The work described in this thesis is part of the research program of the Foundation

for Fundamental Research of Matter (F.O.M.), which is financially supported by the

Netherlands Organization for Pure Research (Z.W.O.).

Introduction 1 References 4" I. The Rough Sphere Model 5

1.1. Conservation Laws 6 1.2. The Collision Model 10

1.3. Invariances 13

1.3.a. Time Reversal Invarianae 13

I.3.b. Spatial Inversion Invarianae 14

I.3.a. Permutation Invarianae 14

1.4. Conclusion 15 References 18

II. The Correlation Functions 19 11.1. Molecular Intermediate Incoherent Scattering Function 20

II.1.a. Rotational Symnetry of the Spheres 21

Il.l.b, Isotropy in Space 22

g

11.2. Approximations for F (t) 25

II.2.a. Expansion of ^(t) '^n Legendre Polynomials 25

11.3. Free Rotator 27 11.4. The Gaussian Approximation 28

11.5. Velocity Autocorrelation Functions 30

References 32 III. The Computer Experiment 33

111.1. Initial Conditions 34

111.2. Equilibrium 34 111.3. Track Generation 35 111.4. Time Correlation Functions 36

IV. 1.a. Series Trunoation 41

IV.l.b, Product Assumption 41

IV.2. The Gaussian Approximation 43

IV.2.a. The Translational Part 44

IV.2.b. The Rotational Part 44

IV.3. Influence of the Roughness Parameter tfi 46 IV.4. Comparison with a Neutron Scattering Experiment 50

IV.5. The Orientation Correlation Functions 53 IV.6. The Velocity Autocorrelation Functions 56

IV.7. Accuracy Estimation 59

References 61

V. Comparison with Stochastic Models 62

V.l. The J-Diffusion Model 62 V.2. The Cumulant Expansion 64 V.3. The Stochastic Liouville Equation 71

References 77 Appendix A 78 Appendix B 90 Summary 99 Samenvatting 101 Acknowledgements 103

Introduction

In this thesis we will discuss computer simulations of a molecular fluid. The purpose of this computer experiment, as in general for

simulations, is twofold viz. o

1 . It provides the ideal experiment for the theory with all parameters perfectly under control cind it contains all the necessary input data needed to compare theory with the numerical results (see chapter V ) .

2 . It can be seen as the ideal theory for the experiment. In practically all cases the theory for a many body system can not be carried out from first principles to a level where reliable predictions for the experiment can be made. With the interactions between the particles as only input one can decide from the machine results and the experimental data in how far this model system describes the real experiment. Also the assumptions which are made in the course of explaining the experimental results can be studied carefully in this way

(see chapter IV).

After the molecular dynamics studies of Alder and Wainwright (1958) for hard sphere systems and Rahman (1964) for Lennard-Jones potentials to describe single atomic fluids, the interest to simulate more complex many body systems on a con5>uter has rapidly increased. One of the possible extensions is to include rotational degrees of freedom to describe

molecular fluids, which is the subject of this thesis.

As one of the objects of this study is the comparison with a real molecular fluid viz. liquid neopentane (C(CH ) ) we have to search for a potential which can describe these spherical molecules. Since neopentane has no dipole- or quadrupolemoment (by reasons of symmetry) it is not possible to use a Stockmayer potential (Lennard-Jones potential plus a dipole-dipole term etc., see Harp and Berne (1970)). A possibility

is to represent every atom by a Lennard-Jones potential (Barojas, J., Levesque, D. and Quentrec, B., 1973) but this is a very time consuming procedure since the molecule consists of 17 atoms. However it may be possible to describe every methyl group by a Lennard-Jones potential reducing the number of potentials to 5. The alternative, which has the advantage of being more tractable for theory, is to take a hard core potential including a mechanism to change the angular velocities during a collision. We have decided in favour of this somewhat un-realistic potential because its simple equations of motion allow us to study a large number of physically different situations for not too much computer time.

The so called rough sphere model used in this thesis has already been proposed by Bryan (1894) and Pidduck (1922). In this model the sign of the relative velocity of the spheres at their point of impact is

reversed. This model has recently been used by Subramanian et.al. (1974, 1975), O'Dell and Berne (1975) and Kestemont (1976) to calculate the

linear and angular velocity autocorrelation functions and the orientational correlation functions. We have extended this model by introducing a

continuously varying roughness parameter, which includes the smooth and completely rough sphere as limiting cases (see also Offerhaus, M.J. 1978). This model has been first proposed by Sather and Dahler (1951) , although in quite a different formulation, and will be described in chapter I.

In addition to the correlation functions mentioned above, we have calculated the intermediate scattering function and some related

correlation functions (see chapter II) . These time correlation functions are the basic quantities which describe the dynamical behaviour of fluids in equilibrium and what is more inportant some of them are measiirable by experiment (Sears, V.F., 1966; Gordon, R.G., 1968).

In chapter III we describe briefly the way our molecular dynamics experiment is performed. For more details on computer experiments we refer the reader to review articles by Wood (1975) and Erpenbeck and Wood (1977).

only a fraction of the total number of curves obtained, which we consider to be representative. These and all the other data are tabulated elsewhere and are available on request. We show the correlation functions described in chapter II and discuss the validity of the various approximations, which can be made, as a function of density, roughness and wave number. The influence of the roughness parameter and the possibility to adjust it so as to explain a neutron scattering experiment on neopentane is analyzed, which serves as an illustration how computer simulations can be used as the ideal flexible theory.

In the last chapter we compare the predictions from some simple stochastic models (for a review see van Kampen (1976)) with the computer results. This is only done for one suitable correlation function, but we do not expect that the conclusions would be different for the other correlation functions. Sometimes we use numerically calculated correlation functions as input data for calculating other more complex correlation functions, illustrating the use of a molecular dynamics experiment as-the ideal experiment with information on all intermediate levels. The conclusion of this chapter is that, although the existing theories give a qualitative agreement with the computer results, one needs a better theory to give a quantitative description of the correlation function.

References

1. Alder, B.J. and Wainwright, T.E. (1958),in Transport Processes in Stat. Mech., Ed. I. Prigogine, Interscience, Brussels.

2. Barojas, J., Levesque, D. and Quentrec, B. (1973), Phys. Rev. A7, 1092.

3. Bryan, G.H. (1894), Report British Association for the Advancement of Science 64, 64.

4. Erpenbeck, J.J. and Wood, W.W. (1977), in Modern Theoretical Chemistry 6, Ed. B.J. Berne, Plenum Press, New York 1977. 5. Gordon, R.G. (1968), in Advances in Magnetic Resonance 3, Ed.

J.S. Waugh, Academic Press, New York.

6. Harp, G.D. and Berne, B.J. (1970), Phys. Rev. A2, 975. 7. van Kampen, N.G. (1976), Physics Reports, 2£, 171. 8. Kestemont, E. (1976), J. Phys. C£, 2651.

9. O'Dell, J. and Berne, B.J. (1975), J. Chem. Phys. 63_, 2376. 10. offerhaus, M.J. (to be published, 1978) Physica.

11. Pidduck, F.B. (1922), Proc. Roy. Soc. AlOl, 101. 12. Rahman, A. (1964), Phys. Rev. 136, A405.

13. sears, V.F. (1966), Can. J. Phys. 44_, 1279, 1299; 45_, 237.

14. Subramanian, G., Levitt, D. and Davis, H.T. (1974), J. Chem. Phys. 60, 591.

15. Subramanian, G. and Davis, H.T. (1975), Phys. Rev. All, 1430. 16. Wood, W.W. (1975),in Fundamental Problems in Statistical Mechanics

I. The Hough Sphere Model

A rough sphere is described by its mass m, diameter a and moment of inertia I. Another important parameter is the loading factor

2

K(=4l/ma ) , characterizing the radially symmetric mass distribution. To describe the motion of such a body we need twelve coordinates viz.

- > • - > • -*

the position r, the orientation u, the linear velocity v and the angular ->•

velocity o). The position and the orientation of the i-th sphere only change between collisions according to

r^(t) = v^(t)

and

u. (t) = to. (t) X u. (t)

where the dot means a time derivative. The velocities only change instant-aneously as the result of a collision. A hard sphere collision is an instcuitaneous one point event. This means that only two spheres can touch each other at the same time and that they then have one point in common. So the hard sphere dynamics is a sequence of binary collisions. To find the velocities after the collision, for a pair (1,2) indicated as:

( V ^ , V ^ , 10^, (1)2) -* ( V j , V ^ , tO^, 0)^)

we need twelve equations. The unprimed quantities denote velocities before the collision and the primed quantities denote velocities after the

I . l . Conservation Laws

Three equations follow from conservation of linear momentum, giving

V = - ( v '

1 ^ 2 -2> (1)

The following six equations follow from conservation of angular momentum and using the concept of a point collision, which applies for hard spheres Let us denote by L. the angular momentum of sphere 1 with respect to the

origin 0. Now we can write for the change of L due to the collision

H

dt q X F = q X dtF Collision Collision Here F is the force and q X F is the torque exerted by sphere 2 on sphere 1 during the collision. The integrated force is equal to the change in linear momentum

dt F = m(v' - V.) Collision

Thus we can write for the change in

„• , - -- . . , , angular momentum for sphere 1:

Fvgure 1. Coll%ston between ^ '^ particles 1 and 2 at the

- > • - > • - > - > • - *

j - j . • j . / ^ / i - H - L , = q X m(v,' - v.)

contact povnt Q. 0 %s 1 1 ^ 1 1

the origin.

and for sphere 2:

L^ - L^ = q X m(v^ - V2) = - q X m(v' - v^) = -(L^ - L^)

This is just conservation of angular momentiim. If we now choose the origin O at the contact point Q we have

5 = 0

S - h = 4 - ^2 = °

From the definition of angular momentum

Lj = r X mv^ + lo).

and from figure (1)

-»•-»• - > - > •

r^ = a/2, r2 = -a/2

we observe that this choice in^lies

L^ - Lj = (ma/2) X (vj - Vj) + I(wJ - u^) = 0

and

L^ - L^ = -(ma/2) X (v^ - V2) + K u ^ - m^) = 0

For the change in angular velocity we thus find using conservation of linear momentum (equation 1)

(o| - (0^ = 0)^ - ^2 = -(nia/2l) X (vj - v^) (

That this result is not a consequence of our particular choice of the origin can be seen if we take q arbitrary

We then have

L^ - Lj = r J X m(v' - Vj) - (ma/2) X (vj[ - v^)

and from the definition of angular momentum

- > • - > • - » • - * • - * • - »•- *•, LJ - L^ = rj X m(v' - v^) + 1(0)^ - u^)

Subtracting these two equations again results in equation (2).

We have now found 9 equations, which arise from conservation laws having their origin in the invariance of the system under a symmetry operation such as translation or rotation of the system or a rotation of each sphere separately.

To proceed we introduce the new set of velocities (g, G, V, sT) defined as, g = v + a/2 X 0) G = (K/K+1)^^-^ (U + (2/Ka^)a X v) (3) V = (v^ + V2)/2 ^ = (U2 - to^)/2 with V = (v^ - V2) 0) = (00^ + ti^)

This is a transformation with Jacobicin unity. From equations (1) and (2) we now have

V = V

^' = Q (4)

- » • -*•

G' = G

So the problem is reduced to finding a collision model for the vector g which is the difference velocity between the colliding points on the two spheres. Only one conservation law viz., energy conservation

2 2 2 2 2 2 2 2

AE = >an(v^ + v^ " '^1 ~ '^2' "*" ''•'•'"l "^ "2 - "^ ~ 1^2' " °

has not been used. Written in the new varicJsles this reads as

AE = AE,, + AE, = 0 _E.. + AE.

AE|| = m/4 (g|[ - g.| ) (5)

AE^ = (mic/4(K+l)) (g^^ - g^)

- * • - » • - > •

Here g.. and g. are the parallel and perpendicular components of g with respect to a, the collision axis. We note in passing that one obvious and non trivial way to satisfy energy conservation is to take

g"

T h i s is the straight forward extension of the smooth sphere collision

- > •

to rough spheres: complete reversal of the difference velocity g in the contact point Q (Chapman, S. and Cowling, T.G., 1970).

1.2. The Collision Model

For an acceptable collision model still another condition has to be fulfilled, which is of a statistical nature viz., equipartition between rotational and translational degrees of freedom

<AE > = <AE > = 0 T R

where

AE^ = 'im(vj^ + v'^ - vj - V2) (6)

AEj^ = isKu)'^ + 03^^ - 1^1 - '"2'

and the bracket denotes an average over a Maxwell-Boltzmann distribution,

< ^ V = '2fk7)' ' 2 ^ ^ '

dv. dv. dtu. do)2 A E ^ e -E/kT

with the energy E given by

E = m(v^ + V2)/2 + 1(00^ + W2)/2

We can easily transform this expression into the new variables because the Jacobian is unity. We then have

and

^^ = |t^'^||'-^?)-<^)'''ng--g)x3}„]

In the integration the last term in AE drops out as it contains G, linearly (assuming that g' does not depend on G. ) and G, appears

quadratically in E. So to guareintee equipartition we have the condition

<g,;^ - g,f> = 0 (7)

which IS fulfilled if we take

(8)

analogous to the smooth sphere collision.

In section (4)we will show that this choice is in fact the only possibility if one wants to keep a linear relation between g' and g.

Returning now to the energy conservation law we see that equation (8) implies

AEj^ = 0

gj; = g^

The most general solution for this equation is a rotation around a over an angle ¥> of g, in the plane perpendicular to a

gj| = cos I/" gj^ + sin f {a X g^) = Ri'(>)g^ (9)

For 1^ = 0 we find again the smooth sphere collision model:

91 = g ,

and for f = u

H

=

-^1

The choice of i^ is to a large extent arbitrary. We will take the magnitude of ¥> (denoted by ¥> ) as a parameter for the roughness

o

,p = S <fi 0 < ¥ > < TT (10)

o o

The sign S of tp, defining the sense of the rotation is taken as

S = sign {a.(a).+n)„)} = sign (o.G)

as will be motivated in the next section. We mention here only that this choice of sign is independent of G. in accordance with the argument leading to equation (7).

We now have twelve equations, fulfilling all conservation laws, from which the new velocities can be calculated. From the definition of g equation (3), we immediately find

(^1 - ^ l ) l = ^ ^ ( ^ ' - ^ ' l

Thus we have from equations (1), (2), (8) and (9)

V' - v^ = - (v- - V2) - > • - > • - » • - » • -y -V -V (j),' - u, = w! - a)„ = -(ma/2l) x (vl - v.) (11) (v' - V ) .a = -(v - V2).a (v| - v^)^ = h —• {co3P-i)g^ + sin¥'(0X gj^') }

In fact ^ may be a function of the relative quantities which describe - - > - - » • - > • - > •

the collision viz. a, (v. - v„) and (o), + 10.), but for calculational simplicity we will assume that f is a constant parameter.

1.3. Invariances

There are still a number of other symmetry properties that have to be fulfilled viz. time reversal, spatial inversion and permutation invariance. In this section we will discuss these in connection with the choice(equation 10) we made for f.

I.3.a. Time Reversal Invarianae

Time reversal symmetry means that the initial and final positions must be the same and the initial and final velocities the opposite if we reverse the velocities halfway. So the following sequence of events must be checked

->• Collision -»-, velocity -»•„ , -*•,. Collision -*•,,,? ->• g ^ g' •: > g" = -g') > g " ' = -g

' inversion '

- * •

For the parallel component g., this is obviously true because two times a sign flip gives the initial situation back. For the perpendicular component we have to check if

-g = g ' " = R(¥'")g"

which i s t r u e b e c a u s e

cuid

s i g n ( a . G " ) ¥ ' = - s i g n ( a . G ) ¥ ' = -V*

I.S.b. Spatial Inversion Invarianae

The second invariance property to be fulfilled is invariance under spatial inversion. Under this symmetry operation the coordinates trans-form according to

(12)

G -* G

and S according to equation (10) S -»• -S.

So we have to find out whether a trajectory starting from a spatially inverted configuration develops in the spatially inverted trajectory. Again for gj this is obviously true. To check this condition for g. we apply a spatial inversion to equation (9). We get

->• -> ~ - -+• - -•

-g.' = -cos(-¥')g. + sin(-¥') (a X g) = -{cos^f g, + sinv(a X g) }

Thus t h i s equation trauisforms i n t o i t s e l f and invariance under s p a t i a l i n v e r s i o n i s f u l f i l l e d .

I.3.a. Permutation Invarianae

The last condition to be checked is invariance under the change of numbering the particles. It means that the development of the trajectory is independent of the choice of numbering the particles. Under a permuta-tion the coordinates transform also according to equapermuta-tion (12), so by the same arguments as given above this condition is fulfilled.

That invariance under permutation symmetry is not a trivial condition can be seen if we take another choice for S, for instance

S = sign{a.(u. - u)_) }

This choice fulfils spatial inversion and time reversal invariance. However, under peinnutation S does not change, thus arguing again as after equation (12) the permutation invariance is not fulfilled.

We now have shown that the choice

S = sign (a .G)

is a proper choice because it fulfils all invariance properties.

1.4. Conclusion

We have obtained a collision model for rough hard spheres based on the assumptions that

o . -••

1 . the relative linear velocity along a (g,, ) changes upon a collision and does not change magnitude.

o -> 2 . the perpendicular relative velocity g. rotates over a constant

angle v' , the sense of rotation depending on the sign S given by equation (10).

discuss law equation (5)

o

To discuss assumption 1 we go back to the energy conservation

«„ in r , ,2 2. K , ,2 2. , AE = 4 {(g,,- -g,,) + ^ (gi - g ^ ) } = 0

If we want to allow for energy transfer (Offerhaus, M.J.» 1974) from

E|| -> E^

a g,,

gj- = R(^)g^ [ 1 - ^ (a^^-Dgpgf]'' (13)

The factor P, which can take the values plus or minus one, is included to obey invariance properties. For time reversal symmetry P ^ -P and for spatial inversion and permutation invariance P ^ P. These equations are rather complicated and have the difficulty that the square root in equation (13) can become negative. Another problem arises if we look to the equipartition condition equation (7)

<g,\^ - gii> = 0

(7)

/ 2P ,, 2 <(a -l)g,, > = 0

This condition shows that a must be a function of at least one of the variables g, G, V and U. For such a model there is no physical evidence and besides that it is less attractive from a computational point of

view. The only simple choice one can make is a = - 1 , which is our collision model.

The second assumption is made because a constant if is the easiest o

computational choice which is allowed. Moreover it is hard to decide what kind of function ^ should be from a molecular viewpoint. The sense of rotation as defined by S = sign(a.G) is the only simple choice which fulfils all the invariance properties the system should have. Other choices for ^ are of course possible for instance

ifi = -n tanh{ (o.G)/to } o

f = -n sin{ (a.G)/u) } o

where u can be regarded as a roughness parameter.

The model as described in the foregoing sections was first proposed by Sather and Dahler (1961) although they did not formulate the model

-*• a

as a rotation of g. , which may be the reason why they did not define the sense of rotation.

Other models, which include a roughness parameter have been proposed by Berne (1977) viz. the rough domain model, the kinetic stick model and the rough screwball. For the first two models molecular dj^amics cal-culations have been performed (Pangali, C.S. and Berne, B.J., 1977). These models have in common, that they are composed of two other models viz. the smooth sphere and the con5)letely rough sphere model and a criterion to decide whether the collision is smooth or rough. For the first model this criterion is in fact a random number generator, while for the second model it is a certain property of the colliding particles. So only on the average there is an intermediate roughness, whereas we have roughness introduced on the microscopic level. The rough screwball model is proposed to describe the dynamics of a chiral molecule in a nonchiral solvent. It is just the rough sphere model where the sense of rotation is defined as the chirality of the molecule.

The rough sphere model has also been used by Offerhaus (1974) to study collision processes between atoms and molecules, but his sign convention S = sign(a.i2) does not fulfil the permutation invariance requirement (Offerhaus, M.J., 1978).

The resemblance between our model and the Sather-Dahler model was pointed out to us by Dr. M.J. Offerhaus.

References

1.1. Berne, B.J. (1977), J. Chem. Phys. 66, 2821.

1.2. Chapman, S. and Cowling, T.G. (1970), The Mathematical Theory of Non-Uniform Gases, Cambridge University Press, Cambridge.

1.3. Offerhaus, M.J. (1974), Proc. of the IXth International Symposium on Rarefied Gas Dynamics, Gottingen, W. Germany.

1.4. Offerhaus, M.J. (1978), preprint to be published in Physica A. 1.5. Pangali, C.S. and Berne, B.J. (1977), J. Chem Phys. 67_, 4571. 1.6. Sather, N.F. and Dahler, J.S. (1961), J. Chem. Phys. 35_, 2029.

II. The Correlation Functions

An important function in the theory of the microscopic behaviour of fluids is the intermediate incoherent scattering function

F ^ t ) ='-Z <eik.5i(t)^-ik.5i(0)^ (^j

k N j^

Here k is a wave vector. N is the total number of particles and the

average <..> is over a canonical ensemble defined as

<f(r)> = dr p(r)f(r)

with

p(r) = exp -BH(r)/ dr exp -eH(r) ,

8 = 1/k T, k = Boltzmemn's constant and T the eibsolute temperature.

a B

r stands for the phase space coordinates of the system. In the case of a fluid of smooth hard spheres the relevant phase space coordinates are (q,,v,),the position of the center of mass and the linear velocity of the i-th sphere respectively. The spatial fourier transform of (1)

G^(?,t) = 3

(2Tr)

-»• -••-»• s dr exp ik.r F (t)

= 1 E <6{? _ (5^(t) - q^(0))}>

is called the Van Hove self correlation function. It is the probability that a particle travels a distance r in a time t. Thus it is a measure for the self diffusion of a particle in a fluid. The fourier transform

s **"

of (1) with respect to the time S (k,u), is called the incoherent scattering function and can be measured by inelastic slow neutron scattering.

For an extension to rough hard spheres one has to include the rotational

->•->->-»

-degrees of freedom and the phase space coordinates become (r.,v.,u.,a).):

->• ->•

the position of the centre of mass r., the linear velocity v., the

-•• ->•

orientation of the sphere u. and the angular velocity u. of the i-th particle.

II.1. Molecular Intermediate Incoherent Scattering Function

To include rotational motion in the scattering function (1) we write

q. (t) = r. (t) + u. (t) la 1 la

where q. (t) is the position of point a on the surface of sphere i at time t. r.(t) is the position of the center of mass of sphere i at time

t and u. (t) is the vector pointing from the center of mass of sphere i to point a on the surface of the sphere. Because of the spherical symmetry of the rough sphere all directions are equivalent and there-fore we can drop the index a. So we write for the scattering function

^s,^. ^ik.(?. (t) - ri(0) ik. (ui(t) - ui (0) ^ F, (t) = < e J- - ^ e -^ ^ >

where the summation over all particles has disappeared because they all contribute equally to the average.

For the sake of convenience we define

r^(t) - r^(0)

and

where the caret means a unit vector: u = u/|u|. With these definitions we have

r-it) =<eii^.5e^!5^-(-(t)-u(0))^ ^^^

In the next sections this formula will be investigated for its

symmetry-properties.

II.1.a. Rotational Symmetry of the Spheres

As explained in chapter I the direction u is an irrelevant vector for the dynamical evolution of the system; so we can perform this part of the ensemble average. To do this we define

G(t) = R(t) .G(0)

and

k(t) = k.R(t)

Here R(t) is a 3 by 3 rotation matrix which rotates u(0) to u(t) Thus ' - "'•""•

k.(G(t) - u(0)) = k.(R(t)-l') .G(0) = (k(t)-k) .u(0))

The average over u results in the following integral

4TT du exp {j- (k(t)-k).u} = J Q ( | lk(t)-k|) = JQ(|2- /2-2k.R(t).k) /2-2k.I

Il.l.b. Isotropy in Space

s "*•

The independence of F (t) on the direction of k can be used to simplify the expression by performing an integration over k. For this purpose we write

Y is a spherical harmonic and the asterix denotes a complex conjugate.

3

We now have for F, (t) k

^k't) '^'^

iW^'^^^'^^^^ ^'^^h ]^ ^l^^^

j^(f /2-2k.lk)>

A rotating vector evolves in time as

k.R(t) = (k.^(t))|(t) + cos ct)(t){ (^(t)x k)X *(t) } + sin <))(t) {^(t)xk}

where (j> (t) is the direction of the rotation cixis at time t and (t> (t) denotes the total angle over which the sphere is rotated. This (f) (t) is the net rotation which is built up from all the successive intercollisional rotations with their different angular velocities and free rotating

times. Thus we can write

1 - k.R(t) .k = (l-cos(t>(t)) (l-k.^^(t))

If we now choose ^(t) as the polar axis we can perform the integration over the azimuth directly yielding

k r * k ^ I (^) = ^ m , 0 ' ^ ^ <i^-*'^^)

Because

we are left with the following expression

F^(t) = Zd)*- ^ j ! ^ <j^(kd)P^(i(t).d)

•iSk

^ ^ " V J I ' ' = ° ^ \ >

j _ ( ^ /2(l-cos(t.(t) (1-cos^e, )>

•"O 2

Now the integral over 9 gives zero for odd I because of the following symmetry property of the Legendre polynomial

Pj^(-x) = (-1) Pj^(x)

So we can write for the scattering function

F^(t) = J o (4n+l)(-l) <J2jkd)P2jd.*(t)L2jz)> ₍₃₎ with L. (z) = ^ 2n 2 de P2n(9) 3o (^•''^^^ and z = ka|sin(ti(t)/2

A power expansion of j gives for L„ (z) 2m

L- (z) = ^ Z(-l)°'_{2n 2 m (2m+l):} ^•x—rr. de P2n(e)(i-e^)"

This integration can be done (Abramowitz, M. and Stegun,I.A., 1965) and gives

T / ^ V „m 2 m

L (z) = E S.T z 2n m=n 2n

\ "

0 1 2 3 4 0 1 0 0 0 0 1 - 3-2 3 - ^ 5 - ^ 0 0 0 2 - 2 3-h-h-' 0 0 3 - 3 - 2 5 - 2 7 - 2 3 - ^ 5 - 2 7 - 2 9 - 1 - 5 - 2 7 - 2 9 - ^ 1 1 - 1 3 - 1 5 - 1 7 - 2 9 - 1 1 1 - 1 1 3 - 1 0

Table I. Values for i'

2n

with

m ^ / _ n ™ 2L (mi)

2n ^ 2 (2m+l) ; (m+n+Js): (m-n):n:(-n-^):

The summation starts at m = n because for n > m l/(m-n)l = 0. So the first term in the series is

22"-'l(2n):

i" = (_1)" l

Jh

2n ^ 2 (2n+l) : (2n+»5) : (-»2-n) ; (4n+2) For the other terras we have the recurrence relation:

C /iC ^= -m/((2m+l) (m-n) (2m+2n+l)) 2n 2n

(see also Table I ) .

TT

For n = 0 one finds L (z) = — H (z)/z where H (z) is the Struve function of order zero (Abramowitz, M. and Stegun, I.A., 1965) . For higher n no such simple relation has been found.

g

II.2. Approximations for F (t)

A frequently made assumption for the scattering function equation (2) is the decoupling between the translational and the rotational part of this function .:;*• t i T k.(G(t)-G(0) s, , ik.d Z F (t) = <e e >

i k d i f k.(G(t)-u(0))

fb <e'-'^-S<e ^ > s Fj^(t) F^(t) (4)

By performing the same integrations as before this can be written as:

F^(t) = <JQ(kd)>

(5)

^k(^> = J%(^)/^>

From equation (3) and equation (5) one can see that this approximation in fact consists of two parts viz.

1 . Truncation of the series expansion after n = 0

F^(t) % j<JQ(kd)H^(z)/z>

This amounts to ignoring the correlations between the directions of the total rotation ^(t) and the total tramslation d ( t ) .

o

2 . Approximation of the average of the product by the product of averages for the functions involving the total displacement d(t) and the

total rotation (j) (t) .

II.2.a. Expansion of PV,(t) in Legendre Polynomials

T3

If we go back to the definition of F, (t)

F«(t) = <ei7k.(G(t)-G(o))^ k

and integrate over the direction of k we are left with:

Ff(t) = < j , ( ^ /2-2u(t).u(0))> k 0 2

This formula can be expanded in Legendre Polynomials in u(t).u(01 as

F^(t) = Z (2S,+ l)j^(ji) Cj^(t) (6)

with

C^(t) = <P^(u(t).u(0))>

The orientation correlation functions C (t) for H = 1,2 can b e measured by respectively infrared absorption and Raman scattering. The orientation correlation functions are independent of the initial orientation u ( 0 ) . Therefore we can perform this part of the ensemble average to obtain

%

expressions in terms of the rotation matrix R(t) defined by

u(t) = R(t)u(0)

For S, = 1 this gives

-»- 1 C, (t) = <u(0).R(t).u(0)> = I <u.R. . (t)u.> = — <TrR>

1 ij 1 i: 3 3 The second orientation correlation function is given a s :

C2(t) = i- [ 3<(u{0) .R(t) .u(0))2> - 1]

Performing the integration gives for the average

<(u(0).R(t).u{0))2> = E <R. . (t)R, „ (t)u.u.u, u > _j • 1] Ki. 1 J k J. kZ = 1 ^ 1 <R, . (t)R, . (t)> (6. .6, „+6., 6. +<5.„6., ) _{15 j_j xj Ti{,' 1] k£ ik ]£ iJ, fk.} , kl 1 2 2 = — {<(TrR) > + 3 + <TrR >}

Thus we have

C2(t) = Y^ {<(TrR)2> + <TrR2> - 2}

2 2 2

Because TrR = 1 + 2 cos 2^ = 4 cos v - 1 = (TrR) - 2 TrR we can write this as

1 2 C„(t) = r <(TrR) - TrR - 1>

2 3

For the higher orientation correlation functions we can proceed in the same way. We find for S. = 3

1 3 2

C (t) = ^ <(TrR) - 2(TrR) - TrR + 1>

Higher order functions will not be studied in this thesis.

II.3. Free Rotator

In the case of free rotators these functions can be calculated explicitly. We will give here the results because it is important to know how they behave in order to compare them with the theoretical approximations. Using that (f(t) = ut we can write

2

dto e-^"''" '2 (2coswt + 1)

S''I =!<-«>= 'g'^'H

Straightforward integration gives

C^(t) = j(l + 2(l-t2/6i)exp - t2/2Bl)

For general H we have for the orientation correlation function

Cj^(t) = 2I73; S (1 - m2t2/Bl)exp(- m2t2/2ei) m=-Jl

and for the rotational scattering function

F ^ t ) = 2 j 2 ( ^ ) E (1 - m2t2/Bl)exp(-m2t2/26l)

•k^ ji ^i'2

m=-li

In the limit t -> <» only the m = 0 term (the component along the rotation axis) survives thus

Cj^(t -> ") = l/(2Jl+l) (7)

and

F^(t ^ •» ) = E j 2 { ^ ) = si(ka)/ka K S, >!• 2

where Si(x) = / dt sint/t is the Sine integral. This result is different from that of the rotational scattering function F (t). Here we have for

t -> CO

„R,^ . ik.u(t) -ik.u(O) .2,ka. Fj^(t ^ 00) = <e ' '> <e > = 3. (—)

That we have only one term here instead of a infinite series as in the free rotator case can be seen from the expansion (see equation 6 ) .

R 9 kry F^(t) = I (2£+l)j^(-)C^(t)

Now because C (t ->• <») = 0 and C (t) = 1, in contrast with the free R 2 ka

rotator limit given by (7) , we have F. (t ->•">)= j_(-r—) .

This function does not vanish in the limit t ^ •" because the orientation space is finite, so that the orientation vector can not diffuse away as the position vector can do.

II.4. The Gaussian Approximation

Another frequently made assumption is the gaussian approximation which may be seen as the first nontrivial term in a cumulant expansion in k.

For smooth spheres where the orientation plays no role this approximation amounts to

2

F^(t) = exp - ^ <(q(t) - q(0))2>

k b

For a rough sphere where

q(t) = ?(t) + u(t)

we have

2

F^(t) = exp - | - <{(?(t)-?(0)) + (u(t)-u(0))}2>

Due to the spherical symmetry the cross product gives a zero contribution thus it automatically satisfies the decoupling approximation

F^(t) = F ^ (t)R^<^(t) where o T G k -> ->• ? Fj^^(t) = exp - ^ <(r(t)-r(0) > and F^=(t) = exp - - ^ ^ (1 - C^(t))

Here we used the following equality 2

<(u(t)-u(0))2> = ~- <2-2u(t).u(0)> Because we have the relation

ft

we can construct the scattering function F (t) from knowledge of the velocity autocorrelation function (see section 5). There is no such

RG

simple description possible for F (t), there we have to know the orientation correlation function C.(t). It is worthwhile to check these approximations because they are often used to understand neutron

scattering data (Boutin, H. and Yip, S. , 1968 ; Agrawal, A.K. and Yip, S.^ 1968).

II.5. Velocity Autocorrelation Functions

Two other important functions are the linear and angular velocity autocorrelation functions defined as:

6 m ->• •*

C^(t) = ^ <v{t).v(o)>

and

C (t) = | i <u(t) .u(0)>

From these functions the translational and rotational diffusion coefficient can be calculated

"x = i

°„=i

dt <v(t).v(0)>

dt <u)(t) .a)(0)>

For the translational degrees of freedom there is a simple connection between C (t) and F^(t) (Egelstaff, P.A., 1967) viz.:

1 *'^k<^' C (t) = em lim — — —

&t-For the rotational degrees of freedom no such a simple relation exists. It is well known that the velocity autocorrelation functions roughly behave as exponentials albeit with substantial deviations for intermediate and long times. In order to get an idea of the decay times T one may estimate for hard spheres these quantities from the initial decay of the velocity autocorrelation functions (Chandler, D., 1974)

^ = - ^ C . ( t ) T. 3t 1

1

1 = V , U)

t=0

For these decay times we find for our collision model

— = —

(1

- V ^ ^ ^ 'c°s

'^-D)

T^ Tg 2 ( K + 1 )

and _{i. _ 1 . l-cosyJ}

\ ^ ^ E 2(K+1)'

For smooth spheres (v' = 0) they reduce to

and

\

^E

i - = o

T

0)

The reciprocal Enskog relaxation time (l/x ) equals E

Here n is the density and g(a) is the pair correlation function in the contact point. For completely rough spheres (.f = ir) we have

^ = ^ ( 1 + - ^ ) and ^ = ^ ( - 1 ^ )

T T„ K + 1 T T K + 1 V E (1) E

References

11.1. Abramowitz, M. and Stegun, I.A. (1965), Handbook of Mathematical Functions.

11.2. Agrawal, A.K. and Yip, S. (1968) Phys. Rev. 171 (263).

11.3. Boutin, H. and Yip, S. (1968), Molecular Spectroscopy with Neutrons, M.I.T. Press, Cambridge, Massachusetts.

11.4. Chandler, D. (1974), J. Chem. Phys. 60, 3500 J. Chem. Phys. 6£, 3508.

11.5. Egelstaff, P.A. (1967), An Introduction to the Liquid State, Academic Press, London.

III. The Computer Experiment

To calculate the correlation functions mentioned in chapter II numerically, we have simulated the dynamics of a system of rough hard

spheres on a computer. From the phase space trajectories that one obtains by solving the equations of motion it is then possible to calculate the equilibrium time correlation functions. This so called molecular dynamics method is reviewed extensively in the literature

(Erpenbeek, J.J. and Wood, W.W., 1977; Wood, W.W., 1975) and therefore we will describe our computer experiment in a global way and only discuss some special features, we used in our calculations, in more detail.

Tracks of 10 collisions (400 mean free times) were generated for systems of 500 particles. Such a relative small number is sufficient if we use periodic boundary conditions (Alder, B.J. and Wainwright, T.E.,

1959; Wood, W.W., 1975) as has been shown by Alder et.al.(Alder, B.J., Gass, D.M. and Wainwright, T.E., 1970) and Wood (Wood, W.W., 1975) by checking the number dependence of the velocity autocorrelation function. The data are stored on magnetic tape and used as input for the calculation of the correlation functions. This computer experiment is performed for 8 densities viz. na"^ = 0.88, 0.79, 0.71, 0.47, 0.28, 0.20, 0.14, 0.10 and 5 values of the roughness parameter f (see chapter I) viz. ^ = 0,

* 3 N 3

IT/4, ir/2, 3TT/4, TT. Here n = na = ^ ci is a dimensionless density, N is the particle number, V is the volume of the box and a is the diameter of the sphere. For these 40 different situations the correlation functions are calculated and for the intermediate scattering function also for 6 values of the dimensionless wave number ka (see chapter II).

In addition calculations are performed to simulate the dynamics of neopentane (C(CH ) ) for the comparison with the neutron scattering experiment. This will be discussed in chapter IV.

111.1. Initial Conditions

Before starting the calculation of the time history of the system we have to generate an equilibrium situation for the particular values of the density (n ) and roughness if) for the system under consideration. First an initial configuration is set up with the spheres on the points of an f.c.c. lattice. The di^iensionless edge of the unit cell is

* 1/3

(4/n ) . Then for every direction a linear and angular dimensionless velocity is assigned to each sphere according to a gaussian distribution with zero mean and with a standard deviation of 1/3 for the linear velocity and 4/(3K) for the angular velocity. For the value of K in the

2

theoretical oriented simulations we have taken K = 4I/ma = 0 . 4 , corresponding to a homogeneous mass distribution over the sphere. This choice exaggerates the effect of the rotational degrees of freedom but for comparison with the literature (O'Dell, J. and Berne, B.J., 1975) K = 0.4 IS a convenient choice. For the neopentane oriented simulations we have as appropriate value K = 0.2. Our choice for the standard deviation means that we have made the velocities dimensionless with V =(/3kT/m).Thus the total dimensionless energy is twice the number of particles and the total dimensionless linear velocity is made zero. These two properties are maintained during the track generation m contrast with the total angular velocity which is made zero too at the beginning but obtains a finite fluctuating value after a certain time. This IS due to the periodic boundary conditions (Wood, W.W., 1975), because every time when a particle leaves the box and an image particle comes m , a sudden change in orbital angular momentum occurs. We have not observed any consequences from this phenomenon for the equilibrium properties of the system.

111.2. Equilibrium

and then check the translational order of the system (O'Dell, J. and Berne, B.J., 1975). This is done by calculating

1 N

P(t) = ^ - .1, [cos(kx.(t)) + cos(ky.(t) + cos(kz.(t))]

JN 1=1 1 1 1

where k = 4ir/a and a is the edge of the unit cube of the f.c.c. lattice. X.(t) is the position in the x direction of particle i at time t. At t = 0 we have P(0) = 1 , while in ecjuilibrium the value of P fluctuates around zero. For certainty we have run the system longer than strictly necessary and then it took 20 mean free times for low densities up to 140 mean free times for the highest density. After such long equilibration times we can safely assume that we have a Maxwell-Boltzmann distribution for the velocities.

III.3. Track Generation

To minimize the computation time we have divided the volume in a number of cubes and then we only look for collisions in the same or adjacent cubes. If between a certain pair of particles a collision can occur the collision time and the particle numbers are stored in a list unless the collision time exceeds a time T which is chosen so small that

max

it is impossible that a collision between a pair of spheres in non-adjacent cubes can occur within this time T . In this way we can not miss a

max

collision although we do not check every possible pair of particles. 3

Experimentally the optimum number of cubes appears to be 7 = 343 for high density and 8 = 5 1 2 for low density. So the number of pairs to

2 3 3 be checked for a collision is reduced from '5N(N-1) to ^N (—) or to

2 3 3

*5N (—) . These are reduction factors of about 0.08 and 0.05. Because the computer program itself is much more intricate, the track generation is only speeded up twice by this method.

After completion of this time table the first collision is performed and the time table including T is shortened by t,, the first collision

^ max •" 1

the list and the particles under consideration are checked for other collisions again m the way as described above. From this renewed list then the next collision is performed and so on until the list is exhausted. Then a new list is constructed. The initial and final velocities and

positions are stored and for every collision the particle numbers, the collision time and the new velocities are stored. Thus we can reconstruct the time history ot the system.

The computer time needed for generating tracks of 10 collisions runs from 60 minutes for high densities to 75 minutes for low densities on an I.B.M. 370/158.

III.4. Time Correlation Functions

The correlation functions

C,(t) = <A(t)A(0)> A

are calculated as (Berne, B.J., 1971)

1 1 N M m m C (kAt) = 7 7 - 2 , E, <A (t')A (tn+kAt)>

A N M n=l m=l n 0 n 0

N IS the number of particles and M is the number of correlation

functions calculated each starting at a different time t . k runs from 0 to K where KAt = T is the maximum time for which we are interested in the correlation function.

As we have stored only the velocities of the colliding particles, we have to calculate the molecular coordinates for every time kAt. This IS a very time consuming part of the computer program especially for the rotational part of the correlation functions. If we now choose the

m m-1

interval t|-, - t„ , for which we start the calculation of a new correlation function smaller than T we need to perform a substantial amount of

computer operations only once, while they are used m the calculation of several correlation functions viz. all those correlation functions

-*-

t

Figure 1. Sketch of the calculation of correlation functions

{see text).

with starting points t which fall in a region T (Bruin, C , 1974). For instance the distance covered by a particle in the time At (see

figure 1 ) needs to be calculated once while v.»e are using it for 2 3 4 the calculation of the correlation functions starting at t , t , t and t^.

In this way one can proceed and calculate different correlation functions at the same time. So we used three different computer programs for the calculation of the correlation functions, one for the linear and angular velocity autocorrelation functions, one for the orientation correlation functions and one for the scattering functions. These programs took roughly 10, 25, 90 minutes of computertime per track.

For the calculation of the relative orientation of the particles we used Caley-Klein parameters (Goldstein, H., 1959) instead of

direction-cosines as used by Berne (O'Dell, J. and Berne, B.J., 1975), because the change in the orientation has to be calculated from matrix multiplication, which is a very time consuming operation. The reduction from 3X3 matrices to 2X2 matrices reduces the computertime for the calculation of the correlation functions by 35%.

The functions L^ (z) (see chapter II), which are needed for the cal-2n

culation of the scattering functions are calculated separately for z = 0, ..., 400 with an increment of 0.01 and stored in a permanent dataset. In the correlation function program the value of the function L (z) was found by linear interpolation. This makes no difference for

the correlation function compared to the result obtained by direct calculation of L (z) and is certainly faster than direct calculation.

References

111.1. Alder, B.J. and Wainwright, T.E. (1959), J. Chem Phys. 3J_, 459. 111.2. Alder, B.J. Gass, D.M. and Wainwright, T.E. (1970), J. Chem Phys.

53_, 3813.

111.3. Berne, B.J. (1971) in Physical Chemistry, An Advanced Treatise, vol. VIII B, Ed. D. Henderson, Academic Press, New York.

111.4. Brum, C. (1974) Physica 72^, 261.

111.5. Erpenbeck, J.J. and Wood, W.W. (1977) in Modern Theoretical Chemistry, 6, Ed. B.J. Berne, Plenum Press, New York 1977. 111.6. Goldstein, H. (1959), Classical Mechanics.

111.7. O'Dell, J. and Berne, B.J. (1975), J. Chem. Phys. 63^/ 2376.

111.8. Wood, W.W. (1975) in Fundamental Problems in Statistical Mechanics III, Ed. E.G.D. Cohen, North-Holland Publishing Co. Amsterdam.

IV. The Molecular Dynamics Results

In this chapter we will discuss our molecular dynamics results. In the first three sections attention is paid to the intermediate scattering function and its approximations (see chapter II). In the fourth section we compare our results for neopentane (C(CH2) ) with a neutron scatter-ing experiment. The last three sections are concerned with the orientation correlation functions, the velocity correlation functions and an estimate of the accuracy of the computer experiment. In all sections the dependence of the correlation functions on the various parameters such as the density

(n ) , roughness (V) and the wave vector (k) will be discussed.

IV.1. The Decoupling Assumption for F (t)

The decoupling assumption for the intermediate scattering function

F^(t) « Fj^(t)F^(t) (1)

T R into a translational (F ) and a rotational (F ) part is practically always made in the literature describing neutron- and light scattering experiments. Only recently some authors (Montgomery, J.A. and Berne, B.J., 1977; Wolyness, P.G. and Deutch, J.M., 1977) have been concerned with the mutual influence of translational and rotational degrees of freedom in connection with the calculation via a hydrodynamic approach of the diffusion tensor for a molecular fluid. Thus it is interesting to check the validity of this approximation (equation 1) by a computer experiment, where the 3 functions involved ccin be calculated separately. In chapter II it was shown that the approximation (equation 1) consists of two parts viz. a series truncation and an assun^ition of statistical independence of the translational cUid rotational scattering function.

IV. 1, a. Series Truncation

In chapter II (equation 3) it is shown how the intermediate scattering function can be written as

F^(t) = ?„(-l)" (4n+l)F"(t) (2) k n=0 k

with

F;;(t) = <J2„(kd)P2^(d.i)L2jz)>

To see to what extent the scattering function can be approximated by the first term in this expansion both the first and second correlation function in this series are calculated. The result for F (t) is shown in table I (appendix A) for several values of the density, the wave vector and the roughness parameter. To visualize the situation we have plotted F (t) and F (t) in figure (1) for a representative set of para-meter values. On the basis of this table and figure (1) we conclude that to a good accuracy we can approximate the scattering function by the first term in the expeinsion.

IV. l.b. Product Assumption

The next assumption to study is the decomposition

FO«F^(t)F«(t)

with (equation II.4, 5)

F^(t) = <JQ(kd)>

F«(t) =<L^(Z)>

Figure 1.

The scattering functions F.,(t)

(1) and -SF^(t) (2) versus

the time in units of the mean

free time (i). The density

in ) IS 0.88, the roughness

('fi) is TT and the wave vector

(ka) is 4.0.

0 8 0 6 0 4 0 2 0 0

\

A

- \

\

-^^"^^^^^r—-J_

1 1 1 1 1 3 , 1

Figure 2.

As figure 1 for the

scatter-ing functions Fv(t) (1),

^(t) (2) and their

differ-ence (3) .

F°(t) - F^(t) Fj(t)

The results are listed in table II (appendix A) for the same parameter values as in table I. Inspection of table II shows that this assumption is valid within the accuracy (see section 7) of our experiment.

The final conclusion of this section is that the total scattering function can be approximated as (equation 1)

F^(t) «= F^(t) F^(t) (1) a l t h o u g h t h e r e i s a n a d d i t i o n a l c o n t r i b u t i o n (up t o 0 . 0 2 ) d u e t o F ( t ) f o r i n t e r m e d i a t e a n d h i g h k - v a l u e s . T h i s r e s u l t i s n o t u n e x p e c t e d b e c a u s e s p h e r e s i n c o n t r a s t t o s p h e r o c y l i n d e r s f o r e x a m p l e d o n o t h a v e a p r e f e r r e d d i r e c t i o n o f r o t a t i o n . T h e c o m p a r i s o n o f t h e s e r e s u l t s w i t h t h e o r y s o f a r d e v e l o p e d ( B e r n e , B . J . a n d M o n t g o m e r y , J . A . , 1 9 7 6 ) s e e m s n o t v e r y m e a n i n g f u l b e c a u s e t h e i r t h e o r e t i c a l p r e d i c t i o n ( f r o m t h e J d i f f u s i o n m o d e l , s e e c h a p t e r V) d i f f e r s a l r e a d y more f r o m t h e n u m e r i c a l r e s u l t s t h a n t h e c o n t r i b u t i o n o f F, ( t ) t o t h e s c a t t e r i n g f u n c t i o n . k I V . 2 . T h e G a u s s i a n A p p r o x i m a t i o n As shown i n s e c t i o n I I . 4 . t h e g a u s s i a n a p p r o x i m a t i o n r e a d s a s F ^ t ) = F ^ ^ ( t ) F ^ ^ ( t ) k k k w i t h TG 2 2 F j ^ ^ ( t ) = e x p ( k < ( r ( t ) - r ( 0 ) ) > / 6 ) 2

F]^^(t) = exp - ^^~- (l-<u(t).G(0)>)

This already includes the decoupling approximation, so we can study the

translational and rotational part separately.

IV. 2. a. The Translational Part

The gaussian approximation for the translational scattering function

F^(t) . F^-.t)

IS a widely used approximation which has been studied carefully by Rahman (1964) and by Nijboer and Rahman (1966) for atoms interacting via a Lennard-Jones potential. For hard spheres Wood (1975) has given the contribution of the next term in the cumulant expansion (of which the first term gives the gaussian approximation). To provide some more numerical results we have given in table III (appendix A) the difference

T TG

(F - F ) . This shows clearly that the gaussian approximation is satis-factory. Only for high k-values at the highest density is there a sig-nificant deviation. The rotational degrees of freedom seem to have no influence on this difference.

IV. 2.b. The Rotational Part

For the rotational part of the scattering function the gaussian approximation

Pf (t) « Ff^(t)

k k

has been studied by Agrawal and Yip (1958) (see also Boutin, H. and Yip, S., 1958) utilizing free rotator correlation functions, but only for small k-values. As is shown in table IV (appendix A) and figure (2)

R RC

the difference F (t)-F (t) is larger than the corresponding quantity for the translational part and it shows up for lower k-values. This can be understooc

and Table I)

R RC

the difference F (t)-F (t) is larger than the corresponding quantity lower k-va

R RC

FJ^(t ^- ") = j2(ka/2)

F^^(t ->•")= exp-(ka) 2/12

For low k-values (ka < 2.0) there i s no difference between these two

expressions and the agreement for a l l values of the density and the

roughness parameter i s good. For the higher k-values (ka=4.0, 8.0) the

R RG

agreement i s poor because the difference between F and F for long

times i s much l a r g e r (-0.057 and 0.031 r e s p e c t i v e l y ) .

Also for s h o r t times and e s p e c i a l l y for low d e n s i t i e s t h e r e i s a

considerable difference (of the order 0 . 1 0 ) . This finds i t s o r i g i n i n

the poor r e p r e s e n t a t i o n of the free r o t a t o r function by the gaussian

approximation. That the r o t a t i o n a l s c a t t e r i n g function for s h o r t times

(up t o 0.7 mean free times) behaves l i k e a free r o t a t o r function can be

seen from figure (2a) (appendix B). I f we make a short time expansion

R RG

for F, (t) and F, (t) in the free r o t a t o r l i m i t , the difference between

_{k k .}

4

these two series appears in the t term because the fourth cumulant 4

starts with a t behaviour.

4

F^(t) « F^'^ exp -— {<(u(t)-u(0))'*> - I <(u(t)-u(0))2>2}

\ <i "^ l o ' H T ^ B T ' • • • '

Tj TIC* 4

So the difference between F (t) - F (t) s t a r t s with a t te'-m

2 2

L. , (ka) t 2

10 * 12 Bl'

4

which is already considerable compared with the t term from the gaussicui

+ L ,(ka)2 2 t2 2

|-(l/(l + 10/(ka)2))

which amounts to 0.12 and 0.17 for ka=4.0 and 8.0 respectively. So the conclusion is that one should certainly not use the gaussian approximation for molecular gases.

IV.3. Influence of the Roughness Parameter ^

In section (1) of this chapter we have already shown that for all practical purposes the total intermediate scattering function can be approximated as

F^(t) «=F°(t) =« F^(t)Fj^(t) ' • ••• '•"' • - •• • '• •- (1)

To understand the influence of the roughness parameter ^ we will

0 T R '•..•-discuss in this section the functions F, (t) , F, (t) and F, (t) for 3

* k k k different densities viz. n = 0,10; 0,47; 0,88 as functions of the rough-ness. The functions are plotted for the most significant values of ka and 5 values of the roughness parameter viz. • - •• ^ '

¥> = (0,1,2,3,4) X 71/4 "'• " ". • • - ' • '

The smooth sphere curves (v' = 0) are easily recognized because the

p

asynptotic behaviour of F (t) for ka ^ 2.0 differs considerably from the rough sphere curves (table I) . The wiggles in the v'=0 case for longer times are due to poor statistics. Since the angular velocities of the particles do not change by collisions, the average is only over the initial velocities i.e. 500 different values. As soon as the angular velocities change in the course of time the statistics are improved. In the following we discuss the roughness influence for the various densities.

ka 0.5 1.0 2.0 4.0 8.0 F^^t) 0.979 0.919 0.708 0.207 0.036

-r<"

0.979 0.920 0.717 0.264 0.005 F (t) (free rotator) 0.986 0.946 0.803 0.440 0.197

Table I. Asymptotic behaviour of the rotational

scattering function, gaussian approximation

and the free rotator.

High density (n = 0.88)

Here the i n t e r e s t i n g range of k-values i s 2.0 < ka < 8.0. The

influence of the roughness parameter i s large as can be seen from

figure (3) . At the value of f = T^/2 for the roughness parameter

the curves are the most s e n s i t i v e to a change in ip, while between

the V' = 3ir/4 and v = TT curves there i s l i t t l e d i f f e r e n c e . This

holds for the t r e i n s l a t i o n a l as well as for t h e r o t a t i o n a l s c a t t e r

-ing function. All the curves for F (t) are monotonically decreas-ing

functions, no free r o t a t o r behaviour appears a t t h i s d e n s i t y .

Intermediate d e n s i t y (n = 0.47)

The range of r e l e v a n t k-values i s 1.0 < ka < 4 . 0 . The influence

of the roughness parameter i s much smaller (figure 1 in appendix B),

but the type of dependence i s the same. From t h i s figure one can

Figure 3a.

The -scattering function FV(t)

for 3 Values of the wave

vector (ka). The numbers

ly 2 ...S indicate the

rough-ness, being (i-l)'X-n/4. The

unumbered curves should be

read in the same order. The

density is n = 0.88.

10 20 30 40 50 60 1.0 0 8 -06 0 4 0.2 0 0

•\xr'*'*^^^^2==-jii^2

\ ^$^^^=::i-!l'^=^

\^ ^"\r\r~^-^

\ ^ ^ ^ \ ^

\ 0 ^ : v

X X ^ ^ ^ - ^ ^ ^ ^ l < cr= 8 ^ ^ ~ ~ ~ ~ ^ ^ ^ ^ ^ ^ ^ ^ ^

^ ^ ^ ^

^ - ^

>

^

1

^ ^ =

1 L I 1 1 1

Figure 3b.

As figure 3a for F^(t),

0 10 20 30 40 50 60

Figure 3c.

As figure 3a for F (t)

also see that there is already free rotator behaviour; this of course is most evident for the ^ = ir/4 curves.

Low density (n = 0.10)

For this density we have plotted the curves (figure 2 in appendix B) for the ka-values 0.5, 1.0 and 2.0. Again the influence of the roughness parameter is small. The dependence on f of F (t) is almost totally due to the translational part of the scattering function, which for this density can be described quite accurately by the gaussian approximation. Thus the influence of f on F (t) can be predicted from the behaviour of the velocity correlation function (see section 6 ) . For this density all curves (for F ) behave like free rotators at least for short times and then fall off slowly to their long time limit

(j2(ka/2)).

In conclusion we can summarize this section with the statement that although there is roughness dependence for low and intermediate density, the most significant effects occur at high density.

IV.4. Comparison with a Neutron Scattering Experiment

In order to check our collision model against reality we compared the molecular dynamics results with a neutron scattering experiment on neopentane (C(CH)) performed by Chr. Steenbergen at I.R.I. Delft (to be published). Neopentane is a molecule with tetrahedral symmetry (one central carbon atom with methyl (CH ) groups at the corners), a so called spherical top for which the moment of inertia tensor can be re-placed by a scalar. Hence for such molecules we expect that a rough sphere description is possible.

The neutron scattering experiment is performed at liquid density (under 1 atm. pressure) at T = 279 K (p = 0.606 gr/cmi ). Besides the corrections for background scattering and the resolution of the measuring apparatus, these measurements are also corrected for the contribution of the internal degrees of freedom, which could be done because the peaks in S(k,w) due to the combined translational-rotational motion and the one caused by the internal motion are quite well separated.

The accuracy varies from ±0.04 (for 0.5-2 psec.) to ±0.02 (for 2-4 p s e c ) . To compare the experimental results with the computer ejcperiment we have to know the equivalent hard sphere diameter for neopentane. It is taken to be (van Loeff, J.J., to be published)

a(hard sphere) = 5.60A

Thus the values of the density (n ) and the loading parameter (K) for our experiment are

n =0.89

K = - ^ = 0.20 ma

To calculate the intermediate scattering function from this computer simulation, we have to correct for the fact that the protons (which

0 8 0 6 0 4 0 2 0 0 ' •< Y . T r

-io

20 30 40

Figure 4a.

The scattering function

F^(t) for completely rough

spheres (f - it, solid lines)

compared with the neutron

scattering data (marks).

The values of the wave are

k - .23, .44, .50, .60, .73

and .S3^ from top to

hottom. The upper horizontal

axis inaiaates the time in

picoseconds,

psec

0 0

"• Igure -'b. ' - "*

4s figu^e 4a for half rough

are the neutron scatterers) are not located on the surface of the sphere. For the distance between the central carbon atom and the hydrogen atoms we have taken (Bierdermann, S., 1972)

C -H = 2.16^ c

With these choices for the various parameters we have calculated the intermediate scattering function as

F^(t) « F°(t)-5 X Fl(t) f o r 3 d i f f e r e n t v a l u e s of t h e roughness p a r a m e t e r v i z . i^ = ( 2 , 3 , 4 ) X Tr/4 and f o r 6 v a l u e s of t h e wave v e c t o r v i z . k = . 2 3 , . 4 4 , . 5 0 , . 6 0 , .73 and . 9 3 S " T h i s c o n ^ u t e r s i m u l a t i o n h a s a l s o been performed f o r 10 c o l l i s i o n s , y i e l d i n g v e r y good a c c u r a c y (see s e c t i o n 7 ) .

The r e s u l t s of t h e s e two t y p e s of e x p e r i m e n t s a r e compared i n f i g u r e ( 4 ) , from which i t can be s e e n t h a t t h e i n f l u e n c e of t h e

rough-2

n e s s p a r a m e t e r i s d e c r e a s e d f o r a s m a l l e r v a l u e of K = 4l/ma . Never-t h e l e s s a r e l a Never-t i v e l y l a r g e e f f e c Never-t , e s p e c i a l l y f o r Never-t h e h i g h e r wave numbers, e x i s t s . We have n o t shown t h e ip = 3ir/4 c a s e b e c a u s e , as a l r e a d y remarked i n s e c t i o n (3) t h e r e i s l i t t l e d i f f e r e n c e between t h e tp = -n and f = 3Tr/4 c u r v e s .

I n t h e c o m p l e t e l y rough sphere c a s e t h e h i g h e s t k - v a l u e ( . 9 3 A ) f i t s t h e t h e e x p e r i m e n t v e r y w e l l , b u t a l l t h e o t h e r k - v a l u e s give t o o h i g h r e s u l t s f o r t h e m o l e c u l a r dynamics c u r v e s . For h a l f r o u g h s p h e r e s t h e agreement i s much b e t t e r a l t h o u g h f o r t i m e s between 1 a n d 2 p i c o -s e c o n d -s t h e m o l e c u l a r dynamic-s r e -s u l t -s -s t i l l d i f f e r t o o much from t h e r e a l e x p e r i m e n t . S i n c e t h i s i s t h e r e g i o n where t h e i n f l u e n c e of t h e r o u g h n e s s p a r a m e t e r i s t h e l a r g e s t and as we have seen i n t h e f o r e g o i n g s e c t i o n Fj^(t) i s t h e most s e n s a t i v e t o changes i n t h e r o u g h n e s s a t a b o u t

somewhat higher value of the roughness parameter ^.

IV.5. The Orientation Correlation Functions

As shown in chapter II the rotational scattering function can be expressed in terms of the orientation correlation functions C (t)

F^(t) =^(2^+l)3,(-)C,(t)

Cj^(t) = <Pj^(G(t) .u(0))>

Hence the time dependence of F (t) is completely expressible in terms of the orientation correlation functions. In practice the summation up to Jl=3 suffices. We have calculated C (t) and C (t) also because of their interest for infrared spectroscopy and Raman scattering. Of these functions C, (t) is the most interesting function because, due to its larger numerical value, the roughness dependence is more pronounced. For the completely rough sphere these functions have already been calculated by C Dell and Berne (1975), although less accurately, from

4

a computer simulation with 10 collisions for 108 particles.

We have plotted the C (t) for £ = 1,2 and 3 (figure 5 and appendix B ) . We have included the )l=3 term to see whether it is possible to predict the behaviour of the rotational scattering function if one knows C.(t) and C (t) already. This is possible for ka < 3.0 as can be seen from the

^ 2 numerical value of the weight factor (7 X 3 (3.0/2) = 0.0056). For

higher values of the wave vector and for relatively short times it is not possible because of the contribution of C (t). As the decay times of the orientation correlation functions C , C and C differ considerably

* R for high density (n = 0.88) we can predict F (t) from knowledge of C (t)

K 1

and C (t) for longer times. This is in contrast to the low density

^

*

case (n = 0.10) where the decay times are almost the same.

Again we can immediately recognize the smooth sphere (1^=0) curves as having a finite asymptotic value