Introduction

Imagine you want to predict whether a drug molecule will bind to a protein, or whether a new alloy will be stronger than steel. Both questions ultimately reduce to quantum mechanics: the electrons in every atom obey the Schrödinger equation, and the answer is written in their collective wavefunction.

The catch is brutal. A system of NN electrons needs a wavefunction Ψ(r1,r2,…,rN)\Psi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N) — a function of 3N3N continuous coordinates. For a modest fifty-electron molecule, storing that function even at coarse resolution is more data than all computers on Earth can hold. The many-body problem is not merely hard; it is exponentially hard.

In 1964, Pierre Hohenberg and Walter Kohn published a theorem that changed everything. The ground-state energy of any quantum system, they proved, is a functional of the electron density n(r)n(\mathbf{r}) — a function of just three coordinates, regardless of how many electrons the system contains. The entire quantum information is encoded in that single three-dimensional cloud. One year later, Kohn and Lu Jeu Sham showed how to use that theorem in practice through a set of single-particle equations. Walter Kohn received the Nobel Prize in Chemistry in 1998 for this work.

Density Functional Theory (DFT) is the result: a method that converts the impossible many-body Schrödinger equation into a tractable set of independent-electron equations, each living in an effective potential shaped by all the others.

Kohn-Sham Orbitals

The demo below solves the Kohn-Sham equations in one dimension for a simple external potential. Each step of the self-consistency loop updates the effective potential from the current electron density, recomputes the lowest orbitals, and rebuilds the density — repeating until nothing changes.

<!-- Kohn-Sham DFT demo: controls and canvas for orbital/density visualization -->
<div class="controls">
  <label>Potential: <select id="potSel">
    <option value="harmonic">Harmonic well</option>
    <option value="double">Double well</option>
    <option value="box">Particle in a box</option>
  </select></label>
  <label>Electrons: <input id="nElec" type="number" min="1" max="6" value="2" style="width:3rem"></label>
  <button id="runBtn" type="button">Run SCF</button>
  <button id="resetBtn" type="button" class="ghost">Reset</button>
</div>
<canvas id="canvas" width="560" height="280" aria-label="Kohn-Sham orbitals and electron density plot"></canvas>
<div id="status" class="status">Choose a potential and number of electrons, then click Run SCF.</div>
<div class="legend">
  <span class="leg-item"><span class="leg-line pot"></span> External potential</span>
  <span class="leg-item"><span class="leg-line dens"></span> Electron density n(x)</span>
  <span class="leg-item"><span class="leg-line orb"></span> KS orbitals</span>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; background: #fff; }
