PyRAMD Scheme: A Protocol for Computing the Infrared Spectra of Polyatomic Molecules Using Ab Initio Molecular Dynamics

1. Introduction

Molecular dynamics (MD) is the method of describing nuclear motions using classical mechanics [1,2,3,4]. If the potential energy surface (PES) for the simulation is taken from the quantum-chemical calculations, such simulations are called ab initio MD (AIMD) [4]. AIMD is a powerful tool for simulating various physical observables of medium to large systems (from tens to a few hundreds of atoms) with non-trivial dynamical behavior [3,4,5,6]. Examples of modeled observables can be radial distribution functions in diffraction experiments [7,8,9,10], dynamical properties of solutions [6,11], and various types of spectra [5,12,13]. Calculating of the spectral response in the case of electronic excitations, such as ultraviolet or X-ray spectroscopy, usually relies on sampling the conformational space of molecular systems [12,14]. In contrast, in the case of the vibrational spectra, i.e., infrared (IR) [4,5,10], Raman [13], vibrational circular dichroism [15], etc., the fluctuation-dissipation theorem (FDT) [16,17] has to be applied.

The bare AIMD simulations, by definition, model the motion of classical particles, thus producing a classical microcanonical ($NVE$) ensemble [18]. To emulate the canonical ensemble ($NVT$ or $NPT$), an artificial external system (thermostat and/or barostat) is added to the MD simulations [18]. As the real nuclei are quantum objects, the absent nuclear quantum effects (NQEs) can also be accounted for in the MD by various approaches [19], such as path-integral MD [20,21], quantum thermostats [22], Wigner sampling [14,23,24], etc. [19]

In this work, we provide a protocol for simulating the IR spectra of molecules with BOMoND software. This MD interface software can take energies and gradients from the quantum-chemical package ORCA 5 [25] or the semi-empirical package xTB [26]. The BOMoND is part of the PyRAMD packet for AIMD-based simulations, which is also capable of metadynamics [27] and MD-based mass-spectra calculations [28]. We will demonstrate examples of applications of BOMoND software, illustrating the various aspects of the MD procedure and subsequent analysis.

2. Methods

The results presented here were obtained using the BOMoND software from the PyRAMD package [29]. This software can be obtained from the corresponding repository (Ref. [29]). The density functional theory (DFT) calculations[30] at the BLYP-D3(BJ)/6-31G,[31,32,33,34] PBE-D3(BJ)/6-31G,[33,34,35] and PBEh-3c[36] and semi-empirical calculations at the GFN2-xTB levels of theory [37], including structure optimization and harmonic frequencies computation, were done using ORCA 5 [25] and xTB software [26], respectively. The BOMoND used the same program suits to get the gradients for the AIMD simulations. The dipole moment autocorrelation functions calculations and the subsequent Fourier transform (FT) [38] and regularized least-squares spectral analysis (rLSSA) [39] were done using the script from the PyRAMD package.

3. Simulation Protocol for the Vibrational Spectra Using the Molecular Dynamics

3.1. General Idea of Calculation of the IR Spectra from MD Trajectories

The MD simulation produces an MD trajectory, $N_{\mathrm{steps}}$ simulations points, which span the dynamics in time t from $t=0$ to $N_{\mathrm{steps}}·\Delta t$, where $\Delta t$ is the time step of the simulation. Therefore, for various computed physical observables ($\mathcal{O}$, e.g., coordinates, velocities, energy, dipole moment, polarizability, etc.), we can tell their values ($\mathcal{O}(n·\Delta t)={\mathcal{O}}_n$) at a n-th given point of simulation ($0\le n<N_{\mathrm{steps}}$), corresponding to the time $t=n·\Delta t$ [18].

As said in the introduction, the vibrational spectra from MD are computed with the help of the FDT [16,17]. We first need to compute the autocorrelation function ($\langle \mathcal{O}\left(0\right)·\mathcal{O}\left(t\right)\rangle$) of the variable $\mathcal{O}$ that is responsible for the response to the external excitation. For the IR spectra, observable is the total dipole moment of the molecule; for the Raman spectra, it is the molecular polarizability tensor; for the power spectrum, these are simply the mass-weighted coordinates, etc. [13]. For a real-valued observable, the autocorrelation function is defined as [4,5,13,40]

