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:

\[\begin{split} | \Phi \rangle = | \phi_{1} \cdots \phi_{N} \rangle = \frac{1}{\sqrt{N!}} \det \begin{pmatrix} \phi_{1}(x_1) & \cdots & \phi_{N}(x_1) \\ \vdots & \ddots & \vdots \\ \phi_{1}(x_N) & \cdots & \phi_{N}(x_N) \\ \end{pmatrix}, \end{split}\]

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:

\[ E_{\mathrm{HF}} = \sum_{i=1} \langle \phi_{i} | -\frac{1}{2}\nabla^{2} + V_{\mathrm{Ne}} | \phi_{i} \rangle + \frac{1}{2} \sum_{ij} (\langle \phi_{i}\phi_{j}| \phi_{i}\phi_{j} \rangle - \langle \phi_{i}\phi_{j}| \phi_{j}\phi_{i} \rangle) = \sum_{i}h_{ii} + \frac{1}{2} \sum_{ij} \langle ij | ij \rangle - \langle ij | ji \rangle. \]

Atomic orbital basis expansion¶

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

\[ \phi^{\sigma}_{p} = \sum_{\mu}^{N_{\mathrm{AO}}} \chi_{\mu}C^{\sigma}_{\mu p} \]

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

\[ E_{\mathrm{HF}} = \sum_{\mu \nu} D_{\mu\nu} h_{\mu\nu} + \frac{1}{2} \sum_{\mu\nu\kappa\lambda} D_{\mu\nu}D_{\kappa\lambda}(\mu\nu|\kappa\lambda) - \frac{a}{2}\sum_{\mu\nu\kappa\lambda} \left[ D^{\alpha}_{\mu\kappa}D^{\alpha}_{\nu\lambda} + D^{\beta}_{\mu\kappa}D^{\beta}_{\nu\lambda} \right] (\mu\nu|\kappa\lambda) \]

with the AO density matrices defined as:

\[ D^{\sigma}_{\mu\nu} = \sum_{i=1}^{N_{\sigma}} C_{\mu i}^{\sigma} C_{\nu i}^{\sigma}, \quad D_{\mu\nu} = D^{\alpha}_{\mu\nu} + D^{\beta}_{\mu\nu}. \]

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:

\[ D^{\alpha}_{\mu\nu} = D^{\beta}_{\mu\nu} = \frac{1}{2} D_{\mu\nu} = \sum_{i=1}^{N/2} C_{\mu i} C_{\nu i} \]

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:

\[ L_{\mathrm{RHF}} = E_{\mathrm{RHF}} - \sum_{ij}\varepsilon_{ij}[\langle i | j \rangle - \delta_{ij} ], \]

and seek one of its stationary states:

\[ \frac{\partial L_{\mathrm{RHF}}}{\partial C_{\theta p}} = 0. \]

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:

\[ \mathbf{F}\mathbf{C} = \mathbf{S}\mathbf{C}\mathbf{\varepsilon}, \]

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

\[ F_{\mu\nu} = h_{\mu\nu} + \sum_{\kappa\lambda} \left[ 2(\mu\nu|\kappa\lambda) - (\mu\kappa|\nu\lambda) \right]D_{\lambda\kappa} = h_{\mu\nu} + 2J_{\mu\nu} - K_{\mu\nu}. \]

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:

Roothaan-Hall self-consistent field iterations

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:

\[ C_{ij} = \sum_{k}A_{ik}B_{kj} \]

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:

\[ \mathbf{X}^{t} \mathbf{S} \mathbf{X} = \mathbf{1} \]

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:

\[ \mathbf{S} = \mathbf{U}\mathbf{\sigma}\mathbf{U}^{t} \]

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:

\[ \mathbf{D}\mathbf{S}\mathbf{D} = \mathbf{D} \]

and normalized to the number of electrons:

\[ N = 2\,\mathrm{tr} \, \mathbf{D}\mathbf{S} = 2 \sum_{\mu\nu}D_{\mu\nu}S_{\nu\mu} \]

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^{[0]} = E_{\mathrm{NN}} + \mathrm{tr}\,\mathbf{h}\mathbf{D}^{[0]} \]

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:

\[ J^{[n]}_{\mu\nu} = \sum_{\kappa\lambda}(\mu\nu|\kappa\lambda)D_{\lambda\kappa}^{[n]} \]

and

\[ K^{[n]}_{\mu\nu} = \sum_{\kappa\lambda}(\mu\kappa|\nu\lambda)D_{\lambda\kappa}^{[n]} \]

Finally, the energy is:

\[ E^{[n]} = E_{\mathrm{NN}} + \mathrm{tr}\,(\mathbf{h} + \mathbf{F}^{[n]})\mathbf{D}^{[n]} \]

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.

The Roothaan-Hall self-consistent field procedure

Contents