Chemical Dynamics Simulations of Intermolecular Energy Transfer

Department of Chemistry and Biochemistry, Texas Tech University, Lubbock, Texas 79409, United States. J. Phys. Chem. A , 2016, 120 (27), pp 5187–519...
1 downloads 14 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCA

Chemical Dynamics Simulations of Intermolecular Energy Transfer: Azulene + N2 Collisions Hyunsik Kim, Amit K. Paul, Subha Pratihar, and William L. Hase* Department of Chemistry and Biochemistry, Texas Tech University, Lubbock, Texas 79409, United States S Supporting Information *

ABSTRACT: Chemical dynamics simulations were performed to investigate collisional energy transfer from highly vibrationally excited azulene (Az*) in a N2 bath. The intermolecular potential between Az and N2, used for the simulations, was determined from MP2/6-31+G* ab initio calculations. Az* is prepared with an 87.5 kcal/mol excitation energy by using quantum microcanonical sampling, including its 95.7 kcal/mol zero-point energy. The average energy of Az* versus time, obtained from the simulations, shows different rates of Az* deactivation depending on the N2 bath density. Using the N2 bath density and Lennard-Jones collision number, the average energy transfer per collision ⟨ΔEc⟩ was obtained for Az* as it is collisionally relaxed. By comparing ⟨ΔEc⟩ versus the bath density, the single collision limiting density was found for energy transfer. The resulting ⟨ΔEc⟩, for an 87.5 kcal/mol excitation energy, is 0.30 ± 0.01 and 0.32 ± 0.01 kcal/mol for harmonic and anharmonic Az potentials, respectively. For comparison, the experimental value is 0.57 ± 0.11 kcal/mol. During Az* relaxation there is no appreciable energy transfer to Az translation and rotation, and the energy transfer is to the N2 bath.

I. INTRODUCTION Understanding collisional intermolecular energy transfer from highly excited molecules is important in unimolecular processes.1−3 A number of experimental studies of energy transfer from highly vibrationally excited molecules have been performed. CS2, SO2, cycloheptatrienes, and toluene were studied by time-resolved UV absorption (UVA) spectroscopy.4−7 Pyrazine with various colliders was studied by infrared fluorescence (IRF).8,9 NO2 was studied by resolved Fourier transform infrared (FTIR) emission spectroscopy.10 Benzene was studied with IRF.11−14 Azulene was studied with IRF,15−19 UVA,20 and time-resolved diode laser absorption spectroscopy.21 Also, to investigate potential energy surface, collisional, and statistical properties for unimolecular collisions, many simulations have been performed for benzene,22−27 toluene,7,28,29 hexafluorobenzene,30−32 and azulene.33−40 Azulene (Az) is ideal molecule to study energy transfer processes experimentally, because it is capable of a one-photon transition16 that produces unambiguous initial energies that may easily be varied. It eliminates complications by not reacting with many bath gases and has very clean photophysics, with very low quantum yields for forming naphthalene or an Az triplet state.16 Experimentally, highly vibrational excited Az is generated by 193−571 nm laser excitation.15−21,41,42 It is first excited electronically to S2 or S1 and then undergoes rapid internal conversion to form highly vibrationally excited ground state Az, that is, S0* or Az*. Intermolecular energy transfer (IET) from Az* is then studied for collisions with bath molecules. A statistical model has been used to interpret this IET.43 © XXXX American Chemical Society

Carefully characterized Az* + N2 IET experiments are those by Hippler and co-workers20 for a laser excitation of 337 nm (30 600 cm−1). In the work presented here, chemical dynamics simulations were used to compare with these experiments. The simulations are performed, as for previous recent studies,31,32,44 using a bath of N2 molecules to average over the random, thermal properties of the N2 + Az* collisions. Different bath densities are considered, and the density is lowered to reach the single collision limit. Both harmonic and anharmonic intramolecular potentials are used for Az* to see how they affect the IET dynamics. Electronic structure calculations were performed to determine the Az* + N2 intermolecular potential, and the effect of this potential on the IET is considered. Of particular interest is the average energy transferred per collision between Az* and N2 in the low-density single collision limit and how it compares with experiment.20

II. POTENTIAL ENERGY To perform an accurate chemical dynamics simulation, it is essential to use a proper potential energy function. The total potential energy for the Az* + N2 bath system is written as Vtot = VAz + VN2 + VAz − N2 + VN2 − N2

(1)

Special Issue: Piergiorgio Casavecchia and Antonio Lagana Festschrift Received: January 27, 2016 Revised: May 13, 2016

A

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

harmonic bends the force constant is the same, and fθ = 0.572 mdyn Å rad −2. The equilibrium values in degrees for the angles θe are C2−C1−H1 = 124.25; C2−C1−C10 = 111.5; C3−C2− H2 = 126.45; C1−C2−C3 = 107.10; C2−C3−C4 = 126.0; C2−C3−C9 = 107.2; C4−C3−C9 = 126.8; C3−C4−C5 = 130.1; C3−C4−H3 = 114.95; C4−C5−C6 = 128.0; C4−C5− H4 = 116.0; C5−C6−C7 = 130.2, and C5−C6−H5 = 114.9. The parameters for harmonic α-bending/wagging are fα = 0.450 mdyn Å rad−2 and αe = 0°; and values for torsion angle: barrier (Vtor in kcal/mol) are C3−C2−H2 vs C10−C1−H1:16.192, C1−C10−H8 vs C3−C9−C8:17.644, C3−C9−C10 vs C7− C8−H7:16.644, C9−C8−H7 vs C6−C7−H6:17.541, C8− C7−H6 vs C5−C6−H5:14.843 and C2−C3−C4 vs C10− C9−C8:59.891. To include anharmonicity in the azulene intramolecular potential, the harmonic stretch potential was replaced by the Morse stretch potential [V(R) = De(1 − eβe(R − Re))2] for all the CC and CH bonds. The bond-energy-bond-order relation (BEBO)45 was used to determine bond energies (De) and beta (βe) parameters for the CC bonds. The De (kcal/mol) and βe (Å−1) are 112.5, 2.32 for C1−C2; 88.125, 2.57 for C2−C3; 128.75, 1.85 for C3−C9; 88.125, 2.63 for C3−C4; 120.625, 2.25 for C4−C5; and 104.375, 2.40 for C5−C6. The values De = 112.2 kcal/mol and βe = 1.79 Å−1 for the C−H bonds are similar to those for benzene.46 Stretch−bend anharmonicity was added to the potential energy surface by attenuating the potentials for the H−C−C and C−C−C harmonic bends, each given by V(θ) = fθ(θ − θ0)2/2 where fθ = fθ0 S(rj) S(rk).47 Here, rj and rk define the bend. The attenuation S(r) is defined as exp(−CΔr2),43 with C = 0.825 Å−2 a standard value,47 Δr = (r − re), and S(r) equal to unity if r < re. The Morse intramolecular potential for N2 is obtained from the literature and the same as used previously.31,32 The parameters are De = 228.3 kcal/mol, βe = 2.699 Å−1, and Re = 1.0945 Å.

