Multireference: CASSCF, CASPT2, NEVPT2¶

Some molecules are not honestly described by a single electron configuration. Stretch a bond toward dissociation, open a diradical, expose a near-degenerate transition-metal d-shell, or ask for an excited state, and the wavefunction becomes an unavoidable mixture of several determinants of comparable weight. Hartree-Fock, built on one determinant, gets these qualitatively wrong, and patching it with a single-reference correlation method on top (MP2, CCSD) inherits the error. The cure is a wavefunction that holds the competing configurations on an equal footing from the start: the complete active space (CAS) family. This page walks the standard ladder, RHF, then CASCI, then orbital-optimized CASSCF, then the CASPT2 and NEVPT2 dynamic-correlation corrections, on real systems with real numbers, and the Theory section explains what each rung buys you and why.

Every method here runs through the same run_job entry point you already know, selected with method= and an active_space=(n_orb, n_elec) tuple. The lower-level vibeqc.solvers.{casci,casscf,nevpt2, caspt2} functions are available when you need finer control (state averaging, custom references); the non-mean-field solvers tutorial covers that API and the active-space contract in depth.

Setup¶

The CAS methods need nothing beyond a molecule and a basis. A CAS calculation is specified by two numbers: how many electrons and how many spatial orbitals you place in the active space. vibe-qc takes them as active_space=(n_orb, n_elec) and partitions the RHF orbitals for you: the lowest (n_total_elec - n_elec) / 2 orbitals become a doubly-occupied frozen core, the next n_orb are active, and the rest are dropped. The frozen core is dressed exactly (its mean field enters the active one-electron integrals and a constant E_core is folded into the energy offset), so the number run_job reports is a real total energy you can put on a curve, not an active-space-only fragment. The active-space section of the solvers tutorial spells out that dressing.

from vibeqc import Atom, Molecule
from vibeqc.runner import run_job

# H2 at equilibrium (bond length 1.4 bohr)
h2 = Molecule([Atom(1, [0.0, 0.0, 0.0]), Atom(1, [0.0, 0.0, 1.4])])

result = run_job(h2, basis="sto-3g", method="casscf", active_space=(2, 2))
print(f"E(CASSCF(2,2)) = {result.energy:.8f} Ha")

E(CASSCF(2,2)) = -1.13727594 Ha

For H2 in a minimal basis the full valence space is CAS(2,2), the bonding \(\sigma_g\) and antibonding \(\sigma_u^*\) orbitals with both electrons, so this CASSCF is the exact (full-CI) STO-3G answer. Every result object exposes .energy (the total energy), .method (a label recording the active space and reference), and .converged; the PT2 and state-averaged variants add .root_energies for the per-state energies. We will read those off as we go.

When comparing dense and selected CI at the exact-selection limit, use identical orbitals, the same active-space Hamiltonian and the same target state. CASSCF also optimizes the orbitals: independent runs can reach different stationary solutions with different energies or iteration histories. Exact CI agreement at fixed orbitals does not establish a common orbital-optimization minimum or the global CASSCF minimum. The regression comparison corrected in #56 uses this fixed-orbital boundary.

The full ladder on N2¶

N2 at its equilibrium bond length is the cleanest illustration of the ladder, because vibe-qc’s CASSCF and CASPT2 energies for it are pinned to an external OpenMolcas reference (recorded in tests/test_solvers_mrpt_parity.py), so you can trust every digit. The honest active space for the nitrogen triple bond is CAS(6,6): the six electrons of the \(\sigma\) and two \(\pi\) bonds in the six valence orbitals \(\sigma_{2p}\), \(\pi_{2p}\times 2\), \(\pi^*_{2p}\times 2\), \(\sigma^*_{2p}\). The 1s cores and the 2s-derived \(\sigma/\sigma^*\) pair are the inactive frozen core.

This script climbs the whole ladder. It runs RHF, CASCI (CI in the frozen HF orbitals), CASSCF (the same CI but with the orbitals relaxed to minimize the CAS energy), then NEVPT2 and CASPT2 on the orbital-optimized CASSCF reference, which is the standard production composition.

from vibeqc import Atom, Molecule
from vibeqc.runner import run_job
from vibeqc.solvers import CASSCFOptions

# N2 at equilibrium, R = 2.074 bohr (1.098 Angstrom)
n2 = Molecule([Atom(7, [0.0, 0.0, 0.0]), Atom(7, [0.0, 0.0, 2.074])])

