Improved Theory of Difference Vibrational Spectroscopy and

Article pubs.acs.org/JCTC

Improved Theory of Difference Vibrational Spectroscopy and Application to Water Tatsuya Joutsuka† and Akihiro Morita*,†,‡ †

Department of Chemistry, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Kyoto 615-8520, Japan

‡

ABSTRACT: Measuring difference spectrum of two closely related systems is a versatile experimental method to highlight the difference in the two systems. However, the computation of a minuscule difference between two spectra by molecular dynamics (MD) simulation is far more challenging than that of each spectrum in terms of statistical convergence. Therefore, we have proposed a theory of difference spectra, which accelerates the calculation of difference spectra by several orders of magnitude [Sakaguchi, S. et al. J. Chem. Phys. 2014, 140, 144109]. The present paper reports our subsequent advances in the computational method to greatly improve its accuracy, numerical stability, and applicability. The present method of computation based on the nonequilibrium MD simulation allows for general molecular models including polarizable ones without sacrificing its computational efficiency. The improved method was applied to polarizable liquid water and yielded difference spectra of infrared, Raman, and vibrational sum frequency generation spectroscopies. The present method enables us to analyze various difference spectra of large molecular systems using MD simulation. has been established for the IR, Raman,11 and SFG spectroscopies.12 Therefore, one might apply the MD simulation of the vibrational spectra to the difference spectra as well, via subtraction of two calculated spectra, in principle. However, such a conventional subtraction method would be computationally too demanding and often infeasible, because the satisfactory convergence of tiny difference spectra would require prohibitively large statistical sampling to eliminate a common huge background. To overcome this difficulty, we have developed a general theory of difference spectra, which allows us to directly calculate the difference spectra themselves in background-free conditions.13,14 The proposed theory represents the difference spectra with two terms, namely the configuration term and the time evolution term. These terms were treated with first-order perturbation with the assumption that the small difference is treated as a perturbation. We have applied the theory to liquid argon and demonstrated that the new theory can reproduce the accurate difference spectra much more efficiently than the conventional subtraction by several orders of magnitude. We think that the present theory will become a promising route to analyze general difference spectra. In the present paper, we substantially improve the stability and applicability of the theory for the purposes of practical applications. The previous theory involved the stability matrix in the time evolution term, which is time-consuming to calculate and often numerically unstable.15,16 It also hinders the use of polarizable molecular models from the computational cost, since the computation of stability matrix would require to solve the coupled-perturbed equations of self-consistent

1. INTRODUCTION Vibrational spectroscopies, such as infrared (IR) absorption, Raman scattering,1 and sum frequency generation (SFG),2 provide powerful means to identify molecular species. These vibrational spectroscopies decompose observed spectra into various bands in a frequency domain. Each band signifies characteristic moiety of molecular structure since molecular vibrations absorb or scatter light by their own resonant frequency. Such selectivity to the molecular species renders the vibrational spectroscopies general probe methods of molecules in diverse areas of chemistry, biology, physics, etc. However, such decomposition into frequency domain is sometimes not sufficient to selectively probe some molecules by the vibrational spectroscopies. Such difficulties typically arise when we are interested in a partial system embedded in a large system, and the vibrational band from the partial system overlaps with that from the rest. If the vibrational bands overlap with each other in the frequency domain, the interpretation of vibrational spectra becomes ambiguous. To overcome these difficulties, difference spectroscopy is an effective and versatile method.3−6 Measuring the difference spectrum of two closely related systems lets us eliminate the unnecessary background and thus highlight some spectral bands of our interest. The difference spectra have been applied to various systems, such as simple solution,3 interface,4,5 and proteins.6 Interpretation of the vibrational spectra is greatly supported by advancement of theoretical analysis with molecular dynamics (MD) simulation.7−10 The MD simulation can derive a vibrational spectrum of a system from its microscopic structure and dynamics and thereby offer valuable insight into the observed/calculated spectral features. The computational method of vibrational spectra by the time correlation formulas © 2016 American Chemical Society

Received: July 13, 2016 Published: September 12, 2016 5026

DOI: 10.1021/acs.jctc.6b00697 J. Chem. Theory Comput. 2016, 12, 5026−5036

Article

Journal of Chemical Theory and Computation polarizations. To resolve the instability problem, first we tried to employ the multiple time step approach of time evolution17 to the stability matrix. The approach of multiple time step surely improved the numerical stability, though does not fully resolve the instability problem. Therefore, we further pursue an alternative route, and adopt the exact, nonperturbative expressions of the configuration and time evolution terms. The exact expressions are more accurate in principle than the first-order perturbation, while these methods require comparable computational costs. Furthermore, the nonperturbative expressions obviate the use of stability matrix and hence eliminate the practical problems associated with the stability matrix. The revised method therefore allows for polarizable models, which greatly expands our applications to various vibrational spectra. Using the polarizable model of water we developed, we apply the present method to the IR and Raman difference spectra of bulk water and to the SFG difference spectra of the air/water interface. The performance of the present calculation is demonstrated in comparison to the experimental difference spectra. The outline of the remainder of this paper is as follows. In Section 2, we present the improved computational scheme to calculate difference spectra. Section 3 discusses the computational accuracy and efficiency of the present method. The proposed computational method is applied to analyze the experimental IR and Raman spectra in bulk water and the SFG spectrum at air/water interface for temperature change in Section 4. Concluding remarks follow in Section 5.

with each other. This is a straightforward method to obtain the difference spectra, though it costs an enormous computational effort to deal with a huge background, as stated in Section 1. This numerical subtraction is called the “numerical method” hereafter, while our proposed method described in the following is called the “analytical method”. The difference in the time correlation function in eq 2 is rewritten using the phase space variables Γ as ΔC(t ) = ⟨B(0)B(t )⟩ − ⟨B(0)B(t )⟩0 =

(3)

Here, ρ and ρ0 are the thermal distribution functions for the perturbed (H) and reference (H0) states, respectively ρ[Γ] =

Z=

∫ d Γexp(−βH[Γ]),

Z0 =

⎡

⎤ ∂H ∂ ∂H ∂ ⎥ − , ∂qj ∂pj ⎥⎦ ⎣ ∂pj ∂qj

∑ ⎢⎢ j

(4)

∫ d Γexp(−βH0[Γ]) ⎡

iL0 =

⎤ ∂H0 ∂ ∂H0 ∂ ⎥ − ∂qj ∂pj ⎥⎦ ⎣ ∂pj ∂qj (5)

∑ ⎢⎢ j

