The purpose of this project is to provide a deeper understanding of Hartree-Fock theory by demonstrating a simple implementation of the self-consistent-field method. The theoretical background can be found in Ch. 3 of the text by Szabo and Ostlund or in the nice set of on-line notes written by David Sherrill. A concise set of instructions for this project may be found here.
We thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project.
The test case used in the following discussion is for a water molecule with a bond-length of 1.1 Å and a bond angle of 104.0o with an STO-3G basis set. The input to the project consists of the nuclear repulsion energy and pre-computed sets of one- and two-electron integrals: overlap integrals, kinetic-energy integrals, nuclear-attraction integrals, electron-electron repulsion integrals.
Read the nuclear repulsion energy from the enuc.dat file.
Read the AO-basis overlap,
and store them in appropriately constructed matrices. Then form the “core Hamiltonian”:
Note that the one-electron integrals provided include only the permutationally unique integrals, but you should store the full matrices for convenience. Note also that the AO indices on the integrals in the files start with “1” rather than “0”.
Read the two-electron repulsion integrals from the eri.dat file. The integrals in this file are provided in Mulliken notation over real AO basis functions:
Hence, the integrals obey the eight-fold permutational symmetry relationships:
and only the permutationally unique integrals are provided in the file, with the restriction that, for each integral, the following relationships hold:
Note that the two-electron integrals may be stored efficiently in a one-dimensional array and the above relationship used to map between given , , , and indices and a “compound index” defined as:
Diagonalize the overlap matrix:
where is the matrix of eigenvectors (columns) and is the diagonal matrix of corresponding eigenvalues.
Build the symmetric orthogonalization matrix using:
where the tilde denotes the matrix transpose.
Form an initial (guess) Fock matrix in the orthonormal AO basis using the core Hamiltonian as a guess:
Diagonalize the Fock matrix:
Note that the matrix contains the initial orbital energies.
Transform the eigenvectors into the original (non-orthogonal) AO basis:
Build the density matrix using the occupied MOs:
where m indexes the columns of the coefficient matrices, and the summation includes only the occupied spatial MOs.
The SCF electronic energy may be computed using the density matrix as:
The total energy is the sum of the electronic energy and the nuclear repulsion energy:
where 0 denotes the initial SCF energy.
Start the SCF iterative procedure by building a new Fock matrix using the previous iteration's density as:
where the double-summation runs over all the AOs and i-1 denotes the density for the last iteration.
Form the new density matrix following the same procedure as in Step #5 above:
Compute the density:
where i denotes the current iteration density.
Compute the new SCF energy as before:
where i denotes the SCF energy for the ith iteration.
Test both the energy and the density for convergence:
If the difference in consecutive SCF energy and the root-mean-squared difference in consecutive densities do not fall below the prescribed thresholds, return to Step #7 and continue from there.
At convergence, the canonical Hartree-Fock MOs are, by definition, eigenfunctions of the Fock operator, viz.
If we multiply on the left by an arbitrary MO and integrate, we obtain:
In other words, the Fock matrix should be diagonal in the MO basis, with the orbital energies as its diagonal elements. We can demonstrate this explicitly using the AO-basis Fock matrix by first re-writing the above expression using the LCAO-MO coefficients:
Use the above equation to transform the Fock matrix from the AO basis to the MO basis and demonstrate that it is indeed diagonal (to within the convergence limits of the SCF iterative procedure).
As discussed in detail in Ch. 3 of the text by Szabo and Ostlund, the calculation of one-electron properties requires density matrix and the relevant property integrals. The electronic contribution to the electric-dipole moment may be computed using,
where the vector notation implies three sets of dipole-moment integrals – one for each Cartesian component of the dipole operator.
Two points to note:
The test cases provided below include the structural information dipole integrals needed to compute the dipole moment.
A Mulliken population analysis (also described in Szabo & Ostlund, Ch. 3) requires the overlap integrals and the electron density, in addition to information about the number of basis functions centered on each atom. The charge on atom A may be computed as:
where the summation is limited to only those basis functions centered on atom A.