feat(chemistry): Boys-Bernardi counterpoise / BSSE correction

Eighth Tier 3 item. Closes the "Counterpoise correction (BSSE)"
row in LIMITATIONS.md — unlocks noncovalent interaction energies
(H-bond, π-stacking, vdW) at HF.

Ghost-atom support (atoms.ts):
- New `ghost?: boolean` flag on `Atom`. When true, the atom
  contributes its basis functions to the AO space but **no nuclear
  charge** (Z=0) and **no electrons**. `moleculeToShellsNuclei`
  honors the flag: shells emitted as normal, Z=0 in Nucleus entry
  so V_ne and V_nn contributions vanish, electron count unchanged.

Counterpoise module (counterpoise.ts, ~150 lines):
- `runCounterpoise(atoms, fragments, basis, opts)` does 1 + 2·N
  HF runs (N fragments): full supermolecule + each fragment in
  the dimer basis (others ghosted) + each fragment in bare basis.
- Returns: supermolecule energy, per-fragment CP + bare energies,
  uncorrected ΔE, CP-corrected ΔE_CP, and BSSE = ΔE_CP − ΔE
  (always ≥ 0 by variational principle).
- Validates fragmentation: every atom in exactly one fragment.

Tests (5 green) — `counterpoise.test.ts`:
- H₂...H₂ STO-3G at 3 Å: BSSE positive, sub-mHa.
- Supermolecule energy from counterpoise matches plain HF.
- Single-fragment case (no ghosts): BSSE = 0, CP == bare exactly.
- Variational: each fragment's ghost-augmented energy ≤ bare.
- API validation: rejects overlapping / missing / out-of-range
  fragment assignments.

Limitations (queued):
- HF only at this entry. MP2 / CCSD counterpoise = mechanical
  follow-up (pass ghost atoms through post-HF correlation step).
- Closed-shell only (assumes RHF). UHF counterpoise similar.

Full suite: 434 passing / 1 skipped (+5 from this commit).

Tier 3 progress: 8 items shipped this session.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: abgnydn

LIMITATIONS.md

@@ -109,7 +109,7 @@ will silently truncate large dispatches.
 
 | missing | impact | roadmap |
 |---|---|---|
-| Counterpoise correction (BSSE) | Can't quote noncovalent interaction energies | Tier 3 — additive on top of HF/CCSD |
+| Counterpoise / BSSE (HF) | **Shipped 2026-05**. `runCounterpoise(atoms, fragments, basis, opts)` in `src/chemistry/counterpoise.ts`. Boys-Bernardi: per-fragment HF in the full dimer basis via ghost atoms (new `ghost?: boolean` flag on `Atom` — basis present, Z=0, no electrons). Returns supermolecule energy, per-fragment energies in dimer + bare basis, uncorrected + CP interaction energies, and `bsseCorrection = ΔE_CP − ΔE ≥ 0`. Verified by `counterpoise.test.ts` on H₂...H₂ STO-3G at 3 Å. MP2 / CCSD counterpoise is a mechanical follow-up (pass ghost atoms through post-HF). Open-shell counterpoise similar. | shipped ✓ (HF only) |
 | Frozen-core in CCSD(T) | **Audited 2026-05** — CCSD(T) frozen-1s on H₂O / CH₄ lies above all-electron by < 30 mHa, (T) stays negative. Verified by `tests/chemistry/frozen-core-audit.test.ts`. | shipped ✓ |
 | Frozen-core in UCCSD | **Audited + fixed 2026-05**. Audit found UCCSD's frozen-core was freezing α-occupied SOs [0, nFrozenSO) instead of α + β of the lowest k spatials (SO-ordering mismatch with ccsdIterate's contiguous-frozen contract). Fix: ccsdIterate now accepts `ReadonlySet<number>` of frozen SO indices; UCCSD constructs the correct interleaved (α-spatial-s, β-spatial-s) set. Verified by `frozen-core-audit.test.ts` — UCCSD frozen-core on closed-shell H₂O now matches RHF-CCSD frozen-core to 1e-6 Ha. | shipped ✓ |
 | Frozen-core in EE-EOM-CCSD | **Audited + shipped 2026-05**. `runEOMCCSD` accepts `nFrozenCore`. Implementation: packed (singles + antisym doubles) basis restricted to occupied indices ≥ 2·nFrozenCore; σ-equation internals unchanged (R_1 / R_2 are zero at frozen indices, so frozen-index contributions to the inner sums vanish automatically). Verified by `frozen-core-audit.test.ts` — H₂O STO-3G frozen-1s lowest 3 EOM excitations are real-positive, ordered, and shift by < 100 mHa from all-electron. Dense + Davidson paths both honor frozen-core. | shipped ✓ |

