Hartree-Fock Approximation

Since you should already be familiar with the Hartree-Fock approximation from the bachelor course on theoretical chemistry, we will only briefly summarize the algorithm here. A derivation of the Hartree-Fock equations has been or will be covered in the master courses on quantum chemistry.

The SCF procedure in the Hartree-Fock approximation is as follows:

1. Specify a molecule and a basis set.
2. Calculate all required molecular integrals (overlap integrals \(\mathbf{S}\), kinetic energy integrals \(\mathbf{T}\), nuclear attraction integrals \(\mathbf{V}_\mathrm{ne}\), electron repulsion integrals \((\mu\nu|\sigma\lambda\)).
3. Calculate the transformation matrix \(\mathbf{X} = \mathbf{S}^{-1/2}\).
4. Obtain a guess for the density matrix \(\mathbf{P}\).
5. Calculate the matrix \(\mathbf{G}\) from the density matrix using \(G_{\mu\nu} = \sum_{\lambda\sigma} P_{\lambda\sigma} [ (\mu\nu|\sigma\lambda) - \frac{1}{2} (\mu\lambda|\sigma\nu)\)].
6. Obtain the Fock matrix \(\mathbf{F} = \mathbf{T} + \mathbf{V}_\mathrm{ne} + \mathbf{G}\).
7. Calculate the transformed Fock matrix \(\mathbf{F}' = \mathbf{X}^\dagger \mathbf{F} \mathbf{X}\)
8. Diagonalize the transformed Fock matrix \(\mathbf{F}'\) to obtain the orbital energies \(\varepsilon_i\) and the transformed MO coefficients \(\mathbf{C}'\).
9. Transform the MO coefficients back to the original basis set \(\mathbf{C} = \mathbf{X} \mathbf{C}'\).
10. Form a new density matrix \(\mathbf{P}\) from \(\mathbf{C}\) using \(P_{\mu\nu} = 2 \sum_{i=1}^{n_{\mathrm{occ}}} C_{\mu i} C_{\nu i}\).
11. Check for convergence. If not, return to step 5 with the new density matrix.
12. If the procedure has converged, calculate the quantities of interest (e.g. total energy, dipole moment, etc.) from the converged solution.

We can implement the Hartree-Fock procedure in a Python class. Since we want to access lots of properties of the molecule, we shall let the HartreeFock class inherit from the Molecule class.

class HartreeFock(Molecule):

    def initialize(self):
        self.nel = np.array([
            self.atomlist[i].atnum for i in range(0, len(self.atomlist))
        ]).sum()
        self.nocc = self.nel // 2
        
        # precalculate integrals 
        self.get_S()
        self.get_T()
        self.get_V()
        self.get_twoel()
        
        # orthogonalize AO
        eigval, eigvec = np.linalg.eigh(self.S)
        self.X = eigvec @ np.diag(1.0 / np.sqrt(eigval))

        # core hamiltonian + initialize density matrix
        self.hcore = self.T + self.Ven
        orb_en, orb = np.linalg.eigh(self.hcore)
        
        c = orb[:, :self.nocc]
        self.p = c @ c.T
    
    def get_fock(self, p):
        g = np.einsum(
            'kl, ijkl -> ij', p, 
            2.0 * self.twoel - self.twoel.transpose(0, 2, 1, 3),
        )
        return self.hcore + g
    
    def run_hf(self, max_iter=100, threshold=1e-6, verbose=5):
        energy_last_iteration = 0.0
        self.p_old = np.copy(self.p)
        for iteration in range(max_iter):
            # calculate Fock-matrix
            f = self.get_fock(self.p)
            # orthogonalize Fock-Matrix
            f_ortho = self.X.T @ f @ self.X
            # diagonalize Fock-Matrix
            eigvals, eigvect = np.linalg.eigh(f_ortho)
            # get new density matrix 
            c = eigvect[:, :self.nocc]
            self.p = c @ c.T
            self.p = self.X @ self.p @ self.X.T
            # calculate energy 
            energy = np.trace((self.hcore + f) @ self.p)
            
            if verbose > 0:
                print(f"Iteration {iteration}, Energy = {energy} Hartree")
            if verbose > 1:
                print(f"MO energies: {eigvals}")
            
            if np.abs(energy - energy_last_iteration) < threshold:
                break
            energy_last_iteration = energy
        
        self.mo_energy = eigvals
        self.mo_coeff = self.X @ eigvect
        self.energy = energy

        return energy

Now we test our HF implementation on the water molecule using the STO-3G basis set.

    from atom import Atom
    
    # Coordinates are in the unit of Angstrom.
    o1 = Atom('O', [ 0.000,  0.000,  0.000], unit='A')
    h1 = Atom('H', [ 0.758,  0.587,  0.000], unit='A')
    h2 = Atom('H', [-0.758,  0.587,  0.000], unit='A')
    
    water = HartreeFock()
    water.set_atomlist([o1, h1, h2])
    water.get_basis('sto-3g')
    
    water.initialize()
    e_scf = water.run_hf()
    print(f"SCF energy: {e_scf} Hartree")

We get the final (electronic) energy of \(E_\mathrm{SCF} = -84.143659\ \mathrm{a.u.}\) To obtain the total HF energy, the nuclear repulsion energy \(E_\mathrm{ne}\) has to be added to it.