Fock matrix diagonalization

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#

The canonical Hartree–Fock equation reads

\[ \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
"""

molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ", ostream=None)

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]