# Fock matrix diagonalization#

import numpy as np
import scipy

np.set_printoptions(precision=4, suppress=True, linewidth=260)

import veloxchem as vlx


## Hartree–Fock equation#

$\mathbf{F C} = \mathbf{S C} \boldsymbol{\varepsilon}$

It represents a generalized eigenvalue equation.

### Reference calculation#

Let us adopt the water molecule and perform an SCF optimization of the ground state.

h2o_xyz = """3

O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""


norb = basis.get_dimensions_of_basis(molecule)
nocc = molecule.number_of_alpha_electrons()

print("Number of contracted basis functions:", norb)
print("Number of doubly occupied molecular orbitals:", nocc)

Number of contracted basis functions: 24
Number of doubly occupied molecular orbitals: 5

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()

scf_results = scf_drv.compute(molecule, basis)


We get the Fock and overlap matrices.

F = scf_results["F_alpha"]
S = scf_results["S"]


And we obtain the orbital energies for later comparisons against those obtained from our own Fock matrix diagonalization.

print(f"Orbital energies of occupied orbitals: {scf_results['E'][:nocc]}")

Orbital energies of occupied orbitals: [-20.5482  -1.3452  -0.7058  -0.5711  -0.4946]


## Solving a generalized eigenvalue problem#

### Standard case#

The Hartree–Fock equation can be solved with a standard SciPy library routine.

epsilon, C = scipy.linalg.eigh(F, S)

print(f"Orbital energies of occupied orbitals: {epsilon[:nocc]}")

Orbital energies of occupied orbitals: [-20.5482  -1.3452  -0.7058  -0.5711  -0.4946]


We note that the orbital energies are in perfect agreement with the reference values.

### Non-standard case#

#### Linear dependence in basis set#

Let us introduce a linear dependence in the atomic orbital (AO) basis by duplicating the oxygen 1s-function. In this extended basis set, we determine the overlap and Fock matrices by duplication of the first row and column.

S_ext = np.zeros((norb + 1, norb + 1))

S_ext[1 : norb + 1, 1 : norb + 1] = S

S_ext[0, 1 : norb + 1] = S[0, :]
S_ext[1 : norb + 1, 0] = S[:, 0]

S_ext[0, 0] = S[0, 0]

F_ext = np.zeros((norb + 1, norb + 1))

F_ext[1 : norb + 1, 1 : norb + 1] = F

F_ext[0, 1 : norb + 1] = F[0, :]
F_ext[1 : norb + 1, 0] = F[:, 0]

F_ext[0, 0] = F[0, 0]


If we now try to diagonalize the Fock matrix in the extended basis set using a standard library routine we will encounter an error due to linear dependence.

try:
epsilon, C = scipy.linalg.eigh(F_ext, S_ext)
except:
print("Error: Linear dependence in the atomic orbital basis!")

Error: Linear dependence in the atomic orbital basis!


#### Diagonalizing the overlap matrix#

To address the issue of linear dependence in the AO basis, we first diagonalize the symmetric overlap matrix by means of a unitary transformation

$\mathbf{U}^{\dagger} \mathbf{S U} = \boldsymbol{\sigma}$

where $$\boldsymbol{\sigma}$$ is a diagonal matrix collecting the eigenvalues.

sigma_ext, U_ext = np.linalg.eigh(S_ext)


We expect one of the eigenvalues to be equal to zero due the linear dependence that we have introduced in the AO-basis.

np.isclose(0.0, sigma_ext)

array([ True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False])


Thereafter, we remove eigenvectors in $$\mathbf{U}$$ associated with zero eigenvalues. In real program implementations a threshold is used to detect cases of near degeneracy.

U = U_ext[:, 1:]
sigma = sigma_ext[1:]


After this operation, the shape of $$\mathbf{U}$$ becomes rectangular.

U.shape

(25, 24)


#### Orthogonalizing the AO basis#

We retrieve an orthonormal AO (OAO) basis by means of a non-unitary transformation matrix $$\mathbf{X}$$ such that

$| \overline{\chi^\mathrm{OAO}} \rangle = | \overline{\chi} \rangle \mathbf{X}$

and

$\mathbf{S}^\mathrm{OAO} = \mathbf{X}^\dagger \mathbf{S X} = \mathbf{I}$

With use of $$\mathbf{U}$$, it is straightforward to construct explicit forms of $$\mathbf{X}$$ and there exist two common choices:

symmetric form

canonical form

$$\mathbf{X} = \mathbf{U} \boldsymbol{\sigma}^{-\frac{1}{2}} \mathbf{U}^\dagger$$

$$\mathbf{X} = \mathbf{U} \boldsymbol{\sigma}^{-\frac{1}{2}}$$

We can readily verify that both expressions satisfy $$\mathbf{X}^\dagger \mathbf{S X} = \mathbf{I}$$.

# canonical transformation
X = np.einsum("ik,k->ik", U, 1 / np.sqrt(sigma))


In the OAO basis, the Hartree–Fock equation takes the form

$\mathbf{F}^\mathrm{OAO} \mathbf{C}^\mathrm{OAO} = \mathbf{C}^\mathrm{OAO} \boldsymbol{\varepsilon}$

where

$\mathbf{F}^\mathrm{OAO} = \mathbf{X}^\dagger \mathbf{F X}$
F_OAO = np.einsum("ki,kl,lj->ij", X, F_ext, X)
epsilon, C_OAO = np.linalg.eigh(F_OAO)


The expression for the associated transformation of MO coefficients is determined from the relation

$| \overline{\phi} \rangle = | \overline{\chi} \rangle \mathbf{C} = | \overline{\chi^\mathrm{OAO}} \rangle \mathbf{X}^{-1} \mathbf{C}$

We identify

$\mathbf{C}^\mathrm{OAO} = \mathbf{X}^{-1} \mathbf{C}$

or

$\mathbf{C} = \mathbf{X C}^\mathrm{OAO}$
C = np.einsum("ik,kj->ij", X, C_OAO)

print(f"Orbital energies of occupied orbitals: {epsilon[:nocc]}")

Orbital energies of occupied orbitals: [-20.5482  -1.3452  -0.7058  -0.5711  -0.4946]


We note that the orbital energies are in perfect agreement with the reference values.