In static quantum chemistry, IR spectra are calculated as described in section 2.4.2. As soon as the normal modes of a molecule are known, the IR intensities can be computed by taking the derivatives of the dipole moment along the normal coordinates. In an MD simulation, these derivatives are not directly accessible. An alternative approach is, however, possible by the following line of thought: If the MD simulation is in equilibrium, all modes are excited at once, so the oscillation of the dipole moment contains contributions from all modes at the
same time. Each of these contributions oscillates at the frequency of the corresponding mode, suggesting to employ a Fourier transform to obtain a frequency-resolved representation of these oscillations. In the same manner as power spectra (see section 3.2.1), these spectra show peaks at all vibrational frequencies of the molecule, but this time, the intensity is proportional to the amplitudes of the dipole oscillations at these frequencies. In linear approximation, the latter are proportional to the derivatives of the dipole moment along the normal coordinates (cf. equation (2.58)), so the result of the Fourier transform is just the IR spectrum of the system.
For a more detailed mathematical insight, the one-dimensional classical harmonic oscillator is considered again. Using the mass-weighted coordinateq = √
mx, its trajectory is given by (see section 2.4.1)
q(t) =√
mx0cos(ω0t+φ), (3.48)
wherex0is the amplitude,ω0is the eigenfrequency, andφis the phase. The dipole momentµis assumed to depend linearly on the coordinate, so the Taylor expansion around the equilibrium position can be truncated after the linear term:
µ(q)=µ0+ ∂µ
∂q
!
0
q, (3.49)
whereµ0is the dipole moment at the equilibrium position. Inserting the trajectory yields for the time-dependent dipole moment
µ(t)=µ0+ ∂µ
∂q
!
0
√mx0cos(ω0t +φ). (3.50)
For power spectra (see section 3.2.1), it was discussed that it is more convenient to use the time derivative of the coordinate instead of the coordinate itself. Transferring this to the dipole moment bears the additional advantage that the equilibrium dipole moment, which does not carry information relevant for the IR intensities, is removed:
µ˙(t) =− ∂µ
∂q
!
0
√mx0ω0sin(ω0t+φ). (3.51)
To remove also the dependence on the initial conditions of the MD simulation in terms of the phaseφ, the autocorrelation function is calculated (see appendix B.1):
Rµ˙(t)= 1 2
∂µ
∂q
!2
0
mx20ω02cos(ω0t). (3.52)
The Fourier transform of this expression can be written as F(Rµ˙)(ω)=
Z ∞
−∞
Rµ˙(t)exp(−iωt)dt = 1 2
∂µ
∂q
!2
0
mx02ω02π (δ(ω−ω0)+δ(ω+ω0)). (3.53) Due to the properties of the delta distribution, the total integral of the spectrum is
Z ∞
−∞
F(Rµ˙)(ω)dω = ∂µ
∂q
!2
0
mx02ω02π. (3.54)
Recalling that the amplitude is related to the average kinetic energy byx0= √
4hEkini/m/ω0 (see section 3.2.1), this is equal to
Z ∞
−∞
F(Rµ˙)(ω)dω= ∂µ
∂q
!2
0
·4πhEkini= ∂µ
∂q
!2
0
·2πkBT. (3.55) This means that the integral of the Fourier transform of the dipole moment autocorrelation in a harmonic oscillator is proportional to the average kinetic energy (or the temperature) and the squared derivative of the dipole moment along the mass-weighted coordinate at the equi-librium position. The latter is the same term as the one appearing in static quantum chemistry (see equation (2.59)). Comparing these two expressions leads to the following definition of the IR spectrum in MD simulations:
A(ω)= NA 24πε0c2kBT
Z ∞
−∞
µ(τ˙ )·µ(τ˙ +t)
τexp(−iωt)dt. (3.56) In this expression, the vector autocorrelation of the dipole moment is introduced to switch over to the general three-dimensional case. The prefactor ensures that—for a harmonic os-cillator with a linear dependence of the dipole moment on the position—the integral of the spectrum is equal to the integral absorption coefficient in static quantum chemistry. For com-parison to experimental data later on, it is more convenient to use the wavenumber-dependent representation
A(ν˜)= NA 12ε0ckBT
Z ∞
−∞
µ(τ˙ )·µ(τ˙ +t)
τexp(−2πicνt˜ )dt, (3.57) which is the one actually calculated in the Travis implementation.
The concept of IR spectra from the Fourier transform of the dipole moment autocorrelation function has been known in the literature for a long time261. The usual derivation starts with Fermi’s golden rule (2.56) and transforms the formula for the IR absorption coefficient from the Schrödinger picture to the Heisenberg picture of quantum mechanics51,148. In the resulting
expressions, the quantum correlation functions are approximated by classical correlation func-tions. Since these do not satisfy the detailed balance condition262, several quantum correction factors have been discussed52,67,75,139,140,263,264, and it has been concluded that the best choice is the factor called “harmonic correction”. Although different variants of the constant prefactor have been applied, this correction factor yields the same frequency dependence as in equa-tion (3.56). Here, this formula is directly derived in an alternative fashion by comparing the results for a harmonic oscillator. This clarifies the close relation of MD and static calculations regarding vibrational spectra.
Replacing the dipole moment by the polarizability, the same derivation can in principle be carried out for Raman spectra. The only point requiring special attention is that the polar-izability is a second-order tensor and the static quantum chemistry formulas contain various combinations of its elements (see section 2.4.3). The general rule is to replace each quadratic term containing polarizability derivatives along the normal coordinates by the corresponding autocorrelation function of the polarizability time derivatives. Choosing the prefactor in such a way that the integral of the MD spectrum agrees with the differential Raman scattering cross section in static quantum chemistry (for a harmonic oscillator with a linear dependence of the polarizability on the position), the polarized frequency-dependent differential Raman scatter-ing cross sections for fixed molecular orientation with respect to the laboratory coordinate system are given by
Ik(ω)= h¯ 64π3ε02c4kBT
(ωin−ω)4 ω
1−exp
−k¯hω
BT
Z ∞
−∞
Dα˙xx(τ)α˙xx(τ+t)E
τ exp(−iωt)dt (3.58) and
I⊥(ω)= h¯ 64π3ε02c4kBT
(ωin−ω)4 ω
1−exp
−k¯hω
BT
Z ∞
−∞
Dα˙xy(τ)α˙xy(τ +t)E
τ exp(−iωt)dt. (3.59) For random molecular orientation with respect to the laboratory coordinate system, the ex-pressions read as
Ik(ω)= h¯ 64π3ε20c4kBT
(ωin−ω)4 ω
1−exp
−khω¯
BT
45a(ω)+4γ(ω)
45 (3.60)
and
I⊥(ω)= h¯ 64π3ε20c4kBT
(ωin−ω)4 ω
1−exp
−khω¯
BT
3γ(ω)
45 (3.61)
with the isotropic contribution a(ω)=Z ∞
−∞
*α˙xx(τ)+α˙yy(τ)+α˙zz(τ) 3
α˙xx(τ +t)+α˙yy(τ+t)+α˙zz(τ+t) 3
+
τ
exp(−iωt)dt (3.62) and the anisotropic contribution
γ(ω)= Z ∞
−∞
"
1 2
D α˙xx(τ)−α˙yy(τ) α˙xx(τ+t)−α˙yy(τ +t) E
τ
+ 1 2
D α˙yy(τ)−α˙zz(τ) α˙yy(τ+t)−α˙zz(τ+t) E
τ
+ 1 2
D α˙zz(τ)−α˙xx(τ) α˙zz(τ +t)−α˙xx(τ +t) E
τ
+3D
α˙xy(τ)α˙xy(τ +t)E
τ +3D
α˙yz(τ)α˙yz(τ +t)E
τ
+3D
α˙zx(τ)α˙zx(τ+t)E
τ
#
exp(−iωt)dt.
(3.63)
Similar to the IR spectrum (3.56), these formulas can be transferred to the wavenumber-dependent representation, which is the one actually calculated in Travis. For fixed molecular orientation with respect to the laboratory coordinate system, this yields
Ik(ν˜)= h 8ε20kBT
(ν˜in−ν˜)4 ν˜
1−exp
−hck ν˜
BT
Z ∞
−∞
Dα˙xx(τ)α˙xx(τ +t)E
τ exp(−2πicνt˜ )dt (3.64) and
I⊥(ν˜)= h 8ε20kBT
(ν˜in−ν˜)4 ν˜
1−exp
−hck ν˜
BT
Z ∞
−∞
Dα˙xy(τ)α˙xy(τ +t)E
τexp(−2πicνt˜ )dt. (3.65) The formulas for random molecular orientation with respect to the laboratory coordinate sys-tem are given by
Ik(ν˜)= h 8ε02kBT
(ν˜in−ν˜)4 ν˜
1−exp
−hck ν˜
BT
45a(ν˜)+4γ(ν˜)
45 (3.66)
and
I⊥(ν˜)= h 8ε02kBT
(ν˜in−ν˜)4 ν˜
1−exp
−khcν˜
BT
3γ(ν˜)
45 (3.67)
with the isotropic contribution a(ν˜)=Z ∞
−∞
*α˙xx(τ)+α˙yy(τ)+α˙zz(τ) 3
α˙xx(τ +t)+α˙yy(τ+t)+α˙zz(τ+t) 3
+
τ
·exp(−2πicνt˜ )dt
(3.68)
and the anisotropic contribution γ(ν˜) =Z ∞
−∞
"
1 2
D α˙xx(τ)−α˙yy(τ) α˙xx(τ +t)−α˙yy(τ+t) E
τ
+1 2
D α˙yy(τ)−α˙zz(τ) α˙yy(τ +t)−α˙zz(τ +t) E
τ
+1 2
D α˙zz(τ)−α˙xx(τ) α˙zz(τ +t)−α˙xx(τ+t) E
τ
+3D
α˙xy(τ)α˙xy(τ+t)E
τ +3D
α˙yz(τ)α˙yz(τ +t)E
τ
+3D
α˙zx(τ)α˙zx(τ +t)E
τ
#
exp(−2πicνt˜ )dt.
(3.69)
The quotientI⊥(ν˜)/Ik(ν˜)provides the depolarization ratioρ(ν˜)as a function of the wavenum-ber. The result, however, is only meaningful in the region of Raman bands. Anywhere else, it is just numerical noise due to the division of two very small numbers.