# CI response#

## Full-CI response#

As we saw in the previous section, for an exact state, we have the formula

$\langle \! \langle \hat{\Omega}; \hat{V}^{\omega} \rangle \! \rangle = - \frac{1}{\hbar} \sum_{n>0} \left[ \frac{\langle 0 | \hat{\Omega} | n \rangle \langle n | \hat{V}^{\omega} | 0 \rangle } {\omega_{n0}-\omega} + \frac{\langle 0 | \hat{V}^{\omega} | n \rangle \langle n | \hat{\Omega} | 0 \rangle } {\omega_{n0}+\omega} \right]$

For example, in the case of a molecule irradiated by light, we are interested in computing the polarizability $$\alpha_{\alpha\beta}(-\omega; \omega)$$, corresponding to the response with $$\hat{\Omega} = \hat{\mu}_\alpha$$ and $$\hat{V}^\omega = - \hat{\mu}_\beta$$, where $$\alpha, \beta \in \{x,y,z\}$$. Looking more precisely at the isotropic component, $$\bar\alpha (-\omega; \omega) = \frac{1}{3} \sum_\alpha \alpha_{\alpha\alpha} (-\omega; \omega)$$, we obtain

$\bar\alpha (-\omega; \omega) = \frac{2}{3} \sum_{n>0} \sum_\alpha \frac{ \omega_{n0} | \langle 0 | \hat{\mu}_\alpha | n \rangle |^2 }{ \omega_{n0}^2 - \omega^2 }$

We now define the oscillator strength $$f_n$$ for a transition $$0 \rightarrow n$$ as

$f_n = \frac{2}{3} \sum_\alpha \omega_{n0} | \langle 0 | \hat{\mu}_\alpha | n \rangle |^2$

such that

$\bar\alpha (-\omega; \omega) = \sum_{n>0} \frac{ f_n }{ \omega_{n0}^2 - \omega^2 }$

Those formulas can be applied directly for a configuration interaction wave function, where it is easy to get the excited state energies simply by obtaining additional roots in the diagonalization of the Hamiltonian. We then only need to compute the matrix element $$\langle 0 | \hat{\mu}_\alpha | n \rangle$$. In second quantization, the $$\alpha$$ component of the dipole operator can be written as

$\hat{\mu}_{\alpha} = \sum_{pq} \mu_{\alpha, pq} \hat{a}_p^\dagger \hat{a}_q$

We thus need to compute the matrix $$\boldsymbol{\gamma}$$, with elements $$\gamma_{pq} =\langle 0 |\hat{a}_p^\dagger \hat{a}_q | n \rangle$$, which is called the transition density matrix, and contract it with the dipole moment integrals.

Let’s apply these for full CI (FCI) of a water molecule.

import veloxchem as vlx
import multipsi as mtp
import numpy as np