def energy(method, **kw):
    return run_job(n2, basis="sto-3g", method=method, **kw).energy

e_rhf    = energy("rhf")
e_casci  = energy("casci",  active_space=(6, 6))
e_casscf = energy("casscf", active_space=(6, 6))

# PT2 on the orbital-optimized CASSCF reference (pass casscf_options to
# switch the reference from CASCI-on-HF-orbitals to a converged CASSCF).
e_nevpt2 = energy("nevpt2", active_space=(6, 6), casscf_options=CASSCFOptions())
e_caspt2 = energy("caspt2", active_space=(6, 6), casscf_options=CASSCFOptions())

print(f"  RHF              {e_rhf:.8f}")
print(f"  CASCI(6,6)       {e_casci:.8f}   (static, frozen HF orbitals)")
print(f"  CASSCF(6,6)      {e_casscf:.8f}   (static, optimized orbitals)")
print(f"  NEVPT2 / CASSCF  {e_nevpt2:.8f}   (+ dynamic, intruder-free)")
print(f"  CASPT2 / CASSCF  {e_caspt2:.8f}   (+ dynamic)")
print()
print(f"  static  correlation (RHF -> CASSCF): {(e_rhf - e_casscf) * 627.509:7.2f} kcal/mol")
print(f"  dynamic correlation (CASSCF -> CASPT2): {(e_casscf - e_caspt2) * 627.509:7.2f} kcal/mol")

This is exactly what you get:

  RHF              -107.49584213
  CASCI(6,6)       -107.62174545   (static, frozen HF orbitals)
  CASSCF(6,6)      -107.63683517   (static, optimized orbitals)
  NEVPT2 / CASSCF  -107.64546059   (+ dynamic, intruder-free)
  CASPT2 / CASSCF  -107.64834004   (+ dynamic)

  static  correlation (RHF -> CASSCF):  88.47 kcal/mol
  dynamic correlation (CASSCF -> CASPT2):  7.22 kcal/mol

Read the ladder top to bottom. Each rung lowers the energy, and each lowering has a distinct physical meaning:

RHF to CASCI (\(-107.4958 \to -107.6217\), about 79 kcal/mol): letting the six active electrons distribute over the bonding and antibonding valence orbitals, instead of forcing them all into the bonds. This is the bulk of the static correlation, the part a single determinant structurally cannot hold.
CASCI to CASSCF (\(-107.6217 \to -107.6368\), about 9 kcal/mol): relaxing the orbitals themselves. The HF orbitals were optimal for the single HF determinant, not for the multi-determinant CAS wavefunction. Re-optimizing them under the CAS density recovers the rest of the static correlation.
CASSCF to NEVPT2 / CASPT2 (about 5 to 7 kcal/mol): the dynamic correlation, the short-range, many-small-contributions avoidance from all the excitations outside the active space. This is what a perturbation theory on top of a good reference is for.

The two external checks: vibe-qc’s CASSCF(6,6) is \(-107.63683517\) Ha against the OpenMolcas RASSCF value of \(-107.63683522\) Ha (agreement to 0.05 microhartree), and the CASPT2-on-CASSCF total is \(-107.64834004\) Ha against OpenMolcas \(-107.64834018\) Ha (0.14 microhartree). Those OpenMolcas numbers were generated out-of-process and recorded as constants in the test suite (the project rules in AGENTS.md disallow importing another quantum-chemistry program into vibe-qc), and the parity is what lets this tutorial state real digits rather than hand-wave at a trend.

Two references, two PT2 totals

By default method="caspt2" and method="nevpt2" build their CAS reference as CASCI on the HF orbitals (matching OpenMolcas RASSCF ... CIonly). Passing casscf_options=CASSCFOptions() switches the reference to an orbital-optimized CASSCF, the standard workflow and the one pinned to OpenMolcas above. The energies differ because the reference differs: on this N2 case the CASCI-referenced CASPT2 is \(-107.64577302\) Ha, the CASSCF-referenced one \(-107.64834004\) Ha. Pick the CASSCF reference unless you have a specific reason to perturb the frozen-orbital CASCI.

Why a single determinant fails: H2 dissociation¶

The ladder above was run at equilibrium, where the static correlation is real but moderate. To see why the multireference machinery is not optional, pull a bond apart and watch RHF break. H2 in a 6-31G basis is the smallest honest example: the full valence active space is again CAS(2,2), the smallest space that dissociates correctly, and RHF is the single-determinant approximation to it. CASSCF(2,2) is not the exact-in-basis curve (6-31G gives H2 four orbitals, so full CI would be CAS(2,4)), but it holds the static correlation that RHF misses. This script scans the bond from equilibrium out toward dissociation and prints the gap.