where Vtot is total potential energy, VAz and VN2 are the intramolecular potentials for Az and N2, respectively, VAz−N2 is the intermolecular potential between Az and N2, whereas VN2−N2 is the intermolecular potential between N2 molecules. In the following, details of this potential are described. A. Intramolecular Potentials. The minimum energy geometry40 of Az is shown in Figure 1. The harmonic Az

Figure 1. Minimum-energy geometry of the Az molecule with atom numbering.

intramolecular potential is obtained from a previous work done by Bernshtein and Oref.40 The harmonic stretching force constants and equilibrium bond lengths ( f R in mdyn Å−1, Re in Å) are 8.433 50, 1.404 for C1−C2; 8.100 94, 1.404 for C2−C3; 6.132 68, 1.490 for C3−C9; 8.501 16, 1.380 for C3−C4; 8.501 16, 1.391 for C4−C5; 8.332 73, 1.385 for C5−C6; and 5.0, 1.084 for C−H stretches. For the H−C−C and C−C−C

Figure 2. Intermolecular potential energies of Az−N2 for eight different orientations. The red ● and blue ■ are for the MP2/6-31+G* and MP2/6311+G* calculations, respectively. The solid black line is the fit to the MP2/6-31+G* calculations using the two-body potential in eq 2. B

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. Parameters for the Azulene-N2 Intermolecular Potential Fit to ab Initio Calculationsa interactionb N−C N−H

A

B

C

D

n

m

32267.81 21476.80 10692.68 18694.06

3.291957 3.229249 3.403010 3.705687

−1864.199 −1411.447 −2693.609 −2591.526

667.9112 13298.14 5981.596 4584.946

7 7 8 8

10 10 11 11

The units of A, B, C, D are kcal/mol, Å−1, kcal Ån/mol, and kcal Åm/mol, respectively. The top and bottom rows of parameters are those for the fits to the MP2/6-31+G* and SCS-MP2/6-311++G** potential energy curves, respectively. bFit with eq 2. a

Figure 3. Intermolecular potential energies of Az−N2 for eight different orientations. The red ● and blue ■ are for the MP2/6-31+G* and SCSMP2/6-311++G** calculations, respectively. The solid black line is the fit to the SCS-MP2/6-311++G** calculations using the two-body potential in eq 2.

fine details of the intermolecular potential, such as potential energy minima for different collisional orientations. In previous work,46,49,50 MP2 was found to give accurate intermolecular potentials for Ar−CF4 and benzene−Na+, and MP2 with the 631+G* and 6-311+G* basis sets were used to calculate Az−N2 potential energy curves for the eight different orientations. The results of these calculations are given in Figure 2, and the two basis sets give similar potential energy curves. All of the MP2/ 6-31+G* points for the eight orientations were simultaneously fit using a combined genetic/nonlinear least-squares algorithm.31,51 The potential energy curves were fit by a sum of N− C and N−H two-body potentials, each given by the modified Buckingham equation

B. Intermolecular Potential. The intermolecular N2−N2 potential was described in ref 31 and is not discussed here. To obtain an intermolecular potential for the Az−N2 interaction, ab initio calculations were performed for eight different orientations of this interaction as depicted in Figure 2. In panels 2a−2d, the N2 molecule approaches the Az plane perpendicularly, directly above a C atom; that is, atoms C2, C3, C4, and C5 in panels 2a, 2b, 2c, and 2d, respectively (see Figure 1 for atom numbers). R on the x-axis is the distance between the C atom and the nearest N atom. For panels 2e and 2f the N2 molecule approaches Az, in-plane along the H−C bond axis of the H5-atom. In 2e the N2 molecule and the H−C bond are collinear, and in 2f the midpoint of the N2 molecule approaches the H−C bond perpendicularly. In 2e and 2f, R is the distance between the H atom and the nearest N atom and the midpoint of the N2 bond, respectively. Panels 2g and 2h are both out-ofplane approaches of N2 toward a C−C bond in Az. Panel 2g is an approach of the midpoint of N2 toward the midpoint of the C5−C6 bond, with the C5−C6 and N−N bonds perpendicular. Finally, panel 2h is a perpendicular approach of the midpoint of the N2 molecule toward the midpoint of the C6−C7 bond, with the C6−C7 and N−N bonds parallel. Though quite accurate intermolecular potentials may be obtained with the computationally demanding CCSD(T) theory,48,49 previous simulations for C6F6 + N2 collisions31 illustrated that the IET dynamics are not strongly sensitive to

V = A e−Br +

C D n + m r r

(2)

Figure 2 shows that this model provides an excellent fit to the ab initio points. The fitting parameters are given in Table 1. The respective MP2/6-31+G* and fitted V0 (kcal/mol), R0 (Å) values for orientation in Figure 2a−2h are −0.525, 4.29 and −0.415,4.46; −0.647, 3.90 and −0.779, 3.71; −0.567, 4.44 and −0.597, 4.44; −0.515, 4.81 and −0.476, 4.96; −0.193, 7.39 and −0.492, 7.20; −0.295, 7.47 and −0.166, 7.47; −0.516, 6.16 and −0.449, 6.16; and −0.755, 4.73 and −0.958, 4.73. Spin-component-scaled (SCS) MP2 corrects some of the shortcomings of MP2 and often gives more accurate C

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A intermolecular potential energy curves.52−54 To compare with the above MP2 calculations, SCS-MP2/6-311++G** potential energy curves were calculated for all eight Az−N2 orientations and are given in Figure 3, where they are compared with the MP2/6-31+G* curves. Overall, the SCS-MP2/6-311++G** curves are less attractive. The solid black line in Figure 3 is the fit to the SCS-MP2/6-311++G** curves with the two-body potential in eq 2 as described above. The parameters for the fit are listed in Table 1. The respective SCS-MP2/6-311++G** and fitted V0 (kcal/mol), R0 (Å) values for orientation in Figure 3a−3h are −0.314, 4.46 and −0.227,4.46; −0.367, 4.09 and −0.467, 3.90; −0.339, 4.60 and −0.376, 4.60; −0.314, 4.96 and −0.297, 4.96; −0.245, 8.58 and −0.636, 8.39; −0.385, 8.53 and −0.106, 8.93; −0.318, 4.70 and −0.312, 4.70; and −0.443, 5.92 and −0.664, 5.92.

