# Advances in chemical physics - Prigogine I.

**Download**(direct link)

**:**

**21**> 22 23 24 25 26 27 .. 51 >> Next

(4) To test various stochastic assumptions for molecular motion that would simplify the N-body problem if they were valid. Molecular dynamics is far superior to experiment for this purpose since it provides much more detailed information on molecular motion than is provided by any experiment or group of experiments.

These computer experiments have provided insight into the microscopic dynamical behavior of real diatomic liquids for both the experimentalist

122

В. J. BERNE AND G. D. HARP

and theoretician alike. Further it is hoped that these studies will motivate more realistic approximate theories of the liquid state and provide “ experimental” data to test these theories against.

The molecular dynamics calculations were carried out in a manner similar to that used by Rahman in his original study of liquid argon.32 A finite number of molecules, N, were assumed to interact pairwise through a given truncated intermolecular pair potential. In addition, the atoms on the same molecule were assumed to interact through a harmonic potential, \Kv{ri — r)2, where Kv is the ground state vibrational force constant for the molecule, rt is the internuclear separation for the rth molecule, and f is the ground state equilibrium internuclear separation. For carbon monoxide, Kv = 1.9020 x 106 dynes/cm and r = 1.1281 A. For nitrogen, Kv = 2.2962 x 106 dynes/cm and r= 1.094 A.45 The harmonic potential was added because the calculations were done in the cartesian coordinates of the atoms forming the molecules. These atoms were originally separated by the equilibrium internuclear distance. They remained separated by this distance to within ~10-4 A throughout the course of the calculations. Therefore the results of these computations are essentially those for systems of rigid rotors.

The center of mass of each molecule was initially placed in a cubic lattice system within a large cube. The length, L, of a side of the cube was (NMjp)1/3 where M was the mass of a molecule and p was the density of the fluid. L was typically 30 A in these calculations. The molecular orientation angles were chosen randomly on a unit sphere. That is, if 0 and ф were the usual molecular axis, polar orientation angles, then ф was assumed to be uniformly distributed between 0 and 2n and cos 0 was assumed to be uniformly distributed between —1 and 1. The relative and center of mass velocity components were chosen by the Von Neuman46 rejection method from Gaussian distributions appropriate to a gas of rigid rotors at some preselected temperature, T. That is, if Vx, Vy, and Vz were the center of mass velocity components and 6 and ф were angular velocities, then the probability of selecting these velocities was given by

B. Method Employed

ON THE CALCULATION OF TIME CORRELATION FUNCTIONS

123

The selected relative and center of mass positions and velocities for each molecule were then transformed to cartesian coordinates to give an initial set of velocities and positions to each of the two atoms making up the molecule. Periodic boundary conditions were imposed.31,32 That is, if (*, y, z) were the coordinates of a particle in the original cube, then there were assumed to be 26 images of this particle with coordinates gotten by adding and subtracting L from each cartesian coordinate of the original particle. A given particle i then interacted not only with every particle j within the original cube but also with all of particle j’s images. Further, if during the course of the calculation a particle passed outside the original cube, then it was replaced by a particle entering the opposite side of the cube and having the same velocity as the particle that left. In other words, the number of molecules in the cube was constant. Hamilton’s equations of motion for the N molecules were then solved numerically using the Runga-Kutta-Gill method with a stepsize in time, At, of 5 x 10"15 s. See Appendix A for a general discussion of the reasons why this numerical method47 and time step were used.

During the course of the calculations the translational and rotational temperatures, TT and TR, respectively, were monitored at each step. These temperatures were defined by the equipartition theorem:

where V{ is the center of mass velocity of the /th molecule,V/ is the velocity of the jth atom on the /th molecule, and Mt and M2 are the masses of the two atoms making up the molecule. The formula for TR assumes explicitly that we were dealing with systems of rigid rotors. This was actually a very good assumption since the vibrational coordinates only contributed ~0.02°K to TR and the variance of TR due to statistical fluctuations was typically 1000 times larger than this contribution. The total kinetic temperature, TK, was then defined by

The initial distribution of positions and velocities was not that of an equilibrium fluid at the preselected temperature T. Therefore, during the initial or equilibration phase of these calculations the monitored temperatures fluctuated wildly. This behavior is illustrated in Figure 1 where TT

**21**> 22 23 24 25 26 27 .. 51 >> Next