Notes on simulating two-dimensional Raman and terahertz-Raman signals with a full molecular dynamics simulation approach.

Recent developments in two-dimensional (2D) THz-Raman and 2D Raman spectroscopies have created the possibility for quantitatively investigating the role of many dynamic and structural aspects of the molecular system. We explain the significant points for properly simulating 2D vibrational spectroscopic studies of intermolecular modes using the full molecular dynamics approach, in particular, regarding the system size, the treatment of the thermostat, and inclusion of an Ewald summation for the induced polarizability. Moreover, using the simulation results for water employing various polarization functions, we elucidate the roles of permanent and induced optical properties in determining the 2D profiles of the signal.

INTRODUCTION

Intermolecular vibrations of molecular liquids and biological material in the frequency range 0–700 cm−1 play an essential role in many chemical and biological processes, because they promote reactive dynamics via interactions through intramolecular modes and because they are active at room temperature. While one-dimensional (1D) IR, THz, and Raman spectroscopic approaches are versatile tools for investigating collective intermolecular motion, using such approaches alone, it is not clear whether the vibrational modes that they investigate are mutually coupled, and whether the width of the vibrational mode peaks that they measure are from an inhomogeneous origin or a homogeneous origin, because these contributions are usually broadened and overlap in the 1D spectra. To elucidate these points, in 1993, the fifth-order two-dimensional (2D) Raman spectroscopy approach was proposed to distinguish the contributions from inhomogeneous broadenings using three sets of Raman pulses. This initiated the development of 2D infrared (IR) spectroscopy for intramolecular vibrational modes. While the possibility for applying 2D Raman spectroscopy to the study of anharmonicity, mode-mode coupling mechanisms, and dephasing processes of intermolecular modes has been explored theoretically, due to the unforeseen cascading effect of light emissions, experimental signals have been obtained only for CS2, benzene, and formamide liquids. However, recently developed single-beam spectrally controlled 2D Raman spectroscopy method overcomes the difficulty of cascading effects and creates a new possibility for measuring intra-molecular interactions of liquids by means of 2D Raman spectroscopy.

While the 2D Raman signals of liquid water has not yet been observed, 2D THz-Raman (or 2D Raman-THz) spectroscopy was introduced in 2012 to measure the water signal. In this method, the cascading effects are suppressed using two THz pulses and one set of Raman pulses. Because 2D THz-Raman utilizes the dipole moment, in addition to the polarizability, the applicability of this spectroscopy is different from that of 2D Raman spectroscopy. Moreover, the information obtained from 2D Raman and 2D THz-Raman spectroscopies can be used in a complementary manner to investigate the fundamental nature of intermolecular interactions.

Theoretically, the linear absorption (1D IR or 1D THz spectroscopy) signal and the third-order Raman spectroscopy (1D Raman spectroscopy) signal are obtained from the linear response functions of the optical observables. There, the main contribution to a signal arises from harmonic vibrational motion as a function of molecular polarizability or the dipole moment. Contrastingly, in 2D Raman and 2D THz-Raman spectroscopic approaches, the three-body nonlinear response function of the molecular polarizability and/or the dipole moment is measured to monitor molecular motion. Because the complex profiles of such 2D signals depend on many dynamic and structural aspects of the molecular system, full molecular dynamics (MD) simulations for the nonlinear response function play important roles in the design of 2D spectroscopy experiments and the analysis of their results, particularly in regard to intermolecular vibrations. In the 2D case, the anharmonicity of potentials and the nonlinearities of the polarizability contribute significantly to the signals. For these reasons, a reliable potential model and an accurate polarizability function are necessary to accurately predict quantitative features of 2D signals. Moreover, the theoretical details of the MD simulations, in particular, regarding the system size, the treatment of the thermostat, and the proper inclusion of an Ewald summation for the induced polarizability, must be carefully examined, because the 2D profile of a signal is also very sensitive to numerical errors. The purpose of this paper is to elucidate the significant points regarding the calculation of 2D Raman and THz-Raman signals using full MD approach while minimizing computational costs.