where qj and pj are the coordinate and momentum of the jth degree of freedom. The expression of ΔC(t) in eq 3 indicates that the difference is attributed to two terms, i.e. the distribution functions (ρ and ρ0) and the time evolution operators (exp(iLt) and exp(iL0t)). Accordingly, ΔC(t) in eq 3 is decomposed into two terms ΔC(t ) = ΔCc(t ) + ΔCt(t )

(6)

where the first term is called the configuration term

(1)

where C(t) = ⟨B(0)B(t)⟩ denotes the arbitrary time correlation function, and the angle bracket means an ensemble average. The physical quantity B depends on the type of spectroscopy, such as dipole moment for IR spectra and polarizability for Raman spectra. C(t) in eq 1 could be a cross-correlation function ⟨B(0)B′(t)⟩ with B′ ≠ B, e.g. in the case of the sum frequency generation spectroscopy.18,19 We deal with the autocorrelation function in the following discussion of this section for simplicity, though the cross-correlation function can be treated in the same manner. Suppose that the reference state described by the Hamiltonian H0 is changed to the new state described by the Hamiltonian H = H0 + ΔH. Then the difference in the line shape functions is given as

ΔCc(t ) =

∫ d Γ(ρ[Γ] − ρ0 [Γ])B[Γ] exp(iL0t )B[Γ]

(7)

and the second term is the time evolution term ΔCt(t ) =

∫ d Γρ[Γ]B[Γ][exp(iLt ) − exp(iL0t )]B[Γ] (8)

The above decomposition eqs 6−8 provides the basis for the present theory of difference spectra. We note that the above formulas are based on the assumption that the spectral difference is attributed to the change in the Hamiltonian from H0 to H. The change in temperature is treated with some modification, as we detailed in ref 13. There are some cases that the physical quantity B is influenced by the perturbation parameter besides the phase space variables Γ, and such cases are treated in the Appendix. It is also straightforward to extend the above formulas to nonHamiltonian dynamics, as long as the equilibrium distribution functions of eq 4 and the Liouville operators of eq 5 are defined.20 2.2. Calculation Scheme. Here we present the improved calculation scheme of the two terms, ΔCc(t) in eq 7 and ΔCt(t)

∫ dt e−iωt ΔC(t ) ∫ dt e−iωt [⟨B(0)B(t )⟩ − ⟨B(0)B(t )⟩0 ]

exp( −βH0[Γ]) Z0

ρ0 [Γ] =

L and L0 in eq 3 are the Liouville operators for the perturbed and reference states

2. THEORY Here we briefly outline the theory of the difference spectra presented in our previous papers13,14 and further improvement in the computational method mentioned above. 2.1. Decomposition Scheme. Vibrational spectra are generally represented by the line shape function I(ω) as a function of frequency ω. It is written by the Fourier transformation of the time correlation function as11 1 I(ω) = ∫ dt e−iωt ⟨B(0)B(t )⟩ = 21π ∫ dt e−iωt C(t ) 2π

1 2π 1 = 2π

exp( −βH[Γ]) , Z

where β = 1/(kBT) is the inverse of the Boltzmann constant kB times temperature T, and Z and Z0 are the partition functions of the two states

iL =

ΔI(ω) =

∫ d Γρ[Γ]B[Γ]exp(iLt )B[Γ] − ∫ d Γρ0 [Γ]B[Γ]exp(iL0t )B[Γ]

(2)

where the subscript 0 means the ensemble average over the reference state (H0). Eq 2 could be evaluated in principle by conducting two different MD simulations and subtract them 5027


Article

Journal of Chemical Theory and Computation

perturbation and on the nonequilibrium MD simulation, will be discussed in Section 3.3. We further performed the above calculations of ΔCc(t) and ΔCt(t) by switching H and H0. Since the difference between H and H0 is mutual, switching the roles of H and H0 leads to equivalent results except that ΔCc(t) and ΔCt(t) change their signs. We call the original route as “forward” and the switched route as “backward”. The 2-fold calculations double the statistical sampling without substantial increase in the computational cost. This is due to the fact that eq 9 for ΔCc(t) involves the ensemble average at the reference state (H0) while eq 8 for ΔCt(t) at the sample state (H). The former ensemble average is taken during equilibrium MD trajectories at the reference state (H0), while the latter is taken during equilibrium MD at the sample state (H). Therefore, the set of MD trajectories by H and H0 can be commonly utilized in the forward and backward routes. The agreement between the two routes is a working criterion to confirm the validity of the present calculation. On the other hand, the deviation may occur by a large perturbation so that the phase-space regions sampled at the two states (H0 and H) do not overlap sufficiently. This is the same problem with the usual free energy calculations.24 In the present paper, we mainly adopted the analytical method without using the first-order perturbation, where the configuration term is calculated using eq 9 and the time evolution term using the nonequilibrium MD scheme of Figure 1. We call the present method the “nonperturbative” method, in contrast to the first-order perturbation formulas previously proposed.13

in eq 8, without resorting to the perturbation expansion. The configuration term ΔCc(t) in eq 7 is rewritten as ΔCc(t ) =

=

∫ dΓ

exp(− βH ) B[Γ] exp(iL0t )B[Γ] − ⟨B(0)B(t )⟩0 Z

∫ dΓ

exp(−βH0) Z0

exp(− β ΔH )B[Γ] exp(iL0t )B[Γ]

∫ dΓ =

exp(−βH0) Z0

exp(− β ΔH )

− ⟨B(0)B(t )⟩0

⟨exp[− β ΔH(0)]B(0)B(t )⟩0 − ⟨B(0)B(t )⟩0 ⟨exp[− β ΔH ]⟩0

(9)

This expression can be evaluated with MD simulation at the reference state (H0). If we expand eq 9 with respect to ΔH with the assumption that ΔH is treated as a perturbation, ΔCc(t) in eq 9 is expressed up to the first order of ΔH as ΔCc(t ) ≈ −β⟨(ΔH(0) − ⟨ΔH ⟩)B(0)B(t )⟩0 + O({ΔH }2 ) (10)

The above first-order expression eq 10 has been proposed and employed in our previous paper.13 We notice that computational cost of eq 9 as well as eq 10 is of the same order to that of the usual time correlation function, ⟨B(0)B(t)⟩0, using the equilibrium MD trajectories at the reference state (H0). On the other hand, the time evolution term ΔCt(t) in eq 8 was calculated via the first-order perturbation, which includes the stability matrix,15,16 in the previous scheme.13 The time evolution of the stability matrix tends to be numerically unstable, as we mentioned in Section 1. This problem becomes particularly serious when fast intramolecular vibrational motions have to be considered. Therefore, we adopt an alternative, numerically more stable way to calculate ΔCt(t) by the nonequilibrium MD simulation. Such nonequilibrium MD simulations have been employed in the calculations of the multidimensional vibrational spectra.21−23 The computational scheme is depicted in Figure 1. During time evolution of an equilibrium MD trajectory at the sample