of atomic positions and velocities are saved at random intervals for the N2 bath. Each of these configurations was selected as the solvent’s initial conditions for 12 trajectories, for which azulene had random initial conditions as specified by the quantum microcanonical sampling. Each of these 8 × 12 = 96 trajectories was then further integrated by keeping azulene fixed and equilibrating the N2 bath for 12 ps. These 96 trajectories, with random initial conditions for both azulene and the N2 bath, comprised an initial trajectory ensemble. For the simulations reported here, ensembles with 48 trajectories were employed. As shown previously,31 an ensemble of 24 trajectories provides semiquantitative results. All of the calculations, including the equilibrations and simulations of energy transfer, were performed with the nearest-neighbor list algorithm.57 This algorithm was used to save calculation time with negligible loss in accuracy. A “cut-off” distance of 18 Å was used for the nearest neighbor list. For the C6F6 + N2 simulations, the single collision limit was attained at a bath density of 40 kg/m3, that is, 35 atm.31,32 However, Az is larger than C6F6, and the Az mean free path within the N2 bath is roughly one-half of that for C6F6. Thus, for the Az + N2 bath system it is expected that a bath density lower than 40 kg/m3 is required to obtain the single collision limit. The Az mean free path is 5.05 Å at a 40 kg/m3 N2 bath density and 10.1 Å at a 20 kg/m3 bath density. For comparison, the experimental Az−N2 collision diameter is 6.07 Å.20 The chemical dynamics simulations were performed for the four bath densities of 10, 20, 40, and 80 kg/m3. The corresponding sides of the “box” for these simulations are 167.2, 132.7, 105.4, and 83.7 Å. For all the densities the trajectories were integrated up to 3.3 ns, with 2 fs of integration step size, using a sixth-order symplectic numerical integration algorithm.69 As described above, an ensemble of 48 trajectories was calculated for each simulation. The majority of the simulations were performed with the harmonic Az intramolecular potential and the Az−N 2 intermolecular potential fit to the MP2/6-31+G* calculations. One simulation was performed with the anharmonic Az intramolecular potential, and another was performed with the SCS-MP2/6-311++G** intermolecular potential for Az−N2, both for the bath density of 20 kg/m3, to compare with the results of the harmonic intramolecular and MP2/6-31+G* intermolecular potential. Overall, the simulations results were found to be insensitive to the intramolecular and intermolecular potentials. The results are discussed in the following section.

III. SIMULATION METHODOLOGY The simulations reported here were performed using a condensed-phase version of the VENUS chemical dynamics computer program.55,56 Azulene (Az) collisional relaxation in a N2 bath was simulated by initially placing randomly excited Az in the center of a cubical box with the surrounding N2 bath thermally equilibrated. Periodic boundary conditions57 were used for this simulation box. Details of the simulation methodology have been described previously31,32 and are briefly outlined in the following. Experimentally Az was excited to the S2 electronic state by 337 nm laser excitation.15−21 This electronic state then undergoes an internal conversion (IC) to the vibrationally excited ground electronic state (S0*), which contains 87.5 kcal/ mol of vibrational energy. To simulate this experiment, a quasiclassical microcanonical sampling algorithm58,59 was used to add this energy to the Az vibrational energy levels, as done in previous work.32 The Az zero-point energy (ZPE) of 95.7 kcal/ mol is explicitly added by this algorithm. The initial translational and rotational energies of azulene were sampled from their Boltzmann distributions at 298 K. To obtain agreement with experiment for classical dynamics simulations, previous work has illustrated the importance of using the quasiclassical method for which ZPE is included in the initial conditions.58−63 For the simulations reported here, previous chemical dynamics simulations of intramolecular vibrational energy redistribution (IVR) are the most demonstrative. IVR rates are obtained in good agreement with experiment and quantum dynamics simulations if ZPE is included in the initial conditions.64−67 Without the addition of ZPE the IVR rate is much too small.68 As described above, for the simulations reported here Az ZPE is included in the initial conditions. For comparison, a simulation is also performed without Az ZPE. The N2 bath was equilibrated to a temperature of 298 K by a molecular dynamics (MD) simulation, as described previously.31,32 This equilibration is performed with the azulene molecule in the middle of the simulation cubical box, with the randomly chosen coordinates and velocities of azulene kept fixed. The MD equilibration was performed using a maximum of seven stages, with each stage consisting of velocity rescaling and monitoring components of the N2 bath, for a maximum total equilibration time of 172 ps. The equilibration was complete when an accurate temperature, radial distribution function, and vibrational, rotational, and center-of-mass translational energies were obtained for the bath. After equilibration, the trajectory was further integrated for 12 ps, where eight sets

IV. RESULTS AND DISCUSSION The simulation results presented in the following were determined using the Az harmonic intramolecular potential and the Az−N2 MP2/6-31+G* intermolecular potential unless stated otherwise. One simulation was performed with the Az anharmonic intramolecular potential, and another was performed with the SCS-MP2/6-311++G** intermolecular potential for Az−N2. A. Average Azulene Energy Versus Time. As shown in Figure 4, for all bath densities the average vibrational energy of Az versus time, with the harmonic intramolecular potential and averaged over 48 trajectories, is well fit by the biexponential expression ⟨E(t )⟩ = (E(0) − E(∞)) × [f1 e−k1t + f2 e−k 2t ] + E(∞) (3) D

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 4. Average vibrational energy of Az, using the Az harmonic intramolecular potential and the MP2/6-31+G* Az−N2 intermolecular potential, vs time for N2 bath densities of 10, 20, 40, and 80 kg/ m3, averaged over 48 trajectories. There are 1000 N2 molecules in the bath. Average energy of Az vs time × N2 bath density (i.e., t × ρ) is plotted in the inset.

