Configuration interaction#
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
"""
molecule = vlx.Molecule.read_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "6-31g")
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
# 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
Show 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.12020540305328
- 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.80967663263307
- 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.00250 0.00263
* State 3
- Energy: -75.72755706157497
- 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.71682482474822
- 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.62861277592332
- 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.00204
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#
In truncated CI, the energy does not just depend on the CI coefficients but also on the molecular orbitals, which makes the truncated CI response more complicated. However, it is interesting to look at the results we obtain by simply applying the above equations. This gives 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']
Show 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 4.7e-27 0:00:00
** Convergence reached in 1 iterations
Final results
-------------
* State 1
- Energy: -75.98387037576919
- 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.6391400584789
- 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.5683315663458
- 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.54906255700254
- 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.47296765017849
- 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.1134021156927
- 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.73519313729659
- 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.64356171507328
- 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.55743398573307
- 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.11438963659847
- 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.80003060665926
- 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.00236 0.00237
* State 3
- Energy: -75.71893849552859
- 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.00197 0.00953 0.01223 0.00311
* State 4
- Energy: -75.70731804012325
- 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.62036297238913
- 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.00187
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:00
2 -75.793118310 -9.5e-02 3.6e-02 0:00:00
3 -75.798908867 -7.2e-03 2.0e-03 0:00:00
4 -75.799401560 -5.8e-04 2.2e-04 0:00:00
5 -75.799440899 -6.4e-05 2.9e-05 0:00:00
6 -75.799446125 -9.5e-06 3.3e-06 0:00:00
7 -75.799446623 -9.0e-07 5.5e-07 0:00:00
8 -75.799446711 -2.5e-07 2.1e-07 0:00:00
9 -75.799446734 -9.2e-08 6.1e-08 0:00:00
10 -75.799446738 -2.0e-08 7.4e-09 0:00:00
11 -75.799446739 -2.1e-09 8.5e-10 0:00:00
** Convergence reached in 11 iterations
Final results
-------------
* State 1
- Energy: -76.12003281964522
- 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.80828932732828
- 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.7262299016006
- 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.71537749268838
- 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.6273041535163
- 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.00200
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.