3. ACCURACY AND EFFICIENCY This section describes the accuracy and efficiency of the present computational method by MD simulation of difference vibrational spectra in liquid water. Section 3.1 describes the MD conditions to evaluate the accuracy and efficiency. Sections 3.2 and 3.3 examine the computational accuracy and efficiency of the difference vibrational spectra of water in various types of perturbation. 3.1. Computational Conditions. The MD simulation of a model liquid water was conducted using 100 water molecules enclosed in a cubic cell of 14.41 Å with three-dimensional periodic boundary conditions. Such a small size is chosen here to allow for the calculation of numerical method, which requires huge statistical sampling. Once we confirm the reliability of the analytical method in comparison to the numerical one, a larger system will be used to investigate experimental difference spectra in Section 4. The water molecules are described with the charge response kernel (CRK) model,25 which incorporates intramolecular vibration and electronic polarization. The intermolecular interactions consist of the Lennard-Jones and electrostatic terms. The Lennard-Jones potential is considered between oxygen sites of water as

Figure 1. Computational scheme for the time evolution term ΔCt(t) in eq 8 by the nonequilibrium MD. The equilibrium trajectory at the sample state (H) is denoted with a solid black line, and the bifurcated trajectories by H0 are denoted with red dashed lines. Every bifurcated trajectory lasts for time Ttcf, during which the difference of the two time evolutions (exp(iLt)B − exp(iL0t)B) is sampled at t (0 < t < Ttcf). This nonequilibrium simulation with the bifurcation is repeated with the interval of Ttcf.

state (H), we bifurcate the trajectory at t = 0 by the two Liouville operators L and L0 and run the two trajectories from the same initial conditions during time Ttcf. Then we obtain the difference in the time evolutions, exp(iLt)B(0) − exp(iL0t) B(0), during the time t (0 < t < Ttcf), which is multiplied by B(0) to yield a statistical sample of ΔCt(t) in eq 8. This bifurcation procedure is repeated every time Ttcf for acquiring the statistical sampling. The performance of the above two calculation schemes of ΔCt(t), based on the first-order

U LJ =

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ ϵ 4 ∑ ⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r r ⎝ ij ⎠ ⎦ i>j ⎣⎝ ij ⎠

(11)

where rij is the distance between i-th and j-th oxygen sites, ϵ = 649.6 J/mol, and σ = 3.205 Å. The cutoff length of the Lennard-Jones interaction is set to rcut = 5.205 Å. The electrostatic interaction is described with fluctuating charges by the CRK model. The long-range electrostatic interactions are 5028


Article

Journal of Chemical Theory and Computation treated by the Ewald method,26 with the Ewald parameter set to 0.604 Å−1 and the reciprocal-space cutoff to 3.49 Å−1. The simulation was performed in the canonical (NVT) ensemble using the Nosé−Hoover thermostat.27,28 The temperature T was set to 300 K unless otherwise noted. The time evolution was treated by the multiple time step method, called RESPA (reversible reference system propagator algorithm),17 where the time steps are 2.0 and 0.1 fs for intermolecular and intramolecular interactions, respectively. In the calculation of the time evolution term, the equilibrium trajectory is bifurcated every 512 fs for the nonequilibrium simulation. All the calculations were performed by an in-house code. 1600 and 160−320 independent initial configurations are generated for the numerical and analytical methods, respectively, and equilibration runs are carried out for 200 ps for each trajectory. Then the sampling runs follow for a sufficiently long time to reach the statistical convergence, as discussed below. The difference vibrational spectra are examined with the following four types of perturbations: (1) Mass (isotope substitution): one H2O molecule is changed to HOD, (2) Dispersion potential: the Lennard-Jones parameter ϵ in eq 11 of one water molecule is increased by 10%, (3) Temperature: the temperature is increased from 300 to 303 K (by 1%), (4) Electrostatic potential: the equilibrium partial charges of one water molecule are increased by 1%. The parameters to be perturbed in four cases are summarized in Table 1. Note that one molecule among 100 in the MD cell

Figure 2 (0). The vibrational spectrum of neat H2O mainly consists of three vibrational bands; the high-frequency band around 3420 cm−1 arises from the OH stretch, the band around 1660 cm−1 arises from the water bend, and the low-frequency band arises from the intermolecular motions.29 The difference spectra for the four cases in Table 1 are shown in Figure 2 (1)−(4). In each panel, the black line shows the result of the nonperturbative analytical method, and the red line shows the result of the numerical method. The blue and green dashed lines are the configuration and time evolution terms, respectively. The forward and backward results are shown in black and red lines, respectively, in the inset of each panel. 3.2.1. Mass Change. Figure 2 (1) shows the calculated difference vibrational spectrum between neat H2O and one HOD in 99 H2O. The relative intensity of the difference spectrum is about 1/200 (≈0.005) of the reference spectrum, which is understood as one OH oscillator among 200 is replaced with OD. This isotope substitution is signified in the obvious negative amplitude in the OH stretch region (around 3400 cm−1) and the positive one in the OD stretch region (around 2490 cm−1). It also causes red shift in both bands of water bending and intermolecular vibration due to the increasing effective mass. The result by our analytical method agrees well with that by the numerical method. The decomposition of the analytical result shows that the time evolution term is dominant in this case, particularly in the intramolecular vibrations. Relative importance of the time evolution term over the configuration term is understandable in the mass change, since the mass change does not affect the equilibrium configuration while it directly affects the nuclear dynamics. The inset shows that the forward and backward calculations fairly well agree with each other, though the deviation is discernible. The deviation between the two routes implies that the isotope substitution may have substantial perturbation on the phase-space sampling, as it may have large influence on the local dynamics near the isotope atom. 3.2.2. Dispersion Potential Change. Figure 2 (2) shows the difference spectrum when the Lennard-Jones parameter ϵ of one water molecule is increased by 10%. The result of the analytical method can reproduce numerically exact spectrum, and the agreement between the forward and backward routes is quite satisfactory in the inset. In this case the configuration term is dominant in the whole frequency range, which is in contrast to the case of the mass change. The configuration term becomes important as the change in the intermolecular potential directly affects the equilibrium configuration. Regarding the band shape of the difference spectrum, it is noteworthy that the OH stretch band shows a blue shift. This indicates that increasing Lennard-Jones ϵ parameter results in a longer intermolecular distances and thus weaker hydrogen bonds. Since σ = 3.205 Å in eq 11 is longer than the average O−O distance of a hydrogen-bonding pair, the Lennard-Jones interaction in the first hydration shell is repulsive in liquid water. On the other hand, the water bending region shows a red shift, which is consistent with the fact that the stretch and bend frequencies tend to anticorrelate with each other.30 The difference spectrum in the libration band shows a positive peak around 400 cm−1 and a negative region in higherfrequency side. This feature means red shift in the libration band due to the weakened hydrogen-bond interactions.29

