Efficiently Calculating Anharmonic Frequencies of Molecular Vibration by Molecular Dynamics Trajectory Analysis.

Two efficient methods, the Eckart frame algorithm and the multiorder derivative algorithm, for vibrational frequency calculation directly based on the raw data of atomic trajectory from the state-of-the-art first-principles molecular dynamics simulation are presented. The Eckart frame approach is robust to retrieve the full set of anharmonic fundamental frequencies of any molecule from the atomic trajectory for a sufficiently long molecular dynamics simulation at a temperature close to 0 K. In addition to the fundamental vibrational frequencies, the multiorder derivative approach is universal for the calculations of vibrational frequencies based on the molecular dynamics result in a wide range of temperatures. The accuracy, efficiency, and applicability of these two methods are demonstrated through several successful examples in calculating the anharmonic fundamental vibrational frequencies of methane, ethylene, water, and cyclobutadiene.

Introduction

The fundamental frequencies of molecular vibration, often simplified as fundamentals, are exclusively depended on the geometrical configuration and elemental composition of the molecule, and therefore, are considered as the molecular fingerprint. Information on fundamental frequencies is important to identify molecular species and to determine the molecular structure. The accurate calculation of fundamental frequencies is essential for a reasonable interpretation of the experimental observation of the molecular vibrational spectrum. The methodological development of the vibrational frequency calculation was initiated in the early part of the 20th century for explaining the experimental vibrational spectra of several small molecules.[1,2] The traditional frequency calculations are typically performed by the force field method with the force field parameters either fitted by experimental data or from ab initio calculations.[3−6] Because the principle of this realization is dependent on the knowledge of molecular normal vibrational modes,[1] it is often referred to as the normal mode analysis (NMA) method. At present, it has become a routine task to perform fundamental frequency calculations through the NMA method based on the second derivatives (SDs) of the adiabatic potential energy at the local energy minimum evaluated by ab initio calculations. In this way, the calculated frequencies are the harmonic frequencies. The algorithms for NMA and harmonic frequency calculations have been integrated in many of the commercial/non-commercial electronic structure software packages such as GAMESS, GAUSSIAN, HyperChem, MOLPRO, and others. However, the anharmonicity of the potential energy surface, the vibrational mode–mode coupling, and the Coriolis coupling are not taken into consideration in this kind of calculation. The accuracy of such calculations is poor,[7−13] and subsequent scaling of the calculated frequencies should be applied.[8,14−18] The scaling factor is usually a single-parameter empirical value between 0.89 and 0.99, and therefore, leads to a certain degree of uncertainty.[13,19−24]

Molecular vibration is the periodic motion of the constituent atoms around their mean positions in a molecule. Each atom oscillates around an equilibrium point in the molecular vibration. The interatomic bonds are non-rigid, may be stretched or compressed, thus have a variable bond-length. The bond angles that are formed between adjacent bonds may also change to cause bond bending, dihedral or torsion motions. This leads to the synergic motions of atoms in the molecule like a spring network, i.e. the intra-molecular atomic vibration at a finite temperature. The geometrical symmetry of a molecule can be described by a specific point group in the group theory, e.g. C2 for water, T for methane, and D6 for benzene.[25] According to the NMA theory, molecular vibration is described by a set of normal vibrational modes. Each of these normal modes belongs to an irreducible representation of a symmetry species in the point group. Each normal mode involves the associated motions of all atoms in the molecule. Based on the harmonic vibration theory, the atomic motion between different vibrational modes must be independent of each other and the potential surface must be a parabola. However, the potential surface is always slightly anharmonic in any real molecular system. This anharmonicity causes the coupling between different vibrational modes and results in anharmonic molecular vibration.

Several first-principles theoretical methods to calculate high resolution vibrational spectra by taking the anharmonic effect into consideration were proposed since the late 20th century, a few of which can sometimes even be compared with the experimental accuracy.[16] These methods are basically rooted from the vibrational self-consistent field (VSCF) approach and vibrational perturbation theory (VPT), including the configuration interaction VSCF,[26] correlation-corrected VSCF,[27,28] vibrational multi-configuration self-consistent field,[29] size-extensive VSCF method,[30] vibrational second-order Moller–Plesset perturbation theory (VMP2 or VPT2),[31−34] vibrational coupled cluster theory,[35,36] vibrational multi-reference coupled cluster theory,[37] etc. These high-accuracy methods are only limited to small molecules except VSCF and VPT2. A fully automated code for an anharmonic vibration study with VPT by the second-order Moller–Plesset perturbation had been developed.[34] The VPT2 method had demonstrated its accuracy for anharmonic vibrational frequency calculations for molecular systems consisting of dozens of atoms.[38−40] A VSCF approximation with the second-order perturbation correction (VSCF-PT2) algorithm[28] had been recently applied to anharmonic vibration calculations for a few large biomolecules with several hundred atoms.[41−43] Nevertheless, both VPT2 and VSCF-PT2 suffer the problem of the resonance effect,[43−45] and several terms in the formalism of the VSCF algorithm show nonphysical size dependence.[30] A more general computational method to solve vibrational Schrödinger equations for larger molecules and solids remains to be developed.

