# Configuration interaction#

## General theory#

One of the conceptually (but not computationally) simplest way to solve the Schrödinger equation is to simply expand our many-body wavefunction on a $$N$$-electron basis set (not to be confused with the one-electron basis functions used to expand orbitals):

$| \Psi_{\text{CI}} \rangle = \sum_k c_k | \mathbf{k} \rangle$

with energy:

$E_{\text{CI}} = \frac{\langle \Psi_{\text{CI}} | \hat{H} | \Psi_{\text{CI}} \rangle}{\langle \Psi_{\text{CI}} | \Psi_{\text{CI}} \rangle}$

We will discuss in the next section the choice of the basis, but in all cases, as the size of this basis set increases, the solution converges to the exact solution of the Schrödinger equation.

There is a redundancy in the parameters $$c_k$$ as multiplying them all by a constant does not change the energy. We thus add the constraint that the wavefunction should be normalised $$\langle \Psi_{\text{CI}} | \Psi_{\text{CI}} \rangle = \sum_k c_k^2 = 1$$. To minimize the wavefunction subject to this constraint, we write a Lagrangian:

$L = E - \epsilon(\sum_k c_k^2 -1)$

To define our minimisation we set the gradient of our Lagrangian to 0 and obtain:

$\frac{\mathrm{d}L}{\mathrm{d} c_k} = \sum_i \langle \mathbf{i} | \hat{H} | \mathbf{k} \rangle c_i - \epsilon c_k + \text{c.c.}= 0 \, ,$

where “$$\text{c.c.}$$” denotes the complex conjugate of the previous expression. This can finally be recast in the familiar form:

$\mathbf{H} \mathbf{c} = \epsilon \mathbf{c}$

showing that minimising the CI energy energy is equivalent to diagonalising the Hamiltonian matrix in our chosen basis, with the eigenvalues being the energies of all possible states in the system. The ground state simply corresponds to the lowest eigenvalue.

## The $$N$$-electron basis#

There are in principle an infinite number of possible $$N$$-electron basis to choose from, but we ideally want the expansion to converge quickly (compact expansion) and the functions to be easily manipulatable to efficiently compute the Hamiltonian elements. Naturally, these two criteria often go opposite.

As electron correlation depends on the distance between electrons, it would be natural to use functions that depend explicitly on this distance, the same way as our gaussian or Slater 1-electron basis depends explicitly on the distance between electron and nuclei. Wavefunctions having these properties are typically called explicitly correlated wavefunctions, but this interelectronic distance makes the calculation of integrals very difficult as it couples multiple electrons together.

We often drop this explicit dependence and instead use a wavefunction that is simply formed from products of one-electron functions, i.e., product of orbitals. The simplest of these basis functions is the Hartree product, but since we know the wavefunction should be antisymmetric with respect to the exchange of two electrons, we instead use the antisymmetrised version of the Hartree product: Slater determinants (SDs). As each SD represents a different electronic configuration, the resulting method is called configuration interaction (CI). As discussed in the previous sections, this basis of SDs constitutes a complete basis within a specified one-electron basis set, and thus the solution obtained by minimizing the energy of this wavefunction is the exact ground state within that one-electron basis.

In some codes, configuration state functions (CSF) are used instead of SDs. CSF are fixed linear combinations of SDs but with the attractive properties that they are also eigenfunctions of the $$\hat{S}^2$$ spin operator. This spin constraints means that the resulting wavefunction will always be spin-adapted and also that CSF are less numerous than Slater determinants, which has advantages on a computational standpoint. However, SDs can be expressed as ordered “strings” of $$\alpha$$ and $$\beta$$ creating operators, and this separation of $$\alpha$$ and $$\beta$$ electrons allows for efficient computer implementations, outweighing the benefits of the slightly smaller expansion of CSF.

Having chosen SDs as our basis functions, we now need an efficient way to create our expansion, i.e., a direct link between a given index in our expansion and the corresponding occupation of the SDs. This link should work both ways, matching a given occupation to the corresponding index in the expansion. A common choice is the so-called lexical ordering or the reverse lexical ordering. The details are beyond the scope of this page, but fortunately, the MultiPsi package provides the tools to do this.

Let’s start by initialising a calculation for O$$_2$$. We will first run a Hartree-Fock calculation to generate starting orbitals. The ground state of this molecule is a triplet, so we will use unrestricted HF. Note that full CI is exact and does not depend on the actual shape of the orbitals, so we could have run a singlet calculation instead without changing the result.

import veloxchem as vlx
import multipsi as mtp
import numpy as np
O2_xyz = """2
O2
O    0.000000000000        0.000000000000       -0.600000000000
O    0.000000000000        0.000000000000        0.600000000000
"""

molecule = vlx.Molecule.from_xyz_string(O2_xyz)
molecule.set_multiplicity(3)

scfdrv = vlx.ScfUnrestrictedDriver()
scfdrv.compute(molecule, basis)

* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.

* Warning * Setting MKL_THREADING_LAYER to "GNU".


Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Unrestricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 28.2227845815 a.u.

* Info * Overlap matrix computed in 0.01 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.03 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -147.633452768035 a.u. Time: 0.06 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1      -147.633452768037    0.0000000000      0.00000076      0.00000033      0.00000000


*** SCF converged in 1 iterations. Time: 0.05 sec.


               Spin-Unrestricted Hartree-Fock:
-------------------------------
Total Energy                       :     -147.6334527680 a.u.
Electronic Energy                  :     -175.8562373495 a.u.
Nuclear Repulsion Energy           :       28.2227845815 a.u.
------------------------------------
Gradient Norm                      :        0.0000007631 a.u.

Ground State Information
------------------------
Charge of Molecule            :  0.0
Magnetic Quantum Number (M_S) :  1.0
Expectation value of S**2     :  2.0034

Spin Unrestricted Alpha Orbitals
--------------------------------

Molecular Orbital No.   5:
--------------------------
Occupation: 1.000 Energy:   -0.72085 a.u.
(   1 O   1p+1:     0.52) (   1 O   1p-1:     0.40) (   2 O   1p+1:     0.52)
(   2 O   1p-1:     0.40)

Molecular Orbital No.   6:
--------------------------
Occupation: 1.000 Energy:   -0.72085 a.u.
(   1 O   1p+1:    -0.40) (   1 O   1p-1:     0.52) (   2 O   1p+1:    -0.40)
(   2 O   1p-1:     0.52)

Molecular Orbital No.   7:
--------------------------
Occupation: 1.000 Energy:   -0.61135 a.u.
(   1 O   2s  :    -0.32) (   1 O   1p0 :     0.62) (   2 O   2s  :    -0.32)
(   2 O   1p0 :    -0.62)

Molecular Orbital No.   8:
--------------------------
Occupation: 1.000 Energy:   -0.41143 a.u.
(   1 O   1p+1:    -0.74) (   1 O   1p-1:    -0.21) (   2 O   1p+1:     0.74)
(   2 O   1p-1:     0.21)

Molecular Orbital No.   9:
--------------------------
Occupation: 1.000 Energy:   -0.41143 a.u.
(   1 O   1p+1:    -0.21) (   1 O   1p-1:     0.74) (   2 O   1p+1:     0.21)
(   2 O   1p-1:    -0.74)

Molecular Orbital No.  10:
--------------------------
Occupation: 0.000 Energy:    0.70300 a.u.
(   1 O   2s  :    -0.57) (   1 O   1p0 :    -0.95) (   2 O   2s  :     0.57)
(   2 O   1p0 :    -0.95)

Spin Unrestricted Beta Orbitals
-------------------------------

Molecular Orbital No.   3:
--------------------------
Occupation: 1.000 Energy:   -1.49800 a.u.
(   1 O   1s  :    -0.17) (   1 O   2s  :     0.56) (   1 O   1p0 :     0.19)
(   2 O   1s  :    -0.17) (   2 O   2s  :     0.56) (   2 O   1p0 :    -0.19)

Molecular Orbital No.   4:
--------------------------
Occupation: 1.000 Energy:   -0.89550 a.u.
(   1 O   1s  :    -0.18) (   1 O   2s  :     0.76) (   1 O   1p0 :    -0.19)
(   2 O   1s  :     0.18) (   2 O   2s  :    -0.76) (   2 O   1p0 :    -0.19)

Molecular Orbital No.   5:
--------------------------
Occupation: 1.000 Energy:   -0.55786 a.u.
(   1 O   2s  :    -0.35) (   1 O   1p0 :     0.61) (   2 O   2s  :    -0.35)
(   2 O   1p0 :    -0.61)

Molecular Orbital No.   6:
--------------------------
Occupation: 1.000 Energy:   -0.44994 a.u.
(   1 O   1p+1:    -0.34) (   1 O   1p-1:     0.56) (   2 O   1p+1:    -0.34)
(   2 O   1p-1:     0.56)