$$\langle \mathcal{O}\left(0\right)·\mathcal{O}\left(t\right)\rangle ={\int }_0^{\infty }\mathcal{O}\left(t^{\prime }\right)\mathcal{O}(t^{\prime }+t)d{t}^{\prime },$$

(1)

and it shows the repeatable periods in the time evolution of the given observable. The numerical calculation of the autocorrelation function for finite time series is available in various packages. For instance, in Python, one could use either the “correlate” method from the NumPy module [41] or the “signal.correlate” method from the SciPy module [42]. By performing the FT of the $\langle \mathcal{O}\left(0\right)·\mathcal{O}\left(t\right)\rangle$ as [4,5,13,40]

$$S\left(\nu \right)={\nu }^2{\int }_0^{\infty }\langle \mathcal{O}\left(0\right)·\mathcal{O}\left(t\right)\rangle ·exp(-i2\pi \nu )dt$$

(2)

we can transfer the response characteristic from the time domain (t) to the frequency domain ($\nu$), thus obtaining the spectrum of a given process. For numerical reasons, a more advantageous computational strategy is to use the velocity of the observable ($\dot{\mathcal{O}}$) to obtain the spectrum [5]. In this case, the expression to be used is[5]

$$S\left(\nu \right)={\int }_0^{\infty }\langle \dot{\mathcal{O}}\left(0\right)·\dot{\mathcal{O}}\left(t\right)\rangle ·exp(-i2\pi \nu )dt.$$

(3)

The resulting spectra from Equations 2 or 3 can be multiplied with the correction factors that account for the statistical properties at the different spectral ranges [5,43].

The FT procedure on the MD results (Equation 2 or 3) is usually done using a fast-FT (FFT) approach [38,44], which provides certain limitations on the frequency discretization of the resulting spectra. The frequency grid increment ($\Delta \nu$) is determined by the total trajectory duration ${\tau }_{\mathrm{tot}}$ through the Nyquist–Shannon–Kotelnikov (NSK) theorem [45,46,47] as $\Delta \nu =1/{\tau }_{\mathrm{tot}}$, while the time step $\Delta t$ determines the upper boundary of the spectrum as ${\nu }_{max}=1/\Delta t$. To increase the formal resolution, extending the trajectory for a given observable by an arbitrary number of steps with zeros is possible, which is called zero-padding [48]. However, such direct extension will cause artifacts in the resulting spectrum; thus, the zero-padding should be combined with either frequency filters or smoothing techniques [39,48].

An alternative procedure to the FFT could be the regularized least-squares spectral analysis (rLSSA). The idea is to transform the time-dependent dataset

$$\{(t_1,y_1),(t_2,y_2),\dots ,(t_N,y_N)\}$$

(4)

composed of N time values $t_k$ ($1\le k\le N$) with corresponding observable values $y_k$ into a frequency domain with M points

$$\{({\nu }_1,f_1),({\nu }_2,f_2),\dots ,({\nu }_M,f_M)\},$$

(5)

where ${\nu }_l$ ($1\le l\le M$) is the frequency and the $f_l$ is the corresponding spectral amplitude value. The rLSSA procedure is derived from the regularized weighted least-squares analysis (rwLSSA) procedure by choosing a unit weights matrix (see Refs. [39,49]). As a result, the rLSSA is just a matrix transformation[39,49]

$$\mathbf{f}=\left|{\Sigma }_{\alpha }{\mathcal{S}}^{†}\mathbf{y}\right|,$$

(6)