Table 1. Parameters To Be Changed in Four Cases of Perturbationa perturbation

(1) mass H/D

(2) dispersion ϵ [J/mol]

(3) temp T [K]

(4) electrostatic qH [au]

original perturbed

HOH HOD

649.6 714.6

300.0 303.0

0.3286 0.3319

a

Note that the parameter(s) of one molecule is changed in case (1), (2), or (4).

is changed in the case (1), (2), or (4). In case (4), the partial charge of the oxygen site of the molecule changes accordingly (qO = −0.6572 → −0.6638) to keep the charge neutrality, while the charge response kernel of the molecule is unchanged. The effects of these perturbations are examined through calculating the difference spectra of velocity autocorrelation function ΔI(ω) =

1 2π

∞

N

∫−∞ dt e−iωt N1 ∑ [⟨vj(0)·vj(t )⟩ − ⟨vj(0)·vj(t )⟩0 ] j=1

(12)

where j denotes the atomic sites of the liquid water, and vj denotes the velocity. We calculate the difference spectra in four cases by both the analytical and numerical methods and assess the computational accuracy and efficiency of the former. Such difference spectra of velocity autocorrelation functions are convenient to the discussion, since they are related to the vibrational density of states and require modest computational costs to allow for the numerical method. 3.2. Accuracy. In this section, we examine the accuracy of our analytical method in comparison to the numerical method. The reference vibrational spectrum of liquid water is shown in 5029


Article


Figure 2. (0) Reference vibrational spectrum (Fourier transform of velocity correlation function) of liquid water. (1)−(4) Difference spectra (eq 12) in cases (1)−(4) of Table 1, respectively. In each panel, the black line shows the nonperturbative analytical method, and the red line shows the numerical method. The result by the analytical method is decomposed into the configuration (blue dashed) and time evolution (green dashed) terms, and the inset of each panel displays the forward (black) and backward (red) calculations. All the intensities are normalized with the maximum intensity of the reference spectrum (0).

3.2.4. Electrostatic Potential Change. Figure 2 (4) shows the difference spectrum when the partial charges of one water molecule are increased by 1%. The OH stretch band shows a red shift, while the bend band shows a blue shift. All the bands in this panel (4) have the opposite signs to those of (2) Lennard-Jones ϵ or (3) temperature. Those opposite signs indicate that the increasing polarity strengthens the hydrogen bonds, in contrast to the increasing temperature or ϵ. We note in this case that the configuration and time evolution terms have comparable contributions. The forward and backward calculations also show an excellent agreement. We have shown in the above four cases that the present analytical method very well reproduces the results of numerical subtraction. Regarding the relative contributions of the configuration and time evolution terms, the four cases show rather different trends. In general, the former reflects the difference in the equilibrium structures between the two states, and the latter depends on the difference in the nuclear dynamics. In many cases including change in potential

3.2.3. Temperature Change. Figure 2 (3) shows the difference spectra when temperature is increased by 1% from 300 to 303 K. The analytical result reproduces the numerical difference spectrum accurately. The time evolution term is much smaller than the configuration term because the temperature change has a minor indirect influence on the time evolution of molecules through the thermostat. Overall, the line shape of difference spectrum is very similar to that by the Lennard-Jones ϵ change since both perturbations weaken the hydrogen-bond interactions. The OH stretch region shows a blue shift with increasing temperature, whereas the bend band shows a red shift, consistent to experimental results of IR and Raman spectroscopies.31−35 The agreement between forward and backward calculations is excellent in this case, apparently better than that in the mass change in (1), though the amplitude of the difference spectrum in panel (3) is much larger than that in (1). The excellent agreement suggests that the structure and dynamics of liquid water cover a similar region of phase space when the temperature slightly increases. 5030


Article


estimated convergence times by the numerical and analytical methods are summarized in four cases (1)−(4) in Table 2, which shows that the analytical method accelerates the statistical convergence by factors ranging from 54.9 to 5650.

functions, the configuration term tends to play a major role, though the time evolution term may have non-negligible contribution. 3.3. Computational Efficiency. In this section, we discuss the statistical convergence of calculated difference spectra by various methods. Section 3.3.1 compares the present analytical method to the numerical one, and Section 3.3.2 compares the present nonperturbative method to the first-order perturbation treatment. 3.3.1. Numerical and Analytical Methods. To compare computational efficiency of various methods, we first define the criterion of statistical convergence of the calculated difference spectra by estimating their error bars. Using the MD trajectories for a total of time Ttot, we divide them into n segments of the simulation time t = Ttot/n each, calculate ΔIk(ω) for each segmented sample k (= 1 ∼ n), and take the sample variance of the calculated spectra by ⎛ T tot ⎞ 1 S ⎜ω , t = ⎟= ⎝ n ⎠ n−1 2

Table 2. Convergence Times To Calculate the Difference Spectra in Cases (1)−(4) of Water by the Numerical and Analytical Nonperturbative Methodsa

a

k=1

(13)

where ΔI (ω) =

1 n

n

∑ ΔIk(ω) k=1

(1) mass

(2) dispersion

(3) temp

(4) electrostatic

8840 2.75 (3210)

12900 10.4 (1240)

80.8 1.47 (54.9)

29400 5.20 (5650)

The ratios of acceleration are shown in parentheses.