from vibeqc import Atom, Molecule
from vibeqc.runner import run_job

def energy(mol, method, **kw):
    return run_job(mol, basis="6-31g", method=method, **kw).energy

print(f" {'R/bohr':>7s} {'E(RHF)':>14s} {'E(CASSCF)':>14s} {'gap/kcal':>9s}")
for R in [1.4, 2.0, 3.0, 4.0, 5.0, 6.0]:
    h2 = Molecule([Atom(1, [0.0, 0.0, 0.0]), Atom(1, [0.0, 0.0, R])])
    e_rhf = energy(h2, "rhf")
    e_cas = energy(h2, "casscf", active_space=(2, 2))
    print(f" {R:>7.2f} {e_rhf:>14.8f} {e_cas:>14.8f} {(e_rhf - e_cas) * 627.509:>9.1f}")

  R/bohr         E(RHF)      E(CASSCF)  gap/kcal
40    -1.12674270    -1.14624789      12.2
00    -1.08381072    -1.11498960      19.6
00    -0.98167091    -1.04404698      39.1
00    -0.90055091    -1.00937107      68.3
00    -0.84318557    -0.99927739      97.9
00    -0.80456891    -0.99705370     120.8

Two things to read off. First, the CASSCF curve dissociates correctly: it flattens toward \(\approx -0.997\) Ha, which is twice the energy of a single hydrogen atom in this basis, the physically correct limit of two neutral H atoms. Second, the RHF curve does not: it keeps climbing, and at 6 bohr it sits 120.8 kcal/mol above the true limit. The error is not a small quantitative miss, it is the wrong shape of the curve, the RHF dissociation catastrophe.

The reason is structural. As the bond stretches, the exact ground state becomes an equal mixture of the closed-shell determinant (\(\sigma_g^2\)) and the doubly-excited one (\(\sigma_u^{*2}\)). RHF is locked into keeping \(\sigma_g\) doubly occupied, which forces in high-energy ionic terms (H\(^+\cdots\)H\(^-\)) that should vanish at dissociation. No single determinant can represent the correctly-mixed state, and so low-order single-reference corrections on top of RHF cannot fix it (MP2 diverges here). CCSD escapes only because H2 has two electrons, where CCSD equals full CI; for many-electron bond breaking it fails as well. CAS(2,2) holds both determinants from the start, which is why its curve is right at every distance. The gap between the two curves is the static correlation, and the fact that it keeps growing as you stretch the bond (toward a finite limit at infinite separation in this basis) is the fingerprint that tells you a single-reference method has no business being used.

At equilibrium (\(R = 1.4\)) the CASSCF(2,2) value \(-1.14624789\) Ha matches the OpenMolcas reference of \(-1.14624788\) Ha in the test suite, the same external pinning as the N2 case.

A water example with CASPT2 parity¶

For a polyatomic that is not bond-breaking, H2O at equilibrium in STO-3G with CAS(4,4) gives a compact ladder whose CASPT2 total is again OpenMolcas-pinned. Here the static correlation is small (water is a well-behaved closed shell), so the rungs are close together and most of the recovered correlation is dynamic, the regime where a single CASPT2/NEVPT2 step does the heavy lifting.

from vibeqc import Atom, Molecule
from vibeqc.runner import run_job

h2o = Molecule([
    Atom(8, [0.0,  0.0,   0.0]),
    Atom(1, [0.0,  1.43, -0.93]),
    Atom(1, [0.0, -1.43, -0.93]),
])

def energy(method, **kw):
    return run_job(h2o, basis="sto-3g", method=method, **kw).energy

print(f"  RHF          {energy('rhf'):.8f}")
print(f"  CASCI(4,4)   {energy('casci',  active_space=(4, 4)):.8f}")
print(f"  CASSCF(4,4)  {energy('casscf', active_space=(4, 4)):.8f}")
print(f"  NEVPT2(4,4)  {energy('nevpt2', active_space=(4, 4)):.8f}")
print(f"  CASPT2(4,4)  {energy('caspt2', active_space=(4, 4)):.8f}")

  RHF          -74.94142025
  CASCI(4,4)   -74.94699163
  CASSCF(4,4)  -74.95274478
  NEVPT2(4,4)  -74.97146663
  CASPT2(4,4)  -74.97302709

