Calculation of IR spectra using VMD IRspec Plug-in

The purpose of this tutorial to illustrate the basic capabilities of VMD IRspec plug-in. On example of the Wannier center dynamics from Car-Parrinello simulation we will calculate IR spectra of water molecule. Details of CPMD simulations are available as a separate tutorial.

First, we will need IONS+CENTERS.xyz (10mb) from that tutorial. It includes trajectory of the ionis and the Wannier centers (labeled as dummy atom X) in XYZ format. After proper visualization our water molecule will looks close to the picture on the right. Where Wannier centers represented by transparent green spheres. The final visualization state water_wannier.vmd file is available.

The upper two Wannier centers represent two lone pairs of oxygen. The other two centers represent bonds between oxygen and hydrogen. Please, note there are only 8 electrons in our system since we use pseudopotentials to describe core electrons. Finally, one can enjoy the dynamics on his/her own.

We need to supply several parameters from our simulations, namely Time between frames (simulation timestep), Temperature, maximum frequency, type of Correction, etc. In order to get the spectra, charges have to be assigned for all Wannier centers and atoms of interest. Despite of classical molecular dynamics there is no topology file available. The guess from CHARMM does not work, because the names of the elements do not match. Therefore we have to prepare charge.dat file by ourselves. It will looks like following:

# Charges for water molecule + Wannier centers O 6.0 H 1.0 H 1.0 X -2.0 X -2.0 X -2.0 X -2.0

Where numbers 6.0 and 1.0 represent effective atomic charges of corresponding “pseudoatoms”. Hence, the value of this charge for oxygen does not correspond to its true atomic charge 8. As we mentioned before, two core 1s electrons represented by pseudopotential. The charge of each Wannier center either -2.0 (two electrons) or -1.0 in case of spin polarized calculations.

Charge file should be loaded via “Utilities” –> “Load name” -> :Charge map from file” menu in the GUI.

You can optionally select Save to File checkbox in order to export spectral data in text format. The IR spectra will be plotted by Multiplot plug-in. One might export spectra into PostScript format or Grace program (on *nix systems) via File menu. Due to several reasons, our spectra is quite coarse and noisy. The water molecule is extremely small system, therefore the thermal fluctuations during molecular dynamics are high. Time scale of the simulation is small (just few picosecond). The fictitios masses distort the forces, shifting some frequencies by as much as few hudreds wavenumbers. Partly, this can be resolved by using Born-Oppenheimer Molecular Dynamics (BOMD), that is far more accurat than CPMD.

The lower the frequencies the longer time you have to sample them. For the frequencies above 1000 wavenumbers the trajectory about 20ps should be sufficient. Furthermore, the broad group of peaks in 0-500 cm^-1 range due to molecular librations and does not contribute to real IR spectra. These motions can be subtracted by special technique [link ??], however this goes beyond the scope of present tutorial.

In the spectra we can immediately identify three modes of water: asymmetric stretch in the region 1400-1800 cm^-1 and symmetric stretch+bend around 3600 wavenumbers. These signals are shifed comparing to the experimental values cm^-1 1595, 3657 and 3756 cm^-1 respectively. (DOI here)

For water molecules it is also possible to use partial charges from e.g. SPC classical water potential (H 0.42; O -0.84). BTW, very nice introduction into availabe water models is given by Martin Chaplin. This approach also gives you a reasonable approximation of the spectrum.