# Spectrum calculation#

This section currently focus on the calculation of absorption spectra, with fluorescence and vibrational effects to added.

import copy

import gator
import matplotlib.pyplot as plt
import multipsi as mtp
import numpy as np
import veloxchem as vlx
from matplotlib import gridspec
from scipy.interpolate import interp1d
from veloxchem.lreigensolver import LinearResponseEigenSolver
from veloxchem.rsplinabscross import LinearAbsorptionCrossSection

# for vlx
silent_ostream = vlx.OutputStream(None)
from mpi4py import MPI

comm = MPI.COMM_WORLD

# au to eV conversion factor
au2ev = 27.211386

def lorentzian(x, y, xmin, xmax, xstep, gamma):
"""

Call: xi,yi = lorentzian(energies, intensities, start energy, end energy, energy step, gamma)
"""
xi = np.arange(xmin, xmax, xstep)
yi = np.zeros(len(xi))
for i in range(len(xi)):
for k in range(len(x)):
yi[i] = yi[i] + y[k] * (gamma / 2.0) / (
(xi[i] - x[k]) ** 2 + (gamma / 2.0) ** 2
)
return xi, yi


## Absorption spectra#

### Spectra from ADC eigenstates#

The first six excited states of water is calculated as:

water_xyz = """
O       0.0000000000     0.0000000000     0.1178336003
H      -0.7595754146    -0.0000000000    -0.4713344012
H       0.7595754146     0.0000000000    -0.4713344012
"""

# Construct structure and basis objects
struct = gator.get_molecule(water_xyz)
basis = gator.get_molecular_basis(struct, "6-31G")

# Perform SCF calculation
scf_gs = gator.run_scf(struct, basis)

# Calculate the 6 lowest eigenstates


The resuls can be printed as a table, and convoluted and plotted using built-in functionalities (which can use several different energy-axis, as shown below).

# Print information on eigenstates

# Plot using built-in functionalities
plt.figure(figsize=(10, 6))
plt.subplot(221)
plt.title("In eV")
plt.subplot(222)
plt.title("In atomic units")
plt.subplot(223)
plt.title("In nm")
plt.subplot(224)
plt.title(r"In cm$^{-1}$")
plt.tight_layout()
plt.show()

+--------------------------------------------------------------+
| adc3 (adc2)                             singlet ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0     0.3127134      8.509364   0.0134     0.951   0.04895  |
|  1     0.3933877      10.70463   0.0000    0.9532   0.04681  |
|  2     0.4056694      11.03883   0.1161    0.9491   0.05091  |
|  3     0.4916579      13.37869   0.1103     0.949   0.05096  |
|  4     0.5655096       15.3883   0.4351    0.9625   0.03751  |
|  5     0.6977152       18.9858   0.2604    0.9537   0.04626  |
+--------------------------------------------------------------+ Note that the nm selection (which is the inverse of the energy) is here not broadened.

The built-in functionality uses a small default broadening and plots the full region which is resolved, while we might want to focus on lower-energy and larger broadenings, as can be done with manual routines:

plt.figure(figsize=(6, 6))
# Convolute using functionalities available in gator and adcc
plt.subplot(211)

# Convoluted using information on excitation energies and oscillator strengths
plt.subplot(212)
plt.title("Custom broadening (0.5 eV)")
x, y = au2ev * adc_res.excitation_energy, adc_res.oscillator_strength
xi, yi = lorentzian(x, y, 7, 14, 0.01, 0.5)
plt.plot(xi, yi)
plt.tight_layout()
plt.show() We further note that high-energy features can often be a result of a discretized continuum region. A larger basis set will flatten out this region, but care should be taken for any analysis of that part of the spectrum.

### Spectra from TDDFT eigenstates#

Calculating the first six eigenstates using TDDFT:

# Prepare molecule and basis objects
basis = vlx.MolecularBasis.read(molecule, "6-31G")

# SCF settings and calculation
scf_drv = vlx.ScfRestrictedDriver(comm, ostream=silent_ostream)
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(molecule, basis)

# resolve four eigenstates
rpa_solver = LinearResponseEigenSolver(comm, ostream=silent_ostream)
rpa_solver.update_settings({"nstates": 6}, method_settings)
rpa_results = rpa_solver.compute(molecule, basis, scf_drv.scf_tensors)


There are currently no built-in functionalities for obtaining the broadened spectra, so we instead construct this from the eigenvalues and oscillator strengths, as well as printing a table with energies (here in atomic units), oscillator strengths, and transition dipole moments:

# Print results as a table
print("Energy [au]  Osc. str.   TM(x)     TM(y)     TM(z)")
for i in np.arange(len(rpa_results["eigenvalues"])):
e, os, x, y, z = (
rpa_results["eigenvalues"][i],
rpa_results["oscillator_strengths"][i],
rpa_results["electric_transition_dipoles"][i],
rpa_results["electric_transition_dipoles"][i],
rpa_results["electric_transition_dipoles"][i],
)
print("   {:.3f}     {:8.5f}  {:8.5f}  {:8.5f}  {:8.5f}".format(e, os, x, y, z))

plt.figure(figsize=(6, 4))
x = au2ev * rpa_results["eigenvalues"]
y = rpa_results["oscillator_strengths"]
xi, yi = lorentzian(x, y, min(x) - 1.0, max(x) + 1.0, 0.01, 0.5)
plt.plot(xi, yi)
plt.tight_layout()
plt.show()

Energy [au]  Osc. str.   TM(x)     TM(y)     TM(z)
0.286      0.01128   0.00000   0.24321   0.00000
0.364      0.09652  -0.00000  -0.00000   0.63080
0.364      0.00000  -0.00000   0.00000  -0.00000
0.454      0.08684  -0.53570  -0.00000  -0.00000
0.540      0.41655  -1.07524   0.00000   0.00000
0.666      0.24377  -0.00000   0.00000  -0.74073 The resulting spectra agrees well with ADC(3), with a common shift in energies of about 0.7 eV:

plt.figure(figsize=(6, 4))
xmin,xmax = 7,17
x = au2ev * rpa_results["eigenvalues"]
y = rpa_results["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin,xmax, 0.01, 0.5)
plt.plot(xi, yi)
x, y = au2ev * adc_res.excitation_energy, adc_res.oscillator_strength
xi, yi = lorentzian(x, y, xmin,xmax, 0.01, 0.5)
plt.plot(xi, yi)
plt.xlim((xmin,xmax))
plt.tight_layout()
plt.show() ### Spectra from CPP-DFT#

Using CPP-DFT, the linear absorption cross-section is resolved over a range of energies, which is chosen as the 7-17 eV from above (with a resolution of 0.1 eV):

# Define spectrum region to be resolved
freqs = np.arange(7.0, 17.0, 0.1) / au2ev
freqs_str = [str(x) for x in freqs]

# Calculate the response
cpp_prop = LinearAbsorptionCrossSection(
{"frequencies": ",".join(freqs_str), "damping": 0.3 / au2ev}, method_settings
)
cpp_prop.init_driver(comm, ostream=silent_ostream)
cpp_prop.compute(molecule, basis, scf_drv.scf_tensors)

# Extract the imaginary part of the complex response function and convert to absorption cross section
sigma = []
for w in freqs:
axx = -cpp_prop.rsp_property["response_functions"][("x", "x", w)].imag
ayy = -cpp_prop.rsp_property["response_functions"][("y", "y", w)].imag
azz = -cpp_prop.rsp_property["response_functions"][("z", "z", w)].imag
alpha_bar = (axx + ayy + azz) / 3.0
sigma.append(4.0 * np.pi * w * alpha_bar / 137.035999)


Plotting the raw output, the raw output with a splined function, and a comparison to the eigenstate results above (here plotted as bars):

# Make figure with panels of 3:1 width
plt.figure(figsize=(9, 8))
gs = gridspec.GridSpec(3, 2, width_ratios=[3, 1])

# Raw results for the full region
plt.subplot(gs)
plt.plot(au2ev * freqs, sigma, "bx-")
plt.legend(("Raw", ""),loc='upper left')
plt.xlim((xmin,xmax))

# Raw results for a zoomed in region
plt.subplot(gs)
plt.plot(au2ev * freqs, sigma, "bx-")
plt.xlim((9.4, 10.4))
plt.ylim((0, 0.50))

# Raw and splined spectra for the full region
plt.subplot(gs)
plt.plot(au2ev * freqs, sigma, "bx")
x = np.arange(min(au2ev * freqs), max(au2ev * freqs), 0.01)
y = interp1d(au2ev * freqs, sigma, kind="cubic")
plt.plot(x, y(x), "r")
plt.legend(("Raw", "Splined"),loc='upper left')
plt.xlim((xmin,xmax))

# Zoomed in raw and splined spectra
plt.subplot(gs)
plt.plot(au2ev * freqs, sigma, "bx")
plt.plot(x, y(x), "r")
plt.xlim((9.4, 10.4))
plt.ylim((0, 0.50))

# Zoomed in raw and splined spectra for the full region
plt.subplot(gs)
x = np.arange(min(au2ev * freqs), max(au2ev * freqs), 0.01)
y = interp1d(au2ev * freqs, sigma, kind="cubic")
plt.plot(x, y(x), "r")
xi = au2ev * rpa_results["eigenvalues"]
yi = rpa_results["oscillator_strengths"]
plt.bar(xi, 3.0*yi,width=0.1, color='k')
plt.legend(("Splined", "Eigenstates"),loc='upper left')
plt.xlim((xmin,xmax))

# Zoomed in raw and splined spectra
plt.subplot(gs)
plt.plot(x, y(x), "r")
plt.bar(xi, 3.0*yi,width=0.25, color='k')
plt.xlim((9.4, 10.4))
plt.ylim((0, 0.50))

plt.show() ### Spectra from MCSCF#

In some cases, it can be preferable to use a multiconfigurational wavefunction to compute excitation energies. This is necessary if, for instance, the molecule is suspected to have strong correlation effects, or if the user is interested in analyzing a conical intersection — something that ADC and TD-DFT often fail to properly describe.

In principle, a MCSCF spectrum can be computed using response theory similarly to ADC and DFT. However, as configuration interaction (CI) can naturally provide not just the lowest but any state, it is possible to produce excited states by simply increasing the number of requested states or “roots” in the CI. However, in traditional MCSCF, the orbitals cannot be simultaneously optimized for each state, and instead we use a set of orbitals that is a compromise between all states: the state-averaged orbitals.

For water in its equilibrium distance, there is no strong correlation, and thus the only orbitals we need to include in the active space are those that can be excited in the UV-visible spectrum. Here, will use a CASSCF with an active space comprising the molecular orbitals formed by the oxygen 2p and hydrogen 1s, corresponding to the 2 $$\sigma$$ and $$\sigma^*$$ and one oxygen lone pair. We could in principle include the other lone pair, formed mostly by the 2s orbital of the oxygen, but its orbital energy is much lower and the orbital is thus not involved in the lowest UV-visible transitions.

The orbitals are conveniently located around the HOMO-LUMO gap, so it is sufficient to request a CAS(6,5) (6 electrons in 5 orbitals) to get the desired active space. First, we calculate the SCF ground state:

# Prepare molecule and basis objects
basis = vlx.MolecularBasis.read(molecule, "6-31G")

# SCF calculation
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.compute(molecule, basis)


Then we resolve the six lowest excited states:

# Active space settings
space = mtp.OrbSpace(molecule, scf_drv.mol_orbs)
space.CAS(6, 5)  # 3 O_2p and 2 H_1s

# CASSCF calculation
mcscf_drv = mtp.McscfDriver()
mcscf_drv.compute(molecule, basis, space, 6)  # 6 excited states

# Transition properties
SI = mtp.InterState()
DipOsc = SI.diposc(molecule, basis, mcscf_drv.CIVecs)

          Active space definition:
------------------------
Number of inactive (occupied) orbitals: 2
Number of active orbitals:              5
Number of virtual orbitals:             6

This is a CASSCF wavefunction: CAS(6,5)

CI expansion:
-------------
Number of determinants:      55

MCSCF Iterations
----------------

Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. |   Time
---------------------------------------------------------------------

        1     -75.617648941     0.0e+00      4.8e-01          1   0:00:00

        2     -75.642903350    -2.5e-02      1.7e-01          1   0:00:00

        3     -75.646125925    -3.2e-03      1.9e-02          1   0:00:00

        4     -75.646212175    -8.6e-05      3.1e-03          1   0:00:00

        5     -75.646218727    -6.6e-06      1.2e-03          1   0:00:00

        6     -75.646219045    -3.2e-07      2.1e-04          1   0:00:00

        7     -75.646219063    -1.8e-08      5.1e-05          1   0:00:00

        8     -75.646219065    -1.9e-09      1.4e-05          1   0:00:00


** Convergence reached in 8 iterations


        Final results
-------------

* State 1
- Energy: -75.98281079982928
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99960 1.99159 1.98486 0.01544 0.00851

* State 2
- Energy: -75.7040497309247
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.00000 1.99728 1.98828 0.01531 0.99912

* State 3
- Energy: -75.62045164342095
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.00000 1.99422 1.99530 0.99784 0.01264

* State 4
- Energy: -75.6022500250873
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99969 0.77478 1.96141 0.04148 1.22264

* State 5
- Energy: -75.5200407815364
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99960 1.01676 1.97752 0.97865 0.02746

* State 6
- Energy: -75.44771140909488
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99895 1.97785 1.01935 0.02395 0.97990


Spin Restricted Orbitals
------------------------

Molecular Orbital No.   1:
--------------------------
Occupation: 2.000 Energy:  -20.68879 a.u.
(   1 O   1s  :    -1.00)

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.45143 a.u.
(   1 O   1s  :     0.22) (   1 O   2s  :    -0.49) (   1 O   3s  :    -0.48)

Molecular Orbital No.   3:
--------------------------
Occupation: 1.816 Energy:   -0.76930 a.u.
(   1 O   1p+1:    -0.56) (   1 O   2p+1:    -0.29) (   2 H   1s  :     0.23)
(   3 H   1s  :    -0.23)

Molecular Orbital No.   4:
--------------------------
Occupation: 1.666 Energy:   -0.52443 a.u.
(   1 O   1p-1:    -0.67) (   1 O   2p-1:    -0.48)

Molecular Orbital No.   5:
--------------------------
Occupation: 1.666 Energy:   -0.59973 a.u.
(   1 O   2s  :    -0.18) (   1 O   3s  :    -0.24) (   1 O   1p0 :    -0.60)
(   1 O   2p0 :    -0.40)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.501 Energy:    0.10225 a.u.
(   1 O   2s  :     0.16) (   1 O   3s  :     1.09) (   1 O   1p0 :    -0.25)
(   1 O   2p0 :    -0.40) (   2 H   2s  :    -0.94) (   3 H   2s  :    -0.94)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.351 Energy:    0.21690 a.u.
(   1 O   1p+1:     0.42) (   1 O   2p+1:     0.69) (   2 H   2s  :     1.24)
(   3 H   2s  :    -1.24)

Molecular Orbital No.   8:
--------------------------
Occupation: 0.000 Energy:    0.99875 a.u.
(   1 O   2p+1:    -0.69) (   2 H   1s  :    -0.97) (   2 H   2s  :     0.59)
(   3 H   1s  :     0.97) (   3 H   2s  :    -0.59)

Molecular Orbital No.   9:
--------------------------
Occupation: 0.000 Energy:    1.06993 a.u.
(   1 O   1p-1:     0.94) (   1 O   2p-1:    -1.05)

Molecular Orbital No.  10:
--------------------------
Occupation: 0.000 Energy:    1.08788 a.u.
(   1 O   1p0 :     0.94) (   1 O   2p0 :    -0.96) (   2 H   1s  :     0.30)
(   2 H   2s  :    -0.41) (   3 H   1s  :     0.30) (   3 H   2s  :    -0.41)



Total MCSCF time: 00:00:00


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2        7.58548         1.220154e-02    2.799323e-02
1       4       10.35559         1.097824e-01    1.281782e-01
1       5       12.59261         1.545527e-01    1.469866e-01
1       6       14.56080         5.993590e-01    3.307423e-01


The resulting states are printed above, or can be printed or plotted:

# Print results as a table
print("Energy [au]  Osc. str.")
for i in np.arange(len(DipOsc["energies"])):
e, os = (
DipOsc["energies"][i],
DipOsc["oscillator_strengths"][i],
)
print("   {:.3f}     {:8.5f}".format(e, os))

plt.figure(figsize=(6, 4))
x = au2ev * DipOsc["energies"]
y = DipOsc["oscillator_strengths"]
xi, yi = lorentzian(x, y, min(x) - 1.0, max(x) + 1.0, 0.01, 0.5)
plt.plot(xi, yi)
plt.tight_layout()
plt.show()

Energy [au]  Osc. str.
0.279      0.01220
0.362      0.00000
0.381      0.10978
0.463      0.15455
0.535      0.59936 Comparing to ADC(3) we see a good agreement in relative features, but the absolute energies and intensities are a bit different:

plt.figure(figsize=(6, 4))
xmin,xmax = 7,17
x, y = au2ev * adc_res.excitation_energy, adc_res.oscillator_strength
xi, yi = lorentzian(x, y, xmin,xmax, 0.01, 0.5)
plt.plot(xi, yi)
# SA-MCSCF
x = au2ev * DipOsc["energies"]
y = DipOsc["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin,xmax, 0.01, 0.5)
plt.plot(xi, yi) 