All information on molecular vibration is encrypted in the trajectory of the atomic motion evolved over time. Unfortunately, it is not feasible to directly record the atomic trajectory by any means of experimentation due to the micro-amplitude and rapid displacement of the atoms in the molecular vibration. Molecular vibration energy is quantized. Molecules can change their vibrational state and interact with the outside world through absorbing or releasing energy that matches their vibrational frequencies. Based on this principle, experimental scientists try to use a variety of ways to measure the amount of energy exchanged by molecular vibration with the outside world to probe the behavior of the molecular vibration. Various spectroscopic techniques were invented owing to these efforts, such as infrared (IR) spectroscopy, ultraviolet spectroscopy,[46] Raman spectroscopy, photon fluorescence spectroscopy, cold ion spectroscopy,[42,47] etc. Because it is hardly possible to conduct a comprehensive study on the molecular vibrational behavior only using a single spectroscopic device in most cases, it is always a tough task for experimental study of the vibrational properties even for many well-known molecules. On the other hand, there is no difficulty to record the atomic trajectory in theoretical simulations. The only foreseeable obstacle might just be the limitation on data storage of the computer. The simplest and most straightforward way to reasonably study the anharmonic behavior in the molecular vibration can be done by molecular dynamics (MD) simulations. MD was first developed by Fermi and co-workers in the mid 50s of the last century,[48] and successfully applied to simulate the elastic collisions between hard spheres with an IBM 704 computer by Alder and Wainwright in 1960.[49] MD is designed to simulate the physical movements of atoms and molecules. The trajectories of atoms and molecules are determined by numerically solving Newton’s equations of motion for the system. The forces between the constituent atoms are calculated either using the empirical interatomic potentials or by directly solving the Schrodinger equation.

Theoretically, the dynamic behavior of atoms in any physical system can be simulated by MD as long as the exact information on the interaction between atoms is known. An empirical interatomic potential function incorporated correctly with non-harmonic and polarization effects is required in the classical MD simulation for vibrational frequency calculations. As long as the potential function is accurate enough, detailed information on molecular vibration can be recorded in the atomic trajectory of a sufficiently long MD simulation. However, such a potential function is very difficult to find for most interested molecules. The first-principles MD (FPMD)[50] or ab initio MD methods have emerged in the last three decades to study the dynamic problems of atomic and molecular systems in quantum mechanical precision. FPMD within the framework of density-functional theory (DFT)[51,52] deals with the interatomic forces by accurate electronic structure calculation, the effects of charge redistribution, potential energy surface anharmonicity, and polarization can be naturally taken into account in such calculations. The state-of-the-art FPMD offers one of the most appropriate choices for molecular vibration simulations.

The earliest application of computer simulation to molecular vibrational spectra was attempted by Noid and co-workers for anharmonic two- and three-dimensional model systems in 1977.[53] Adams et al. studied the vibration of icosahedral C60 by a simplified quantum-molecular-dynamics simulation.[54] The thermal vibration of hydrogen on the Si(111) surface was simulated by Car–Parrinello MD (CPMD); the vibration frequencies of Si–H stretching and bending, and Si–Si–H angle bending modes were calculated by the Fourier transformation of the simulated atomic trajectories near the surface.[55] Wathelet et al. tested the possibility to calculate the vibrational frequencies of diatomic molecules by the first-principles CPMD simulation.[56] The IR absorption spectrum of water was calculated by Fourier transformations of the velocity autocorrelation function[57] and the dipole moment autocorrelation function[58,59] from MD simulations, respectively. The hydrogen molecular vibration in crystalline silane under high-pressure was simulated by ab initio canonical ensemble MD at 300 K.[60] The phenomenon of vibrational spectral diffusion in liquid methanol was captured by the CPMD simulation.[61] The Fourier transform of atom-velocity and bond-length autocorrelation functions from the CPMD simulation was carried out to study the molecular vibration of methane in clathrate hydrate.[62] Analysis on the electronic excitation state from the simulated trajectory of the Born–Oppenheimer MD (BOMD) simulation was realized.[63] On the whole, the reported MD studies are mostly focused on some specific aspects of molecular vibration. The calculation and analysis methods adopted are limited only to certain special purposes. The feasibility of accurately calculating the anharmonic frequencies of all vibrational modes of a molecule directly by the MD trajectory has not been verified. A universal, easy-to-use method for this purpose has yet to be developed.

The Fourier transform provides a primary theoretical tool for spectral analysis. A Fourier transform maps a time series g(t) into a frequency spectrum G(ν) using