This paper is organized as follows. In Sec. II, we explain the methodology for simulating 2D Raman and 2D THz-Raman signals with full MD simulations. In Sec. III, we investigate the size dependence of simulations, the artifacts of the thermostat, sensitivity of the signals to the choice of the force field, and the importance of employing an Ewald sum of the induced polarizability function in calculating 2D Raman signals. In Sec. IV, we explain the role of the force field and the sensitivity of the signals to the choice of the polarizability functions, particularly in calculations of 2D THz-Raman signals for liquid water. Section V is devoted to concluding remarks.

FULL MD APPROACH FOR 2D RAMAN AND 2D THz-RAMAN SPECTROSCOPY

The optical observable in 2D spectroscopy is expressed in terms of the three-body response function as where , and can be the dipole moment, , or the polarizability, , of the molecules expressed as functions of the molecular position coordinates, . The classical mechanical expression for the response functions can be obtained by replacing the commutator and operators with the Poisson bracket and c-number observables as The response function in the classical limit is then expressed as

Equilibrium MD approach

Three full MD approaches have been developed to this time for evaluating the above response function. The first approach is based on equilibrium MD simulations. We evaluate the outer Poison bracket in terms of the time derivative of the observable, employing the relation, , with the molecular Hamiltonian , as follows: Here, β is the inverse temperature divided by the Boltzmann constant, and is the time derivative of the observable , defined as . Then, we obtain the response function in terms of the stability matrix as We can calculate this response function using the time and coordinate derivatives of or evaluated from the molecular trajectories and that are obtained from the equilibrium MD simulations. While the equilibrium MD approach is convenient for analyzing the 2D profile of a signal as the contribution of the anharmonicity and the nonlinear polarizability, this requires a great deal of CPU time and memory due to the computational intensiveness of the treatment of the equation of motion used to calculate the Jacobian (or stability block) matrix, . Moreover, because of the stability matrix, the convergence of the signal becomes very slow, particularly for a large molecular system.

Non-equilibrium finite field MD approach

The second approach that we consider, the non-equilibrium MD (NEMD) approach, was developed to evaluate double Poisson brackets from non-equilibrium trajectories in the case that there exist multiple external perturbations. In this approach, the outer Poison bracket in Eq. (3) is evaluated from the relation , where is the expectation value corresponding to calculated from the trajectories subjected to weak perturbations , with the external field acting on . Then, the Poisson bracket is also evaluated from , where is defined in the same manner as . This approach does not involve a convergence problem of the stability matrix. However, it is too computationally intensive, because we have to repeat the calculations four times in order to account for the configurations of the external perturbations.

Equilibrium-non-equilibrium hybrid MD approach

The third approach that we consider, the equilibrium-non-equilibrium hybrid approach, was developed to take advantage of the merits of both the equilibrium and non-equilibrium approaches. In the hybrid approach, the outer Poison bracket is evaluated in terms of the time derivative of the observable, while the inter-Poison bracket is evaluated using the NEMD approach. Then, the response function is expressed as In this approach, we first obtain the time derivative of the dipole moment or polarizability from the equilibrium trajectories at time t = −t1. Next, we evaluate the dipole moments and or the polarizabilities and at time t = t2 from the non-equilibrium trajectories, which are generated by a perturbation at time t = 0, either or , resulting from the external electric field of the jth pulse acting on the dipole moment or the polarizability .

SIMULATING 2D RAMAN SIGNALS

Single-beam spectrally controlled 2D Raman spectroscopy, which was developed recently, allows us to obtain clear 2D spectra of molecular liquids in the THz region by suppressing cascading effects. This motivated us to carry out MD simulations of 2D Raman spectroscopy, in addition to 2D THz-Raman spectroscopy. Here, we discuss the important points in carrying out full MD simulations of 2D Raman spectroscopy, specifically with regard to the conditions of the simulations and the verification of force fields and polarizability functions.

