CASCI / CASSCF / NEVPT2 / CASPT2 — Implementation Handover

Last updated: 2026-05-31 (CASSCF added; rebuilt + validated after the 2026-05-30 audit) Branch: main Status:

  • CASCI — fixed, validated vs PySCF, ungated.

  • CASSCF is orbital-optimized CASCI, validated vs PySCF mcscf.CASSCF, ungated.

  • NEVPT2 (strongly-contracted) — fixed, validated vs PySCF, ungated.

  • CASPT2 (strongly-contracted) — implemented, validated in its limits, gated (VIBEQC_EXPERIMENTAL_MULTIREF=1) pending an external reference.

What was wrong (audit 2026-05-30) and how it was fixed

These three methods landed in 0c8f0fcf as scaffolding on the v0.21 roadmap path. The audit found all three broken; they have since been rebuilt:

Method

Was

Now

CASCI

crashed for >1 same-spin e⁻; cored ≈ −104 Ha (below FCI)

rebuilt on the validated build_hamiltonian_matrix_unrestricted engine + correct frozen-core dressing; matches PySCF mcscf.CASCI

NEVPT2

3/8 classes; wrong-sign E_corr

strongly-contracted NEVPT2 via the uniform perturber-projection definition; matches PySCF mrpt.NEVPT per class

CASPT2

returned 0 for no-core

strongly-contracted CASPT2 (generalized-Fock H₀); MP2 limit exact, gated

Key design decisions

  • CASCI is built on the proven FCI engine. The only CASCI-specific physics is the frozen-core dressing (_casci.py::_frozen_core_dressing): E_core = ⟨core|H|core⟩ and h̃_pq = h_pq + Σ_c (2 g_{pcqc} − g_{pccq}) (physicist’s g). Validated to Δ=0 against an explicit CAS-determinant FCI over the full orbital space, and to ~1e-12 against PySCF.

  • RDMs (_rdm.py) are spin-summed 1/2/3-RDM builders in PySCF convention, validated elementwise vs fci.direct_spin1.make_rdm123 (~1e-15) and by the energy self-test energy_from_rdms == eigenvalue.

  • NEVPT2/CASPT2 share one engine (_mrpt.py). The strongly-contracted PT2 energy is computed from its definition — one perturber |V_g⟩ per external (core-hole / virtual-particle) pattern, E = Σ_g −⟨V_g|V_g⟩ / (⟨V_g|H₀|V_g⟩/⟨V_g|V_g⟩ − E₀) — using a small spin-orbital second-quantization engine, so all spin/antisymmetry and combinatorial factors are handled by construction (no per-class formulas). NEVPT2 uses Dyall’s H₀; CASPT2 uses the generalized-Fock H₀. Orbitals are semicanonicalized (core + virtual blocks of the generalized Fock diagonalized); the active CI vector is invariant under that rotation.

  • CASSCF optimizes the orbitals on top of CASCI (_casscf.py). Two-step macro-iteration: solve CASCI, build the non-symmetric generalized Fock and the orbital gradient g_pq = 2(F_qp − F_pq) over the non-redundant (inactive↔active / inactive↔virtual / active↔virtual) rotations, then take a backtracked augmented-Hessian step whose Hessian comes from central finite differences of that gradient. The gradient is exact (re-solving the CI in every basis removes the CI-response term) and is validated elementwise against finite differences of the trusted CASCI energy (~1e-9). The backtracking line search makes the energy monotone, so it cannot diverge. Like every CASSCF the surface is non-convex: the optimizer converges to the stationary point in the basin of the starting (HF) orbitals, which need not be the global minimum (PySCF’s default start can land in a different basin). Basis-independent correctness is pinned by the variational sandwich E_FCI ≤ E_CASSCF ≤ E_CASCI, a vanishing orbital gradient, and exact rotation-invariance of the full-space CASCI energy. NEVPT2/CASPT2 can run on the CASSCF reference via CASSCFResult.{cas, h1e_cas, h2e_cas}.

Validation (PySCF out-of-process, CLAUDE.md §10)

System

Quantity

vs reference

H₂O/STO-3G CAS(2,2)/(4,4), H₂/6-31g, LiH

CASCI

PySCF mcscf.CASCI, Δ ≤ 1e-9

H₂/6-31g, LiH, HF, Be (unique-minimum)

CASSCF

PySCF mcscf.CASSCF, Δ ≤ 1e-7

H₂/6-31g

CASSCF orbital gradient

finite diff of CASCI energy, Δ ≤ 1e-6

H₂O/STO-3G CAS(2,2)

CASSCF sandwich

E_FCI ≤ E_CASSCF ≤ E_CASCI

H₂/6-31g, H₂O/STO-3G CAS(2,2)/(4,4)