The MD trajectory consists of discrete time series of atomic coordinates. The high-efficient fast Fourier transform (FFT) algorithm is usually used for such discrete Fourier transform. All information about the vibrational frequencies should be obtained in principle through the FFT of a sufficiently long MD trajectory in molecular vibrational simulations. However, the problem that complicates the matter is that the equilibrium positions of atoms in vibration are constantly changing because of the overall rotating, random tilting, and drifting motions of the molecule in a vacuum. The regular atomic vibrational signal and the irregular random signal are mixed together in the MD trajectory. This problem will cause serious noise in the spectrogram by the direct FFT of the raw MD trajectory. As a typical example, the influence of noise on the obtained spectrogram of methane is shown in Figure . Figure a presents the raw atomic trajectory for the y coordinate of a hydrogen atom in the molecule between 5 000 000 and 5 200 000 steps in the FPMD simulation. The inset of Figure a shows the enlarged trajectory of the small rectangular area to show that the faint vibrational signal is buried in the strong noises. The direct FFT of this raw trajectory is given in Figure b. It is seen that the extremely strong noise frequencies occupy the low-frequency end on the left side of the figure. Still, the fundamental frequencies of the normal vibration modes ν1 and ν3 in the vicinity of 3000 cm–1 are distinguishable, and the other two fundamental frequencies, ν2 and ν4, can also be identified in the inset enlarged three hundred times. The pure vibrational trajectory between 5 090 000 and 5 100 000 MD steps of this atom after noise reduction using the Eckart frame algorithm is shown in Figure c which is retrieved from the inset area of the raw MD trajectory in Figure a. Figure a shows that the maximum amplitude of the raw MD trajectory is about 0.6 Å. Figure c shows that the maximum amplitude of the pure vibrational signal is less than 0.002 Å. Therefore, the noise is a few hundred times stronger than the pure molecular vibration in the raw MD trajectory. Despite the serious noise, the information on molecular vibration is completely kept in the raw MD trajectory as shown in the inset of Figure a and in the Fourier transform spectrogram of Figure b. Evidently, the calculation of highly accurate molecular vibrational frequencies from the MD trajectory depends on powerful noise reduction and signal enhancement techniques.

Figure 1

Spectrogram calculation based on the MD trajectory. (a) The raw MD trajectory of the y-coordinate of a hydrogen atom in methane by FPMD simulations. The trajectory in the small rectangular frame area is enlarged in the inset which shows that the faint vibrational signal is buried in strong noise. (b) FFT of the raw MD trajectory in (a). (c) The pure vibrational trajectory of the same atom after noise reduction using the Eckart frame algorithm.

Although the MD simulation is one of the most natural and direct methods to study molecular vibration theoretically, it is usually still rather difficult to directly retrieve the vibrational information from the raw data of the MD trajectory because of the serious noise interference. The great potential of MD simulation in molecular vibration research has not been fully demonstrated. Therefore, it is highly necessary to find efficient numerical algorithms for MD trajectory analysis. The purpose of the present study is to develop accurate, robust, and easy-to-use methods for retrieving molecular vibrational information based on the result of FPMD simulations.

Algorithms for Molecular Vibrational Frequency Calculation by MD Trajectory Analysis

Two algorithms, the Eckart frame algorithm and the multiorder derivative (MOD) algorithm, will be developed for retrieving the vibrational information from the raw data of the MD trajectory through noise reduction and signal enhancement in the following two subsections.

Calculation of Molecular Vibrational Frequencies using the Eckart Frame Algorithm

It is seen from Figure that there are serious noises in the MD trajectory resulted from the rotation, the random drift, and tilting motions of the molecule as a whole. In order to eliminate these noises, the algorithm for vibrational frequency calculation using the Eckart molecule-fixed frame includes the following steps.

First of all, the random drift motion can be eliminated by the translational transformation of the coordinate system through always setting the mass center of the molecule as the origin of the coordinate system in each MD step. For a molecule with n atoms assigned by coordinates r⃗ and masses m (i = 1...n), its mass center is given by

Secondly, a molecule-fixed Eckart frame coordinate system is adopted to remove the effect of global tilting and rotation of the molecule. Eckart proposed a moving Cartesian coordinate system to separate the rotational and vibrational behaviors of polyatomic molecules to the maximum extent possible, i.e., the Eckart frame.[64] The instantaneous orientation of the Eckart frame changes with the instantaneous coordinate vectors of the atoms in the molecule at the moment. All atoms look seem to be undergoing the vibration motion around their equilibrium positions in the Eckart frame. The condition of the Eckart frame is[64−66]Here, r⃗e is the initial coordinate vector of the ith atom at equilibriuma⃗1, a⃗2, and a⃗3 are the Cartesian basic vectors of the real space R3. The instantaneous coordinate vector r⃗(t) of the ith atom in real space R3 is