Unless otherwise noted, the MD simulations were carried out as follows. The 2D Raman signals were obtained using the hybrid MD approach described in Sec. II C. Such a signal is given by where Δt is the time step used in integrating the equations of motion and is the polarizability with the equilibrium trajectories at time −t1. For comparison, we also calculated the 1D Raman spectra using the expression In the cases of formamide and carbon disulfide (CS2), we introduced a harmonic quantum correction factor and evaluated the signal from . The MD simulations were carried out with 108 molecules for the CS2 liquid and 64 molecules for the other liquids. The simulations were carried out in a cubic box with periodic boundary conditions. The interaction potentials were cut off smoothly at a distance equal to half the box length using a switching function, and the long-range Coulomb interactions were calculated with the Ewald summation. The intramolecular geometries were kept rigid throughout the simulations, using the RATTLE algorithm. The equations of motion were integrated using the velocity-Verlet algorithm with Δt = 2.5 fs for water and Δt = 5.0 fs for the other molecular liquids. The system volume and total energy were fixed (in accordance with a microcanonical simulation) after the completion of the isothermal simulations that were carried out in equilibration with a Nose-Hoover chain thermostat. The volume was chosen to reproduce experimental densities: 0.997 g/cm3 for water, 1.120 g/cm3 for formamide, 0.815 g/cm3 for formaldehyde, and 1.270 g/cm3 for CS2. The temperature was set to 300 K for water and formamide, 255 K for formaldehyde, and 270 K for the CS2. To estimate the polarizability of each liquid, we employed the dipole-induced-dipole (DID) model with an Ewald summation using the permanent molecular polarizability, which was determined from the Huiszoon polarizability for water, the atomic polarizability for formamide and formaldehyde, and experimental data for CS2. In the hybrid MD approach, the NEMD part of the calculation was carried out with the Raman laser fields E = 4.0 V/Å for water, E = 1.0 V/Å for formamide and formaldehyde, and E = 2.0 V/Å for CS2. These signals were obtained by averaging over 106 initial configurations.

Size dependence of the simulations

Because our objective is to describe fast intermolecular modes, which arise from short-range intermolecular interactions, it is not necessary to carry out large-scale simulations with many molecules. In general, 64 molecules are sufficient to obtain a reliable signal. To illustrate this point, we employed 2D Raman signals for water calculated with various system sizes. In these computations, the interactions between the molecules were modeled using the TIP4P/2005 potential, and the full-order (see Appendix A 1) DID models were employed to evaluate the polarizability using an Ewald summation (see Sec. III D). The computational results are presented in Fig. 1 with (a) 32, (b) 64, (c) 108, and (d) 216 water molecules. No clear size dependence is observed for water molecules. Although there is a possibility that the size dependence varies among the types of molecules, the present results suggest that a system size of 32 molecules is sufficient to elucidate the qualitative properties of 2D signals, at least for liquid water.

FIG. 1.

2D Raman signals of the zzzzzz tensor elements for water with (a) 32, (b) 64, (c) 108, and (d) 216 molecules. The red and blue shadings represent positive and negative signals, respectively. The signal intensities are normalized with respect to the absolute value of the peak signal intensities.

Microcanonical (NVE) and canonical (NVT) simulations

In MD simulations, a thermostat is often employed to maintain the system temperature (canonical simulation). Because 2D spectroscopy is very sensitive to molecular motion, however, the thermostat may alter the 2D profiles of the signal. To demonstrate this point, we calculated and compared these 2D Raman signals for CS2 liquid by using the Nose-Hoover chain thermostat, i.e., a canonical ensemble (NVT), and compared these results with those obtained from microcanonical (NVE) simulations. These simulations were carried out with the Lennard-Jones (LJ) model using the polarizability described by the full-order DID model (see Appendix A 1). To evaluate the induced polarizability, we employed the Ewald summation. Figures 2(a)–2(c) illustrate the effects of the thermostat in the 1D and 2D Raman spectra. The broadened peaks near 50 cm−1 in Fig. 2(a) reflect the presence of intermolecular vibrations. While the 1D Raman signals in Fig. 2(a) are similar in both cases, there is a difference along the t2 axis in the thermostatted case for the 2D Raman signals in Figs. 2(b) and (c). The elongation of the negative signal along the t2 direction arises from the anharmonicity of the potential. This indicates that the Nose-Hoover chain thermostat acts as an undesirable source of anharmonicity for the molecular dynamics. Because the Nose-Hoover chain thermostat was formulated to study thermal equilibrium states, the dynamics obtained using this thermostat are not necessarily accurate. The sensitivity of the 2D measurements may reveal such inaccuracy. Thus, it is most prudent to simulate the 2D spectrum using the NVE ensemble without a thermostat.

FIG. 2.