where $\mathbf{y}=(y_1,\dots ,y_N)$ is the N-dimensional vector of the time-domain values, $\mathbf{f}=(f_1,\dots ,f_M)$ is the M-dimensional vector of the frequency-domain spectrum amplitude values, $\mathcal{S}$ is the $N\times M$ matrix of elements $S_{kl}=exp(-\mathfrak{i}2\pi {\nu }_l{t}_k)$ with $\mathfrak{i}$ being imaginary unit, and ${\mathcal{S}}^{†}$ is the conjugate transpose of matrix $\mathcal{S}$. The matrix ${\Sigma }_{\alpha }$ of size $M\times M$ is defined as ${\Sigma }_{\alpha }=\alpha \mathcal{E}+{\mathcal{S}}^{†}\mathcal{S}$ with $\mathcal{E}$ being a unit matrix of the size $M\times M$ and the regularization parameter $\alpha$ provided by the regularization criterion[39,49]

$$\alpha =tr\left({\mathcal{S}}^{†}\mathcal{S}\right)·{\left(M+{\displaystyle \frac{tr\left({\mathcal{S}}^{†}\mathcal{S}\right)·\left({\mathbf{y}}^{†}\mathbf{y}\right)}{M}}\right)}^{-1}.$$

(7)

To smooth the data even further, the parameter $\alpha$ can be increased compared to that from Equation 7.

To illustrate that, we can take the vibrational spectrum of the carbon dioxide (${\mathrm{CO}}_2$) shown in Figure 1. A single $NVT$-MD trajectory at the GFN2-xTB level of theory was obtained for a temperature of 300 K using the Berendsen thermostat and Maxwell-Boltzmann sampling of initial conditions. The trajectory was 1 ps in length with a time step of 0.5 fs, and half of the trajectory was taken as the equilibration phase. The NSK theorem provides the frequency resolution from the dataset of total length ${\tau }_{\mathrm{tot}}=0.5$ ps to be $\Delta \nu ={\left(c·{\tau }_{\mathrm{tot}}\right)}^{-1}=66.8$ cm⁻¹, where $c=0.02998$ cm/ps is the speed of light, therefore upon applying the FFT procedure without zero-padding with frequency filtering we obtain quite coarse spectra. By applying the rLSSA approach to the same dataset with the requested frequency increment of $\Delta \nu =1$ cm⁻¹, from the Equation 7 we get $\alpha =1000$ and the more detailed spectrum (see Figure 1). Note that the rLSSA procedure cannot provide new spectroscopic details; it just provides a smoother representation of the same dataset.

3.2. Cheapening Simulations by Using Large Integration Steps with Frequency Correction

The usual recommendation for computing the molecules’ vibrational spectra is to use as small time step $\Delta t$ as possible. The reason for that is the artificial blue shift from the numerical integration error, which is especially prominent for the high-frequency vibrational bands [50,51]. However, it is possible to correct this behavior by using an analytical formula that compensates for such a shift. In the case of the Verlet algorithm and its equivalent methods (velocity Verlet and leapfrog) [18], one needs to replace the observed frequencies from the FT ($\nu$) with the corrected ones (${\nu }_{\mathrm{corr}}$) using the equation [40,52]

$${\nu }_{\mathrm{corr}}={\displaystyle \frac{\sqrt{2·\left(1-cos(2\pi ·\Delta t·\nu )\right)}}{2\pi ·\Delta t}}=\nu ·sinc(\pi ·\Delta t·\nu ),$$

(8)

where $sinc\left(x\right)=sin\left(x\right)/x$. This correction is applicable if the NSK-like limit $\Delta t\le {\left(\pi {\nu }_{max}\right)}^{-1}$ is satisfied, where ${\nu }_{max}$ is the maximal vibrational frequency in the system. In the case of the higher-order methods, similar correction schemes can be derived; see an example in the Appendix A.

To illustrate this artificial blue shift of vibrational bands and the effectiveness of the correction in Equation 8, we ran four sets of ten $NVE$-MD trajectories of methane (${\mathrm{CH}}_4$) at GFN2-xTB level of theory with initial conditions sampled from the Maxwell-Boltzmann distribution at 300 K, each trajectory was 0.5 ps in length. The sets were differing by the time step $\Delta t$, namely with $\Delta t=0.1,0.5,1.0,1.5$ fs. The mean vibrational spectra for each set of trajectories and their frequency-corrected counterparts are shown in Figure 2. As one can see, the larger time steps lead to a noticeable shift of the C–H stretching band. In the reference spectrum with $\Delta t=0.1$ fs such peak is located at approximately 3015 cm⁻¹, whilst at simulations with $\Delta t=1.5$ fs this peak appears at approximately 3150 cm⁻¹. However, rescaling the frequency axis using Equation 8 restores the positions of the vibrational bands to their correct origin.