Figure 5. Comparison of the average vibrational energy of Az vs time for N2 bath density of 20 kg/m3 for harmonic and anharmonic intramolecular potential as discussed in Section II.A, averaged over 48 trajectories. The other trajectory conditions are the same as those in Figure 4.

E(0) is the initial energy of excited Az, E(∞) is the Az final energy, f1 + f 2 = 1, and k1 and k2 are rate constants. These fitting parameters are presented in Table 2. In Figure 5 the average energy of Az versus time is compared for simulations with harmonic and anharmonic Az intramolecular potentials, for a N2 bath density of 20 kg/m3. As seen from the figure changing the intramolecular potential has only a slight effect on the energy transfer rate. The curve for the anharmonic potential is also fit by eq 3, and the parameters are listed in Table 2. B. Single Collision Limit. As discussed in previous work,31,32 collision energy transfer dynamics become independent of the bath density (i.e., pressure) in the single collision limit. In this limit, the relaxation rate constants k1 and k2 are proportional to the bath density ρ, with proportionality constants C1 and C2 defined by k1 = ρ × C1, k2 = ρ × C2. The criterion for the single collision limit is that C1 and C2 become independent of ρ. As shown in Table 2 this behavior is found at ρ = 20 kg/m3, where C1 is the same within three significant figures as that for the lower ρ of 10 kg/m3, and the C2 of 0.000 0207 is nearly the same as the 0.000 0208 value at 10 kg/ m3. Similar identical C1 and nearly identical C2, at different ρ, were used to identify the single collision limit for the C6F6 + N2 simulations.32 The approach to the single collision limit is also illustrated by plots of ⟨E(t)⟩ versus t × ρ for different bath densities.32 Such curves are presented for all the densities in the Figure 4 inset.

For the previous C6F6 + N2 simulations,31,32 the single collision limit is found at ρ = 40 kg/m3, where the mean free path is 12.5 Å. In contrast, for the Az + N2 simulations the single collision limit is found for the lower ρ = 20 kg/m3, with a similar mean free path of 10.1 Å. The Az + N2 mean free path is 5.05 Å for ρ = 40 kg/m3 and too small to obtain the single collision limit. The Az−N2 collision diameter is 6.07 Å.20 The C1 and C2 values for Az + N2 at 40 kg/m3 are similar to those at 80 kg/m3, but not identical as those above for ρ of 10 and 20 kg/m3, where the energy transfer dynamics are in the single collision limit. For the C6F6 + N2 simulation,31 similar C1 values were found at ρ of 201 and 324 kg/m3, where the dynamics are not in the single collision limit. The dependence of C1 and C2 on ρ in the nonsingle collision limit is uncertain, but the previous simulations for C6F6 + N2 indicate that they decrease with increase in ρ in this limit, which is found here for ρ of 40 and 80 kg/m3. In future work it would be of interest to consider nonbinary collision models for interpreting C1 and C2 values at high bath densities. C. Average Energy Transfer per Collision. In the single collision limit, experimental and simulation values of the average energy transfer per collision, ⟨ΔEc⟩, may be compared. In this limit ⟨ΔEc⟩ is given by ⟨ΔEc⟩ =

d⟨E⟩ 1 × dt ω

(4)

Table 2. Parameters for Fits to ⟨E(t)⟩ for Different N2 Bath Densitiesa ρ (kg/m3)b

E(∞) (kcal/mol)

f1

f2

k1

k2

C1

C2

10 20 20 (anhar) 40 80

31.0 33.0 33.2 33.0 33.4

0.0120 0.0130 0.0137 0.0140 0.0170

0.988 0.987 0.986 0.986 0.983

0.001 18 0.002 35 0.002 52 0.003 00 0.005 90

0.000 208 0.000 414 0.000 433 0.000 761 0.001 45

0.000 118 0.000 118 0.000 126 0.000 0750 0.000 0738

0.000 0208 0.000 0207 0.000 0217 0.000 0190 0.000 0181

a The fits are to eq 3. The sum of f1 and f 2 is set to unity in the fitting, and the k values are in unit of inverse picoseconds. C1 and C2 are calculated by k1/ρ and k2/ρ, respectively. All of the simulations are performed with the harmonic azulene intramolecular potential, except for the one simulation at 20 kg/m3, where the anhramonic intramolecular potential was used. bThe pressures for the lowest and highest densities of 10 and 80 kg/m3 are 8.75 and 70.0 atm, respectively.

E

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

energies were monitored during the simulation. The Az translational and rotational energies versus time are shown in Figure 7 upper and lower panels, respectively, and it is seen that

where ω is collision frequency and d E is energy transfer per dt unit time. The experimental collision frequency has been given previously31 and is obtained from ⎛ T ⎞−1/2 * ω = 4.415 × 107 × σA2− N p⎜ ⎟ × Ω12 ⎝ μ⎠

(5)

with −1 ⎡ ⎛ ⎞⎤ * = ⎢0.636 + 0.567 × log⎜ kT ⎟⎥ Ω12 ⎢⎣ ⎝ ϵ12 ⎠⎥⎦

(6)

Here σA‑N = (σA + σN)/2 and ϵ12 (= ϵ1ϵ2 ) are the LennardJones parameters for azulene + N2, p is the pressure in torr, the temperature T is in Kelvin, μ is the reduced mass in amu, and Ω is the collision integral, with σA = 6.61 Å, ϵA = 523 K,20 σN = 3.74 Å, and ϵ2 = 82 K.7 From the ⟨E(t)⟩ in Figure 3, with parameters in Table 2, the energy transfer per collision (ΔEc) is obtained. These data are plotted in Figure 6 for the bath densities ρ of 10, 20, and 40 kg/

Figure 7. Center of mass translational (top) and rotational (bottom) energies of azulene molecule vs simulation time, averaged over 48 trajectories.

Figure 6. Average energy transferred per collision vs average energy of azulene for the N2 bath densities of 10, 20, and 40 kg/m3. (inset) ⟨ΔEc⟩ vs [⟨E⟩ − E(∞)] [see eq 3]. The trajectory conditions are the same as those in Figure 4.