The 1D and 2D Raman signals for CS2 calculated with NVE and NVT simulations. We also tested a different force field in the case of the NVE simulations. Panels (b)-(d) represent the zzzz tensor elements obtained using (b) the LJ model in the NVE, (c) the LJ model in the NVT simulations, and (d) the LJ + Coulomb model in the NVE simulations. The red and blue shadings represent positive and negative signals, respectively. Panel (a) displays 1D Raman signals for the cases (b)-(d). The signal intensities are normalized with respect to the absolute value of the peak signal intensities.

Choosing the force field

As was shown in a formamide case, the profiles of 2D Raman signals are very sensitive to the choice of the force field. Here, we demonstrate this point using two types of force fields for CS2 while keeping the molecular polarizability function fixed. The first model consists of Lennard-Jones interactions only (LJ model), in order to simulate intermolecular vibrational spectra, whereas the second model also includes Coulomb interactions (LJ + Coulomb model), in order to account for the structural information obtained from neutron and X-ray scattering experiments. While almost all full MD simulations for 2D Raman spectroscopy of CS2 have been carried out using the first model, here we examine the validity of results obtained using a more reliable potential that includes the Coulomb interactions. In both cases, full-order DID models were employed to evaluate the polarizability using the Ewald summation.

While the LJ and LJ + Coulomb results are similar in the 1D case considered in Fig. 2(a), we observe a difference in the 2D Raman case considered in Figs. 2(b) and 2(d). Because the anharmonicity of the potential is assumed to be large in the LJ + Coulomb case, the negative signal along the t2 axis does not decay even at 1 ps, while the elongation vanishes within 500 fs in the LJ case. It should be noted that the LJ + Coulomb potential is more reliable than the LJ potential, because the LJ + Coulomb potential reproduces the structure of the CS2 liquid more accurately. The present results indicate that an accurate force field is necessary in order to obtain the correct profile of a 2D Raman signal even in the case of a nonpolar molecular system, due to the sensitivity of the nonlinear response function. For this reason, the accuracy of the potential should be evaluated using the experimental results.

Ewald summation of the induced polarizability: Long-range effects

By nature, the 2D Raman profile is very sensitive to the functional form of the polarizability. For this reason, the signal profile is significantly affected by the cutoff of the polarizability function that is introduced to reduce the computational intensiveness of the simulation. To demonstrate this point, we calculated the 2D Raman signals for water, formamide, and formaldehyde described by the TIP4P/2005 potential, the modified T potential, and the 4-site model, respectively, with and without the Ewald summation for the induced polarizability. The DID model was employed to evaluate the total polarizability in all of the simulations. In the induced polarizability, we include the dipole-dipole interaction term defined by Eq. (A4), which includes the contributions from the infinite periodic images while taking into account the changes in the distances between all of the molecular pairs. Without using the Ewald summation, however, the long-range effect is ignored.

Figure 3 displays the computational results of the 2D Raman signals for (i) water, (ii) formamide, and (iii) formaldehyde. The 2D Raman signals displayed in Fig. 3(a) were obtained using the Ewald summation, while the 2D Raman signals displayed in Fig. 3(b) were calculated by cutting off the dipole-dipole interaction at half of the box length without using the Ewald summation. All of the simulations whose results are presented here were carried out with 64 molecules. Note that we examined the size dependence of the results by considering systems with 64, 108, and 216 molecules without using the Ewald summation. We found that there were only minor differences among the signals obtained from these small systems. While the 2D Raman signals were calculated from the same trajectories, they are significantly different for formamide and formaldehyde. By contrast, the 2D Raman signals for water exhibit only a slight difference near t1 = t2 = 0 fs, with the negative signals being weaker when computed using the Ewald summation. The cause of these differences is analyzed in Appendix B.

FIG. 3.

The 2D Raman signals of the zzzzzz tensor elements calculated (a) with and (b) without the Ewald summation of the induced polarizability for (i) water, (ii) formamide, and (iii) formaldehyde.

The results discussed above indicate that, although the effects depend on the type of molecule, the long-range contributions to the polarizability must be computed from the dipole-dipole interactions using the Ewald summation in the 2D Raman case.

SIMULATING 2D THZ-RAMAN SIGNALS

