Plotting of Spectra

This tutorial explains how to produce various types of spectra through the Dalton Project platform.

Note, that the quantum chemistry methods used in this tutorial may not be suitable for the particular problem that you want to solve. Moreover, the basis sets used here are minimal in order to allow you to progress fast through the tutorial and should under no circumstances be used in production calculations.

UV/vis Spectrum

To create a UV/vis spectrum, we need to obtain excitation energies and oscillator strengths. First, we will need to specify our molecule and basis set:

import matplotlib.pyplot as plt

import daltonproject as dp

molecule = dp.Molecule(input_file='water.xyz')
basis = dp.Basis(basis='cc-pVDZ')

The excitation energies and oscillator strengths can be obtained with any method of choice, for example HF, CASSCF, or CC, but for this example, we will obtain them using TD-DFT with the CAM-B3LYP functional. In this case, we include the six lowest excitation energies and associated oscillator strengths.

dft = dp.QCMethod('DFT', 'CAMB3LYP')
prop = dp.Property(excitation_energies=True)
prop.excitation_energies(states=6)
result = dp.dalton.compute(molecule, basis, dft, prop)
print('Excitation energies =', result.excitation_energies)
# Excitation energies = [ 7.708981 10.457775 11.631572 11.928512 14.928328 15.733272]

print('Oscillator strengths =', result.oscillator_strengths)
# Oscillator strengths = [0.0255 0.     0.1324 0.0883 0.0787 0.171 ]

To visualize the resulting spectrum, we can use the built-in plotting functionality from the spectrum module.

ax = dp.spectrum.plot_one_photon_spectrum(result, color='k')
plt.savefig('1pa.svg')

And the resulting spectrum looks like:

../../_images/1pa.svg

The plotting function returns a new matplotlib Axes object by default. An existing Axes object can also be passed to the plotting function with the ax= keyword argument. The frequencies at which the spectrum is by default selected based on the excitation energies, but can also be specified manually with the frequencies= keyword argument. Optional standard plotting arguments, such as color or line-style can be passed as additional keyword arguments, for example with color='k' to set a black line color.

The complete Python script for this example is:

import matplotlib.pyplot as plt

import daltonproject as dp

molecule = dp.Molecule(input_file='water.xyz')
basis = dp.Basis(basis='cc-pVDZ')

dft = dp.QCMethod('DFT', 'CAMB3LYP')
prop = dp.Property(excitation_energies=True)
prop.excitation_energies(states=6)
result = dp.dalton.compute(molecule, basis, dft, prop)
print('Excitation energies =', result.excitation_energies)
# Excitation energies = [ 7.708981 10.457775 11.631572 11.928512 14.928328 15.733272]

print('Oscillator strengths =', result.oscillator_strengths)
# Oscillator strengths = [0.0255 0.     0.1324 0.0883 0.0787 0.171 ]

ax = dp.spectrum.plot_one_photon_spectrum(result, color='k')
plt.savefig('1pa.svg')

The plot is generated by making a convolution of the excitation energies by the following equation:

\[\varepsilon(\omega) = \frac{e^2\pi^2N_A}{\ln(10)2\pi\epsilon_0nm_ec}\sum_i \frac{\omega f_i}{\omega_i}g(\omega,\omega_i,\gamma_i)\]

here \(N_a\) is Avogadro’s constant, \(e\) is the elementary charge, \(m_e\) is the electron mass, \(c\) is the speed of light, \(\epsilon_0\) is the vacuum permittivity, \(n\) is the refractive index (approximated to be one), \(f_i\) is the calculated oscillator strength of the \(i\) th transition, and \(\omega_i\) is the corresponding transition angular frequency.

It can be noted that another form is seen often in the literature:

\[\varepsilon(\omega) = \frac{e^2\pi^2N_A}{\ln(10)2\pi\epsilon_0nm_ec}\sum_i f_ig(\omega,\omega_i,\gamma_i)\]

This form just corresponds to approximating \(\frac{\omega f_i} {\omega_i} \approx f_i\).

Two-Photon Absorption Spectrum

To create a two-photon absorption spectrum, we need to calculate excitation energies and two-photon strengths. First, we specify the molecule and basis set:

import matplotlib.pyplot as plt

import daltonproject as dp

molecule = dp.Molecule(input_file='water.xyz')
basis = dp.Basis(basis='cc-pVDZ')