NEVPT2 (all 8 classes)

PySCF mrpt.NEVPT per class, Δ ≤ 5e-9

H₂/6-31g

CASPT2 empty-active (MP2 limit)

PySCF mp.MP2, Δ ≤ 1e-9

H₂/6-31g

CASPT2 full-space

0 (exact)

Tests: tests/test_solvers_mrpt_parity.py.

CASPT2 — why it stays gated

PySCF ships no CASPT2, so the active-space terms have no local external reference. The implemented variant is strongly-contracted CASPT2 with the generalized-Fock H₀ and no IPEA/level shift; the production reference codes (OpenMolcas, Molpro) use the internally-contracted variant with an IPEA shift, so their numbers will differ even where both are “CASPT2”. Note that ORCA is not a CASPT2 reference: ORCA’s multireference PT2 is NEVPT2, not CASPT2 (CASPT2 is the Roos/Andersson method native to OpenMolcas/Molpro), so ORCA can only cross-check NEVPT2 / CASSCF / CASCI here, not CASPT2. Per CLAUDE.md §7 the active-space energy is not presented as trustworthy: it is gated behind VIBEQC_EXPERIMENTAL_MULTIREF=1. To un-gate: validate against OpenMolcas/Molpro (or implement the internally-contracted variant) and flip the require_experimental_multiref("CASPT2") call in _mrpt.py::caspt2.

Files

File

Role

python/vibeqc/solvers/_rdm.py

spin-summed 1/2/3-RDM builders (new)

python/vibeqc/solvers/_casci.py

CASCI on the FCI engine + frozen-core dressing

python/vibeqc/solvers/_casscf.py

CASSCF: orbital optimization on CASCI (new)

python/vibeqc/solvers/_mrpt.py

spin-orbital engine + SC-NEVPT2 + SC-CASPT2

python/vibeqc/solvers/_common.py

require_experimental_multiref (CASPT2 gate)

python/vibeqc/runner.py

method="casci"/"casscf"/"nevpt2"/"caspt2" dispatch

python/vibeqc/memory.py

casci/casscf/nevpt2/caspt2 memory estimators

python/vibeqc/output/citations/database.toml

Roos 1980 (CASCI+CASSCF) / Angeli 2001 / Andersson 1990 + routes

tests/test_solvers_mrpt_parity.py

parity + limit + gate tests

Performance / scope

The spin-orbital engine scales steeply with the orbital count, so NEVPT2 / CASPT2 are for small active spaces + small bases (consistent with the rest of vibeqc.solvers). CASCI scales as FCI in the active space. CASSCF rebuilds the orbital Hessian by finite differences of the analytic gradient each macro-iteration (O(n_pairs) CASCI solves per iteration), which is fine for the small active spaces this engine targets; it is not tuned for large CAS.

Next steps

Internally-contracted CASPT2 + IPEA is the active path to replace the strongly-contracted variant and un-gate. Build order, foundation first:

  1. 4-RDM (make_rdm1234 in _rdm.py): DONE 2026-06-01, validated vs PySCF make_rdm1234 to ~1e-14 (tests/test_solvers_rdm.py).

  2. Per-class overlap metric S_uv and the generalized-Fock H0 in the internally-contracted first-order interacting space (uses RDMs through 4th order). New machinery.

  3. Metric-orthonormalized linear solve for the per-class E(2); validate the MP2 and full-space limits first, then class-by-class against OpenMolcas.

  4. Real / imaginary level shifts (Roos-Andersson 1995 / Forsberg-Malmqvist 1997) and the IPEA shift (Ghigo 2004, default 0.25) drop in as diagonal terms on the orthonormalized denominators.

  5. Validation reference: OpenMolcas, now installed and verified (driven as a subprocess per CLAUDE.md section 10; ORCA cannot reference CASPT2, its MR-PT2 is NEVPT2). Cross-check confirmed on H2/6-31g CAS(2,2): vibe-qc CASSCF matches OpenMolcas RASSCF to ~1e-8 (same geometry/basis/active space, so the comparison is apples-to-apples). On that system the current strongly- contracted CASPT2 sits +3.6 mHa above OpenMolcas internally-contracted CASPT2 (IPEAshift=0.0); that gap is the correlation IC-CASPT2 must recover to un-gate. Validate IPEA=0 first, then the OpenMolcas default 0.25.

The strongly-contracted perturber engine in _mrpt.py is NOT reusable for the internally-contracted path (SC materializes one perturber per external pattern; IC solves a metric linear system). The RDMs (now through 4th order), the generalized Fock, and the semicanonicalization ARE reusable.

Also open: optimize the spin-orbital engine for larger CAS; an analytic CASSCF orbital Hessian; CASSCF state-averaging / excited states / open-shell.