Using the hybrid approach discussed in Sec. II C, the response functions for the 2D Raman-THz-THz (RTT), THz-Raman-THz (TRT), and THz-THz-Raman (TTR) spectroscopic approaches are calculated as respectively. Here, in addition to the 1D Raman spectrum defined by Eq. (8), we also evaluated 1D THz spectrum expressed as where the pre-factor is a harmonic quantum correction factor that must be applied to the classical calculations. It should be noted that the types of information that we can obtain from the 2D Raman signal and each of the three THz-Raman signals are different, due to the role of the nonlinear polarizability. We can separate the inhomogeneous and anharmonic contributions to the signal more clearly in the cases of 2D THz-Raman spectroscopy than in the case of 2D Raman spectroscopy, because the inhomogeneous contribution arises from the TRT pulse configuration, while the anharmonic contribution arises from the RTT configuration.

The conditions of the 2D Raman simulation discussed in Sec. III, regarding the system size, usage of a thermostat, and the sensitivity to the force field, can be applied to the present 2D THz-Raman case as well. However, the 2D THz-Raman signals are insensitive to the use of the Ewald summation for the induced polarizability, because the 2D THz-Raman response functions account for only one Raman process that is sensitive to the induced polarizability contribution. In this section, we explain the significant points involved in the simulation of 2D THz-Raman signals specifically for the case of water, focusing on the polarizability and potential model, because presently, water is of primary interest in experimental investigations employing 2D THz-Raman spectroscopy.

Unless otherwise noted, the MD simulations considered here were carried out for water under the same conditions as in the 2D Raman case. The NEMD part of the calculation involved in the TRT response was carried out with a Raman pulse intensity of E = 4.0 V/Å, whereas those involved in the RTT and TTR responses were carried out with a THz pulse intensity of E = 0.1 V/Å. The 2D signals were obtained by averaging over 106 initial configurations.

Choosing the permanent polarizability

Here, we first demonstrate the sensitivity of the 2D THz-Raman signals to the choice of the permanent polarizability. We compared the 2D profiles obtained using the Huiszoon permanent polarizability and the coupled-cluster single- and double-excitation (CCSD) permanent polarizability (see Appendix A 4), while the TIP4P/2005 potential was used for the force field in both cases. The total dipole moment and polarizability were calculated using the DID model with the Ewald summation.

The parameter values for the Huiszoon polarizability were set to α = 1.626 Å3, α= 1.495 Å3, and α = 1.286 Å3, while those for the CCSD polarizability were set to α = 1.442 Å3, α = 1.375 Å3, and α = 1.321 Å3. Here, the X axis is defined as that connecting the hydrogen atoms, the Y axis lies along the bisector of the H–O–H angle, and the Z axis is perpendicular to the XY plane.

The computational results for the 2D RTT, 2D TRT, and 2D TTR signals with the Huiszoon permanent polarizability and the CCSD permanent polarizability are displayed in Fig. 4. As explained in Appendix C, the differences among the 2D profiles in the cases of the Huiszoon polarizability and the CCSD polarizability arise from the librational and translational modes, which are activated by the permanent and induced optical properties, respectively. We found that the anisotropy of the permanent polarizability in the Huiszoon case is larger than that in the CCSD case. For this reason, in the Huiszoon case, the signal intensity from the permanent optical properties was stronger, whereas in the CCSD case, the signal intensity from the induced optical properties was stronger.

FIG. 4.

The 2D profiles of the zzzz tensor elements calculated with (a) the Huiszoon polarizability and (b) the CCSD polarizability for the (i) 2D RTT, (ii) 2D TRT, and (iii) 2D TTR signals of water.

Choosing the polarizability functions

Next, we elucidate the role of the polarizability functions in 2D THz-Raman signals. While the inter-molecular interaction potentials were modeled by the TIP4P/2005 potential, we employed the full-order DID, atomic site dipole-induced-dipole (ASDID), and charge-flow dipole-induced-dipole (CFDID) polarizability function models to calculate the total dipole moment and polarizability with the Ewald summation. The definitions of these polarizability functions are presented in Appendix A. We display the computational results for 2D THz-Raman signals in Fig. 5. In the 2D RTT results, the first negative peak, close to t1 = t2 = 0, which arises from the librational motion, is prominent in the DID case, while it cannot be observed in the ASDID and CFDID cases. The 2D profiles of the RTT signals are similar, whereas those of the 2D TRT and TTR signals for these three cases differ significantly. This is due to the fact that the 2D TRT and TTR signals are more sensitive to the polarizability function, because the main contribution to the signal comes from the nonlinear element of the polarizability function. More specifically, the sensitivities of the 2D signals are due to the induced polarizabilities, as illustrated in Appendix D.

