Molecular Hartree-Fock

Compute the restricted-Hartree-Fock energy of the water molecule in the 6-31G* basis. This page is also the right place to learn what vibe-qc actually solves under the hood, the theory section below spells out the Roothaan-Hall equations, the Fock matrix, and Pulay’s DIIS acceleration.

Use the installed engine environment and save the script in a run directory outside the checkout. Coordinates constructed with Atom are in bohr; energies below are in hartree. See running a script if this is your first calculation.

from vibeqc import Atom, Molecule, BasisSet, run_rhf

# Water, illustrative fixed geometry in bohr (not an optimized structure).
mol = Molecule([
    Atom(8, [0.0,  0.0,  0.0]),
    Atom(1, [0.0,  1.43, -0.98]),
    Atom(1, [0.0, -1.43, -0.98]),
])
basis = BasisSet(mol, "6-31g*")
result = run_rhf(mol, basis)

print(f"E_RHF       = {result.energy:.10f}")
print(f"converged   = {result.converged}")
print(f"iterations  = {result.n_iter}")
print(f"MO energies = {result.mo_energies}")
if not result.converged:
    raise RuntimeError("RHF did not converge; inspect result.scf_trace")

Illustrative output from an earlier build follows. Treat the iteration count and rounded orbital list as orientation, not a regression fixture for every current build or set of options. Compare energies only at identical geometry, basis and convergence settings:

E_RHF       = -76.0066778869
converged   = True
iterations  = 9
MO energies = [-20.54  -1.36  -0.74  -0.56  -0.50   0.22   0.31   ...]