The instantaneous coordinate vector r⃗*(t) of the ith atom in the Eckart frame at time t is given bye⃗1(t), e⃗2(t), and e⃗3(t) are the three unitary Eckart vectors[65] that define the instantaneous Eckart frame. The relationship between r⃗*(t) and r⃗(t) is

For the sake of brevity, the time dependence of the atomic coordinate will be hidden hereafter. U is the direction-cosine matrix for the coordinate transformation from the R3 real space to Eckart space. It is expressed in quaternion Q = (q1, q2, q3, q4) algebra by[67]with |Q|2 = 1. Through minimizing the sum of mass-weighted squared deviations of the atomic coordinates of the displaced and the equilibrium configurations, it was deduced that Q is the eigenvector with the minimal eigenvalue of the following symmetrical matrix C(66)

The matrix elements of C are calculated usingHere, the ξ and ζ variables are defined as

The physics behind this algorithm is that the Eckart frame rotates together with the global rotation of the molecule; thus, the atoms seem to engage only in the vibrational motion in the instantaneous Eckart coordinate system. In this way, the details of molecular vibration are completely kept in the time series r⃗* (i = 1...n), and all information on the global molecule rotation is kept in the time series of the Eckart frame, e⃗(t) (k = 1...3).

Finally, it is a simple task by the FFT of these atomic trajectories in the Eckart frame to calculate the vibrational frequencies of the molecule as soon as the time series of r⃗* (i = 1...n) is ready. The spectrograms for each of the three coordinate components of the ith atom are calculated usingwhere Δt is the timestep in the MD simulation, N represents the total MD steps, and i in the exponent is the imaginary unit. The spectrogram of the ith atom is given by

Calculation of Molecular Vibrational Frequencies using the MOD Algorithm

The Eckart frame algorithm presented in the former subsection enables one to effectively separate vibration from rotational motion through intensely compressing the rotational effect. Because the fundamental vibrational modes of a molecule are only related to the pure vibrational motion of atoms, the method is especially suitable for fundamental frequency calculations. However, the algorithm may fail at high temperature when the strict Eckart frame is hard to find due to the serious molecular geometry distortion.[66] The Eckart frame algorithm is also not applicable to analyze the rovibrational spectrum caused by the interactive coupling between the rotational and vibrational motions. An alternate solution will be presented in this subsection.

Vibration is the periodic motion of an object around its equilibrium position. In addition to the periodic change of the atomic coordinate, the velocity, acceleration, and all the MODs of the atomic coordinate trajectory also periodically vary with the same frequency as the atomic coordinate in molecular vibration. Unlike the constant drift of the atomic equilibrium position in the MD simulation of an isolated molecule, the equilibrium points of velocity, acceleration, and any MODs of the atomic coordinate trajectory are always zero. Hence, the problem of the equilibrium point drift of the vibratory variable in molecular vibration can be solved effectively. For this reason, the FFT spectrum of the MOD of the MD trajectory should provide more accurate information on the vibrational frequency than the atomic coordinate trajectory. Mathematically, the MD trajectory is the superposition of the vibrational motions, non-vibrational motions, and random noises. The coordinates of all atoms in the molecule can be written aswhere v⃗(t) is the contribution of the vibrational motion to r⃗(t). Its three components can be written asM is the total number of vibrational modes, ω is the angular frequency of vibrational mode j, and φ is the phase difference of the mode. A strong anharmonic vibration gives rise to frequency multiplication by overtones o(t), and combination bands c(t), upon the fundamental vibrational motions to the system. μ⃗(t) includes all non-vibrational motions and noises. This item represents the random thermal fluctuation of the system, its lack of periodicity. Because NVE ensemble is adopted, no periodic noise source is involved in the present FPMD simulation. Thermal fluctuation leads to an irregular drift and tilt of the molecule in a short period of time. All of these aperiodic noises can be described using a polynomial expansion function of time t

Based on eqs –17, a noise reduction algorithm by the successive derivative of the atomic coordinate trajectory with respect to t is suggestedeqs and 19 show that the derivative operation will change π/2 in the phase angle and enhance the amplitude by 2πω to the periodic part of the derivative function. All information on the vibration frequencies will be retained. As for those non-vibrational signals and noises, the first derivative of r⃗(t) will remove the signal by transient random jumps of the molecule. The SD eliminates the effect of global molecule translation at a steady speed, and the third derivative removes the irregular accelerating and tilting motions, etc. It is seen that the net effect of derivative operations will lead to the significant enhancement of the vibrational signal/noise ratio in the derivative functions. Therefore, a much clean molecular vibrational spectrogram can be recorded through the Fourier transform of these MOD functions of the MD trajectory.

The numerical five-point central difference equations[68] were used for MOD operation in this study

The vibrational spectrogram is calculated using

By this algorithm, all information related to the vibrational frequencies will be preserved.