src/chemistry/atoms.ts

@@ -42,6 +42,13 @@ export interface Atom {
   readonly symbol: AtomSymbol;
   /** Position in Ångströms (will be converted to Bohr). */
   readonly pos: readonly [number, number, number];
+  /** If true, this atom contributes its basis functions to the calculation
+   *  but **no nuclear charge** (Z=0) and **no electrons**. Used by
+   *  Boys-Bernardi counterpoise correction to compute a fragment's energy
+   *  in the full dimer basis — the "ghost" augments the AO space without
+   *  changing the molecular Hamiltonian's nuclear / electron count.
+   *  Default false. */
+  readonly ghost?: boolean;
 }
 
 /** Atomic number for each supported atom. */
@@ -258,8 +265,14 @@ export function moleculeToShellsNuclei(
     const newShells = atomShells(a.symbol, pos_bohr, basis);
     shells.push(...newShells);
     for (let s = 0; s < newShells.length; s++) shellAtomIdx.push(ai);
-    nuclei.push({ Z: Z_FOR[a.symbol], pos: pos_bohr });
-    nElectrons += N_ELECTRONS_FOR[a.symbol];
+    // Ghost atoms contribute their AO basis but no nuclear charge and
+    // no electrons (Boys-Bernardi counterpoise). Push a Z=0 nucleus
+    // entry so the per-atom index stays aligned with shellAtomIdx but
+    // V_ne contributions from this site vanish, and V_nn pair terms
+    // involving a ghost are zero.
+    const isGhost = a.ghost === true;
+    nuclei.push({ Z: isGhost ? 0 : Z_FOR[a.symbol], pos: pos_bohr });
+    if (!isGhost) nElectrons += N_ELECTRONS_FOR[a.symbol];
   }
   return { shells, nuclei, nElectrons, shellAtomIdx };
 }

src/chemistry/counterpoise.ts