3.3. Simplified Wigner Sampling for Generating Initial Conditions

To include the NQEs in the MD simulations, we can use the approach of simplified Wigner sampling (SWS) introduced in Ref. [14]. The sampling relies on generating simultaneously a displacement along the coordinate and momentum axis in a given direction. The displacements are generated using a Gaussian distribution with standard deviations[14]

$${\sigma }_x^2={\displaystyle \frac{\hslash {\tau }_{\mathrm{SWS}}}{2m}}$$

(9)

for coordinate and

$${\sigma }_T^2=m{k}_{\mathrm{B}}{T}_{\mathrm{eff}}$$

(10)

for momentum. Here, $\hslash =1.05\times {10}^{-34}$ J$·$s is the Planck constant, m is the mass of the nucleus, $k_{\mathrm{B}}=1.38\times {10}^{-23}$ J/K is the Boltzmann constant, ${\tau }_{\mathrm{SWS}}$ is a free parameter with the dimension of time, and $T_{\mathrm{eff}}$ is the effective temperature defined as

$$T_{\mathrm{eff}}(T,{\tau }_{\mathrm{SWS}})=T+{\displaystyle \frac{\hslash }{2k_{\mathrm{B}}{\tau }_{\mathrm{SWS}}}},$$

(11)

where $T\ge 0$ is the temperature of the system.

The main problem of the SWS is the choice of the free parameter ${\tau }_{\mathrm{SWS}}$. In work [14], it was proposed to perform several MD simulations with different values of ${\tau }_{\mathrm{SWS}}$ and then take the parameter value that minimizes the total energy or the mean temperature of the simulation. In some sense, this approach relates to the variational principle, adjusting the underlying Wigner distribution to fit the PES as well as possible. However, such an approach is quite slow and inefficient. Thus, a cheaper criterion for the choice of ${\tau }_{\mathrm{SWS}}$ is required.

To provide an alternative criterion, we may use the harmonic approximation as the guiding principle, as single harmonic frequency calculation with analytical Hessian generally takes less computational time than running multiple MD trajectories. The total kinetic energy of the system is defined as [18]

$$\mathrm{KE}=\sum _{k=1}^{N_{\mathrm{at}}}{\displaystyle \frac{{\mathbf{p}}_k^2}{2m_k}},$$

(12)

where ${\mathbf{p}}_k$ is the momentum vector of k-th atom, $m_k$ is the mass of the k-th atom, and $N_{\mathrm{at}}$ is the total number of atoms in the system. In the case of the distribution sampled at $T=0$ with SWS, the mean value of ${\mathbf{p}}_k^2$ is given as $\langle {\mathbf{p}}_k^2\rangle =3{\sigma }_0^2=3\hslash m_k/\left(2\tau \right)$ (see Equation 10). This gives the total mean kinetic energy of $\langle \mathrm{KE}\rangle =3N_{\mathrm{at}}\hslash m_k/\left(4\tau \right)$, while for a given degree of freedom, the mean kinetic energy is $\langle {\mathrm{KE}}_1\rangle =\langle \mathrm{KE}\rangle /\left(3N_{\mathrm{at}}\right)=\hslash /\left(4\tau \right)$. In the harmonic approximation, for a given l-th mode ($l=1,2,\dots ,N_{\mathrm{f}}$, where $N_{\mathrm{f}}$ is the number of vibrational degrees of freedom) with vibrational frequency ${\nu }_l$, the kinetic energy is ${\mathrm{KE}}_l=h{\nu }_l/4$ [53], and the mean kinetic energy over all of vibrational modes is