.controls { display: flex; flex-wrap: wrap; gap: .5rem; align-items: center; margin-bottom: .6rem; }
label { font-size: .88rem; display: flex; align-items: center; gap: .3rem; }
select, input[type=number] { font: inherit; border: 1px solid #bbb; border-radius: 6px; padding: .2rem .4rem; }
button { font: 600 14px system-ui, sans-serif; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
canvas { display: block; border: 1px solid #dde3ea; border-radius: 10px; max-width: 100%; }
.status { font-size: .92rem; font-weight: 600; min-height: 1.4em; margin: .4rem 0; color: #555; }
.status.ok { color: #0a7d33; }
.status.err { color: #c92f3c; }
.legend { display: flex; gap: 1.2rem; flex-wrap: wrap; font-size: .82rem; margin-top: .25rem; }
.leg-item { display: flex; align-items: center; gap: .35rem; }
.leg-line { display: inline-block; width: 24px; height: 3px; border-radius: 2px; }
.leg-line.pot { background: #e63946; }
.leg-line.dens { background: #1d7fd4; }
.leg-line.orb { background: #2d9c6f; }

// Kohn-Sham DFT solver: find the ground-state electron density of a 1D system
// Hohenberg-Kohn theorem: the ground-state energy is a functional of density n(x).
// Self-consistency loop: update veff from density, solve KS equations, rebuild density.

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));

const W = 560, H = 280;
const NPTS = 200;
const dx = 1.0 / (NPTS - 1); // uniform grid spacing over [0,1]

// External potential definitions (scaled for visibility on canvas)
function makePotential(name) {
  const v = new Float64Array(NPTS);
  for (let i = 0; i < NPTS; i++) {
    const x = i * dx - 0.5; // x in [-0.5, 0.5]
    if (name === 'harmonic') v[i] = 60 * x * x;
    else if (name === 'double') v[i] = 60 * (x * x - 0.0625) * (x * x - 0.0625) * 16 - 1;
    else v[i] = (Math.abs(x) > 0.45) ? 80 : 0; // infinite-wall box potential (large V outside walls)
  }
  return v;
}

// Build the tridiagonal Hamiltonian matrix using finite differences
// H = -d^2/dx^2 / 2 + veff(x)  (atomic units, hbar=m=e=1)
function buildH(veff) {
  // Finite-difference kinetic energy: T_ii = 2/dx^2, T_{i,i+1} = -1/dx^2
  const t = 1.0 / (dx * dx); // kinetic energy coefficient 1/dx^2
  const n = NPTS;
  const diag = new Float64Array(n);
  const off = new Float64Array(n - 1);
  for (let i = 0; i < n; i++) diag[i] = 2 * t + veff[i];
  for (let i = 0; i < n - 1; i++) off[i] = -t;
  return { diag, off };
}

// Find the nStates lowest eigenvectors by power iteration with deflation
// Suitable for small grids; a production solver would use LAPACK dsygvd.
function lowestEigenvectors(H, nStates) {
  const n = H.diag.length;
  const vecs = [];
  const vals = [];
  // working copy of diagonal for deflation (only diagonal is modified)
  const deflated = H.diag.slice();
  const offD = H.off.slice();

  for (let s = 0; s < nStates; s++) {
    // initial guess: sine basis function for state s+1
    let v = new Float64Array(n);
    for (let i = 0; i < n; i++) v[i] = Math.sin(Math.PI * (s + 1) * (i + 1) / (n + 1));
    normalize(v);

    // power iteration: repeatedly apply H, then orthogonalize and normalize
    for (let iter = 0; iter < 80; iter++) {
      let w = matvec(deflated, offD, v);
      // Gram-Schmidt: subtract projections onto already-converged states
      for (let k = 0; k < s; k++) {
        const dot = dotv(w, vecs[k]);
        for (let i = 0; i < n; i++) w[i] -= dot * vecs[k][i];
      }
      const lam = dotv(w, v) / dotv(v, v);
      normalize(w);
      v = w;
      if (iter > 20 && Math.abs(lam - (vals[s] || 0)) < 1e-7) break;
      vals[s] = lam;
    }
    vecs.push(v);
    vals.push(dotv(matvec(deflated, offD, v), v));
    // deflation: large positive diagonal shift removes this state from next search
    for (let i = 0; i < n; i++) deflated[i] += 200 * v[i] * v[i];
  }
  return { vecs, vals };
}

function matvec(diag, off, v) {
  const n = diag.length;
  const w = new Float64Array(n);
  for (let i = 0; i < n; i++) {
    w[i] = diag[i] * v[i];
    if (i > 0) w[i] += off[i - 1] * v[i - 1];
    if (i < n - 1) w[i] += off[i] * v[i + 1];
  }
  return w;
}
function dotv(a, b) { let s = 0; for (let i = 0; i < a.length; i++) s += a[i] * b[i]; return s; }
function normalize(v) {
  let s = 0; for (let i = 0; i < v.length; i++) s += v[i] * v[i];
  s = Math.sqrt(s * dx); for (let i = 0; i < v.length; i++) v[i] /= s;
}

// Self-consistent field (SCF) loop
function runSCF(potName, nElec) {
  const vext = makePotential(potName);
  let density = new Float64Array(NPTS);
  // uniform initial density: nElec electrons spread evenly over the grid
  for (let i = 0; i < NPTS; i++) density[i] = nElec * dx;
  let orbs = [], energies = [];
  let converged = false;
  const history = [];

  for (let scfIter = 0; scfIter < 60; scfIter++) {
    // Build effective potential: veff = vext + v_Hartree (for nElec > 1)
    // Exchange-correlation omitted for numerical stability on this coarse grid.
    const veff = new Float64Array(NPTS);
    if (nElec > 1) {
      // Hartree potential: direct electron-electron Coulomb repulsion
      // Softened 1D Coulomb: 1/|r| replaced by 1/(|r|+0.01) to avoid divergence
      for (let i = 0; i < NPTS; i++) {
        let vh = 0;
        for (let j = 0; j < NPTS; j++) {
          const r = Math.abs(i - j) * dx + 0.01; // softened Coulomb: avoids 1/0 divergence in 1D
          vh += density[j] / r * dx;
        }
        veff[i] = vext[i] + vh;
      }
    } else {
      // single electron: no Hartree term (no electron-electron repulsion)
      for (let i = 0; i < NPTS; i++) veff[i] = vext[i];
    }

    const H = buildH(veff);
    const { vecs, vals } = lowestEigenvectors(H, nElec);
    orbs = vecs;
    energies = vals;

    // new density = sum of |orbital_i(x)|^2 over all occupied states
    const newDensity = new Float64Array(NPTS);
    for (let s = 0; s < nElec; s++) {
      for (let i = 0; i < NPTS; i++) newDensity[i] += vecs[s][i] * vecs[s][i];
    }

    // linear mixing: blend new and old density to damp oscillations
    const alpha = 0.5; // alpha=0.5: equal blend of new and old density
    let diff = 0;
    for (let i = 0; i < NPTS; i++) {
      diff += Math.abs(newDensity[i] - density[i]);
      density[i] = alpha * newDensity[i] + (1 - alpha) * density[i];
    }
    history.push({ density: density.slice(), orbs: orbs.map(o => o.slice()), energies: energies.slice() });
    if (diff * dx < 1e-5 && scfIter > 2) { converged = true; break; }
  }
  return { history, converged, nIter: history.length };
}

// Canvas drawing routines
let animFrame = null;
function draw(canvas, history, step) {
  const ctx = canvas.getContext('2d');
  ctx.clearRect(0, 0, W, H);
  const PAD = { l: 30, r: 10, t: 15, b: 20 };
  const pw = W - PAD.l - PAD.r, ph = H - PAD.t - PAD.b;

  const { density, orbs } = history[step];
  const vext = currentVext;

  // compute scale factors for the three plotted curves
  const maxV = Math.max(...vext) || 1;
  const minV = Math.min(...vext);
  const maxD = Math.max(...density) || 1;

  function px(i) { return PAD.l + (i / (NPTS - 1)) * pw; }
  function pyV(val) { return PAD.t + ph - ((val - minV) / (maxV - minV + 1e-9)) * ph * 0.7; }
  function pyD(d) { return PAD.t + ph - (d / maxD) * ph * 0.85; }
  function pyOrb(o) { return PAD.t + ph * 0.5 - o * ph * 0.25; }

  // draw external potential as dashed red line
  ctx.strokeStyle = '#e63946'; ctx.lineWidth = 1.8; ctx.setLineDash([5, 3]);
  ctx.beginPath();
  for (let i = 0; i < NPTS; i++) {
    i === 0 ? ctx.moveTo(px(i), pyV(vext[i])) : ctx.lineTo(px(i), pyV(vext[i]));
  }
  ctx.stroke(); ctx.setLineDash([]);

  // draw each Kohn-Sham orbital in green shades (semi-transparent)
  const orbColors = ['#2d9c6f', '#6dbb9a', '#a8d9bb', '#c9e8d4', '#e2f3ea', '#f0f9f4'];
  for (let s = 0; s < orbs.length; s++) {
    ctx.strokeStyle = orbColors[s % orbColors.length]; ctx.lineWidth = 1.5; ctx.globalAlpha = 0.65;
    ctx.beginPath();
    for (let i = 0; i < NPTS; i++) {
      const y = pyOrb(orbs[s][i]);
      i === 0 ? ctx.moveTo(px(i), y) : ctx.lineTo(px(i), y);
    }
    ctx.stroke(); ctx.globalAlpha = 1;
  }

  // draw electron density as solid blue line
  ctx.strokeStyle = '#1d7fd4'; ctx.lineWidth = 2.4;
  ctx.beginPath();
  for (let i = 0; i < NPTS; i++) {
    i === 0 ? ctx.moveTo(px(i), pyD(density[i])) : ctx.lineTo(px(i), pyD(density[i]));
  }
  ctx.stroke();

  // draw SCF iteration counter in top-left corner
  ctx.fillStyle = '#555'; ctx.font = '12px system-ui, sans-serif';
  ctx.fillText(i18nFmt('SCF step %n% / %total%', { n: step + 1, total: history.length }), PAD.l + 4, PAD.t + 14);
}

let currentVext = makePotential('harmonic');
let currentHistory = null;

function setStatus(msg, cls) {
  const el = document.getElementById('status');
  el.textContent = msg; el.className = 'status ' + (cls || '');
}

document.getElementById('runBtn').addEventListener('click', () => {
  const potName = document.getElementById('potSel').value;
  const nElec = parseInt(document.getElementById('nElec').value, 10) || 2;
  if (nElec < 1 || nElec > 6) { setStatus('Number of electrons must be between 1 and 6.', 'err'); return; }
  currentVext = makePotential(potName);
  setStatus('Running self-consistency loop...', '');
  if (animFrame) cancelAnimationFrame(animFrame);
  // defer to allow the browser to repaint the status before the SCF loop
  setTimeout(() => {
    const { history, converged, nIter } = runSCF(potName, nElec);
    currentHistory = history;
    const canvas = document.getElementById('canvas');
    let step = 0;
    function animate() {
      draw(canvas, history, step);
      if (step < history.length - 1) { step++; animFrame = requestAnimationFrame(animate); }
      else {
        const msg = converged
          ? i18nFmt('Converged in %n% SCF iteration(s).', { n: nIter })
          : i18nFmt('Did not converge after %n% iterations — try fewer electrons.', { n: nIter });
        setStatus(msg, converged ? 'ok' : 'err');
      }
    }
    animate();
  }, 10);
});

document.getElementById('resetBtn').addEventListener('click', () => {
  if (animFrame) cancelAnimationFrame(animFrame);
  currentHistory = null;
  const canvas = document.getElementById('canvas');
  const ctx = canvas.getContext('2d');
  ctx.clearRect(0, 0, W, H);
  setStatus('Choose a potential and number of electrons, then click Run SCF.', '');
});

setStatus('Choose a potential and number of electrons, then click Run SCF.', '');
Result
Console

Watch how the electron density (solid curve) converges as the loop iterates. The density is simply the sum of the squares of the occupied orbitals: n(x)=∑i∣ϕi(x)∣2n(x) = \sum_{i} |\phi_i(x)|^2. This is the quantity that — according to the Hohenberg-Kohn theorem — contains all the information needed to determine the ground-state energy, without ever touching a many-body wavefunction.

The Real Complexity

DFT is celebrated for what it avoids, but the story has a sting in its tail.

  • Exponential to cubic. The many-body Schrödinger equation scales as O(eN)O(e^{N}) in the worst case. Kohn-Sham DFT scales as O(N3)O(N^{3}) in standard implementations (diagonalizing an N×NN \times N matrix each iteration), making thousands-of-atom simulations routine. Linear-scaling DFT variants exist for large systems.
  • The exchange-correlation problem. The Hohenberg-Kohn theorem guarantees the exact functional Exc[n]E_{xc}[n] exists, but no one knows what it is. Every practical DFT calculation uses an approximation: Local Density Approximation (LDA), Generalized Gradient Approximation (GGA), hybrid functionals mixing in exact exchange, and so on. This is the central open problem of the field — sometimes called climbing Jacob's ladder of approximations.
  • v-representability. Not every density n(r)n(\mathbf{r}) can be produced by some external potential v(r)v(\mathbf{r}) — a subtle mathematical gap between the theorem and its use. Most practical densities are fine, but the gap matters in formal developments.
  • Strongly correlated systems. When electron-electron interactions dominate — Mott insulators, high-temperature superconductors, heavy-fermion materials — standard DFT approximations fail badly. Capturing strong correlation is an active research frontier connecting DFT to quantum computing.

The result is a method that is exact in principle and approximate in practice, with accuracy that is often remarkable and occasionally wrong in important ways.

Where It Matters

DFT is the most widely used method in computational science today — the workhorse behind:

  • Drug and protein design: computing binding energies and reaction barriers for candidate molecules before synthesizing them in the lab. DFT screens millions of compounds in silico, cutting the cost of early-stage drug discovery.
  • Materials science: predicting crystal structures, band gaps, and mechanical properties of alloys, ceramics, and two-dimensional materials like graphene and MoS₂.
  • Catalysis: understanding how molecules adsorb on metal surfaces and which sites lower activation barriers — essential for designing better industrial catalysts and electrocatalysts for green hydrogen production.
  • Battery and solar cell research: calculating ion migration barriers in electrolytes, and the electronic structure of photovoltaic absorbers to optimize efficiency.
  • Semiconductor physics: modeling dopant behavior, defect levels, and interface properties in silicon, III-V semiconductors, and emerging two-dimensional materials for next-generation transistors.

The connection to computational complexity is not accidental: DFT is effective precisely because it trades an exponentially hard problem for a polynomial one — at the cost of an approximation whose quality must always be validated.

Conclusion

Density Functional Theory is one of the great intellectual shortcuts of modern science. By swapping an exponentially large wavefunction for a three-dimensional electron density, it made the quantum mechanics of real materials computable — and earned Walter Kohn a Nobel Prize.

The catch is honest: the exact exchange-correlation functional is unknown, and the approximations we use work brilliantly in many regimes and fail quietly in others. The field advances by climbing Jacob's ladder, developing better functionals, and connecting DFT to machine learning and quantum computing for the hardest correlated systems.

Every time a chemist computes a reaction barrier, or an engineer predicts a band gap, they are betting on that foundational theorem: the ground-state energy is a functional of the density, and the density alone is enough.

Share this article

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

Comments

Loading comments...

https://www.kipuhub.com/en/article/density-functional-theory/Content licensed under CC BY-NC 4.0.