Absorption#

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

import gator
import matplotlib.pyplot as plt
import multipsi as mtp
import numpy as np
import py3Dmol as p3d
import veloxchem as vlx
from matplotlib import gridspec
from scipy.interpolate import interp1d

# au to eV conversion factor
au2ev = 27.211386

We here consider the UV/vis spectrum of gaseous water, with molecular structure:

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

viewer = p3d.view(width=400, height=300)
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})
viewer.addModel("3\n" + water_mol_str)
viewer.setStyle({"stick": {}})
viewer.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

TDDFT#

Calculating the first six eigenstates using TDDFT:

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

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

# resolve four eigenstates
rpa_solver = vlx.lreigensolver.LinearResponseEigenSolver()
rpa_solver.update_settings({"nstates": 6}, method_settings)
rpa_results = rpa_solver.compute(molecule, basis, scf_results)

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:

Note

Broadening and \(\sigma\) extraction routines are currently being implemented.

def lorentzian(x, y, xmin, xmax, xstep, 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


# 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][0],
        rpa_results["electric_transition_dipoles"][i][1],
        rpa_results["electric_transition_dipoles"][i][2],
    )
    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.00000  -0.24321
   0.364      0.09652  -0.00000   0.63080  -0.00000
   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.74073  -0.00000
../../_images/c79859dcdded8a6ddbafdf6b6cc15f78162031c8e4185211726beaf0b0e8c7ff.png

CPP-DFT#

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

Note

The frequency specification is currently being rewritten.

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

# Calculate the response
cpp_drv = vlx.rsplinabscross.LinearAbsorptionCrossSection(
    {"frequencies": ",".join(freqs_str), "damping": 0.3 / au2ev}, method_settings
)
cpp_drv.init_driver()
cpp_results = cpp_drv.compute(molecule, basis, scf_results)

# Extract the imaginary part of the complex response function and convert to absorption cross section
sigma = []
for w in freqs:
    axx = -cpp_drv.rsp_property["response_functions"][("x", "x", w)].imag
    ayy = -cpp_drv.rsp_property["response_functions"][("y", "y", w)].imag
    azz = -cpp_drv.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[0])
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[1])
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[2])
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[3])
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[4])
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", "States"), loc="upper left")
plt.xlim((xmin, xmax))

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

plt.show()
../../_images/dbee3386027a83ac8c4d20e89f63a328756335f4f0092a04f2c3a60429b46b48.png

ADC#

The first six excited states of water is calculated as:

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

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

# Calculate the 6 lowest eigenstates
adc_results = gator.run_adc(molecule, basis, scf_gs, method="adc3", singlets=6)

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).

Note

The adc3 (adc2) designation from adc_results.describe() means that energies are calculated at an ADC(3) level, while properties are given at an ADC(2) level. This is sometimes referred to as ADC(3/2), as well.

# Print information on eigenstates
print(adc_results.describe())

# Plot using built-in functionalities
plt.figure(figsize=(10, 6))
plt.subplot(221)
plt.title("In eV")
adc_results.plot_spectrum(xaxis="eV")

plt.subplot(222)
plt.title("In atomic units")
adc_results.plot_spectrum(xaxis="au")

plt.subplot(223)
plt.title("In nm")
adc_results.plot_spectrum(xaxis="nm", broadening=None)

