# The Roothaan-Hall self-consistent field procedure¶

## Theory refresher¶

In Hartree-Fock (HF) theory, we approximate the molecular electronic wavefunction as a single Slater determinant:

with \(N = N_{\alpha} + N_{\beta}\) is the number of electrons and \(\phi_{p}\) the *molecular spin-orbitals* (MOs). From the Slater-Condon rules, the HF energy is:

### Atomic orbital basis expansion¶

The MOs will be expanded in a fixed, localized, *nonorthogonal* atomic orbital (AO) basis:

where \(\sigma\) is the spin of the MO: either \(\alpha\) or \(\beta\). The energy expression becomes:

with the AO density matrices defined as:

We will develop our HF code for *closed-shell molecules*, \(N_{\alpha} = N_{\beta} = \frac{N}{2}\), and adopt a *spin-restricted formulation*, where \(\alpha\) and \(\beta\) spin-orbitals share the same spatial part:

### Variational optimization¶

In order to obtain the best possible single-determinant representation in the given AO basis, we optimize the expansion coefficients to obtain a stationary point of the restricted Hatree-Fock energy under the constraint that the MOs remain orthonormal.
We set up a *Lagrangian*:

and seek one of its stationary states:

Through some algebra and remembering that the orbitals are *invariant* to unitary transformation *within* the occupied and virtual spaces, this is equivalent to the generalized eigenvalue problem:

with \(\mathbf{F}\) the *Fock matrix*:

The eigenvalue problem is *nonlinear*: the Fock matrix depends on its eigenvectors through the density matrix.

## The Roothaan-Hall algorithm¶

At a glance, the algorithm for solving the **Roothaan-Hall equations** is as follows:

We will start our implementation by creating a molecule and assigning it a basis set.

```
import veloxchem as vlx
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
"""
mol = vlx.Molecule.from_xyz_string(h2o_xyz)
print(mol.get_string())
basis = vlx.MolecularBasis.read(mol, "sto-3g")
print(basis.get_string(mol))
```

Next, we generate all necessary one-electron integrals: overlap, kinetic, and nuclear potential. Note that we *immediately* convert all the integrals to NumPy arrays using the `to_numpy`

function.

```
ovldrv = vlx.OverlapIntegralsDriver()
ovl = ovldrv.compute(mol, basis)
S = ovl.to_numpy()
```

```
kindrv = vlx.KineticEnergyIntegralsDriver()
kin = kindrv.compute(mol, basis)
T = kin.to_numpy()
```

```
naidrv = vlx.NuclearPotentialIntegralsDriver()
nai = naidrv.compute(mol, basis)
V = -nai.to_numpy()
```

Finally, we generate the ERI 4-index tensor. When using the `compute_in_mem`

method, VeloxChem will return a NumPy array directly.

```
eridrv = vlx.ElectronRepulsionIntegralsDriver()
eri = eridrv.compute_in_mem(mol, basis)
```

### Matrix and tensor operations with NumPy¶

We will leverage the NumPy library for all our linear algebra needs. Before moving on to the implementation of the actual algorithm, we will review the functions that we will use.

```
import numpy as np
```

For matrix multiplications and, in general, tensor contractions, it is very convenient to use np.einsum. This function accepts a string, definining the operation to perform and a variadic pack of operands. The string uses the Einstein summation convention. For example, the matrix-matrix multiplication:

is performed with:

```
C = np.einsum("ik,kj->ij", A, B)
```

NumPy includes functions for linear algebra, such as decompositions and direct linear solvers. These are in the `linalg`

submodule. For example, `np.linalg.eigh`

performs the full diagonalization of a Hermitian matrix and return its eigenvalues *and* its eigenvectors:

```
v, w = np.linalg.eigh(A)
```

We will also need to take *slices* of matrices and tensors and advanced indexing in NumPy will help with this task. For example, to select the columns between indices \(i\) and \(j\) in a matrix:

```
slice = U[:, i:j]
```

the single `:`

means to select all rows.

## Orthogonalization of the AO basis¶

The Roothaan-Hall equations are a generalized eigenvalue problem, due to the nonorthogonality of the AO basis. In order to use the `np.linalg.eigh`

function from NumPy, we first need to obtain the *orthogonal* AO (OAO) basis, that is, find a transformation matrix \(\mathbf{X}\) such that:

In the absence of linear dependencies in the AO basis, the overlap matrix is symmetric (and positive-definite) It is thus diagonalizable in a basis of its eigenvectors:

with \(\mathbf{\sigma} = \mathrm{diag}(\sigma_{0}, \cdots, \sigma_{n})\) the diagonal matrix of eigenvalues and \(\mathbf{U}\) its eigenvectors. Any function of the overlap matrix can be defined in terms of its eigenvalues and eigenvectors. This gives us two choices for the orthogonalization:

**Symmetric (Löwdin)**, with \(\mathbf{X} = \mathbf{U}\mathbf{\sigma}^{-\frac{1}{2}}\mathbf{U}^{t}\).**Canonical**, with \(\mathbf{X} = \mathbf{U}\mathbf{\sigma}^{-\frac{1}{2}}\).

In both cases, we achieve a factorization of the overlap matrix as \(\mathbf{S} = \left(\mathbf{X}^{t}\mathbf{X}\right)^{-1}\).

```
import numpy as np
sigma, U = np.linalg.eigh(S)
# TODO form the X matrix for *symmetric* orthogonalization
X =
# check that the orthogonalizer is correct
# TODO re-form S from the orthogonalizer matrix
test_S =
np.testing.assert_allclose(test_S, S, atol=1e-10)
```

## Initial guess¶

We use the very simple core guess: the initial MO coefficients are the eigenvectors of the core Hamiltonian: \(h = T + V\). The algoritm goes as follows:

Transform the core Hamiltonian \(h\) to OAO basis.

Diagonalize it.

Backtransform the coefficients to AO basis.

Form the density matrix using \(D_{\mu\nu} = \sum_{i=1}^{N/2} C_{\mu i} C_{\nu i}\).

```
N_O = mol.number_of_electrons() // 2
# TODO form the core Hamiltonian
h =
# TODO transform to orthogonal AO basis using the orthogonalizer X
h_OAO =
# TODO diagonalize
_, C_OAO =
# TODO backtransform coefficients from OAO basis to AO basis
C =
# TODO form density matrix from occupied slice of coefficients in AO basis
D =
```

We can check that the density matrix so produced is idempotent:

and normalized to the number of electrons:

These checks will ensure that our implementation produces sane density matrices!

```
# TODO check that the density matrix is normalized to the number of electrons
np.testing.assert_allclose(???, mol.number_of_electrons(), atol=1e-10)
# TODO check that the density matrix is idempotent
np.testing.assert_allclose(???, D, atol=1e-10)
```

The initial energy is:

```
E_NN = mol.nuclear_repulsion_energy()
E = np.einsum("pq,qp->", h, D) + E_NN
print(f"Initial SCF energy: {E:20.12f}")
```

## The iterative procedure¶

With a guess density at hand, we can start the Roothaan-Hall self-consistent iterations. We will run until a convergence criterion is met, `conv_thresh`

, or until we exceed the maximum number of iterations, `max_iter`

.
We choose the *norm of the difference between two successive density matrices*, \(\| \mathbf{D}^{[n+1]} - \mathbf{D}^{[n]} \|\), as convergence criterion.

```
max_iter = 50
conv_thresh = 1e-6
```

At each iteration, we form the *Coulomb*, \(J\), and *exchange*, \(K\), matrices from a contraction of the density matrix and the pre-computed ERI tensor:

and

Finally, the energy is:

```
print("# it. SCF energy (E_h) |Delta_E| Delta_D")
print(f" 0 {E:20.12f} N/A N/A")
for it in range(1, max_iter+1):
# TODO build Fock matrix in AO basis
F =
# TODO transform Fock matrix to OAO basis using the orthogonalizer X
F_OAO =
# TODO diagonalize F_OAO
_, C_OAO =
# TODO backtransform coefficients to AO basis
C =
# TODO form density matrix from occupied slice of coefficients
D_ =
# evaluate energy
E_ = np.einsum("pq,qp->", h+F, D_) + E_NN
# TODO compute convergence metric
Delta_D =
# print iteration stats
print(f" {it:>2d} {E_:20.12f} {np.abs(E_ - E):7.5e} {Delta_D:7.5e}")
# update current density and current energy
D = D_
E = E_
# TODO check if convergence threshold was met
# TODO fail if maximum number of iterations was exceeded without achieving convergence
```

## Testing our implementation¶

Let’s compare the result of our own SCF with that obtainable from VeloxChem.

```
scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(mol, basis)
ref_E = scfdrv.get_scf_energy()
```

```
print(f"Difference with reference result: {E-ref_E:5.8e}")
```

## References¶

Lehtola, S.; Blockhuys, F.; Van Alsenoy, C. An Overview of Self-Consistent Field Calculations Within Finite Basis Sets. Molecules 2020, 25 (5). https://doi.org/10.3390/molecules25051218.