What just happened

  1. Building the basis. BasisSet(mol, "6-31g*") looks up the 6-31G* contraction (Ditchfield, R., Hehre, W. J., and Pople, J. A. (1971). *Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules.* Journal of Chemical Physics, **54**(2), 724–728. doi:[10.1063/1.1674902](https://doi.org/10.1063/1.1674902), Hariharan, P. C. and Pople, J. A. (1973). *The influence of polarization functions on molecular orbital hydrogenation energies.* Theoretica Chimica Acta, **28**(3), 213–222. doi:[10.1007/BF00533485](https://doi.org/10.1007/BF00533485)) through the bundled basis library. Libint (Valeev, Edward F. *Libint: A library for the evaluation of molecular integrals of many-body operators over Gaussian functions.* [Software]. <https://github.com/evaleev/libint>) places its shells on each atomic center.

  2. SCF from an AUTO seed. RHFOptions.initial_guess defaults to AUTO. For this ordinary closed-shell light molecule it selects PATOM, a SAD density followed by an in-field repolarization step. The default accelerator is EDIIS_DIIS, which switches from energy-DIIS to commutator-DIIS near convergence; it is not plain Pulay DIIS throughout. Convergence tolerances default to conv_tol_energy = 1e-8 Ha and conv_tol_grad = 1e-6. See initial guesses and SCF convergence.

  3. Output. The RHFResult carries the total energy, converged MO coefficients, orbital energies, density matrix, and the full SCF trace for post-analysis.

Tuning convergence

Most of the time the defaults are fine. For tight benchmarks you may want:

# doc-audit: skip - continues with mol and basis from the first block
from vibeqc import InitialGuess, RHFOptions

opts = RHFOptions()
opts.conv_tol_energy = 1e-12
opts.conv_tol_grad = 1e-10
opts.use_diis = True
opts.initial_guess = InitialGuess.AUTO  # retain the route-aware starting policy

result = run_rhf(mol, basis, opts)
if not result.converged:
    raise RuntimeError("The tighter RHF calculation did not converge")

Tightening the stopping tolerances does not remove basis-set or model error. If comparing guesses, change only initial_guess in a separate run and keep all other settings fixed. Convergence to a stationary solution alone does not establish the global electronic ground state.

Computing gradients

For geometry optimization or force evaluation, vibe-qc exposes analytic HF gradients:

# doc-audit: skip - uses the converged result and matching mol/basis above
from vibeqc import compute_gradient

if not result.converged:
    raise RuntimeError("Use a converged RHF result before computing its gradient")
grad = compute_gradient(mol, basis, result)   # (n_atoms, 3), Hartree / bohr

grad is an energy derivative, not a force: forces are -grad. Keep the molecule and basis identical to those used for result; changing coordinates requires a new SCF calculation before evaluating its gradient.

For this molecular RHF example the gradient is analytic; accuracy tracks SCF convergence (a few 1e-8 Ha/bohr at default tolerances: 3 to 4e-8 Ha/bohr against a tightly converged SCF for water at HF/6-31G*, vibe-qc 0.17.17 development build).

Theory

The sections below trace what the driver above actually solves: the single-determinant ansatz, the Roothaan-Hall equations it leads to, the self-consistent iteration that solves them, the DIIS acceleration that makes it converge fast, and the accuracy ceiling that comes with a single determinant.

The variational ansatz

Hartree-Fock approximates the electronic ground-state wavefunction as a single Slater determinant built from \(N\) orthonormal spin-orbitals \(\phi_i(\mathbf{r}, \sigma)\). The antisymmetric product

\[ \Psi_\text{HF}(\mathbf{r}_1 \sigma_1, \dots, \mathbf{r}_N \sigma_N) = \frac{1}{\sqrt{N!}} \det\!\bigl[\phi_i(\mathbf{r}_j, \sigma_j)\bigr] \]

encodes Pauli antisymmetry exactly. Minimizing within the chosen determinant class and finite basis gives the Hartree-Fock variational problem; an SCF iteration can still reach a local stationary solution rather than its global minimum. RHF also restricts paired electrons to a common spatial orbital. Electron correlation beyond that determinant is missing; methods such as MP2 and coupled cluster address it explicitly.

The Roothaan-Hall equations

Expand each orbital in a finite set of \(N_b\) atomic-orbital basis functions \(\{\chi_\mu\}\):

\[ \phi_i(\mathbf{r}) = \sum_\mu C_{\mu i} \, \chi_\mu(\mathbf{r}). \]

The variational minimisation of the HF energy with respect to the \(C_{\mu i}\) coefficients (subject to orthonormality \(\mathbf{C}^\dagger \mathbf{S} \mathbf{C} = \mathbf{I}\)) yields the Roothaan-Hall equations (C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951), doi:10.1103/RevModPhys.23.69; G. G. Hall, Proc. R. Soc. Lond. A 205, 541 (1951), doi:10.1098/rspa.1951.0048), a generalised eigenvalue problem in the AO basis:

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

where \(\mathbf{S}_{\mu\nu} = \braket{\chi_\mu | \chi_\nu}\) is the overlap matrix, \(\boldsymbol{\varepsilon}\) is the diagonal matrix of orbital energies, and \(\mathbf{F}\) is the Fock matrix:

\[ F_{\mu\nu} = h_{\mu\nu} + \sum_{\lambda\sigma} P_{\lambda\sigma} \Bigl[ (\mu\nu \,|\, \lambda\sigma) - \tfrac{1}{2} \, (\mu\sigma \,|\, \lambda\nu) \Bigr]. \]

Here \(h_{\mu\nu}\) is the one-electron core Hamiltonian (kinetic + nuclear-attraction), \(P_{\lambda\sigma} = 2 \sum_{i \in \text{occ}} C_{\lambda i} C_{\sigma i}\) is the one-particle density matrix in the closed-shell (RHF) case, and \((\mu\nu|\lambda\sigma)\) are the two-electron repulsion integrals in chemists’ notation. The two terms inside the brackets are the classical Coulomb \(J\) and quantum-mechanical exchange \(K\) contributions.

Self-consistency via iteration

\(\mathbf{F}\) depends on \(\mathbf{P}\), which depends on \(\mathbf{C}\), which is the eigenvector we’re trying to find, so the equations are non-linear. The standard fix is the self-consistent-field (SCF) iteration:

  1. Construct an initial density \(\mathbf{P}^{(0)}\) using the selected guess. AUTO chooses PATOM for this example, built from SAD (Van Lenthe, J. H., Zwaans, R., et al. (2006). *Starting SCF calculations by superposition of atomic densities.* Journal of Computational Chemistry, **27**(8), 926–932. doi:[10.1002/jcc.20393](https://doi.org/10.1002/jcc.20393)) plus an in-field step. Explicit unsupported guesses are refused; they are not silently replaced by HCORE.

  2. Build \(\mathbf{F}^{(n)}\) from \(\mathbf{P}^{(n)}\).

  3. Solve the generalised eigenvalue problem for \(\mathbf{C}^{(n+1)}\) and \(\boldsymbol{\varepsilon}^{(n+1)}\).

  4. Form \(\mathbf{P}^{(n+1)}\) from the occupied columns of \(\mathbf{C}^{(n+1)}\).

  5. Stop when \(|E^{(n+1)} - E^{(n)}| < \texttt{conv\_tol\_energy}\) and \(\|\mathbf{F}\mathbf{P}\mathbf{S} - \mathbf{S}\mathbf{P}\mathbf{F}\|_\text{F} < \texttt{conv\_tol\_grad}\). This is the AO stationarity residual used by the driver. The overlap factors matter in a nonorthogonal basis.

DIIS acceleration

Fixed-point iteration of the above recipe typically oscillates and converges slowly. Pulay’s Direct Inversion in the Iterative Subspace (DIIS) (Pulay, Péter (1980). *Convergence acceleration of iterative sequences. The case of SCF iteration.* Chemical Physics Letters, **73**(2), 393–398. doi:[10.1016/0009-2614(80)80396-4](https://doi.org/10.1016/0009-2614(80)80396-4), Pulay, Péter (1982). *Improved SCF convergence acceleration.* Journal of Computational Chemistry, **3**(4), 556–560. doi:[10.1002/jcc.540030413](https://doi.org/10.1002/jcc.540030413)) extrapolates the Fock matrix instead of using \(\mathbf{F}^{(n)}\) directly:

\[ \mathbf{F}_\text{DIIS} = \sum_{i} c_i \, \mathbf{F}^{(i)}, \qquad \text{minimize } \Bigl\| \sum_i c_i \, \mathbf{e}^{(i)} \Bigr\|^2 \text{ subject to } \sum_i c_i = 1, \]

where $\mathbf{e}^{(i)} = \mathbf{F}^{(i)}\mathbf{P}^{(i)}\mathbf{S}

  • \mathbf{S}\mathbf{P}^{(i)}\mathbf{F}^{(i)}\( is the AO stationarity error matrix at iteration \)i$, flattened for DIIS, the choice introduced in the 1982 paper. Near the solution DIIS converges rapidly in practice (superlinearly, not quadratically); at default settings vibe-qc takes 8-12 iterations on most closed-shell molecules, where naïve mixing would take dozens.

Where HF helps, where it hurts

HF is exact for one-electron systems (H, He⁺, Li²⁺, …) and is the starting point for the whole correlated-wavefunction hierarchy (MP2, Møller, Christian and Plesset, Milton S. (1934). *Note on an Approximation Treatment for Many-Electron Systems.* Physical Review, **46**(7), 618–622. doi:[10.1103/PhysRev.46.618](https://doi.org/10.1103/PhysRev.46.618); CCSD, Purvis, George D. III and Bartlett, Rodney J. (1982). *A full coupled-cluster singles and doubles model: The inclusion of disconnected triples.* Journal of Chemical Physics, **76**(4), 1910–1918. doi:[10.1063/1.443164](https://doi.org/10.1063/1.443164); CCSD(T), Raghavachari, Krishnan, Trucks, Gary W., et al. (1989). *A fifth-order perturbation comparison of electron correlation theories.* Chemical Physics Letters, **157**(6), 479–483. doi:[10.1016/S0009-2614(89)87395-6](https://doi.org/10.1016/S0009-2614(89)87395-6)). It has a hard ceiling on accuracy because the single determinant cannot represent the instantaneous electron-electron correlation, dynamic correlation is missing entirely. For benchmark-quality numbers you graduate to Kohn-Sham DFT (W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965), doi:10.1103/PhysRev.140.A1133; see Molecular DFT (RKS)), MP2, or coupled cluster.

Resources

~70 MB peak RAM, <0.1 s on one core (Apple M2 baseline) for the H₂O / 6-31G* example. The Fock build is \(\mathcal{O}(N_\text{bf}^4)\) in basis size, cc-pVTZ on the same H₂O is ~3 s, ~250 MB.

References

Foundation papers and a canonical textbook.

  • Roothaan-Hall equations. C. C. J. Roothaan, “New developments in molecular orbital theory,” Rev. Mod. Phys. 23, 69 (1951); G. G. Hall, “The molecular orbital theory of chemical valency. VIII,” Proc. R. Soc. Lond. A 205, 541 (1951).

  • Slater determinant. J. C. Slater, “The theory of complex spectra,” Phys. Rev. 34, 1293 (1929).

  • Self-consistent field method. D. R. Hartree, “The wave mechanics of an atom with a non-Coulomb central field,” Proc. Cambridge Philos. Soc. 24, 89 (1928); V. A. Fock, “Näherungs- methode zur Lösung des quantenmechanischen Mehrkörperproblems,” Z. Phys. 61, 126 (1930).

  • DIIS extrapolation. P. Pulay, “Convergence acceleration of iterative sequences. The case of SCF iteration,” Chem. Phys. Lett. 73, 393 (1980); P. Pulay, “Improved SCF convergence acceleration,” J. Comput. Chem. 3, 556 (1982).

  • Textbook. A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover (1996). Chapter 3 is the clearest walk through HF; chapter 5 covers MP2.

Next