$${\mathrm{KE}}_{\mathrm{h}}=\sum _{l=1}^{N_{\mathrm{f}}}{\displaystyle \frac{h{\nu }_l}{4}}={\displaystyle \frac{N_{\mathrm{f}}h\langle \nu \rangle }{4}},$$

(13)

where $\langle \nu \rangle =\left({\sum }_{l=1}^{N_{\mathrm{f}}}{\nu }_l\right)/N_{\mathrm{f}}$ is the mean vibrational frequency of the system. By taking $\langle {\mathrm{KE}}_1\rangle ={\mathrm{KE}}_{\mathrm{h}}/N_{\mathrm{f}}$, we arrive at the following expression

$${\tau }_{\mathrm{SWS}},\mathrm{h}={\displaystyle \frac{1}{2\pi \langle \nu \rangle }}={\left(2\pi \stackrel{\langle \nu \rangle }{\overbrace{{\displaystyle \frac{1}{N_{\mathrm{f}}}}\sum _{l=1}^{N_{\mathrm{f}}}{\nu }_l}}\right)}^{-1}.$$

(14)

Equation 14 provides a criterion for the choice of $\tau$ from the harmonic vibrational frequencies of the system. To show its applicability, we ran a set of $NVE$-MD calculations for molecules XH ($\mathrm{X}=\mathrm{F},\mathrm{Cl},\mathrm{Br},\mathrm{I}$), XH₂ ($\mathrm{X}=\mathrm{O},\mathrm{S},\mathrm{Se},\mathrm{Te}$), XH₃ ($\mathrm{X}=\mathrm{N},\mathrm{P},\mathrm{As},\mathrm{Sb}$), and XH₄ ($\mathrm{X}=\mathrm{C},\mathrm{Si},\mathrm{Ge},\mathrm{Sn}$) at the GFN2-xTB level of theory. The MD trajectories were obtained with 1 fs time step for 0.5 ps in total. The ${\tau }_{\mathrm{SWS}}$ for SWS initial conditions was scanned from 1 to 10 fs with 1 fs step, and the MD simulation for a given ${\tau }_{\mathrm{SWS}}$ contained ten individual trajectories. The optimal ${\tau }_{\mathrm{SWS}}$ was chosen as the one minimizing the mean total temperature of the MD simulation. The results are given in Figure 3. As one can see, both criteria correlate with each other (Pearson correlation coefficient is 0.71). Therefore, we can apply the Equation Figure 3 for other molecular systems.

3.4. Thermostats Incorporating Simplified Wigner Sampling

SWS approach, despite being quite crude, has its merits. Unlike the standard Wigner sampling, which can be used only at the beginning of the MD simulation, as the displacements are connected to the equilibrium reference structure, the SWS allows a re-sampling of the displacements during the simulation. This property makes the SWS compatible with canonical ensemble simulations, particularly with Andersen [54] and Berendsen [55] thermostats. Here, we will describe how to do it for both of these thermostats.

The canonical Andersen thermostat works by reassigning the velocity (or of one of the three components) of a randomly chosen atom according to the Maxwell-Boltzmann distributions at random points of MD simulation. The probability of reassignment is given as $p_{\mathrm{A}}=\Delta t/{\tau }_{\mathrm{A}}$, where ${\tau }_{\mathrm{A}}$ is the time of the reassignment. At each point of time, a random number $0\le p_{\mathrm{trial}}<1$ is generated from the uniform distribution, and if condition $p_{\mathrm{trial}}<p_{\mathrm{A}}$ is fulfilled, then the velocity reassignment procedure is initiated. Typically, ${\tau }_{\mathrm{A}}$ should be greater than the characteristic periods of motions in the system, so the Andersen thermostat does not significantly disrupt the dynamics of the molecule. To combine the Andersen thermostat with the SWS, one has to replace the Maxwell-Boltzmann distribution resampling stage with the SWS routine, sampling both the momentum along given degrees of freedom and the displacement of the coordinate.

The Berendsen thermostat relies on the soft resampling of the velocities at each step of the MD simulation by multiplying them with a scale factor [18,55]

