Spectrum Function (TDDFT) Giving Error

I’m trying to reproduce the example here

import numpy as np

import psi4

from psi4.driver.procrouting.response.scf_response import tdscf_excitations
from psi4.driver.p4util import spectrum

psi4.core.set_output_file("moxy.out")

moxy = psi4.geometry("""0 1
C  0.152133 -0.035800  0.485797
C -1.039475  0.615938 -0.061249
C  1.507144  0.097806 -0.148460
O -0.828215 -0.788248 -0.239431
H  0.153725 -0.249258  1.552136
H -1.863178  0.881921  0.593333
H -0.949807  1.214210 -0.962771
H  2.076806 -0.826189 -0.036671
H  2.074465  0.901788  0.325106
H  1.414895  0.315852 -1.212218
""", name="(S)-methyloxirane")

psi4.set_options({
    'save_jk': True,
})

e, wfn = psi4.energy("HF/cc-pvdz", return_wfn=True, molecule=moxy)
res = tdscf_excitations(wfn, states=8, triplets="also")

# get poles and residues to plot OPA and ECD spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENT"])**2 for r in res]
ecd_residues = [r["LENGTH-GAUGE ROTATORY STRENGTH"] for r in res]

opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")

But get the error:

Traceback (most recent call last):
  File "/data/Apps/anaconda3/envs/psi4env/bin/psi4", line 333, in <module>
    exec(content)
  File "<string>", line 42, in <module>
  File "<string>", line 42, in <listcomp>

KeyError: 'LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENT'


Printing out the relevant lines from the Psithon --> Python processed input file:
        'save_jk': True,
    })
    e, wfn = psi4.energy("HF/cc-pvdz", return_wfn=True, molecule=moxy)
    res = tdscf_excitations(wfn, states=8, triplets="also")
    poles = [r["EXCITATION ENERGY"] for r in res]
--> opa_residues = [np.linalg.norm(r["LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENT"])**2 for r in res]
    ecd_residues = [r["LENGTH-GAUGE ROTATORY STRENGTH"] for r in res]
    opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
    ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")

!---------------------------------------------------!
!                                                   !
!  'LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENT' !
!                                                   !
!---------------------------------------------------!

It looks like some labels were changed since that example was written. Try changing these:

LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENTELECTRIC DIPOLE TRANSITION MOMENT (LEN)

LENGTH-GAUGE ROTATORY STRENGTHROTATORY STRENGTH (LEN)

That worked – thanks.

Although I am not getting the error anymore when running the commands, I don’t see the output produced from:

# get poles and residues to plot OPA and ECD spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["ELECTRIC DIPOLE TRANSITION MOMENT (LEN)"])**2 for r in res]
ecd_residues = [r["ROTATORY STRENGTH (LEN)"] for r in res]

opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")

The last output I get are from ‘Contributing excitations and de-excitations’ → For example:

Contributing excitations and de-excitations
Only contributions with coefficients > 1.00e-01 will be printed:

Excited State    1 (1 A):   0.24630 au   185.00 nm f = 0.0434
  Sums of squares: Xssq =  1.000490e+00; Yssq =  4.902097e-04; Xssq - Yssq =  1.000000e+00
     5 ->  6    0.991925 (98.391%)
     5 ->  8    0.117839 ( 1.389%)

My full input file (using water as example):

molecule mol {
0 1
O 0.00000000 0.00000000 0.00000000
H 0.00000000 -0.79069000 0.61221700
H 0.00000000 0.79069000 0.61221700

symmetry c1
}

set scf_type df
set basis aug-cc-pVTZ
set reference rks
set s_tolerance 1e-9
set_num_threads(62)

set save_jk true

import numpy as np
import psi4
from psi4.driver.procrouting.response.scf_response import tdscf_excitations
from psi4.driver.p4util import spectrum

e, wfn = energy('b3lyp-d3bj', return_wfn=True)
res = tdscf_excitations(wfn, states=20)

# get poles and residues to plot OPA and ECD spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["ELECTRIC DIPOLE TRANSITION MOMENT (LEN)"])**2 for r in res]
ecd_residues = [r["ROTATORY STRENGTH (LEN)"] for r in res]

opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")

mol.print_out()
print_variables()

The spectrum function returns a dictionary with different peak types as keys (“convolution” for broadened, “sticks” for sticks), and each peak type has an “x” and “y” key. In the TDSCF example you linked, they say they used the Altiar library to plot the spectra they show. I’m not familiar with this library, but if you want to use matplotlib in your input file to save an image of an ECD spectrum, it would look something like:

import matplotlib.pyplot as plt
ecd_spectrum = spectrum(
    poles=poles,
    residues=ecd_residues,
    kind="ECD",
    gamma=0.01,
    out_units="nm",
)
plt.plot(ecd_spectrum["convolution"]["x"], ecd_spectrum["convolution"]["y"])
plt.savefig("ecd.png")

You could also write the x and y data in CSV files so you can plot them in Excel or something.

Although you might need the flexibility of using your own code to process spectra, I’ll note that there are software packages that can plot the data in the output file without having to include anything after the tdscf_excitations function is called (though you should still be able to put things there). I think Avogadro can do this, but I also help develop some other stuff for it called AaronTools and SEQCROW.