We see from this table that the numerical method yields a particularly short convergence time (80.8 ns) in case (3), because the perturbation is relatively large among the four cases. This is evidenced in the ordinates of the difference spectra in Figure 2, where the relative intensity in case (3) is plotted in the order of 0.01 while those in the other cases are ∼0.001. The analytical method becomes more advantageous when the perturbation is smaller. 3.3.2. Nonperturbative and First-Order Methods. Next, we discuss the efficiency of the nonperturbative analytical method in comparison to the method based on the first-order perturbation. The former is adopted in this paper, while the latter was proposed and employed in our previous paper.13 We evaluated the convergence times by the nonperturbative and first-order analytical methods in the same manner as in Section 3.3.1. We also analyzed the convergence behaviors of the configuration term and time evolution term separately and estimated the convergence time for each term by both methods. The comparison of the convergence times by the two analytical methods for the configuration term, time evolution term, and total (configuration + time evolution) are summarized in Table 3.

n

∑ (ΔIk(ω) − ΔI(ω))2

perturbation numerical/ns analytical/ns (ratio)

(14)

The statistical uncertainty in the calculated difference spectra is inferred from the Student confidence interval36 δ(ω, t) = tn−1S(ω, t)/√n, where tn−1 is the confidence coefficient at 95% confidence. The convergence criterion is given with a sufficiently small δ (ω, t) relative to the spectral amplitude of ΔI(ω) . Figure 3 displays examples of the convergence behavior in case (4) electrostatic potential change by various methods. The

Table 3. Convergence Times by the Nonperturbative and First-Order Analytical Methods for the Configuration Term (CT), the Time Evolution Term (TET), and the Total (CT +TET) (1) mass perturbation nonperturbative/ns first-order/ns

Figure 3. Log−log plot of the relative uncertainty of the spectral amplitude δ(ωmax,t)/|ΔI(ωmax)| to the simulation time t in case (4) electrostatic potential change. Red symbols denote the numerical method, black symbols denote nonperturbative analytical, and blue symbols denote the first-order analytical. The line of each color shows the linear regression with a slope of −1/2. The convergence criterion δ/|ΔI |= 0.1 is shown with the dashed line.

total

CT

(2) dispersion TET

total

CT

TET

2.75 187 2.74 6.24 126 2.88 (3) temperature

10.4 12.1 41.2 12.4 14.9 18.6 (4) electrostatic

perturbation

total

CT

TET

total

CT

TET

nonperturbative/ns first-order/ns

1.47 1.58

1.48 1.28

2.54 0.840

5.20 1.64

12.5 13.0

1.42 0.122

Table 3 indicate that the performance of convergence is more or less equivalent between the nonperturbative and firstorder analytical methods. Such equivalent performance is rather surprising, particularly for the time evolution term. During the MD sampling of the time evolution term B(0)(exp(iLt) − exp(iL0t))B(0), the origin t = 0 of the time correlation function is taken at the interval of Ttcf (= 512 fs) by the nonequilibrium MD scheme as shown in Figure 1, while the origin is taken every time step Δt (= 2.0 fs) in the scheme of first-order analytical method like the calculation of a usual time correlation

relative uncertainty δ(ω,t)/|ΔI(ω)| is monitored at the frequency of maximum absolute amplitude, ω = ωmax, in the difference spectrum ΔI(ω) (i.e. ωmax ≈ 3490 cm−1 in Figure 2 (4)). The logarithm of relative uncertainty, log10 [δ(ωmax,t)/ |ΔI(ωmax)|], is represented with a linear regression line with a slope of −1/2 in each method. We estimate the convergence time t from the regression line so that log10 [δ(ωmax,t)/ |ΔI(ωmax)|] = −1.0 (or δ(ωmax,t)/|ΔI(ωmax)| = 0.1). The 5031


Article


Figure 4. Samples of the time evolution term in eq 15 in case (2) dispersion change. The results by the nonperturbative formula are shown in black, and those by the first-order formula in red. The ordinate is scaled in the atomic unit. (a) The shaded area by each color represents the region that 1000 trajectories cover. An exemplary trajectory by each formula is shown by the thick line, and the standard deviation at t = 220 fs is shown by the error bar. (b) Averaged time evolution term among 1000 samples using the nonperturbative (black) and first-order (red) formulas.

function.26 Therefore, the first-order formula is expected to outperform the nonperturbative method in terms of the sampling efficiency, though the sampling efficiency is not necessarily inversely proportional to the sampling interval. In spite of the apparent advantage in the sampling efficiency, there is a problem in the first-order method to hinder the convergence. The first-order method includes successive multiplications of stability matrix to trace the time evolution, which is prone to be numerically unstable.15,16 The statistical convergence of the time evolution term is examined by investigating individual samples of the velocity autocorrelation function δct(t ) =

1 N

at air/water interface. The computations of difference spectra are carried out with the nonperturbative analytical method. 4.1. IR and Raman for Temperature Change. In this subsection, we calcualte the IR and Raman difference spectra in bulk water for the temperature change and compare the calculated results with experiments.31−35,37−39 IR and Raman Formulas. We first define the IR line shape function and the absorption spectrum. The line shape function IIR(ω) is given by eq 1 with B replaced by dipole moment of the measured system M I IR (ω) =

N

∫ dt e−iωt ⟨M(0)·M(t )⟩

(16)

The absorption coefficient α(ω) is written as

∑ vj·[exp(iLt ) − exp(iL0t )]vj j=1

1 2π

(15)

α(ω)n(ω) =

Note that the ensemble average of δct(t) at the sample state (H) yields the time evolution term ΔCt(t) = ⟨δct(t)⟩ in eq 8. Figure 4a displays 1000 samples of δct(t) in eq 15 (black) and its first-order approximation (red). This panel shows that the first-order results of 1000 trajectories are much more scattered than the corresponding exact results. The standard deviation at t = 220 fs (error bars) by the first-order method is 3.5 times larger than by the exact formula, and the difference in the standard deviation keeps increasing as time progresses. Figure 4b plots the average of the 1000 samples by the nonperturbative (black) and first-order (red) formulas in a longer time scale. The first-order result shows a divergent behavior even after the average of 1000 samples, whereas the nonperturbative result shows no such behavior. In practical calculations of the Fourier transformation of ΔCt(t), this divergent behavior is suppressed by multiplying a damping function in the long time limit. Nevertheless, such huge unphysical fluctuations in the first-order method have an effect on the convergence efficiency. Thus, the nonperturbative method is preferred to the first-order method in future applications since the former guarantees a better accuracy in principle and, at the same time, holds a similar speed-up to calculate difference spectra.

4π 2ω 2 IR I (ω) 3VckBT

(17)

where n(ω) is the refractive index, V is the volume of the system, and c is the speed of light. In this formula, we assume the harmonic quantum correction factor to replace the quantum time correlation function with the classical one.40 The expressions of the Raman spectra are given as follows. The line shape functions of polarized and depolarized Raman scattering are Raman Ipol (ω) =

