Skip to content

Broaden and plot IR spectra

Harmonic vibrational calculations produce discrete mode frequencies and intensities. Gaussian broadening turns those sticks into a smooth illustrative spectrum. This tutorial covers both Plot.gaussian_broadening() and the older Plot.ir() helper, including its current limitations.

The examples use invented frequencies so that they run without an external quantum-chemistry calculation. For real data, read a supported Gaussian output using the molecular file readers, then inspect the state's Freq and Ints arrays.

Make a broadened spectrum

gaussian_broadening(state, broaden=1, resolution=1) reads state.Freq and state.Ints, then sets state.IR to a two-column array containing wavenumber and broadened intensity. It returns None; it modifies the supplied state rather than the Plot.

import numpy as np
import matplotlib.pyplot as plt
from matmalib import State, Plot

state = State(
    Freq=np.array([980.0, 1050.0, 1080.0, 1250.0, 1450.0, 1650.0]),
    Ints=np.array([0.25, 1.0, 0.65, 0.4, 0.8, 0.55]),
)
plot = Plot()
plot.gaussian_broadening(state, broaden=12, resolution=1)
wavenumber, intensity = state.IR.T

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(wavenumber, intensity, color="#2563eb", label="Gaussian broadened")
ax.vlines(state.Freq, 0, state.Ints, color="#d97706", alpha=0.65, label="Modes")
ax.set(xlim=(900, 1750), xlabel="Wavenumber (cm⁻¹)",
       ylabel="Intensity (arbitrary units)", title="Six synthetic vibrational modes")
ax.legend()
fig.tight_layout()
print("IR array shape:", state.IR.shape)
print("Sampling interval:", state.IR[1, 0] - state.IR[0, 0], "cm⁻¹")

Each mode contributes

intensity_i × exp(-0.5 × ((wavenumber - frequency_i) / sigma)²).

Here broaden is the Gaussian standard deviation sigma in cm⁻¹, not its full width at half maximum. FWHM is approximately 2.355 * sigma. The function currently applies int(broaden), so use a positive integer; values between 0 and 1 become zero and cause invalid calculations. Intensities are peak coefficients; the Gaussian is not normalized to preserve area when its width changes.

The frequency grid always covers 0–4000 cm⁻¹, with int(4000 / resolution) + 1 points. resolution is approximately the spacing in cm⁻¹, not points per wavenumber. The endpoints are included, so a resolution that does not divide 4000 exactly has a slightly different actual step. Use a finite positive value, such as 1 or 2.

The output is not normalized, frequency-scaled, or converted to a transmittance spectrum. Frequencies and intensities should be numeric one-dimensional arrays of equal length; the implementation uses zip, so mismatched arrays silently drop the unmatched tail.

Explore peak merging and normalization

Broadening can make nearby modes appear as one feature. Compare widths while normalizing each spectrum to its own maximum:

import numpy as np
import matplotlib.pyplot as plt
from matmalib import Plot, State

state = State(Freq=np.array([1000.0, 1030.0, 1120.0]),
              Ints=np.array([1.0, 0.7, 0.8]))
plot = Plot()
fig, ax = plt.subplots(figsize=(8, 3))
for width in (4, 12, 25):
    plot.gaussian_broadening(state, broaden=width, resolution=1)
    x, y = state.IR.T
    ax.plot(x, y / y.max(), label=f"σ = {width} cm⁻¹")
ax.set(xlim=(950, 1200), xlabel="Wavenumber (cm⁻¹)",
       ylabel="Normalized intensity", title="Linewidth changes apparent peak structure")
ax.legend()
fig.tight_layout()

At the smallest width, the modes near 1000 and 1030 cm⁻¹ remain separate. At larger widths their envelopes overlap. Try moving the second mode from 1030 to 1060 to see how separation competes with linewidth.

Each call overwrites state.IR. Copy it with spectrum = state.IR.copy() if you need to retain several raw spectra. Normalizing each curve is convenient for shape comparison but removes differences in absolute intensity.

Work with a calculated state

Run this locally after replacing the path with your Gaussian output:

import matplotlib.pyplot as plt
from matmalib import State, Plot

state = State.from_file("calculation/input.log", software="gaussian")
print(state.Freq)
print(state.Ints)

plot = Plot()
plot.gaussian_broadening(state, broaden=10, resolution=1)
fig, ax = plt.subplots()
ax.plot(state.IR[:, 0], state.IR[:, 1])
ax.set(xlim=(900, 1800), xlabel="Wavenumber (cm⁻¹)", ylabel="IR intensity")
fig.savefig("calculated-spectrum.svg", bbox_inches="tight")

Check that the output actually contains a frequency calculation and that Freq and Ints exist before plotting. Negative frequencies and modes outside 0–4000 cm⁻¹ are not validated by the broadening helper; inspect them as part of interpreting the calculation.

If your analysis calls for a frequency scaling factor, apply it explicitly to a copied frequency array before broadening. There is no automatic scaling-factor parameter in either helper.

The legacy ir() convenience method

ir(state, xmin=900, xmax=1700, normal_modes=False) creates a normalized black spectrum with a colored filled region and an optional stick overlay. It unconditionally calls gaussian_broadening(state, resolution=0.001, broaden=2), overwriting any existing state.IR and creating 4,000,001 samples. This can be expensive, so the following example is intended for a local Python session rather than the browser.

The default Plot palette contains hex strings, but ir() expects mutable RGB lists and adds 0.25 to their channels. Supply a dark RGB palette immediately before calling it:

import numpy as np
import matplotlib.pyplot as plt
from matmalib import Plot, State

state = State(Freq=np.array([1000.0, 1250.0, 1500.0]),
              Ints=np.array([0.5, 1.0, 0.7]))
plot = Plot()
plot.set_colors([[0.1, 0.2, 0.4]])
plot.ir(state, xmin=900, xmax=1700, normal_modes=True)
fig, ax = plt.gcf(), plt.gca()
ax.set_xlabel("Wavenumber (cm⁻¹)")
ax.set_ylabel("Normalized intensity + offset")
fig.savefig("legacy-ir.png", dpi=150, bbox_inches="tight")

The current helper has these specific behaviors:

  • It mutates the input color lists in place and does not clip RGB values. Reset the palette before every call; initial channel values should stay at or below 0.75.
  • It normalizes by the maximum within the requested region, adds a 0.05 vertical offset, and hides y tick labels. A window with no nonzero spectrum cannot be normalized safely.
  • Its width is fixed to sigma=2 cm⁻¹; it offers no broadening, resolution, experimental-data, or frequency-scaling argument despite historical docstring wording.
  • It creates a Matplotlib figure but does not return it or set plot.fig and plot.ax. Use plt.gcf() and plt.gca() immediately afterward.

For control over linewidth, grid size, axis orientation, intensity scaling, or multiple spectra, use gaussian_broadening() with Matplotlib as in the interactive examples. See plot styling and export for the shared configuration tools.