The AaronTools python module has a command line utility for plotting UV/vis spectra called plotUVVis.py. For example:

plotUVVis.py test_tdscf.out --plot-type uv-vis-velocity --ranges 150-450 --output test_uvvis.png

will plot the spectrum in the 150-450nm region. If you don’t use --output, it will show you the spectrum instead. Run with just --help or look on the GitHub wiki to see all the options the script has. There is an installation guide elsewhere on the wiki, but simplest option for most is just running python -m pip install AaronTools on the command line.

I also have a plugin for UCSF ChimeraX called SEQCROW. This plugin will allow you to open Psi4 (and other) output files in ChimeraX and plot spectra with the tool under ‘Tools → AaronTools → UV/Vis Spectrum’ on ribbon menu. If you aren’t Boltzmann-weighting the spectra, you will have to set the “energy for weighting” option to “electronic” and don’t add extra conformers.

You can install ChimeraX from their website. They have installers for Windows, MacOS, and various Linux distros. To install the SEQCROW plugin, go to ‘Tools → More Tools…’ on the ribbon menu, find SEQCROW and click install.

Thanks, this helps a lot. I’ll check out those packages.

Following up; some fundamental questions:

When printing both ECD and OPA →

poles = [r["EXCITATION ENERGY"] for r in res]

opa_residues = [np.linalg.norm(r["ELECTRIC DIPOLE TRANSITION MOMENT (LEN)"])**2 for r in res]
opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
plt.plot(opa_spectrum["convolution"]["x"], opa_spectrum["convolution"]["y"])
plt.savefig("opa.png")

ecd_residues = [r["ROTATORY STRENGTH (LEN)"] for r in res]
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")
plt.plot(ecd_spectrum["convolution"]["x"], ecd_spectrum["convolution"]["y"])
plt.savefig("ecd.png")

Q1: The second one (ecd.png) has both of the spectra in it - how do I write both, but keep them separate?

Q2: What are you using to write the CSV (do you import csv, for example). Could show an example?

thanks again!

Using plt.clf() between calling plt.plot will clear the figure so you don’t get both datasets on it.

For saving the data to a CSV file, you can use NumPy’s savetxt function:

import numpy as np
# get data arrays to save
data_x = ecd_spectrum["convolution"]["x"]
data_y = ecd_spectrum["convolution"]["y"]
# transpose because numpy treats (array_1, array_2) as a 2 x N matrix
# and I want a N x 2 matrix
np.savetxt("ecd.csv", np.transpose((data_x, data_y)), delimiter=",")

This worked – thanks again.

In attempt to write the raw data for later processing, I did:

data_x_opa_raw = opa_spectrum["sticks"]["x"]
data_y_opa_raw = opa_spectrum["sticks"]["y"]
# transpose because numpy treats (array_1, array_2) as a 2 x N matrix
# and I want a N x 2 matrix
np.savetxt("opa_raw.csv", np.transpose((data_x_opa_raw, data_y_opa_raw)), delimiter=",")

but get and error message:

KeyError: 'x'

Printing out the relevant lines from the Psithon --> Python processed input file:
    plt.plot(opa_spectrum["convolution"]["x"], opa_spectrum["convolution"]["y"])
    plt.savefig("opa.png")
    data_x_opa = opa_spectrum["convolution"]["x"]
    data_y_opa = opa_spectrum["convolution"]["y"]
    np.savetxt("opa.csv", np.transpose((data_x_opa, data_y_opa)), delimiter=",")
--> data_x_opa_raw = opa_spectrum["sticks"]["x"]
    data_y_opa_raw = opa_spectrum["sticks"]["y"]
    np.savetxt("opa_raw.csv", np.transpose((data_x_opa_raw, data_y_opa_raw)), delimiter=",")
    ecd_residues = [r["ROTATORY STRENGTH (LEN)"] for r in res]
    ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")
    plt.clf()

Is ‘sticks’ not the right term?

Sorry, I assumed sticks had the same “x” and “y” keys without checking. The keys for sticks are “poles” and “residues”. These are basically your inputs to spectrum, but converted to the units you request.

Also, if you want to plot them with matplotlib the hlines and vlines functions are better for this. They just plot vertical lines and horizontal lines, respectively. So you can use vlines to plot vertical lines for each excitation from 0 to the corresponding intensity. With hlines, you can plot the baseline x=0. Here’s an example for ECD:

plt.clf()
plt.vlines(
    ecd_spectrum["sticks"]["poles"],
    np.zeros(len(poles)),
    ecd_spectrum["sticks"]["residues"]
)
# scaling the min. and max. x range by a little bit
# because I think it looks better that way
plt.hlines(
    [0, 0],
    0.9 * min(ecd_spectrum["sticks"]["poles"]),
    1.1 * max(ecd_spectrum["sticks"]["poles"])
)
plt.savefig("ecd_sticks.png")

thanks! Now, If I wanted to write it to file with matplotlib?

You can use plt.savefig, like before.

There’s a lot you can do with matplotlib. I don’t use it often enough to warrant committing much of it to memory, but there are quite a few tutorials out there to help you get plots to look however you want them to. The matplotlib website is a good jumping-off point: Tutorials — Matplotlib 3.4.3 documentation

Cool – I used the format you provided before. And thanks for the link.