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.
sakibmatin/calc-ir-spectra-from-lammps
This code shows how to calculate an IR spectrum from a LAMMPS MD simulation.
PythonMIT