@@ -0,0 +1,171 @@
+// ─────────────────────────────────────────────────────────────
+// counterpoise.ts — Boys-Bernardi counterpoise correction for
+// the basis-set superposition error (BSSE) in supermolecular
+// interaction energies.
+//
+// Why BSSE exists:
+//   When you compute ΔE = E(AB) − E(A) − E(B) for a dimer A···B,
+//   in the dimer calculation A's electrons can "borrow" B's basis
+//   functions (and vice-versa), giving A a better effective basis
+//   than it has alone. The monomer A calculation uses only A's
+//   own basis, so its energy is described in a worse basis. Result:
+//   the dimer looks artificially more stable than it should.
+//
+// Boys-Bernardi fix (CP):
+//   ΔE_CP = E(AB) − E(A in dimer basis) − E(B in dimer basis)
+//
+//   where "X in dimer basis" means: run HF on X's real atoms,
+//   keeping every OTHER atom's basis functions present at their
+//   normal positions but with NO nuclear charge (Z=0) and NO
+//   electrons. The ghost atom provides AO basis only.
+//
+// API:
+//   `runCounterpoise(atoms, fragments, basis, opts)`
+//   - `atoms`: the full supermolecule.
+//   - `fragments`: array of arrays-of-atom-indices. Each fragment
+//     is one monomer. The union must cover every atom index
+//     exactly once.
+//   - Returns `CounterpoiseResult` with:
+//       supermolecule energy,
+//       per-fragment energy in the full dimer basis (CP),
+//       per-fragment energy in the bare monomer basis (no CP),
+//       uncorrected interaction energy ΔE,
+//       counterpoise-corrected ΔE_CP,
+//       BSSE = ΔE − ΔE_CP   (positive when the dimer was
+//       artificially stabilized by basis-set borrowing).
+//
+// Limitations:
+//   - HF only at this entry point. MP2 / CCSD counterpoise is a
+//     mechanical follow-up (same pattern — pass ghost atoms
+//     through to the post-HF correlation step).
+//   - Assumes every fragment is closed-shell. Open-shell (UHF)
+//     counterpoise is also mechanical to add.
+//   - Per-fragment fragmentation must be a strict partition; we
+//     don't validate physical sensibleness (you can put atoms
+//     in arbitrary groups, but the BSSE interpretation only
+//     makes sense for non-bonded fragments).
+// ─────────────────────────────────────────────────────────────
+
+import type { Atom, BasisName } from "./atoms.js";
+import { moleculeToShellsNuclei } from "./atoms.js";
+import { computeMolecularIntegrals } from "./cg-molecular.js";
+import { runRHFSCF, type HFOpts } from "./hf-scf.js";
+
+export interface CounterpoiseFragment {
+  /** Indices into the `atoms` array that belong to this fragment. */
+  readonly atomIndices: readonly number[];
+}
+
+export interface CounterpoiseResult {
+  /** E(AB) — full supermolecule with all real atoms, Hartree. */
+  readonly supermoleculeEnergy: number;
+  /** Per-fragment energies in the full dimer basis (ghost-augmented). */
+  readonly fragmentEnergiesCP: readonly number[];
+  /** Per-fragment energies in the bare monomer basis (no ghosts). */
+  readonly fragmentEnergiesBare: readonly number[];
+  /** Uncorrected ΔE = E(AB) − Σ E(fragment bare). */
+  readonly interactionEnergy: number;
+  /** Counterpoise-corrected ΔE_CP = E(AB) − Σ E(fragment in dimer basis). */
+  readonly interactionEnergyCP: number;
+  /** BSSE = ΔE_CP − ΔE = Σ_k (E_bare(k) − E_CP(k)). Always ≥ 0 by
+   *  the variational principle: ghost-augmented monomer energies
+   *  are always ≤ bare-basis energies, so CP raises the interaction
+   *  energy (makes the dimer look less attractive / more repulsive).
+   *  In small basis sets (STO-3G, cc-pVDZ) this can be 10-50% of the
+   *  true interaction energy for hydrogen-bonded systems. */
+  readonly bsseCorrection: number;
+  /** True iff every HF run (supermolecule + 2·N fragment runs)
+   *  converged. If false, the energies are returned but should
+   *  not be trusted. */
+  readonly allConverged: boolean;
+}
+
+/**
+ * Run the Boys-Bernardi counterpoise correction on a closed-shell
+ * supermolecule partitioned into fragments. Performs (1 + 2·N) HF
+ * runs where N is the number of fragments:
+ *   1× E(AB) on the full supermolecule
+ *   N× E(fragment in dimer basis) — fragment real, others ghost
+ *   N× E(fragment bare) — fragment real, others absent
+ */
+export function runCounterpoise(
+  atoms: readonly Atom[],
+  fragments: readonly CounterpoiseFragment[],
+  basis: BasisName = "sto-3g",
+  hfOpts: HFOpts = {},
+): CounterpoiseResult {
+  // ── Validate fragmentation: every atom in exactly one fragment. ──
+  const claimed = new Array<boolean>(atoms.length).fill(false);
+  for (const f of fragments) {
+    for (const idx of f.atomIndices) {
+      if (idx < 0 || idx >= atoms.length) {
+        throw new Error(`runCounterpoise: fragment atom index ${idx} out of range [0, ${atoms.length})`);
+      }
+      if (claimed[idx]) {
+        throw new Error(`runCounterpoise: atom ${idx} appears in multiple fragments`);
+      }
+      claimed[idx] = true;
+    }
+  }
+  for (let i = 0; i < atoms.length; i++) {
+    if (!claimed[i]) {
+      throw new Error(`runCounterpoise: atom ${i} is not assigned to any fragment`);
+    }
+  }
+
+  let allConverged = true;
+
+  // ── Supermolecule HF. ──
+  const superSh = moleculeToShellsNuclei(atoms, basis);
+  const superInt = computeMolecularIntegrals(superSh.shells, superSh.nuclei);
+  const superHF = runRHFSCF(superInt, superSh.nElectrons, hfOpts);
+  if (!superHF.converged) allConverged = false;
+
+  // ── Per-fragment runs. ──
+  const fragmentEnergiesCP: number[] = [];
+  const fragmentEnergiesBare: number[] = [];
+
+  for (const f of fragments) {
+    const fragmentSet = new Set(f.atomIndices);
+    // a) Fragment in dimer basis: all atoms present, but only this
+    //    fragment's atoms keep their nuclei + electrons. Others are
+    //    ghosts (Z=0, no electrons, basis functions present).
+    const dimerBasisAtoms: Atom[] = atoms.map((a, i) =>
+      fragmentSet.has(i) ? a : { ...a, ghost: true },
+    );
+    const cpSh = moleculeToShellsNuclei(dimerBasisAtoms, basis);
+    const cpInt = computeMolecularIntegrals(cpSh.shells, cpSh.nuclei);
+    const cpHF = runRHFSCF(cpInt, cpSh.nElectrons, hfOpts);
+    if (!cpHF.converged) allConverged = false;
+    fragmentEnergiesCP.push(cpHF.energy);
+
+    // b) Fragment in bare monomer basis: only this fragment's atoms exist.
+    const bareAtoms = atoms.filter((_, i) => fragmentSet.has(i));
+    const bareSh = moleculeToShellsNuclei(bareAtoms, basis);
+    const bareInt = computeMolecularIntegrals(bareSh.shells, bareSh.nuclei);
+    const bareHF = runRHFSCF(bareInt, bareSh.nElectrons, hfOpts);
+    if (!bareHF.converged) allConverged = false;
+    fragmentEnergiesBare.push(bareHF.energy);
+  }
+
+  let sumCP = 0;
+  let sumBare = 0;
+  for (let k = 0; k < fragments.length; k++) {
+    sumCP += fragmentEnergiesCP[k]!;
+    sumBare += fragmentEnergiesBare[k]!;
+  }
+  const interactionEnergy = superHF.energy - sumBare;
+  const interactionEnergyCP = superHF.energy - sumCP;
+  // BSSE = ΔE_CP − ΔE = Σ (E_bare − E_CP). Variationally ≥ 0.
+  const bsseCorrection = interactionEnergyCP - interactionEnergy;
+
+  return {
+    supermoleculeEnergy: superHF.energy,
+    fragmentEnergiesCP,
+    fragmentEnergiesBare,
+    interactionEnergy,
+    interactionEnergyCP,
+    bsseCorrection,
+    allConverged,
+  };
+}

