Introduction

Every atom heavier than hydrogen poses an impossible problem. You cannot write down the exact quantum state of two electrons, let alone the 92 in a uranium atom, because the electrons do not move independently — each one pushes and pulls every other through the Coulomb repulsion, and that entanglement makes the wavefunction a function of all 3N3N coordinates at once. For NN electrons sampled on even a coarse grid of 10 points per coordinate, you need 103N10^{3N} numbers. Ten electrons: a trillion trillion entries.

Douglas Hartree (1928) and Vladimir Fock (1930) found a way out. Instead of tracking every electron's exact effect on every other, replace all those interactions with a single average field — each electron feels not the others' positions but their smeared-out average charge cloud. Solve for one electron in that field, use the result to update the field, and repeat until nothing changes. That is the self-consistent field (SCF) loop, and it is still the engine of computational chemistry today.

The approximation is bold: it ignores electron correlation — the tendency of electrons to avoid each other more than the average field predicts. But it captures the lion's share of the energy, scales polynomially with system size, and provides the launching pad for every higher-accuracy method in the field.

Watch the SCF Converge

The SCF loop starts with a guess for the average field, computes new orbitals, updates the field, and repeats until the energy barely changes between cycles. The demo below mimics this with a simplified model: the orbital energy ε\varepsilon at each iteration depends on the previous iteration's density.

<!-- SCF convergence demo: adjust electrons and coupling, then run or step through iterations -->
<p class="hint">Adjust the number of electrons and coupling strength, then press <b>Run SCF</b> to watch the mean-field energy converge iteration by iteration — or step through one cycle at a time.</p>
<div class="controls">
  <label>Electrons (N): <input type="range" id="nElec" min="2" max="10" value="4" step="1"> <span id="nElecVal">4</span></label>
  <label>Coupling (g): <input type="range" id="coupling" min="0.1" max="2.0" value="0.8" step="0.1"> <span id="couplingVal">0.8</span></label>
</div>
<div class="chart-wrap">
  <canvas id="chart" width="480" height="200"></canvas>
</div>
<div class="status" id="status">Ready — press Run SCF or Step.</div>
<div class="btns">
  <button id="btnRun" type="button">Run SCF</button>
  <button id="btnStep" type="button">Step</button>
  <button id="btnReset" type="button" class="ghost">Reset</button>
</div>

