918
J . Phys. Chem. 1990, 94, 918-923
A Breakdown of Equilibrium Statistical Mechanics? W. P. Keirsteadt and Kent R . Wilson*,* Department of Chemistry, B-039, and Institute f o r Nonlinear Science, University of California, San Diego, La Jolla. California 92093 (Received: March 31, 1989; In Final Form: July 27, 1989)
We consider the recent paper “Ultrafast Dynamics of a Quasi-Dissociative Diatomic Molecule in Solution”by Zhu and Robinson, in which the results of molecular dynamics simulations and their interpretation by the authors seem to contradict basic principles of classical equilibrium statistical mechanics. We discuss both the results expected from equilibrium statistical mechanics and the results of tests by computer simulation. Two types of simulations are presented: (i) a “canonical”case with frequent connections to a heat bath, and (ii) a “microcanonical”case which repeats a critical portion of Zhu and Robinson’s molecular dynamics. We pay particular attention to the question of whether the numerical data can justify the rejection of equilibrium statistical mechanical expectations.
Introduction In a recent paper’ in this Journal, Zhu and Robinson present results of molecular dynamics (MD) trajectory calculations for a model of a chemical reaction in solution, a quasi-dissociative diatomic molecule in a rare gas solvent. They conclude that for appropriately chosen potential parameters neither the velocity along the line of centers of the two diatoms nor the velocities of the nearby solvent atoms are Maxwell-Boltzmann (MB) distributed. Furthermore, they find that this noncanonical distribution is a function of the distance between the two diatomic atoms (the reaction coordinate). From their findings of local nonequilibrium and noncanonical behavior, they draw a number of conclusions as to the microscopic nature of chemical reactions in solution. Zhu and Robinson have reached related conclusions in studies of cis-trans i s ~ m e r i z a t i o n , ~but - ~ the quasi-dissociative system which they have recently presented is a particularly clear and simple case which avoids such complexities as angle-dependent internal coordinates and the care which must be taken to use properly canonically related coordinate and momentum representations. Their results are ones which we find surprising, as they seem to contradict the predictions of classical equilibrium statistical mechanics. Consider a reacting system in thermal equilibrium with a solvent where reactions are caused by equilibrium fluctuations. We expect that, whenever the Hamiltonian is of the form p2/2m plus terms not depending on p , where m is independent of coordinate, the distribution function for p will be proportional to the Boltzmann factor exp(-@p2/2m), regardless of any restriction on the sampling range involving the other coordinates or momenta of the system. In the Appendix, we present this argument in greater detail for the Zhu-Robinson quasi-diatomic system and show that we expect from equilibrium statistical mechanics that both the velocities of the nearby solvent atoms and the relative velocity of the diatomic along the reaction coordinate will be indeed MB distributed, independent of the distance between the diatoms. Although we cannot believe for the above reason that Zhu and Robinson’s results reflect a basic property of chemical reactions in nature, it is of interest to investigate their importance for MD simulations, which are only approximations to nature. Thus, we present both what we consider a better simulation approach (in the sense that it is closer to experimental reality) and a repetition of a critical portion of Zhu and Robinson’s numerical work. We analyze the resulting trajectories with respect to the questions raised above, paying particular attention to the statistical uncertainties in the data as well as to the effects of initial state preparation and simulation technique. The outline of this paper is as follows. We first describe the quasi-diatomic system and the MD methods used in our simulations. We discuss the sampling and statistical techniques used
’
Department of Chemistry and Institute for Nonlinear Science. *Department of Chemistry.
0022-3654/90/2094-0918%02.50/0
TABLE I: Potential Parameters for the Simulations 1490
1.293
4.117
4
TABLE II: Run Specific Parameters run
type‘
1-20 21
c
m
elapsed time, ns
time step, fs
density, g/cm3
T,b K
1.5 20.14
5 5
1.3447 1.3441
296.587 295.977
a c = “canonical” run, m = “microcanonical”run. bThe temperature is determined from the kinetic energy averaged over the entire simulation.
to generate and analyze the trajectory data. We then consider the results from the viewpoint of whether they are consistent with classical equilibrium statistical mechanics and in particular whether the velocity distributions of the diatomic along the reaction coordinate and of the nearby solvent atoms are MB distributed. Finally, we try to express clearly some basic distinctions among the concepts of equilibrium/nonequilibrium and adiabatic/ nonadiabatic as applied to reactive trajectories in the hope of avoiding future confusion.
Simulation Details We study the molecular dynamics for Zhu and Robinson’s system I1 for which they find the greatest deviation from MB expectations. This system consists of 106 argon-like solvent atoms and a diatomic molecule which has nitrogen atomic masses, but with a special quasi-dissociative potential. The argon-argon and argon-nitrogen interactions are Lennard-Jones potentials of the form =
(
:)I2
-
($1
where rij is the distance between atoms i and j , while the diatomic potential is given by a symmetric double well where the first minimum corresponds to a bound state and the second to a quasi-dissociativestate. The precise form of the diatomic potential is U ( r ) = Uo[l - A ( r - ro)2+ B(r - ro)2n] (2) where A = n / [ ( n- l ) b 2 ] , B = A / [n b2n-2 1, 6 = Ir2 - rol, and ro = ( r l + r 2 ) / 2 . Here, rl and r2 > rl are the locations of the two potential well minima. The specific parameters e and D used for ( 1 ) Zhu, (2) Zhu, (3) Zhu, (4) Zhu, (5) Zhu,
S.-B.; Robinson, G. W. J . Phys. Chem. 1989, 93, 164. S.-B.;Lee, J.; Robinson, G. W. J . Phys. Chem. 1988,92,2401. S.-B.;Lee, J.; Robinson, G. W. J . Chem. Phys. 1988,88, 7088. S.-B.; Robinson, G.W. Chem. Phys. Lett. 1988, 153, 539. S.-B.;Robinson, G . W. Proc. 3rd Inr. Supercomput. Insr. 1988,
I, 300.
0 1990 American Chemical Society
A Breakdown of Equilibrium Statistical Mechanics? r 5.0 4.0
0.0 1 ,
1.0 I
.
I
(4 3.0
2.0
The Journal of Physical Chemistry, Vol. 94, No. 2, 1990 919
1
t iR
4.0 5.0 . , . , . , . I 0.6
( a ) i,
1
0.5
0.4
0.3 D
0.2 0.1
0.0
1-1 $1.00
0'95 0.90
u
1.05 Y
1.00
0.45
0.40
1 '
0.0
1.0
2.0
r
3.0
4.0
5.0
(4
Figure 1. Output from the "canonical" simulation, plotted as a function of the reaction coordinate r. (a) The dotted line is the gas-phase diatomic potential energy U ( r ) and the circles are numerically calculated probabilities p ( r ) to be in bins along the reaction coordinate r. The asymmetry between the populations of the wells is a result of solvent forces. (b) The , by average kinetic energy in the reaction coordinatep , 2 / 2 ~normalized k B T / 2 . (c) The average kinetic energy of an inner core solvent atom, normalized by 3 k B T / 2 . (d) The fraction of configurations with negative reaction coordinate velocity. The error bars are h one standard deviation, as discussed in the text.
the Lennard-Jones potentials, and rl, rz, n, and U, for the diatomic potential, are summarized in Table I, and the diatomic potential is shown as the dashed curve in Figure la. The run specific parameters are given in Table 11. For our computer simulations, we use periodic truncated octahedral boundary conditions and the leapfrog Verlet integration alg~rithm.~.'A time step of 5 fs is used, which gives excellent energy conservation even for our very long runs. We handle the potential cutoff problem by feathering the potentials smoothly (6) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (7) Swope, W. C.; Anderson, H. C.; Berens, P. H.; Wilson, K. R. J . Chem. Phys. 1982, 76, 637.
to zero at the boundary edges,'*8 thus avoiding force discontinuities. Our trajectories are generated from 20 seed states prepared in the following way. First, the atoms are randomly placed in a simulation box of an appropriate size to give the desired overall density, the only constraint being that no atom be too close to any neighboring atom. This initial configuration is then minimized in potential energy by using a conjugate-gradient m e t h ~ d . ~ Finally, each system is equilibrated to a temperature of 298 K by choosing an initial random set of velocities corresponding to a MB distribution and then integrating constant temperature dynamics'O for 500 fs, rerandomizing the velocities according to a MB distribution and running for another 500 fs, and so on. This process of randomizing/integrating is repeated 20 times, for a total of 10 ps of equilibration for each seed. We analyze our results in a manner similar to Zhu and Robinson, that is, by studying the behavior of relevant quantities as a function of reaction coordinate r, the distance between the two diatoms. The reaction coordinate is binned in 11 boxes of equal length over the range 0.7635-4.6465 A (in order that there be bins centered at the two minima and at the maximum of the diatomic potential), and distributions and averages in each bin are then calculated from those configurations where the reaction coordinate is in the bin range. We collect and analyze the distribution of the momentum along the reaction coordinate, and the kinetic energies of the solvent atoms, both inner core (those with distances from the diatomic less than the average solventsolvent separation 4.1 17 A) and outer core (those solvent atoms not in the inner core). We compute as a function of reaction coordinate (i) the average kinetic energy corresponding to the momentum along the reaction coordinate, (ii) the inner core average kinetic energies, and (iii) the fraction of configurations f-where the velocity along the reaction coordinate, v, is negative (the molecule is contracting). An important addition in our work, as compared to the study of Zhu and Robinson, is a consideration of the statistical significance of the numerical results. The kinetic energy and velocity sign averages presented in the figures are taken over all the time steps accumulated in relevant runs. As pointed out to us by Ciccotti and Hynes,]' because the MD data are computed at short time intervals compared to the decay of the relevant velocity autocorrelations, there are strong serial correlations in the data which make the effective independent sample size much smaller than the full data set. To correct for these correlations in our measurements, we employ the method outlined in Allen and Tildesley6 to compute the standard errors in the mean. The data for any given bin is split into m blocks, the size b of which is large compared to some correlation time for the system (which we take to be about 150 fs, corresponding approximately to the l / e time for decay of the solvent velocity autocorrelation function). The relevant averages (A),, are calculated within each block, and then the standard deviation " ( ( A ) b )in these block averages is found. When b is large compared to the correlation time, 02(( A ) b ) l / b , and hence scaling the block standard error by b'lz should give a quantity approximately independent of b. The standard , ) the overall measured average (A)m,,,where m error u ( ( A ) ~ " ,in 1 , is then u ( ( A ) b ) / m ' / z This . quantity should be fairly independent of b and m if there are a sufficient number of counts to analyze. In our analysis, we typically choose m = 25, although we have also examined results with a greater number of blocks. The calculated errors for the bins with the fewest number of counts can vary by as much as 13% as a function of m, although the bins with large numbers of counts will have errors smaller than 1%. (Since the bins in regions of high potential energy, such as near the barrier maximum, have small populations, their standard errors
-
-
(8) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J . Chem. Phys. 1987, 86, 1356. (9) Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling, W. T. Numerical Recipes in C; Cambridge University Press: Cambridge, 1988; p 317. (10) Nosi, S . Mol. Phys. 1984, 52, 255. ( 1 1) Ciccotti, G.; Hynes, .I. T. University of Colorado, Boulder, CO, private communication.
920
Keirstead and Wilson
The Journal of Physical Chemistry, Vol. 94, No. 2, 1990
in the mean will be considerably larger than for bins near the potential wells.)
Numerical Results: Canonical Case We have chosen two different methods for carrying out the MD simulations. For the work described in this section, we choose 20 independently prepared initial seed configurations, as described in the previous section. From each equilibrated seed, we choose an initial set of velocities randomly selected from a MB velocity distribution at temperature 298 K and run constant energy dynamics for 1.5 ns, using a time step of 5 fs and collecting data every 10 fs (for a total elapsed time of 30 ns). Every 100 ps the velocities are rerandomized according to the MB distribution. A bin average is then taken over all configurations in all runs where the reaction coordinate lies within the given bin. Because we choose the initial velocity distribution canonically, and because we periodically rerandomize the velocities according to the canonical MB distribution to approximately take into account the heat bath of the rest of the solution, we refer to this case as the “canonical” one. We find that barrier crossings occur on average about every 20 ps and that the time for any given transition is around 20 fs. In Figure la, we present the distribution function p ( r ) for the position r of the diatomic along the reaction coordinate superimposed on the diatomic potential U(r)as a function of r . In Figure 1 b, we show the normalized average kinetic energy in the reaction coordinate, (p:/2p),/(kBT/2) (in the notation of the Appendix), versus r. In Figure IC, we show the normalized average inner core atomic kinetic energy versus r , normalized by a factor 3kBT/2. In Figure Id, we present the velocity sign fractionf versus r . The solid lines in each of these figures represent the expected canonical result based on the theoretical calculations shown in the Appendix, and the error bars (A one standard deviation) are obtained as described in the previous section. As can be seen in Figure I , our results are consistent with MB expectations within the computed uncertainties. The regions with the poorest agreement with MB are seen to be precisely those visited least often. If one ran a longer simulation and connected with the heat bath adequately, the results should be in even better agreement with MB. We chose not to try to improve our statistics further because of the large computer time requirements. Comparing our results with those of Zhu and Robinson, we find our results are considerably closer to the MB ones expected for a canonical ensemble. In their Figure 10, they find bins in which the value of the normalized kinetic energy in the reaction coordinate differs from unity by as much as +22%. Our observed distribution as shown in Figure l b is considerably flatter and deviates from unity by no more than +7%. In addition, the basic shape of our distribution does not agree with theirs, evidence that the observed fluctuations are random in nature and do not reflect any underlying physical property of the system. I n their Figure 2, they find that the normalized inner core solvent kinetic energy is about 6% below the canonical value near the first potential minimum, and monotonically increases to about 20% above the canonical value near the second potential minimum. In Figure I C , we find no such monotonic trend and observe deviations of at most 2%. Furthermore, our outer core solvent results are nearly identical with the inner core ones. I n our Figure Id, we see that f-deviates significantly from one-half only in the central three bins near the barrier maximum, where our statistics are poorest, and the magnitude of the deviation is at most 1 o/r low. This result is in dramatic contrast to Zhu and Robinson’s observations for the velocity sign ratio S, which is simplyf-/(I -f-) (their Figure 9). They find that S deviates from unit) by as much as 60%. Furthermore, they observe a particular shape to the distribution in which S follows the negative sign of the barrier force. We observe no such trend. The central question for this simulation is whether or not the M B distribution for the reaction coordinate and nearby solvent atom velocities can be rejected based on the numerical evidence found. Zhu and Robinson claim that the shape of the velocity distribution is not Gaussian and present evidence in terms of a
f
0.00 -0.50
t c
-1 00 - 1 50
- 2 00
, 01
1
1
I
l
l
!
,
I
-
I
5 1 0 2030 5 0 7080 9 0 9 5
99
9 9 9 9 9 99
Percent
Normal probability plots for the reaction coordinate velocity for the “canonical”simulation, which for a Maxwell-Boltzmann distribution would give a straight line. The solid line is data from the bin centered at the first potential minimum ( r = 1.293 A), the dashed line from the bin centered at the maximum ( r = 2.705 A), and the dotted line from the bin centered at the second minimum ( r = 4.1 17 A). In each case we plot 1000 points chosen at equally spaced intervals from the full data set for the bin. Figure 2.
computed ratio of moments ( V * ) ~ / ( U ~ ) We . suspect that this is not a particularly robust statistical test in the presence of inevitable uncertainties at the tails of the velocity distribution due to low sampling probability. Instead we present a more direct test of the agreement with the theoretical distribution function by producing a normal probability plot.12 In Figure 2, we show such a plot for the data from three of the bins-the two bins centered at the potential minima, and the one bin centered at the potential maximum. The expected MB distribution would give a straight line. We see that our data points lie very close to a straight line, except near the tails. The deviations at large absolute values of velocity are not unexpected since it is difficult to get good sampling for these unlikely large velocity fluctuations in a reasonable length simulation run.
Numerical Results: Microcanonical Case We have argued thus far that neither experimental systems under typical experimental situations nor suitably designed MD simulations-such as our “canonical case” above-would give velocity distributions contradicting with statistical significance the MB result. The question remains, however, as to whether or not other forms of MD simulation would give the same answer. A particularly simple MD technique is to perform a single long run using a constant energy integration technique. In this method, one is sampling a constant energy surface in phase space and is hence directly simulating a microcanonical ensemble (we refer to this case as “microcanonical”). Although the results of any simulation should be independent of which ensemble one uses in the thermodynamic limit of large number of particles, perhaps \,,henone does microcanonical simulations using systems with a limited number of atoms there are significant deviations from the expected results for the canonical ensemble. Because the total system energy must be constant in the microcanonical simulation, any fluctuation in energy in one part of the system must be compensated for in the other parts. In the Zhu-Robinson system, for example, when the diatomic molecule is at the barrier top, the extra 3 kcal/mol of potential energy required for this spatial configuration will take away an equal amount from the total of the potential and kinetic energies in the other degrees of freedom. In our simulations, the total available energy K + V - Vminis around 160 kcal/mol, and hence such a move by the diatomic can take a nonnegligible fraction of the total available energy away from the other degrees of freedom in the system. It is then conceivable that the solvent or diatomic velocity ( 1 2) Rice, J. A. Mathematical Statistics and Data Analysis; Wadsworth & Brooks/Cole Advanced Books & Software: Pacific Grove, CA, 1988; p 292.
The Journal of Physical Chemistry, Vol. 94, No. 2, 1990 921
A Breakdown of Equilibrium Statistical Mechanics? 1.0
0.0
r (A, 2.0 3.0
4.0
5.0
0.6
4.00
I
3.00
1
1
,
I
/
/
/
,
,
I
/
/
/
I
I
'1
-r = 1.293 A - - r = 2.705 A
0.5
4.0
0.4
3.0
0.3 2.0
14
0.2 - 1 .oo
1 .o
0.1 -2.00
0.0
0.0
-3.00
cr"' I .
~
.01
,1
1
5 1 0 2030 50 7080 9 0 9 5
99
99.999.99
Percent
0.901
'
"
"
"
Figure 4. Normal probability plots for the reaction coordinate velocity for the "microcanonical" simulation, which for a Maxwell-Boltzmann distribution would give a straight line. The solid line is data from the bin centered at the first potential minimum ( r = 1.293 A), the dashed line from the bin centered at the maximum (r = 2.705 A), and the dotted line from the bin centered at the second minimum (r = 4.1 17 A). In each case we plot 1OOO points chosen at equally spaced intervals from the full data set for the bin.
' " 3
1.05
0.55 I c
.
energy dynamics for 20.14 ns, using a 5-fs time step and collecting data every 10 fs. In Figures 3 and 4, we show the microcanonical versions of the same quantities shown in Figures 1 and 2. The data are generated and analyzed in the same manner as in the canonical case. Examination of the plots in Figures 3 and 4 reveals that the magnitude of the deviations from the MB velocity distributions are comparable to those found in the canonical simulation. In Figure 3c, however, one can clearly see that the average inner core solvent kinetic energy varies systematically with reaction coordinate, in precisely the way one would expect: when the diatomic is at the potential maximum, the inner core solvent energy is reduced; when the diatomic is at either of the potential minima, the solvent energy is increased. The outer core results are nearly identical. We may try to quantitatively account for this systematic deviation using the following simple equilibrium argument.13 When the diatomic is at reaction coordinate r, it takes energy -6E 61/(r) out of the solvent, in which 6U(r) is the excess potential energy removed. The corresponding change in the temperature (and hence total kinetic energy) of the solvent should be given by 6T = 6E/C,, where C, is the heat capacity at constant volume. We calculate C, from our microcanonical trajectory data using the following relationship6J4
-
0.50
0.45 0.40
0.0
1.0
2.0
3.0
(4
4.0
5.0
Figure 3. Output from the "microcanonical" simulation, plotted as a function of the reaction coordinate r. (a) The dotted line is the gas-phase diatomic potential energy U(r)and the circles are numerically calculated probabilities p ( r ) to be in bins along the reaction coordinate r. The asymmetry between the populations of the wells is a result of solvent forces. (b) The average kinetic energy in the reaction coordinate p > / 2 p , normalized by k 6 T / 2 . (c) The average kinetic energy of an inner core solvent atom, normalized by 3k6T/2. The dashed curve is the theoretical calculation discussed in the text. (d) The fraction of configurations with negative reaction coordinate velocity. The error bars are i one standard deviation, as discussed in the text. distributions could be measurably different when the diatomic is at the barrier top than when it is at the barrier minimum, provided statistical noise does not wash out the effect. In the "canonical" simulation we describe above, we overcome the limitation of a small number of degrees of freedom by periodically introducing a heat bath which effectively couples the system to an infinite energy reservoir. To investigate the consequences of constant energy dynamics for a system with a small number of particles, we carry out the following simulation. We take one initial seed (prepared as above) and generate its initial velocity distribution choosing atom velocities randomly from the MB distribution. We then integrate constant
The change in temperature 6T may be related to the average kinetic energy of an inner core solvent atom, and hence we may calculate as a function of r the deviation from unity of the normalized value of this quantity. In Figure 3c, the dashed line represents the results of this calculation. The fit is quite close and leads one to suspect that this equilibrium picture is the correct interpretation of the observed systematic deviations. In summary, then, we find that even a constant energy (microcanonical) simulation gives results in good agreement both with physical expectations and with computer simulations which attempt to introduce the heat bath more directly. The systematic deviations from MB as a function of reaction coordinate in the average solvent kinetic energy for the microcanonical case can be accounted for by the above equilibrium argument, and the remaining statistical errors are rather small. Certainly, we find no valid reason to reject the results expected from equilibrium (13) Wheeler, J. University of California, San Diego, CA, private communication. (14) Lebowitz, J . L.; Percus, J. K.;Verlet, L. Phys. Reu. 1967, 153, 250.
922
The Journal of Physical Chemistry, Vol. 94, No. 2, 1990
statistical mechanics based upon our numerical simulations.
Conclusion Zhu and Robinson have presented molecular dynamics (MD) calculations for a simple model of a chemical reaction in solution, consisting of a quasi-dissociative diatomic molecule in a raregas-like solvent. They interpret their results as showing that both the momentum along the reaction coordinate and the momenta of the nearby solvent atoms are non-Maxwell-Boltzmann (MB) distributed and that this noncanonical distribution is a function of reaction coordinate. We have repeated their MD calculations for the particular case where they note the greatest deviations from canonical behavior, and have supplemented them with MD calculations which we believe better reflect experimental reality. We conclude, in contrast to their paper, that the sampled velocities both along the reaction coordinate and in the nearby solvent atoms are consistent with the expectations of classical equilibrium statistical mechanical analysis. For our "canonical" MD simulations, the results agree with MB distributions within statistical uncertainties. For our "microcanonical" calculations, the observed systematic deviations from the MB distribution can be accounted for by equilibrium analysis and are in any case rather small. Thus we find no substantiation for the view that nonequilibrium effects on velocity distributions (in particular effects dependent on reaction coordinate position) are present in this model for a chemical reaction in solution. With respect to other points raised in the paper by Zhu and Robinson, which we do not wish to discuss at length here, we suspect it may be useful in trying to think clearly about reactions in solution to try to carefully distinguish among the concepts of equilibrium/nonequilibrium versus adiabatic/nonadiabatic. Lack of clarity with respect to such concepts may lead to unnecessary confusion and miscommunication. Our view is if one takes a reacting system in a heat bath which has come to equilibrium, for which (i) the Hamiltonian can be written in terms of a proper conjugate set of momenta and coordinates (one of which being the reaction coordinate), (ii) the Hamiltonian separates into a kinetic energy term which is a function only of the square of a particular momentum and other terms which are not a function of that momentum, and (iii) the mass associated with this momentum is constant and independent of coordinate, then one will find that the particular momentum is distributed in a normal equilibrium canonical MB distribution, independent of reaction coordinate position. This does not imply that these distributions will remain the same as equilibrium ones if a subset selected to be within a certain reaction coordinate range is later allowed to evolve in time. Then the specialness or nonequilibrium aspect of the reaction coordinate (which was selected for that sample to be in a particular interval and which was thus not distributed in an equilibrium manner) can spread into other coordinates and momenta and they in time can also become distributed in a nonequilibrium manner. Whether the other degrees of freedom do or do not become distributed in a nonequilibrium manner as the selected sample evolves and the reaction coordinate distribution spreads out from the originally selected (r, r dr) interval is a question as to whether the other degrees of freedom do or do not adiabatically follow the change in the reaction coordinate r. Such adiabaticity or nonadiabaticity of the other coordinates of the reacting system with respect to motion along the reaction coordinate involves the detailed coupling among the degrees of freedom and thus depends on the particular Hamiltonian and conditions such as the temperature and pressure of the reacting system. We should therefore indeed be surprised if we find that the other degrees of freedom of such an otherwise equilibrated reacting system were immediately distributed in a noncanonical nonequilibrium manner just because we selected a sample only within a certain range of the reaction coordinate, but we should not be at all surprised if a decent interval after selection the momenta in the sample have evolved away from canonical equilibrium distributions.
Keirstead and Wilson discussions, and the National Science Foundation for supporting this research. Appendix For completeness, we include a detailed equilibrium statistical mechanics calculation of the quantities discussed in the text. Consider the diatomic molecule and solvent atoms as being in a closed container in thermal contact with an external heat bath which maintains the system at temperature T. Denote by rl and r2 the coordinates of the two diatomic atoms, and by xi the coordinates of the solvent atoms. Let m be the mass of one of the diatomic atoms and m, be the mass of one of the solvent atoms. The Lagrangian for the system is L = Y2m(il2+ i,2)+ C'/2m,x: I
Here U is the internal diatomic potential energy and V,,, is the internal solvent potential energy plus the interaction energy with the diatomic. If we replace the diatomic Cartesian coordinates by spherical coordinates R = (r, + rz)/2 and r = ri - r2, and convert to the Hamiltonian formalism, we find that P? Psz H=-+-+ 2~ 2pr2
PQZ 2prZ sinZ8
E"' i
2m,
PZ + -2M +
+ U ( r ) + Vnb(r,8,$,R,x, ,...,xN) (5)
where p = m / 2 and M = 2m. The conjugate momenta in the above equation are related to the velocities in the Lagrangian formalism by p, = p i , ps = pr*%,p s = p r2 sinZ84, P = M R , xi = m,ki. The probability that the system is in any given microstate is proportional to the Boltzmann factor e-flH.The probability density that the system has a particular value of pr and r, regardless of the values of the other coordinates and momenta, is proportional to the integral of the Boltzmann factor over the unconstrained degrees of freedom. Applying this result to the above Hamiltonian gives prob @,,r) = C(r)e-flP:/Zpe-flu(r)
(6)
where C ( r ) is a function solely of the reaction coordinate r a n d is related to the potential of mean force produced by the average influence of the solvent. Let A@,) be some function of p, which is independent of the other coordinates and momenta. Then the average of A@,) over an ensemble restricted to those microstates having reaction coordinate r is *jl:Abr)
rob @r*r) dpr
(A@,)),= Jm
Prob @ J )dPr
+
Acknowledgment. We thank J. T. (Casey) Hynes, Raphy Levine, John Rice, Wilse Robinson, and John Wheeler for helpful
So, if A@,) is independent of r, then the ensemble average of A@,) restricted to those microstates with reaction coordinate r is independent of r. Applying this result to two of the quantities of interest in this paper, we find
f- =
no. of configurations with p r < 0 = 1/2 total no. of configurations
(9)
We may also consider the averages of the kinetic energies of atoms in the inner core solvent (the outer core solvent is handled
J . Phys. Chem. 1990, 94, 923-929 similarly, and the results are the same). Consider a function B(lri) which is independent of the other coordinates and momenta. The average of B ( r i ) over an ensemble restricted to those microstates where the coordinate xi satisfies min (Ixi - rll, Ixi - r21) < d is given by
S_: B(ri) ~ r o b l c
( T i ) d3ri
(B(ri))lc =
where problc (ni) is the probability density that the solvent atom has momentum lri and is in the inner core. Because the ri term in the Hamiltonian separates out, problc (ni)is proportional to exp(-Plr,2/2ms), and the average of B ( r i ) has no position dependence, either on the reaction coordinate r or on other solvent atom coordinates xi. For the particular case where B is the kinetic energy, this gives
(10)
~ r o b l c(Ti)d 3 r i
923
( ~ ; / 2 m , ) ~=~ 3kBT/2
(1 1)
The same result holds for the outer core solvent.
Accurate Thermodynamic Properties of the Six Isotopomers of Diatomic Hydrogen Robert J. Le Roy,* Steven G. Chapman, and Frederick R. W. McCourt Guelph- Waterloo Centre for Graduate Work in Chemistry, University of Waterloo, Waterloo, Ontario N 2 L 3G1, Canada (Received: April 20, 1989; In Final Form: August 1 1989) ~
Tabulations of the thermodynamic properties of the six isotopomers of diatomic hydrogen have been generated from a knowledge of the ab initio nonadiabatic eigenvalues for all bound and quasibound levels of the ground electronic state. The new results for H,, HD, and D2 should be. more accurate than the existing JANAF values generated from an eigenvalue spectrum obtained empirically by interpolating over and extrapolating beyond the available spectroscopic data. For the tritium isotopomers, the present results comprise the first comprehensivetabulation of the thermodynamic properties. Comparisons between ortho, para, equilibrium, and "frozen" (or "normal") results for the homonuclear systems are also presented.
I. Introduction Although diatomic hydrogen is already one of the most thoroughly characterized molecular species, it remains topical because of its importance in many physical and chemical processes and because of its central role as a meeting ground for theory and experiment. It was one of the first systems for which detailed statistical mechanical calculations were employed to generate comprehensive tables of thermodynamic properties.' The landmark tabulations of Wooley et al.' were based on level energies obtained empirically by interpolating over and extrapolating beyond the available spectroscopic data, and many subsequent tabulations were generated in that same manner.2 For diatomic hydrogen, however, this approach was superseded in 1973 when Kosloff et aL3 reported more accurate thermodynamic property calculations based on a complete set of purely theoretical level energies. This was the first example of an accurate determination of macroscopic properties directly from first principle^.^ Moreover, it appears to have been the first such calculation to take proper account of the metastable nature of quasibound levels by integrating over the continua associated with finite line widths. However, as in the calculations of Wooley et al.,' Kosloff et aL3 considered only the three common isotopomers, HZ,HD, and D2. In recent years, a b initio calculations of the properties of molecular hydrogen have continued to improve in both accuracy and range. In particular, there have been significant improvements in both the clamped-nuclei (simple Born-Oppenheimer) potential5 and its adiabatic corrections.6 Effective adiabatic eigenvalues obtained on combining these functions with older estimates of the relativistic and radiative correction^'^ have been corrected by using (1) Wooley, H. W.; Scott, R. B.; Brickwedde, F. G. J . Res. Nail. Bur. Stand. 1948, 41, 379.
(2) Chase, M. W., Jr.; Curnutt, J. L.; Downey, J. R., Jr.; McDonald, R. A.; Syverud, A . N.; Valenzuela, E. A. J . Phys. Chem. Ref. Data 1982, 1 1 , 695. (3) Kosloff, R.; Levine, R. D.; Bernstein, R. B. Mol. Phys. 1974, 27, 981. Erratum: Ibid. 1976. 31. 323. (4) Unfortunately,'the detailed results of Kosloff et al.'seem to have been overlooked in the 1982 revisions of the JANAF thermochemical tables2 (5) Kobs, W.; Szalewicz, K.; Monkhorst, H. J . Chem. Phys. 1986, 84, 3278. ( 6 ) Wolniewicz, L. J . Chem. Phys. 1983, 78, 6173. (7) Kolos, W.; Wolniewicz, L. J . Chem. Phys. 1964, 41, 3663. (8) Bishop, D. M.; Cheung, L. M. J . Chem. Phys. 1978,69, 1881.
0022-3654/90/2094-0923$02.50/0
a generalized representation of the nonadiabatic level shifts, yielding accurate nonadiabatic energies for all bound and quasibound vibration-rotation levels of the six isotopomers of ground-state hydrogen.'OJ1 Since certain of these energies differ significantly from those used previously, thermodynamic properties calculated from them might be expected to be improved significantly over those reported earlier. These considerations, taken together with the fact that there existed no comprehensive tabulations of thermodynamic properties for the tritium isotopomers, motivated our work. The present paper reports new calculations of the thermodynamic properties of all six isotopomers of the ground electronic state of diatomic hydrogen. In the following, section I1 describes the input data and outlines how our calculations have been performed. The results are then presented in section 111, together with an examination of the differences between the present results and those of ref 3, and a discussion of questions such as the effect on the results of the improvements in the level energies and the effect of the finite width of the quasibound levels. 11. Methodology
A. The Input Energies. The level energies used in the present calculations are the nonadiabatic eigenvalues reported by Schwartz and Le R O Y . ' ~ JThey ~ were obtained by solving the radial Schrodinger equation for the best existing potential energy curves and adding nonadiabatic level shifts to the resulting eigenvalues. The latter corrections were based on directly calculated nonadiabatic level shifts for H2, HD, and D26 but were generalized both to allow extrapolation to the tritium isotopomers and to take centrifugal distortion effects into account." As a result, the level energies calculated for HT, DT, and T2 are believed to have the same accuracy as those for H,, HD, and D,, and unlike the input energies used by Kosloff et al.,3 the accuracy of these level energies does not deteriorate with increasing rotational energy. The maximum errors in the calculated level energies due to the computational procedures used to generate them are believed to be less than 0.02 cm-I, and comparisons with the extensive (but (9) Garcia, J. D. Phys. Rev. 1966, 147, 66. (10) Schwartz, C.; Le Roy, R. J. J . Mol. Spectrosc. 1987, 121, 420. (1 1) Le Roy, R. J.; Schwartz, C. University of Waterloo Chemical Physics Research Report CP-301R, 1987.
0 1990 American Chemical Society