plt.subplot(224)
plt.title(r"In cm$^{-1}$")
adc_results.plot_spectrum(xaxis="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  |
+--------------------------------------------------------------+
../../_images/3beb430beea29004a5da6a5dabfa740e6e711c93c59d4aa709940911ce39802e.png

Fix: broaden spectra expressed in wavelength.

Note that some high-energy features can be the 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.

CPP-ADC#

To be added.

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 TDDFT 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, we will use a CASSCF with an active space comprising the molecular orbitals formed by the oxygen 2p and hydrogen 1s, corresponding to the two \(\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
molecule = vlx.Molecule.read_molecule_string(water_mol_str)
basis = vlx.MolecularBasis.read(molecule, "6-31G")

# SCF calculation
scf_drv = vlx.ScfRestrictedDriver()
scf_results = 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_results = mcscf_drv.compute(molecule, basis, space, 6)

# Transition properties
SI = mtp.InterState()
DipOsc = SI.compute(molecule, basis, mcscf_results)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Self-Consistent Field Driver
                         ====================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          State-averaged calculation
            Number of states      : 6
            Equal-weights 
          Max. iterations         : 50
          BFGS window             : 5
          Convergence thresholds:
            - Energy change       : 1e-08
            - Gradient norm       : 0.0001
                                                                                                                          
          Integrals in memory
                                                                                                                          
          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
        9     -75.646219065    -5.3e-11      2.7e-06          1   0:00:00
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -75.98281062496628
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98486 1.99960 1.99160 0.00851 0.01544
                                                                                                                          
* State 2
- Energy: -75.70404926196595
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98828 1.00000 1.99728 0.99912 0.01531
                                                                                                                          
* State 3
- Energy: -75.6204513864733
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99530 1.00000 1.99422 0.01264 0.99784
                                                                                                                          
* State 4
- Energy: -75.6022504197667
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.96141 1.99969 1.22263 0.77478 0.04148
                                                                                                                          
* State 5
- Energy: -75.52004110975075
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.97752 1.99960 1.01676 0.02746 0.97865
                                                                                                                          
* State 6
- Energy: -75.44771158728817
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.01935 1.99895 1.97785 0.97990 0.02395
                                                                                                                          
                                                 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   1p0 :    -0.67) (   1 O   2p0 :    -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   1p-1:     0.60)                              
               (   1 O   2p-1:     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   1p-1:     0.25)                              
               (   1 O   2p-1:     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   1p0 :    -0.94) (   1 O   2p0 :     1.05)                                                        
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :         0.000000 a.u.         0.000000 Debye                                     
                                   Y   :        -0.848594 a.u.        -2.156911 Debye                                     
                                   Z   :         0.000000 a.u.         0.000000 Debye                                     
                                 Total :         0.848594 a.u.         2.156911 Debye                                     
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 2                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :         0.088075 a.u.         0.223865 Debye                                     
                                   Z   :        -0.000000 a.u.        -0.000000 Debye                                     
                                 Total :         0.088075 a.u.         0.223865 Debye                                     
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 3                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :         0.000000 a.u.         0.000000 Debye                                     
                                   Y   :        -0.019117 a.u.        -0.048591 Debye                                     
                                   Z   :         0.000000 a.u.         0.000000 Debye                                     
                                 Total :         0.019117 a.u.         0.048591 Debye                                     
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 4                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :         0.127030 a.u.         0.322878 Debye                                     
                                   Z   :        -0.000000 a.u.        -0.000000 Debye                                     
                                 Total :         0.127030 a.u.         0.322878 Debye                                     
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 5                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :         0.000000 a.u.         0.000000 Debye                                     
                                   Y   :         0.138716 a.u.         0.352582 Debye                                     
                                   Z   :        -0.000000 a.u.        -0.000000 Debye                                     
                                 Total :         0.138716 a.u.         0.352582 Debye                                     
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 6                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :        -0.255728 a.u.        -0.649995 Debye                                     
                                   Z   :        -0.000000 a.u.        -0.000000 Debye                                     
                                 Total :         0.255728 a.u.         0.649995 Debye                                     
                                                                                                                          
                                                                                                                          
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.220177e-02    2.799370e-02
     1       4       10.35557         1.097829e-01    1.281789e-01
     1       5       12.59260         1.545515e-01    1.469851e-01
     1       6       14.56079         5.993553e-01    3.307390e-01
                                                                                                                          
List of rotatory strengths greater than 1e-10
                                                                                                                          
  From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)

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
../../_images/bbf9e19e204b51b99318ed1cdd8bae2adbc83eb1fab80bda348aed68c1e85b7e.png

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))

# ADC(3)
x, y = au2ev * adc_results.excitation_energy, adc_results.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)
plt.legend(("ADC(3)", "SA-MCSCF"))
plt.xlim((xmin, xmax))
plt.tight_layout()
plt.show()
../../_images/2f049ed84647d9727eb32f4d0fa678bdfcb397c6cc9b21d9bfd3cbc9817682f1.png