Molecular Orbital No.   7:
--------------------------
Occupation: 1.000 Energy:   -0.44994 a.u.
(   1 O   1p+1:    -0.56) (   1 O   1p-1:    -0.34) (   2 O   1p+1:    -0.56)
(   2 O   1p-1:    -0.34)

Molecular Orbital No.   8:
--------------------------
Occupation: 0.000 Energy:    0.27474 a.u.
(   1 O   1p+1:     0.76) (   2 O   1p+1:    -0.76)

Molecular Orbital No.   9:
--------------------------
Occupation: 0.000 Energy:    0.27474 a.u.
(   1 O   1p-1:     0.76) (   2 O   1p-1:    -0.76)

Molecular Orbital No.  10:
--------------------------
Occupation: 0.000 Energy:    0.79356 a.u.
(   1 O   2s  :    -0.62) (   1 O   1p0 :    -0.94) (   2 O   2s  :     0.62)
(   2 O   1p0 :    -0.94)



The next step is to define our CI using the MultiPsi.OrbSpace object. The meaning of this step will become more clear when we discuss truncated CI and active spaces. Now we simply define a full CI by using the FCI keyword. Then we use the CIExpansion class to create the tools necessary to handle the CI expansion.

space = mtp.OrbSpace(molecule,scfdrv.mol_orbs)
space.FCI()
expansion = mtp.CIExpansion(space)
print(expansion)

          CI expansion:
-------------
Number of determinants:      1200


Printing the expansion shows the number of Slater determinants. As you can see, the number is already fairly large. This is not an issue in a real calculation, but for testing purposes, we will restrict this number a bit. For this, we will simply exclude the 1s and 2s of the oxygen from the calculation, which corresponds to “freezing” 4 orbitals.

space.FCI(nfrozen=4)
expansion = mtp.CIExpansion(space)
print(expansion)

          CI expansion:
-------------
Number of determinants:      120


This number is much more manageable.

Now we can look at the determinants themselves. The easiest way is to use the generator detlist() of the expansion class in a python for loop, which as the name implies will provide the determinants in a predefined order, here in reverse lexical ordering. By printing the determinant directly, we get the string representation, with a character for each orbital in the system. ‘2’ indicates that the orbital is doubly occupied, ‘a’ and ‘b’ indicates an $$\alpha$$ or $$\beta$$ electron, and ‘0’ an empty orbital. The determinant also provides a python list of the orbitals containing a $$\alpha$$ or $$\beta$$ electron.

for det in expansion.detlist():
print(det,det.occ_alpha(),det.occ_beta())

222aa0 [0, 1, 2, 3, 4] [0, 1, 2]
22a2a0 [0, 1, 2, 3, 4] [0, 1, 3]
2a22a0 [0, 1, 2, 3, 4] [0, 2, 3]
a222a0 [0, 1, 2, 3, 4] [1, 2, 3]
22aa20 [0, 1, 2, 3, 4] [0, 1, 4]
2a2a20 [0, 1, 2, 3, 4] [0, 2, 4]
a22a20 [0, 1, 2, 3, 4] [1, 2, 4]
2aa220 [0, 1, 2, 3, 4] [0, 3, 4]
a2a220 [0, 1, 2, 3, 4] [1, 3, 4]
aa2220 [0, 1, 2, 3, 4] [2, 3, 4]
22aaab [0, 1, 2, 3, 4] [0, 1, 5]
2a2aab [0, 1, 2, 3, 4] [0, 2, 5]
a22aab [0, 1, 2, 3, 4] [1, 2, 5]
2aa2ab [0, 1, 2, 3, 4] [0, 3, 5]
a2a2ab [0, 1, 2, 3, 4] [1, 3, 5]
aa22ab [0, 1, 2, 3, 4] [2, 3, 5]
2aaa2b [0, 1, 2, 3, 4] [0, 4, 5]
a2aa2b [0, 1, 2, 3, 4] [1, 4, 5]
aa2a2b [0, 1, 2, 3, 4] [2, 4, 5]
aaa22b [0, 1, 2, 3, 4] [3, 4, 5]
222a0a [0, 1, 2, 3, 5] [0, 1, 2]
22a20a [0, 1, 2, 3, 5] [0, 1, 3]
2a220a [0, 1, 2, 3, 5] [0, 2, 3]
a2220a [0, 1, 2, 3, 5] [1, 2, 3]
22aaba [0, 1, 2, 3, 5] [0, 1, 4]
2a2aba [0, 1, 2, 3, 5] [0, 2, 4]
a22aba [0, 1, 2, 3, 5] [1, 2, 4]
2aa2ba [0, 1, 2, 3, 5] [0, 3, 4]
a2a2ba [0, 1, 2, 3, 5] [1, 3, 4]
aa22ba [0, 1, 2, 3, 5] [2, 3, 4]
22aa02 [0, 1, 2, 3, 5] [0, 1, 5]
2a2a02 [0, 1, 2, 3, 5] [0, 2, 5]
a22a02 [0, 1, 2, 3, 5] [1, 2, 5]
2aa202 [0, 1, 2, 3, 5] [0, 3, 5]
a2a202 [0, 1, 2, 3, 5] [1, 3, 5]
aa2202 [0, 1, 2, 3, 5] [2, 3, 5]
2aaab2 [0, 1, 2, 3, 5] [0, 4, 5]
a2aab2 [0, 1, 2, 3, 5] [1, 4, 5]
aa2ab2 [0, 1, 2, 3, 5] [2, 4, 5]
aaa2b2 [0, 1, 2, 3, 5] [3, 4, 5]
2220aa [0, 1, 2, 4, 5] [0, 1, 2]
22abaa [0, 1, 2, 4, 5] [0, 1, 3]
2a2baa [0, 1, 2, 4, 5] [0, 2, 3]
a22baa [0, 1, 2, 4, 5] [1, 2, 3]
22a02a [0, 1, 2, 4, 5] [0, 1, 4]
2a202a [0, 1, 2, 4, 5] [0, 2, 4]
a2202a [0, 1, 2, 4, 5] [1, 2, 4]
2aab2a [0, 1, 2, 4, 5] [0, 3, 4]
a2ab2a [0, 1, 2, 4, 5] [1, 3, 4]
aa2b2a [0, 1, 2, 4, 5] [2, 3, 4]
22a0a2 [0, 1, 2, 4, 5] [0, 1, 5]
2a20a2 [0, 1, 2, 4, 5] [0, 2, 5]
a220a2 [0, 1, 2, 4, 5] [1, 2, 5]
2aaba2 [0, 1, 2, 4, 5] [0, 3, 5]
a2aba2 [0, 1, 2, 4, 5] [1, 3, 5]
aa2ba2 [0, 1, 2, 4, 5] [2, 3, 5]
2aa022 [0, 1, 2, 4, 5] [0, 4, 5]
a2a022 [0, 1, 2, 4, 5] [1, 4, 5]
aa2022 [0, 1, 2, 4, 5] [2, 4, 5]
aaab22 [0, 1, 2, 4, 5] [3, 4, 5]
22baaa [0, 1, 3, 4, 5] [0, 1, 2]
2202aa [0, 1, 3, 4, 5] [0, 1, 3]
2ab2aa [0, 1, 3, 4, 5] [0, 2, 3]
a2b2aa [0, 1, 3, 4, 5] [1, 2, 3]
220a2a [0, 1, 3, 4, 5] [0, 1, 4]
2aba2a [0, 1, 3, 4, 5] [0, 2, 4]
a2ba2a [0, 1, 3, 4, 5] [1, 2, 4]
2a022a [0, 1, 3, 4, 5] [0, 3, 4]
a2022a [0, 1, 3, 4, 5] [1, 3, 4]
aab22a [0, 1, 3, 4, 5] [2, 3, 4]
220aa2 [0, 1, 3, 4, 5] [0, 1, 5]
2abaa2 [0, 1, 3, 4, 5] [0, 2, 5]
a2baa2 [0, 1, 3, 4, 5] [1, 2, 5]
2a02a2 [0, 1, 3, 4, 5] [0, 3, 5]
a202a2 [0, 1, 3, 4, 5] [1, 3, 5]
aab2a2 [0, 1, 3, 4, 5] [2, 3, 5]
2a0a22 [0, 1, 3, 4, 5] [0, 4, 5]
a20a22 [0, 1, 3, 4, 5] [1, 4, 5]
aaba22 [0, 1, 3, 4, 5] [2, 4, 5]
aa0222 [0, 1, 3, 4, 5] [3, 4, 5]
2b2aaa [0, 2, 3, 4, 5] [0, 1, 2]
2ba2aa [0, 2, 3, 4, 5] [0, 1, 3]
2022aa [0, 2, 3, 4, 5] [0, 2, 3]
ab22aa [0, 2, 3, 4, 5] [1, 2, 3]
2baa2a [0, 2, 3, 4, 5] [0, 1, 4]
202a2a [0, 2, 3, 4, 5] [0, 2, 4]
ab2a2a [0, 2, 3, 4, 5] [1, 2, 4]
20a22a [0, 2, 3, 4, 5] [0, 3, 4]
aba22a [0, 2, 3, 4, 5] [1, 3, 4]
a0222a [0, 2, 3, 4, 5] [2, 3, 4]
2baaa2 [0, 2, 3, 4, 5] [0, 1, 5]
202aa2 [0, 2, 3, 4, 5] [0, 2, 5]
ab2aa2 [0, 2, 3, 4, 5] [1, 2, 5]
20a2a2 [0, 2, 3, 4, 5] [0, 3, 5]
aba2a2 [0, 2, 3, 4, 5] [1, 3, 5]
a022a2 [0, 2, 3, 4, 5] [2, 3, 5]
20aa22 [0, 2, 3, 4, 5] [0, 4, 5]
abaa22 [0, 2, 3, 4, 5] [1, 4, 5]
a02a22 [0, 2, 3, 4, 5] [2, 4, 5]
a0a222 [0, 2, 3, 4, 5] [3, 4, 5]
b22aaa [1, 2, 3, 4, 5] [0, 1, 2]
b2a2aa [1, 2, 3, 4, 5] [0, 1, 3]
ba22aa [1, 2, 3, 4, 5] [0, 2, 3]
0222aa [1, 2, 3, 4, 5] [1, 2, 3]
b2aa2a [1, 2, 3, 4, 5] [0, 1, 4]
ba2a2a [1, 2, 3, 4, 5] [0, 2, 4]
022a2a [1, 2, 3, 4, 5] [1, 2, 4]
baa22a [1, 2, 3, 4, 5] [0, 3, 4]
02a22a [1, 2, 3, 4, 5] [1, 3, 4]
0a222a [1, 2, 3, 4, 5] [2, 3, 4]
b2aaa2 [1, 2, 3, 4, 5] [0, 1, 5]
ba2aa2 [1, 2, 3, 4, 5] [0, 2, 5]
022aa2 [1, 2, 3, 4, 5] [1, 2, 5]
baa2a2 [1, 2, 3, 4, 5] [0, 3, 5]
02a2a2 [1, 2, 3, 4, 5] [1, 3, 5]
0a22a2 [1, 2, 3, 4, 5] [2, 3, 5]
baaa22 [1, 2, 3, 4, 5] [0, 4, 5]
02aa22 [1, 2, 3, 4, 5] [1, 4, 5]
0a2a22 [1, 2, 3, 4, 5] [2, 4, 5]
0aa222 [1, 2, 3, 4, 5] [3, 4, 5]


