13
$\begingroup$

Does anyone know how to export a calculated Raman spectrum (my calculation was performed in Q-Chem) as a data-file, which can be used to plot the profile of a Gaussian-shape spectrum (in origin, excel, etc.), rathe than centroid data points?

$\endgroup$
1
  • 3
    $\begingroup$ There's also QMForge which runs on cclib. $\endgroup$ May 19 '20 at 19:33
14
$\begingroup$

I wrote the following simple Python function that allows one to use any distribution from scipy.stats for peak broadening (the scale parameter determines how broad peaks will be):

import numpy as np
from scipy.stats import cauchy, norm


def broaden_spectrum(x, x0, y0, distribution="gaussian", scale=1.0, 
                     fit_points=True, *args, **kwargs):
    if distribution in {"gaussian", "norm"}:
        distribution = norm
    elif distribution in {"lorentzian", "cauchy"}:
        distribution = cauchy

    s = np.sum([yp * distribution.pdf(x, xp, scale=scale, *args, **kwargs)
                for xp, yp in zip(x0, y0)], axis=0)

    if fit_points:
        s_max = np.max(s)
        if s_max == 0.0:
            s_max = 1.0
        return s * np.max(y0) / s_max
    return s

As Geoff Hutchison mentioned, you can use cclib to read your logfile:

import cclib


data = cclib.ccopen("data/2-propanol.out").parse()
x = np.linspace(data.vibfreqs.min() - 100.,
                data.vibfreqs.max() + 100., num=1000)
y = broaden_spectrum(x, data.vibfreqs, data.vibirs, scale=40.)

Plotting is easy:

import matplotlib.pyplot as plt


plt.plot(x, y)
plt.gca().invert_xaxis()
plt.gca().invert_yaxis()
plt.ylabel("Absorbance (arb. unit.)")
plt.xlabel("Wavenumber (cm$^{-1}$)")
plt.title("Predicted IR spectrum of 2-propanol at revPBE-D4-gCP/def2-SVP")
plt.savefig("ir-spectrum.png")

Predicted IR spectrum of 2-propanol at revPBE-D4-gCP/def2-SVP

This is easily customizable and the broaden_spectrum function works for all kinds of spectra.

$\endgroup$
1
  • 3
    $\begingroup$ +1 Great answer! Haven't seen you for a while. Welcome back! Thank you for your contributions! Notice also that we now have syntax coloring! $\endgroup$ May 20 '20 at 2:14
11
$\begingroup$

Many thanks to Felipe, the code is very useful! I know this is just a simple editing, but this is how I adjusted part of the code for Raman spectra:

import cclib


data = cclib.ccopen("molecule.out").parse()
x = np.linspace(data.vibfreqs.min() - 100.,
                data.vibfreqs.max() + 100., num=1000)
y = broaden_spectrum(x, data.vibfreqs, data.vibramans, scale=10.)

And also this will save a .csv which can be used for further data processing:

import matplotlib.pyplot as plt
import pandas as pd


plt.plot(x, y)
plt.ylabel("Intensity (a.u.)")
plt.xlabel("Raman Shift (cm$^{-1}$)")
plt.title("Predicted Raman spectrum of molecule at B3LYP/6-311++G**")
plt.savefig("Raman_gauss.pdf", dpi=300, bbox_inches='tight')


dataframe = pd.DataFrame({'x': x,
                          'y':y})
dataframe.to_csv("Raman_spectrum.csv", index=False)

I'm sure this has multiple options for customization, but I though it might be useful.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.