/* Chart and controls layout for the SCF demo */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .88rem; color: #444; margin: 0 0 .7rem; line-height: 1.5; }
.controls { display: flex; flex-direction: column; gap: .4rem; margin-bottom: .7rem; }
.controls label { font-size: .88rem; display: flex; align-items: center; gap: .5rem; }
.controls input[type=range] { flex: 1; max-width: 200px; accent-color: #1d3557; }
.controls span { font-weight: 700; font-size: .9rem; min-width: 2.5em; color: #1d3557; }
.chart-wrap { background: #f4f7fa; border: 1px solid #cdd9e3; border-radius: 10px; padding: .5rem; margin-bottom: .6rem; }
canvas { display: block; width: 100%; height: auto; }
.status { font-size: .95rem; font-weight: 600; min-height: 1.4em; margin: .3rem 0; }
.status.ok { color: #0a7d33; }
.status.bad { color: #c92f3c; }
.status.info { color: #1d3557; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: not-allowed; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* Simplified SCF model: one orbital energy level epsilon filled by N electrons. */
/* epsilon(rho) = -1 - g*rho  where rho = N (electron density, mean-field approximation). */
/* Total energy E = N*epsilon - 0.5*g*N^2*rho (double-counting correction for mean field). */

const canvas = document.getElementById('chart');
const ctx = canvas.getContext('2d');
const statusEl = document.getElementById('status');
const nElecInput = document.getElementById('nElec');
const couplingInput = document.getElementById('coupling');
const nElecVal = document.getElementById('nElecVal');
const couplingVal = document.getElementById('couplingVal');

const MAX_ITER = 40;
const TOL = 1e-5;
const DAMP = 0.5; // linear mixing (damping) factor to aid convergence for large coupling

let energyHistory = [];
let iteration = 0;
let prevEnergy = null;
let running = false;
let animId = null;

function getParams() {
  return {
    N: parseInt(nElecInput.value, 10),
    g: parseFloat(couplingInput.value)
  };
}

/* Perform one SCF energy evaluation given current density */
function scfEnergy(density, g, N) {
  /* mean-field orbital energy: kinetic (-1) plus electron-electron repulsion (g*rho) */
  const epsilon = -1.0 - g * density;
  /* total energy with double-counting correction */
  const E = N * epsilon - 0.5 * g * N * N * density;
  return { epsilon, E, newDensity: N };
}

function reset() {
  energyHistory = [];
  iteration = 0;
  prevEnergy = null;
  running = false;
  if (animId) { cancelAnimationFrame(animId); animId = null; }
  setStatus('Ready — press Run SCF or Step.', '');
  drawChart();
  document.getElementById('btnRun').disabled = false;
  document.getElementById('btnStep').disabled = false;
}

function setStatus(t, cls) {
  statusEl.textContent = t;
  statusEl.className = 'status ' + cls;
}

function step() {
  const { N, g } = getParams();
  /* start from half-filled density on first iteration */
  const density = iteration === 0 ? N * 0.5 : N;
  const damped = iteration === 0 ? density : DAMP * N + (1 - DAMP) * (prevEnergy !== null ? N : N);
  const { E } = scfEnergy(damped, g, N);
  /* apply linear mixing to smooth convergence */
  const dampedE = prevEnergy !== null
    ? DAMP * E + (1 - DAMP) * prevEnergy
    : E;
  energyHistory.push(dampedE);
  iteration++;
  const converged = prevEnergy !== null && Math.abs(dampedE - prevEnergy) < TOL;
  prevEnergy = dampedE;
  drawChart();
  if (converged || iteration >= MAX_ITER) {
    running = false;
    document.getElementById('btnRun').disabled = false;
    document.getElementById('btnStep').disabled = false;
    if (converged) {
      setStatus(i18nFmt('Converged in %iter% iteration(s). Final energy: %E% a.u.', { iter: iteration, E: dampedE.toFixed(4) }), 'ok');
    } else {
      setStatus(i18nFmt('Did not converge in %max% iterations. Try reducing coupling.', { max: MAX_ITER }), 'bad');
    }
    return false;
  }
  setStatus(i18nFmt('Iteration %iter% — energy: %E% a.u.', { iter: iteration, E: dampedE.toFixed(4) }), 'info');
  return true;
}

function runAll() {
  running = true;
  document.getElementById('btnRun').disabled = true;
  document.getElementById('btnStep').disabled = true;
  function tick() {
    if (!running) return;
    const cont = step();
    if (cont) animId = requestAnimationFrame(tick);
  }
  animId = requestAnimationFrame(tick);
}

/* Draw the energy vs. iteration chart */
function drawChart() {
  const W = canvas.width, H = canvas.height;
  ctx.clearRect(0, 0, W, H);
  /* draw axes */
  ctx.strokeStyle = '#cdd9e3';
  ctx.lineWidth = 1;
  ctx.beginPath(); ctx.moveTo(40, 10); ctx.lineTo(40, H - 30); ctx.lineTo(W - 10, H - 30); ctx.stroke();
  /* axis labels */
  ctx.fillStyle = '#888';
  ctx.font = '11px system-ui';
  ctx.fillText('Iteration', W / 2 - 20, H - 8);
  ctx.save();
  ctx.translate(12, H / 2 + 20);
  ctx.rotate(-Math.PI / 2);
  ctx.fillText('Energy (a.u.)', 0, 0);
  ctx.restore();
  if (energyHistory.length < 2) return;
  const minE = Math.min(...energyHistory);
  const maxE = Math.max(...energyHistory);
  const rangeE = maxE - minE || 1;
  const toX = (i) => 40 + (i / (MAX_ITER - 1)) * (W - 55);
  const toY = (e) => H - 30 - ((e - minE) / rangeE) * (H - 45);
  /* draw energy curve */
  ctx.strokeStyle = '#1d3557';
  ctx.lineWidth = 2;
  ctx.lineJoin = 'round';
  ctx.beginPath();
  energyHistory.forEach((e, i) => {
    if (i === 0) ctx.moveTo(toX(i), toY(e));
    else ctx.lineTo(toX(i), toY(e));
  });
  ctx.stroke();
  /* highlight current point in red */
  energyHistory.forEach((e, i) => {
    ctx.fillStyle = i === energyHistory.length - 1 ? '#e63946' : '#1d3557';
    ctx.beginPath();
    ctx.arc(toX(i), toY(e), i === energyHistory.length - 1 ? 5 : 3, 0, 2 * Math.PI);
    ctx.fill();
  });
}

nElecInput.addEventListener('input', () => { nElecVal.textContent = nElecInput.value; reset(); });
couplingInput.addEventListener('input', () => { couplingVal.textContent = parseFloat(couplingInput.value).toFixed(1); reset(); });
document.getElementById('btnRun').addEventListener('click', runAll);
document.getElementById('btnStep').addEventListener('click', () => { if (!running) step(); });
document.getElementById('btnReset').addEventListener('click', reset);
drawChart();
Result
Console

Notice that a weak coupling (small electron count or interaction strength) converges in just a few steps. A strong coupling takes many more — and in real quantum-chemistry codes, convergence failure is a genuine headache that requires damping, level shifting, or more sophisticated algorithms like DIIS.

The Real Complexity

What does it cost to run Hartree-Fock — and why does it matter?

  • Exact solution (Full CI): the many-electron Schrödinger equation in a basis of KK functions has a Hilbert space of (KN)\binom{K}{N} Slater determinants. For moderate systems this is billions to trillions of states — exponential in NN.
  • Hartree-Fock: each SCF cycle requires computing the Fock matrix, whose most expensive part is the two-electron repulsion integrals. Naively there are O(K4)O(K^4) of them, giving a quartic scaling with basis size. Modern algorithms with screening can reach O(K2.5)O(K^{2.5}) or better for large sparse systems.
  • The correlation gap: Hartree-Fock misses roughly 1% of the total energy — the correlation energy. That sounds small, but chemical bonds involve energy differences of 0.1–1%, so the error matters. Post-HF methods (MP2, CCSD, CCSD(T)) recover correlation at costs of O(K5)O(K^5), O(K6)O(K^6), and O(K7)O(K^7) respectively.
  • Koopmans' theorem (1934): a striking exactness within the approximation — the negative of each occupied orbital energy equals the ionization potential for removing that electron, to first order. An approximation that forgets correlation yet predicts ionization energies well.

The Hartree-Fock method was proved exact in the limit of infinite basis size for its own variational problem by Hartree and Fock themselves. Its limitation is not numerical but conceptual: the ansatz forces the wavefunction into a single Slater determinant, which cannot represent entangled many-electron states. Related ideas appear in quantum simulation, where representing even a modest entangled state on a classical computer costs exponential memory.

Where It Matters

Hartree-Fock is the foundation on which nearly all of computational chemistry is built:

  • Drug discovery: HF and its post-HF descendants compute binding energies, charge distributions and reaction paths for candidate molecules. Every major pharmaceutical company runs these calculations.
  • Materials design: predicting band gaps, magnetic moments and phase stability of solids starts with a mean-field ground state — often density functional theory (DFT), which inherits the SCF loop directly from HF.
  • Post-HF hierarchy: Møller-Plesset perturbation theory (MP2), coupled-cluster (CCSD, CCSD(T) — often called the "gold standard"), and multireference methods all use the Hartree-Fock reference state as their starting point.
  • Density functional theory (DFT): Kohn and Sham (1965) recast the many-electron problem as an SCF problem over single-particle orbitals, keeping the mean-field structure of HF while replacing the exchange with an exchange-correlation functional. DFT now dominates materials science and large-molecule chemistry.
  • Teaching quantum complexity: HF is a vivid illustration of how an exponentially hard problem (quantum simulation) can be brought to polynomial cost by a well-chosen approximation — with a known, bounded error.

From aspirin to semiconductor design, the self-consistent field loop that Hartree and Fock devised nearly a century ago remains the workhorse of molecular science.

Conclusion

Hartree-Fock theory offers a lesson that echoes far beyond chemistry: when the exact problem is exponentially hard, a well-designed approximation that captures most of the physics can make the problem polynomial. The cost is the correlation energy — roughly 1% of the total — but the gain is a method that runs on a laptop for molecules with hundreds of atoms.

The SCF loop is a fixed-point iteration: start with any reasonable guess, update the field, repeat. If it converges, you have found a self-consistent solution where each electron's orbital is optimal given the average field of all the others. If it diverges, you need better numerics — a reminder that tractability and convergence are not the same thing.

Every time a chemist runs a DFT calculation or a physicist simulates a crystal, they are running a descendant of the SCF idea that Hartree and Fock worked out nearly a century ago. The exponential wall of quantum simulation has not moved — but the mean field gives us a polynomial ladder that reaches most of what we need.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/hartree-fock/Content licensed under CC BY-NC 4.0.