Infrared spectra
This tutorial shows how to calculate infrared spectra with MLatom using a commonly used way: calculate frequencies and infrared (IR) intensities of static molecules.
Please check out how to calculate vibrational spectra from molecular dynamics trajectories.
Citation
When you use the features described here, please cite:
Yi-Fan Hou, Cheng Wang, Pavlo O. Dral. Accurate and Affordable Simulation of Molecular Infrared Spectra with AIQM Models. J. Phys. Chem. A 2025, ASAP. DOI: 10.1021/acs.jpca.5c00146. Preprint on ChemRxiv: https://doi.org/10.26434/chemrxiv-2024-r5prt (2024-11-08).
Theory
In the harmonic approximation, the frequencies of a stationary point can be calculated by diagonalizing its Hessian matrix (second derivatives of energy with respect to coordinates). The corresponding IR intensities need dipole derivatives (first derivatives of dipole moment with respect to coordinates):
\(I_{n} = \sum \limits _{t=1} ^{3} (\frac {\partial p_{t}} {\partial Q_{n}}) ^{2}\) ,
where \(t\) corresponds to each Cartesian dimension, \(p\) is the dipole moment of the molecule, \(Q_{n}\) is the \(n\) th normal coordinate. Usually we get dipole derivatives with respect to each atomic coordinate, so the above equation can be rewritten as:
\(I_{n} = \sum \limits _{t=1} ^{3} (\sum \limits _{x} \frac {\partial p_{t}} {\partial q_{x}} \frac {\partial q_{x}} {\partial Q_{n}}) ^{2}\) ,
where \(q_{x}\) is the component of atomic coordinates.
Example
We present a simple example of calculating harmonic frequencies and IR intensities of ethanol, whose geometry is provided here ethanol_init.xyz. We need to optimize its geometry first and then calculate its frequencies as frequency calculations are only reasonable for stationary points (e.g., local minima and transition states).
Example with input files
First we need to optimize the geometry of ethanol. Below is the input file for this task. You can refer to geometry optimization tutorial for more details.
geomopt
B3LYP/6-31G*
xyzfile=ethanol_init.xyz
optxyz=ethanol_opt.xyz
After getting the geometry of the optimized ethanol (ethanol_opt.xyz) we can calculate its frequencies and IR intenities. For IR intenisties, the keyword ir is needed. You can refer to frequencies tutorial for information of other keywords.
ir
B3LYP/6-31G*
xyzfile=ethanol_opt.xyz
You should be able to find in the MLatom output the vibrational analysis and IR intensities.
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants IR intensities
(cm^-1) (AMU) (mDyne/A) (km/mol)
1 248.4747 1.1524 0.0419 32.3833
2 301.0663 1.0696 0.0571 99.1615
3 417.5299 2.6381 0.2710 11.2144
4 830.0330 1.0800 0.4384 0.0152
5 912.5116 2.2005 1.0796 9.7442
6 1042.3161 2.1512 1.3770 56.5354
7 1125.4816 2.2612 1.6876 23.7752
8 1193.9827 1.4953 1.2559 4.7774
9 1290.5694 1.2508 1.2275 87.6053
10 1310.5341 1.1222 1.1356 0.0903
11 1427.9291 1.2321 1.4802 0.4122
12 1480.9785 1.4590 1.8854 16.5845
13 1513.8300 1.0411 1.4058 4.4003
14 1532.0798 1.0458 1.4464 3.0858
15 1561.1233 1.0910 1.5666 2.6189
16 2980.0716 1.0552 5.5213 73.9417
17 3003.5497 1.1085 5.8921 78.0445
18 3058.1727 1.0352 5.7041 15.7969
19 3129.0297 1.1026 6.3606 31.5280
20 3133.5337 1.1038 6.3856 34.2645
21 3746.8800 1.0665 8.8215 10.4090
On the XACS platform, we provide our own implementation of calculating IR intensities. IR intensities are available when using PySCF for HF and DFT, which are based on the pyscf.prop module.
There are several ways to visualize the vibrations:
At the end of the calculations, MLatom will dump the
freq_gaussian{mol_index}.logfile, which uses Gaussian-style output format for frequencies. You can open it in ChemCraft.You can also dump text output file using Gaussian-style format for frequencies from the Python API:
mymolecule.dump(filename='gauss.log', format='gaussian').
Linear scaling of frequencies
It is a common practice to scale the frequencies linearly to make the calculated IR spectra match better with the experiments. We provide a keyword scaling to deal with it. When the user specifies its value, MLatom will output both unscaled and scaled frequencies. For AIQM1 and AIQM2 methods, the scaling is activated by default and their frequencies will be scaled by 0.957 and 0.962, respectively.
Now we will show an example of calculating the IR spectrum of ethanol with AIQM2. First, we need to optimize its geometry, which is provided here: ethanol_init.xyz.
The input file should look like:
geomopt
AIQM2
xyzfile=ethanol_init.xyz
optprog=ase
optxyz=ethanol_opt.xyz
Then we can calculate the frequencies and IR intensities:
ir
AIQM2
xyzfile=ethanol_opt.xyz
As scaling factor is applied by default for AIQM2, the above input file is equivalent to:
ir
AIQM2
xyzfile=ethanol_opt.xyz
scaling=0.962
You will find the following things in your output file:
==============================================================================
Vibration analysis for molecule 1
==============================================================================
Multiplicity: 1
This is a nonlinear molecule
Mode Frequencies Reduced masses Force Constants IR intensities
(cm^-1) (AMU) (mDyne/A) (km/mol)
1 233.0042 1.1460 0.0367 24.3148
2 284.2353 1.0716 0.0510 136.0188
3 420.9724 2.6476 0.2765 12.7217
4 825.6527 1.0829 0.4350 2.5855
5 913.6812 2.2319 1.0978 6.3449
6 1055.6910 2.1786 1.4305 39.9877
7 1115.0663 2.2802 1.6704 37.0315
8 1187.5787 1.5239 1.2663 13.1134
9 1280.7458 1.2320 1.1907 103.6870
10 1313.8324 1.1081 1.1270 0.9832
11 1409.4878 1.2342 1.4447 2.0347
12 1453.2597 1.4434 1.7961 18.2637
13 1493.2140 1.0392 1.3652 1.7485
14 1512.6421 1.0510 1.4169 4.6961
15 1531.5126 1.0895 1.5056 1.8487
16 3009.1176 1.0544 5.6250 69.1294
17 3044.8103 1.0348 5.6522 14.5552
18 3047.3488 1.1081 6.0628 74.9335
19 3121.0591 1.1002 6.3146 38.3573
20 3130.4621 1.1021 6.3636 33.9217
21 3847.8329 1.0663 9.3014 6.0155
Scale frequencies linearly
Scaling factor: 0.962
Mode Frequencies Reduced masses Force Constants IR intensities
(cm^-1) (AMU) (mDyne/A) (km/mol)
1 224.1500 1.1460 0.0367 24.3148
2 273.4344 1.0716 0.0510 136.0188
3 404.9755 2.6476 0.2765 12.7217
4 794.2779 1.0829 0.4350 2.5855
5 878.9613 2.2319 1.0978 6.3449
6 1015.5748 2.1786 1.4305 39.9877
7 1072.6938 2.2802 1.6704 37.0315
8 1142.4507 1.5239 1.2663 13.1134
9 1232.0775 1.2320 1.1907 103.6870
10 1263.9068 1.1081 1.1270 0.9832
11 1355.9273 1.2342 1.4447 2.0347
12 1398.0358 1.4434 1.7961 18.2637
13 1436.4719 1.0392 1.3652 1.7485
14 1455.1617 1.0510 1.4169 4.6961
15 1473.3151 1.0895 1.5056 1.8487
16 2894.7711 1.0544 5.6250 69.1294
17 2929.1075 1.0348 5.6522 14.5552
18 2931.5496 1.1081 6.0628 74.9335
19 3002.4588 1.1002 6.3146 38.3573
20 3011.5045 1.1021 6.3636 33.9217
21 3701.6153 1.0663 9.3014 6.0155
You will see two blocks of vibration analysis. Unscaled frequencies come first, then the scaled frequencies. MLatom will also dump the scaled_freq_gaussian{mol_index}.log file that saves the scaled frequencies in Gaussian-style output format.
Example with Python
Infrared spectra can be calculated and analyzed with publication-ready quality figures. In a Jupyter notebook, the spectra and vibrations can be visualized at the same time, i.e., we offer a slider to choose which normal model to show and at the same time highlight the corresponding IR band in the spectrum. You can download the Jupyter notebook and required files ethanol_init.xyz (initial geometry guess for ethanol molecule) and exp.txt (experimental spectrum for comparison, from the NIST database):
import mlatom as ml
from mlatom.spectra import ir
# Visualize the calculated spectrum
mol = ml.molecule.from_xyz_file('ethanol_init.xyz')
aiqm2 = ml.methods(method='AIQM2')
_ = ml.optimize_geometry(molecule=mol, model=aiqm2, program='geometric')
_ = ml.freq(molecule=mol, model=aiqm2, ir=True, program='pyscf')
mol.dump('ethanol_freq.json', format='json')
#mol.load('ethanol_freq.json',format='json')
print('Optimized geometry')
print(mol.get_xyz_string())
print()
print("Mode Frequencies Reduced masses Force Constants")
print(" (cm^-1) (AMU) (mDyne/A)")
for ii in range(len(mol.frequencies)):
print("%d %13.4f %13.4f %13.4f"%(ii,mol.frequencies[ii],mol.reduced_masses[ii],mol.force_constants[ii]))
Optimized geometry
9
C -3.8757888541437 -1.2622761408403 0.0093531161968
H -3.5147169293322 -2.2897467272791 -0.0115720789819
H -3.5012865398627 -0.7487568663630 -0.8753184958339
H -4.9661466783656 -1.2732515032080 -0.0275690943305
C -3.3998580115024 -0.5638274393319 1.2625112591570
H -3.7703918347260 0.4660630794194 1.2842713105245
H -3.7793994352488 -1.0810871776526 2.1495613143141
O -1.9771043502805 -0.5811675365947 1.2466673219004
H -1.6564204641551 -0.1380559319389 2.0329922668624
Mode Frequencies Reduced masses Force Constants
(cm^-1) (AMU) (mDyne/A)
0 232.1385 1.1451 0.0364
1 285.0563 1.0721 0.0513
2 420.9774 2.6478 0.2765
3 824.8721 1.0830 0.4342
4 913.3403 2.2298 1.0959
5 1055.4557 2.1873 1.4356
6 1114.9687 2.2775 1.6681
7 1187.2235 1.5246 1.2661
8 1280.8234 1.2307 1.1895
9 1313.2426 1.1081 1.1260
10 1409.6175 1.2341 1.4447
11 1453.3303 1.4433 1.7961
12 1493.2445 1.0392 1.3652
13 1512.6419 1.0510 1.4169
14 1531.6462 1.0895 1.5060
15 3011.6033 1.0544 5.6346
16 3044.6061 1.0349 5.6519
17 3049.1645 1.1080 6.0694
18 3120.6948 1.1003 6.3132
19 3129.9131 1.1021 6.3613
20 3847.5977 1.0663 9.3002
spectrum = ir.lorentzian(molecule=mol,fwhm=30)
spectrum.plot('IR.png')
linespec = ir(frequencies=mol.frequencies,infrared_intensities=mol.infrared_intensities)
# Check vibration analysis
# Check thermochemistry results
print(f"Zero-point vibrational energy: {mol.ZPE} Hartree")
print(f"Enthalpy at 298 K: {mol.H} Hartree")
print(f"Gibbs Free energy at 298 K: {mol.G} Hartree")
Zero-point vibrational energy: 0.08026476279214442 Hartree Enthalpy at 298 K: -154.80626341289735 Hartree Gibbs Free energy at 298 K: -154.83694374633512 Hartree
# Visualize the experimental spectrum
exp = ml.spectra.ir.load('exp.txt',format='txt')
exp.plot('exp.png')
# Plot both spectra together
ml.spectra.plot_ir(spectra=[exp],linespectra=[linespec],
molecule=mol,
normalize='average',
labels=['experiment','theory'],
lorentzian=True,
scaling_factor=0.958,
plotstart=500,plotend=4000)
Theoretical spectrum is scaled by 0.958
# Visualize vibrations and highlight corresponding peaks
exp = ml.spectra.ir.load('exp.txt',format='txt')
ml.spectra.plot_ir(spectra=[exp],linespectra=[linespec],
molecule=mol,
normalize='average',
labels=['experiment','theory'],
lorentzian=True,
band_width_slider=True,
peak_highlight_slider=True,
scaling_factor=0.958,
plotstart=500,plotend=4000)
interactive(children=(FloatSlider(value=30.0, description='width/cm^-1:', max=50.0, min=10.0, readout_format='…
You can directly view the vibrations in
Jupyter notebookby usingmymolecule.view(normal_mode=...), as shown below for an examplemolvibr.jsonfile:
import mlatom as ml
# let's load the molecule with the required properties
molvibr = ml.molecule()
molvibr.load(filename='molvibr.json')
# choose the normal mode to visualize
normal_mode = 1
# let's view it
molvibr.view(normal_mode=normal_mode)
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
# we can also check the information for this normal mode
print(f'Frequency: {molvibr.frequencies[normal_mode]:.2f} cm^-1')
print(f'Normal mode in Angstrom')
for atom in molvibr.atoms:
disp = atom.normal_modes[normal_mode]
print(f" {disp[0]:8.3f} {disp[1]:8.3f} {disp[2]:8.3f}")
Frequency: 288.70 cm^-1
Normal mode in Angstrom
-0.000 0.000 -0.025
-0.085 0.168 -0.141
0.085 -0.168 -0.141
0.000 0.000 0.184
-0.000 -0.000 -0.003
0.018 -0.021 -0.014
-0.018 0.021 -0.014
0.000 0.000 0.085
-0.000 -0.000 -0.890
# the non-displaced geometry for reference
print(molvibr.get_xyz_string())
9 C -1.2121068000000 -0.2252488500000 0.0000164100000 H -1.2701429600000 -0.8599818400000 -0.8855315500000 H -1.2701264500000 -0.8599834500000 0.8855643400000 H -2.0710636400000 0.4504289000000 0.0000252100000 C 0.0804149500000 0.5556635800000 0.0000041500000 H 0.1395635900000 1.1970970800000 -0.8856250600000 H 0.1395848700000 1.1970854800000 0.8856411200000 O 1.1428329300000 -0.3971737900000 -0.0000106600000 H 1.9796722200000 0.0702558200000 -0.0001121300000
# non-displaced geometry
molvibr.view()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
Any questions or suggestions?
If you have further questions, criticism, and suggestions, we would be happy to receive them in English or Chinese via email, Slack (preferred), or WeChat (please send an email to request to add you to the XACS user support group).