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}

(1)

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()
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_alpha'][:nocc]}")

Orbital energies of occupied orbitals: [-20.54818972  -1.345205    -0.70584505  -0.57108597  -0.49456798]

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.5481899   -1.34520506  -0.7058451   -0.57108599  -0.49456809]

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}

(2)

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}

(3)

and

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

(4)

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}

(5)

where

\mathbf{F}^\mathrm{OAO} = \mathbf{X}^\dagger \mathbf{F X}

(6)

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}

(7)

We identify

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

(8)

\mathbf{C} = \mathbf{X C}^\mathrm{OAO}

(9)

C = np.einsum("ik,kj->ij", X, C_OAO)

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

Orbital energies of occupied orbitals: [-20.5481899   -1.34520506  -0.7058451   -0.57108599  -0.49456809]

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