tests/chemistry/counterpoise.test.ts

@@ -0,0 +1,148 @@
+// Counterpoise / BSSE — Tier 3 (Boys-Bernardi).
+//
+// Pass bars:
+//   1. Sign: BSSE correction is POSITIVE for non-bonded fragments
+//      in any incomplete basis (dimer's borrowed basis artificially
+//      stabilizes; CP removes the stabilization, making ΔE_CP > ΔE).
+//   2. Ghost-only sanity: a fragment computed with its real atoms +
+//      no ghosts must match the same fragment computed normally
+//      (counterpoise reduces to plain HF when there's nothing
+//      ghosted).
+//   3. Self-consistency: the supermolecule energy returned by
+//      counterpoise must equal the plain HF energy on the same atoms.
+//   4. Magnitude bounded: STO-3G BSSE on H₂...H₂ at 3 Å should be
+//      sub-mHa (small basis on a non-interacting pair).
+
+import { describe, expect, test } from "vitest";
+import { runCounterpoise } from "../../src/chemistry/counterpoise.js";
+import { computeMolecularIntegrals } from "../../src/chemistry/cg-molecular.js";
+import { moleculeToShellsNuclei, type Atom } from "../../src/chemistry/atoms.js";
+import { runRHFSCF } from "../../src/chemistry/hf-scf.js";
+
+describe("Counterpoise / BSSE", () => {
+  test("H₂...H₂ STO-3G at 3 Å: BSSE positive, sub-mHa, all converged", () => {
+    // Two H₂ molecules end-to-end, 3 Å between the closest H atoms.
+    // STO-3G is a poor basis → BSSE is the *only* "binding" you see
+    // for two H₂ molecules far apart.
+    const atoms: Atom[] = [
+      { symbol: "H", pos: [0, 0, 0]      },
+      { symbol: "H", pos: [0, 0, 0.7414] },
+      { symbol: "H", pos: [0, 0, 3.7414] },  // 3 Å gap
+      { symbol: "H", pos: [0, 0, 4.4828] },
+    ];
+    const cp = runCounterpoise(
+      atoms,
+      [
+        { atomIndices: [0, 1] },
+        { atomIndices: [2, 3] },
+      ],
+      "sto-3g",
+      { useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8 },
+    );
+
+    expect(cp.allConverged).toBe(true);
+    // BSSE = ΔE_CP − ΔE = Σ(E_bare − E_CP) ≥ 0 by the variational
+    // principle (ghost-augmented monomer is always ≤ bare).
+    expect(cp.bsseCorrection).toBeGreaterThan(0);
+    // Magnitude sub-mHa on STO-3G H₂...H₂ at 3 Å.
+    expect(cp.bsseCorrection).toBeLessThan(1e-3);
+    // CP-corrected ΔE > uncorrected ΔE always (CP raises the
+    // interaction energy by exactly bsseCorrection).
+    expect(cp.interactionEnergyCP).toBeGreaterThan(cp.interactionEnergy);
+    expect(cp.interactionEnergyCP - cp.interactionEnergy).toBeCloseTo(cp.bsseCorrection, 12);
+  });
+
+  test("supermolecule energy returned by counterpoise matches plain HF", () => {
+    const atoms: Atom[] = [
+      { symbol: "H", pos: [0, 0, 0]      },
+      { symbol: "H", pos: [0, 0, 0.7414] },
+      { symbol: "H", pos: [0, 0, 3.0]    },
+      { symbol: "H", pos: [0, 0, 3.7414] },
+    ];
+    // Plain HF reference.
+    const { shells, nuclei, nElectrons } = moleculeToShellsNuclei(atoms, "sto-3g");
+    const integrals = computeMolecularIntegrals(shells, nuclei);
+    const plain = runRHFSCF(integrals, nElectrons, {
+      useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8,
+    });
+    expect(plain.converged).toBe(true);
+
+    const cp = runCounterpoise(
+      atoms,
+      [{ atomIndices: [0, 1] }, { atomIndices: [2, 3] }],
+      "sto-3g",
+      { useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8 },
+    );
+
+    expect(Math.abs(cp.supermoleculeEnergy - plain.energy)).toBeLessThan(1e-9);
+  });
+
+  test("single-fragment counterpoise reduces to plain HF (no ghosts)", () => {
+    // If you put every atom in ONE fragment, there are no ghosts —
+    // the "fragment in dimer basis" run is identical to the bare
+    // monomer run, and BSSE = 0.
+    const atoms: Atom[] = [
+      { symbol: "H", pos: [0, 0, 0]      },
+      { symbol: "H", pos: [0, 0, 0.7414] },
+    ];
+    const cp = runCounterpoise(
+      atoms,
+      [{ atomIndices: [0, 1] }],
+      "sto-3g",
+      { useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8 },
+    );
+    expect(cp.allConverged).toBe(true);
+    // With only one fragment: ΔE = ΔE_CP, so BSSE = 0.
+    expect(cp.bsseCorrection).toBeLessThan(1e-9);
+    // Fragment-with-ghosts == fragment-bare when there are no other atoms.
+    expect(Math.abs(cp.fragmentEnergiesCP[0]! - cp.fragmentEnergiesBare[0]!)).toBeLessThan(1e-9);
+  });
+
+  test("ghost-augmented fragment energy is lower (variational basis-set expansion)", () => {
+    // E(fragment in dimer basis) ≤ E(fragment bare) by the
+    // variational principle — adding basis functions can only
+    // lower (or equal) the SCF energy.
+    const atoms: Atom[] = [
+      { symbol: "H", pos: [0, 0, 0]      },
+      { symbol: "H", pos: [0, 0, 0.7414] },
+      { symbol: "H", pos: [0, 0, 3.0]    },
+      { symbol: "H", pos: [0, 0, 3.7414] },
+    ];
+    const cp = runCounterpoise(
+      atoms,
+      [{ atomIndices: [0, 1] }, { atomIndices: [2, 3] }],
+      "sto-3g",
+      { useDIIS: true, maxIter: 200, energyTol: 1e-10, densityTol: 1e-8 },
+    );
+    expect(cp.allConverged).toBe(true);
+    for (let k = 0; k < 2; k++) {
+      // ghost-augmented ≤ bare (variational; tiny shift in STO-3G).
+      expect(cp.fragmentEnergiesCP[k]!).toBeLessThanOrEqual(cp.fragmentEnergiesBare[k]! + 1e-12);
+    }
+  });
+
+  test("validates fragmentation: throws on overlap + on missing atoms", () => {
+    const atoms: Atom[] = [
+      { symbol: "H", pos: [0, 0, 0]      },
+      { symbol: "H", pos: [0, 0, 0.7414] },
+    ];
+    // Overlap.
+    expect(() => runCounterpoise(
+      atoms,
+      [{ atomIndices: [0, 1] }, { atomIndices: [1] }],
+      "sto-3g",
+    )).toThrow(/multiple fragments/);
+    // Missing.
+    expect(() => runCounterpoise(
+      atoms,
+      [{ atomIndices: [0] }],
+      "sto-3g",
+    )).toThrow(/not assigned/);
+    // Out of range.
+    expect(() => runCounterpoise(
+      atoms,
+      [{ atomIndices: [0, 1, 5] }],
+      "sto-3g",
+    )).toThrow(/out of range/);
+  });
+});