Details of the FPMD Simulation

All FPMD simulations in this study were carried out using the CP2K code[69] which is an open-source software package for atomic and molecular simulations. The quickstep algorithm[70] by Gaussian and plane wave formalism, using a dual basis of atom centered Gaussian orbitals and plane waves for regular grids, is implemented for solving DFT Kohn–Sham equations in the code. FPMD simulations were performed using the quickstep module of CP2K in the scheme of BOMD. The double-zeta valence polarized Gaussian basis sets optimized for the Goedecker–Teter–Hutter pseudopotentials[71] and the Pade exchange–correlation functional in the local density approximation (LDA) were chosen for C and H elements in methane, cyclobutadiene, and ethylene calculations. The triple-zeta valence Gaussian basis sets with two sets of polarization functions optimized for GTH pseudopotentials and the Pade LDA exchange–correlation functional were used for H and O elements in H2O molecule calculations. A 280 Ry cutoff for the finest grid level of plane waves was used for GTH pseudopotentials. Orbital transformation methods by the conjugate gradient and direct inversion in the iterative subspace were adopted for wave function optimization in molecular geometry optimization and MD simulation, respectively.[70]

The supercell modeling method is adopted to build the structure model for the FPMD simulation; the molecule is placed in the center of a vacuum cubic box. The edges of these boxes are set to 20 Å for ethylene and cyclobutadiene, 15 Å for water, and 10 Å for methane, respectively. Both periodic and aperiodic boundary conditions were tested. The Kohn–Sham equation under aperiodic boundary conditions was solved using the wavelet algorithm. The decoupling technique[76] was adopted to eliminate the effect of the periodic charge images in the evaluation of the Coulomb interaction for the infinitely replicated periodic system. The vibrational frequencies calculated using these two boundary conditions with the above model settings are almost the same; thus, the reliability of these models for simulating an isolated molecule is verified. Hereafter, all the reported results were simulated using the periodic boundary conditions. As a large displacement from atom’s equilibrium position will cause an overtone phenomenon, the MD simulation at extreme low temperature is preferable to calculate the fundamental frequencies of molecular vibration. The microcanonical ensemble (NVE ensemble), in which the total energy, particle number, and volume of the system remain unchanged during the FPMD simulation, was adopted throughout this study. It was suggested that a timestep shorter than 0.5 fs is essential for accurate spectrum calculations by the MD simulation.[77] A variety of timesteps from 0.01 to 0.5 fs were tested for the system stability during the FPMD simulation with the NVE ensemble in this work. It is found that a timestep of 0.1 fs will keep the system energy stable up to a few million steps when the system temperature is close to 0 K, which will be enough for fundamental frequency calculation of small molecules. Before MD simulation, the ground-state equilibrium geometry of the molecule was obtained using the geometry optimization routine in the CP2K code. The convergence criterion of force is set to 1.0 × 10–4 hartree/Bohr. The optimized equilibrium geometries of the molecules in this work are listed in Table . It is seen that the optimized configurations agree well with the experimental reports. Then, the NVE FPMD simulation was set up from the optimized ground-state configuration by slightly displacing the atoms from their equilibrium positions according to the initial temperature. No atomic velocity or temperature rescale was performed throughout the simulation. This realization is essential for simulating the intrinsic vibrational behavior of a molecule under the conditions without outside interference. Starting from 0 K, the final stabilized temperatures of our simulated systems are 0.16, 0.17, 0.35, and 0.12 K for water, ethylene, cyclobutadiene, and methane, respectively, in the fundamental frequency calculations. The closer the optimized molecule configuration is to the ground state, the lower the stabilized system temperature. The molecular vibrational frequencies are calculated using the atomic coordinate trajectories after the system temperature was stabilized. The total time of the MD simulation for spectrum calculations is 0.29, 0.29, 0.24, and 1.80 ns for water, ethylene, cyclobutadiene, and methane, respectively, these correspond to 2.4–18.0 million MD steps.

Table 1

Optimized Molecular Geometries for the FPMD Simulation (Length and Angle Units are in Å and Degree, Respectively)

molecule	geometry	optimized	experimental	error (%)
cyclobutadiene	rCC1	1.347	1.380[72]	–2.4
	rCC2	1.571	1.560[72]	0.7
	rCH	1.094
	α1	134.62
	α2	135.38
ethylene	rCC	1.343	1.335[73]	0.6
	rCH	1.098	1.090[73]	0.7
	αHCC	121.51	121.7[73]	–0.2
methane	rCH	1.102	1.086[74]	1.5
	αHCH	109.47
water	rOH	0.971	0.959[75]	1.3
	αHOH	104.79	103.9[75]	0.9

Verification of the Numerical Algorithms for Calculating Molecular Vibrational Frequencies