$$s=1+{\displaystyle \frac{\Delta t}{{\tau }_{\mathrm{B}}}}·\left({\displaystyle \frac{T_{\mathrm{d}}}{T\left(t\right)}}-1\right),$$

(15)

where ${\tau }_{\mathrm{B}}$ is the relaxation time, a free parameter of the Berendsen thermostat, $T_{\mathrm{d}}$ is the desired temperature of the MD simulation and $T\left(t\right)$ is the instant temperature of the nuclei in the MD simulations, defined as $T\left(t\right)=2·\mathrm{KE}\left(t\right)/N_{\mathrm{f}}$ with $\mathrm{KE}$ being the total kinetic energy of the molecule at time t, and $N_{\mathrm{f}}$ is the number of degrees of freedom (i.e., $N_{\mathrm{f}}=3\times N_{\mathrm{at}}-N_{\mathrm{constr}}$, where $N_{\mathrm{at}}$ is the number of atoms in the system and $N_{\mathrm{constr}}$ is number of constraints). By replacing the $T_{\mathrm{d}}$ with the effective desired temperature according to Equation 11, we also make the Berendsen thermostat compatible with the SWS sampling.

3.5. Scaling of the vibrational spectra from the molecular dynamics

Scaling the harmonic vibrational frequencies and zero-point vibrational energies is a simple yet powerful tool for increasing the accuracy of IR spectra and thermodynamic properties calculations [56]. In the case of vibrational frequencies, the multiplication by a precomputed scale factor for a given method accounts for the systematic errors due to the quantum-chemical approximation quality and the absence of the anharmonic effects [56,57]. Therefore, it seems like a reasonable suggestion to also allow for scaling of the anharmonic spectra from the MD simulations.

Here, we propose the procedure for obtaining the scale factors ($\gamma$) for IR spectra from the MD and also provide a set of precomputed values for three quantum-chemical approximations: BLYP-D3(BJ)/6-31G, PBE-D3(BJ)/6-31G, and PBEh-3c. For this, we have taken six molecules: water (${\mathrm{H}}_2\mathrm{O}$), ammonia (${\mathrm{NH}}_3$), methane (${\mathrm{CH}}_4$), ethane (${\mathrm{C}}_2{\mathrm{H}}_6$), methylamine (${\mathrm{CH}}_3{\mathrm{NH}}_2$), and methanol (${\mathrm{CH}}_3\mathrm{OH}$). For each of these molecules at each of the quantum-chemical methods, three trajectories were computed using the combination of Andersen and Berendsen thermostats with ${\tau }_{\mathrm{A}}=100$ fs and ${\tau }_{\mathrm{B}}=50$ fs and the SWS sampling with the ${\tau }_{\mathrm{SWS}}$ taken according to criterion from Equation 14. Each trajectory was at $T=0$ K with the time step of $\Delta t=1$ fs and the total duration of 1 ps. The spectra were computed using the FFT procedure without ignoring the equilibration phase and with the application of the frequency correction from Equation 8. The experimental gas-phase IR spectra of the same molecular species were taken from the NIST Chemistry WebBook [58].

It is not straightforward to compare theoretical spectra with a given frequency increment with their experimental counterparts since they have different frequency increments. Moreover, when the frequency axis is scaled for the spectra from the MD, the frequency increment will also change. Therefore, the following procedure was applied to determine the scale factor of the MD spectra for a given molecule. Let us assume that the MD spectrum is given as a set of M points $\{({\nu }_1^{\left(\mathrm{MD}\right)},I_1^{\left(\mathrm{MD}\right)}),({\nu }_2^{\left(\mathrm{MD}\right)},$$I_2^{\left(\mathrm{MD}\right)}),\dots ({\nu }_M^{\left(\mathrm{MD}\right)},I_M^{\left(\mathrm{MD}\right)})\}$ and the experimental spectrum is a set of N points $\{({\nu }_1^{\left(\exp \right)},I_1^{\left(\exp \right)}),({\nu }_2^{\left(\exp \right)},$$I_2^{\left(\exp \right)}),\dots ({\nu }_N^{\left(\exp \right)},I_N^{\left(\exp \right)})\}$, where ${\nu }_i^{\left(X\right)}$ and $I_i^{\left(X\right)}$ are the i-th frequency and intensity obtained with method “X”. To remove the rotational substructure of the bands for small molecules, such as water or methane, the experimental data were also smoothed by a convolution with the Gaussian function. The frequency scale factor $\gamma$ is obtained by maximizing the following expression

