5682
J. Phys. Chem. A 2010, 114, 5682–5688
Vibrational Energy Relaxation Rates via the Linearized Semiclassical Method without Force Derivatives Francisco X. Va´zquez, Irina Navrotskaya, and Eitan Geva* Department of Chemistry, UniVersity of Michigan, Ann Arbor, Michigan 48109-1055 ReceiVed: February 3, 2010; ReVised Manuscript ReceiVed: March 23, 2010
A new computational scheme for calculating vibrational energy relaxation rate constants is presented that is based on applying the linearized semiclassical approximation to the symmetrized force-force correlation function. Unlike the previous scheme of Geva and Shi [Shi, Q.; Geva, E. J. Phys. Chem A 2003, 107, 9059], which was based on applying the linearized semiclassical approximation to the standard force-force correlation function, the new scheme does not involve a power expansion of the initial force in terms of the Wigner transform integration variable ∆ and as a result is more accurate and does not require force derivatives as input. The main disadvantages of the new scheme are its slower convergence rate in comparison to the Shi-Geva scheme and its ambiguity with respect to the choice of configuration around which the local harmonic approximation is performed. The computational feasibility and accuracy of the new approach is demonstrated on a number of benchmark models, including the highly challenging case of nonpolar diatomic liquids. It is observed that performing the local harmonic approximation around the configuration used for computing the initial force yields the best agreement with experiment. I. Introduction Vibrational energy relaxation (VER) is the process by which an excited vibrational mode releases its excess energy into other intramolecular and/or intermolecular degrees of freedom (DOF). VER is prevalent in many systems of fundamental, technological, and biological importance and plays a central role in determining chemical reactivity.1-49 Recent theoretical and computational studies of VER have been mostly based on the Landau-Teller formula,15,50,51 which gives the VER rate constant in terms of the Fourier transform (FT), at the vibrational transition frequency, of the quantum-mechanical autocorrelation function of the fluctuating force exerted on the relaxing mode by the other DOF. Importantly, replacing the quantum-mechanical force-force correlation function (FFCF) by its classical counterpart can only be justified in cases where the vibrational transition frequency is significantly smaller than kBT/p. Indeed, discrepancies by many orders of magnitude have been reported between experimentally measured VER rates and those calculated based on purely classical molecular dynamics (MD) simulations when this condition is not met.7,52-67 However, the exact calculation of real-time quantum-mechanical correlation functions for general anharmonic many-body systems remains far beyond the reach of currently available computer resources.68 A popular approach for dealing with this difficulty in the case of VER calls for multiplying the classical VER rate constant by a frequency-dependent quantum correction factor (QCF).7,52-67 A variety of different approximate QCFs have been proposed in the literature. However, the choice of QCF is somewhat arbitrary, and estimates obtained from different QCFs can differ by orders of magnitude. Thus, the development of more rigorous methods for computing VER rate constants is highly desirable. In a series of recent papers,69-73 we have proposed a more rigorous approach for calculating VER rate constants, which is based on the linearized semiclassical (LSC) method. The LSC method is based on linearizing the forward-backward action * To whom correspondence should be addressed.
in the path-integral expression for the quantum-mechanical FFCF with respect to the difference between the forward and backward paths.74 It should be noted that the same approximation can be derived in several other ways, including linearization of the forward-backward action in the semiclassical initial value representation approximation for the correlation function,75-81 and starting from the Wigner representation formalism.82 The resulting LSC approximation for a real-time quantum-mechanical correlation function of the form: ˆ ˆ CAB(t) ) Tr(AˆeiHt/pBˆe-iHt/p)
(1)
is given by LSC CAB (t) )
1 (2πp)N
∫ dQ0 ∫ dP0 AW(Q0, P0)BW(Qt(Cl), Pt(Cl)) (2)
where, N is the number of DOF, Q0 ) (Q0(1), · · · , Q0(N)) and P0 ) (P0(1), · · · , P0(N)) are the corresponding coordinates and momenta,
AW(Q, P) )
〈
∫ d∆ e-iP∆/p Q + ∆2 |Aˆ|Q - ∆2
〉
(3)
) is the Wigner transform of the operator Aˆ,83 and Q(Cl) t (Cl) Q(Cl) ) P(Cl) t (Q0, P0) and Pt t (Q0, P0) are propagated classically with the initial conditions Q0 and P0. The LSC approximation is known to be exact at t ) 0, at the classical limit, and for harmonic systems. At the same time, it also provides a convenient starting point for introducing computationally feasible schemes to calculate quantum-mechanical time correlation functions. The main disadvantage of the LSC approximation is that it can only capture quantum effects at short times.77 However, it should be noted that in
10.1021/jp1010499 2010 American Chemical Society Published on Web 04/22/2010
Vibrational Energy Relaxation Rates
J. Phys. Chem. A, Vol. 114, No. 18, 2010 5683
condensed phase systems in general, and in the case of highfrequency VER in particular, the quantities of interest are often dominated by the behavior of short-lived correlation functions at relatively short times. In practice, applying the LSC approximation, eq 2, requires the calculation of the phase-space integrals underlying the Wigner transforms, eq 3. The numerical calculation of the integral in eq 3 can be extremely difficult in practice because of the oscillatory phase factor, e-iP∆/p, in the integrand. In the case of the standard FFCF, Aˆ and Bˆ in eq 2 correspond to ˆ ˆ δFˆe-βH/Z and δFˆ, respectively, where δFˆ ) Fˆ - Tr[e-βHFˆ]/Z, Fˆ corresponds to the force exerted on the relaxing mode by the bath, which consists of all the other DOF except for the relaxing ˆ is the free-bath Hamiltonian, and Z ) Tr[e-βHˆ] is the mode, H free-bath quantum partition function. Neglecting centrifugal forces, which are often found to be of only secondary importance, so that Fˆ is a function of only the ˆ )]W(Q(Cl) bath coordinates, one finds that [F(Q t ) reduces into its ). Hence, the only remaining classical counterpart, F(Q(Cl) t computational challenge has to do with calculating the following Wigner transform:
ˆ )e-βHˆ]W(Q0, P0) ) [δF(Q
〈
〉
| |
∫ d∆ e-iP · ∆/p Q0 + ∆2 e-βH Q0 - ∆2 δF(Q0 + ∆2 ) 0
ˆ
(4)
Within the original LHA-LSC method,69 one proceeds by introducing the following two additional approximations: ˆ ˆ (1) Evaluating 〈Q0 + (∆/2)|e-βH|Q0 - (∆/2)〉/〈Q0|e-βH|Q0〉 within the framework of the local harmonic approximation (LHA), which amounts to expanding the Hamiltonian to second order around Q0. (2) Expanding F(Q0 + (∆/2)) to second order in terms of ∆ around Q0. The resulting Gaussian integral over ∆ can then be solved analytically to yield the LHA-LSC approximation for the FFCF:
CLHA-LSC(t) )
∫
ˆ
〈Q0 |e-βH |Q0〉 Z
∫
N
∏
(
1 (j) 2 j)1 R πp (Cl) [δF(Q0) + D(Q0, Pn,0)]δF(Qt )
dQ0
dPn,0
)
1/2
[
exp -
(j) 2 (Pn,0 )
p2R(j)
]
(5) Here, Pn ) Pn(Q0) ) (Pn(1)(Q0), ..., Pn(N)(Q0)), where {Pn(k)(Q0) ) ∑Nl)1Tl, k(M(l))-1/2P(l)} are the normal mode momenta that emerge from diagonalizing the Hessian matrix underlying the quadratic expansion of the bath potential energy around Q ) Q0 and
[
(j)
(j)
Ω (Q0) βpΩ (Q0) R(j) ) R(j)(Q0) ) coth p 2
]
(6)
where {(Ω(k))2(Q0)} are the eigenvalues of the Hessian matrix. The term D(Q0, Pn,0) is given by:69 N
D(Q0, Pn,0) ) -i
(k) F˜k′Pn,0
∑ pR(k)
k)1
N
+
′′ F˜k,k
∑ 4R(k)
k)1
-
N
(k) (l) ′′ Pn,0 F˜k,l Pn,0
k,l)1
2p2R(k)R(l) (7)
∑
and its calculation requires as input the first and second derivatives of the force with respect to the bath coordinates, that is, N
F˜k′ )
∑ (M(l))-1/2Tl,kFl′; l)1
N
N
∑ ∑ (M(i)M(j))-1/2Ti,lTj,kFi,j′′
′′ ) F˜k,l
(8)
i)1 j)1
where
Fk′ )
∂F ∂Q(k)
|
Q)Q0
;
′′ ) Fk,l
∂2F ∂Q(k)∂Q(l)
|
Q)Q0
(9)
Quantum effects enter the expression for CLHA-LSC(t), eq 5, in two ways: (1) Nonclassical initial sampling of the bath coordinates and momenta. (2) The nonclassical term D(Q0, Pn,0) which can be shown to vanish in the classical limit. Both effects turned out to be important for the applications considered in refs 69-73. Importantly, the nonclassical term D(Q0, Pn,0), which can be traced back to the first- and secondorder terms in the above-mentioned expansion of F(Q0 + ∆/2) in powers of ∆, was found to play an important role in enhancing the VER rate. Although the calculation of the first and second derivatives of the bath-induced force with respect to the bath coordinates, which is required as input for calculating D(Q0, Pn,0), is in principle straightforward, in practice it can become increasingly more tedious and cumbersome with the increasing complexity of the force fields. Furthermore, unlike the second derivatives required for calculating the Hessian matrix underlying the LHA, the above-mentioned force derivatives are usually not part of the output of standard packages for performing MD simulations. In this paper we develop an alternative scheme for calculating VER rate constants that, while still based on LSC and LHA, avoids the above-mentioned expansion of the bath-induced force in terms of ∆ and as a result is more accurate and does not require force derivatives as input. The remainder of this paper is organized as follows. We outline the theory of VER and derive the new force-derivativefree approach in Section II, test the new approach on several benchmark systems, including liquid oxygen and nitrogen, in Section III, and conclude in Section IV. II. Theory Consider the following general quantum-mechanical Hamiltonian of a vibrational mode linearly coupled to a bath: 2 ˆ tot ) pˆ + V(qˆ) + H 2µ
N
ˆ (j)
2
) ∑ (P (j) 2M
ˆ ) - qˆF(Q ˆ) + V(Q
(10)
j)1
Here, qˆ, pˆ, µ, and V(qˆ) are the relaxing mode coordinate, momentum, reduced mass, and bath-free vibrational potential, ˆ ) are the coordinates, ˆ , Pˆ, {M(1), ..., M(N)} and V(Q respectively; Q momenta, masses, and potential energy of the bath DOF, ˆ ) is the potential force exerted by the bath respectively; and F(Q on the relaxing mode. The standard expression of the population relaxation rate constant between the first-excited and ground vibrational states
5684
J. Phys. Chem. A, Vol. 114, No. 18, 2010
Va´zquez et al.
is usually given in terms of the FT, at the transition frequency, of the standard FFCF:15,50,51
k10
(11)
〈 ∫ 〈 | | | | 〉 〈 | | | | 〉 〈 | | | | 〉
)
Here, ω10 is the transition frequency and
∫ dQ δF(Q′) Q0 + ∆2 e-βH/2 Q′〉〈Q′ e-βH/2 Q0 - ∆2 ˆ
(12)
ˆ
ˆ
ˆ
ˆ ∆ -βHˆ/2 ∆ e Q′〉〈Q′ e-βH/2 Q0 2 2 ˆ
ˆ
Q0 e-βH/2 Q′〉〈Q′ e-βH/2 Q0
(19)
is the FT of the FFCF
C(t) )
The ratio in the integrand,
ˆ ˆ ˆ 1 Tr[δFˆe-βHeiHt/pδFˆe-iHt/p] Z
(13)
ˆ ) is the free bath Hamiltonian. ˆ ) ΣNj)1(Pˆ(j))2/(2M(j)) + V(Q where H The population relaxation rate constant k10 can also be cast in terms of a variety of alternative expressions that are completely equivalent to eq 11. One such expression, which will turn out particularly useful for our purpose, can be based on the symmetrized FFCF given by:
Cs(t) )
ˆ ˆ ˆ ˆ 1 Tr[e-βH/2δFˆe-βH/2eiHt/pδFˆe-iHt/p] Z
(14)
The fact that the FTs of Cs(t) and C(t) are related in a simple manner, namely
〈
(15)
1 eβpω10/2C˜s(ω10) 2µpω10
(16)
In the next step, we attempt to evaluate Cs(t) within the framework of the LSC and LHA approximations. We first note that within the LSC approximation (see eq 2):
CsLSC(t) ) 1 (2πp)N
∫ dQ0 ∫ dP0 [e-βH/2δFˆe-βH/2]W(Q0, P0)δF(Qt(Cl)) ˆ
ˆ
(17) Clearly, the main challenge in calculating CsLSC(t) lies in evaluating the following Wigner transform: ˆ
ˆ
[eβH/2δFˆe-βH/2]W(Q0, P0) )
〈
|
|
∫ d∆ e-iP ∆/p Q0 + ∆2 e-βH/2δFˆe-βH/2 Q0 - ∆2 0
ˆ
ˆ
〉
(18)
Using the closure relation ∫dQ′|Q′〉〈Q′| ) 1ˆ we can rewrite the matrix element in the integrand in the following from:
| |
ˆ ∆ -βHˆ/2 ∆ e Q'〉〈Q' e-βH/2 Q0 2 2 ˆ /2 -βH
Q0 e
| | 〉 ˆ /2 -βH
Q'〉〈Q' e
〉
(20)
Q0
can then be evaluated within the LHA. To this end, we expand the bath potential energy to second order around the yet to be specified point Q ) Q*, followed by analytically solving the resulting Gaussian integral over ∆, to obtain the LHA-LSC approximation of Cs(t):
CsLHA-LSC(t)
∫ dPn ∏ j)1
then allows us to obtain the following alternative expression for the population relaxation rate constant:
| | 〈 | |
Q0 +
N
C˜(ω) ) eβpω/2C˜s(ω)
k10 )
〉
dQ δF(Q′) Q0 e-βH/2 Q′〉〈Q′ e-βH/2 Q0
Q0 +
∫-∞∞ dt eiωtC(t)
| 〉 | | | |
|
∆ -βHˆ/2 ˆ -βHˆ/2 ∆ e δFe Q0 2 2
Q0 +
)
1 ˜ ) C(ω10) 2µpω10
C˜(ω) )
〈
(
ˆ
ˆ
〈Q0 |e-βH/2 |Q'〉〈Q'|e-βH/2 |Q0〉 ) dQ0 dQ' Z (j) 2 (P ) 1/2 1 n exp 2 (j) δF(Q')δF(Qt(Cl)[Q0, P0]) (j) 2 R πp pR (21)
∫
∫
)
[ ]
Here, Pn ) Pn(Q*) ) (Pn(1)(Q*), ..., Pn(N)(Q*)) are the normal mode momenta that emerge from diagonalizing the Hessian matrix underlying the quadratic expansion of the bath potential energy around Q ) Q* and {R(j)} are as in eq 6, where {(Ω(k))2} are the eigenvalues of the Hessian matrix. Equation 21 represents the main theoretical result of this paper. Similarly to eq 5, it reduces to the classical FFCF in the classical limit. However, it should be noted that the highfrequency FT of the standard and symmetrized FFCFs is expected to exhibit nonclassical behavior and may therefore differ quite significantly from one another. It should also be noted that in deriving eq 21 we did not resort to expanding δF(Q′) to second order around Q0 (or any other point). In this respect, eq 21 is in fact less approximate than eq 5. At the same time calculating eq 21 does not require force derivatives as input! One consequence of this is that the D(Q0, Pn,0) term, eq 7, which played such an important role in accounting for the nonclassical behavior of the standard FFCF, is no longer present. Instead, the nonclassical behavior of the symmetrized FFCF is accounted for by the following attributes: (1) Nonclassical sampling of bath coordinates and momenta (similar to that in eq 5). (2) The initial force is not calculated at the initial position used to generate the classical trajectory leading to the force at a later time t. (3) The factor eβpω10/2 (see eq 16), which actually coincides with the so-called Schofield QCF.65
Vibrational Energy Relaxation Rates It should also be noted that, up to this point, we have not specified the point Q* around which the LHA is performed. This choice will obviously affect the actual frequencies of the instantaneous normal modes and thereby the values of R(j), which will in turn affect the values of the initial momenta that are (t). Below, we will consider sampled and the resulting CLHA-LSC s three possible and rather natural choices for Q*, namely: (1) Q* ) Q0, which corresponds to the initial configuration for the classical trajectory that generates the force at time t, δF(Q(Cl) t ). (2) Q* ) Q′, which corresponds to the configuration used for calculating the initial force, δF(Q′). (3) Q* ) Qc, which correspond to the centroid configuration of the cyclic imaginary time path integral. It should be noted that the distinction between the first two choices is absent from eq 5 since in this equation the initial force is calculated at the initial configuration that is used as an initial condition for the classical trajectory. As will be shown below, at least for the applications that we have considered, performing the LHA around Q′ appear to yield results that are similar to these obtained via eq 5 and also compare well with experiment. It should be noted that a similar calculation of the symmetrized flux-side quantum-mechanical correlation function within the LSC and LHA approximations was recently presented by Liu and Miller in the context of calculating reaction rate constants.84 Although the calculation follows a similar procedure in both cases, it is worth noting that the actual challenges involved with calculating reaction rate constants are somewhat different than these associated with calculating VER rate constants (e.g., the reaction rate constant can be expressed in terms of a correlation function that involves quantities that only depend on the reaction coordinate and does not require calculating its high-frequency FT). Nevertheless, in both cases one needs to make a choice regarding the bath configuration around which the LHA is performed. The authors of ref 84 have chosen to perform the LHA around Q* ) (Q0 + Q′)/2, which is similar, but not identical, to the above-mentioned choice in Q* ) Qc. The calculations of CsLHA-LSC(t) reported below were based on eq 21 and were carried out following the algorithm outlined below: (1) Perform an imaginary-time path integral moleculardynamics or Monte Carlo simulation (PIMD and PIMC, respectively),85,86 and use it to sample the initial configuration, Q0, and the configuration at which the initial force is calculated, Q′. To this end, it should be noted that within the context of a PIMD/PIMC simulation, each DOF is represented by a cyclic necklace of P beads labeled 0, 1, 2, ..., P - 1. Assuming that P is even, Q0 is identified with the configuration of the beads labeled 0, and Q′ is identified with the configuration of the beads labeled P/2. (2) Perform a LHA around either Q′ or Q0 or Qc, find the normal-mode frequencies, {Ω(k)}, and corresponding transformation matrix, {Tl, k}, and use it to calculate {R(k)} and sample the initial (normal-mode) momenta, {Pn(k)}. (3) Calculate Q(Cl) via a classical MD simulation for each t sampled initial configuration Q0 and normal mode momenta Pn,0, and time correlate δF(Q(Cl) t ) with δF(Q′). (4) Repeat steps 1-3 until reaching the desired convergence. III. Applications A. Exponential Coupling to a Harmonic Bath. The first model that we will consider involves a bath consisting of uncoupled harmonic oscillators of different frequencies,
J. Phys. Chem. A, Vol. 114, No. 18, 2010 5685 N
ˆ ) H
∑ j)1
(
(Pˆ(j))2 1 ˆ (j))2 + M(j)(ω(j))2(Q (j) 2 2M
)
(22)
and a force that is exponential in the bath coordinates,
[
ˆ ) ) exp Fˆ(Q
N
∑ c(j) j)1
2M(j)ω(j) ˆ (j) Q p
]
(23)
The exact quantum-mechanical FFCF can be obtained analytically for this model and is given by:64,87,88
C(t) ) eB(0)(eB(t) - 1)
(24)
where
B(t) )
∫0∞ dω Γ(ω){[n(ω) + 1]e-iωt + n(ω)eiωt} (25) N
Γ(ω) )
∑ (c(k))2δ(ω - ω(k))
(26)
k)1
and n(ω) ) [exp(βpω) - 1]-1. It should be noted that the LSC approximation is exact when the system is harmonic. Thus, for this model, eq 21 actually coincides with the exact symmetrized FFCF. To this end, we note that for the bath harmonic Hamiltonian in eq 22 one can show that: ˆ ˆ 1 〈Q |e-βH/2 |Q'〉〈Q'|e-βH/2 |Q0〉 ) Z 0
exp[-R(j)(Q0(j))2 + (Q'(j))2] exp
∏ ( Mπpω N
[
j)1 (j)
(j)
(j)
)
×
]
2 M (ω(j)) Q0(j)Q'(j) p sinh(βpω(j) /2)
(27)
Substituting this result into eq 21 then yields:
[
CsLHA-LSC(t) ) exp
N
∑ (c(j))2 j)1
{ ( ) coth
βpω(j) + 2
} ]
cos(ω(j)t) -1 sinh(βpω(j) /2)
(28)
The exact standard FFCF, eq 24, can then be obtained from eq 28 by using the identity C(t) ) Cs(t + (iβp/2)). This should be contrasted with eq 5, which involves the additional expansion of the force δF(Q0 + (∆/2)) around Q0, to second order in terms of ∆. Indeed, in this case it can be shown that69
[
CLHA-LSC(t) ) eBR(0) (eBR(t) - 1) + ieBR(t)BI(t) 1 BR(t) 2 e BI (t) 2
]
(29)
where BR(t) and BI(t) are the real and imaginary parts of B(t), respectively, which clearly differs from the exact result, eq 24.
5686
J. Phys. Chem. A, Vol. 114, No. 18, 2010
Va´zquez et al.
Figure 1. The exact C˜(ω) in the case of a harmonic bath and a force that is exponential in the bath coordinates. Note that in this case C˜(ω) ˜ LHA-LSC(ω), C ˜ Cl(ω), ˜ LHA-LSC (ω)eβpω/2. Also shown are C coincides with C s and eβpω/2C˜Cl(ω) (Schofield QCF).
C˜(ω) as obtained via eqs 24 and 29 are compared to the corresponding fully classical result in Figure 1. The calculations were performed with the following spectral density,
Γ(ω) ) 2λ
( )
ωR ω2 exp ωR+1 ω2c c
(30)
and for the following values of the parameters: λ ) 0.20, R ) 3, and βpωc ) 4.0. Figure 1 shows that although eq 5 is in rather good agreement with the exact results for a wide range of frequencies, discrepancies start appearing with increasing frequency, which can be attributed to the additional approximation embodied in the above-mentioned expansion of δF(Q0 + (∆/2)) to second order with respect to ∆. At the same time, eq 21 is seen to coincide with the exact result at all frequencies, thereby testifying to its less-approximate nature. B. Breathing Sphere Model. We next consider the VER of a spherically symmetric solute (so-called breathing sphere) in a monatomic solvent.42,69,89,90 Calculations were performed on a two-dimensional liquid and under the assumption that the solvent atoms and the solute have the same mass and interact via identical pair potentials of the Lennard-Jones type (for more details see ref 69). The calculations reported below were performed using a square simulation cell containing 81 atoms at a reduced density and temperature of F* ) 0.70 and T* ) 0.68, respectively. Periodic boundary conditions and a potential cutoff at 3σ have been employed. PIMD simulations were performed with 16 beads per atom. The classical equations of motion were solved via the velocity Verlet algorithm, and the FFCF was obtained by averaging over 30 000 trajectories of length 8.0 ps each. In Figure 2a, we compare CsLHA-LSC(t) as obtained via eq 21 with the real part of CLHA-LSC(t) as obtained via eq 5, and the corresponding classical FFCF, CCl(t). The LHA underlying the calculation of CsLHA-LSC(t) was performed around the abovementioned three different configurations, namely Q0, Q′, and Qc. It should be noted that the initial value and rate of decay of CsLHA-LSC(t) are smaller than these of CCl(t). This should be contrasted with the initial value and rate of decay of the real part of CLHA-LSC(t) which are larger in comparison to their classical counterparts. As a result, C˜sLHA-LSC(ω) < C˜Cl(ω) < C˜LHA-LSC(ω) throughout the entire range of frequencies. Thus, in the case of CsLHA-LSC(t) the combined effect of nonclassical initial sampling and the fact that the initial force is calculated at Q ) Q′, rather than at Q0, is to diminish the value of C˜sLHA-LSC(ω) relative to its classical counterpart. In contrast, in the case of CLHA-LSC(t), the combined effect of the very same
nonclassical sampling and a nonclassical initial force that includes the D(Q0, Pn,0) term turns out to enhance C˜LHA-LSC(ω) relative to its classical counterpart. However, it should be ˜ LHA-LSC ˜ (ω) ) eβpω/2C (ω) (see eq 16). Hence, remembered that C s a more meaningful comparison is between C˜LHA-LSC(ω) and eβpω/2C˜sLHA-LSC(ω). Indeed, as we show in Figure 2b, eβpω/2C˜sLHA-LSC(ω) is comparable to C˜LHA-LSC(ω), regardless of the choice of LHA. It should also be noted that eβpω/2C˜Cl(ω) is ˜ LHA-LSC ˜ LHA-LSC(ω) and eβpω/2C (ω). significantly larger than both C s It should also be noted that the Schofield QCF significantly overestimates of the VER rate (see Figure 2b). An interesting observation is related to the numerical ˜ LHA-LSC(ω). As ˜ LHA-LSC (ω) in comparison to C convergence of C s is well-known, in the absence of resonance with other vibrations, C˜(ω) is expected to decay asymptotically with frequency in an exponential manner. As a result, it becomes increasingly more difficult to average out the statistical noise accompanying any real-life simulation, which is needed in order to calculate the increasingly small value of C˜(ω) at high frequencies. In fact, it is a common practice to obtain C˜(ω) at high frequencies by extrapolating the exponential gap law, which usually emerges at low frequencies, to much higher frequencies.87,91 Now, because C˜LHA-LSC(ω) ) eβpω/2C˜sLHA-LSC(ω), its high frequency value is significantly smaller than C˜(ω), which in turn makes it even more difficult to obtain a converged result. Thus, one expects that obtaining a converged estimate of C˜s(ω) at a given frequency will be computationally more demanding than obtaining a converged estimate of C˜(ω) at the same frequency. However, in practice, we found that it is still possible to obtain converged estimates of C˜s(ω) on a wide enough range of frequencies for the exponential gap law to emerge and serve as basis for extrapolation to higher frequencies. Another interesting observation is the dependence of C˜s(ω) on the LHA. Although the results of different types of LHA are relatively similar, there are discernible differences that become increasingly larger upon extrapolation to the high frequency domain. In practice, we have observed that, at least for the systems considered here, the best agreement with the estimates obtained via the original LHA-LSC method, as well with experiment when available (see next section), were obtained when the LHA was performed around Q′, that is around the configuration used for calculating the initial force. This can be explained by the fact that performing the LHA around Q′ leads to the most nonclassical sampling of momenta, thereby maximizing the quantum effect, which in this case corresponds to lowering the value of C˜s(ω). C. Neat Liquid Oxygen and Nitrogen. Finally, we consider the VER of the homonuclear diatomic molecules O2 and N2 in the corresponding neat liquid. Calculations were performed on a three-dimensional liquid at 77K, with densities of 22.64 nm-3 and 17.37 nm-3 for O2 and N2, respectively, using site-site pair potentials of the Lennard-Jones type (for more details see refs 70 and 71). All calculations were preformed with 108 molecules contained in a cubic cell with periodic boundary conditions. PIMD simulations were performed with 16 beads per atom. Each FFCF was averaged over 50 000 trajectories, (t), the calculations each of length 40 ps. In the case of CLHA-LSC s were performed using the above-mentioned three types of LHA around Q′, Q0, and Qc. In Figures 3a and 4a, we compare CsLHA-LSC(t) as obtained via eq 21 with the real part of CLHA-LSC(t) as obtained via eq 5, and the corresponding classical CCl(t), in the cases of liquid oxygen and liquid nitrogen, respectively. The LHA underlying (t) was performed around the abovethe calculation of CLHA-LSC s
Vibrational Energy Relaxation Rates
J. Phys. Chem. A, Vol. 114, No. 18, 2010 5687
Figure 2. (a) A comparison of CsLHA-LSC(t) (with the LHA performed around either Q′ or Q0 or Qc), Re [CLHA-LSC(t)] and CCl(t) for the breathing sphere model; (b) A semilog plot of C˜sLHA-LSC(ω)eβpω/2 (with the LHA performed around either Q′ or Q0 or Qc) C˜LHA-LSC(ω), C˜Cl(ω), and eβpω/2C˜Cl(ω) (Schofield QCF) for the breathing sphere model.
Figure 3. (a) A comparison of CsLHA-LSC(t) (with the LHA performed around Q′), Re [CLHA-LSC(t)] and CCl(t) for neat liquid oxygen; (b) A semilog plot of C˜sLHA-LSC(ω) (with the LHA performed around Q′), C˜LHA-LSC(ω), C˜Cl(ω), and eβpω/2C˜Cl(ω) (Schofield QCF) for neat liquid oxygen. Solid lines were obtained from the simulation, and dashed lines correspond to extrapolations. The experimental value is indicated by an asterisk (*) at the vibrational transition frequency of oxygen (1553 cm-1). The insert shows C˜sLHA-LSC(ω) on a narrower range of frequencies.
˜ LHA-LSC(ω) ) e-βpω/2C ˜ LHA-LSC due to the fact that C (ω) becomes s increasingly smaller than C˜(ω) with increasing frequency. (3) The best agreement with the predictions of the LHALSC method as well as with experiment (see table) are obtained when the LHA is performed around Q′. (4) The agreement of the prediction obtained via eq 21 with experiment in the case of oxygen is somewhat inferior to that obtained via eq 5 in the case of oxygen, and vice versa in the case of nitrogen. The differences between the two schemes can be attributed to the fact that they originate from applying the LSC and LHA approximations to different forms of the correlation functions (i.e., standard vs symmetrized). However, it is important to note that both predictions are far superior to the corresponding classical predictions that underestimate the VER rates by many orders of magnitude. (5) The Schofield QCF alone significantly overestimates the VER rate. However, the agreement with experiment can be significantly improved by replacing C˜Cl(ω) by C˜sLHA-LSC(ω). IV. Concluding Remarks
Figure 4. (a) A comparison of CsLHA-LSC(t) (with the LHA performed around Q′), Re [CLHA-LSC(t)] and CCl(t) for neat liquid nitrogen; (b) A semilog plot of C˜sLHA-LSC(ω) (with the LHA performed around Q′), C˜LHA-LSC(ω), C˜Cl(ω), and eβpω/2C˜Cl(ω) (Schofield QCF) for neat liquid nitrogen. Solid lines were obtained from the simulation and dashed lines correspond to extrapolations. The experimental value is indicated by an asterisk (*) at the vibrational transition frequency of nitrogen (2327 cm-1). The insert shows C˜sLHA-LSC(ω) on a narrower range of frequencies.
mentioned three different configurations, namely Q0, Q′, and Qc. The corresponding frequency domain results are shown in Figures 3b and 4b. The results are seen to follow trends similar to these observed in the case of the breathing sphere model, namely: (1) The initial value and rate of decay of CsLHA-LSC(t) are smaller than these of CCl(t), which are in turn smaller than these ˜ LHA-LSC ˜ Cl(ω) (ω) < C of the real part of CLHA-LSC(t). As a result, C s LHA-LSC ˜ (ω) throughout the entire range of frequencies.