there are no significant V → T or V → R energy transfers for Az. On the one hand, for all of the Az−N2 simulation densities, the Az rotational and translational energies are never over 3.6 and 2.1 kcal/mol, respectively. On the other hand, as shown in Figure 8, the energy transferred from Az is used to increase the translational (top) and rotational (middle) energies of the N2 bath. There is ∼50% more energy transferred to translation than rotation, consistent with equipartitioning of energy to these degree of freedoms. Because of the high N2 vibrational frequency, energy transfer to this degree of freedom is negligible, and the vibrational energy of the N2 bath fluctuates near its 298 K equilibrium value, which is ∼600 kcal/mol for the 1000 molecules (bottom). The energy transfer pathways for the Az−N2 system are the same as those found previously32 for the C6F6−N2 system. E. Comparison of Classical and Quasiclassical Simulations with and without Azulene Zero-Point Energy. A chemical dynamics simulation was also performed for N2 + Az using classical microcanonical normal mode sampling for Az, for which only the excitation energy of 87.5 kcal/mol is added to Az without its ZPE. The calculation is done at the single collision limiting density of 20 kg/m3 of the N2 bath. A ⟨ΔEc⟩

m3. For ρ of 10 and 20 kg/m3, the ⟨ΔEc⟩ versus ⟨E ⟩ plots are almost identical. However, there is a significant difference between the plots for the 20 and 40 kg/m3 bath densities. This observation, along with the C1 and C2 parameters for the different densities in Table 2, shows that the single collision limit is not attained at ρ of 40 kg/m3. As discussed previously,32 slopes of plots of ⟨ΔEc⟩ versus [⟨E⟩ − E(∞)] become identical at or below ρ for the single collision limit, as depicted in the Figure 6 inset for the current simulations. ⟨ΔEc⟩ for the harmonic azulene intramolecular potential is 0.30 ± 0.01 kcal/mol for ⟨E ⟩ of 183.2 kcal/mol, that is, at the excitation energy of 87.5 kcal/mol with ZPE. This value is somewhat smaller than the experimental value of 0.57 ± 0.11 kcal/mol.20 ⟨ΔEc⟩ at this energy for the anharmonic azulene intramolecular potential is only slight larger and 0.32 ± 0.01 kcal/mol. D. Energy Transfer Pathways. The excitation energy of vibrationally excited Az can be transferred to translational and rotational energy of Az, as well as translational, rotational, and vibrational energy of the N2 bath. To investigate this, these F

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the harmonic intramolecular potential for Az, but with the Az− N2 intermolecular potential fit to the higher level SCS-MP2/6311++G** calculations. As shown in the plot of ⟨E(t)⟩ versus t in the Supporting Information (Figure S1), the IET dynamics with this intermolecular potential are statistically the same as that for the MP2/6-31+G* intermolecular potential, used for the simulation results presented above. This insensitivity of the IET to fine details of the intermolecular potential is the same result found for the C6F6 + N2 IET.31

V. SUMMARY In this article, gas-phase energy transfer from highly vibrationally excited azulene (Az) to a N2 bath is studied with chemical dynamics simulations. Intermolecular potentials between Az and N2 were determined by performing single-point electronic structure calculations for eight Az−N2 orientations, with both MP2/6-31+G* and SCS-MP2/6-311++G** theory, and then simultaneously fitting the points by a sum of C−N and H−N two-body potentials as given by eq 2. The same N2 potential and N2−N2 and intermolecular potential were used as described and used previously.31 A harmonic Az intramolecular potential was obtained from previous work done by Bernshtein and Oref.41 An anharmonic Az intramolecular potential was developed by including anharmonicity in all the stretch and bend potential energy terms. Chemical dynamics simulations, for collisional deactivation of highly vibrationally excited Az by the N2 bath, that is, Az* + N2, were initiated with a total Az* vibrational energy of 183.2 kcal/mol (i.e., 87.5 kcal/mol excitation energy and a 95.7 kcal/ mol ZPE) to compare with experiment.15−21 Simulations for the harmonic Az potential and MP2/6-31+G* Az−N 2 intermolecular potential were performed for a N2 bath temperature of 298 K and bath densities of 10, 20, 40, and 80 kg/m3; that is, pressures ranging from 9 to 70 atm were considered. The average energy of Az* versus time was fit by the biexponential function in eq 3, and the best fit values for the f1, f 2, k1, and k2 parameters are listed in Table 2. Values for C1 = k1/ρ and C2 = k2/ρ, where ρ is the bath density, become independent of pressure (or density) in the single collision limit. These parameters are given in Table 2 for the different densities, and it is seen that the single collision limit is attained at the bath density of 20 kg/m3. This is a lower density than required for the previous C6F6 + N2 simulations,31,32 as a result of the larger collision diameter for Az as compared to C6F6. For the 20 kg/m3 bath density, simulations were compared for the harmonic and anharmonic azulene intramolecular potential. As shown in Figure 5, the anharmonic potential has only a small effect on the energy transfer dynamics. Simulations at the 20 kg/m3 bath density were also performed with purely classical initial conditions for Az, without ZPE, and using the SCS-MP2/ 6-311++G** Az−N2 intermolecular potential, and they give statistically the same results as found with quasiclassical sampling for Az and the MP2/6-31+G* Az−N2 intermolecular potential. The average energy transferred from Az* per collision, ⟨ΔEc⟩, was determined from experiment at an excitation energy of 87.5 kcal/mol and is ⟨ΔEc⟩ = 0.57 ± 0.11 kcal/mol.20 The value determined from the simulation is slightly lower, but still in good agreement with experiment; that is, ⟨ΔEc⟩ is 0.30 ± 0.01 and 0.32 ± 0.01 kcal/mol for the harmonic and anharmonic Az potentials, respectively.

Figure 8. Center of mass translational (top), rotational (middle), and vibrational (bottom) energy of 1000 N2 bath molecules vs simulation time, averaged over 48 trajectories.