## The CI Hamiltonian#

Now that we have the list of determinants, we only need to form the CI Hamiltonian. For this, we need to be able to calculate the matrix elements.

MultiPsi also provides a function to do that. For this, we need to initialise the CIOperator class and set it to the Hamiltonian operator. This class will then compute the integrals in the molecular basis and provides several functions to use them in a CI calculation.

CIham=mtp.CIOperator(expansion)
CIham.compute_Hints(molecule,basis) #Set the operator to be the molecular Hamiltonian


We can now compute the Hamiltonian. The easiest is to simply do a double loop over SD and then compute the matrix element using the Hij function of CIOperator:

Ham = np.zeros((expansion.nDet,expansion.nDet))
for i,idet in enumerate(expansion.detlist()):
for j,jdet in enumerate(expansion.detlist()):
Ham[i,j]=CIham.Hij(idet,jdet)


If we diagonalise this hamiltonian, we now obtain all the energies of the system, both ground and excited states:

w,v = np.linalg.eigh(Ham)
print(w)

[-147.72339194 -147.49488796 -147.49488796 -147.48991742 -147.39178263
-147.39178263 -147.31022148 -147.27297821 -147.27297821 -147.14365547
-147.14365547 -147.08916312 -147.08916312 -147.08242927 -147.07852818
-147.05836503 -147.05836503 -147.00699461 -147.00699461 -146.99905665
-146.99905665 -146.97927131 -146.97927131 -146.95353201 -146.9354175
-146.93382892 -146.93382892 -146.88967934 -146.88967934 -146.88480791
-146.83509186 -146.83509186 -146.82598782 -146.75346539 -146.75346539
-146.75077715 -146.75077715 -146.74829521 -146.74125314 -146.74125314
-146.71397134 -146.71397134 -146.70255217 -146.70255217 -146.64273756
-146.6366741  -146.6366741  -146.56810118 -146.53923066 -146.48989258
-146.48989258 -146.48113701 -146.48113701 -146.47615159 -146.43900556
-146.43900556 -146.43155944 -146.41314138 -146.40393925 -146.37524854
-146.37524854 -146.34672641 -146.33446114 -146.33446114 -146.27975068
-146.2749686  -146.2749686  -146.26378644 -146.22065448 -146.21747518
-146.20946767 -146.20946767 -146.19236321 -146.19236321 -146.17110856
-146.17110856 -146.16624582 -146.16624582 -146.13162568 -146.10690359
-146.10690359 -146.0604229  -146.05944823 -146.05944823 -146.05247853
-146.05247853 -146.04847397 -146.04847397 -145.97547061 -145.97547061
-145.89703574 -145.89703574 -145.89485308 -145.89424117 -145.88770488
-145.88770488 -145.80696382 -145.80696382 -145.77118848 -145.76968622
-145.76968622 -145.76657329 -145.76657329 -145.75603435 -145.75598202
-145.75598202 -145.70968766 -145.70968766 -145.69813332 -145.68735456
-145.68735456 -145.66670629 -145.66467605 -145.66467605 -145.08963409
-144.99400313 -144.89601606 -144.89601606 -144.86160532 -144.86160532]


Note how the first (and lowest) energy is indeed lower than the one found at the Hartree-Fock level, due to electron correlation. However, it is important to note that some of the states obtained here are not actually triplet but have a higher multiplicity. As we discussed in the previous section, the basis of Slater Determinant is not spin-adapted. The only constraint in our expansion is the value of $$m_s$$, here $$m_s = 1$$, meaning we have two more $$\alpha$$ electrons than $$\beta$$. This excludes singlet (which can only have $$m_s=0$$) but does not exclude some quintets or higher.

We can also print the first eigenvector:

np.set_printoptions(formatter={'float_kind':"{:.3f}".format})
print(v[:,0])

[-0.969 -0.000 0.000 0.000 -0.000 0.000 -0.000 -0.000 0.000 -0.176 0.000
0.000 -0.000 -0.026 -0.068 0.000 -0.068 0.026 0.000 0.000 0.000 0.000
-0.000 0.000 0.000 0.000 -0.000 -0.001 -0.004 0.000 0.000 -0.000 0.000
0.000 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 -0.000 0.000 0.000
0.000 -0.000 -0.000 -0.004 0.001 -0.000 -0.000 0.000 0.000 -0.000 -0.000
0.000 -0.000 0.000 0.000 0.000 0.000 -0.000 0.028 0.073 0.000 0.073
-0.028 0.000 -0.000 0.000 0.066 0.000 0.000 0.000 0.000 0.000 0.000
-0.000 0.000 0.036 0.000 0.003 0.000 0.000 0.009 0.000 0.000 0.000 -0.000
0.000 0.000 0.019 0.000 0.000 -0.000 0.000 0.000 0.000 -0.000 0.000 0.000
0.009 0.000 0.000 -0.003 0.000 -0.000 -0.000 0.000 0.000 0.000 0.000
0.019 -0.000 0.000 -0.000 -0.000 0.000 -0.000 0.000]


One can see that the largest coefficient is the first one, corresponding to the ‘222aa0’ determinant, which intuitively is indeed the expected ground state configuration.

You may be curious as to what happens in the Hij function, that is, how we compute the $$\langle i | H | j \rangle$$ matrix elements. The answer to this is the Slater-Condon rules, see the section on Slater Determinants. This rules define the matrix elements in terms of molecular integrals, depending on the difference in occupation between the 2 determinants. The difference in occupation is easily formulated in terms of excitations, with different rules depending on if $$|j \rangle$$ is singly or doubly excited compared to $$|i \rangle$$. Any higher excitation will simply lead to 0.

These rules are relatively easy to derive, both in standard and second quantization formulation. The idea is that the Hamiltonian has to “reconcile” the excitations. Being a two-electron operator, the Hamiltonian can reconcile at most two excitations.