The NEVPT2 and CASPT2 lines here use the default CASCI-on-HF reference (no casscf_options), which is what OpenMolcas reproduces with its CIonly keyword. vibe-qc’s CASPT2 total \(-74.97302709\) Ha matches the recorded OpenMolcas &CASPT2 value of \(-74.97302712\) Ha to 0.03 microhartree. Note that the dynamic-correlation step (CASCI to CASPT2, about 16 kcal/mol, since this CASPT2 is CASCI-referenced) is much larger than the static step (RHF to CASSCF, about 7 kcal/mol) here: at a closed-shell equilibrium geometry, most of the missing correlation is dynamic, and the active space is doing comparatively little. That is the opposite balance from stretched N2, and recognizing which regime you are in is half the skill of using these methods.

Theory¶

The sections below build the picture from the ground: what an active space is, the difference between static and dynamic correlation, why one determinant is not enough, what CASSCF’s orbital optimization adds over CASCI, and how the two PT2 corrections recover the rest, including the intruder-state difference that distinguishes CASPT2 from NEVPT2.

Static versus dynamic correlation¶

Electron correlation, the energy a single Slater determinant misses, comes in two physically distinct flavors, and the whole CAS workflow is organized around treating them separately.

Static (or strong, or non-dynamic) correlation appears when two or more electron configurations have comparable weight in the true wavefunction. It is a near-degeneracy effect: bond-breaking (the bonding and antibonding configurations cross), diradicals (two near-degenerate singly-occupied orbitals), transition-metal complexes (partially-filled d-shells), and many excited states. Static correlation is large, it is qualitative, and no single determinant can capture it, because the very thing that makes it static is that several determinants matter at once.

Dynamic correlation is the short-range, instantaneous avoidance of electrons, the Coulomb hole. It is the sum of an enormous number of small contributions from excitations into the full virtual space. It is present in every system, including closed-shell molecules at equilibrium where there is no static correlation at all. Dynamic correlation is what MP2 and coupled cluster recover on top of a single good reference.

The division of labor in the CAS family follows this split exactly: a CAS wavefunction handles static correlation rigorously (by including all configurations within a chosen active space), and a perturbation theory on top of it (CASPT2, NEVPT2) recovers the dynamic correlation from outside the active space. The art is choosing an active space large enough to capture the static correlation but small enough to be affordable.

The active space¶

A complete-active-space wavefunction partitions the molecular orbitals into three sets:

Inactive (core): doubly occupied in every determinant. These orbitals are frozen, they contribute a constant energy and a mean field, but no correlation.
Active: the orbitals over which the wavefunction is allowed to correlate freely. The CAS ansatz includes all Slater determinants obtained by distributing the active electrons among the active orbitals in all symmetry-allowed ways, a full CI within the active space.
Virtual (secondary): empty in every determinant.