The obtained spectrograms of water, methane, ethylene, and cyclobutadiene in this work are presented in Figures –5, respectively. The vibrational frequency calculations were conducted by applying eq to an arbitrarily selected hydrogen atom in these molecules. The quantitative comparisons of the present calculated fundamental frequencies with the published experimental results are given in Tables –5. The fundamental frequencies in the “Eckart calc.” column are calculated using the Eckart frame algorithm, the data in the “Error” column gives the relative error between the present calculated and the experimental results. It is seen from Tables and 3 that the calculated fundamental frequencies of water and methane match well with the experimental measurements. Table shows that the frequencies of ν1 and ν3, for the stretching vibrations along the C–H bond of methane, are slightly overestimated; meanwhile, ν2 and ν4 for the waging vibrations are underestimated, which means that the present DFT description somewhat enhanced the tension–compression strength of the C–H bond but weakened the bending strength between C–H bonds. The calculated fundamental frequencies of ethylene match closely with the infra-red spectrum analysis by Arnett and Crawford[78] in Table . The vibration modes of ω1, ω9, and ω11 are typical stretching vibrations along the C–H bond. The calculated frequencies of these modes are all higher than those obtained in the experimental measurements, which is similar to the enhancement of C–H tension–compression strength in methane. The mode of ω2 is for the stretching vibration between the two carbon atoms in the middle of ethylene, and ω10 is for the waging vibration between the neighboring C–H bonds of the two carbons, respectively. The calculated frequencies of these two modes reveal that the tension–compression enhancing and bending weakening are also true for this molecule. The fundamental frequencies of cyclobutadiene by the SD of MD trajectories, using the same MD trajectories as for the Eckart frame algorithm, are presented in the “SD calc.” column in Table for a quantitative comparison with the results from the Eckart frame algorithm. It is seen that both the Eckart frame and MOD algorithms have the same accuracy for the fundamental frequency calculation of molecular vibration, but the latter approach is simpler in concept and much easy to implement.

Figure 2

FFT spectrogram of water using the Eckart frame algorithm. The spectrogram was calculated for an arbitrarily selected hydrogen atom in the molecule.

Figure 5

FFT spectrogram of cyclobutadiene using the Eckart frame algorithm. The spectrogram was calculated for an arbitrarily selected hydrogen atom in the molecule.

Table 2

Calculated and Experimental Fundamental Frequencies of Water (Unit in cm–1)

	expt.	Eckart cal.	error (%)
ω1	3755.92[79]	3819.3	1.7
ω2	3657.05[79]	3711.1	1.5
ω3	1594.75[80]	1560.6	–2.1

Table 5

Calculated and Experimental Fundamental Frequencies of Cyclobutadiene (Unit in cm–1)a

no.	sym.	exp.	Eckart calc.	SD calc.	err. (%)	active	action
1	b2g	531[83]	505.7	505.3	–4.8	Raman	CH OPLA deform
2	b3u	569[86]	556.5	556.9	–2.1	IR	CH OPLA
3	b2u	719[86]	709.9	710.0	–1.4	IR	ring deform
4	b3g	723[83]	827.8	827.6	14.5	Raman	ring deform
5	ag	989[83]	935.0	935.1	–5.5	Raman	CH deform
6	b1u	1028[86]	1002.1	1002.2	–2.5	IR	CH deform
7	ag	1059[83]	1093.7	1093.6	3.2	Raman	CH deform
8	b3g		1121.9	1121.8		Raman	C–C stretch
9	b2u	1242[86]	1198.8	1198.8	–3.5	IR	C–H deform
10	ag	1678[83]	1611.0	1611.1	–4.0	Raman	C=C stretch
11	b3g	3093[83]	3121.4	3121.3	0.9	Raman	CH stretch
12	b2u	3105[84,86]	3139.3	3139.3	1.0	IR	CH stretch
13	b1u	3124[85,86]	3159.7	3159.9	1.1	IR	CH stretch
14	ag	3140[83]	3172.6	3172.7	1.0	Raman	CH stretch

The data in the SD calc. column are calculated using the second derivative of the MD trajectory.

Table 3

Calculated and Experimental Fundamental Frequencies of Methane (Unit in cm–1)

band	ν1	ν2	ν3	ν4
symmetry	A1	E	F2	F2
Eckart calc.	2969.8	1488.1	3097.7	1253.6
exp.	2917[81]	1533.6[81]	3019.7[81]	1310.8[82]
error (%)	1.8	–3.0	2.6	–4.4

Table 4

Calculated and Experimental Fundamental Frequencies of Ethylene (Unit in cm–1)

band	ω1	ω2	ω3	ω4	ω5	ω6
symmetry	A1g	A1g	A1g	A1u	B1g	B1g
Eckart calc.	3076.5	1661.6	1317.9	1017.4	3142.2	1162.0
exp.[78]	3019.3	1623.3	1342.4	1027.0	3272.3	1236.0
error (%)	1.9	2.4	–1.8	–0.9	–4.0	–6.0