As an alternative to using state-averaging, we can also use linear response CASSCF. One advantage in this case is that we do not need to add the oxygen lone pair in the active space, since the response includes also excitation outside of the active space. Let’s try it:

# Active space settings
space = mtp.OrbSpace(molecule, scf_drv.mol_orbs)

# The list of active orbitals, orbital 4 is not here since it is the lone pair
space.cas_orbitals([2, 3, 5, 6])

# CASSCF calculation
mcscf_drv = mtp.McscfDriver()
mcscf_results = mcscf_drv.compute(molecule, basis, space)  # Only ground state

# Response calculation
mcrpa = mtp.Mclr_EigenSolver()
rspdict = mcrpa.compute(molecule, basis, mcscf_drv, 5)
Hide code cell output
                                                                                                                          
                          Multi-Configurational Self-Consistent Field Driver
                         ====================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          State-specific calculation
          Max. iterations         : 50
          BFGS window             : 5
          Convergence thresholds:
            - Energy change       : 1e-08
            - Gradient norm       : 0.0001
                                                                                                                          
          Integrals in memory
                                                                                                                          
          Active space definition:
          ------------------------
Number of inactive (occupied) orbitals: 3
Number of active orbitals:              4
Number of virtual orbitals:             6

    This is a CASSCF wavefunction: CAS(4,4)

          CI expansion:
          -------------
Number of determinants:      21


                                                                                                                          
                   MCSCF Iterations
     ---------------------------------------------------------------------
                                                                                                                          
     Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. |   Time
     ---------------------------------------------------------------------
        1     -75.991165249     0.0e+00      3.9e-02          1   0:00:00
        2     -76.005855741    -1.5e-02      6.6e-02          1   0:00:00
        3     -76.031621117    -2.6e-02      3.4e-02          1   0:00:00
        4     -76.034371736    -2.8e-03      4.0e-02          1   0:00:00
        5     -76.036920247    -2.5e-03      1.6e-02          1   0:00:00
        6     -76.037128966    -2.1e-04      7.3e-03          1   0:00:00
        7     -76.037231330    -1.0e-04      5.0e-03          1   0:00:00
        8     -76.037270318    -3.9e-05      3.3e-03          1   0:00:00
        9     -76.037273265    -2.9e-06      1.5e-03          1   0:00:00
       10     -76.037273978    -7.1e-07      4.1e-04          1   0:00:00
       11     -76.037274020    -4.2e-08      2.2e-04          1   0:00:00
       12     -76.037274031    -1.1e-08      6.1e-05          1   0:00:00
       13     -76.037274032    -1.3e-09      1.7e-05          1   0:00:00
** Convergence reached in 13 iterations
       14     -76.037274032    -2.0e-10      5.0e-06          1   0:00:00
                                                                                                                          
        Final results
        -------------
                                                                                                                          
                                                                                                                          
