390
J. Phys. Chem. B 2008, 112, 390-398
Vibrational Energy Relaxation of the Bend Fundamental of Dilute Water in Liquid Chloroform and d-Chloroform† Y.-S. Lin, S. G. Ramesh, J. M. Shorb, E. L. Sibert III, and J. L. Skinner* Theoretical Chemistry Institute and Department of Chemistry, UniVersity of Wisconsin, Madison, Wisconsin 53706 ReceiVed: July 19, 2007; In Final Form: September 17, 2007
The population lifetimes of the bend fundamental of dilute water in liquid chloroform (8.5 ps) and d-chloroform (28.5 ps) display an interesting solvent isotope effect. As the lowest excited vibrational state of the molecule, the water bend fundamental relaxes directly to the ground state with about 1600 cm-1 of energy released to the other degrees of freedom. The strong solvent isotope effect along with the large energy gap indicates the participation of solvent vibrational modes in this vibrational energy relaxation process. We calculate the vibrational energy relaxation rates of the water bend in chloroform and d-chloroform using the LandauTeller formula with a new potential model developed and parametrized self-consistently to describe the chloroform-water interaction. The computed values are in reasonable agreement with the experimental results, and the trend for the isotope effect is correct. It is found that energy transfer to the solvent vibrations does indeed play an important role. Nevertheless, no single dominant solvent accepting mode can be identified; the relaxation appears to involve both the bend and the C-Cl stretches, and frequency changes of all of these modes upon deuteration contribute to the observed solvent isotope effect.
I. Introduction An ensemble of molecules in thermal equilibrium has a Boltzmann distribution of its vibrational eigenstate populations. This distribution can be perturbed during the course of a chemical reaction or by interaction with external fields. The evolution of the resulting nonequilibrium vibrational energy distribution toward equilibrium is the problem of vibrational energy relaxation (VER).1-7 Understanding VER is important for understanding the rates, mechanisms, and pathways of chemical reactions. In particular, assessing the feasibility of mode-selective chemistry8-11 relies to a large extent on VER ratessif they are too large, then the energy deposited (for example, by a laser) will not remain in the selected modes long enough for reaction to occur. As such, VER has long been the object of intense theoretical and experimental study. The study of VER is particularly illuminating in small solute molecules with well-separated vibrational levels. In this case individual eigenstates can be excited, and the evolution of population to other individual eigenstates can be monitored. In a VER process, energy is typically transferred from solute vibrations to solvent modes, and so the question naturally arises as to the nature of the accepting solvent modes. In the simplest situation of a dilute diatomic solute in an atomic solvent (for example, dilute oxygen in liquid argon12,13), with the exception of the possible excitation of solute rotation, the accepting modes are translational. Next most complicated would be a mixture of a dilute diatomic in a solvent of a different diatomic molecule. If the vibrational frequency of the solvent was substantially higher than that of the solute, so that vibrational excitation of the solvent would be very unlikely, then the accepting modes for the solute VER process would include rotations and translations. More generally, for a polyatomic solvent, transla†
Part of the “James T. (Casey) Hynes Festschrift”. * Author to whom correspondence should be addressed. E-mail:
[email protected].
tions, rotations, and vibrations can all act as accepting modes, and it is of great interest to understand circumstances when one relaxation “channel” dominates over the others. A nice example of a relatively simple VER system involves excitation of the asymmetric stretch of azide anion in H2O and in D2O.14-17 It appears18 that the first step in the VER process is population transfer to the symmetric stretch of azide, with the concomitant deposition of some 700 cm-1 of energy into the solvent. From the observed solvent isotope effect (the lifetime is twice as long in D2O as it is in H2O)15-17 and from theoretical calculations,18 it seems clear that the primary relaxation channel involves excitation of the hindered rotations (librations) of nearby water molecules. Another interesting example involves the VER of the cyanide fundamental in D2O.19,20 In these experiments different isotopes of the cyanide solute are studied, with corresponding fundamental frequencies varying from 2004 to 2079 cm-1. Because cyanide is a diatomic, the only available relaxation pathway is directly to the ground state. If the accepting modes were low-frequency rotations and translations, then energy-gap-law considerations would indicate that the VER rate would decrease as the fundamental frequency increases. In fact the opposite is observed experimentally, prompting Hamm et al. to conclude19 that D2O stretch modes, centered at about 2400 cm-1 but with a low-frequency tail extending to below 2000 cm-1, act as energy acceptors. A final illustrative example is VER of iodine in the series of solvents xenon, carbon tetrachloride, chloroform, and methylene chloride.21 Experimentally, the VER rates increase as one progresses through this series in the order listed, as intermolecular vibration-vibration (V-V) energy transfer channels become operative and the characteristic frequencies of solvent vibrational, translational, and rotational modes change. VER has also been studied in liquid water22-31 and in isotopic mixtures of HOD in H2O and D2O.32-34 For the problem of
10.1021/jp075682s CCC: $40.75 © 2008 American Chemical Society Published on Web 11/29/2007
VER of Dilute H2O in Chloroform
J. Phys. Chem. B, Vol. 112, No. 2, 2008 391
TABLE 1: Fundamental Frequencies of Chloroform and d-Chloroform from Experiment46 and from the Sibert-Rey Model39 a frequency (cm-1) symmetry mode A1 A1 A1 E E E a
ν1 ν2 ν3 ν4 ν5 ν6
expt (H)
model (H)
exp (D)
model (D)
C-H(D) stretch 3019 C-Cl stretch 668 Cl-C-Cl bend 366 Cl-C-H(D) bend 1216 C-Cl stretch 761 Cl-C-Cl bend 262
3151 681 377 1240 786 268
2256 651 367 908 738 262
2321 661 374 923 759 267
description
The E modes are doubly degenerate.
HOD in heavy water, it is generally believed that upon excitation of the OH stretch the first step in the relaxation is dominated by transition to the HOD bend overtone,35-38 with the concomitant release of several hundred wave numbers of energy, which is accepted by low-frequency solvent modes. More complicated examples involve VER in neat chloroform and bromoform following CH fundamental excitation.39-44 Because of the dense vibrational level structure of these molecules, relaxation occurs to a number of different levels in the originally excited molecule, and a number of combinations of vibrational, rotational, and translational solvent modes serve as energy acceptors. In this paper we build on these last two systems by considering VER of the bend fundamental of dilute water in chloroform and d-chloroform. This is an interesting system for experimental and theoretical study because while there is only one relaxation pathway from the water bend fundamental (directly to the ground state), there are several possibilities for accepting the roughly 1600 cm-1 of energy released. Experimentally it was found that the VER rate increases by over a factor of 3 in going from d-chloroform to chloroform.45 Deuteration of the solvent has little effect on its translational and rotational motions, while the solvent vibrational level structure changes considerably. The strong solvent isotope effect along with the large energy gap clearly indicates a V-V mechanism for this VER process, although the precise nature of the accepting vibrational modes is at this point still unclear.45 Chloroform has nine vibrational modes, with frequencies as shown in Table 1.46 Deuteration shifts the fundamental frequencies, most significantly of the C-H(D) stretch and the bend. The fundamental closest to the 1601 cm-1 of the water bend is that of ν4, at 1216 or 908 cm-1, with an energy mismatch ∆E of about 400 or 700 cm-1 for chloroform or d-chloroform, respectively. This mode is therefore a likely candidate as an acceptor, with the remaining energy difference made up by lowfrequency solvent translational or rotational modes. Other candidates include combinations of solvent vibrations, leading to smaller values of ∆E. It is generally understood that a smaller ∆E will typically lead to a larger rate. However, another relevant quantity is ∆V, the total number of solute and solvent vibrational quanta destroyed or created during the relaxation process. A larger ∆V corresponds to a higher-order term in the interaction potential, generally resulting in a slower rate.10,47,48 While combinations of the solvent fundamentals can form numerous accepting states having fairly small ∆E, they have larger ∆V. Thus the effects of ∆E and ∆V tend to compete with each other. Theoretical approaches to calculating VER rates and to elucidating VER intramolecular pathways and intermolecular channels2,13,20,35-40,49-61 have for the most part invoked the system-bath-coupling (SBC) method.2 The solute vibrations are separated from the other degrees of freedom, and the Hamil-
tonian is partitioned accordingly into the system, which involves the solute vibrations, the bath, which includes the translation and rotation of the solute and all the motions of the solvent, and the coupling between system and bath. The transition rates between pairs of system eigenstates are given by Fermi’s Golden Rule as derived from first-order time-dependent perturbation theory in the coupling. In a simple VER problem, the latter can often be approximated as the product of the solute vibrational coordinate of interest and the force acting on it. In the more useful time-correlation function (TCF) formalism the transition rates are described by the Landau-Teller (LT) equation.1,62 The LT rate is proportional to the Fourier transform, evaluated at the relevant transition frequency, of the quantum-mechanical force-force TCF. Due to the difficulty in evaluating a quantum TCF in condensed phases, typically a molecular dynamics (MD) simulation is performed to extract the Fourier transform of the classical TCF, which is then multiplied by an appropriate quantum correction factor (QCF).63-67 The perturbative LT treatment of VER is valid only if the lifetime of the initially excited state is long compared to the bath correlation time, which is typically on the order of 2 ps; otherwise one should appeal to higher-order perturbation theories or the fluctuating LT approach.68,69 Several options are available for determining the force on the vibrational coordinate within the LT framework. In the conventional SBC approach, the interaction potential is described by an empirical force field, typically with Lennard-Jones (LJ) parameters and partial charges assigned to each atomic site. Atomic forces are dictated by the potential gradients, which are then projected onto the solute vibrational coordinate to generate the required force. In the extended system-bath-coupling (ESBC) approach the atomic charges of the solute are allowed to depend on the vibrational coordinate, therefore leading to an additional contribution to the force. This variation of the solute charges with vibrational coordinate has been found to play an essential role in the relaxation of the asymmetric stretch of azide in water.14,18 Recently a new optimized system-bath-coupling (OSBC) approach using a parametrized electronic structure method to calculate the force has been proposed and applied to the VER of azide in water; it is found that the OSBC approach gives theoretical lifetimes in good agreement with experiment.18 In what follows we perform theoretical calculations of the VER rate for the water bend fundamental in chloroform and d-chloroform based on the LT theory. The experimental lifetimes are about an order of magnitude larger than the bath correlation time, suggesting that treating this VER problem by first-order perturbation theory should be adequate. We first adopt the conventional SBC approach for the force calculations, using existing empirical force fields for the water and the solvent, together with Lorentz-Berthelot combining rules.39,70-72 Our VER rates using this approach, however, are not in good agreement with experiment, indicating that either the empirical force fields or the method of calculating the forces on the bend coordinate need to be improved. To this end we develop and parametrize a new intermolecular chloroform-water force field, using ab initio calculations on appropriately chosen chloroformwater dimers, within a novel self-consistent framework. The computed relaxation rates with this new potential, but with the traditional SBC method for calculating forces on the vibrational coordinate, are in reasonable agreement with the experimental results, and the trend for the solvent isotope effect is correct. Therefore we did not make further efforts to improve upon the SBC calculation of forces on the vibrational coordinate. We found that V-V channels contribute significantly to the water
392 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Lin et al.
bend relaxation in both chloroform and d-chloroform and that there is no single solvent vibrational state acting as a dominant energy acceptor.
The SBC method for VER calculations starts by partitioning the total Hamiltonian into system, bath, and coupling terms: H ˆ )H ˆS + H ˆB + H ˆ ′(b,q). The single solute vibrational coordinate, q, in this VER study of the water bend fundamental is simply the bending coordinate. The water stretches are not considered here because their fundamentals are roughly 2000 cm-1 higher than the bend, making the upward transitions very unlikely. The bath coordinates b involve the translations, rotations, and vibrations of the solvent along with the translation and rotation of the solute. A. Water Bend Vibration. The bend vibration of the water molecule is approximated by a one-dimensional harmonic rigid bender, with the system Hamiltonian
H ˆS ) -
H ˆ )H ˆS + H ˆB + H ˆ′
Tˆ i + Tˆ TN + Tˆ RN + W(b,0) ∑ i)1
(2)
and the harmonic frequency ω10 is taken to be 1601 cm-1, the experimental value of the water bend fundamental in chloroform.45 B. System-Bath Decomposition. The Hamiltonian for N 1 chloroform and 1 water molecules can generally be written as a sum of kinetic and potential energy terms
H ˆ ′ = -Fq
(3)
where Tˆ i is the kinetic energy operator for the ith molecule, with the Nth one being for water, and U(b,q) represents the total potential energy, which depends on both the bath coordinates b and the solute vibrational coordinate q. The kinetic energy operator for the water molecule is separated into translational, rotational, and vibrational terms
Tˆ N ) Tˆ TN + Tˆ RN -
2
2
p ∂ 2µr02 ∂q2
(4)
N-1
Wi(b) + V(b,q) ∑ i)1
(5)
(11)
where Wi(b) is the intramolecular potential of the ith solvent (chloroform) molecule, defined by
Wi(b) )
1
9
9
∑ ∑ FmnSimSin
2 m)1 n)1
(12)
where Sim are the symmetry coordinates of molecule i and the force constants Fmn are taken from the FF2 force field developed by Sibert and Rey.39 The fundamental frequencies based on this model are given in Table 1. For both chloroform isotopes the agreement with experiment is good enough for our purposes. Note that the model frequencies are harmonic, while the experimental frequencies are for the (anharmonic) fundamentals, which is one reason why the latter are smaller. According to eqs 6 and 9, the potential energy for the bath dynamics is then N-1
W(b,0) )
Wi(b) + V(b,0) ∑ i)1
(13)
while F is given by
F)-
∂V(b,q) | ∂q q)0
(14)
The intermolecular interaction term V(b,q) includes atom-atom LJ and Coulomb contributions and explicitly depends on q because the positions of the water atoms and hence the distances between water and chloroform atoms are functions of the bending coordinate. Using the chain rule, we can write
The total potential energy is written as
1 U(b,q) ) µr02ω102q2 2
(10)
Thus the bath potential energy is W(b,0), and the force F is a nonlinear function of (only) the bath coordinates. C. Solvent Vibrations and Calculation of the Force. W(b,q) in eq 5 consists of the vibrational potential energy of the solvent molecules and the intermolecular interactions
N
Tˆ i + U(b,q) ∑ i)1
(9)
and
(1)
where r0 is the equilibrium O-H bond length, the bending coordinate q is defined by the difference between the H-O-H bond angle θ and its equilibrium value q ≡ θ - θ0, the reduced mass µ is related to the atomic masses and the equilibrium bond angle by73
2 2 1 ) + [1 - cos(θ0)] µ mH mO
N-1
H ˆB )
W(b,q) )
2
p ∂ 1 + µr02ω102q2 2 2 2 2µr0 ∂q
H ˆ )
(8)
where H ˆ S is given in eq 1
II. Theoretical Considerations
2
Now the Hamiltonian can be rewritten as
F)
∑j
∂rj |q)0 ∂q
fj‚
(15)
Performing a Maclaurin expansion in q of the latter term in eq 5 leads to
where j ) 1, 2, 3 for the three water atoms, rj is the position of the jth atom, and
W(b,q) ) W(b,0) - Fq + ‚‚‚
fj ) -∇jV(b,0)
(6)
with
F≡-
∂W(b,q) | ∂q q)0
(7)
(16)
is the total force on atom j from the bath variables during the course of the simulation. The derivatives of the water atomic positions with respect to q are most easily evaluated in the molecule-fixed frame, defined
VER of Dilute H2O in Chloroform
J. Phys. Chem. B, Vol. 112, No. 2, 2008 393
TABLE 2: Potential Parameters for the DH Chloroform Model,71,79 the SPC/E Water Model,72 and the New Chloroform-Water Potentiala DH chloroform
SPC/E water
new
geometry
RCH RCCl θHCCl
1.100 1.758 107.6
ROH θHOH
1.000 109.47
ROH θHOH
0.95718 104.523
charges
QC QH QCl
0.179 0.082 -0.087
QO QH
-0.8476 0.4238
QO QH
-0.6592 0.3296
LJ parameters
σC σH σCl C H Cl
σO σH O H
3.40 2.2 3.44 4.26 × 102 8.31 × 101 1.25 × 103
3.166 0.0 6.49 × 102 0.0
CCO AHO AHH AClO AClH
4.665 × 106 2.818 × 107 2.168 × 107 6.877 × 109 9.718 × 106
a Bond distances R and LJ σ are in angstroms, bond angles θ are in degrees, charges Q are in e, values are in J/mol, CCO is in Å6 J/mol, and A are in Å12 J/mol.
as follows: The z-axis bisects the H-O-H angle and points from the center of mass to the oxygen, the x-axis is in the molecular plane, and the y-axis is perpendicular to the molecular plane. The water atoms are labeled such that the hydrogen atom with a positive x-component is 1, the other hydrogen is 2, and the oxygen is 3. Defining
cjk )
( | ) ∂rj ∂q
to be the kth component (k ) 1, 2, 3, corresponding to x, y, z in the molecule-fixed frame) of the derivative for the jth atom, we find that (keeping the center of mass fixed and the molecule nonrotating)73
(
)
where
R ) r0 cos(θ0/2)/2
(18)
β ) mOr0 sin(θ0/2)/2(2mH + mO)
(19)
γ ) mHr0 sin(θ0/2)(2mH + mO)
(20)
Finally, F can then be written as
F)
fjkcjk ∑ jk
(21)
where fjk is the kth component (in the molecule-fixed frame) of fj. D. LT Formalism and Quantum Correction Factors. The transition rate from the bend fundamental to the vibrational ground state is, within the TCF formulation of first-order timedependent perturbation theory, given by the LT formula
ˆ (ω10)/p2 k10 ) C
(22)
with
C ˆ (ω) ) q102
k10 = Q(ω10)C ˆ cl(ω10)/p2
(17)
q)0 k
R 0 β c ) -R 0 β 0 0 -γ
average (with F(t) ) eiHˆ Bt/pF e-iHˆ Bt/p), and the brackets denote the quantum statistical equilibrium average with respect to the bath. The quantum force-force TCF can be approximated by its classical counterpart, whose Fourier transform is then multiplied by a QCF (discussed below) to ensure that detailed balance is obeyed, so that
∫-∞∞ dt eiωt〈δF(t)δF(0)〉
(23)
In the above q10 ) 〈1|q|0〉 ) xp/2µω10r02, with |0〉 and |1〉 being the ground and the first excited eigenstates of the system, δF(t) ) F(t) - 〈F〉 is the difference between the force and its
(24)
where
C ˆ cl(ω) ) q102
∫-∞∞ dt eiωt〈δF(t)δF(0)〉cl
(25)
Various forms for QCFs have been proposed.2,63,66,74-77 A common choice is the harmonic QCF, QH(ω), which is exact if the bath is harmonic and the force is linear in the bath coordinates.63,66 For the VER of the water bend at room temperature, QH(ω10) = 8. If the relaxation mechanism can be identified, then other QCFs may be more appropriate.78 For example, for a vibrational mode relaxing by a vibrationtranslation/rotation V-T/R channel, involving the excitation of multiple quanta of low-frequency modes, the harmonic/Schofield QCF, QHS(ω10), may be a better choice.78 Thus, if the water bend relaxes through this multiphonon fashion, then the classical force-force TCF should be multiplied by QHS(ω10) = 20. E. Simulation Details. Molecular dynamics (MD) simulations of 255 chloroform molecules plus a single water are performed at 295 K to extract the classical force-force TCFs. Cubic periodic boundary conditions are applied, where the box size of is chosen such that the density corresponds to the experimental value of 1.484 g/cm3 for chloroform at 295 K.79 One of the major challenges in computing the LT rates involves finding a satisfactory model potential W(b,0). Among the significant number of extant chloroform71,79-82 and water72,83-90 force fields, the chloroform potential parametrized by Dietz and Heinzinger (DH)71,79-82 and the SPC/E water model72 (see Table 2 for the parameters) were chosen for their abilities to give nice overall agreement with liquid-phase experimental properties. They were applied to the present system with the use of the Lorentz-Berthelot combining rules70 to determine the intermolecular LJ interactions between chloroform and water. As discussed below, we also perform simulations with a new potential model developed herein. During the simulations, the water geometry is fixed using the SHAKE algorithm91 such that the q ) 0 condition of the bath Hamiltonian in eq 9 is fulfilled. As described in section II.C, the truncated intramolecular force field by Sibert and Rey39
394 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Lin et al.
Figure 2. Normalized distributions of F and Fa for 2000 dimers sampled from an MD simulation of 255 flexible DH chloroform and 1 SPC/E water molecules.
Figure 1. Spectral densities for SPC/E water in DH chloroform.
is implemented to describe solvent vibrations. Intermolecular interactions, to which the shifted-force potential method with a cutoff radius of 15 Å is applied, include site-site LJ and Coulomb contributions. The leapfrog algorithm70 is used to integrate the equations of motion with a time step of 0.2 fs. The conventional SBC forces on the water bending coordinate, as described in eq 21, are calculated in the production runs of approximately 8 million time steps (on the order of 1 ns). The flexible solvent model implemented in the MD simulations allows for VER contributions from intermolecular energy transfer. To help analyze the VER mechanism, we also performed simulations with a rigid solvent, in which case the solvent geometry is constrained by SHAKE. III. Results and Discussion A. VER Results for the DH-SPC/E Potential. A classical MD simulation, with the bath Hamiltonian in eq 9, yields a force trajectory that is used to calculate C ˆ cl(ω) defined by eq 25. The latter (divided by p2) is plotted in Figure 1 for both chloroform and d-chloroform, using both rigid and flexible solvent models. These results are for the DH-SPC/E model described above. As shown in the figure, when the solvent is rigid the spectrum displays a sloping plateau from 200 to 500 cm-1 due to librational motions, with a subsequent approximately exponential decay consistent with the energy-gap law. The low-frequency motions of chloroform and d-chloroform are expected to be very similar, and indeed, the spectral densities for rigid chloroform and rigid d-chloroform are nearly indistinguishable. When vibrations are included in the simulation model, for chloroform there is an increase in the magnitude of the spectral density at the frequencies of the stretches and bend (Table 1). Similar features are noticed in the case of flexible d-chloroform, with vibrational frequencies red-shifted due to the deuteration. In the case of flexible chloroform solvent, when evaluated at ω10 ) 1601 cm-1, C ˆ cl(ω)/p2 (the fully classical approximation to k10) yields 0.000361 ps-1, which is more than 300 times smaller than the experimental value of 0.12 ps-1. For the rigid solvent, the rate is even smaller, by about 2 orders of magnitude, showing the pronounced effect of including solvent vibrations. To estimate the magnitude of solvent quantum effects, the rate for the flexible solvent should be multiplied by Q(ω10), as
discussed earlier. Different choices for the QCF lead to different estimates for the VER rate. Probably a reasonable upper bound comes from choosing QHS(ω10) ≈ 20, which gives a rate that is still 15 times smaller than the experimental result. Thus we must conclude that our calculations of the spectral density are not sufficiently accurate. Assuming that the harmonic approximation to the system Hamiltonian is adequate, it follows that the error lies in the high-frequency Fourier transform of the classical force-force TCF. Thus either the SBC approach to calculating the force is not sufficiently accurate, our approximation to the bath Hamiltonian is not adequate, or both. To examine this further, we can write 〈δF(t)δF(0)〉cl ) 〈δF2〉clψ(t), where ψ(t) is the normalized classical force-force TCF. Therefore C ˆ cl(ω10) ) q102〈δF2〉clψ ˆ (ω10). It would appear that either 〈δF2〉cl or ψ ˆ (ω10) (or both) are too small. The former involves equilibrium fluctuations of the force, Boltzmann weighted by the bath potential W(b,0), while the latter involves high-frequency components resulting from short-time dynamics of F as generated by the bath potential. Thus we must have significant errors associated with the solvent vibrational or intermolecular components of the bath Hamiltonian or with the calculations of the force. To clarify this situation, we will make use of ab initio calculations. To this end from an MD simulation with the flexible solvent, 2000 chloroform-water dimers were selected under the criterion that C-O distances were smaller than 6 Å. For each dimer single-point interaction energies Ea were calculated at the B3LYP level with the aug-cc-pVDZ basis set,92 counterpoise-corrected for basis set superposition error. We also calculated the force on the bend coordinate from ab initio theory. By specifying the atoms of chloroform as ghost atoms, intramolecular contributions to the force were computed and subtracted to obtain the ab initio force Fa. For example, in Figure 2, we compare the distribution of ab initio forces Fa for the 2000 dimers with the distribution of forces F using the SBC approach and the DH-SPC/E potential. The distributions are qualitatively similar. The standard deviation of F (2.30 × 10-21 J) is smaller than that of Fa (3.34 × 10-21 J) by about 30%, which is insufficient to account for the discrepancy between theory and experiment described above. Therefore it seems likely that the error lies in ψ ˆ (ω10). This could in turn result from problems in either our description of the intramolecular solvent vibrations or the empirical DH-SPC/E model for the intermolecular interactions.
VER of Dilute H2O in Chloroform
J. Phys. Chem. B, Vol. 112, No. 2, 2008 395
Figure 3. Binding energies from empirical potential, E, vs those from ab initio calculations, Ea, for dimers sampled from an MD simulation of 255 flexible DH chloroform and 1 SPC/E water molecules (200 of the 2000 data points are shown). The diagonal line indicates perfect correlation.
Figure 4. Binding energies from empirical potential, E, vs those from ab initio calculations, Ea, for dimers sampled from an MD simulation of 255 flexible DH chloroform and 1 water molecules using the new potential (200 of the 2000 data points are shown). The diagonal line indicates perfect correlation.
To investigate the latter, in Figure 3 we show a scatter plot of the empirical interaction energies E from the DH-SPC/E force field versus the ab initio energies Ea for a subset of the dimers. A root-mean-square (rms) deviation of 3.62 kJ/mol, more than kT at room temperature, is found between E and Ea. Furthermore the empirical binding energies are systematically too low. (Note that we are not suggesting any major flaws with the DH potential, but rather, with the Lorentz-Berthelot combining rules used for determining the chloroform-water interaction and to a lesser extent the appropriateness of the SPC/E geometry and charges.) Therefore, we decided to develop a new model potential for the chloroform-water interaction (retaining the DH model for the chloroform-chloroform interactions). B. New Water-Chloroform Potential. In the SPC/E model, the intramolecular geometry of water is modified from that of the gas phase to account for average changes in structure in the liquid state. For dilute water in chloroform, however, among other things the small vibrational frequency shifts from the gasphase values indicate that the properties of the solute are more similar to those of gas-phase water than to those of bulk water.45,46 Thus the development of the chloroform-water model potential starts with adopting the gas-phase geometry for water and adjusting the atomic partial charges to reproduce the experimental gas-phase dipole moment93,94 (Table 2). Although the chloroform-water interaction potential will still consist of LJ and Coulomb contributions, the Lorentz-Berthelot combining rules for the LJ interactions will be discarded. It then takes two parameters to describe the A/R12 repulsion and -C/R6 attraction between each type of atomic pair, incurring a total of 12 parameters to define the LJ interaction between chloroform and water! To reduce the number of parameters, only five types of interactions are considered: C-O attraction, H(chloroform)-O repulsion, H(chloroform)-H(water) repulsion, Cl-O repulsion and Cl-H(water) repulsion. The C-O attraction is designed to represent the overall induced dipole-induced dipole interaction between the two molecules, while the four repulsions between outlying atoms arise from steric interactions. It was determined that including other interactions did not significantly improve the results. In general, parametrization of an intermolecular potential is often done by fitting to ab initio binding energies of certain “typical” dimer geometries. Instead of preselecting these dimer geometries, herein dimer configurations are sampled from MD simulations. First, an MD simulation was performed for 255
flexible DH chloroform molecules and 1 rigid water molecule with an O-H bond length of 0.95718 Å, an H-O-H bond angle of 104.523°, and QH ) 0.3296e (Table 2).93,94 The LJ interaction between chloroform and water in this simulation was still determined by the DH and SPC/E LJ parameters with the implementation of the mixing rules. Using the same procedure described earlier in the section, 2000 dimer configurations were sampled, and their ab initio binding energies Ea were calculated. The five LJ parameters were optimized by fitting the empirical binding energy E to Ea for each dimer with the help of a genetic algorithm.95,96 The new potential with these LJ parameters is now different from that used in the MD simulation, and in that sense the new potential is inconsistent with the configurations used in its fit. To make the potential and configurations selfconsistent, we then ran a new MD simulation using this updated set of LJ parameters, and another 2000 dimer configurations were sampled. By fitting E to Ea, a new set of parameters with a better rms deviation were obtained. This process was iterated until convergence was reached (that is, the LJ parameters no longer changed appreciably). The final optimized parameters are given in Table 2. To demonstrate the improvement in the potential, 2000 dimers were sampled from the MD simulation with the new model, and the E-Ea comparison is plotted in Figure 4. By comparing Figure 4 with Figure 3, a notable improvement can be observed (rms deviation is 1.41 kJ/mol, in contrast with 3.62 kJ/mol), and the deviations are much less systematic. We also wanted to make sure that this new potential produces an accurate distribution of forces, and so in Figure 5 we show the new F and Fa distributions for these dimers. F and Fa again share rather similar standard deviations (2.45 × 10-21 J for F and 2.76 × 10-21 J for Fa). Interestingly the F distribution of the new potential does not reveal an appreciable difference from that of the old potential (Figures 2 and 5). C. VER Results for the New Potential. Armed with this new potential, we now repeat the MD simulations and spectral density calculations for both chloroform and d-chloroform solvents. Results for chloroform and d-chloroform are plotted in Figures 6a and 6b, respectively. The spectral densities shown in Figures 1 and 6 all start at approximately the same value and decay similarly until about 500 cm-1. After that, however, the spectral densities for the new potential decay significantly more slowly, confirming that the new intermolecular potential generates different short-time dynamics of the forces. (Note that
396 J. Phys. Chem. B, Vol. 112, No. 2, 2008
Lin et al. TABLE 3: Some Possible Solvent Accepting States chloroform
Figure 5. Normalized distributions of F and Fa for 2000 dimers sampled from an MD simulation of 255 flexible DH chloroform and 1 water molecules using the new potential.
Figure 6. Spectral densities for new potential model: (a) chloroform and (b) d-chloroform.
the two figures have different scales on the y-axes.) For flexible chloroform, when evaluated at 1601 cm-1, C ˆ cl(ω10)/p2 leads to 0.00735 ps-1; for flexible d-chloroform C ˆ cl(ω10)/p2 ) 0.00293 ps-1. As discussed in sections I and II, k10 can be estimated from the product of C ˆ cl(ω10)/p2 and an appropriate QCF. If the VER process involved the excitation of one quantum of a nearresonant solvent vibration, which is approximately harmonic, then the model of linear coupling to a harmonic bath would be appropriate, indicating the use of the harmonic QCF.63,78 For other more complicated scenarios involving excitation of several quanta of solvent vibrations (see below), then this model becomes less appropriate but is still often a reasonably good approximation. In any case, for the moment we use the harmonic QCF, whose value for the water bend relaxation is QH(ω10) ) 7.82. With this choice T1 ) 1/k10 ) 17.4 and 43.6 ps for chlorofrom and d-chloroform, respectively, with a ratio of about 2.5 between the lifetimes in the two solvents. The computed lifetimes are approximately twice as long as the experimental
state
frequency (cm-1)
ν4 2ν2 ν2 + ν4 ν2 + ν5 ν3 + ν4 ν3 + ν5 2ν4 ν4 + ν 5 ν4 + ν6 2ν5
1216 2 × 668 668 + 1216 668 + 761 366 + 1216 366 + 761 2 × 1216 1216 + 761 1216 + 262 2 × 761
d-chloroform
∆E (cm-1)
frequency (cm-1)
∆E (cm-1)
385 265 -283 172 19 474 -831 -376 123 79
908 2 × 651 651 + 908 651 + 738 367 + 908 367 + 738 2 × 908 908 + 738 908 + 262 2 × 738
693 299 42 212 326 496 -215 -45 431 125
results of T1 ) 8.5 ps in chloroform and T1 ) 28.5 ps in d-chloroform. The theoretical isotope ratio is in the same direction as experiment and with roughly the same magnitude (about 3.4 for experiment). These theoretical results, while still not quantitative, are a significant improvement over those from the DH-SPC/E model. D. Solvent Accepting Modes. To begin to understand the nature of the accepting modes, for both isotopes we can compare the spectral densities at 1601 cm-1 for rigid and flexible solvents, as shown in Figure 6. In chloroform the flexible spectral density is more than a factor of 3 higher than the rigid spectral density, and in d-chloroform this ratio is more like a factor of 2. This shows that in both cases V-V processes make a significant contribution but also that V-T/R channels are active as well. We can use this information to make a slightly more informed estimate of the VER rate by decomposing the classical spectral density into two terms, one for the V-T/R channel (this comes from the rigid solvent) and one for the V-V channel (from the difference between the flexible and the rigid spectral densities). To estimate the total VER rate, we multiply each of these contributions by its own QCF:37 A reasonable approach is to use QHS for the V-T/R channel and QH for the V-V channel.78 This approach gives T1 ) 12.8 and 24.7 ps for the lifetimes in chloroform and d-chloroform, respectively. These lifetimes, expected to be somewhat more accurate than above, are indeed in slightly better agreement with experiment, although the isotope ratio of 1.9 is somewhat too small. To understand which chloroform vibrations comprise the accepting modes (for the V-V channel), we perform another simulation where the solvent bond lengths are fixed by the SHAKE algorithm. One sees in Figure 6 that when evaluated at ω10 these results reduce the spectral densities from their flexible solvent values but not all the way to the rigid solvent values. This indicates that the C-Cl stretching modes (ν2 and ν5 in Table 1) are involved in the relaxation process, but that other solvent vibrational modes also participate. As discussed in the Introduction, an efficient channel is likely to have a rather small total vibrational quantum number change ∆V. In addition, the energy mismatch between the solute vibration and the accepting solvent vibrational state should be small enough for the low-frequency translations and rotations to compensate. Regarding the latter, from the spectral densities for both chloroform and d-chloroform, we see that the edge of the low-frequency plateau occurs at about 500 cm-1. Thus, for example, one finds that for chloroform there is one possible accepting state with ∆V ) 2 and ∆E e 500 cm-1, eight possible states with ∆V ) 3 and ∆E e 500 cm-1, and many more such possible states with ∆V g 4. For d-chloroform, there are no states with ∆V ) 2 (and with ∆E e 500 cm-1) and nine states with ∆V ) 3. The candidates with ∆V e 3, together with their frequencies and energy mismatches, are summarized in Table
VER of Dilute H2O in Chloroform
J. Phys. Chem. B, Vol. 112, No. 2, 2008 397 correlation time of this bond vector, for H217O, can also be measured with quadrupolar relaxation experiments,98 but to the best of our knowledge this has not been done for water in chloroform. It has, however, been measured recently in carbon tetrachloride;99 this correlation time is 91 fs. For D2O (O-D bond vector), the correlation time is over a factor of 2 longer in chloroform than in carbon tetrachloride.97 Assuming the same is true for H217O (out-of-plane bond vector), this would again be consistent with the conclusion that our theoretical orientation times in chloroform are somewhat too short. IV. Concluding Remarks
Figure 7. Second-rank orientational TCF for water in chloroform. P2 is the Legendre polynomial of order 2; the unit vector e bisects the H-O-H bond angle.
3. Constraining the C-Cl bond lengths is tantamount to eliminating all states involving ν2 and ν5. From the table, we see that all of the remaining states involve ν4, the bend. Our conclusion then is that this VER process involves the excitation of both C-Cl stretching and bending modes. Moreover, from examining the spectral densities at 1601 cm-1 in Figure 6 for both chloroform and d-chloroform, we conclude that contributions from both the C-Cl stretches and the Cl-C-H(D) bends are increased in chloroform, implying that the observed isotope effect is related to frequency changes of all of these modes. It is possible that states with ∆V ) 4 or higher contribute to the relaxation process; a similar analysis for all states with ∆V ) 4 and ∆E e 500 cm-1 leads to the identical conclusions. An alternative complementary approach for identifying the nature of the accepting modes comes from instantaneousnormal-mode analysis.59-61 Diagonalization of the massweighted Hessian matrix of the potential energy yields the frequencies and eigenvectors of the instantaneous normal modes for a given configuration. The influence spectrum for the VER process of interest can then be computed using these instantaneous normal modes, and it is possible to determine the contribution to the relaxation process from specific solvent modes. It would be interesting to see if such an analysis confirms our mechanistic conclusions. E. Water Rotational Dynamics. In addition to measuring VER rates, Seifert et al. measured the anisotropy decay for the water bend fundamental in both chloroform and d-chloroform.45 Specifically, they measured the long-time decay of the secondrank orientation TCF for the unit vector bisecting the H-O-H bond angle (the direction of the bend transition dipole), finding identical relaxation times of 2.7 ( 0.7 ps for both solvents. Our calculated TCF (also indistinguishable for the two solvents) is shown in Figure 7 and has a long-time decay of 1.2 ( 0.1 ps, which is roughly a factor of 2 faster than experiment. We can also make a connection with NMR quadrupolar relaxation experiments on dilute D2O in chloroform.97 These experiments measure the (integrated) correlation time of the second-rank orientation TCF for the O-D bond vector, giving a value of 230 fs. From our simulations of H2O in chloroform, we obtain 76 fs for the same quantity (but for the O-H bond vector). To make a direct comparison with experiment, we reran our simulations with D2O, finding a slightly longer value of 83 fs. Thus our theoretical correlation time is again over a factor of 2 too short. We also calculated the correlation time for the out-of-plane bond vector, which for H2O is 102 fs. The
Vibrational energy relaxation rates of the water bend fundamental in chloroform and d-chloroform are calculated using the Landau-Teller formula. The empirical chloroform force field developed by Dietz and Heinzinger and the SPC/E water model are used, together with the standard Lorentz-Berthelot combining rules, to determine the intermolecular solvent-solute interaction, and we used the intramolecular force field of Sibert and Rey to describe the solvent vibrations. The required force on the water bend coordinate was calculated from the empirical intermolecular potential in the standard manner. It is found that this treatment of the chloroform-water interaction does not lead to lifetimes in good agreement with experiment. We then developed a new chloroform-water intermolecular potential model using a novel scheme, whereby model potential energies were fit to those from ab initio calculations for dimers obtained from an MD simulation, and the process is iterated by using the new potential parameters to generate a new sample of dimers until self-consistency is achieved. This new potential gives relaxation rates in reasonable agreement with the experimental results and also gives a correct trend for the solvent isotope effect. Our analysis shows that both V-T/R and V-V channels contribute to the relaxation in this system. In terms of the specific vibrational accepting modes of the solvent, by performing simulations where the solvent is rigid and also where only the C-Cl stretches are frozen, we conclude that both C-Cl stretches and Cl-C-H(D) bends function as solvent accepting modes and the changes in the frequencies of all of these modes contribute to the observed solvent isotope effect. Thus we find that even for this relatively simple solvent there are a number of possible solvent accepting modes, and many of these channels are operative. For VER of other solutes in other polyatomic solvents with five or more atoms, the same conclusion is likely to be true. VER lifetimes for the water stretch in chloroform and d-chloroform have also been measured experimentally; the lifetimes are 72 and 89 ps, respectively.100 This process is of course substantially more complicated than that of the water bend due to the presence of lower-frequency vibrational states of the water solute. Recently the relaxation pathway has been determined experimentally: The OH stretch relaxes via the bend overtone and bend fundamental to the ground state.101 With a generalization of the system Hamiltonian to include the water stretches, one could perform similar theoretical calculations of the water stretch VER process, which, among other things, would provide an additional test of the new water-chloroform potential developed herein. Acknowledgment. The authors are grateful for support from the National Science Foundation (Grant Nos. CHE-0446666 and CHE-0615165) and from an American Chemical Society Petroleum Research Foundation grant.
398 J. Phys. Chem. B, Vol. 112, No. 2, 2008 References and Notes (1) Owrutsky, J. C.; Raftery, D.; Hochstrasser, R. M. Annu. ReV. Phys. Chem. 1994, 45, 519. (2) Oxtoby, D. W. AdV. Chem. Phys. 1981, 47 (Part 2), 487. (3) Fujisaki, H.; Bu, L.; Straub, J. E. AdV. Chem. Phys. 2005, 130 (Part B), 179. (4) Harris, C. B.; Smith, D. E.; Russell, D. J. Chem. ReV. 1990, 90, 481. (5) Stratt, R. M.; Maroncelli, M. J. Phys. Chem. 1996, 100, 12981. (6) Okazaki, S. AdV. Chem. Phys. 2001, 118, 191. (7) Rey, R.; Møller, K. B.; Hynes, J. T. Chem. ReV. 2004, 104, 1915. (8) Crim, F. F. J. Phys. Chem. 1996, 100, 12725. (9) Crim, F. F. Acc. Chem. Res. 1999, 32, 877. (10) Elles, C. G.; Crim, F. F. Annu. ReV. Phys. Chem. 2006, 57, 273. (11) Zare, R. N. Science 1998, 279, 1875. (12) Faltermeier, B.; Protz, R.; Maier, M.; Werner, E. Chem. Phys. Lett. 1980, 74, 425. (13) Everitt, K. F.; Skinner, J. L. J. Chem. Phys. 1999, 110, 4467. (14) Morita, A.; Kato, S. J. Chem. Phys. 1998, 109, 5511. (15) Li, M.; Owrutsky, J.; Sarisky, M.; Culver, J. P.; Yodh, A.; Hochstrasser, R. M. J. Chem. Phys. 1993, 98, 5499. (16) Hamm, P.; Lim, M.; Hochstrasser, R. M. Phys. ReV. Lett. 1998, 81, 5326. (17) Zhong, Q.; Baronavski, A. P.; Owrutsky, J. C. J. Chem. Phys. 2003, 118, 7074. (18) Li, S.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2006, 125, 244507. (19) Hamm, P.; Lim, M.; Hochstrasser, R. M. J. Chem. Phys. 1997, 107, 10523. (20) Rey, R.; Hynes, J. T. J. Chem. Phys. 1998, 108, 142. (21) Harris, A. L.; Berg, M.; Harris, C. B. J. Chem. Phys. 1986, 84, 788. (22) Ashihara, S.; Huse, N.; Espagne, A.; Nibbering, E. T. J.; Elsaesser, T. J. Phys. Chem. A 2007, 111, 743. (23) Wang, Z.; Pang, Y.; Dlott, D. D. J. Phys. Chem. A 2007, 111, 3196. (24) Lindner, J.; Vo¨hringer, P.; Pshenichnikov, M. S.; Cringus, D.; Wiersma, D. A.; Mostovoy, M. Chem. Phys. Lett. 2006, 421, 329. (25) Ashihara, S.; Huse, N.; Espagne, A.; Nibbering, E. T. J.; Elsaesser, T. Chem. Phys. Lett. 2006, 424, 66. (26) Huse, N.; Ashihara, S.; Nibbering, E. T. J.; Elsaesser, T. Chem. Phys. Lett. 2005, 404, 389. (27) Bakker, H. J.; Lock, A. J.; Madsen, D. Chem. Phys. Lett. 2004, 384, 236. (28) Pakoulev, A.; Wang, Z.; Pang, Y.; Dlott, D. D. Chem. Phys. Lett. 2003, 380, 404. (29) Lock, A. J.; Bakker, H. J. J. Chem. Phys. 2002, 117, 1708. (30) Pakoulev, A.; Wang, Z.; Dlott, D. D. Chem. Phys. Lett. 2003, 371, 594. (31) Wang, Z.; Pakoulev, A.; Pang, Y.; Dlott, D. D. Chem. Phys. Lett. 2003, 378, 281. (32) Woutersen, S.; Emmerichs, U.; Nienhuys, H.-K.; Bakker, H. J. Phys. ReV. Lett. 1998, 81, 1106. (33) Gale, G. M.; Gallot, G.; Lascoux, N. Chem. Phys. Lett. 1999, 311, 123. (34) Schwarzer, D.; Lindner, J.; Vo¨hringer, P. J. Chem. Phys. 2005, 123, 161105. (35) Rey, R.; Hynes, J. T. J. Chem. Phys. 1996, 104, 2356. (36) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 5827. (37) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 1623. (38) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 119, 3840. (39) Sibert, E. L.; Rey, R. J. Chem. Phys. 2002, 116, 237. (40) Ramesh, S. G.; Sibert, E. L. J. Chem. Phys. 2006, 124, 234501. (41) Ramesh, S. G.; Sibert, E. L. J. Chem. Phys. 2006, 125, 244512. (42) Ramesh, S. G.; Sibert, E. L. J. Chem. Phys. 2006, 125, 244513. (43) Graener, H.; Zu¨rl, R.; Hofmann, M. J. Phys. Chem. B 1997, 101, 1745. (44) Seifert, G.; Zu¨rl, R.; Patzlaff, T.; Graener, H. J. Chem. Phys. 2000, 112, 6349. (45) Seifert, G.; Patzlaff, T.; Graener, H. J. Chem. Phys. 2004, 120, 8866. (46) Herzberg, G. Molecular Spectra and Molecular Structure. II. Infrared and Raman Spectra of Polyatomic Molecules; Van Nostrand: New York, 1950. (47) Gruebele, M.; Bigwood, R. Int. ReV. Phys. Chem. 1998, 17, 91. (48) Gruebele, M.; Wolynes, P. G. Acc. Chem. Res. 2004, 37, 261. (49) Shiga, M.; Okazaki, S. J. Chem. Phys. 1999, 111, 5390. (50) Nesbitt, D. J.; Hynes, J. T. J. Chem. Phys. 1982, 77, 2130. (51) Miller, D. W.; Adelman, S. A. J. Chem. Phys. 2002, 117, 2672. (52) Li, S.; Thompson, W. H. J. Phys. Chem. A 2003, 107, 8696. (53) Li, S.; Thompson, W. H. Chem. Phys. Lett. 2004, 383, 326. (54) Gulmen, T. S.; Sibert, E. L. J. Phys. Chem. A 2004, 108, 2389. (55) Gulmen, T. S.; Sibert, E. L. J Phys. Chem. A 2005, 109, 5777. (56) Ka, B. J.; Geva, E. J. Phys. Chem. A 2006, 110, 13131. (57) Leitner, D. M. Phys. ReV. Lett. 2001, 87, 188102.
Lin et al. (58) Fujisaki, H.; Straub, J. E. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6726. (59) Ladanyi, B. M.; Stratt, R. M. J. Phys. Chem. A 1998, 102, 1068. (60) Deng, Y.; Ladanyi, B. M.; Stratt, R. M. J. Chem. Phys. 2002, 117, 10752. (61) Graham, P. B.; Matus, K. J.; Stratt, R. M. J. Chem. Phys. 2004, 121, 5348. (62) Zwanzig, R. J. Chem. Phys. 1961, 34, 1931. (63) Bader, J. S.; Berne, B. J. J. Chem. Phys. 1994, 100, 8359. (64) Rostkier-Edelstein, D.; Graf, P.; Nitzan, A. J. Chem. Phys. 1997, 107, 10470. (65) Egorov, S. A.; Skinner, J. L. Chem. Phys. Lett. 1998, 293, 469. (66) Egorov, S. A.; Everitt, K. F.; Skinner, J. L. J. Phys. Chem. A 1999, 103, 9494. (67) Egorov, S. A.; Rabani, E.; Berne, B. J. J. Phys. Chem. B 1999, 103, 10978. (68) Bakker, H. J. J. Chem. Phys. 2004, 121, 10088. (69) Gulmen, T. S.; Sibert, E. L. J. Chem. Phys. 2005, 123, 204508. (70) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U. K., 1987. (71) Dietz, W.; Heinzinger, K. Ber. Bunsen-Ges. Phys. Chem. 1984, 88, 543. (72) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (73) Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations; Dover: New York, 1980. (74) Berens, P. H.; White, S. R.; Wilson, K. R. J. Chem. Phys. 1981, 75, 515. (75) Schofield, P. Phys. ReV. Lett. 1960, 4, 239. (76) Egelstaff, P. A. AdV. Phys. 1962, 11, 203. (77) Berne, B. J.; Jortner, J.; Gordon, R. J. Chem. Phys. 1967, 47, 1600. (78) Skinner, J. L.; Park, K. J. Phys. Chem. B 2001, 105, 6716. (79) Dietz, W.; Heinzinger, K. Ber. Bunsen-Ges. Phys. Chem. 1985, 89, 968. (80) Evans, M. W. AdV. Mol. Relax. Interact. Processes 1982, 24, 123. (81) Kovacs, H.; Kowalewski, J.; Laaksonen, A. J. Phys. Chem. 1990, 94, 7378. (82) Chang, T.-M.; Dang, L. X.; Peterson, K. A. J. Phys. Chem. B 1997, 101, 3413. (83) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (84) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (85) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (86) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (87) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (88) Stern, H. A.; Rittner, F.; Berne, B. J.; Friesner, R. A. J. Chem. Phys. 2001, 115, 2237. (89) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (90) Soper, A. K. J. Phys.: Condens. Matter 2005, 17, 3273. (91) Ciccotti, G.; Ryckaert, J.-P. Comput. Phys. Rep. 1986, 4, 345. (92) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (93) Benedict, W. S.; Gailar, N.; Plyler, E. K. J. Chem. Phys. 1956, 24, 1139. (94) Lide, D. R. Handbook of Chemistry and Physics, 73rd ed.; CRC Press: Boca Raton, FL, 1992. (95) Michalewicz, Z. Genetic Algorithms + Data Structures ) EVolution Programs, 3rd ed.; Springer-Verlag: Berlin, 1996. (96) Carroll, D. L. GA DriVer; CU Aerospace: Urbana, IL, 2001. (97) Nakahara, M.; Wakai, C. J. Chem. Phys. 1992, 97, 4413. (98) Ropp, J.; Lawrence, C.; Farrar, T. C.; Skinner, J. L. J. Am. Chem. Soc. 2001, 123, 8047. (99) Goodnough, J. A.; Goodrich, L.; Farrar, T. C. Preprint, 2007. (100) Graener, H.; Seifert, G. J. Chem. Phys. 1993, 98, 36. (101) Seifert, G.; Patzlaff, T.; Graener, H. J. Chem. Phys. 2006, 125, 154506.