value of 0.27 ± 0.02 kcal/mol is obtained from this calculation for the 87.5 kcal/mol excitation energy. The simulations were performed with the Az harmonic intramolecular potential, and the Az−N2 intermolecular potential was determined from the MP2/6-31+G* calculations. Using these potentials and a quasiclassical simulation, which includes ZPE as described above, a similar value of 0.30 ± 0.01 kcal/mol was obtained for ⟨ΔEc⟩. The experimental value is 0.57 ± 0.11 kcal/mol. F. Comparison of the Energy Transfer Dynamics for MP2/6-31+G* and SCS-MP2/6-311++G** Az−N2 Intermolecular Potentials. The above simulations were performed with an Az−N2 intermolecular potential fit to MP2/6-31+G* calculations. To determine if the IET dynamics depend on the nature of the Az−N2 intermolecular potential, simulations were performed at ρ = 20 kg/m3 using quasiclassical sampling and G

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

■

Vibrational energy transfer from excited Az is primarily to the translational and rotational degrees of freedom of the N2 bath, with negligible energy transfer from Az vibration to either Az rotation or translation. The energy transfer to translation of the bath is ∼50% larger than to bath rotation, which is consistent with equipartition energy transfer to bath translation and rotation. These dynamics for energy transfer from vibrationally excited Az to the N2 bath are the same as those found previously for the C6F6−N2 system.31,32 The “rate constants” k1 and k2 in eq 3 are proportional to the bath density; that is, ki = Ciρ. As discussed above, in the lowdensity, single collision limit the Ci values become constant and independent of ρ. Also of interest is to consider the Ci at high densities, ρ, and compare the Ci found here for azulene−N2 with those reported previously for C6F6−N2. The C1 and C2 for azulene−N2 and C6F6−N2 are compared in Figure 9, upper and

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b00893. The comparison of the averaged vibrational energy of Az versus time for N2 bath density of 20 kg/m3 between MP2/6-31+G* and SCS-MP2/6-311++G** intermolecular potential. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The research reported here is based upon work supported by the Air Force Office of Scientific Research, BRI Grant No. FA 9550-12-1-0443, and the Robert A. Welch Foundation under Grant No. D-0005. Support was also provided by the High Performance Computing Center at Texas Tech Univ, under the direction of P. W. Smith. Parts of the computations were also performed on Robinson, a general computer cluster of the Dept. of Chemistry and Biochemistry, Texas Tech Univ, purchased by the NSF CRIF-MU Grant No. CHE-0840493.

■

REFERENCES

(1) Tardy, D. C.; Rabinovitch, B. S. Intermolecular Vibrational Energy Transfer in Thermal Unimolecular Systems. Chem. Rev. 1977, 77, 369−408. (2) Quack, M.; Troe, J. Gas Kinetics and Energy Transfer; The Chemical Society: London, U.K., 1977, Vol: 2. (3) Oref, I.; Tardy, D. C. Energy Transfer in Highly Excited Large Polyatomic Molecules. Chem. Rev. 1990, 90, 1407−1445. (4) Dove, J. E.; Hippler, H.; Troe, J. Direct Study of Energy Transfer of Vibrationally Highly Excited CS2 Molecule. J. Chem. Phys. 1985, 82, 1907−1919. (5) Heymann, M.; Hippler, H.; Nahr, D.; Plach, H. J.; Troe, J. UV Absorption Study of Collisional Energy Transfer in Vibrationally Highly Excited Sulfur Dioxide Molecules. J. Phys. Chem. 1988, 92, 5507−5514. (6) Hippler, H.; Troe, J.; Wendelken, H. J. Collisional Deactivation of Vibrationally Highly Excited Polyatomic Molecules. III. Direct Observations for Substituted Cycloheptatrienes. J. Chem. Phys. 1983, 78, 6718−6724. (7) Hippler, H.; Troe, J.; Wendelken, H. J. UV Absorption Spectra of Vibrationally Highly Excited Toluene Molecules. J. Chem. Phys. 1983, 78, 5351−5357. (8) Miller, L. A.; Barker, J. R. Collisional Deactivation of Highly Vibrationally Excited Pyrazine. J. Chem. Phys. 1996, 105, 1383−1391. (9) Miller, L. A.; Cook, C. D.; Barker, J. R. Temperature Effects in the Collisional Deactivation of Highly Vibrationally Excited Pyrazine by Unexcited Pyrazine. J. Chem. Phys. 1996, 105, 3012−3018. (10) Hartland, G. V.; Qin, D.; Dai, H. L. Collisional Deactivation of Highly Vibrationally Excited NO2 Monitored by Time-Resolved Fourier Transform Infrared Emission Spectroscopy. J. Chem. Phys. 1994, 100, 7832−7835. (11) Toselli, B. M.; Barker, J. R. Quantum Effects in Large Molecule Collisional Energy Transfer? Chem. Phys. Lett. 1990, 174, 304−308. (12) Yerram, M. L.; Brenner, J. D.; King, K. D.; Barker, J. R. Collisional Deactivation of Highly Vibrationally Excited Benzene Pumped at 248 nm. J. Phys. Chem. 1990, 94, 6341−6350. (13) Toselli, B. M.; Barker, J. R. Isotope Effects in the Vibrational Deactivation of Large Molecules. J. Chem. Phys. 1992, 97, 1809−1817.

Figure 9. C1 (= k1/ρ) and C2 (= k2/ρ) values for C6F6 + N2 (see refs 31 and 32), and azulene + N2 simulations for all N2 bath densities of respective calculations. The y-axis is presented in logarithmic scale.