FIG. 5.

The zzzz tensor elements of (i) 2D RTT, (ii) 2D TRT, and (iii) 2D TTR signals for water calculated using the (a) DID, (b) ASDID, and (c) CFDID models, respectively.

Finally, we examine the sensitivity of the 2D THz-Raman signals to the choice of the force field by comparing the results obtained using the TIP4P/2005 (Ref. 61) and POLI2VS (Ref. 67) potentials. The TIP4P/2005 model is a simple point-charge model that can properly simulate several macroscopic thermodynamic properties of water, while the POLI2VS model is a polarizable water model developed for vibrational spectroscopies that can simulate a wide range of vibrational modes, from low-frequency intermolecular modes to high-frequency intramolecular modes. The POLI2VS model is a flexible potential model. However, in our simulations, we fixed the intramolecular geometries by setting the O–H bond length to 0.977 Å and the H–O–H bend angle to 105.14° in order to accelerate the simulations. To elucidate the sensitivity to the choice of the force field alone, we employed the CFDID polarizability function developed for the flexible POLI2VS model to calculate the total dipole moment and the polarizability, as explained in Appendix A 3 and Ref. 67 in the cases of the TIP4P/2005 and rigid-POLI2VS potentials, respectively.

The 2D RTT, 2D TRT, and 2D TTR signals calculated using the TIP4P/2005 model and the POLI2VS model are presented in Fig. 6. It should be noted that the CFDID polarizability function for the TIP4P/2005 and the POLI2VS model exhibits similar behavior. Thus, we conclude that the difference between the 2D THz-Raman signals in the two cases arises from the dynamical aspects of the potential models. Because the largest contribution to the signal is from the nonlinear polarizability rather than the anharmonic motion of the molecules in the 2D TRT and TTR cases, these 2D profiles are similar. By contrast, because the main contribution to the 2D RTT signals is from the anharmonicity motion of the molecules, these signals differ significantly.

FIG. 6.

The zzzz tensor elements of the (i) 2D RTT, (ii) 2D TRT, and (iii) 2D TTR signals calculated using (a) the TIP4P/2005 and (b) the POLI2VS potential models.

CONCLUSION

In this paper, we elucidated the important points involved in full MD simulations of 2D Raman and THz-Raman spectroscopic approaches. We found that the 2D signals obtained in such approaches are sensitive to the nonlinearity of the polarizability and the anharmonicity of the force fields. For the purpose of obtaining accurate 2D profiles of the signals, the proper choice of the polarizability functions is more important than that of the force field, particularly in the 2D Raman, 2D TRT, and 2D TTR cases. Although we used a small system size, we found that even in this case, in order to obtain accurate signals, we had to employ the Ewald summation for both the force field and the induced polarizability function. Although a rigorous force field model is necessary for obtaining quantitatively accurate results, the use of such a model does not guarantee that we will obtain the correct 2D profile unless we also use the right polarizability function.

It should be noted that conventional MD methods employ both a force field and polarizability functions in a rather empirical manner. For this reason, the verification of the calculated results is not easy. Although computationally intensive, the ab initio MD approach may be more useful for calculating 2D spectra, because it utilizes not only a force field but also a dipole moment and a polarizability based on the electronic state of the molecules. Even in this approach, however, many assumptions and approximations are involved. Thus, verification must be carried out by comparing the 2D profiles obtained with simulations and experiments.

In this paper, we restricted our investigation to 2D Raman and 2D THz-Raman spectroscopy, but most of the points discussed in this paper also apply to the full MD simulation of 2D IR spectroscopy. For the purpose of studying high-frequency intra-molecular modes, in addition to inter-molecular modes, which is a primary application of 2D IR spectroscopy, a quantum treatment of the atomic motion within molecules is also important. In any case, simulating multidimensional vibrational spectroscopy of molecular liquids is a stringent test to verify the accuracy of MD simulations.