FFT spectrogram of methane using the Eckart frame algorithm. The spectrogram was calculated for an arbitrarily selected hydrogen atom in the molecule.

FFT spectrogram of ethylene using the Eckart frame algorithm. The spectrogram was calculated for an arbitrarily selected hydrogen atom in the molecule.

In general, it should be enough to obtain the full set of fundamental frequencies by performing FFT of the MD trajectory for any one of the atoms with the Eckart frame algorithm because of the overall synergic movement of all atoms in the molecule. However, exceptions may occur for some molecules with a high structural symmetry. For example, the carbon atom is immobile for the vibration modes with A1 and E symmetries, but all hydrogen atoms are active for any vibration mode in methane. The trajectory of the lightest atom in the molecule is assumed to carry the most detailed information on the molecular vibration. Figures –5 are obtained using the Eckart frame algorithm with the MD trajectory of an arbitrarily selected hydrogen atom in these molecules.

Figure presents an example to demonstrate the effect of the MOD noise-reduction algorithm for the vibrational spectrum calculation of ethylene. The MD simulation was initialized in the NVE ensemble at 10 K in this case. The system temperature was finally stabilized at 6.77 K. The thick black lines are the raw MD trajectory of the x-coordinate for a carbon atom in ethylene, and the thin blue lines show the first to fourth derivative curves of the raw MD trajectory in sequence from Figure a–d. The amplitudes of the curves for the first to fourth derivatives are enlarged 35.0, 116.0, 233.0, and 437.0 times, respectively, to facilitate direct visual comparison with the raw MD trajectory. The insets in these figures are the FFT spectrograms of the MODs of the raw trajectory. It is seen from Figure a that the vibrational signal is greatly enhanced against the noise by the first derivative operation, but the noise interference is still cleanly seen in the first derivative curve. The influence of noise is also reflected in the left lower-end of the inset spectrogram near zero cm–1. The noise interference is almost completely removed in the SD curve and in the corresponding spectrogram as shown in Figure b. By comparing the derivative curves and the inset spectrograms from Figure a–d, it is seen that the multiorder derivation has a remarkable effect on the enhancement of the periodical signals, especially the high frequency signals. As is shown in eqs and 19, the derivation operation does not cause any loss of vibrational frequency information other than changing the amplitude of the vibrational signal. The MOD algorithm has a special advantage in the spectral calculation of the weak high-frequency vibrations, for example, the polyad structures of methane in the region from several thousands to tens of thousands cm–1.[87−90] Because the initial system temperature of the MD simulation is 10 K in this example, the spectrograms shown in the insets of Figure are actually the rovibrational spectra of ethylene, which are reflected from the broadening of the peaks around the fundamental frequencies.

Figure 6

Noise reduction of the MD trajectory using the MOD algorithm. The thick black lines are the raw MD trajectory of the x-coordinate for a carbon atom in ethylene, and the thin blue lines show the first to fourth derivative curves of the raw MD trajectory of this atom in sequence from (a–d). The amplitudes of the curves for the first to fourth derivatives are enlarged 35.0, 116.0, 233.0, and 437.0 times, respectively. The obtained spectrogram is given in the inset of each figure. The undifferentiated spectrogram is shown in the upper inset of (a).

Discussion

Molecular vibration is intrinsically anharmonic because of the anharmonicity of the potential surface in real molecules. FPMD can simulate the vibrational process of real molecules without approximation. The simulation takes into account the anharmonic effect of molecular vibration in the most natural way. The respective scopes of application of the vibrational frequency calculation methods developed in this study, as well as the considerations in the application, are discussed below.

Performance Comparison between the Eckart Frame and MOD Algorithms

The Eckart frame algorithm of the vibrational frequency calculation developed in this work enables the separation of the molecular vibration from molecular rotational and random motions to the maximum extent possible. This method is most suitable for calculating the anharmonic fundamental frequencies, which can be conveniently realized to any molecule with the atomic coordinate trajectory simulated by FPMD at a temperature close to 0 K. It is a robust method, and is very suitable for the construction of large scale molecular databases based on the idea of materials genome research.[91−93] The MOD algorithm is inspired by the analysis that the MD atomic coordinate trajectory comprises periodic sinusoidal vibrational signs, aperiodic non-vibrational signals, and random noises. The derivative operation specializes in compressing aperiodic signals and intensifying the periodic signals. The higher the derivative order, the more significant the effect is. A visual comparison of these two algorithms is presented in Figure . The FPMD simulations of an isolated methane molecule were performed for three stabilized temperatures Ts at 0.41, 9.70, and 198.02 K, respectively. It is seen from the left column spectrograms, obtained using the Eckart frame algorithm, in Figure a,c,e that there are very sharp peaks at the positions of the fundamental frequencies, but the rovibrational, combination bands, and overtone features are seriously distorted because of the suppression effect of the algorithm on the rotational motion. In contrast, all rovibrational, combination bands, and overtones are gradually resolved with the increase of temperature in the spectrograms calculated using the MOD algorithm as shown in the right column FFT spectrograms of Figure . It is seen that the fundamental frequencies can always be determined by the local strongest peaks in the recorded spectrogram at various simulated temperatures with the Eckart frame method. As the simulation temperatures rise, it becomes increasingly difficult to distinguish the fundamental frequencies in the obtained spectrogram using the MOD algorithm. However, the latter spectrogram contains more information on the vibrational property of the real molecule.