* State 1
- Energy: -76.03727403227734
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.97794 1.97490 0.02386 0.02330
                                                                                                                          
                                                 Spin Restricted Orbitals                                                 
                                                 ------------------------                                                 
                                                                                                                          
               Molecular Orbital No.   1:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:  -20.55476 a.u.                                                                  
               (   1 O   1s  :     1.00)                                                                                  
                                                                                                                          
               Molecular Orbital No.   2:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -1.07550 a.u.                                                                  
               (   1 O   1s  :     0.22) (   1 O   2s  :    -0.49) (   1 O   3s  :    -0.58)                              
               (   1 O   1p-1:    -0.24) (   1 O   2p-1:    -0.19)                                                        
                                                                                                                          
               Molecular Orbital No.   3:                                                                                 
               --------------------------                                                                                 
               Occupation: 2.000 Energy:   -0.50212 a.u.                                                                  
               (   1 O   1p0 :    -0.64) (   1 O   2p0 :    -0.51)                                                        
                                                                                                                          
               Molecular Orbital No.   4:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.978 Energy:   -0.83333 a.u.                                                                  
               (   1 O   1p-1:     0.51) (   1 O   2p-1:     0.35) (   2 H   1s  :    -0.21)                              
               (   3 H   1s  :    -0.21)                                                                                  
                                                                                                                          
               Molecular Orbital No.   5:                                                                                 
               --------------------------                                                                                 
               Occupation: 1.975 Energy:   -0.70219 a.u.                                                                  
               (   1 O   1p+1:    -0.51) (   1 O   2p+1:    -0.25) (   2 H   1s  :     0.27)                              
               (   3 H   1s  :    -0.27)                                                                                  
                                                                                                                          
               Molecular Orbital No.   6:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.024 Energy:    0.75293 a.u.                                                                  
               (   1 O   1p+1:     0.87) (   2 H   1s  :     0.46) (   2 H   2s  :     0.30)                              
               (   3 H   1s  :    -0.46) (   3 H   2s  :    -0.30)                                                        
                                                                                                                          
               Molecular Orbital No.   7:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.023 Energy:    0.73231 a.u.                                                                  
               (   1 O   2s  :    -0.49) (   1 O   3s  :    -0.17) (   1 O   1p-1:     0.78)                              
               (   1 O   2p-1:    -0.18) (   2 H   1s  :     0.43) (   2 H   2s  :     0.20)                              
               (   3 H   1s  :     0.43) (   3 H   2s  :     0.20)                                                        
                                                                                                                          
               Molecular Orbital No.   8:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.67761 a.u.                                                                  
               (   1 O   2s  :    -0.18) (   1 O   3s  :     1.21) (   1 O   1p-1:     0.41)                              
               (   1 O   2p-1:    -0.76) (   2 H   1s  :     0.40) (   2 H   2s  :    -1.13)                              
               (   3 H   1s  :     0.40) (   3 H   2s  :    -1.13)                                                        
                                                                                                                          
               Molecular Orbital No.   9:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    0.78417 a.u.                                                                  
               (   1 O   1p+1:    -0.27) (   1 O   2p+1:     0.72) (   2 H   1s  :    -0.62)                              
               (   2 H   2s  :     1.60) (   3 H   1s  :     0.62) (   3 H   2s  :    -1.60)                              
                                                                                                                          
               Molecular Orbital No.  10:                                                                                 
               --------------------------                                                                                 
               Occupation: 0.000 Energy:    1.16377 a.u.                                                                  
               (   1 O   1p0 :     0.96) (   1 O   2p0 :    -1.04)                                                        
                                                                                                                          
                                                                                                                          
                                                Dipole moment for state 1                                                 
                                               ---------------------------                                                
                                                                                                                          
                                   X   :        -0.000000 a.u.        -0.000000 Debye                                     
                                   Y   :        -0.983137 a.u.        -2.498886 Debye                                     
                                   Z   :        -0.000000 a.u.        -0.000000 Debye                                     
                                 Total :         0.983137 a.u.         2.498886 Debye                                     
                                                                                                                          
                                                                                                                          
Total MCSCF time: 00:00:00
                                                                                                                          
                          Multi-Configurational Linear response Eigenvalue Solver
                         =========================================================
                                                                                                                          
        ╭────────────────────────────────────╮
        │          Driver settings           │
        ╰────────────────────────────────────╯
          Number of states        : 5
                                                                                                                          
          Max. iterations         : 50
          Residual norm threshold : 1e-08
                                                                                                                          
          Integrals in memory
          Trial vector settings:
            - Linear dependence threshold :1e-10
            - Min. trial vectors          :20
            - Max. trial vectors          :60
                                                                                                                          
          CI parameters:      20
          Orbital parameters: 54
                                                                                                                          
                                                                                                                          
