J. Phys. Chem. 1996, 100, 9495-9499
9495
Molecular Dynamics Simulations of Vibrational Energy Distribution in Vibrational Cooling and Heating Hackjin Kim* Department of Chemistry, College of Natural Sciences, Chungnam National UniVersity, Taejon, 305-764, Korea
Youngdo Won Department of Chemistry, College of Natural Sciences, Hanyang UniVersity, Seoul, 133-791, Korea ReceiVed: October 2, 1995; In Final Form: February 9, 1996X
Temporal evolution of vibrational energy distribution during vibrational cooling and heating in a molecular cluster is simulated using techniques of molecular dynamics. The cluster is composed of hundreds of naphthalene molecules with fused atoms for C-H groups. In vibrational cooling, a hot molecule (300 K) is allowed to cool in a cold cluster (10 K). In vibrational heating, a cold molecule (60 K) is heated in a hot cluster (300 K). Normal vibrations of the model molecule are analyzed, and instantaneous amplitudes of the vibrational modes are calculated from trajectories. Vibrational energy transfer processes are discussed. In vibrational cooling, the out-of-plane vibrations play important roles in intramolecular transfer processes and the optically active in-plane vibrations are mainly involved in intermolecular resonant transfer processes. Intramolecular transfer processes through out-of-plane vibrations are significant in vibrational heating.
Introduction Molecular vibrations play important roles in various dynamics of condensed phases through interactions with phonons or other molecular motions. Thermal reactions are initiated by the excitation of molecular vibrations and the vibrational modes work as the reaction coordinates.1 Most of the reaction energy, which is initially located in the vibrational modes of product molecules, transfers to nearby molecules to provoke the spatial propagation of chemical reactions in condensed phases.2 Shockinduced chemistry of energetic materials results from shock energy transfer to phonons to molecular vibrations.3-5 Vibrational relaxation of a molecule adsorbed on a surface causes thermal desorption of molecules through energy transfer from vibrations to phonons.6 Molecular vibrations also serve as intermediate states in many processes such as laser ablation, in which electronic changes into translational energy.7 However, microscopic understanding of vibrational energy transfers in condensed phases is still poor. Knowledge of vibrational dynamics is helpful to the investigation of complicated processes in real systems such as energetic materials, molecular crystals, polymers, or liquids. Molecular dynamics (MD) techniques, which have been widely used to study the microscopic structures of macromolecules,8 have not been employed extensively to investigate dynamic processes. Good potential surfaces are required for realistic simulations of dynamic processes. Some examples of dynamic processes simulated using MD techniques have been reported: propagation of heat waves in crystalline Fe,9 vibrational cooling of hot porphyrin embedded in a protein,10 detonation chemistry of a molecular solid at nanoscale resolution,11 evolution of hot spots in rapidly compressed crystals,12 and vibrational cooling, heating, and nanoscale heat conduction in crystalline molecular clusters.13,14 We use MD techniques to simulate the vibrational energy distribution inside a molecule during mechanical energy transfer in molecular clusters: vibrational cooling (VC) of a hot molecule (300 K) in a cold cluster (10 K), and vibrational X
Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(95)02916-9 CCC: $12.00
heating (VH) of a cold molecule (60 K) in a hot cluster (300 K). The dynamics of VC and VH are not simple reverse processes because the cluster temperatures differ. The energy distribution among the molecular vibrations is analyzed by computing the vibrational amplitude of each mode from atomic trajectories. VC and VH in the molecular clusters were previously investigated by using the same techniques.14 However, in those studies, only total vibrational energy was calculated by subtracting translational and librational energy from total kinetic energy without analyzing the distribution of energy among the vibrational modes. Evolution of the vibrational energy distribution in VC and VH yields the detailed information about the features of the vibrational dynamics in condensed phases. Computational Methods Details of the model system and the computational procedure, except the calculation of the vibrational mode amplitude, have been discussed previously.13,14 A brief description of the system is given below. A schematic diagram of the model cluster system is shown in Figure 1. The cluster with radius ∼20 Å contains 235 naphthalene molecules which are modeled with fused atoms for C-H groups. Since the C-H stretch modes have frequencies higher than 3000 cm-1, the use of fused atoms should not affect the vibrational dynamics unless the cluster temperature is too high. The initial coordinates of the cluster are taken from the experimental crystal structure of naphthalene.15 Close contacts in van der Waals radii are removed by 30 steps of ABNR energy minimization. A boundary potential,8,14 given by a function of the distance from the origin, is employed to maintain the structure of the cluster at a high temperature. Details of the potential parameters are given in refs 13 and 14. The Chemistry at Harvard Macromolecular Mechanics (CHARMM) version 22g416 is used for the computation of atomic trajectories, the normal mode analysis of the model molecule, and the calculation of the mode amplitudes. In the VC simulations, the energy-minimized cluster is heated to 10 K and then equilibrated for 20 ps. The center molecule is then heated and equilibrated at 300 K with all other molecules © 1996 American Chemical Society
9496 J. Phys. Chem., Vol. 100, No. 22, 1996
Kim and Won TABLE 1: Normal Vibrations of the Model Moleculea symmetry
frequency (cm-1)b
Ag
489 (1) 820 (2) 1141 (3) 1539 (4) 1713 (4) 560 (1) 857 (2) 1258 (3) 1755 (4)
symmetry
frequency (cm-1)b
In Plane
B1g
Au B2g Figure 1. Schematic diagram of the model cluster. Inside the reaction zone of normal dynamics, 235 naphthalene molecules are contained. The molecules at the perimeter between 16 and 20 Å are subject to scaled Langevin dynamics and an infinite boundary potential is imposed at 20 Å. The molecule at the center of the cluster is initially hotter (300 K) than the surroundings (10 K) in VC. In VH, the center molecule is initially cooler (60 K) than the surrounding (300 K).
284 (1) 1089 (3) 521 (1)
B2u
B3u
Out of Plane B3g B1u
420 (1) 746 (2) 1261 (3) 1420 (4) 708 (2) 1063 (3) 1440 (4) 1609 (4)
911 (2) 1322 (3) 143 (1) 688 (2)
a The modes of Bu symmetries are optically active. b The numbers in the parentheses indicate the mode group of Figure 6.
kept fixed. Constraints are then lifted and the center molecule is allowed to cool to the final temperature of ∼12 K. In the VH simulations, the center molecule is kept at 60 K and the other molecules are heated and equilibrated at 300 K. After the constraints on the center molecule are removed, it is allowed to heat up to the final temperature, which is close to 300 K. The boundary potential shown in Figure 1 is kept during the dynamics of both VC and VH. Translational, librational, and vibrational energies of the center molecule are calculated from the amplitudes of the normal modes, which are determined from atomic trajectories and results of the normal-mode analysis. Since atomic trajectories are computed in Cartesian coordinates, the normal-mode analysis by the standard FG matrix method17 is carried out in Cartesian coordinates instead of the internal coordinates. Details of the mathematical scheme for calculation of the mode amplitudes are given in the appendix. Results and Discussion Normal Vibration Analysis. Like the real naphthalene molecule, the model naphthalene molecule of fused atoms belongs to D2h group. The symmetries of the 24 vibrational modes of the model molecule are 5Ag + 4B1g + 2B2g + 2B3g + 2Au + 2B1u + 4B2u + 4B3u. The results of the normal vibration analysis are summarized in Table 1. The vibrational frequencies of the model molecule agree with the ones previously determined by a Fourier transform of a velocity autocorrelation function.14 The in-plane vibrational modes have the symmetries of Ag, B1g, B2u, and B3u, and the out-of-plane modes have the symmetries of Au, B2g, B3g, and B1u. The modes of the Bu symmetries are optically active. For the real naphthalene molecule, the in-plane vibrations have the symmetries of Ag, B3g, B1u, and B2u, and the out-of-plane modes have the symmetries of Au, B1g, B2g, and B3u.18 The difference in the mode symmetry between the model and real molecule results from the lack of wagging motions of H atoms in the model molecule. Evolution of Kinetic Energy in VC and VH. The amplitudes of all 30 modes including translations and librations are calculated from atomic trajectories and normal-mode analysis. Figure 2 shows the kinetic terms of translational, librational, and vibrational motions, determined from the mode amplitudes for the center molecule in VC and VH. Data presented in this work are the average of five runs. In Figure 2, the kinetic energy
Figure 2. Kinetic amplitude change of translational, librational, and vibrational motions of the center molecule in (a) VC and (b) VH.
changes just as in the previous simulations14 where translational, librational, and total kinetic energies were calculated using the atomic trajectories. Identical evolution of the kinetic energy confirms that the two calculation methods for kinetic energy are equivalent. Slower equilibration of the vibrations than translations and librations is well-known in the dynamics of molecular crystals.2,13,14,19 The initial vibrational energy of the center molecule in VH shown in Figure 2b indicates that there was energy leaking from the cluster to the center molecule during the equilibration period for VH. The intended temperature for the center molecule was 10 K but the initial temperature of the molecule is about 60 K. However, the initial temperature of the molecule is not critical in the study of mechanical energy dynamics. Evolution of the Vibrational Energy Distribution in VC and VH. VC and VH are composed of many individual stateto-state vibrational relaxation events mediated by cubic anharmonicity.20 Energy-transfer processes can be classified into three types: single-site, intermolecular resonant, and intermolecular nonresonant, as shown in Figure 3.14 The singlesite processes involve energy transfer between vibrations on the same molecule, or energy transfer between a vibration and
Simulations of Vibrational Energy Distribution
Figure 3. Schematic diagram of vibrational energy transfer. The hatched region denotes collective phonons; the horizontal lines, internal vibrations. The center molecule of the cluster is initially set to different temperature from the surrounding. Processes 1 and 3 involve phononassisted intra- and intermolecular transfer, which accelerate at high temperature. Process 2 involves two-site resonant energy transfer, which deactivates due to modulation of energy levels by phonons at high temperature.
Figure 4. Evolution of the vibrational energy distribution of the center molecule in (a) VC and (b) VH. The average of 1 ps data around the specified time is plotted. The atomic trajectories are collected at intervals of 0.05 ps. At equilibrium, the kinetic energy of all vibrations should be equal in the classical limit.
phonons at a single site. The intermolecular resonant processes require neighboring identical molecules. The intermolecular nonresonant processes involve a phonon-assisted transfer to a vibration on another molecule. The rates of intermolecular resonant processes decrease at high temperature due to modulation of energy levels by phonons, while the rates of the other phonon-assisted processes increase at high temperature. The rates of the intermolecular resonant processes decrease by nearly 2 orders of magnitude between 10 and 300 K.21 Figure 4 shows the evolution of the vibrational energy distribution of the center molecule in VC and VH. The intensity of the spectra of Figure 4 represents the kinetic energy of the mode. The spectra are the average of 1 ps data around the specified time. The atomic trajectories were collected at intervals of 0.05 ps. If only the single-site processes were effective, the high-frequency modes would lose energy faster than the low-frequency modes in VC and would be heated after the low-frequency modes in VH. The greater energy content of the lower frequency modes renders the energy distribution asymmetric as in a Boltzmann distribution. If only the intermolecular processes were effective and the rates of energy
J. Phys. Chem., Vol. 100, No. 22, 1996 9497
Figure 5. The standard deviation of 24 mode energies of the center molecule during VC and VH. In VC, the standard deviation decreases with the progress of dynamics. The constant deviation in VH indicates that the vibrations are in thermal equilibrium.
transfer were independent of the mode frequencies, the energy distribution would change while maintaining the pattern of equilibrium. If the single-site and the intermolecular transfer processes were equally effective, the rates of the energy changes would depend more or less on the mode frequency, and the distribution of energy would change in a complex manner. However, it is difficult to study the mechanisms of mechanical energy transfer directly from the evolution of the energy distribution shown in Figure 4 because of the fluctuations of mode energy. The effect of energy fluctuations should be considered in the analysis of the evolution of the energy distribution. The changes of the standard deviation of 24 vibrational mode energies during VC and VH are plotted in Figure 5. The standard deviation is almost constant in VH while the deviation decreases with the progress of VC. A similar trend of deviations is observed in Figure 2 where the deviation of five different runs is marked. The fluctuations of mode energy which are related to vibrational energy as well as phonon energy, increase at higher temperature. In the early stage of VC, the fluctuations are large because the center molecule is hot. As the center molecule cools in VC, the fluctuations decrease. In VH, translational and librational motions of the center molecule reach equilibrium with the hot cluster molecules within a few picoseconds. Therefore, fluctuations at high temperatures are observed in VH. The constant deviation of the vibrational mode energies implies that the vibrations of the center molecule are in thermal equilibrium in VH. Figure 6 shows the evolution of the scaled energy of the mode groups (6 modes) in different frequency regions: 143-600 cm-1, 600-1000 cm-1, 1000-1400 cm-1, and 1400-1755 cm-1. The scaled energy of the mode group i, θgroup,i(t) is calculated by
θgroup,i(t) ) [Agroup,i(t) - Aall(t)]/Sall(t)
(1)
where A and S denote the average and the standard deviation of the vibrational mode energies, and the subscript of “all” corresponds to the values for all 24 modes. When the average energy of a mode group is larger than the average of 24 modes, the scaled energy is positive and the deviation from the average is expressed in terms of the standard deviation of all mode energies. The initial energy spread in VC as shown in Figure 6a indicates that the vibrations of the center molecule are in equilibrium at the onset of the dynamics. Localization of energy in the mode group 1 and less energy content in the mode groups of higher frequencies suggest that the single-site processes are dominant and that the intermolecular processes are inefficient in VC. The energy accumulation into the low-frequency modes by the single-site processes diminishes with the progress of VC
9498 J. Phys. Chem., Vol. 100, No. 22, 1996
Figure 6. Evolution of the scaled energy of the vibrational modes in four different frequency regions in (a) VC and (b) VH. The vibrational modes are grouped in the order of vibrational frequency: 143-600 cm-1 (1), 600-1000 cm-1 (2), 1000-1400 cm-1 (3), and 1400-1800 cm-1 (4).
because the vibrational relaxation rate is proportional to the molecular temperature in the classical limit. As mentioned in the previous section, there was energy leaking from the cluster to the center molecule during the equilibration period for VH. Although the initial energy spread of VH as shown in Figure 6b looks as if the vibrations of the center molecule were not in thermal equilibrium, they are in equilibrium as shown in Figure 5. The difference between VC and VH arises mainly from the cluster temperature, that is, the phonon population. Note that localization of the energy at the early stage of VH is not in the order of the mode frequency. It is not clear why the modes of group 2 are heated preferentially. However, these results indicate that the intermolecular nonresonant processes are more effective than the single-site processes in VH. The intermolecular resonant processes almost stop at high temperature.14 The single-site processes which accelerate at a high temperature are thought to contribute to the redistribution of the vibrational energy rather than the net heating of the vibrational modes. Mode Characteristics in Vibrational Energy Transfer. The mode activities in the energy transfer are related to the relaxation times of the modes which depend on the involved anharmonicities, the density of states, and the phonon population.22 The vibrational relaxation times cannot be directly employed for the interpretation of the simulation results; however, the overall vibrational cooling processes are considerably affected by the vibrational relaxation times.19,23 It has been experimentally observed for the real naphthalene crystal that the Ag modes usually have much longer relaxation times than the out-of-plane modes and the optically active in-plane modes.24 Evolution of the scaled energy of four mode groups of different characters is shown in Figure 7. The 24 vibrational modes are classified into Ag vibrations (five modes), optically inactive in-plane vibrations (four modes), optically active inplane vibrations (eight modes), and out-of-plane vibrations (seven modes). The mode frequencies of each group are widely distributed as given in Table 1. Equation 1 is used to calculate
Kim and Won
Figure 7. Changes of the scaled energy of the vibrational modes of different characters in (a) VC and (b) VH. The same scaling process as in Figure 6 is applied. The vibrational modes are classified into four groups: Ag vibrations (Ag, five modes), optically active in-plane vibrations (act, eight modes), optically inactive in-plane vibrations (inac four modes), and out-of-plane vibrations (out, seven modes).
the scaled energy of each group of Figure 7. Classification of the mode groups based on the mode character shows more localization of energy rather than classification based on the mode frequency. The scaled energy of the Ag modes which maintain their average without large fluctuations in VC is explained by weak interactions of the Ag modes with other vibrations. The relaxation times of the out-of-plane modes and the optically active in-plane modes are similar24 but the energy contents of these modes are quite different in VC. The outof-plane modes have energy greater than the average, but the optically active in-plane modes do not. These results suggest that the out-of-plane modes are mainly involved in the singlesite processes and that the optically active in-plane modes lose energy through the intermolecular transfer processes. The vibrational modes mainly involved in the intermolecular processes have less energy than the modes primarily involved in the single-site processes in VC. Dipole-dipole interactions of the optically active in-plane modes may be considered to enhance the efficiencies of these modes for intermolecular resonant processes, which are effective in VC.14 At a high temperature, a number of phonons reduce the effect of anharmonicity on the vibrational relaxation time; therefore, the mode characters related to the vibrational relaxation time becomes less important in VH. The effect of a large phonon population is confirmed from the behaviors of the Ag modes. Although the Ag modes have relatively small anharmonicities, the energy of the Ag modes changes in a manner similar to that of the in-plane modes. Localization of energy in the out-ofplane modes, which is quite remarkable and is maintained for a long period in VH, suggests that these modes are involved in different energy-transfer processes in VH from those of other modes. As in VC where the out-of-plane modes are primarily involved in the single-site processes, these modes are also heated by the single-site processes in VH. The significance of singlesite processes in VH agrees with the experimental results on vibrational heating in condensed molecular systems.25
Simulations of Vibrational Energy Distribution
J. Phys. Chem., Vol. 100, No. 22, 1996 9499
Acknowledgment. We appreciate the helpful comments of the referees. This work was supported by grants from the Cray R & D program and the Korean Science and Engineering Foundation.
of the mode amplitudes. The molecule is translated and rotated so that the principal axes of the molecule are coincident with the reference molecule. The amplitudes of normal modes are then computed by transforming the atomic trajectories.
Appendix
References and Notes
The normal-mode analysis using the FG matrix method gives the vibrational frequencies, the actual forms of normal vibrations, and the transform matrix connecting the normal vibrations and the general molecular coordinates. The matrix of the normal coordinates including translations and rotations, Q, is related to the general coordinate, q:
(1) McDonald, J. D. Annu. ReV. Phys. Chem. 1979, 30, 29. (2) (a) Dlott, D. D.; Fayer, M. D. J. Chem. Phys. 1990, 92, 3798. (b) Tokamakoff, A.; Fayer, M. D.; Dlott, D. D. J. Phys. Chem. 1993, 97, 1901. (3) Fried, L. E.; Ruggiero, A. L. J. Phys. Chem. 1994, 98, 9786. (4) Sakata, J.; Wight, C. A. J. Phys. Chem. 1995, 99, 6548. (5) (a) Lee, I.-Y. S.; Hill, J. R.; Dlott, D. D. J. Appl. Phys. 1994, 74, 4975. (b) Chen, S.; Hong, X.; Hill, J. R.; Dlott, D. D. J. Phys. Chem. 1995, 99, 4525. (6) Benjamin, I.; Reinhart, W. P. J. Chem. Phys. 1989, 90, 7535. (7) (a) Kim, H.; Zyung, T.; Postlewaite, J. P.; Dlott, D. D. Appl. Phys. Lett. 1989, 54, 2274. (b) Wen, X.; Tolbert, W. A.; Dlott, D. D. J. Chem. Phys. 1993, 99, 4140. (c) Hare, D. E.; Dlott, D. D. Appl. Phys. Lett. 1994, 64, 715. (8) Brooks, C. L., III; Karplus, M. L.; Pettitt, B. M. Proteins: A Theoretical PerspectiVe of Dynamics, Structure and Thermodynamics; Adv. Chem. Phys. Vol. LXXI; Wiley: New York, 1988. (9) Tsai, D. H.; Trevino, S. F. J. Chem. Phys. 1983, 79, 1684. (10) Henry, E. R.; Eaton, W. A.; Hochstrasser, R. M. Porc. Natl. Acad. Sci. U.S.A. 1986, 83, 8982. (11) Brenner, D. W.; Robertson, D. H.; Elert, M. L.; White, C. T. Phys. ReV. Lett. 1993, 70, 2174. (12) Tsai, D. H.; Armstrong, R. W. J. Phys. Chem. 1994, 98, 10997. (13) Kim, H.; Dlott, D. D. J. Chem. Phys. 1990, 94, 8203. (14) Kim, H.; Dlott, D. D.; Won, Y. J. Chem. Phys. 1995, 102, 5480. (15) Cruickshank, D. W. Acta. Crystallogr. 1957, 10, 504. (16) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Kaplus, M. J. Comput. Chem. 1983, 4, 187. (17) Wilson, E. B., Jr.; Decius, J. C.; Cross, P. C. Molecular Vibrations; McGraw-Hill: New York, 1955. (18) Sellers, H.; Pulay, P.; Boggs, J. E. J. Am. Chem. Soc. 1985, 107, 6487. (19) Kim, H.; Dlott, D. D. J. Chem. Phys. 1990, 93, 1695. (20) (a) Velsko, S. P.; Hochstrasser, R. M. J. Chim. Phys. (Paris) 1985, 82, 153. (b) Velsko, S. P.; Hochstrasser, R. M. J. Phys. Chem. 1985, 89, 2240. (21) Grover, M. K.; Silbey, R. J. Chem. Phys. 1970, 52, 2099. (22) Califano, S.; Schettino, V.; Neto, N. Lattice Dynamics of Molecular Crystals; Springer: Berlin, 1981. (23) (a) Hill, J. R.; Dlott, D. D. J. Chem. Phys. 1988, 89, 830. (b) Hill, J. R.; Dlott, D. D. J. Chem. Phys. 1988, 89, 842. (24) (a) Schosser, C. L.; Dlott, D. D. J. Chem. Phys. 1984, 80, 1394. (b) Hill, J. R.; Chronister, E. L.; Chang, T.-C.; Kim, H.; Postlewaite, J. C.; Dlott, D. D. J. Chem. Phys. 1988, 88, 949. (c) Hill, J. R.; Chronister, E. L.; Chang, T.-C.; Kim, H.; Postlewaite, J. C.; Dlott, D. D. J. Chem. Phys. 1988, 88, 2361. (25) Chen, S.; Lee, I.-Y. S.; Tolbert, W. A.; Wen, X.; Dlott, D. D. J. Phys. Chem. 1992, 96, 7178.
Q ) Lq and
dQ/dt ) Lv where Q an q are matrices of (3n × 1) and L is the transform matrix of (3n × 3n). Here n is the number of atoms in the molecule. The matrix of (3n × 1) v represents the atomic velocities in the general coordinate system. When the unit vectors are used for the general coordinates, the resultant normal coordinates matrix and the transform matrix are normalized. The atomic trajectories, qobs and vobs, calculated in the MD simulation, are related to the amplitudes of the normal coordinates. That is
A Q ) Lqobs and
A dQ/dt ) Lvobs where A is a diagonal matrix of (3n × 3n) representing the amplitudes of the normal coordinates. The kinetic energy corresponds to [A dQ/dt]2/2. Since the transform matrix is determined by the normal-mode analysis, the amplitudes of the vibrational modes are calculated by transformation of the atomic trajectories. One problem in applying the above method to the calculation of mode amplitudes is that dynamically moving molecules are randomly oriented in space. Therefore, the molecule under analysis should be adjusted into the reference axis system which is employed for the normal-mode analysis before the calculation
JP952916B
