leucinw / calc-ir-spectra-from-lammps Goto Github PK

To cite: https://zenodo.org/badge/DOI/10.5281/zenodo.154672.svg

This code shows how to calculate an IR spectrum from a LAMMPS MD simulation.
The LAMMPS input file in.infrared outputs the instantaneous net dipole moment
of the system versus time. The python script calc-ir-spectra.py calculates the
dipole moment's autocorrelation function and then performs a Fourier transform
to obtain the IR lineshape. The spectra itself is then calculated by
multiplying the lineshape by an electromagnetic field dependent prefactor.

in.infrared:
Here, I give the LAMMPS calculation that's needed to output the file
dipole.txt, which is the instantaneous value of the dipole moment. This
calculation is given under the comment, '#calculate dipole moment vector'. The
calculation will work as long as your atoms have charge and the net charge of
the system is 0. One thing to consider changing is the '2' in the 'fix
printdipole...' line, which is how often the dipole moment is printed. At the
moment, I print every 1 fs, (which is 2 times the timestep of 0.5 fs). I give
comments regarding how frequently you need to print and how long you need to
run to get the spectrum you want within in.infrared. The simulation is of
MOF-5, using the Dubbeldam force field (doi.org/10.1002/ange.200700218)

calc-ir-spectra.py
Here, I take dipole.txt and I calculate the IR spectrum. After the code
calculates this autocorrelation function, it prints it to the file
autocorr.txt. If you want to rerun calc-ir-spectra.py and don't want to repeat
the calculation of the autocorrelation function, you can set the variable
'autocorrelation_option' to 2, which will then load in a previously-calculated
autocorr.txt. The only things you may want to consider changing are in the
'Inputs' section, where you should change the temperature and the wavenumber at
which you want to cut off the graphs.

In the directory 'output-files', I give the files generated by:
1. Running the LAMMPS code
2. Running calc-ir-spectra.py on the LAMMPS code output
3. Plotting the simulated spectra against experiment
(doi.org/10.1039/B603664C). Note that regarding the seemingly spurious
simulated peak around 3000 cm^-1 that is due to CH stretching, the
experimentalists stated, "We have to note that it is not obvious why the
aromatic CH stretching modes of the framework bdc ligand in pure (and solvent
free) MOF-5 exhibit exceptionally weak intensities and are thus not visible in
our spectra (Fig. 2a), but this matches reported MOF-5 FT-IR data."

As for whether you want to look at the spectra without or with the quantum
correction, that's up to you. See doi.org/10.1021/jp034788u,
doi.org/10.1063/1.441739, and doi.org/10.1080/00268978500102801 for
descriptions of this. If you are comparing to an experimental spectra, you
probably want to use it, though several equally-arbitrary choices of quantum
corrections exist. I used the one recommended by doi.org/10.1021/jp034788u.

A very incomplete list of some good literature on calculating IR spectra:
doi.org/10.1063/1.441739
doi.org/10.1063/1.461069
doi.org/10.1021/jp0014925
doi.org/10.1021/jp034788u
doi.org/10.1039/c3cp44302g

This code has been shown to work using LAMMPS (14 May 2016), Python 2.7.11 and
3.5.1, NumPy 1.10.4 and 1.11.0, Matplotlib 1.5.1, and SciPy 0.17.0 and 0.17.1.