The excitation energies and two-photon strengths can be calculated with many different methods, and for this example we will use TD-DFT with the CAM-B3LYP functional. In this case, we include the three lowest excitation energies and associated oscillator strengths.

dft = dp.QCMethod('DFT', 'CAMB3LYP')
prop = dp.Property(two_photon_absorption=True)
prop.two_photon_absorption(states=3)
result = dp.dalton.compute(molecule, basis, dft, prop)
print('Excitation energies =', result.excitation_energies)
# Excitation energies = [ 7.70898022 10.45777647 11.63157226]

print('Two-photon cross-sections =', result.two_photon_cross_sections)
# Two-photon cross-sections = [0.139 0.776 0.585]

To visualize the resulting spectrum, we can use the built-in plotting functionality from the spectrum module.

ax = dp.spectrum.plot_two_photon_spectrum(result, color='k')
plt.savefig('2pa.svg')

And the resulting spectrum looks like:

../../_images/2pa.svg

The complete Python script for this example is:

import matplotlib.pyplot as plt

import daltonproject as dp

molecule = dp.Molecule(input_file='water.xyz')
basis = dp.Basis(basis='cc-pVDZ')

dft = dp.QCMethod('DFT', 'CAMB3LYP')
prop = dp.Property(two_photon_absorption=True)
prop.two_photon_absorption(states=3)
result = dp.dalton.compute(molecule, basis, dft, prop)
print('Excitation energies =', result.excitation_energies)
# Excitation energies = [ 7.70898022 10.45777647 11.63157226]

print('Two-photon cross-sections =', result.two_photon_cross_sections)
# Two-photon cross-sections = [0.139 0.776 0.585]

ax = dp.spectrum.plot_two_photon_spectrum(result, color='k')
plt.savefig('2pa.svg')

X-ray Absorption Spectroscopy using Coupled-Cluster methods

Using CC methods, you have the option to employ the core-valence separation approximation to produce X-ray absorption spectra. Suppose we want to obtain such a spectrum for acrolein using the CC2 method:

import matplotlib.pyplot as plt

import daltonproject as dp

acrolein = dp.Molecule(input_file='acrolein.xyz')
basis = dp.Basis('STO-3G')

cc = dp.QCMethod('CC2')
exc = dp.Property(excitation_energies={'states': 6, 'cvseparation': [2, 3, 4]})

Here we choose to include six transitions, restricting the excitation channels to specific core orbitals specified in the cvseparation list.

The cvseparation list specifies the active orbitals. In general, they are specified for each irreducible representation (irrep). Suppose we have cvseparation=[[1, 2, 3, 4], [2, 3], [3, 4], [1, 2, 3, 4, 5, 6]]. This would read: from the first irrep use orbitals number 1, 2, 3, and 4; from the second use orbitals number 2 and 3; from the third use orbitals number 3 and 4; and from the fourth irrep use orbitals number 1 to 6. Note that this requires prior analysis of the orbitals using the same symmetry as in excitation energy calculations in order to identify the orbitals of interest.

If we want to include excitations from the lowest lying s-like orbitals of the carbons in acrolein, we thus enter cvseparation=[2, 3, 4], as shown in the example above, skipping the first orbital which belongs to the oxygen.

Dalton Project will by default use any available CPU cores, and since Dalton is only uses MPI parallelization, it will set the number of MPI processes equal to number of CPU cores. However, CC in Dalton is not parallelized so we need to override the default settings:

settings = dp.ComputeSettings(mpi_num_procs=1)

Now we can run the calculation:

result = dp.dalton.compute(
    molecule=acrolein,
    basis=basis,
    qc_method=cc,
    properties=exc,
    compute_settings=settings,
)

and plot it:

print('Excitation energies =', result.excitation_energies)
# Excitation energies = [288.23439 288.37014 289.96145 293.82649 294.58115 296.95772]
../../_images/cvs_cc.svg

The entire script to produce the plots is given below.

import matplotlib.pyplot as plt

import daltonproject as dp

acrolein = dp.Molecule(input_file='acrolein.xyz')
basis = dp.Basis('STO-3G')

cc = dp.QCMethod('CC2')
exc = dp.Property(excitation_energies={'states': 6, 'cvseparation': [2, 3, 4]})

settings = dp.ComputeSettings(mpi_num_procs=1)

result = dp.dalton.compute(
    molecule=acrolein,
    basis=basis,
    qc_method=cc,
    properties=exc,
    compute_settings=settings,
)