Figure 7

Temperature dependence of the FFT spectrograms obtained for an arbitrarily selected hydrogen atom in methane. The spectrograms in the left (a,c,e) and right (b,d,f) columns are recorded using the Eckart frame algorithm and the SD algorithm, respectively. The features of rovibrational bands (d), combination bands, and overtones (f) are emerged successively as temperature increases in the spectrograms by the latter algorithm.

Correlation between Calculation Accuracy and the Simulation Time

For a dynamic system experiencing a sufficient long running time, almost all points in its phase space will be visited and all dynamic properties of the system will be included in system’s evolutionary trajectory according to the Ergodic theory in statistical physics. There are two reasons for using a sufficiently long-time trajectory in vibrational frequency calculations. Firstly, the resolution Δν (unit in cm–1) of the FFT spectrogram is depended on the total MD steps N, the MD timestep Δt, and the light velocity c by

It is seen that the longer the MD trajectory the higher the frequency resolution. For example, the MD trajectories of 330 000 and 3.3 million steps bring about 1.0 and 0.1 cm–1 frequency resolutions, respectively, when the MD timestep is 0.1 fs. Secondly, the MD simulation in this study starts with a small perturbation to the equilibrium molecular configuration. It is most likely that only a few vibrational modes are excited at the initial MD steps. The other modes can be gradually excited through modes coupling during the simulation. A sufficiently long-time trajectory is essential to include the information on all vibrational modes. At higher simulated temperatures, this process can be accelerated, but the spectrogram will become too complex to analyse.

Target Atom Selection in the Case of Macromolecules

In the examples of the former section, only an arbitrarily selected hydrogen atom is used for the FFT spectrogram calculations. This target atom choice method in the vibrational spectrogram calculation is feasible for small, single branched, and covalently bonded molecular systems. But, it might be failure for large organic molecules and biomolecular systems that are composed of multiple functional groups with many complex multi-branched, multi-stranded, and weak hydrogen bonding structures in which the interactions between different functional groups or branches are very weak. In such a case, the selection of a light atom from each of the functional groups or branches for molecular vibration analysis will show better results.

Intensity Comparison between the Calculated and Experimental Spectra

The geometrical configuration and elemental composition are uniquely correlated with a specific set of fundamental frequencies. However, the experimental spectrum intensity is closely related to the experimental methods and test conditions. For example, the intensities of the a1g and e2g vibration modes of benzene can be detected by Raman spectroscopy, but the intensities of these modes are zero in IR spectroscopy. On the contrary, the intensities of the three benzene e1u modes can only be observed by IR spectroscopy rather than the Raman spectroscopy. Temperature is also a critical factor of the experimental spectral intensity. Therefore, in addition to the atomic trajectory, additional information on polarizability or dipole moment is needed to reproduce the experimental intensities of IR and Raman spectra. Because the purpose of the present study is to develop general methods for the molecular vibrational frequency calculation only based on the MD trajectory, the calculated spectrogram intensity is not suitable for direct comparison with the relevant experimental results.

Information on fundamental vibrational frequencies is important for detection, identification, and characterization of organic and inorganic molecules. Molecular vibrational frequencies carry abundant and subtle information of the static and dynamic properties of the molecule. Theoretically, the successful molecular vibrational simulation can be used as an important criterion for verifying the accuracy and validity of various first-principles theories and methods. Such research works will surely push the development of the first-principles calculation to a new height.

Conclusions

In this paper, two algorithms for calculating molecular vibrational frequencies based on the raw data of the atomic coordinate trajectory by the MD simulation are presented. The Eckart frame approach enables the separation of the vibrational and rotational motions of atoms to the greatest extent. It is robust to calculate the full set of fundamental frequencies of any molecules by taking the anharmonic effect into consideration using the atomic trajectory of a sufficiently long time MD simulation near 0 K. The MOD approach fully takes advantage of the enhancement effect of derivative operation on periodic signals. The higher the order, the better the effect of the high-frequency enhancement is. This approach is generally applicable to the calculation of molecular vibrational frequencies based on the MD trajectory in a wide temperature range. Several successful examples are presented to demonstrate the applicability of these two algorithms and the potential of the state-of-the-art FPMD simulation in the study of molecular vibrational properties.