1 2π

∫ dt e−iωt ⟨A̅ (0)A̅ (t )⟩

(18)

Raman Idepol (ω) =

1 2π

∫ dt e−iωt ⟨Tr[B(0)B(t )]⟩

(19)

where A and B are the isotropic and anisotropic parts of the polarizability tensor A, A = Tr[A]/3, and B = A − Tr[A]I. The differential cross sections for polarized and depolarized Raman scattering are d 2σpol

Raman ω − ω ⎤4 β ℏωIpol (ω) ⎡ = ⎢2πn 0 ⎥⎦ ⎣ dωd Ω c 1 − e−βℏω

(20)

d 2σdepol

Raman ω − ω ⎤4 β ℏωIdepol (ω) ⎡ = ⎢2πn 0 ⎥⎦ ⎣ dωd Ω c 1 − e−βℏω

4. ANALYSIS OF VIBRATIONAL SPECTRA OF WATER The present method is applied to difference spectra with temperature change in comparison to the experimental spectra. Section 4.1 deals with IR and Raman difference spectra in bulk water, and Section 4.2 deals with the SFG difference spectrum

(21)

where ω0 = 488 nm is assumed as the frequency of the incident radiation. To compare with the experimental results, we calculate the polarized (VV) and depolarized (VH) Raman spectra defined by11 5032


Article


Figure 5. (a) Calculated IR difference spectrum Δ[α(ω) n(ω)] of liquid water for the temperature change from 25 to 50 °C and (b) the corresponding experimental spectrum.34,41,42 (c) Calculated unpolarized Raman difference spectrum Δ[IVV(ω) + IVH(ω)] from 21 to 44 °C and (d) the corresponding experimental spectrum.38 The intensity of each difference spectrum is normalized with the maximum of the initial spectrum. 2 ⎛ d 2σ ⎞ d 2σpol 2 d σdepol + IVV = ⎜ ⎟ = dωd Ω 15 dωd Ω ⎝ dωd Ω ⎠

(22)

2 ⎛ d 2σ ⎞ 1 d σdepol IVH = ⎜ ⎟ = ⎝ dωd Ω ⎠⊥ 10 dωd Ω

(23)

experimental spectrum which is calibrated with the refractive indices n(ω).41,42 Results. Figure 5 compares the calculated difference spectra to the experimental ones. Panel (a) shows the calculated IR difference spectrum Δ[α(ω)n(ω)] in liquid water for the temperature change from 25 °C (298 K) to 50 °C (323 K), and (b) is the corresponding experimental IR difference spectrum.34,41,42 Panel (c) shows the calculated unpolarized Raman difference spectrum Δ[IVV(ω) + IVH(ω)] for the temperature change from 21 °C (294 K) to 44 °C (317 K), and (d) is the experimental one.38 Each difference spectrum is normalized so that the maximum intensity of the initial spectrum is unity. The calculated difference spectra (a, c) well reproduce the experimental ones (b, d). Those line shapes indicate the blue shift of the OH stretch band by the elevated temperature, which is consistent to the argument in Figure 2 (3) in Section 3.2.3. A similar blue shift was also observed in previous experimental IR34,37,43 and Raman38,39 spectra. 4.2. SFG at Air/Water Interface for Temperature Change. χ(2) Formula. Difference SFG spectra can be calculated in the same manner as described above. The second-order nonlinear susceptibility χ(2) consists of the vibrationally resonant term χR and the nonresonant term χNR, χ(2) = χR + χNR, and the former is represented by the following time correlation function between the polarizability A and the dipole moment M

respectively. MD Conditions. To analyze the experimental difference spectra of temperature change, we calculated the difference IR spectrum from T = 25 to 50 °C and Raman spectrum from T = 21 to 44 °C and compared them with corresponding experimental spectra.34,38 The MD calculations of such large temperature jumps are carried out using intermediate states to reduce the magnitude of the perturbation. The temperature jump in the IR difference spectrum is divided into 25 to 35 °C and 35 to 50 °C, and that in the Raman spectrum the jump is divided into 21 to 31 °C and 31 to 44 °C. The difference spectra from the initial and final states are obtained by summing the segmented difference spectra. In the present MD simulation we expand the box size to L = 24.66 Å containing 500 water molecules. The cutoff length of nonbonded interactions is set to 10.0 Å, and the Ewald coefficient is set to 0.373 Å−1. We prepare 1600 different initial configurations, and generate their MD trajectories in parallel. For each trajectory equilibration is carried out for 20 ps and the production follows for 100 ps. During the calculation of time evolution term, Ttcf in Figure 1 is set to 4096 fs. The other MD conditions are the same as in Section 3. The time correlation functions thus obtained are multiplied with a Gaussian smoothing filter, exp(−t2/τ2) with τ = 916 fs, to reduce computational noise. The calculated IR spectrum IIR(ω) is converted to α(ω)n(ω) in eq 17 and compared to the

R χpqr =

iω kBT

∫0

∞

dt exp(iωt )⟨A pq(t )Mr (0)⟩

(24)

in electronically nonresonant conditions,18,19 where ω is the IR frequency. The latter term χNR is regarded as a constant over the range of IR frequency. 5033


Article


6 are reasonable from related experimental and computational study on the temperature jump.44 Self-Correlation. Finally, we examine the difference spectrum of self-correlation approximation. As the polarizability A and dipole moment M in eq 24 are decomposed into molecular contributions, A = ∑kαk and M = ∑kμk, the time correlation function in eq 24 can be accordingly decomposed into the selfand cross correlations

The difference spectrum of eq 24 is calculated by treating the cross correlation function ⟨A(t)M(0)⟩, as mentioned in Section 2. We calculate difference SFG spectra of χ(2) of water surface. In the following discussion, we deal with the imaginary spectrum of χ(2), Im [χ(2)]. The difference spectra of Im[χ(2)] are determined by the resonant term, i.e. ΔIm [χ(2)] = ΔIm[χR], since χNR is a real quantity in electronically nonresonant conditions. MD Conditions. Computational details of difference SFG spectra follow the previous work of SFG calculation of air/ water interface.25 A slab of 500 water molecules is enclosed in an MD cell of Lx × Ly × Lz = 30 Å × 30 Å × 150 Å under the three-dimensional periodic boundary conditions. 1600 initial configurations are prepared and equilibrated for 80 ps before production for 720 ps. Thus, the total sampling time to calculate the SFG difference spectra is 1.15 μs. The computational procedure for the difference SFG spectrum follows that of difference IR and Raman spectra in Section 4.1. Here we calculate the ssp-polarized SFG difference spectrum for the temperature change from 300 to 303 K. ssp denotes the polarization combination of s-polarized SFG, s-polarized visible, p-polarized IR lights. Results of Difference SFG Spectrum. Figure 6 displays the calculated SFG difference spectrum ΔIm[χ(2)]. We see three