$$C\left(\gamma \right)={\left({\mathbf{I}}^{\left(\mathrm{MD}\right)}\right)}^{†}\mathcal{D}\left(\gamma \right){\mathbf{I}}^{\left(\exp \right)}\to max,$$

(16)

where ${\mathbf{I}}^{\left(\mathrm{MD}\right)}$ and ${\mathbf{I}}^{\left(\exp \right)}$ are the vectors composed of intensity values of N and M dimensions each, and the $\mathcal{D}\left(\gamma \right)$ is the $M\times N$ matrix composed of elements

$$D_{ij}\left(\gamma \right)=exp\left(-{\displaystyle \frac{{\left(\gamma ·{\nu }_i^{\left(\mathrm{MD}\right)}-{\nu }_j^{\left(\exp \right)}\right)}^2}{2{\sigma }^2}}\right)$$

(17)

with $1\le i\le M$, $1\le j\le N$, and $\sigma =min(\gamma ·\Delta {\nu }^{\left(\mathrm{MD}\right)},\Delta {\nu }^{\left(\exp \right)})$ with $\Delta \nu$ being the frequency increments in the respective datasets. The scale factor $\gamma$ was searched for in the interval $0.8\le \gamma \le 1.2$.

Figure 4 shows an example of the result of the procedure described above. In this Figure, the unscaled and scaled spectra of methane obtained from MD simulations at the PBEh-3c level of theory are shown, as well as the corresponding raw and smoothed experimental data from the NIST Chemistry Webbook. As one can see, the scaling indeed improves the positions of the vibrational bands in the spectra. The scale factors obtained from all the training set molecules were averaged to produce the estimated scale factor for a given level of theory. The results are given in Table 1. In the cases of the BLYP and PBE-based spectra, the ammonia had to be removed from the dataset as the most intense low-frequency vibrational band, corresponding to the valence band vibrations, was too far off from the experimental position.

4. Discussion

In the previous section (Section 3), the protocol for the IR spectra simulations from the AIMD using the PyRAMD was introduced. In an orderly fashion, it looks as follows.

• First, we need to optimize the structure of the molecule at the given level of theory and compute harmonic vibrational frequencies. Then, using Equation 14, we can calculate the ${\tau }_{\mathrm{SWS}}$ parameter to define the SWS sampling routine.
• Then, we can set the Berendsen and Andersen thermostats for simultaneous usage in the $NVT$-MD simulation. The combination of the two acts as a friction and random force in more sophisticated thermostats, such as the Langevin-based models [59] (including the color-noise generalized Langevin equation [60]) and the Bussi-Donadio-Parrinello thermostat [61]. This requires setting the two free parameters: relaxation time ${\tau }_{\mathrm{B}}$ for the Berendsen thermostat and the resampling time ${\tau }_{\mathrm{A}}$ for the Andersen thermostat. The SWS compatibility is assured by using the effective temperature (Equation 11) for the Berendsen thermostat. In the case of the Andersen thermostat, the Maxwell-Bolztmann resampling is replaced with the SWS procedure.
• Then, a single or a few MD trajectories are collected with reasonably large time steps. The choice criterion is dictated by the integration method and corresponding frequency correction (Equation 8, see also Appendix A). In the cases of Verlet, velocity Verlet, and leapfrog integration schemes, the limit is given as $\Delta t\le {\left(\pi {\nu }_{\max }\right)}^{-1}$, where ${\nu }_{\max }$ is the maximal vibrational frequency of the system. If we take the H-F stretching frequency in hydrofluoric acid (${\nu }_{\mathrm{HF}}=4138$ cm⁻¹[62]), we get the maximal allowed time step of $\Delta t=2.6$ fs. Therefore, the time steps of around 1 fs are possible for most chemical systems. The total dipole moment of the molecular system is stored at every time step of the MD simulation.
• After the collection of the trajectory, the vibrational spectrum is computed as the FT of the dipole moment (Equation 2) or its velocity (Equation 3) autocorrelation function (Equation 1). The initial part of the trajectory is usually disregarded as the equilibration phase. The frequency resolution of the FT is given as $\Delta \nu =1/{\tau }_{\mathrm{tot}}$, where ${\tau }_{\mathrm{tot}}$ is the total duration of the trajectory (without the equilibration phase). An alternative way to transfer the autocorrelation function from the time domain into the frequency domain with arbitrary frequency increment is the rLSSA routine (Equation 6), albeit this procedure is much more computationally expensive than FFT, thus it makes sense to use it only for short ($\sim {10}^3$ steps) trajectories.
• Finally, the frequency correction (Equation 8) is applied, by transforming the frequency axis. Afterward, a tabulated scale factor for the corrected spectrum can be applied to account for the systematic errors in the quantum-chemical approximation.