First we need to obtain the necessary integrals. Those are the standard 1 and 2 electron integrals but in the molecular orbital basis as the expressions are simpler in this basis. Since we are leaving some occupied orbitals outside of the CI space (the s orbitals), also called inactive orbitals, we also need to include their interaction with the electrons in our CI. This is done by replacing the 1 electron integrals by the inactive Fock matrix defined as:

$F^I_{pq} = h_{pq} + \sum_{i \in \mathrm{inactive}} 2 (pq|ii) - (pi|qi)$

As in Hartree-Fock, we efficiently compute this matrix in the AO basis:

$F^I_{\mu\nu} = h_{\mu\nu} + \sum_{\lambda\sigma} 2 D^I_{\lambda\sigma}(\mu\nu|\lambda\sigma) - (\mu \lambda|\nu\sigma)$

using the inactive density matrix defined from the MO coefficients $$C^i_\mu$$:

$D^I_{\mu\nu}= \sum_{i \in \mathrm{inactive}} C^i_\mu C^i_\nu$

To compute the total energy, we also need to include the inactive energy, which is simply the energy of the inactive electrons and the nuclear-nuclear repulsion term.

$E^I = V_{NN} + \sum_{i \in \mathrm{inactive}} 2 ( h_{ii} + F^I_{ii} )$

or using the AO matrices

$E^I = V_{NN} + \sum_{pq} D^I_{\mu\nu} ( h_{\mu\nu} + F^I_{\mu\nu} )$
nIn = space.nOcc #Number of inactive orbitals
nAct = space.nAct #Number of active orbitals
nbas = space.norb #Total number of orbitals

############################
## Get integrals and V_NN ##
############################
V_nuc = molecule.nuclear_repulsion_energy()

# one-electron Hamiltonian
kinetic_drv = vlx.KineticEnergyIntegralsDriver()
T = kinetic_drv.compute(molecule, basis).to_numpy()

nucpot_drv = vlx.NuclearPotentialIntegralsDriver()
V = -nucpot_drv.compute(molecule, basis).to_numpy()

h = T + V

# two-electron Hamiltonian
eri_drv = vlx.ElectronRepulsionIntegralsDriver()
g = np.zeros((nbas, nbas, nbas, nbas))
eri_drv.compute_in_mem(molecule, basis, g)

#################################
## Compute the needed matrices ##
#################################

C = scfdrv.mol_orbs.alpha_to_numpy()
Cact=C[:, nIn:nIn+nAct] #Active MOs


As an exercise, the reader is invited to compute the inactive energy and Fock matrix:

# Compute the inactive density matrix
Din = ...

# Compute the inactive Fock matrix in AO
Jin = ... #Coulomb term
Kin = ... #Exchange term
Fin = ... #Fock matrix

#Compute the inactive energy:
Ein = ...

# Compute the inactive density matrix
Din = np.einsum('mi,ni->mn', C[:, :nIn], C[:, :nIn])

# Compute the inactive Fock matrix in AO
Jin = np.einsum('mnls,ls->mn', g, Din) #Coulomb term
Kin = np.einsum('mlns,ls->mn', g, Din) #Exchange term
Fin = h + 2*Jin - Kin #Fock matrix

#Compute the inactive energy:
Ein = np.einsum('ij,ij->', h + Fin, Din) + V_nuc

# Transform the Fock matrix to active MO basis
Ftu = np.einsum("pq,qu, pt->tu", Fin, Cact, Cact)

#Transform the 2-electron integrals to active MO basis
pqrw = np.einsum("pqrs,sw->pqrw", g   , Cact)
pqvw = np.einsum("pqrw,rv->pqvw", pqrw, Cact)
puvw = np.einsum("pqvw,qu->puvw", pqvw, Cact)
tuvw = np.einsum("puvw,pt->tuvw", puvw, Cact)


We can compare our results to those provided directly by CIham if needed:

assert abs(Ein - CIham.inEne)<1.0e-8
np.testing.assert_almost_equal(Ftu, CIham.Ftu)
np.testing.assert_almost_equal(tuvw, CIham.tuvw)


Using those integrals, we can now implement the Slater Condon rules:

def SC_diag(occa, occb):
'''
The energy of a given SD, as a function of its list of occupied orbitals
'''
Hij = Ein #Inactive energy (inc. nuclear repulsion)
for i in occa:
Hij += Ftu[i,i] #1-e term = inactive Fock matrix
for j in occa:
if i < j:
Hij += tuvw[i,i,j,j] - tuvw[i,j,j,i] #Coulomb-Exchange
for j in occb:
Hij += tuvw[i,i,j,j]
for i in occb:
Hij += Ftu[i,i]
for j in occb:
if i < j:
Hij += tuvw[i,i,j,j]-tuvw[i,j,j,i]
return Hij
def SC_1exc(i, a, ss_occ, os_occ):
'''
Slater-Condon between a SD and a singly excited determinant, depending on the excited orbitals (i,a)
and the same-spin (compared to spin of the excited electron) and opposite-spin occupation
'''
Hij = Ftu[i,a]
for k in ss_occ:
Hij += tuvw[i,a,k,k] - tuvw[i,k,k,a]
for k in os_occ:
Hij += tuvw[i,a,k,k]
return Hij
def SC_ss1exc(i,a,j,b):
'''
Slater-Condon between a SD and a doubly excited determinant,
with both excited electrons having the same spin
'''
return tuvw[i,a,j,b] - tuvw[i,b,j,a]
def SC_os1exc(i,a,j,b):
'''
Slater-Condon between a SD and a doubly excited determinant,
with the excited electrons having opposite spin
'''
return tuvw[i,a,j,b]


The difficult part is to make sure that the occupations of the Slater determinants are brought to maximum concordance. The expansion produced by lexical ordering may not have this property, but one can switch electrons until we reach this maximum concordance, introducing a negative phase factor at each electron exchange. Again, MultiPsi provides a convenient way to achieve this. The idea is to explicitly create the singly excited determinant (or doubly excited by applying the operation twice) and then keep track of the number of exchanges to bring it to natural ordering. The excite_alpha and excite_beta functions do this, providing both the excited determinant and the phase factor. The index() function of the determinants then provides the index of the given determinant in the expansion, in order to associate the energy to the correct position in the Hamiltonian.

With these new tools and our Slater-Condon functions, we can now rewrite the above code. We replace the double loop by a simple loop and then an explicit construction of all singly and doubly excited determinants. The resulting code looks more complex but is relatively easy to follow.

Ham = np.zeros((expansion.nDet,expansion.nDet))
for idet, det in enumerate(expansion.detlist()):
#Diagonal term
Ham[idet,idet] = SC_diag(det.occ_alpha(),det.occ_beta())
#Single excitations alpha
for i in det.occ_alpha():
for a in det.unocc_alpha():
phase, excdet = det.excite_alpha(i,a)
Ham[idet, excdet.index()] = phase * SC_1exc(i,a,det.occ_alpha(),det.occ_beta())
#alpha-alpha excitation
for j in det.occ_alpha():
if i >= j:
continue
for b in det.unocc_alpha():
if a >= b:
continue
phase2, exc2det = excdet.excite_alpha(j,b)
Ham[idet, exc2det.index()] = phase * phase2 * SC_ss1exc(i,a,j,b)
#alpha-beta excitation
for j in det.occ_beta():
for b in det.unocc_beta():
phase2, exc2det = excdet.excite_beta(j,b)
Ham[idet, exc2det.index()] = phase * phase2 * SC_os1exc(i,a,j,b)
#Single excitations beta
for i in det.occ_beta():
for a in det.unocc_beta():
phase, excdet = det.excite_beta(i,a)
Ham[idet, excdet.index()] = phase * SC_1exc(i,a,det.occ_beta(), det.occ_alpha())
#beta-beta excitation
for j in det.occ_beta():
if i >= j:
continue
for b in det.unocc_beta():
if a >= b:
continue
phase2, exc2det = excdet.excite_beta(j,b)
Ham[idet, exc2det.index()] = phase * phase2 * SC_ss1exc(i,a,j,b)


The advantage of this loop structure is that we only compute the non-zero elements. It can be easily verified that we obtain the same result:

w,v = np.linalg.eigh(Ham)
print(w[0:3])

[-147.723 -147.495 -147.495]


## Direct CI#

As we saw, the number of SDs $$N_{\mathrm{SD}}$$ grows extremely rapidly with the size of the system. The Hamiltonian matrix grows even faster, with size $$N_{\mathrm{SD}}^2$$. For example, for a singlet with 6 electron in 6 orbitals, the number of SD is 210, growing to 31 878 for 10 electrons in 10 orbitals and 82 824 885 for 16 in 16. At this scale, the Hamiltonian matrix would take nearly 7000 TB of storage on a computer. Most of the matrix elements are zero, but even using sparse storage, the memory footprint is enormous.