lower panels, respectively, where it is seen that that the azulene−N2 values are approximately an order of magnitude lower. The Ci decrease, from their values for the single collision limit, as the bath density is increased. Important questions to address in future research is the origin (or origins) of this decrease and if the Ci reach decreased limiting values as the bath density is further increased, for example, to pressures higher than the critical point pressure of 33.5 atm. Additional future work includes probing why the simulation ⟨ΔEc⟩ is somewhat lower than the experimental value and to study collisional energy transfer for additional gas and liquid phase systems. H

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (14) Brenner, J. D.; Erinjeri, J. P.; Barker, J. R. Population Distributions in the Vibrational Deactivation of Benzene and Benzene-d6. First and Second Moments Derived from Two-Color Infrared Fluorescence Measurements. Chem. Phys. 1993, 175, 99−111. (15) Rossi, M. J.; Barker, J. R. Infrared Fluorescence and Collisional Energy Transfer Parameters for Vibrationally Excited Azulene*(So): Dependence on Internal Energy (Evib). Chem. Phys. Lett. 1982, 85, 21− 26. (16) Barker, J. R. Direct Measurements of Energy-Transfer Involving Large Molecules in the Electronic Ground State. J. Phys. Chem. 1984, 88, 11−18. (17) Rossi, M. J.; Pladziewicz, J. R.; Barker, J. R. Energy-Dependent Energy Transfer: Deactivation of Azulene (S0, Evib) by 17 Collider Gases. J. Chem. Phys. 1983, 78, 6695−6708. (18) Shi, J.; Bernfeld, D.; Barker, J. R. Energy Dependence of Infrared Emission from Azulene C-H Stretching Vibrations. J. Chem. Phys. 1988, 88, 6211−6218. (19) Shi, J.; Barker, J. R. Energy-Dependent Collisional Deactivation of Vibrationally Excited Azulene. J. Chem. Phys. 1988, 88, 6219−6227. (20) Hippler, H.; Lindemann, L.; Troe, J. Collisional Energy Transfer of Vibrationally Highly Excited Molecules. V. UV Absorption Study of Azulene. J. Chem. Phys. 1985, 83, 3906−3912. (21) Jalenak, W.; Weston, R. E., Jr.; Sears, T. J.; Flynn, G. W. Energy Transfer from Highly Vibrationally Excited Azulene and Azulene-d8 to Carbon Dioxide. J. Chem. Phys. 1988, 89, 2015−2022. (22) Bernshtein, V.; Oref, I. Intermolecular Energy Transfer Probabilities from Trajectory Calculations: A New Approach. J. Chem. Phys. 1998, 108, 3543−3553. (23) Bernshtein, V.; Oref, I. Termolecular Collisions between Benzene and Ar. J. Chem. Phys. 2003, 118, 10611−10622. (24) Lenzer, T.; Luther, K.; Troe, J.; Gilbert, R. G.; Lim, K. F. Trajectory Simulations of Collisional Energy Transfer in Highly Excited Benzene and Hexafluorobenzene. J. Chem. Phys. 1995, 103, 626−641. (25) Lenzer, T.; Luther, K. Intermolecular Potential Effects in Trajectory Calculations of Collisions between Large Highly Excited Molecules and Noble Gases. J. Chem. Phys. 1996, 105, 10944−10953. (26) Bernshtein, V.; Oref, I. Collisional Energy Transfer between Ar and Normal and Vibrationally and Rotationally Frozen Internally Excited Benzene-Trajectory Calculations. J. Chem. Phys. 1997, 106, 7080−7089. (27) Kable, S. H.; Knight, A. E. W. Semiempirical Model of Vibrational Relaxation for Estimating Absolute Rate Coefficients. J. Phys. Chem. A 2003, 107, 10813−10825. (28) Lim, K. F. Quasiclassical Trajectory Study of Collisional Energy Transfer in Toluene Systems. I. Argon Bath Gas: Energy Dependence and Isotope Effects. J. Chem. Phys. 1994, 100, 7385−7399. (29) Bernshtein, V.; Oref, I. Trajectory Calculations of Relative Center of Mass Velocities in Collisions between Ar and Toluene. J. Chem. Phys. 1996, 104, 1958−1965. (30) Gilbert, R. G. Mechanism and Models for Collisional Energy Transfer in Highly Excited Large Polyatomic Molecules. Aust. J. Chem. 1995, 48, 1787−1817. (31) Paul, A. K.; Kohale, S. C.; Pratihar, S.; Sun, R.; North, S. W.; Hase, W. L. A Unified Model for Simulating Liquid and Gas Phase, Intermolecular Energy Transfer: N2 + C6F6 Collisions. J. Chem. Phys. 2014, 140, 194103. (32) Paul, A. K.; Kohale, S. C.; Hase, W. L. Bath Model for N2 + C6F6 Gas-Phase Collisions. Details of the Intermolecular Energy Transfer Dynamics. J. Phys. Chem. C 2015, 119, 14683−14691. (33) Lim, K. F.; Gilbert, R. G. The Apriori Calculation of Collisional Energy Transfer in Highly Vibrationally Excited Molecules: The Biased Random Walk Model. J. Chem. Phys. 1986, 84, 6129−6140. (34) Lim, K. F.; Gilbert, R. G. Modeling Collisional Energy Transfer in Highly Excited Molecules. J. Chem. Phys. 1990, 92, 1819−1830. (35) Clarke, D. L.; Oref, I.; Gilbert, R. G.; Lim, K. F. Collisional Energy Transfer in Highly Excited Molecules: Calculations of the Dependence on Temperature and Internal, Rotational, and Translational energy. J. Chem. Phys. 1992, 96, 5983−5998.