To illustrate this procedure in action, we will consider a case of protonated methane (${\mathrm{CH}}_5^{+}$), with the vibrational spectrum obtained in Ref. [63]. For that, three $NVT$-MD trajectories at PBEh-3c level of theory with Andersen (${\tau }_{\mathrm{A}}=150$ fs) and Berendsen (${\tau }_{\mathrm{B}}=50$ fs) thermostats were computed. Each trajectory was 3 ps in length with 1 fs time step, and the ${\tau }_{\mathrm{SWS}}=4.2$ fs was taken from the harmonic frequencies at the same level of theory. After the MD simulations ended, for each of the trajectories, an rLSSA-computed spectrum was obtained with the first 1 ps of each trajectory disregarded for the equilibration phase. The spectra were smoothed by a convolution with Gaussian function with full width at half maximum (FWHM) of 50 cm⁻¹. For this, before performing the rLSSA procedure, the dipole velocity autocorrelation function was multiplied by a Fourier image of the frequency smoothing Gaussian, i.e., with another Gaussian distribution in the time domain with FWHM inversely proportional to the one in the frequency domain. Finally, the frequency correction and scale factor were applied. However, in addition to that, the total intensity in the spectrum was corrected by multiplying with the function of the form[14]

$$f\left(\nu \right)=\left\{\begin{array}{c}0,\mathrm{if}\nu <D,\hfill \\ 1-{\displaystyle \frac{D}{\nu }},\mathrm{if}\nu \ge D,\hfill \end{array}\right.$$

(18)

where $\nu$ is the frequency in cm⁻¹ and $D=383$ cm⁻¹ is the reaction energy needed for the signal to be detected. This type of correction was introduced in Ref. [14] and is not specific to the MD spectra but to the action spectroscopic methods in general. As we can see from Figure 5, the MD routine already produces a decent-looking spectrum. However, the application of the frequency correction, scale factor, and intensity correction (Equation 18) brings the theoretical spectrum quite close to the experimental one.

5. Conclusions

We have presented a general approach for computing the IR spectra of polyatomic systems using the AIMD approach implemented in the PyRAMD software. The approach invokes the following essential components: 1) inclusion of the NQEs by the usage of the SWS approach coupled to the thermostats; 2) frequency correction, allowing for larger time steps in the MD simulations; 3) pre-parametrized scaling factors, similar to those used for the harmonic frequency calculations. The synergistic effect of these components was demonstrated using the test case of the protonated methane. Despite this scheme being inferior to the more direct approaches, such as a combination of path integral MD for inclusion of the NQEs and of the smaller time steps in combination with highly accurate methods to obtain the PES, it was demonstrated that it could produce reasonable results on a fraction of computational time. Therefore, the presented algorithm can be recommended to provide theoretical counterparts to the experimental measurements.