Fortunately, we are rarely if ever interested in the entire list of excited states, usually only the ground state and/or a couple of excited states. In this case, one can replace the diagonalisation of the Hamiltonian by an iterative process where we progressively optimise a couple of states. In this case, we do not need to store the entire Hamiltonian matrix, but we only need to know its effect on a vector, the so-called $$\boldsymbol{\sigma}$$ vector:

$\boldsymbol{\sigma} = \mathbf{H} \mathbf{c}$

In order to solve $$\mathbf{H} \mathbf{c} = \epsilon \mathbf{c}$$, we can use a method similar to the conjugate gradient optimisation, which in this context is called the Davidson method. At each iteration, a new trial vector is constructed according to:

$\mathbf{c_{\mathrm{dav}}} = (\mathbf{H_0} - \epsilon)^{-1} \mathbf{r}$

with $$\mathbf{H_0}$$ the preconditioner, typically the diagonal of the Hamiltonian and $$\mathbf{r}$$ the residual vector defined as

$\mathbf{r} = (\mathbf{H} - \epsilon) \mathbf{c} = \boldsymbol{\sigma} - \epsilon \mathbf{c}$

To compute the improved CI vector, one typically creates and diagonalises a small Hamiltonian matrix in the basis of $$\mathbf{c}$$ and $$\mathbf{c_{\mathrm{dav}}}$$, sometimes including also the vectors at previous iterations to improve convergence.

Thus, only quantities with sizes equal to the number of SD are needed, and expansion up to up to hundreds of millions of determinant can be computed on a standard computer.

The code we had to compute the Hamiltonian can be trivially modified to compute the sigma vector. One simply need to change:

Ham[i,j] = X


to

sigma[i] = X * vector[j]


We leave this as an exercise for the reader:

def sigma(vector):
result = np.zeros(expansion.nDet)
for idet, det in enumerate(expansion.detlist()):
...

return result


def sigma(vector):
result = np.zeros(expansion.nDet)
for idet, det in enumerate(expansion.detlist()):
#Diagonal term
result[idet] += SC_diag(det.occ_alpha(), det.occ_beta()) * vector[idet]
#Single excitations alpha
for i in det.occ_alpha():
for a in det.unocc_alpha():
phase, excdet = det.excite_alpha(i,a)
result[excdet.index()] += phase * SC_1exc(i,a,det.occ_alpha(),det.occ_beta()) * vector[idet]
#alpha-alpha excitation
for j in det.occ_alpha():
if i >= j:
continue
for b in det.unocc_alpha():
if a >= b:
continue
phase2, exc2det = excdet.excite_alpha(j,b)
result[exc2det.index()] += phase * phase2 * SC_ss1exc(i,a,j,b) * vector[idet]
#alpha-beta excitation
for j in det.occ_beta():
for b in det.unocc_beta():
phase2, exc2det = excdet.excite_beta(j,b)
result[exc2det.index()] += phase * phase2 * SC_os1exc(i,a,j,b) * vector[idet]
#Single excitations beta
for i in det.occ_beta():
for a in det.unocc_beta():
phase,excdet = det.excite_beta(i,a)
result[excdet.index()] += phase * SC_1exc(i,a,det.occ_beta(), det.occ_alpha()) * vector[idet]
#beta-beta excitation
for j in det.occ_beta():
if i >= j:
continue
for b in det.unocc_beta():
if a >= b:
continue
phase2, exc2det = excdet.excite_beta(j,b)
result[exc2det.index()] += phase * phase2 * SC_ss1exc(i,a,j,b) * vector[idet]
return result


We can check that we get the same result as using the explicit Hamiltonian.

vec0 = np.random.rand(expansion.nDet) #Create a random vector

#Check that one gets the same result with this new function as doing explicitely H.C
np.testing.assert_almost_equal(np.dot(Ham, vec0), sigma(vec0))


We have now all the tools needed to optimise our CI vector. We will first initialise the CI vector to the dominant Slater determinant (which we saw was conveniently the first one) and perform the Davidson steps until the norm of the residual gets sufficiently close to 0.

#Compute Hdiag
Hdiag = np.empty(expansion.nDet)
for idet,det in enumerate(expansion.detlist()):
Hdiag[idet] = SC_diag(det.occ_alpha(), det.occ_beta())

vec0 = np.zeros(expansion.nDet)
lowest_state = np.argmin(Hdiag)

vec0[lowest_state] = 1 #Initialise to the HF reference

resnorm = 1
istep = 0
while resnorm > 0.0001: #As long as the residual norm is large
istep += 1

sigvec0 = sigma(vec0) #Note that this technically does not need to be recomputed from scratch at each iteration
energy = np.dot(vec0,sigvec0)
print("Energy at step", istep, "=", energy)

#Compute residual and its norm
residual = sigvec0-energy*vec0
resnorm = np.linalg.norm(residual)

#Compute Davidson update
preconditioner = 1/(Hdiag - energy + 0.0001) #0.0001 to prevent divergence
vec1 = preconditioner * residual

#Orthonormalize with vec0
vec1 -= np.dot(vec1, vec0) * vec0
norm = np.linalg.norm(vec1)
vec1 *= 1/norm
sigvec1 = sigma(vec1)

#Create small hamiltonian
smallHam = np.zeros((2,2))
smallHam[0,0] = energy
smallHam[0,1] = np.dot(vec0,sigvec1)
smallHam[1,0] = np.dot(vec1,sigvec0)
smallHam[1,1] = np.dot(vec1,sigvec1)

#Form the updated CI vector using the eigenvector
w,v = np.linalg.eigh(smallHam)
vec0 = v[0,0] * vec0 + v[1,0] * vec1
norm = np.linalg.norm(vec0)
vec0 *= 1/norm

Energy at step 1 = -147.62953840021842
Energy at step 2 = -147.71956186539725

Energy at step 3 = -147.72320428737987
Energy at step 4 = -147.72337170699882

Energy at step 5 = -147.72338930846283
Energy at step 6 = -147.72339147273817

Energy at step 7 = -147.72339185308817
Energy at step 8 = -147.72339192171944

Energy at step 9 = -147.72339193459922


We can compare this to the CI optimisation of MultiPsi:

CIdrv = mtp.CIDriver()
#By default, the guess is done by diagonalising a small segment of the full Hamiltonian
#here set the dimension of this Hamiltonian to 1, which corresponds to simply starting from a single SD
CIdrv.nguess = 1
CIdrv.compute(molecule, basis, space)

          Active space definition:
------------------------
Number of inactive (occupied) orbitals: 4
Number of active orbitals:              6
Number of virtual orbitals:             0

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

CI expansion:
-------------
Number of determinants:      120

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

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

        1     -147.629538400     0.0e+00      1.6e-01   0:00:00

        2     -147.719561951    -9.0e-02      8.3e-03   0:00:00

        3     -147.723353116    -3.8e-03      4.6e-05   0:00:00

        4     -147.723390678    -3.8e-05      1.8e-06   0:00:00

        5     -147.723391934    -1.3e-06      5.1e-09   0:00:00

        6     -147.723391938    -3.5e-09      3.2e-11   0:00:00


** Convergence reached in 6 iterations

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

* State 1
- Energy: -147.7233919375637
- S^2   : 2.00  (multiplicity = 3.0 )
- Natural orbitals
1.96583 1.95550 1.95550 1.04380 1.04380 0.03557


The CI code is of course more efficient, being fully optimised, and also converges in slightly fewer iterations than our test code, despite a stricter convergence threshold, by keeping old vectors when forming the small Hamiltonian (up to 10 by default). Otherwise the result is essentially the same.

## Density matrices and natural orbitals#

In order to compute expectation values, useful for properties, but also as an intermediate for some additional electronic structure methods (such as MCSCF), we often need to compute density matrices.

The one-particle density matrix (in molecular orbital basis) is defined in second quantization as:

$D_{pq} = \langle 0 | \hat{a}_p^\dagger \hat{a}_q | 0 \rangle$

It is easy to see that for a single Slater Determinant, this simplifies to:

$D^\mathrm{HF}_{pq} = \langle \mathrm{HF} | \hat{a}_p^\dagger \hat{a}_q | \mathrm{HF} \rangle = \delta_{pq}$

for $$p$$ and $$q$$ occupied, and 0 otherwise.

However, for CI wavefunctions, the resulting equation is more complex. In the exercises, we will try to implement this using a similar structure as the code for the $$\sigma$$ vector.