R χpqr =

=

iω kBT

∫0

∞

⎛ ⎞⎛ ⎞ k dt exp(iωt ) ⎜⎜∑ αpq (t )⎟⎟⎜⎜∑ μrl (0)⎟⎟ ⎝ k ⎠⎝ l ⎠

∞ iω k ∑ dt exp(iωt )⟨αpq (t )μrk (0)⟩ kBT k 0 ∞ iω k + ∑ dt exp(iωt )⟨αpq (t )μrl (0)⟩ kBT k ≠ l 0

∫

∫

self cross = χpqr + χpqr

(25)

The self-correlation χ , is considered as an approximation of χR in eq 25, by neglecting intermolecular correlations.45 The dashed line in Figure 6 shows the SFG difference spectrum calculated using the self-correlation, ΔIm[χself]. It fairly well reproduces the solid line in the higher-frequency region (3500 cm−1 ≲ ω), whereas the agreement becomes worse in the lower frequency region, where the dynamics of intermolecular correlation becomes more substantial.9,46,47 Since the selfcorrelation spectrum χself is easier to calculate than χR and becomes a useful approximation, it can be utilized to discuss the difference spectrum for a qualitative purpose. self

5. CONCLUSION Difference spectra are widely utilized in spectroscopic experiments to selectively detect a part of the observed system, though MD calculation of such tiny difference spectra has been quite challenging in terms of the computational cost. Therefore, we recently proposed a new general theory to calculate difference spectra in background-free conditions,13,14 which drastically facilitated the calculation of such tiny difference spectra. The present work significantly sophisticates the method of calculation by using the exact, nonperturbative formula of difference spectra and the nonequilibrium MD scheme. The refined treatment resolved the numerical problem of instability in calculating stability matrices and also allows for polarizable molecular simulations without sacrificing computational efficiency. The present theory with these refinements becomes a robust and practically useful means to analyze various experimental difference spectra by MD simulation. We demonstrated its performance by applying it to the IR and Raman difference spectra in liquid water and the SFG difference spectra at the air/water interface with the flexible and polarizable CRK model with considerable success. The proposed analytical method is generally applicable to various difference spectra as long as the spectra are described with time correlation functions. The present method is also generally applied to various methods of computing time correlation functions. For example, the use of quantum MD simulations such as the centroid48 and ring-polymer MD simulations49 is straightforward. In the calculation of the time correlation function, use of frequency mapping7 to obtain instantaneous value of B(t) is also possible. Based on the theory

Figure 6. ssp-polarized SFG difference spectra (ΔIm [χ(2)]) for the temperature change from 300 to 303 K (solid line). The dashed line shows the spectrum using the self-correlation (ΔIm[χself]). (Inset) original SFG spectrum (Im[χ(2)]). The amplitudes of the spectra are normalized by the maximum intensity of the SFG spectrum.

bands in this difference spectrum, (i) a positive one at ω ≲ 3350 cm−1, (ii) a negative one at 3350 cm−1 ≲ ω ≲ 3550 cm−1, and (iii) a positive one at 3550 cm−1 ≲ ω. These features are interpreted in reference to the original Im[χ(2)] spectrum in the inset. The original Im[χ(2)] spectrum shows a negative broad band at 3000−3600 cm−1 and a positive sharp band at about 3700 cm−1. The former is assigned to the hydrogen-bonded OH, and the latter is assigned to the free OH at the surface. Therefore, the two features (i, ii) of the difference spectrum indicate a blue shift of the former, negative Im[χ(2)] band. The positive feature (iii) of the difference spectrum with a dip at ω ≃ 3700 cm−1 indicates that the positive Im[χ(2)] band at 3700 cm−1 is broadened. We could not provide the corresponding SFG difference spectrum of this system from 300 to 303 K by the numerical method, as it demands too much computational cost. The analytical method we proposed enables us to calculate this tiny difference spectrum with a feasible cost of computation. We can argue that the spectral features in Figure 5034


Article