(36) Clarke, D. L.; Gilbert, R. G. Collisional Energy Transfer in Highly Excited Molecules: Deuteration Effects. J. Phys. Chem. 1992, 96, 8450−8453. (37) Heidelbach, C.; Fedchenia, I. I.; Schwarzer, D.; Schroeder, J. Molecular-Dynamics Simulation of Collisional Energy Transfer from Vibrationally Highly Excited Azulene in Compressed CO2. J. Chem. Phys. 1998, 108, 10152−10161. (38) Heidelbach, C.; Vikhrenko, V. S.; Schwarzer, D.; Schroeder, J. Molecular Dynamics Simulation of Vibrational Relaxation of Highly Excited Molecules in Fluids. II. Nonequilibrium Simulation of Azulene in CO2 and Xe. J. Chem. Phys. 1999, 110, 5286−5299. (39) Lim, K. F.; Gilbert, R. G. Trajectory Simulations of Collisional Energy Transfer of Highly Vibrationally Excited Azulene. J. Phys. Chem. 1990, 94, 77−84. (40) Bernshtein, V.; Oref, I. Energy Transfer between Azulene and Krypton: Comparison between Experiment and Computation. J. Chem. Phys. 2006, 125, 133105. (41) Barker, J. R.; Golden, R. E. Temperature-Dependent Energy Transfer: Direct Experiments Using Azulene. J. Phys. Chem. 1984, 88, 1012−1017. (42) Damm, M.; Deckert, F.; Hippler, H.; Troe, J. Isomerization and Collisional Deactivation of Highly Vibrationally Excited Azulene Molecules after UV Excitation at 248 and 193 nm. J. Phys. Chem. 1991, 95, 2005−2009. (43) Nilsson, C.; Nordholm, S. Modeling Energy Transfer in Molecular Collisions: Statistical Theory versus Experiment for Highly Excited Toluene and Azulene. J. Chem. Phys. 2003, 119, 11212−11220. (44) Rivera-Rivera, L. A.; Wagner, A. F.; Sewell, T. D.; Thompson, D. L. Pressure Effects on the Relaxation of an Excited Nitromethane Molecule in an Argon Bath. J. Chem. Phys. 2015, 142, 014303. (45) Johnston, H. S. Gas Phase Reaction Rate Theory; The Ronald Press Company: New York, 1966. (46) Kolakkandy, S.; Paul, A. K.; Pratihar, S.; Kohale, S. C.; Barnes, G. L.; Wang, H.; Hase, W. L. Energy and Temperature Dependent Dissociation of the Na+(Benzene)1,2 Clusters: Importance of Anharmonicity. J. Chem. Phys. 2015, 142, 044306. (47) Lu, D. H.; Hase, W. L. Classical Trajectory Calculation of the Benzene Overtone Spectra. J. Phys. Chem. 1988, 92, 3217−3225. (48) Sinnokrot, M. O.; Sherrill, C. D. Highly Accurate Coupled Cluster Potential Energy Curves for the Benzene Dimer: Sandwich, TShaped, and Parallel-Displaced Configurations. J. Phys. Chem. A 2004, 108, 10200−10207. (49) Vayner, G.; Alexeev, Y.; Wang, J.; Windus, T.; Hase, W. L. Ab Initio and Analytic Intermolecular Potentials for Ar-CF4. J. Phys. Chem. A 2006, 110, 3174−3178. (50) Paul, A. K.; Kolakkandy, S.; Hase, W. L. Dynamics of Na+(Benzene) + Benzene Association and Ensuing Na+(Benzene)2* Dissociation. J. Phys. Chem. A 2015, 119, 7894−7904. (51) Marques, J. M. C.; Prudente, V. F.; Pereira, F. B.; Almeida, M. M.; Maniero, A. M.; Fellows, C. E. A New Genetic Algorithm to be used in the Direct Fit of Potential Energy Curves to Ab Initio and Spectroscopic Data. J. Phys. B: At., Mol. Opt. Phys. 2008, 41, 085103. (52) Grimme, S. Improved Second-Order Møller-Plesset Perturbation Theory by Separate Scaling of Parallel- and Antiparallel-Spin Pair Correlation Energies. J. Chem. Phys. 2003, 118, 9095−9102. (53) Grimme, S.; Diedrich, C.; Korth, M. The Importance of Interand Intramolecular van der Waals Interactions in Organic Reactions: the Dimerization of Anthracene Revisited. Angew. Chem., Int. Ed. 2006, 45, 625−629. (54) Grimme, S. Accurate Calculation of the Heats of Formation for Large Main Group Compounds with Spin-Component Scaled MP2Methods. J. Phys. Chem. A 2005, 109, 3067−3077. (55) Hase, W. L.; Duchovic, R. J.; Hu, X.; Komornicki, A.; Lim, K. F.; Lu, D. -h.; Peslherbe, G. H.; Swamy, K. N.; Vande Linde, S. R.; Varandas, A.; et al. VENUS96: A General Chemical Dynamics Computer Program; Texas Tech University: Lubbock, TX, 2005. (56) Hu, X.; Hase, W. L.; Pirraglia, T. Vectorization of the General Monte Carlo Classical Trajectory Program VENUS. J. Comput. Chem. 1991, 12, 1014−1024. I

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (57) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (58) Peslherbe, G. H.; Wang, H.; Hase, W. L. Monte Carlo Sampling for Classical Trajectory Simulations. Adv. Chem. Phys. 1999, 105, 171− 201. (59) Park, K.; Engelkemier, J.; Persico, M.; Manikandan, P.; Hase, W. L. Algorithms for Sampling a Quantum Microcanoical Ensemble of Harmonic Oscillators at Potential Minima and Conical Intersections. J. Phys. Chem. A 2011, 115, 6603−6609. (60) Swamy, K. N.; Hase, W. L. A Quasiclassical Trajectory Calculation of the H + C2H4 → C2H5 Bimolecular Rate Constant. J. Phys. Chem. 1983, 87, 4715−4720. (61) Wang, X.; Bowman, J. M. Zero-point Energy is Needed in Molecular Dynamics Calculations to Access the Saddle Point for H + HCN → H2CN* and cis/trans-HCNH* on a New Potential Energy Surface. J. Chem. Theory Comput. 2013, 9, 901−908. (62) Han, Y.-C.; Bowman, J. M. Reactant Zero-point Energy is Needed to Access the Saddle Point in Molecular Dynamics Calculations of the Association Reaction H + C2D2 → C2D2H. Chem. Phys. Lett. 2013, 556, 39−43. (63) Cotton, S. J.; Miller, W. H. The Symmetrical Quasi-Classical Model for Electronically Non-Adiabatic Processes Applied to Energy Transfer Dynamics in Site-Exciton Models of Light-Harvesting Complexes. J. Chem. Theory Comput. 2016, 12, 983−991. (64) Sibert, E. L., III; Hynes, J. T.; Reinhardt, W. P. Classical Dynamics of Highly Excited CH and CD Overtones in Benzene and Perdeuterobenzene. J. Chem. Phys. 1984, 81, 1135−1144. (65) Lu, D.-h.; Hase, W. L. Classical Trajectory Calculation of the Benzene Overtone Spectra. J. Phys. Chem. 1988, 92, 3217−3225. (66) Wyatt, R. E.; Iung, C.; Leforestier, C. Quantum Dynamics of Overtone Relaxation in Benzene. II. 16 Mode Model for Relaxation from CH(v = 3). J. Chem. Phys. 1992, 97, 3477−3486. (67) Stock, G. Classical Simulation of Quantum Energy Flow in Biomolecules. Phys. Rev. Lett. 2009, 102, 118301. (68) Lu, D.-h.; Hase, W. L. Classical Mechanics of Intramolecular Vibrational Energy Flow in Benzene. V. Effect of Zero Point Energy Motion. J. Chem. Phys. 1989, 91, 7490−7497. (69) Schlier, C.; Seiter, A. High-Order Symplectic Integration: An Assessment. Comput. Phys. Commun. 2000, 130, 176−189.

J

DOI: 10.1021/acs.jpca.6b00893 J. Phys. Chem. A XXXX, XXX, XXX−XXX