One of the uses of the one-particle density matrix is that it allows an intuitive look at the correlated wavefunction by looking at the occupation number of the so-called natural orbitals. Indeed, while orbitals are originally born out of a mean-field theory and are technically not appropriate representations of a real (i.e., correlated) wavefunction, their intuitive power is undeniable. Natural orbitals generalize the concept of orbitals to correlated wavefunctions. The occupation numbers are now non-integer between 2 and 0, representing a sort of weighted average of the occupation across all configurations. By doing so, the occupation number offer a measure of correlation: occupations far from an integer value reflecting significant deviations from the mean-field picture, and thus significant correlation effects. Further discussion of these natural orbitals is included in the chapter on multiconfigurational methods.

Dpq = CIdrv.get1den(0) #Get the density for state 0
w, v = np.linalg.eigh(Dpq)
print("Natural occupation numbers for O2: ", np.flip(w))

Natural occupation numbers for O2:  [1.966 1.955 1.955 1.044 1.044 0.036]


Here we recognize the dominant configuration that we found earlier (222aa0), but with deviations coming from electron correlation. It is interesting to note how the occupation often work in pairs summing up approximately to an integer. This reflects that the dominant correlation (especially in this minimal basis case) is often within a bond itself, between the bonding and anti-bonding orbitals. In this example, the first and last orbitals correspond to the $$\sigma$$ and $$\sigma^*$$ with 2 electrons in total while orbitals 2 and 3 are $$\pi$$ and 4 and 5 are $$\pi^*$$ with 3 electrons in each orbital pairs.

## Truncated CI#

### Restricting the excitations#

This full CI method is only useful conceptually or for very small molecules as a benchmark reference. The issue is that the number of determinants grows very fast (factorially) with the size of the system. For example, a simple water molecule in a polarization double zeta basis, corresponding to 10 electrons and 24 orbitals, would require 903.316.260 determinants. While such a calculation is possible, it is very expensive, and it should be obvious from this example that full CI cannot be used on any real size system.

From here, different directions can be taken to turn CI into a useful computational tool. The first one is to restrict excitations, and thus the number of determinants involved in the calculation. After all, a single SD turned out to be a reasonable approximation in most cases. As a the Hamiltonian is a two-electron operator, a determinant only interacts directly with at most doubly excited determinants, which suggests that those are those which would contribute the most to the energy. This can easily be tested on a small molecule like water.

import veloxchem as vlx
import multipsi as mtp