(2) Shen, Y. R. In Proceedings of International School of Physics “Enrico Fermi; Hansch, T. W., Inguscio, M., Eds.; North Holland: Amsterdam, 1994; Vol. CXX; pp 139−165. (3) Davis, J. G.; Gierszal, K. P.; Wang, P.; Ben-Amotz, D. Nature 2012, 491, 582−585. (4) Yamakata, A.; Soeta, E.; Ishiyama, T.; Osawa, M.; Morita, A. J. Am. Chem. Soc. 2013, 135, 15033−15039. (5) Wiebalck, S.; Kozuch, J.; Forbrig, E.; Tzschucke, C. C.; Jeuken, L. J. C.; Hildebrandt, P. J. Phys. Chem. B 2016, 120, 2249−2256. (6) Furutani, Y.; Shimizu, H.; Asai, Y.; Fukuda, T.; Oiki, S.; Kandori, H. J. Phys. Chem. Lett. 2012, 3, 3806−3810. (7) Skinner, J. L.; Auer, B. M.; Lin, Y.-S. Adv. Chem. Phys. 2009, 142, 59−103. (8) Torii, H. J. Phys. Chem. B 2007, 111, 5434−5444. (9) Buch, V.; Tarbuck, T.; Richmond, G. L.; Groenzin, H.; Li, I.; Shultz, M. J. J. Chem. Phys. 2007, 127, 204710. (10) Laage, D.; Stirnemann, G.; Sterpone, F.; Rey, R.; Hynes, J. T. Annu. Rev. Phys. Chem. 2011, 62, 395−416. (11) McQuarrie, D. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (12) Morita, A. J. Phys. Chem. B 2006, 110, 3158−3163. (13) Sakaguchi, S.; Ishiyama, T.; Morita, A. J. Chem. Phys. 2014, 140, 144109. (14) Sakaguchi, S.; Ishiyama, T.; Morita, A. J. Chem. Phys. 2014, 141, 149901. (15) Meyer, H. J. Chem. Phys. 1986, 84, 3147−3161. (16) Kay, K. G. J. Chem. Phys. 1994, 101, 2250−2260. (17) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990−2001. (18) Morita, A.; Hynes, J. T. J. Phys. Chem. B 2002, 106, 673−685. (19) Morita, A.; Ishiyama, T. Phys. Chem. Chem. Phys. 2008, 10, 5801−5816. (20) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. J. Chem. Phys. 2001, 115, 1678−1702. (21) Jansen, T. l. C.; Duppen, K.; Snijders, J. G. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 134206. (22) Yagasaki, T.; Saito, S. J. Chem. Phys. 2008, 128, 154521. (23) Hasegawa, T.; Tanimura, Y. J. Chem. Phys. 2006, 125, 074512. (24) Frenkel, D.; Smit, B. Understanding Molecular Simulation (Second ed.); Academic Press: San Diego, 2002. (25) Ishiyama, T.; Morita, A. J. Chem. Phys. 2009, 131, 244714. (26) Allen, M. P., Tildesley, D. J. Computer Simulation of Liquids; University Press: New York, 1989. (27) Nosé, S. J. Chem. Phys. 1984, 81, 511−519. (28) Hoover, W. G. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (29) Ashihara, S.; Huse, N.; Espagne, A.; Nibbering, E. T. J.; Elsaesser, T. J. Phys. Chem. A 2007, 111, 743−746. (30) Falk, M. Spectrochim. Acta, Part A 1984, 40, 43−48. (31) Ashihara, S.; Fujioka, S.; Shibuya, K. Chem. Phys. Lett. 2011, 502, 57−62. (32) Šašić, S.; Segtnan, V. H.; Ozaki, Y. J. Phys. Chem. A 2002, 106, 760−766. (33) Libnau, F. O.; Kvalheim, O. M.; Christy, A. A.; Toft, J. Vib. Spectrosc. 1994, 7, 243−254. (34) Venyaminov, S. Y.; Prendergast, F. G. Anal. Biochem. 1997, 248, 234−245. (35) Walrafen, G. E.; Hokmabadi, M. S.; Yang, W. H. J. Phys. Chem. 1988, 92, 2433−2438. (36) Ingram, J. Introductory Statistics; Cummings Publishing Company: Menlo Park, CA, 1974. (37) Larouche, P.; Max, J.-J.; Chapados, C. J. Chem. Phys. 2008, 129, 064503. (38) Walrafen, G. E.; Fisher, M. R.; Hokmabadi, M. S.; Yang, W. J. Chem. Phys. 1986, 85, 6970−6982. (39) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Cohen, R. C.; Geissler, P. L.; Saykally, R. J. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 14171− 14174. (40) Bader, J. S.; Berne, B. J. J. Chem. Phys. 1994, 100, 8359−8366.

and computational methods established above, we are now in the process of expanding the application of the analysis of difference spectra.

■

APPENDIX: MASS CHANGE In calculating a difference spectrum, the physical quantity B in the time correlation function may depend on the perturbation parameter besides the phase space variables. A typical example is seen in the velocity correlation function perturbed by mass change, where B is the velocity of a particle, v. Since the velocity v = p/m depends on the mass m besides the momentum p, the mass change affects the B itself. In such a case, the physical quantity changes from B0 = v0 = p/m0 to B = v = p/m with the mass change from m0 to m. In this section, we extend the theory of difference spectra in such cases. Accordingly, the difference in the time correlation function in eq 2 is rewritten as ΔC(t ) = ⟨B(0)B(t )⟩ − ⟨B0 (0)B0 (t )⟩0 =

∫ d Γρ[Γ]B[Γ] exp(iLt )B[Γ] − ∫ d Γρ0 [Γ]B0 [Γ] exp(iL0t )B0 [Γ] (26)

where ρ, ρ0, L, and L0 are given in eq 4 and eq 5. ΔC(t) in eq 26 is also decomposed into the configuration term ΔCc(t) and the time evolution term ΔCt(t) in eq 6, where ΔCc(t ) =

∫ d Γ(ρ[Γ] − ρ0 [Γ])B0[Γ] exp(iL0t )B0[Γ] (27)

ΔCt(t ) =

∫ d Γρ[Γ]{B[Γ] exp(iLt )B[Γ] − B0[Γ] exp(iL0t )B0[Γ]} (28)

The calculation of eq 27 is analogous to that of eq 9, as B in eq 9 is replaced with B0 in eq 27 during the MD simulation at the reference state (H0). The calculation of eq 28 is also straightforward using the nonequilibrium MD depicted in Figure 1. When the MD trajectory at the sample state (H) is bifurcated, B is switched to B0 in the MD at the reference state (H0) and the difference in the time correlation, B exp(iLt)B − B0 exp(iL0t)B0, is sampled. The auxiliary term in the first-order perturbation in ref 13 is not necessary in the nonequilibrium MD scheme.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Grants-in-Aid for Scientific Research (Nos. JP25104003, JP26288003) by the Japan Society for the Promotion of Science (JSPS) and Ministry of Education, Culture, Sports and Technology (MEXT), Japan. The computations were carried out using the supercomputers in Research Center for Computational Science, Okazaki, Japan.

■

REFERENCES

(1) Advances in Infrared and Raman Spectroscopy; Clark, R. J. H., Hester, R. E., Eds.; Wiley: Heyden, 1975−1985; Vol. 1−12. 5035


Article

Journal of Chemical Theory and Computation (41) Bertie, J. E.; Lan, Z. Appl. Spectrosc. 1996, 50, 1047−1057. (42) Pinkley, L. W.; Sethna, P. P.; Williams, D. J. Opt. Soc. Am. 1977, 67, 494−499. (43) Lock, A. J.; Bakker, H. J. J. Chem. Phys. 2002, 117, 1708−1713. (44) Nagata, Y.; Hasegawa, T.; Backus, E. H. G.; Usui, K.; Yoshimune, S.; Ohto, T.; Bonn, M. Phys. Chem. Chem. Phys. 2015, 17, 23559−23564. (45) Ishiyama, T.; Morita, A. J. Phys. Chem. A 2007, 111, 9277−9285. (46) Imoto, S.; Xantheas, S. S.; Saito, S. J. Chem. Phys. 2013, 138, 054506. (47) Auer, B. M.; Skinner, J. L. J. Phys. Chem. B 2009, 113, 4125− 4130. (48) Voth, G. A. Adv. Chem. Phys. 1996, 93, 135−218. (49) Habershon, S.; Manolopoulos, D. E.; Markland, T. E., III; Miller, T. F. Annu. Rev. Phys. Chem. 2013, 64, 387−413.

5036


Improved Theory of Difference Vibrational Spectroscopy and

Recommend Documents