mol_str = """3

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

scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

# Compute the 5 lowest states of water using full CI

space=mtp.OrbSpace(molecule,scf_drv.mol_orbs)
# Full CI with frozen 1s orbital
space.fci(n_frozen=1)

nstates=5
cidrv=mtp.CIDriver()
ci_results = cidrv.compute(molecule,basis,space,nstates)
CIVecs = cidrv.vecs

Hide code cell output
                          Configuration Interaction Driver
==================================

Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

This is a CASSCF wavefunction: CAS(8,12)

CI expansion:
-------------
Number of determinants:      122760

╭────────────────────────────────────╮
│          Driver settings           │
╰────────────────────────────────────╯
Max. iterations         : 40
Initial diagonalization : 200
Max subspace size       : 50
Convergence thresholds:
- Energy change       : 1e-08
- Residual square norm: 1e-08
Standard Davidson step



CI Iterations
-------------

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

        1     -75.708185001     0.0e+00      3.6e-01   0:00:00

        2     -75.793127822    -9.5e-02      3.7e-02   0:00:00

        3     -75.799861180    -7.3e-03      3.6e-03   0:00:00

        4     -75.800508543    -7.7e-04      3.9e-04   0:00:00

        5     -75.800565380    -8.7e-05      7.3e-05   0:00:00

        6     -75.800574227    -1.6e-05      8.3e-06   0:00:00

        7     -75.800575131    -1.7e-06      1.2e-06   0:00:00

        8     -75.800575290    -3.6e-07      4.2e-07   0:00:00

        9     -75.800575329    -1.4e-07      1.7e-07   0:00:00

       10     -75.800575338    -4.2e-08      2.9e-08   0:00:00

       11     -75.800575340    -5.8e-09      3.9e-09   0:00:00


** Convergence reached in 11 iterations

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

* State 1
- Energy: -76.12020540305299
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98825 1.96801 1.97145 1.98065 0.00050 0.02821 0.00310 0.01813 0.02666 0.01218 0.00063 0.00222

* State 2
- Energy: -75.80967663263324
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98999 1.96500 1.97826 0.99416 0.99612 0.03529 0.00360 0.00679 0.01598 0.00969 0.00360 0.00263

* State 3
- Energy: -75.72755706157507
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98934 1.97697 1.96858 0.99423 0.03496 0.99438 0.00381 0.00671 0.00220 0.01152 0.01394 0.00335

* State 4
- Energy: -75.71682482474854
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98622 1.95656 1.14202 1.98084 0.84833 0.04126 0.00283 0.01796 0.01163 0.00636 0.00360 0.00240

* State 5
- Energy: -75.62861277592349
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98447 1.95661 1.01098 1.98207 0.04117 0.97797 0.00370 0.01636 0.00441 0.00643 0.01380 0.00441

au2ev = 27.211386

# Get the molecular orbitals within the active space
C=space._get_active_mos()

# Compute dipole integrals in AO basis
dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)
dipole=[dipole_mats.x_to_numpy(),dipole_mats.y_to_numpy(),dipole_mats.z_to_numpy()]

# Transform to MO basis
dipole_mo=[]
for icomp in range(0,3):
dipole_mo.append(np.einsum('pq,pt,qu->tu', dipole[icomp], C, C))

# Initialize the CIOperator class to compute the transition densities.
expansion=mtp.CIExpansion(space)
DenDriver=mtp.CIOperator(expansion)

# Compute all 0->n transitions
for n in range(1,nstates):
dE=CIVecs[n].energy-CIVecs[0].energy
tden=DenDriver.get1den(CIVecs[0],CIVecs[n]) #Transition density matrix
dx=np.tensordot(tden,dipole_mo[0]) #<0|x|n>
dy=np.tensordot(tden,dipole_mo[1]) #<0|y|n>
dz=np.tensordot(tden,dipole_mo[2]) #<0|z|n>
F=2/3*dE*(dx*dx+dy*dy+dz*dz)
print(f"Excitation 0->{n} energy: {dE*au2ev:5.3f} eV oscillator strength: {F:.5f}")

Excitation 0->1 energy: 8.450 eV oscillator strength: 0.01319
Excitation 0->2 energy: 10.685 eV oscillator strength: 0.00000
Excitation 0->3 energy: 10.977 eV oscillator strength: 0.11535
Excitation 0->4 energy: 13.377 eV oscillator strength: 0.11499


We thus obtain 3 visible transitions and one dark state (near-zero dipole oscillator strength). We obtain the exact same result by using the built-in class of multipsi:

SI=mtp.InterState()
DipOsc=SI.compute(molecule, basis, ci_results)


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2        8.44992         1.318937e-02    4.625297e-02
1       4       10.97654         1.153538e-01    1.749710e-01
1       5       13.37692         1.149882e-01    1.223962e-01


List of rotatory strengths greater than 1e-10

From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)


## Truncated CI response#

While they do not formally qualify as “exact state”, the response equations for truncated CI are the same as the exact case above. One can thus, in principle, obtain a hierarchy of method from CIS to full CI. However, the excitation energies do not necessarily improve from order to order. Let us demonstrate this on the water example.

nstates=5
# nstates-1 transitions for CIS, CISD, CISDT and CISDTQ and FCI
Energies=np.empty((5,nstates-1))
Energies[4,:]=au2ev*DipOsc['energies'] # Save the FCI result

#CIS to CISDTQ
space=mtp.OrbSpace(molecule,scf_drv.mol_orbs)
cidrv=mtp.CIDriver()
SI=mtp.InterState()

for exc in range(1,5):
space.ci(exc,n_frozen=1)
ci_results = cidrv.compute(molecule,basis,space, nstates)
DipOsc=SI.compute(molecule,basis,ci_results)
Energies[exc-1,:]=au2ev*DipOsc['energies']

Hide code cell output
                          Configuration Interaction Driver
==================================

Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           4               7                8
2          12               8                8

CI expansion:
-------------
Number of determinants:      33

╭────────────────────────────────────╮
│          Driver settings           │
╰────────────────────────────────────╯
Solved by explicit diagonalization



CI Iterations
-------------

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

        1     -75.642674442     0.0e+00      3.1e-27   0:00:00


** Convergence reached in 1 iterations

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

* State 1
- Energy: -75.98387037576913
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
2.00000 2.00000 2.00000 2.00000 -0.00000 -0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

* State 2
- Energy: -75.63914005847884
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
2.00000 2.00000 1.99997 1.00003 0.99997 0.00000 -0.00000 0.00003 0.00000 -0.00000 -0.00000 -0.00000

* State 3
- Energy: -75.56833156634573
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
2.00000 1.99997 2.00000 1.00003 0.00000 0.99997 -0.00000 0.00003 0.00000 -0.00000 -0.00000 0.00000

* State 4
- Energy: -75.54906255700259
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99959 1.98314 1.02031 1.99696 0.97969 0.01686 0.00000 0.00304 0.00000 -0.00000 -0.00000 0.00041

* State 5
- Energy: -75.47296765017845
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99995 1.97129 1.02876 2.00000 0.02871 0.97124 0.00000 0.00000 -0.00000 0.00000 0.00005 0.00000


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2        9.38059         1.468873e-02    4.091187e-02
1       4       11.83172         1.208787e-01    1.355741e-01
1       5       13.90237         1.045445e-01    8.623014e-02


List of rotatory strengths greater than 1e-10

From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)

                          Configuration Interaction Driver
==================================

Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           4               6                8
2          12               8                8

CI expansion:
-------------
Number of determinants:      729

╭────────────────────────────────────╮
│          Driver settings           │
╰────────────────────────────────────╯
Max. iterations         : 40
Initial diagonalization : 200
Max subspace size       : 50
Convergence thresholds:
- Energy change       : 1e-08
- Residual square norm: 1e-08
Olsen step



CI Iterations
-------------

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

        1     -75.708067912     0.0e+00      3.5e-01   0:00:00

        2     -75.739528252    -9.2e-02      9.7e-03   0:00:00

        3     -75.740704223    -3.8e-03      2.2e-04   0:00:00

        4     -75.740732582    -8.3e-05      5.2e-06   0:00:00

        5     -75.740733616    -1.9e-06      1.8e-07   0:00:00

        6     -75.740733655    -1.1e-07      1.2e-08   0:00:00

        7     -75.740733657    -6.5e-09      6.1e-10   0:00:00


** Convergence reached in 7 iterations

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

* State 1
- Energy: -76.11340211569267
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98980 1.97263 1.97571 1.98356 0.00034 0.02409 0.00268 0.01544 0.02279 0.01054 0.00046 0.00195

* State 2
- Energy: -75.7351931372965
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99830 1.99110 1.99600 0.99907 0.99382 0.01175 0.00118 0.00147 0.00444 0.00171 0.00035 0.00082

* State 3
- Energy: -75.65407733183052
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99830 1.99450 1.99144 0.99900 0.01390 0.99166 0.00117 0.00153 0.00080 0.00149 0.00575 0.00047

* State 4
- Energy: -75.64356171507329
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99571 1.97554 1.12840 1.99582 0.86456 0.02646 0.00142 0.00379 0.00507 0.00180 0.00051 0.00091

* State 5
- Energy: -75.55743398573291
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99559 1.96489 1.02434 1.99775 0.03602 0.96813 0.00182 0.00191 0.00134 0.00164 0.00618 0.00039


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2       10.29159         1.548866e-02    3.658376e-02
1       4       12.78501         1.247451e-01    1.355482e-01
1       5       15.12866         1.136810e-01    9.274229e-02


List of rotatory strengths greater than 1e-10

From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)

                          Configuration Interaction Driver
==================================

Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           4               5                8
2          12               8                8

CI expansion:
-------------
Number of determinants:      6329

╭────────────────────────────────────╮
│          Driver settings           │
╰────────────────────────────────────╯
Max. iterations         : 40
Initial diagonalization : 200
Max subspace size       : 50
Convergence thresholds:
- Energy change       : 1e-08
- Residual square norm: 1e-08
Olsen step



CI Iterations
-------------

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

        1     -75.708185001     0.0e+00      3.5e-01   0:00:00

        2     -75.787071882    -9.3e-02      1.4e-02   0:00:00

        3     -75.791912542    -5.3e-03      1.2e-03   0:00:00

        4     -75.792184181    -3.9e-04      8.9e-05   0:00:00

        5     -75.792206026    -4.2e-05      1.1e-05   0:00:00

        6     -75.792207797    -4.0e-06      8.3e-07   0:00:00

        7     -75.792207936    -3.5e-07      7.9e-08   0:00:00

        8     -75.792207948    -2.7e-08      7.6e-09   0:00:00

        9     -75.792207950    -7.4e-09      6.3e-09   0:00:00


** Convergence reached in 9 iterations

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

* State 1
- Energy: -76.11438963659843
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98960 1.97168 1.97493 1.98327 0.00047 0.02478 0.00289 0.01562 0.02338 0.01071 0.00060 0.00208

* State 2
- Energy: -75.8000306066592
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99159 1.97002 1.98187 0.99537 0.99592 0.03029 0.00318 0.00554 0.01330 0.00820 0.00318 0.00237

* State 3
- Energy: -75.71893849552852
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99111 1.97998 1.97319 0.99537 0.03033 0.99427 0.00339 0.00552 0.03033 0.00953 0.01223 0.00311

* State 4
- Energy: -75.70731804012318
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98804 1.96044 1.12926 1.98430 0.86217 0.03726 0.00263 0.01467 0.01016 0.00542 0.00348 0.00217

* State 5
- Energy: -75.620362972389
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.98665 1.95751 1.01386 1.98540 0.04016 0.97604 0.00341 0.01331 0.00418 0.00547 0.01214 0.00418


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2        8.55414         1.328951e-02    4.580239e-02
1       4       11.07698         1.154300e-01    1.705756e-01
1       5       13.44315         1.119484e-01    1.180816e-01


List of rotatory strengths greater than 1e-10

From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)

                          Configuration Interaction Driver
==================================

Active space definition:
------------------------
Number of inactive (occupied) orbitals: 1
Number of active orbitals:              12
Number of virtual orbitals:             0

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           4               4                8
2          12               8                8

CI expansion:
-------------
Number of determinants:      27763

╭────────────────────────────────────╮
│          Driver settings           │
╰────────────────────────────────────╯
Max. iterations         : 40
Initial diagonalization : 200
Max subspace size       : 50
Convergence thresholds:
- Energy change       : 1e-08
- Residual square norm: 1e-08
Olsen step



CI Iterations
-------------

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

        1     -75.708185001     0.0e+00      3.6e-01   0:00:01

        2     -75.793118310    -9.5e-02      3.6e-02   0:00:01

        3     -75.798908867    -7.2e-03      2.0e-03   0:00:01

        4     -75.799401560    -5.8e-04      2.2e-04   0:00:01

        5     -75.799440899    -6.4e-05      2.9e-05   0:00:01

        6     -75.799446125    -9.5e-06      3.3e-06   0:00:01

        7     -75.799446623    -9.0e-07      5.5e-07   0:00:01

        8     -75.799446711    -2.5e-07      2.1e-07   0:00:01

        9     -75.799446734    -9.2e-08      6.1e-08   0:00:01

       10     -75.799446738    -2.0e-08      7.4e-09   0:00:01

       11     -75.799446739    -2.1e-09      8.5e-10   0:00:01


** Convergence reached in 11 iterations

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

* State 1
- Energy: -76.12003281964496
- S^2   : -0.02  (multiplicity = 1.0 )
- Natural orbitals
1.98832 1.96820 1.97163 1.98078 0.00049 0.02805 0.00308 0.01801 0.02650 0.01210 0.00062 0.00221

* State 2
- Energy: -75.80828932732835
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.99041 1.96630 1.97922 0.99448 0.99598 0.03407 0.00351 0.00646 0.01523 0.00935 0.00243 0.00257

* State 3
- Energy: -75.72622990160049
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98982 1.97785 1.96971 0.99453 0.03391 0.99421 0.00372 0.00639 0.00215 0.01098 0.01350 0.00324

* State 4
- Energy: -75.71537749268829
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98672 1.95763 1.14165 1.98172 0.84888 0.04031 0.00280 0.01710 0.01128 0.00608 0.00347 0.00235

* State 5
- Energy: -75.62730415351614
- S^2   : -0.03  (multiplicity = 0.9 )
- Natural orbitals
1.98504 1.95735 1.01140 1.98291 0.04058 0.97774 0.00364 0.01556 0.00427 0.00615 0.01336 0.00427


List of oscillator strengths greather than 1e-10

From     to       Energy (eV)    Oscillator strength (length and velocity)
1       2        8.48297         1.323123e-02    4.607713e-02
1       4       11.01123         1.154562e-01    1.737591e-01
1       5       13.40783         1.150580e-01    1.218071e-01


List of rotatory strengths greater than 1e-10

From     to       Energy (eV)    Rot. strength (a.u. and 10^-40 cgs)

import matplotlib.pyplot as plt

plt.figure(figsize=(6,4))
x = np.array(range(1,6))
plt.plot(x,Energies[:,0], label='0->1')
plt.plot(x,Energies[:,1], label='0->2')
plt.plot(x,Energies[:,2], label='0->3')
plt.plot(x,Energies[:,3], label='0->4')

plt.title('Convergence of excitation energies')
plt.xticks([1, 2, 3, 4, 5], ['CIS', 'CISD', 'CISDT', 'CISDTQ', 'FCI'])
plt.ylabel("Energies (eV)")
plt.legend()
plt.tight_layout(); plt.show()

Energies[0,:]

array([ 9.38058973, 11.30738694, 11.83172339, 13.90237127])


As we can see, we do not have a monotonous convergence towards the full CI result. In particular, CISD is actually worse than CIS. The reason for this is that CIS has a more even-handed treatment of the ground and excited states than CISD. In CIS, due to the Brillouin theorem, the ground state does not improve with the inclusion of single excitations, but the single excitations are needed to generate the dominant configuration for the (singly) excited states. Both ground and excited states are thus treated at a mean-field (HF-like) level.

However, in CISD, the ground state contains double excitation, and thus dynamical correlation, while the excited states would need up to triple excitations to both generate the singly-excited dominant configuration and correlating double excitations on top. There is an unbalance, the ground state being described more accurately than the excited states, leading to too high excitation energies. Adding the triple excitations bring correlation also for excited states, which brings the result closer to the final FCI result, and in this simple case is essentially converged.

For this reason, while CIS is sometimes used to compute excitation energies (and as we can see in the next section, is actually related to HF response), the other truncated CIs are not.