A CAS(\(n\), \(m\)) calculation places \(m\) electrons in \(n\) active spatial orbitals. The number of active determinants grows combinatorially, so active spaces are kept small (typically up to CAS(14,14) or so for conventional determinant CI, Knowles, Peter J. and Handy, Nicholas C. (1984). *A new determinant-based full configuration interaction method.* Chemical Physics Letters, **111**(4), 315–321. doi:[10.1016/0009-2614(84)85513-X](https://doi.org/10.1016/0009-2614(84)85513-X); larger needs DMRG, White, Steven R. (1992). *Density matrix formulation for quantum renormalization groups.* Physical Review Letters, **69**(19), 2863–2866. doi:[10.1103/PhysRevLett.69.2863](https://doi.org/10.1103/PhysRevLett.69.2863), Chan, Garnet Kin-Lic and Head-Gordon, Martin (2002). *Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group.* Journal of Chemical Physics, **116**(11), 4462–4476. doi:[10.1063/1.1449459](https://doi.org/10.1063/1.1449459)). Choosing the active space is the one genuinely non-black-box decision in a multireference calculation: it must contain every orbital involved in the near-degeneracy you care about. For a breaking bond, that is the bonding/antibonding pair; for the N2 triple bond, all three bonding orbitals and their three antibonds; for a \(\pi\) system, the full set of \(\pi\) and \(\pi^*\) orbitals.

CASCI: full CI in a chosen window¶

The cheapest CAS method is CASCI (Roos, Björn O., Taylor, Peter R., and Siegbahn, Per E. M. (1980). *A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach.* Chemical Physics, **48**(2), 157–173. doi:[10.1016/0301-0104(80)80045-0](https://doi.org/10.1016/0301-0104(80)80045-0)): take a fixed set of orbitals (here, the converged RHF orbitals), partition them into inactive/active/virtual, and diagonalize the Hamiltonian in the full determinant space of the active electrons. The inactive orbitals are frozen but dressed: they contribute an effective one-electron potential to the active orbitals plus a constant \(E_\text{core}\), so the energy CASCI reports is a real total energy. CASCI is variational (\(E_\text{CASCI} \le E_\text{HF}\) when the HF determinant lies inside the active space, and \(E_\text{CASCI} \ge E_\text{FCI}\) always), and it captures the static correlation that lives in the active space. What it does not do is improve the orbitals, it inherits whatever orbitals you hand it.

CASSCF: optimizing the orbitals too¶

The HF orbitals are optimal for the single HF determinant. They are not optimal for a multi-determinant CAS wavefunction, the density has changed, so the best orbitals have changed. CASSCF (complete active space self-consistent field, Roos, Björn O., Taylor, Peter R., and Siegbahn, Per E. M. (1980). *A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach.* Chemical Physics, **48**(2), 157–173. doi:[10.1016/0301-0104(80)80045-0](https://doi.org/10.1016/0301-0104(80)80045-0)) closes that loop: it optimizes the CI coefficients and the orbital shapes simultaneously, minimizing the CAS energy with respect to both. Formally it iterates two coupled steps, a CI solve in the current orbitals and an orbital rotation \(\mathbf{C} \to \mathbf{C}\,e^{\boldsymbol\kappa}\) driven by the orbital gradient, until the gradient vanishes.

The payoff is that CASSCF gives the best possible CAS-quality wavefunction for a given active space, and a reference whose orbitals are self-consistent with the correlated density, which matters when you then put perturbation theory on top. The cost is that the orbital optimization is non-linear and can land in different stationary points from different starting orbitals (the N2 and H2 cases here have unique minima; water’s CAS(4,4) is famously basin-rich, which is why the test suite asserts only a variational bound for it rather than exact OpenMolcas equality). In the ladders above, the CASCI-to-CASSCF drop is precisely the static correlation that frozen HF orbitals leave on the table.

CASPT2 and NEVPT2: the dynamic-correlation correction¶

A CASSCF wavefunction is an excellent reference but still misses dynamic correlation, the excitations into the virtual space it never touches. Both CASPT2 and NEVPT2 recover that dynamic correlation by second-order perturbation theory built on the CAS reference, in direct analogy to how 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)) corrects a single HF determinant. The first-order wavefunction spans the excitations out of the active space (internal-to-virtual, internal-to-active, and so on), and the second-order energy is the resulting correction. Both reduce to MP2 exactly when the active space is empty, and both give zero correction when the active space is the full space (nothing left outside to correlate), limits the vibe-qc test suite checks explicitly.

The two methods differ in their choice of zeroth-order Hamiltonian \(\hat H_0\), the operator whose eigenvalues become the energy denominators of the perturbation series, and that difference has a practical consequence:

CASPT2 (Andersson, Malmqvist, Roos; Andersson, Kerstin, Malmqvist, Per-Åke, et al. (1990). *Second-order perturbation theory with a CASSCF reference function.* Journal of Physical Chemistry, **94**(14), 5483–5488. doi:[10.1021/j100377a012](https://doi.org/10.1021/j100377a012), Andersson, Kerstin, Malmqvist, Per-Åke, and Roos, Björn O. (1992). *Second-order perturbation theory with a complete active space self-consistent field reference function.* Journal of Chemical Physics, **96**(2), 1218–1226. doi:[10.1063/1.462209](https://doi.org/10.1063/1.462209)) uses a generalized Fock operator as \(\hat H_0\). It is the established workhorse, but its \(\hat H_0\) is not guaranteed to keep the perturber states above the reference in energy. When a perturber accidentally falls near the reference energy, its denominator goes to zero and the second-order energy diverges, the intruder-state problem. The standard remedies are a real or imaginary level shift that pushes the offending denominator away from zero; vibe-qc exposes an imaginary shift (Forsberg, Niclas and Malmqvist, Per-Åke (1997). *Multiconfiguration perturbation theory with imaginary level shift.* Chemical Physics Letters, **274**(1–3), 196–204. doi:[10.1016/S0009-2614(97)00669-6](https://doi.org/10.1016/S0009-2614(97)00669-6)) through caspt2_options=CASPT2Options(imaginary=...).
NEVPT2 (Angeli, Cimiraglia, Malrieu; Angeli, Celestino, Cimiraglia, Renzo, et al. (2001). *Introduction of n-electron valence states for multireference perturbation theory.* Journal of Chemical Physics, **114**(23), 10252–10264. doi:[10.1063/1.1361246](https://doi.org/10.1063/1.1361246)) uses Dyall’s \(\hat H_0\) (K. G. Dyall, J. Chem. Phys. 102, 4909 (1995), doi:10.1063/1.469539), which retains the two-electron interaction within the active space. This is more expensive per term but structurally intruder-free: the perturber denominators are bounded away from zero by construction, so NEVPT2 needs no level shift and behaves smoothly even on difficult references. That robustness is its main selling point over CASPT2.

On well-behaved references the two agree closely, in the N2 ladder above, NEVPT2/CASSCF and CASPT2/CASSCF differ by under 2 kcal/mol, a typical spread. NEVPT2 is the safer default when you suspect intruder trouble; CASPT2 remains a common, widely benchmarked choice for excited states and spectroscopy. vibe-qc ships internally-contracted CASPT2 un-gated and validated against OpenMolcas to within a few microhartree, including the canonical eight-class single-state IPEA contraction (Ghigo, Giovanni, Roos, Björn O., and Malmqvist, Per-Åke (2004). *A modified definition of the zeroth-order Hamiltonian in multiconfigurational perturbation theory (CASPT2).* Chemical Physics Letters, **396**(1–3), 142–149. doi:[10.1016/j.cplett.2004.08.032](https://doi.org/10.1016/j.cplett.2004.08.032), Nishimoto, Yoshio (2023). *Analytic First-Order Derivatives of CASPT2 with IPEA Shift.* Journal of Chemical Physics, **158**(17), 174112. doi:[10.1063/5.0147611](https://doi.org/10.1063/5.0147611)). The strongly-contracted CASPT2 variant remains gated behind VIBEQC_EXPERIMENTAL_MULTIREF=1.

Status: gated

The strongly-contracted CASPT2 variant (CASPT2Options(variant="sc"), the 1990 single-state formulation) is recognised but fails closed. Without the opt-in it raises RuntimeError: CASPT2 (strongly-contracted) is gated: it is not sub-microhartree reproducible against OpenMolcas, naming the missing validation, and runs only after you set the environment variable VIBEQC_EXPERIMENTAL_MULTIREF=1. With the gate lifted, H2/6-31G CAS(2,2) at 1.6 bohr gives -1.145309 Ha, against -1.145330 Ha for the default internally contracted variant. Only this variant is gated: the default method="caspt2" (internally contracted) and method="nevpt2" (strongly contracted by construction, validated against PySCF) are production. See Methods: theory and status.

Choosing a method¶

A short decision guide, since the ladder offers several stopping points:

If you are at or near a closed-shell equilibrium and just want the correlation energy, you do not need a multireference method at all, MP2 or CCSD on the RHF reference is cheaper and entirely appropriate (see the non-HF solvers tutorial).
If more than one configuration matters, bond-breaking, diradicals, transition-metal d-shells, excited states, reach for the CAS family. Use CASCI for a quick look in fixed orbitals, CASSCF for a proper variational reference, and add NEVPT2 or CASPT2 when you need the dynamic correlation on top.
Between the two PT2s: prefer NEVPT2 when intruder states are a risk (stretched geometries, dense excited-state manifolds) because it needs no shift; prefer CASPT2 when you want to match the large body of CASPT2 benchmark literature, and reach for its imaginary shift if a denominator misbehaves.

Nuclear gradients and geometry optimization¶

A converged state-specific CASSCF wavefunction has a well-defined analytic nuclear gradient, and vibe-qc’s is validated for that envelope: closed shell, one root. Outside it (state-averaged CASSCF with nroots > 1, open-shell references, or compute_wz="numerical") the analytic vector is not checked against finite differences, and the geometry optimizers fall back to central finite differences of the full energy. The PT2 methods are narrower still, as described after the optimisation example below.

run_job attaches the gradient to the result whenever a CASSCF single point converges:

from vibeqc.solvers import CASSCFOptions

result = run_job(mol, basis="cc-pvdz", method="casscf", active_space=(6, 6))
forces = -result.gradient            # (n_atoms, 3), Hartree/bohr

The single-state gradient is the complete analytic derivative of the CASSCF energy: a converged state-specific CASSCF is stationary in every orbital rotation and CI coefficient, so no response (“z-vector”) term exists, and the analytic vector matches a central finite difference of the energy to ~2e-7 Ha/bohr (H2/6-31G and H2/cc-pVDZ CAS(2,2)) and to 1.1e-6 on N2/cc-pVDZ CAS(6,6), where the residual is the finite difference’s own truncation error. Re-checked on a vibe-qc 0.17.10 development build: H2/6-31G CAS(2,2) at 1.8 bohr gives 0.06979236 Ha/bohr analytically against 0.06979230 by central differences, and N2/STO-3G CAS(6,6) at 2.3 bohr gives 0.0446371 against 0.0446365. A finite-difference cross-check is available as CASSCFOptions(compute_wz="numerical") (expensive: the CASSCF is re-converged at every displaced geometry).

CASSCFOptions(compute_wz=True), which once selected an experimental “W^z” CP-MCSCF correction, is now an accepted no-op alias of the analytic gradient and emits a FutureWarning: that correction was a phantom term for a variational energy and produced spurious gradient components on linear molecules with p/d functions (GitLab #516). State-averaged runs (nroots > 1 with weights) expose the state-averaged analytic approximation, which is checked for finiteness and translational invariance but not against finite differences.

For geometry optimization, call optimize_molecule. Inside the validated envelope it walks the analytic gradient; outside it, central finite differences of the full CASSCF energy. It returns the trajectory of frames:

from vibeqc import Atom, Molecule
from vibeqc.molecular_optimize import optimize_molecule
from vibeqc.solvers import CASSCFOptions

h2 = Molecule([Atom(1, [0, 0, 0]), Atom(1, [0, 0, 1.5])])   # stretched, bohr
opt = optimize_molecule(
    h2, "6-31g", method="casscf", active_space=(2, 2),
    casscf_options=CASSCFOptions(), conv_tol_grad=1e-2, max_iter=20,
)
print(opt.energy, opt.trajectory_frames[-1])     # R = 1.5 relaxes toward ~1.4 bohr

run_job(..., optimize=True) makes the same choice. For the H2/6-31G example above, run through run_job from a 1.5 bohr start, the ASE, native and brent optimizer backends all relax to E = -1.14633791 Ha. The regression adjudicator for the analytic gradient is examples/regression/casscf_gradient_fd_reproducer.py.

Gradients and optimisation for the PT2 methods¶

The analytic PT2 gradient (Celani, Paolo and Werner, Hans-Joachim (2003). *Analytical energy gradients for internally contracted second-order multireference perturbation theory.* Journal of Chemical Physics, **119**(10), 5044–5057. doi:[10.1063/1.1597672](https://doi.org/10.1063/1.1597672)) covers only a narrow envelope: state-specific, closed-shell, unshifted, explicit-engine IC-CASPT2 on a converged CASSCF reference (variant="ic", ipea=0, imaginary=0, gas phase, with compute_corr_grad=True on the options); NEVPT2 gets a relaxed finite-difference gradient under the same switch. Everything else is finite differences of the full energy, which is a preview capability. With the default compute_corr_grad=False the returned gradient is the bare CASSCF gradient, and vibe-qc warns that it is not a PT2 force.

run_job(..., optimize=True) for fci, nevpt2, caspt2 and mrci therefore refuses the default optimizer backend: the ASE backend has no correlated-surface calculator for them and would silently relax the mean-field surface while reporting the correlated method. It raises ValueError: run_job(method='caspt2', optimize=True) requires optimizer_backend='native' or 'brent' (GitLab #357). Pass optimizer_backend="native" or "brent" to walk the correlated surface by finite differences of the full energy, and expect minutes, not seconds, even for a small molecule:

result = run_job(
    h2, basis="6-31g", method="caspt2", active_space=(2, 2),
    casscf_options=CASSCFOptions(),
    optimize=True, optimizer_backend="native",
)

Natural orbitals and QVF wavefunction export¶

A converged CASSCF/CASCI result has no ordinary canonical mean-field MOs to export, but it does have natural orbitals, the active-space 1-RDM’s eigenvectors, with fractional occupation numbers that are exactly the multireference character the calculation was run to find. run_job computes these automatically and, with the default output_qvf=True, embeds them as the QVF archive’s wavefunction section, so vibe-view can resample any natural orbital on demand just as it would a mean-field MO:

result = run_job(
    mol, basis="cc-pvdz", method="casscf", active_space=(6, 6),
    output="n2_casscf",   # -> n2_casscf.out, n2_casscf.qvf
)

vibe-view open n2_casscf.qvf   # browse the natural orbitals + occupations

No extra flag is needed: output_qvf already defaults to True. This is the same packaging path mean-field routes use for their MOs; CASSCF/CASCI results are detected because they populate natural_orbitals / natural_occupations on the returned SolverResult where a mean-field result would populate mo_energies instead.

For a second worked example (H2/6-31G/CAS(2,2), the textbook two-configuration case), see examples/vibe_view/showcase_natural_orbitals.py, which regenerates the fixture pinning vibe-view’s HONO/LUNO path against computed (not invented) data.

References¶

The papers below define the methods used on this page. vibe-qc emits the matching citations automatically in the references block of any run that uses them (the routing lives in python/vibeqc/output/citations/database.toml).

CASCI / CASSCF. B. O. Roos, P. R. Taylor, and P. E. M. Siegbahn, “A complete active space SCF method (CASSCF) using a density matrix formulated super-CI approach,” Chem. Phys. 48, 157 (1980), doi:10.1016/0301-0104(80)80045-0. The original density-matrix super-CI CASSCF; vibe-qc cites it for both method="casci" and method="casscf" and as the CAS reference for the PT2 methods.
CASPT2. K. Andersson, P.-Å. Malmqvist, and B. O. Roos, “Second-order perturbation theory with a complete active space self-consistent field reference function,” J. Chem. Phys. 96, 1218 (1992), doi:10.1063/1.462209. The internally-contracted CASPT2 that is vibe-qc’s default method="caspt2". The original single-state formulation is K. Andersson, P.-Å. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, J. Phys. Chem. 94, 5483 (1990), doi:10.1021/j100377a012 (vibe-qc’s gated variant="sc").
NEVPT2. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, “Introduction of n-electron valence states for multireference perturbation theory,” J. Chem. Phys. 114, 10252 (2001), doi:10.1063/1.1361246. Defines NEVPT2 with Dyall’s intruder-free zeroth-order Hamiltonian; vibe-qc implements the strongly-contracted variant as method="nevpt2".
IPEA shift (CASPT2). G. Ghigo, B. O. Roos, and P.-Å. Malmqvist, “A modified definition of the zeroth-order Hamiltonian in multiconfigurational perturbation theory (CASPT2),” Chem. Phys. Lett. 396, 142 (2004). The empirical level shift available via caspt2_options=CASPT2Options(ipea=...). The canonical internally contracted treatment follows Y. Nishimoto, “Analytic First-Order Derivatives of CASPT2 with IPEA Shift,” J. Chem. Phys. 158, 174112 (2023), doi:10.1063/5.0147611.
Textbook. A. Szabo and N. S. Ostlund, Modern Quantum Chemistry, Dover (1996), §3.8 walks the H2 dissociation failure of RHF that motivates the whole multireference picture. For CASSCF and multireference perturbation theory specifically, T. Helgaker, P. Jorgensen, and J. Olsen, Molecular Electronic-Structure Theory, Wiley (2000), chapters 12 and 14, give the complete derivations.

Next¶

Once the ladder makes sense, these are the natural follow-ups, deeper API control, runnable scripts, the single-reference methods for the easy regime, and convergence help when the orbital optimization misbehaves.

The non-mean-field solvers tutorial covers the lower-level vibeqc.solvers API, the active-space contract in full, state-averaged CASSCF, and the other determinant-based solvers (Selected-CI, DMRG, v2RDM) you can use as CAS references.
Runnable scripts for every method live in examples/wavefunction/; 01_casscf_h2o.py is the CASCI / CASSCF / PT2 / state-averaged walkthrough this page distills, and n2_dissociation_active_space.py scans the breaking N2 triple bond with RHF, CASCI(6,6), and CASSCF(6,6); 02_casscf_gradient.py and 03_casscf_geometry_optimization.py cover the analytic nuclear gradient and geometry optimization.
For the single-reference correlation methods appropriate near equilibrium (MP2, CCSD), see the molecular HF tutorial and the post-HF entries in the API reference.
If your CASSCF will not converge or lands in the wrong basin, the SCF convergence guide covers the orbital-optimization knobs.