Initialization time: 00:00:00
                                                                                                                          
                                                                                                                          
        MC-RPA Iterations
        -----------------
                                                                                                                          
     Iter. | Average Energy | Grad. Norm | Converged | CI+Orb vec. | Time
     -----------------------------------------------------------------------
        1       0.711592304      2.9e-01      0/  5     5 +   5      0:00:00    
        2       0.532277685      1.1e-01      0/  5     5 +  15      0:00:00    
        3       0.472274167      4.4e-02      1/  5     7 +  23      0:00:00    
        4       0.455154316      6.4e-03      1/  5     7 +  31      0:00:00    
        5       0.453380481      3.6e-03      2/  5    13 +  32      0:00:00    
        6       0.452031444      7.0e-04      2/  5    13 +  38      0:00:00    
        7       0.451970822      3.6e-04      4/  5    15 +  41      0:00:00    
        8       0.451896428      5.0e-05      4/  5    15 +  43      0:00:00    
        9       0.451894976      2.1e-05      4/  5    15 +  45      0:00:00    
       10       0.451891023      1.0e-06      4/  5    17 +  45      0:00:00    
       11       0.451890969      8.2e-08      4/  5    17 +  47      0:00:00    
       12       0.451890957      6.7e-09      5/  5    19 +  47      0:00:00    
                                                                                                                          
** Convergence reached in 12 iterations
                                                                                                                          
                                                                                                                          
                                                                                                                          
                                                 Oscillator strength (au)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3305189      8.99388          0.0389564        0.0121777
          2     0.4084659     11.11492          0.0000000        0.0000000
          3     0.4307432     11.72112          0.1738586        0.1284586
          4     0.5110794     13.90718          0.1450281        0.1453586
          5     0.5786474     15.74580          0.3465044        0.4062187
                                                                                                                          
                                                                                                                          
                                                 Rotatory strength (cgs)
 Transition   Energy (Ha)    Energy (eV)    Velocity gauge     Length gauge 
 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
          1     0.3305189      8.99388         -0.0000000       -0.0000000
          2     0.4084659     11.11492         -0.0000000       -0.0000000
          3     0.4307432     11.72112          0.0000000        0.0000000
          4     0.5110794     13.90718         -0.0000000       -0.0000000
          5     0.5786474     15.74580          0.0000000        0.0000000
                                                                                                                          
Total MC-RPA time: 00:00:00
plt.figure(figsize=(6, 4))
xmin, xmax = 7, 17

# 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)

# MCSCF response
x = au2ev * rspdict["eigenvalues"]
y = rspdict["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin, xmax, 0.01, 0.5)

plt.plot(xi, yi)
plt.legend(("SA-MCSCF", "MCSCF response"))
plt.xlim((xmin, xmax))
plt.tight_layout()
plt.show()
../../_images/49d9f96c2b59b317128a077ea7b7517d73a5ac7e5366310d154639354ea9b4e6.png

Comparison of spectra#

The spectra from the four different approaches can be plotted as:

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_figheight(4)
fig.set_figwidth(10)

xmin, xmax = 7, 17

# ADC(3)
x, y = au2ev * adc_results.excitation_energy, adc_results.oscillator_strength
xi, yi = lorentzian(x, y, xmin, xmax, 0.01, 0.5)
ax1.plot(xi, yi)
ax2.plot(xi, yi)

# TDDFT
x = au2ev * rpa_results["eigenvalues"]
y = rpa_results["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin, xmax, 0.01, 0.5)
ax1.plot(xi, yi)
ax2.plot(xi, yi + 1)

# SA-MCSCF
x = au2ev * DipOsc["energies"]
y = DipOsc["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin, xmax, 0.01, 0.5)
ax1.plot(xi, yi)
ax2.plot(xi, yi + 2)

# MCSCF response
x = au2ev * rspdict["eigenvalues"]
y = rspdict["oscillator_strengths"]
xi, yi = lorentzian(x, y, xmin, xmax, 0.01, 0.5)
ax1.plot(xi, yi)
ax2.plot(xi, yi + 3)

ax1.legend(("ADC(3)", "TDDFT", "SA-MCSCF", "MCSCF response"))
ax1.set_xlim((xmin, xmax))
ax2.set_xlim((xmin, xmax))
plt.tight_layout()
plt.show()
../../_images/77976783ffea9f5d23c3217e0df079dea3722f145c6c8abfe8c1e794428e071e.png

Plotted either on top of each other, or with a vertical set-off. We note that the features are, in general, similar, with differences mainly taking the form of absolute energy shifts, as well as noticeably more intense features for SA-MCSCF.