2967
J . Phys. Chem. 1991, 95, 2967-2970 "antisymmetric" N - C 4 stretch that is strongly weighted toward the N-C single-bond motion. The (inactive) amide I band is then primarily C 4 stretch with a smaller amount of C-N component with an in-phase relative displacement. Electronic excitation is expected to select the antisymmetric form in the sense that any change in the C - 0 bond length will be opposite to the change in the C-N bond length. Hydrogen bonding or the effect of a high dielectric constant is expected to decrease the C 4 force constant and raise the C-N force constant. This will increase the coupling and change the resulting eigenvectors in a space spanned by the C=O and C-N stretching motions. The result will be a new distribution of intensity which will not necessarily include both modes of motion. The magnitude of this effect is very difficult to estimate although it should be amenable to theoretical treatments that include the effects of solvation.'*
Conclusions The major conclusions of this work are that a consistent picture (18) Blair, J. T.; Westbrook, J. D.; Levy, R. M.; Krogh-Jespersen, K. Simple Models for Solvation Effects on Electronic Transition Energies: Formaldehyde and Water. Chem. Phys. Leu. 1989, 154, 531. Blair, J. T.; Krogh-Jespersen, K.; Levy, R. M. Solvent Effects on Optical Absorption Spectra: The 'A, 'A2 Transition of Formaldehyde in Water. J . Am. Chem. SOC.1989, 111, 6984. (19) The conventional notation for the amide vibrations indicates the mode of the N-deuterated form with a prime.
-
of the electronic excitation of NMA can be constructed that is in agreement with the absorption spectrum and the enhancement of Raman transitions in resonance but that this picture must be adjusted when considering this simple peptide in the vapor phase or in solution in non-hydrogen-bonding solvents relative to the particularly simple picture that appears to apply to aqueous solution. The change in geometry associated with electronic excitation appears, for the deuterated peptide in an aqueous environment, to be along a single mode of vibrational motion of the ground electronic state. The shape of the absorption spectrum, consistent with this simple picture, does not permit significant displacements along other degrees of freedom except possibly those with very low vibrational frequency. The possibilities within which readjustment of this simple picture must be made to accommodate the gas-phase and acetonitrile solution data are outlined. The most important variation is probably the change in the form of the ground electronic state normal modes with isotopic substitution and solvation.
Acknowledgment. The research was sponsored by NIH GM32626 and NSF CHE8816698. We thank Suzanne Hudson for writing the programs used to calculate the resonance Raman and absorption spectra. Registry No. NMA, 79-16-3; 14390-96-6; CHICN, 75-05-8.
D2,7782-39-0; 'IC,14762-74-4; "N,
Ultrafast Dynamics of a Quasi-Dissociative Diatomic Molecule in Solution. 2 S.-B. Zhu and G.W. Robinson* SubPicosecond and Quantum Radiation Laboratory, P.O.Box 4260, Texas Tech University, Lubbock, Texas 79409 (Received: July 3, 1990)
To answer the recent criticism of Keirstead and Wilson concerning our previous work, we perform a new molecular dynamics simulation of a quasi-dissociative diatomic molecule in the condensed phase. In accord with our earlier study, a noncanonical behavior of the low-mass test particle subjected to a strong systematic force is again observed. Deviations from the Maxwellian velocity distribution are found to be larger than the statistical uncertainty. This gives us confidence that the results are reliable and that they reflect the true nature of certain types of many-body systems.
I. Introduction Recently in this journal, Keirstead and Wilson' have published a paper entitled "A Breakdown of Equilibrium Statistical Mechanics?" In that paper they reproduce the conditions of an earlier calculation by Zhu and Robinson.2 They assert that the reported noncanonical behavior contradicts the predictions of classical equilibrium statistical mechanics, and by carefully examining the statistical errors, they conclude that the results in ref 2 can be attributed to poor statistics. While Keirstead and Wilson' criticized the breaking of the principles of statistical mechanics, they forget that the basic assumption in statistical mechanics is no more than an assumption of the "absence of interference between phenomena", that is, the assumption of "randomne~s",~This is the foundation for understanding any stochastic phen~menon.~Therefore, observed deviations from canonical behavior do not necessarily violate statistical mechanics. On the contrary, such deviations may provide an example where 'the development of statistical mechanics has reached the point of delving into its rig in".^ ~
(1) Keirstead. W. P.; Wilson,
K. R. J . Phys. Chem. 1990, 94, 918.
(2) Zhu,S.-B.;Robinson, G. W. J . Phys. Chem. 1989, 93, 164. (3) Ma, S.-K. Stutisricul Mechunics; World Scientific: Philadelphia, 1985. (4) Cramer, H.Murhemaricul Merhods of Srurisrics; Princeton University Press: Cambridge, 1946.
0022-365419112095-2967$02.50/0
The Maxwellian velocity distribution was initially derived for spatially homogeneous gases5 and later was found to be suitable for describing Brownian motion6 and a number of other physical In fact, over the past century, the MaxwellBoltzmann "normal distribution" has been successfully applied to so many problems that most persons, ignoring its origin and assumptions, are accustomed to treat it as an established law, or a "by definition" equilibrium property, rather than merely a mathematical formulation. The basic assumption for the probability distribution of a group of test particles to follow the Maxwell-Boltzmann formulation rests on the "central limit theorem", which requires that these tests particles are in contact with an infinite number of heat bath particles in the time duration of interest and further that the (5) Sears, F. W. An Introduction to Thermodynumics, The Kinetic Theory of Guses and Sratisricul Mechanics, 2nd ed.; Addison-Wesley: Reading, MA, 1953; Chapter 12. (6) Fowler, R.H.Srurisrical Mechanics, 2nd ed.; Cambridge University Press: Cambridge, 1936. (7) Chapman, S.;Cowling, T. G. The Muthemutical Theory of NonUniform Gases; Cambridge University Press: Cambridge, 1970. (8) Gray, C. G.; Gubbins, K.E. Theory of Molecular Fluids; Clarendon: Oxford, U. K., 1984; Vol. 1. (9) Tolman, R. C. The Principles of Srurisrical Mechunics; Clarendon: Oxford, U.K., 1938.
0 1991 American Chemical Society
2968 The Journal of Physical Chemistry, Vol. 95, No.8, 1991 motions of these test particles are mutually in depend en^^ These conditions are approximately satisfied for Brownian particles but are not for the small test particles studied in our work. The basic assumption of the Brownian motion theory is that both the size and mass of the test particles are considerably larger than that of the host particles.1° A number of concepts and rules break down when this assumption is no longer va1id.l' When the small test particles are furthermore subjected to strong systematic forces, their motions are not purely chaotic and independent. Each test particle communicates only with a few cage particles during the time scale of its state changes. Therefore, such a system may simply be beyond the range of applicability of the canonical en~emble.~ When the motion of a test particle is rapid compared with the thermal motion of its surroundings, the friction experienced by the test particle becomes nonlinearly dependent upon its velccity.12 This is because of the slow response of the heat bath.I3 In a nonlinear system, fluctuations can produce order.'* They can propagate unstable dynamics to other regimes and transfer information from short- to long-time behavior.l5 These effects cause a breakdown of the normal fluctuation-dissipation Such systems obviously deviate from Brownian motion behavior. However, in spite of their possible importance in chemical reactions and other microscopic rate phenomena, they have never been thoroughly studied. It is well-known that the configurational distribution of the solute does not follow the Boltzmann factor if one insists on using a "bare" intramolecular potential term.ls2J+21 In the Brownian motion case, the distribution is modified if a "static" contribution to the system-bath interaction is included.22 However, when a strong systematic force is involved and the motion becomes rapid, a dynamical contribution should also be taken into account.l'1*20Jl Since conformational rearrangements in liquids involve the displacement of solvent molecules, the ability of the solvent to respond to such changes plays an important role in determining the equilibrium configurational distribution of the solute. In view of all these complicating effects, in particular the possible mixing between the configurational and dynamical variables, how can Keirstead and Wilson' be certain that the velocity distribution for light test particles subjected to steep potentials remains of the same form as in a Brownian motion process? In ref 2, because of limited computational accessibility,we could only simulate systems with rather low barriers. As a consequence, deviations from the canonical distribution were not very significant. In addition, because of strong entropy effects, which tend to keep the molecule near its inner minimum in the quasi-dissociative diatomic molecule example, the statistical uncertainty is large. Because of these inherent problems, we later published a series of papers on other systems. Three of these papers investigated the velocity distribution functions of test particles subjected to intense corrugated force fields.23 A noncanonical behavior systematically emerged as the test particle mass or the forces to which it is subjected were systematically varied. We have also carried out several molecular dynamics studies of small particles (IO) Einstein, A. Ann. Phys. 1905, 17, 549. ( I I ) Zhu, S.-B.;Singh, S.; Robinson, G. W. Phys. Rev. 1989, A40, 1109.
(12) Zhu, S.-B.; Lee, J.; Robinson, G. W. Chem. Phys. 1990, 42, 3374. ( 1 3) Landau, L. D.; Lifshitz, E.M. Fluid Mechanics; Pergamon: London, 1959. (14) Suzuki, M. Adu. Chem. Phys. 1981,16, 195. (15) Coffey, W.T.; Evans, M.W.; Grigolini, P. Molecular Diffusion and Spectra; Wiley-lnterscience: New York, 1984; Chapters 7 and 8. (16) Ferrario, M.;Grigolini, P.; Leoncini, M.; Pardi, L.; Tani, A. Mol. Phys. 1984, 53, 1251. (17) Ferrario. M.; Grigolini, P.; Tani, A.; Vallauri, R.; Zambon, B. Adu. Chem. Phys. 1985.62, 225. (18) Zhu, S.-B. Phys. Rev. A 1990, A42, 3374. (19) Fixman, M. J. Chem. Phys. 1978,69, 1527. (20) Zhu, S.-8.;Robinson, G. W . Proc. 3rd Inr. Supercompuf. Inst. 1988,
I, 300. ( 2 1 ) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (22) Lindenberg, K.;Seshadri, V. Physica 1981, l09A, 483. Lindenberg, K.;Cortes, E. Physica 1984, 126A, 489. (23) Zhu, S.-B.;Lee, J.; Robinson, G. W. Chem. Phys. Lerr. 1989, 163, 328; 1990, 169, 355; 1990, 170, 368.
Zhu and Robinson TABLE I: Molecular Panmeters# r, = 1.91 8, A = 61.28 (kcal/mol)/A2 B = 315.00 (kcal/mol)/A4 cccl,/ks 3 2 1 K
uCCI,
= 5.881
A
ed/ke 91.5 K ud = 3.025 A
'The subscript d is used to represent the atoms of the diatomic molecule. Ud,CC14 = ( U C C I ~+ ud)/2. Cd.CCI, = (CCCI~CP. 45
4
I
I
I I
I I
1.1
1.5
2 .'3
1.9
2.1
Figure 1. Configurational distribution (solid curve) and the gas-phase intramolecular barriier potential (dashed curve), in units of 100 cal/mol. In all figures, the total height of the error bars represents two standard deviations.
in a heavy solidlike environment where the same type of results were obtained.24 Both of these calculations had much improved statistics. In a more direct response to the criticisms of Keirstead and Wilson' concerning the quasi-dissociative diatomic molecule problem,* in this work we report some new molecular dynamics simulations on this system. 11. Description of the Molecular System As in the earlier work? the system to be studied consists of one hypothetical quasi-dissociative diatomic molecule immersed in a bath containing 106 CCI4-like solvent atoms. The atoms interact through 12-6 Lennard-Jones potentials. In addition, the two atoms of the homonuclear diatomic molecule are linked by a symmetric double-well potential A B U(r) = Uo- -2( r - ro)' -4( r - r0)4
+
were Uo = 2.98 kcal/mol. Table I lists the other parameters employed for this potential function. The inner minimum corresponds to a bound state with a bond length of 1.47 A, while the outer minimum, which represents a quasi-dissociated state, lies at 2.35 A. In this simulation, the two minima have been separated from each other by a rather short distance so that large barrier forces can be produced without having to employ a computationally prohibitive barrier height. The diatomic molecule is also very light in mass, each of its atoms having a mass such that the ratio with the solvent mass is only 0.05. This mass differential, which is greater than in the earlier paper,2 will also help to promote the noncanonical behavior. In order to avoid possible computational errors induced when a strong barrier force is employed, we adopt a very short time increment of 0.1 11 fs for integrating Newton's equations of motion. After reaching equilibrium, 36 million time steps, corresponding to an elapsed time of 3.996 ns, are used to collect the data. This single long trajectory is different from one of the simulation methods of Keirstead and Wilson,l where velocities were rerandomized according to the Maxwell-Boltzmann distribution every 500 fs. This time duration is not long enough when compared with the rotational relaxation time of the diatomic molecule or with its lifetime in the inner minimum. The frequent perturbation of stopping and starting the computer run destroys the real nature (24) Zhu, S.-B.;Lee, J.; Robinson, G. W. Chem. Phys. Left. 1989, 191, 249; J . Fusion Energy 1990, 9, 465.
The Journal of Physical Chemistry, Vol. 95, No. 8, 1991 2969
Quasi-Dissociative Diatomic Molecule in Solution
’.-
I
T
h1.4
1.2 -
C
.-0 1.1 0 E 1.0 --
1.2
c,
Ir 0.9
0.9 -
i 1.3
1.6
1.9
2.2
2.5
0.8, 1.3
I
S
I
I
of the system and is by no means a “better simulation approach”. Determination of error bars in computer simulations is often a difficult problem. Here we adopt a method suggested by Allen and TildesleyZSto estimate standard deviations. The entire simulation is split into nine blocks. Each lasts for 444 ps and is viewed as an independent measure. In order to make certain that the block size is large enough so that the calculated standard deviation is approximately independent of the manner of the division, we recalculated the standard deviation by splitting the simulation into three blocks. The results from these two evaluations are very close to one another. 111. Numerical Results Figure 1 depicts the probability distribution of the bond length of the diatomic moelcule in the bath, compared with the free intramolecular potential energy function. The asymmetry in the condensed-phase population distribution with respect to the barrier top is a consequence of an entropy effect: the molecular state with the short bond length occupies a smaller hydrodynamical volume and consequently is much more stablesz6 Because of this strong entropy effect, the diatomic molecule tends to stay near its inner minimum most of the time. From time to time, it gains a sufficiently high kinetic energy along the reaction coordinate to cross the barrier into the outer minimum. However, because of rotations and other “nonproductive” degrees of freedom, the barrier crossing is always only minimally successful. This is another entropy effect, which tends to produce poor statistics for the more highly excited configurations in this type of problem. Because of the slow energy exchange rate between the low-mass diatomic test particles and the bath particles, the molecule, once excited and before returning to the more stable state, is unable to release all the extra kinetic energy it has gained. As a result, the kinetic energy of the molecular vibration does not uniformly distribute itself along the reaction coordinate. This is similar to the results reported in ref 24. Because of this, on average, the diatomic molecule with the longer bond length is hotter than the molecule with the shorter bond length. As plotted in Figure 2, this variation in the kinetic energy is much larger than the expected standard deviations. In conjunction with Figure 1, the position dependence of the vibrational kinetic energy in Figure 2 shows that the effective Hamiltonian of the diatomic molecule includes cross terms of its coordinate and momentum. This breaks down the condition required for the validity of the law of equipartition of en erg^,^,^ seemingly an important basis for the rejection of our earlier results by Keirstead and Wilson (see Appendix of ref I). A similar but less significant trend can be found for the nearest neighbors, defined as those solvent particles whose instantaneous distance from one of the solute atoms is less than 8.82 A. Under this assumption, the average number of nearest neighbors is 17. This
-
(25) Alien, M. P.; Tildesley, D. J. Compurer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (26) Flory, P. J. J . Chem. Phys. 1942, 10, 51.
,
9
I
9
2 i
2.2
Distance
Distance Figure 2. Configurational dependence of the reduced kinetic energy: asterisks, diatomic molecule (along the reaction coordinate); circles, nearest-neighbor solvent (see text).
T
1.9
1.6
Figure 3. Configurational dependence of S. 3.8
A.
2.61 1.3
I
8
I
I
I
I
8
1.9
1.6
8
I
2.2
8
8
5
Distance Figure 4. Configurational dependence of G.
I
Figure 5. Velocity distribution functions: solid curve, Maxwellian; solid symbols, diatomic molecule (along the reaction coordinate); circles,
nearest-neighbor solvent (see text). definition is somewhat arbitrary, and the curve in Figure 2 certainly relies upon this definition. Local departure from the Gaussian behavior can be characterized by the ratios G = ( U ~ ) ~ / ( U and ~ ) S = n+/n-, where n+ and n- represent respectively the numbers of diatomic molecules having positive (tending to increase the bond length) and negative velocities. Figure 3 depicts the position dependence of S. In accordance with Figure 9 of ref 2, the ratio oscillates, following the variation of the barrier force. Illustrated in Figure 4 is the variation of G along the reaction coordinate. This ratio equals 3.0 for a purely Gaussian velocity distribution. Since the ratio contains the fourth-order moment, which is strongly dependent on the sparsely populated high-energy tail of the distribution function, the error bars are generally much larger than those measured for other properties. Nevertheless, the variations of G are consistent with a velocity distribution function favoring slower particles in the inner well than in the outer well. In Figure 5 we compare the distribution functions for (1) the velocity component of the diatomic molecule along the reaction coordinate, (2) the velocity of the nearest solvent atoms (
