- Email: [email protected]
Non-Maxwell distributions in equilibrated fluids. III
Volume 170, number 4
13 July 1990
CHEMICAL PHYSICS LETTERS
Non-Maxwell distributions in equilibrated fluids. III S.-B. Zhu, J. Lee and G.W. Robinson SubPicosecondandQuantum Radialion Laboratory, P.0. Box 4260, Texas Tech University,Lubbock, TX 79409, USA Received 27 March 1990
We examine by molecular dynamics methods the equilibrium velocity distribution of a finite sized system without using periodic boundary conditions. Extremely long runs with very short integration time steps were performed. The system is two dimensional with 196 host particles and 4 small test particles of low mass, the latter being subjected to an intense two dimensional corrugated force field. The smooth nonMaxwellian distributions obtained for the test particles show a remarkable consistency with those in papers I and I1 of this series.
1. Introduction In a recent series of molecular dynamics (MD) [ l-71, we have discovered that the velocity distribution function of test particles, which are small both in mass and size compared with those of heat bath particles and are concomitantly subjected to steep potential fields, may deviate from Maxwellian significantly. We believe this is a fundamental feature of ultrafast dynamical processes where the nonlinear properties of friction are of importance [ 8,9]. Though our investigations have covered quite a broad range, including liquid and crystalline lattice heat baths and test particles with internal degrees of freedom or those subjected to strong external forces, the results are always consistent. The equilibrium velocity distribution of the test particles is always nonMaxwellian. Theoretical support for non-Maxwellian velocity distributions will be supplied in a future paper [9]. There is also accumulating experimental evidence indicating a strange behavior of hydrogen atoms in certain environments [ 10-l 61. Such experiments seem closely related to our findings. Our results have evoked much controversy, summarized by the following three questions: ( 1) Is the integration time step appropriately chosen so that it does not induce numerical errors? (2 ) Does the periodic boundary condition affect the nature of the distribution function? (3) Has equilibrium really
been established? Here we provide further evidence in support of our previous work.
studies
368
2. Finite-sized system A simple way to test the influence of periodic boundary conditions is to consider a sample with finite volume. In order to avoid possible edge effects, while still using a computational effort of realistic length, we study a two-dimensional system containing 196 host particles and 4 test particles. These particles interact with one another through LennardJones (LJ) 12-6 potentials. The LJ parameters are given in table 1. To create the necessary nonlinear dynamics, the four test particles are further subjected to an X, y corrugated potential surface having the form 4
UC=& 1 (cos&,+sinKy,), i=, where K= 56.3/u,,,, aM being the W diameter of a host particle; X, and yi designate the coordinates of the ith test particle, taking the origin to be the center of the sample box. The reduced barrier height h is ZU,/knT. In addition, a two-dimensional external potential
u,= ; A,(xy+
Y:')
t=I
Elsevier Science Publishers B.V. (North-Holland)
Volume 170,number 4
13July 1990
CHEMICALPHYSICSLETTERS
Table 1 Potential parameters
2-D c’ IIa d, IIb *’ VII d’ VIII d’ a’
%M=
0.579 0.579 0.579 0.579 0.579
0.26 0.26 0.26 0.26 0.26
0.9714 0.9714 0.9714 0.9714 0.9714
4 4 4 4 4
0.366 0.366 0.366 5.856 0.023
0.0005 0.005 0.0025 0.005 0.005
16.92 3.60 13.01 3.84 4.80
:(%+u,wM); htM=JcG.
b, N is the total number of timesteps in millions; but note that there are 8 x more data per timestep in the present computation than in ref. [6]. ‘) Present work: two-dimensional system, two-dimensional corrugated potential, d, Four fully independent runs from ref. [6]: three-dimensional system, one-dimensional corrugated potential. The distribution functions for eachof thesesystemsare depictedindividuallyin ref. [ 61. The four data sets have been averaged in this paper and plotted (circles) in fig. 2 with their standard deviations.
is applied to every particle to diminish its probability of evaporation. In eq. (2) Xi=X,/6.7o,w and Y,=y,/6.7aM represent the reduced coordinates of the ith particle. The choice of the form and magnitude of this potential is quite arbitrary. We choose A,= 16 for a test particle and AM= 1 for a host particle. Using the parameters listed in table 1, we are able to obtain a uniformly distributed system confined within a square of 13.4x 13.4&. Throughout the paper dimensions are expressed in reduced units, taking A4 and the LJ parameters u,,, and tM as measures of mass, length and energy, respectively. The equations of motion are integrated by using a modified version of the Verlet algorithm [ 171. The initial velocities of the particles, obeying a Gaussian distribution, are selected from a random number generator. In order to prevent possible numerical errors caused from an inappropriate choice of integration time step, we adopt here a computationally achievable small time step, At=O.O005r where T=u,(M/e,)“2, i.e. only one fifth to one tenth of that employed in previous work on similar systems. See ref. [ 61 and table 1.
3. Statistics A fully equilibrated system can be prepared by a preliminary integration of the equations of motion using 104At. Information about the numerical ac-
curacy of the results and the approach to equilibrium can be obtained by examining fluctuations of the total energy of the system and fluctuations in the mean square velocity of the test and host particles in the aftermath of this equilibration procedure. These fluctuations all lie within 0.05%. In order to answer the question whether the system is really in equilibrium and whether sufficient statistics have been accumulated, we perform an extremely long run. Intermediate results are extracted after every 1.2X 106 time steps starting from t = 6.12 x 106At. These results can be compared with each other to make certain that a steady-state distribution is actually reached. Altogether, 1.692~ 10’ time steps are carried out to accumulate data. This is at least three to four orders of magnitude longer than the time required for the system, initially far from equilibrium, to reach an equilibrium state. In view of the fact that four test particles are investigated and that two velocity components are taken into account, the quality of the statistics, compared with earlier results [4,6], in increased by a factor of eight. Utilization of the two-dimensional corrugated potential surface also avoids any possible superperiod problem of one-dimensional nonlinear weighted strings of the type first studied by Fermi, Pasta and Ulam [ 18,191. Our previous calculations [4,6] on three-dimensional systems with a one-dimensional corrugated potential were criticized as being “too one-dimensional” to acquire a sensibly rapid approach to equilibrium. The results to be re369
Volume 170,
number4
ported below, showing nearly identical distribution functions from the two types of systems when the parameters (m/M, aJaM, h, T) are the same, fully negates this criticism.
Graphically illustrated in fig. 1 are the distribution functions of the velocity components of the test particles (asterisks) and those host particles (circles) whose distance from one of the test particles is less than 1.5a,. We call these host particles the “cage particles”. Because of the dynamical motion, the cage includes different host particles at different instants of time. The curve symbolized by asterisks is averaged from 1.3536~ 1OSsamples, while the circles are from 6.0066~ 10’ samples. For comparison, the Maxwell velocity distribution is also plotted in fig. 1. The nearby host particles are seen to exhibit a nearly Maxwellian behavior, indicating that these particles are fairly well mixed with other host particles. In agreement with our previous work, we see from fig. 1 an unmistakable noncanonical behavior of the low-mass test particles. In fact, as emphasized in fig. 2, the distribution function is nearly identical to the ones obtained earlier for a thermodynamically sim-
Velocity
Fig. 1. Distributionsof the velocitycomponents.Solid curve: Maxwell distribution; circles: distribution function for the “cage particles”, asterisks: distribution function for the test particles. The reduced velocity here and in fig. 2 is defined as V= [ $1 (kJ/ P) P2 , where u is the velocity component for the test or host particles and p is the respective mass, m or M. The verticle axis is
scaledby percentage. 370
2.7
c2.1
0 .c ii ‘i
Gl.5 .n
4. Results
Reduced
13July 1990
CHEMICALPHYSICSLETTERS
0.9,
-
-
\.* 0.0
I
I
I
,
0.5
Reduced
I
I
I
,
1 .o
,
,
,
,
1.5
Velocity
Fig. 2. Distributions of the velocity components. Solid curve: Maxwell distribution; asterisks: distribution function for the test particles (this work); circles: average of distribution functions for systems IIa, IIb, VII and VIII of ref. [ 61 (see also table 1). The error bars are standard deviations among these four sets of data, assuming the distribution functions are independent oft,/ CM.
ilar three-dimensional system with a one-dimensional corrugated potential. The observations here, combined with those of ref. [ 61, provide evidence that (1) that distribution functions from the four ji& independent runs in ref. [ 61 average, within their standard deviations, to the same non-Maxwellian distribution function as obtained here, though the two-dimensional system needs longer time to reach equilibrium; (2) the use of periodic boundary conditions does not appreciably change the nature of the average dynamics of the test particle; (3 ) the time step employed in the previous work is not too wide; and (4) the distribution function is insensitive to the test-host interaction strength \/E,tM. An important criterion of an equilibrated system is the average kinetic energy. The present calculation gives (KM) =0.99999kBT for host particles and (Km) = 1.00027kJ for test particles, where k, indicates the Boltzmann constant, and /r,i”/~,= 0.97 14 represents the reduced temperature of the heat bath. Thus, the entire system must be in thermal equilibrium and the numerical accuracy ofthe computation is satisfactory. To ensure that the observed distribution really represents a steady-state quantity, we plot in fig. 3 the normalized probabilities, within a small range ( V, V+ A V) of velocity space, for test particle velocities as a function of the number of samples used in the statistical average. Here, V defines the reduced ve-
Volume 170,number 4
CHEMICALPHYSICSLETTERS
6z 0 .-Y 0 .-
: o_
2 0
: -6-
* *
-12
40
*
*
*
* * 8 I I I I I I I I I 8, I I 60 “if %?ples
*
* I 08 120
* 140
No.
Fig. 3. Normalized probability deviations Y( V) of test particles within a small range (I’, T/t AV) of velocity space. Circles: Y’(0); asterisks: ul( 1.19); Ac’=O.O3125.The x-axis is to be multiplied by I 06,the y-axis by 1O-‘.
locity component normalized by the one-dimensional mean square velocity Jk,Tlm of the test particles. Using PO(V) to represent the probability calculated from the Maxwell velocity distribution, we define Y( V) = 1 -P( V)/Po( V) to measure the deviations of the distribution from Maxwellian. The fluctuations of this quantity throughout the final 10.8 million time steps are seen to be small compared with their magnitudes. Moreover, as the number of samples increases, this quantity has tended to a constant that differs significantly from zero. Therefore, the distribution is stable, and the system is in a steady state. These results reply to the criticism that the system has not reached an equilibrium state or that the simulation time is not sufficiently long.
5. Conclusions The present paper is aimed at answering a number of questions raised concerning previous investigations. In order to get rid of the conventional periodic boundary conditions, we have studied here a system with finite size. The doubt about the length of the time step is dismissed by reducing At by a factor of 5-10 (see table 1). In doing so, we have obtained a distribution function nearly identical to the average (fig. 2) of that obtained from four fully independent runs for thermodynamically similar systems (same m Jh4, u,,,/uM, h, T), where the At, the dimensional-
13July 1990
ity and the boundary conditions were radically different. The results shown in fig. 2 thus indicate that the velocity distribution functions, though decidedly non-Maxwellian, are independent of dimensionality. The results also suggest that they are insensitive to the test-host interaction energy CM. To obtain better statistics, we have performed a very long molecular dynamics simulation. In addition, four test particles and two velocity components were considered, further improving the statistics by a factor of eight. The smoother distribution function obtained here is undoubtably a result of these improvements, but may also be partly a result of using the shorter timestep. Fluctuations of the probabilities that test particles remain within a small range of velocity space during the simulation time have also been investigated. The smallness of these fluctuations and the approach of the probabilities to fixed limiting values provide evidence that the system has reached its equilibrium state on a time scale much shorter than that employed to accumulate data. Additional evidence for the system being in equilibrium is given by the fact that the average kinetic energies of the test and host particles are equal and the total energy of the system is conserved to within 0.05% throughout the run. In view of the above results and discussion, we conclude that the observed noncanonical behavior is a real characteristic of equilibrated systems in which ultrafast dynamical processes are taking place. It is not a consequence of computationally introduced errors or poor statistics [20]. The effect of these distorted velocity distributions on velocity-dependent cross sections in chemistry and physics will radically change current views about large classes of thermally activated condensed phase reaction processes.
Acknowledgement Financial support at the SPQR Laboratory been shared by the Robert A. Welch Foundation 0005, 30% and D- 1094, 5%)) the State of Texas vanced Research Program (1306, 47%) and Pittsburgh Supercomputing Center ( 18%).
has (DAdthe
371
Volume 170, number 4
CHEMICAL PHYSICS LETTERS
References [ 11 S.-B. Zhu, J. Lee and G.W. Robinson, I. Chem. Phys. 88 (1988) 7088. [2] S.-B. Zhu and G.W. Robinson, Chem. Phys. Letters 153 (1988) 539. [ 31 S.-B. Zhu and G.W. Robinson, J. Phys. Chem. 93 (1989)
164. [4] S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 161 (1989) 249. [ 51 S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 163 (1989) 328. [ 61 S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 169 (1990) 355. [7] S.-B. Zhu, J. Lee and G.W. Robinson, I. Phys. Chem., submitted for publication. [ 81 S.-B. Zhu, J. Lee and G.W. Robinson, Velocity Dependence of Friction, Chem. Phys., submitted for publication. [ 91 S.-B. Zhu, Velocity Distributions in Nonlinear Systems, Phys. Rev. A, submitted for publication.
372
13 July 1990
[ IO] G. Comsa, R. David and B.-J. Schumacher, Surface Sci. 95 (1980) L210. [I l] G. Comsa and R. David, Surface Sci. 117 ( 1982) 77. [ 121 S. Estreicher, Phys. Rev. B 36 ( 1987) 9 122. [ 131 S.J. Pearton, J.W. Corbett and T.S. Shi, Appl. Phys. A 43 (1987) 153. [ 141 M. Fleischmann, S. Pons and M. Hawkins, J. Electroanal. Chem. 261 (1989) 301. [ 151 SE. Jones, E.P. Palmer, J.B. Czirr, D.L. Decker, G.L. Jensen, J.M. Thome, S.F. Taylor and J. Rafelski, Nature 338 (1989) 737. [ 161 R.J. Beuhler, G. Friedlander and L. Friedman, Phys. Rev. Letters 63 (1989) 1292. [ 171 W.C. Swope, H.C. Anderson, P.H. Berens and K.R. Wilson, J. Chem. Phys. 76 (1989) 637. [ 181 E. Fermi, J.R. Pasta and SM. Ulam, in: Collected works of Enrico Fermi, Vol. 2 (University of Chicago Press, Chicago, 1965). [ 191 J.I. Tuck and M.T. Menzel, Advan. Math. 9 ( 1972) 399. [ 201 W.P. Keirsteadand K.R. Wilson, J. Phys. Chem. 94 ( 1990) 918.
13 July 1990
CHEMICAL PHYSICS LETTERS
Non-Maxwell distributions in equilibrated fluids. III S.-B. Zhu, J. Lee and G.W. Robinson SubPicosecondandQuantum Radialion Laboratory, P.0. Box 4260, Texas Tech University,Lubbock, TX 79409, USA Received 27 March 1990
We examine by molecular dynamics methods the equilibrium velocity distribution of a finite sized system without using periodic boundary conditions. Extremely long runs with very short integration time steps were performed. The system is two dimensional with 196 host particles and 4 small test particles of low mass, the latter being subjected to an intense two dimensional corrugated force field. The smooth nonMaxwellian distributions obtained for the test particles show a remarkable consistency with those in papers I and I1 of this series.
1. Introduction In a recent series of molecular dynamics (MD) [ l-71, we have discovered that the velocity distribution function of test particles, which are small both in mass and size compared with those of heat bath particles and are concomitantly subjected to steep potential fields, may deviate from Maxwellian significantly. We believe this is a fundamental feature of ultrafast dynamical processes where the nonlinear properties of friction are of importance [ 8,9]. Though our investigations have covered quite a broad range, including liquid and crystalline lattice heat baths and test particles with internal degrees of freedom or those subjected to strong external forces, the results are always consistent. The equilibrium velocity distribution of the test particles is always nonMaxwellian. Theoretical support for non-Maxwellian velocity distributions will be supplied in a future paper [9]. There is also accumulating experimental evidence indicating a strange behavior of hydrogen atoms in certain environments [ 10-l 61. Such experiments seem closely related to our findings. Our results have evoked much controversy, summarized by the following three questions: ( 1) Is the integration time step appropriately chosen so that it does not induce numerical errors? (2 ) Does the periodic boundary condition affect the nature of the distribution function? (3) Has equilibrium really
been established? Here we provide further evidence in support of our previous work.
studies
368
2. Finite-sized system A simple way to test the influence of periodic boundary conditions is to consider a sample with finite volume. In order to avoid possible edge effects, while still using a computational effort of realistic length, we study a two-dimensional system containing 196 host particles and 4 test particles. These particles interact with one another through LennardJones (LJ) 12-6 potentials. The LJ parameters are given in table 1. To create the necessary nonlinear dynamics, the four test particles are further subjected to an X, y corrugated potential surface having the form 4
UC=& 1 (cos&,+sinKy,), i=, where K= 56.3/u,,,, aM being the W diameter of a host particle; X, and yi designate the coordinates of the ith test particle, taking the origin to be the center of the sample box. The reduced barrier height h is ZU,/knT. In addition, a two-dimensional external potential
u,= ; A,(xy+
Y:')
t=I
Elsevier Science Publishers B.V. (North-Holland)
Volume 170,number 4
13July 1990
CHEMICALPHYSICSLETTERS
Table 1 Potential parameters
2-D c’ IIa d, IIb *’ VII d’ VIII d’ a’
%M=
0.579 0.579 0.579 0.579 0.579
0.26 0.26 0.26 0.26 0.26
0.9714 0.9714 0.9714 0.9714 0.9714
4 4 4 4 4
0.366 0.366 0.366 5.856 0.023
0.0005 0.005 0.0025 0.005 0.005
16.92 3.60 13.01 3.84 4.80
:(%+u,wM); htM=JcG.
b, N is the total number of timesteps in millions; but note that there are 8 x more data per timestep in the present computation than in ref. [6]. ‘) Present work: two-dimensional system, two-dimensional corrugated potential, d, Four fully independent runs from ref. [6]: three-dimensional system, one-dimensional corrugated potential. The distribution functions for eachof thesesystemsare depictedindividuallyin ref. [ 61. The four data sets have been averaged in this paper and plotted (circles) in fig. 2 with their standard deviations.
is applied to every particle to diminish its probability of evaporation. In eq. (2) Xi=X,/6.7o,w and Y,=y,/6.7aM represent the reduced coordinates of the ith particle. The choice of the form and magnitude of this potential is quite arbitrary. We choose A,= 16 for a test particle and AM= 1 for a host particle. Using the parameters listed in table 1, we are able to obtain a uniformly distributed system confined within a square of 13.4x 13.4&. Throughout the paper dimensions are expressed in reduced units, taking A4 and the LJ parameters u,,, and tM as measures of mass, length and energy, respectively. The equations of motion are integrated by using a modified version of the Verlet algorithm [ 171. The initial velocities of the particles, obeying a Gaussian distribution, are selected from a random number generator. In order to prevent possible numerical errors caused from an inappropriate choice of integration time step, we adopt here a computationally achievable small time step, At=O.O005r where T=u,(M/e,)“2, i.e. only one fifth to one tenth of that employed in previous work on similar systems. See ref. [ 61 and table 1.
3. Statistics A fully equilibrated system can be prepared by a preliminary integration of the equations of motion using 104At. Information about the numerical ac-
curacy of the results and the approach to equilibrium can be obtained by examining fluctuations of the total energy of the system and fluctuations in the mean square velocity of the test and host particles in the aftermath of this equilibration procedure. These fluctuations all lie within 0.05%. In order to answer the question whether the system is really in equilibrium and whether sufficient statistics have been accumulated, we perform an extremely long run. Intermediate results are extracted after every 1.2X 106 time steps starting from t = 6.12 x 106At. These results can be compared with each other to make certain that a steady-state distribution is actually reached. Altogether, 1.692~ 10’ time steps are carried out to accumulate data. This is at least three to four orders of magnitude longer than the time required for the system, initially far from equilibrium, to reach an equilibrium state. In view of the fact that four test particles are investigated and that two velocity components are taken into account, the quality of the statistics, compared with earlier results [4,6], in increased by a factor of eight. Utilization of the two-dimensional corrugated potential surface also avoids any possible superperiod problem of one-dimensional nonlinear weighted strings of the type first studied by Fermi, Pasta and Ulam [ 18,191. Our previous calculations [4,6] on three-dimensional systems with a one-dimensional corrugated potential were criticized as being “too one-dimensional” to acquire a sensibly rapid approach to equilibrium. The results to be re369
Volume 170,
number4
ported below, showing nearly identical distribution functions from the two types of systems when the parameters (m/M, aJaM, h, T) are the same, fully negates this criticism.
Graphically illustrated in fig. 1 are the distribution functions of the velocity components of the test particles (asterisks) and those host particles (circles) whose distance from one of the test particles is less than 1.5a,. We call these host particles the “cage particles”. Because of the dynamical motion, the cage includes different host particles at different instants of time. The curve symbolized by asterisks is averaged from 1.3536~ 1OSsamples, while the circles are from 6.0066~ 10’ samples. For comparison, the Maxwell velocity distribution is also plotted in fig. 1. The nearby host particles are seen to exhibit a nearly Maxwellian behavior, indicating that these particles are fairly well mixed with other host particles. In agreement with our previous work, we see from fig. 1 an unmistakable noncanonical behavior of the low-mass test particles. In fact, as emphasized in fig. 2, the distribution function is nearly identical to the ones obtained earlier for a thermodynamically sim-
Velocity
Fig. 1. Distributionsof the velocitycomponents.Solid curve: Maxwell distribution; circles: distribution function for the “cage particles”, asterisks: distribution function for the test particles. The reduced velocity here and in fig. 2 is defined as V= [ $1 (kJ/ P) P2 , where u is the velocity component for the test or host particles and p is the respective mass, m or M. The verticle axis is
scaledby percentage. 370
2.7
c2.1
0 .c ii ‘i
Gl.5 .n
4. Results
Reduced
13July 1990
CHEMICALPHYSICSLETTERS
0.9,
-
-
\.* 0.0
I
I
I
,
0.5
Reduced
I
I
I
,
1 .o
,
,
,
,
1.5
Velocity
Fig. 2. Distributions of the velocity components. Solid curve: Maxwell distribution; asterisks: distribution function for the test particles (this work); circles: average of distribution functions for systems IIa, IIb, VII and VIII of ref. [ 61 (see also table 1). The error bars are standard deviations among these four sets of data, assuming the distribution functions are independent oft,/ CM.
ilar three-dimensional system with a one-dimensional corrugated potential. The observations here, combined with those of ref. [ 61, provide evidence that (1) that distribution functions from the four ji& independent runs in ref. [ 61 average, within their standard deviations, to the same non-Maxwellian distribution function as obtained here, though the two-dimensional system needs longer time to reach equilibrium; (2) the use of periodic boundary conditions does not appreciably change the nature of the average dynamics of the test particle; (3 ) the time step employed in the previous work is not too wide; and (4) the distribution function is insensitive to the test-host interaction strength \/E,tM. An important criterion of an equilibrated system is the average kinetic energy. The present calculation gives (KM) =0.99999kBT for host particles and (Km) = 1.00027kJ for test particles, where k, indicates the Boltzmann constant, and /r,i”/~,= 0.97 14 represents the reduced temperature of the heat bath. Thus, the entire system must be in thermal equilibrium and the numerical accuracy ofthe computation is satisfactory. To ensure that the observed distribution really represents a steady-state quantity, we plot in fig. 3 the normalized probabilities, within a small range ( V, V+ A V) of velocity space, for test particle velocities as a function of the number of samples used in the statistical average. Here, V defines the reduced ve-
Volume 170,number 4
CHEMICALPHYSICSLETTERS
6z 0 .-Y 0 .-
: o_
2 0
: -6-
* *
-12
40
*
*
*
* * 8 I I I I I I I I I 8, I I 60 “if %?ples
*
* I 08 120
* 140
No.
Fig. 3. Normalized probability deviations Y( V) of test particles within a small range (I’, T/t AV) of velocity space. Circles: Y’(0); asterisks: ul( 1.19); Ac’=O.O3125.The x-axis is to be multiplied by I 06,the y-axis by 1O-‘.
locity component normalized by the one-dimensional mean square velocity Jk,Tlm of the test particles. Using PO(V) to represent the probability calculated from the Maxwell velocity distribution, we define Y( V) = 1 -P( V)/Po( V) to measure the deviations of the distribution from Maxwellian. The fluctuations of this quantity throughout the final 10.8 million time steps are seen to be small compared with their magnitudes. Moreover, as the number of samples increases, this quantity has tended to a constant that differs significantly from zero. Therefore, the distribution is stable, and the system is in a steady state. These results reply to the criticism that the system has not reached an equilibrium state or that the simulation time is not sufficiently long.
5. Conclusions The present paper is aimed at answering a number of questions raised concerning previous investigations. In order to get rid of the conventional periodic boundary conditions, we have studied here a system with finite size. The doubt about the length of the time step is dismissed by reducing At by a factor of 5-10 (see table 1). In doing so, we have obtained a distribution function nearly identical to the average (fig. 2) of that obtained from four fully independent runs for thermodynamically similar systems (same m Jh4, u,,,/uM, h, T), where the At, the dimensional-
13July 1990
ity and the boundary conditions were radically different. The results shown in fig. 2 thus indicate that the velocity distribution functions, though decidedly non-Maxwellian, are independent of dimensionality. The results also suggest that they are insensitive to the test-host interaction energy CM. To obtain better statistics, we have performed a very long molecular dynamics simulation. In addition, four test particles and two velocity components were considered, further improving the statistics by a factor of eight. The smoother distribution function obtained here is undoubtably a result of these improvements, but may also be partly a result of using the shorter timestep. Fluctuations of the probabilities that test particles remain within a small range of velocity space during the simulation time have also been investigated. The smallness of these fluctuations and the approach of the probabilities to fixed limiting values provide evidence that the system has reached its equilibrium state on a time scale much shorter than that employed to accumulate data. Additional evidence for the system being in equilibrium is given by the fact that the average kinetic energies of the test and host particles are equal and the total energy of the system is conserved to within 0.05% throughout the run. In view of the above results and discussion, we conclude that the observed noncanonical behavior is a real characteristic of equilibrated systems in which ultrafast dynamical processes are taking place. It is not a consequence of computationally introduced errors or poor statistics [20]. The effect of these distorted velocity distributions on velocity-dependent cross sections in chemistry and physics will radically change current views about large classes of thermally activated condensed phase reaction processes.
Acknowledgement Financial support at the SPQR Laboratory been shared by the Robert A. Welch Foundation 0005, 30% and D- 1094, 5%)) the State of Texas vanced Research Program (1306, 47%) and Pittsburgh Supercomputing Center ( 18%).
has (DAdthe
371
Volume 170, number 4
CHEMICAL PHYSICS LETTERS
References [ 11 S.-B. Zhu, J. Lee and G.W. Robinson, I. Chem. Phys. 88 (1988) 7088. [2] S.-B. Zhu and G.W. Robinson, Chem. Phys. Letters 153 (1988) 539. [ 31 S.-B. Zhu and G.W. Robinson, J. Phys. Chem. 93 (1989)
164. [4] S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 161 (1989) 249. [ 51 S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 163 (1989) 328. [ 61 S.-B. Zhu, J. Lee and G.W. Robinson, Chem. Phys. Letters 169 (1990) 355. [7] S.-B. Zhu, J. Lee and G.W. Robinson, I. Phys. Chem., submitted for publication. [ 81 S.-B. Zhu, J. Lee and G.W. Robinson, Velocity Dependence of Friction, Chem. Phys., submitted for publication. [ 91 S.-B. Zhu, Velocity Distributions in Nonlinear Systems, Phys. Rev. A, submitted for publication.
372
13 July 1990
[ IO] G. Comsa, R. David and B.-J. Schumacher, Surface Sci. 95 (1980) L210. [I l] G. Comsa and R. David, Surface Sci. 117 ( 1982) 77. [ 121 S. Estreicher, Phys. Rev. B 36 ( 1987) 9 122. [ 131 S.J. Pearton, J.W. Corbett and T.S. Shi, Appl. Phys. A 43 (1987) 153. [ 141 M. Fleischmann, S. Pons and M. Hawkins, J. Electroanal. Chem. 261 (1989) 301. [ 151 SE. Jones, E.P. Palmer, J.B. Czirr, D.L. Decker, G.L. Jensen, J.M. Thome, S.F. Taylor and J. Rafelski, Nature 338 (1989) 737. [ 161 R.J. Beuhler, G. Friedlander and L. Friedman, Phys. Rev. Letters 63 (1989) 1292. [ 171 W.C. Swope, H.C. Anderson, P.H. Berens and K.R. Wilson, J. Chem. Phys. 76 (1989) 637. [ 181 E. Fermi, J.R. Pasta and SM. Ulam, in: Collected works of Enrico Fermi, Vol. 2 (University of Chicago Press, Chicago, 1965). [ 191 J.I. Tuck and M.T. Menzel, Advan. Math. 9 ( 1972) 399. [ 201 W.P. Keirsteadand K.R. Wilson, J. Phys. Chem. 94 ( 1990) 918.