print('Excitation energies =', result.excitation_energies)
# Excitation energies = [288.23439 288.37014 289.96145 293.82649 294.58115 296.95772]

Vibrational IR and/or Raman Spectrum

To simulate vibrational IR and Raman spectra, we need vibrational frequencies and corresponding IR/Raman intensities. The vibrational frequencies are obtained from the molecular Hessian, i.e., the matrix containing second- derivatives with respect to nuclear displacements, while the dipole and polarizability gradients are needed to compute IR and Raman intensities, respectively. This is a two-step procedure where the basic properties are calculated first and subsequently used to derive the final quantities needed to plot the spectra. To proceed, we specify the molecule and level of theory:

import matplotlib.pyplot as plt

import daltonproject as dp

hf = dp.QCMethod(qc_method='HF')
basis = dp.Basis(basis='pcseg-0')
h2o2 = dp.Molecule(input_file='H2O2.xyz')

The properties should be evaluated at the equilibrium geometry. Therefore, we start by optimizing the geometry using LSDalton:

geo_opt = dp.Property(geometry_optimization=True)
opt_output = dp.lsdalton.compute(h2o2, basis, hf, geo_opt)
h2o2.coordinates = opt_output.final_geometry

Using the optimized geometry, we can proceed to calculate the Hessian, needed to derive vibrational frequencies, and dipole and polarizability gradients, from which the IR and Raman intensities are obtained. Raman intensities also depend on the frequency of the incident light. We may choose to compute the polarizability at this frequency (specified in hartree). The default is to compute static polarizability (i.e., a frequency of 0.0 hartree). Here we assume a light source with a wavelength of 633 nm corresponding to 0.07198 hartree:

props = dp.Property(hessian=True, dipole_gradients=True, polarizability_gradients={'frequencies': [0.07198]})
prop_output = dp.lsdalton.compute(h2o2, basis, hf, props)

After the properties have been calculated, a vibrational analysis can be performed, which involves mass-weighting and diagonalizing the Hessian as well as projecting out translational and rotational degrees of freedom. Through this, we obtain the frequencies and the transformation matrix, allowing us to transform quantities from Cartesian coordinates to normal coordinates. If supplied, the dipole and polarizability gradients are immediately transformed and used to compute IR and Raman intensities, respectively:

vib_ana = dp.vibrational_analysis(molecule=h2o2,
                                  hessian=prop_output.hessian,
                                  dipole_gradients=prop_output.dipole_gradients,
                                  polarizability_gradients=prop_output.polarizability_gradients)

The results can be used to plot spectra through the spectrum module. Since it uses Matplotlib, we can take advantage of its great flexibility. Here we are interested in comparing the IR and Raman spectra and therefore want to plot both in the same spectrum, which can be done as follows:

fig, ax1 = plt.subplots()
ax1 = dp.spectrum.plot_ir_spectrum(vib_ana, ax=ax1, color='blue', label='IR')
ax2 = ax1.twinx()
ax2 = dp.spectrum.plot_raman_spectrum(vib_ana, ax=ax2, color='red', label='Raman')
fig.legend(loc='center')
fig.savefig('ir_raman.svg')
../../_images/ir_raman.svg

The complete Python script for this example is:

import matplotlib.pyplot as plt

import daltonproject as dp

hf = dp.QCMethod(qc_method='HF')
basis = dp.Basis(basis='pcseg-0')
h2o2 = dp.Molecule(input_file='H2O2.xyz')

geo_opt = dp.Property(geometry_optimization=True)
opt_output = dp.lsdalton.compute(h2o2, basis, hf, geo_opt)
h2o2.coordinates = opt_output.final_geometry

props = dp.Property(hessian=True, dipole_gradients=True, polarizability_gradients={'frequencies': [0.07198]})
prop_output = dp.lsdalton.compute(h2o2, basis, hf, props)

vib_ana = dp.vibrational_analysis(molecule=h2o2,
                                  hessian=prop_output.hessian,
                                  dipole_gradients=prop_output.dipole_gradients,
                                  polarizability_gradients=prop_output.polarizability_gradients)

fig, ax1 = plt.subplots()
ax1 = dp.spectrum.plot_ir_spectrum(vib_ana, ax=ax1, color='blue', label='IR')
ax2 = ax1.twinx()
ax2 = dp.spectrum.plot_raman_spectrum(vib_ana, ax=ax2, color='red', label='Raman')
fig.legend(loc='center')
fig.savefig('ir_raman.svg')