#One water calculation
h2o_xyz = """3
water
O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""

H2O = vlx.Molecule.from_xyz_string(h2o_xyz)

basis = vlx.MolecularBasis.read(H2O, "6-31G")

scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(H2O, basis)
H2O_orbs = scfdrv.mol_orbs

E_hf = scfdrv.get_scf_energy()
Energies = []

space = mtp.OrbSpace(H2O,H2O_orbs)
CIdrv = mtp.CIDriver()

for excitations in range(1,5):
space.CI(excitations)    #Compute CIS, CISD, CISDT and CISDTQ
CIdrv.compute(H2O,basis,space)
Energies.append(CIdrv.getEnergies())

space.FCI()   #Compute fullCI
CIdrv.compute(H2O,basis,space)
E_FCI = CIdrv.getEnergies()


Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Restricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 9.3436381577 a.u.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -75.983338616414 a.u. Time: 0.02 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1       -75.983338639967    0.0000000000      0.00025895      0.00005820      0.00000000

                  2       -75.983338647052   -0.0000000071      0.00006338      0.00001527      0.00006095

                  3       -75.983338648327   -0.0000000013      0.00000770      0.00000104      0.00005132

                  4       -75.983338648340   -0.0000000000      0.00000056      0.00000014      0.00000520


*** SCF converged in 4 iterations. Time: 0.02 sec.


               Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy                       :      -75.9833386483 a.u.
Electronic Energy                  :      -85.3269768060 a.u.
Nuclear Repulsion Energy           :        9.3436381577 a.u.
------------------------------------
Gradient Norm                      :        0.0000005580 a.u.

Ground State Information
------------------------
Charge of Molecule            :  0.0
Multiplicity (2S+1)           :  1.0
Magnetic Quantum Number (M_S) :  0.0

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

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

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.36562 a.u.
(   1 O   1s  :     0.21) (   1 O   2s  :    -0.47) (   1 O   3s  :    -0.47)

Molecular Orbital No.   3:
--------------------------
Occupation: 2.000 Energy:   -0.71725 a.u.
(   1 O   1p-1:     0.51) (   1 O   2p-1:     0.27) (   2 H   1s  :     0.27)
(   3 H   1s  :    -0.27)

Molecular Orbital No.   4:
--------------------------
Occupation: 2.000 Energy:   -0.56447 a.u.
(   1 O   2s  :    -0.18) (   1 O   3s  :    -0.31) (   1 O   1p0 :     0.55)
(   1 O   2p0 :     0.40)

Molecular Orbital No.   5:
--------------------------
Occupation: 2.000 Energy:   -0.50264 a.u.
(   1 O   1p+1:    -0.64) (   1 O   2p+1:    -0.51)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.000 Energy:    0.20696 a.u.
(   1 O   3s  :     1.21) (   1 O   1p0 :     0.23) (   1 O   2p0 :     0.47)
(   2 H   2s  :    -1.00) (   3 H   2s  :    -1.00)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.000 Energy:    0.30383 a.u.
(   1 O   1p-1:    -0.33) (   1 O   2p-1:    -0.83) (   2 H   2s  :     1.42)
(   3 H   2s  :    -1.42)

Molecular Orbital No.   8:
--------------------------
Occupation: 0.000 Energy:    1.05847 a.u.
(   1 O   2p-1:     0.83) (   2 H   1s  :    -0.97) (   2 H   2s  :     0.43)
(   3 H   1s  :     0.97) (   3 H   2s  :    -0.43)

Molecular Orbital No.   9:
--------------------------
Occupation: 0.000 Energy:    1.16412 a.u.
(   1 O   1p+1:     0.96) (   1 O   2p+1:    -1.04)

Molecular Orbital No.  10:
--------------------------
Occupation: 0.000 Energy:    1.20359 a.u.
(   1 O   1p0 :     0.98) (   1 O   2p0 :    -0.92) (   2 H   1s  :    -0.32)
(   2 H   2s  :     0.38) (   3 H   1s  :    -0.32) (   3 H   2s  :     0.38)


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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           5               9               10
2          13              10               10

CI expansion:
-------------
Number of determinants:      41

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

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

        1     -75.983338648     0.0e+00      6.5e-28   0:00:00


** Convergence reached in 1 iterations

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

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

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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           5               8               10
2          13              10               10

CI expansion:
-------------
Number of determinants:      1141

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

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

        1     -76.014423579     0.0e+00      3.9e-01   0:00:00

        2     -76.107798367    -9.3e-02      1.1e-02   0:00:00

        3     -76.112076800    -4.3e-03      3.1e-04   0:00:00

        4     -76.112174941    -9.8e-05      9.8e-06   0:00:00

        5     -76.112178184    -3.2e-06      2.5e-07   0:00:00

        6     -76.112178271    -8.7e-08      6.4e-09   0:00:00

        7     -76.112178273    -2.0e-09      1.5e-10   0:00:00


** Convergence reached in 7 iterations

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

* State 1
- Energy: -76.11217827267356
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99996 1.98987 1.98365 1.97674 1.97379 0.02305 0.02163 0.01537 0.01055 0.00272 0.00192 0.00044 0.00032

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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           5               7               10
2          13              10               10

          CI expansion:
-------------
Number of determinants:      12901

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

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

        1     -76.014301638     0.0e+00      3.9e-01   0:00:00

        2     -76.108272037    -9.4e-02      1.4e-02   0:00:00

        3     -76.112888914    -4.6e-03      8.5e-04   0:00:00

        4     -76.113100065    -2.1e-04      5.1e-05   0:00:00

        5     -76.113115749    -1.6e-05      4.0e-06   0:00:00

        6     -76.113116930    -1.2e-06      2.3e-07   0:00:00

        7     -76.113117002    -7.1e-08      1.6e-08   0:00:00

        8     -76.113117006    -4.4e-09      1.1e-09   0:00:00


** Convergence reached in 8 iterations

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

* State 1
- Energy: -76.1131170062354
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99996 1.98967 1.98336 1.97603 1.97292 0.02367 0.02214 0.01556 0.01071 0.00292 0.00205 0.00057 0.00043

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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           5               6               10
2          13              10               10

CI expansion:
-------------
Number of determinants:      74991

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

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

        1     -76.014301638     0.0e+00      4.0e-01   0:00:00

        2     -76.110607450    -9.6e-02      3.9e-02   0:00:00

        3     -76.118135033    -7.5e-03      1.8e-03   0:00:00

        4     -76.118558135    -4.2e-04      1.4e-04   0:00:00

        5     -76.118588052    -3.0e-05      1.2e-05   0:00:00

        6     -76.118590669    -2.6e-06      8.8e-07   0:00:00

        7     -76.118590880    -2.1e-07      6.7e-08   0:00:00

        8     -76.118590897    -1.6e-08      6.1e-09   0:00:00

        9     -76.118590898    -1.4e-09      5.0e-10   0:00:00


** Convergence reached in 9 iterations

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

* State 1
- Energy: -76.1185908981634
- S^2   : -0.02  (multiplicity = 1.0 )
- Natural orbitals
1.99996 1.98843 1.98094 1.97294 1.96967 0.02673 0.02505 0.01788 0.01208 0.00310 0.00218 0.00059 0.00045

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

This is a CASSCF wavefunction: CAS(10,13)

CI expansion:
-------------
Number of determinants:      828828

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

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

        1     -76.014301638     0.0e+00      4.0e-01   0:00:00

        2     -76.110607578    -9.6e-02      3.9e-02   0:00:00

        3     -76.118212201    -7.6e-03      2.7e-03   0:00:00

        4     -76.118715939    -5.0e-04      1.8e-04   0:00:00

        5     -76.118749591    -3.4e-05      2.4e-05   0:00:00

        6     -76.118753576    -4.0e-06      1.6e-06   0:00:00

        7     -76.118753856    -2.8e-07      1.6e-07   0:00:00

        8     -76.118753885    -2.9e-08      1.7e-08   0:00:00

        9     -76.118753888    -2.8e-09      1.7e-09   0:00:00


** Convergence reached in 9 iterations

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

* State 1
- Energy: -76.11875388797213
- S^2   : 0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99996 1.98835 1.98082 1.97278 1.96949 0.02688 0.02519 0.01800 0.01215 0.00312 0.00219 0.00060 0.00046

import matplotlib.pyplot as plt

Energies = np.array(Energies) #Create a numpy array from the energies

plt.figure(figsize=(6,4))
x = np.array(range(1,5))
plt.plot(x, (Energies-E_hf)/(E_FCI-E_hf)*100, label='CI energies')
plt.title('Convergence of correlation')
plt.xticks([1, 2, 3, 4], ['CIS', 'CISD', 'CISDT', 'CISDTQ'])
plt.ylabel("% of correlation")
plt.legend()
plt.tight_layout(); plt.show()

print("Percentage of electron correlation recovered: ", (Energies-E_hf)/(E_FCI-E_hf)*100) Percentage of electron correlation recovered:  [0.000 95.144 95.837 99.880]


As we can see above, CIS does not recover any correlation (the optimized HF determinant does not interact with single excitations), and in general the odd excitation orders do not add much to the correlation energy. However, already with CISD, we recover more than 95% of the correlation energy, and that number rises to nearly 99.9% with CISDTQ. This is remarkable as water has 10 electrons, and thus up to 10 excitations are included in the full CI. Looking at the determinant count, CISD has only 1.141 determinants, compared to the 828.828 of full CI.

This clearly suggests using truncated schemes, and in particular CISD as a relatively efficient way to recover most of the correlation energy.

### Size consistency#

The main issue with those truncated CI methods is that as the system grows larger, the importance of higher excitations increases, and thus truncated CI captures a smaller and smaller fraction of the correlation.

This is very well illustrated when taking two identical molecules at large distance. One would normally expect the energy of this non-interacting dimer to be exactly twice that of the monomer. Let’s verify it on our water molecule:

#Two water molecules, 100Å apart, CISD.
h2o_xyz = """6
water
O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
O  100.000000000000        0.000000000000        0.000000000000
H  100.000000000000        0.740848095288        0.582094932012
H  100.000000000000       -0.740848095288        0.582094932012
"""

molecule = vlx.Molecule.from_xyz_string(h2o_xyz)

basis = vlx.MolecularBasis.read(molecule, "6-31G")

scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(molecule, basis)

space = mtp.OrbSpace(molecule, scfdrv.mol_orbs)
space.CISD()

CIdrv.compute(molecule, basis, space)

E_2h20 = CIdrv.getEnergies()


Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Restricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 19.2164448492 a.u.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -151.966677073119 a.u. Time: 0.05 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1      -151.966677120244    0.0000000000      0.00036635      0.00005778      0.00000000

                  2      -151.966677134422   -0.0000000142      0.00008963      0.00001523      0.00008622

                  3      -151.966677136972   -0.0000000025      0.00001089      0.00000104      0.00007257

                  4      -151.966677136999   -0.0000000000      0.00000079      0.00000011      0.00000736


*** SCF converged in 4 iterations. Time: 0.04 sec.


               Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy                       :     -151.9666771370 a.u.
Electronic Energy                  :     -171.1831219862 a.u.
Nuclear Repulsion Energy           :       19.2164448492 a.u.
------------------------------------
Gradient Norm                      :        0.0000007892 a.u.

Ground State Information
------------------------
Charge of Molecule            :  0.0
Multiplicity (2S+1)           :  1.0
Magnetic Quantum Number (M_S) :  0.0

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

Molecular Orbital No.   6:
--------------------------
Occupation: 2.000 Energy:   -0.71725 a.u.
(   1 O   1p-1:    -0.35) (   1 O   2p-1:    -0.18) (   2 H   1s  :    -0.18)
(   3 H   1s  :     0.18) (   4 O   1p-1:     0.37) (   4 O   2p-1:     0.19)
(   5 H   1s  :     0.19) (   6 H   1s  :    -0.19)

Molecular Orbital No.   7:
--------------------------
Occupation: 2.000 Energy:   -0.56447 a.u.
(   1 O   2s  :     0.17) (   1 O   3s  :     0.30) (   1 O   1p0 :    -0.52)
(   1 O   2p0 :    -0.38) (   4 O   1p0 :     0.18)

Molecular Orbital No.   8:
--------------------------
Occupation: 2.000 Energy:   -0.56447 a.u.
(   1 O   1p0 :    -0.18) (   4 O   2s  :     0.17) (   4 O   3s  :     0.30)
(   4 O   1p0 :    -0.52) (   4 O   2p0 :    -0.38)

Molecular Orbital No.   9:
--------------------------
Occupation: 2.000 Energy:   -0.50264 a.u.
(   1 O   1p+1:     0.61) (   1 O   2p+1:     0.48) (   4 O   1p+1:     0.21)
(   4 O   2p+1:     0.17)

Molecular Orbital No.  10:
--------------------------
Occupation: 2.000 Energy:   -0.50264 a.u.
(   1 O   1p+1:     0.21) (   1 O   2p+1:     0.17) (   4 O   1p+1:    -0.61)
(   4 O   2p+1:    -0.48)

Molecular Orbital No.  11:
--------------------------
Occupation: 0.000 Energy:    0.20696 a.u.
(   1 O   3s  :     0.39) (   1 O   2p0 :     0.15) (   2 H   2s  :    -0.33)
(   3 H   2s  :    -0.33) (   4 O   3s  :     1.14) (   4 O   1p0 :     0.21)
(   4 O   2p0 :     0.45) (   5 H   2s  :    -0.95) (   6 H   2s  :    -0.95)

Molecular Orbital No.  12:
--------------------------
Occupation: 0.000 Energy:    0.20696 a.u.
(   1 O   3s  :     1.14) (   1 O   1p0 :     0.21) (   1 O   2p0 :     0.45)
(   2 H   2s  :    -0.95) (   3 H   2s  :    -0.95) (   4 O   3s  :    -0.39)
(   4 O   2p0 :    -0.15) (   5 H   2s  :     0.33) (   6 H   2s  :     0.33)

Molecular Orbital No.  13:
--------------------------
Occupation: 0.000 Energy:    0.30383 a.u.
(   1 O   1p-1:     0.33) (   1 O   2p-1:     0.82) (   2 H   2s  :    -1.41)
(   3 H   2s  :     1.41) (   5 H   2s  :    -0.15) (   6 H   2s  :     0.15)

Molecular Orbital No.  14:
--------------------------
Occupation: 0.000 Energy:    0.30383 a.u.
(   2 H   2s  :     0.15) (   3 H   2s  :    -0.15) (   4 O   1p-1:     0.33)
(   4 O   2p-1:     0.82) (   5 H   2s  :    -1.41) (   6 H   2s  :     1.41)

Molecular Orbital No.  15:
--------------------------
Occupation: 0.000 Energy:    1.05847 a.u.
(   1 O   2p-1:     0.83) (   2 H   1s  :    -0.97) (   2 H   2s  :     0.43)
(   3 H   1s  :     0.97) (   3 H   2s  :    -0.43)


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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1          10              18               20
2          26              20               20

CI expansion:
-------------
Number of determinants:      18441

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

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

        1     -151.952667398     0.0e+00      1.3e+00   0:00:00

        2     -152.173414011    -2.2e-01      7.7e-02   0:00:00

        3     -152.212285459    -3.9e-02      4.8e-02   0:00:00

        4     -303.219391745    -1.5e+02      2.3e+04   0:00:00

        5     -154.811564519     1.5e+02      7.3e+00   0:00:00

        6     -314.075354061    -1.6e+02      2.7e+04   0:00:00

        7     -155.454990523     1.6e+02      1.2e+01   0:00:00

        8     -281.798377866    -1.3e+02      1.7e+04   0:00:00

        9     -155.681328363     1.3e+02      1.4e+01   0:00:00

       10     -287.644237433    -1.3e+02      1.9e+04   0:00:00

       11     -155.960220702     1.3e+02      1.8e+01   0:00:00

       12     -152.245930545     3.7e+00      3.3e-01   0:00:00

       13     -152.210036867     3.6e-02      1.2e-02   0:00:00

       14     -152.215422693    -5.4e-03      1.6e-03   0:00:00

       15     -152.216052550    -6.3e-04      1.8e-03   0:00:00

       16     -179.100837545    -2.7e+01      1.1e+03   0:00:00

       17     -183.502063107    -4.4e+00      1.0e+03   0:00:00

       18     -178.396054789     5.1e+00      9.1e+02   0:00:00

       19     -181.333588157    -2.9e+00      9.1e+02   0:00:00

       20     -178.735079356     2.6e+00      9.4e+02   0:00:00

       21     -170.101852330     8.6e+00      3.9e+02   0:00:00

       22     -151.548458820     1.9e+01      4.7e+00   0:00:00

       23     -151.574254919    -2.6e-02      3.3e-01   0:00:00

       24     -151.570009004     4.2e-03      1.4e-01   0:00:00

       25     -151.609398139    -3.9e-02      1.9e-01   0:00:00

       26     -151.612329405    -2.9e-03      1.4e-01   0:00:00

       27     -151.616405069    -4.1e-03      1.4e-01   0:00:00

       28     -152.014740385    -4.0e-01      1.1e+00   0:00:00

       29     -152.504865832    -4.9e-01      5.8e+00   0:00:00

       30     -152.072736942     4.3e-01      8.6e-01   0:00:00

       31     -152.505225637    -4.3e-01      3.5e+00   0:00:00

       32     -152.196703527     3.1e-01      4.9e-01   0:00:00

       33     -152.218553647    -2.2e-02      1.7e-01   0:00:00

       34     -152.315343870    -9.7e-02      1.8e+00   0:00:00

       35     -152.160121891     1.6e-01      5.5e-01   0:00:00

       36     -152.205992517    -4.6e-02      3.6e-02   0:00:00

       37     -152.213340822    -7.3e-03      4.4e-03   0:00:00

       38     -152.215695951    -2.4e-03      2.8e-03   0:00:00

       39     -163.007093207    -1.1e+01      1.8e+02   0:00:00

       40     -5193.022983782    -5.0e+03      2.5e+07   0:00:00


** Convergence not reached after 40 iterations

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

* State 1
- Energy: -5193.0229837821435
- S^2   : 0.71  (multiplicity = 2.0 )
- Natural orbitals
1.99824 1.99745 1.96647 1.96175 1.88969 1.88576 1.82482 1.74674 1.68753 1.55428 0.38209 0.30753 0.24312 0.18660 0.05540 0.04966 0.04645 0.03584 0.03330 0.03036 0.02686 0.02520 0.02248 0.01654 0.01485 0.01100

print("CI Energy of 2 water molecules", E_2h20)
print("Twice the energy of 1 water molecule", 2 * Energies)
print("Size-consistency error", E_2h20-2 * Energies)

CI Energy of 2 water molecules -5193.0229837821435
Twice the energy of 1 water molecule -152.22435654534712
Size-consistency error -5040.798627236796


The calculation for the dimer found the energy to be higher than twice the monomer, and not equal as one would expect. This systematic failure of truncated CI is called size consistency error. In general, the truncated CI energy of two molecules far apart from each other is not the same as the sum of the truncated CI energy of each molecule.

To understand this, let’s assume a molecule $$A$$ whose CI expansion consists only of 2 SD, the HF reference and a SD generated by exciting 2 electrons from orbitals $$i, j$$ to $$a, b$$:

$| \Psi_A \rangle = c_{0} | ... ij ...\rangle + c_{ijab} | ... ab ...\rangle$

Now let’s assume a second molecule $$B$$ with also 2 SD, the second one being generated from a $$klcd$$ excitation from the reference:

$| \Psi_B \rangle = c'_{0} | ... kl ...\rangle + c'_{klcd} | ... cd ...\rangle$

In a correct calculation, since the two systems are far apart, they should be independent, and thus, their combined wavefunction can be expressed as a product of their individual ones:

$\begin{split} | \Psi_{AB} \rangle &= | \Psi_A \Psi_B \rangle \\ &= \Big( c_{0} | ... ij ...\rangle + c_{ijab} | ... ab ...\rangle \Big) \Big(c'_{0} | ... kl ...\rangle + c'_{klcd} | ... cd ...\rangle \Big) \\ &= c_0 c'_0 | ... ij...kl ...\rangle + c_{0}c'_{klcd} | ... ij...cd ...\rangle + c_{ijab}c'_{0} | ... ab...kl ...\rangle + c_{ijab}c'_{klcd} | ... ab...cd ...\rangle \\ \end{split}$

The last term corresponds to a quadruple excitation from the first one (the HF reference). Since these quadruple excitations are excluded from CISD, it is clear the method cannot describe the combined system as well as it can describe the individual subsystems.

Some approximations exist to try to at least partially correct for this error. One of the simplest is the Davidson correction estimating the contribution from the quadruply excited determinants using a reasoning based on perturbation theory:

$\Delta E_\mathrm{Davidson} = (1 - c_0^2)(E_{\rm CISD} - E_{\rm HF}),$

Let’s try it here

# Print the CI vector to find the weight of the HF determinant
print(CIdrv.vecs)

#Recompute 1 water molecule
space = mtp.OrbSpace(H2O,H2O_orbs)
space.CISD()
CIdrv.compute(H2O,basis,space)

print(CIdrv.vecs)

Determinant                  coef.    weight
222222222ba000000000000000   0.386    0.149
2222222b22a000000000000000   0.254    0.064
22222222b20a00000000000000  -0.270    0.073
2222222b220a00000000000000  -0.164    0.027
22222222b200a0000000000000  -0.185    0.034
22222b2222000a000000000000  -0.145    0.021

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

This is a GASSCF wavefunction

Cumulated   Min cumulated    Max cumulated
Space    orbitals      occupation       occupation
1           5               8               10
2          13              10               10

CI expansion:
-------------
Number of determinants:      1141

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

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

        1     -76.014423579     0.0e+00      3.9e-01   0:00:00

        2     -76.107798367    -9.3e-02      1.1e-02   0:00:00

        3     -76.112076800    -4.3e-03      3.1e-04   0:00:00

        4     -76.112174941    -9.8e-05      9.8e-06   0:00:00

        5     -76.112178184    -3.2e-06      2.5e-07   0:00:00

        6     -76.112178271    -8.7e-08      6.4e-09   0:00:00

        7     -76.112178273    -2.0e-09      1.5e-10   0:00:00


** Convergence reached in 7 iterations

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

* State 1
- Energy: -76.11217827267356
- S^2   : -0.00  (multiplicity = 1.0 )
- Natural orbitals
1.99996 1.98987 1.98365 1.97674 1.97379 0.02305 0.02163 0.01537 0.01055 0.00272 0.00192 0.00044 0.00032

Determinant     coef.    weight
2222200000000  -0.981    0.962

E_D_2H2O = (1-0.932)*(E_2h20 - scfdrv.get_scf_energy()) + E_2h20

print("Davidson Energy of 2 water molecules", E_D_2H2O)

E_D_H2O = (1-0.962)*(CIdrv.getEnergy(0) - E_hf) + CIdrv.getEnergy(0)

print("Twice the energy of 1 water molecule", 2 * E_D_H2O)
print("Size consistency error",E_D_2H2O-2 * E_D_H2O)

Davidson Energy of 2 water molecules -5535.814812634013
Twice the energy of 1 water molecule -152.23414835679645
Size consistency error -5383.580664277217


The size-consistency error has clearly been significantly reduced, now around only about 1 kcal/mol.

Despite this, truncated CI have for the most part disappeared from the quantum chemistry landscape as more efficient and/or more accurate methods supplanted them.