Non-mean-field solvers — Selected-CI, DMRG, v2RDM, and transcorrelated methods

vibe-qc v0.9+ ships a family of non-mean-field wavefunction solvers. These go beyond Hartree–Fock and DFT by operating directly on the many-electron wavefunction — trading the single- determinant picture for systematically improvable representations.

All methods share a common foundation: a Hamiltonian object carrying one- and two-electron integrals in an orthonormal basis, and a solve(hamiltonian, options) → SolverResult protocol. You can run them through the low-level solver API or through the familiar run_job entry point — whichever fits your workflow.

This is an active development area. The Python/NumPy backends are teaching tools and correctness references for small-to-medium active spaces. Production backends (C++ / DMRG with quantum-number symmetries) are on the roadmap. But the APIs you learn here are the ones that will persist.

Quick start

The simplest possible non-mean-field calculation — Selected-CI on H₂O / STO-3G — through run_job:

import vibeqc as vq

mol = vq.Molecule([
    vq.Atom(8, [ 0.0,  0.00,  0.00]),
    vq.Atom(1, [ 0.0,  1.43, -0.98]),
    vq.Atom(1, [ 0.0, -1.43, -0.98]),
])

vq.run_job(
    mol,
    basis="sto-3g",
    method="selected_ci",
    selected_ci_options=vq.SelectedCIOptions(target_size=100, verbose=1),
    output="h2o_selci",
)

That’s it. vibe-qc runs RHF to get the MO basis, builds the Hamiltonian, runs Selected-CI with iterative determinant selection, and writes the result to h2o_selci.out. The whole call takes a few tens of milliseconds for this system.

For the low-level API (more control, returns a SolverResult object directly):

from vibeqc.solvers import (
    Hamiltonian,
    SelectedCIOptions,
    build_hamiltonian_mo,
    get_hf_orbital_provider,
    solve_selected_ci,
)
import vibeqc as vq

mol = vq.Molecule([
    vq.Atom(8, [ 0.0,  0.00,  0.00]),
    vq.Atom(1, [ 0.0,  1.43, -0.98]),
    vq.Atom(1, [ 0.0, -1.43, -0.98]),
])
basis = vq.BasisSet(mol, "sto-3g")

# 1. Get orthonormal orbitals (HF canonical)
C = get_hf_orbital_provider(mol, basis)

# 2. Build MO-basis Hamiltonian
H = build_hamiltonian_mo(mol, basis, C)

# 3. Solve
opts = SelectedCIOptions(target_size=100, verbose=1)
result = solve_selected_ci(H, opts)

print(f"E(Selected-CI) = {result.energy:.8f} Ha")
print(f"ndet           = {len(result.ci_labels)}")

What are non-mean-field methods?

In Hartree–Fock and DFT, the many-electron wavefunction is approximated as a single Slater determinant — each electron moves in the average field of the others. This is a mean-field approximation. It captures ~99% of the total electronic energy but misses the remaining ~1% of electron correlation: the instantaneous, correlated motion of electrons avoiding each other.

Non-mean-field methods recover this correlation by expanding the wavefunction as a linear combination of many determinants:

\[|\Psi\rangle = \sum_I c_I |D_I\rangle\]

where each \(|D_I\rangle\) is a Slater determinant (a specific occupation pattern of the molecular orbitals). The coefficients \(c_I\) are determined variationally by diagonalising the Hamiltonian in the determinant basis. Different methods differ in how they choose which determinants to include and how they solve for the coefficients — but they all share this many-determinant ansatz.

Available methods

Method

run_job keyword

What it does

Selected-CI

"selected_ci"

Iteratively builds a compact determinant space by adding the most important singles and doubles, then diagonalises. Optional PT2 correction.

Full CI (manual)

(low-level API)

Diagonalises the Hamiltonian over the complete determinant space — exact within the one-electron basis.

DMRG

"dmrg"

Matrix-product-state wavefunction optimised by two-site sweeps. Bond-dimension controls accuracy.

v2RDM

"v2rdm"

Minimises the energy directly as a functional of the 2-RDM, under N-representability constraints (P, Q, G).

Transcorrelated-CI

"transcorrelated_ci"

Applies a similarity transformation \(\tilde{H} = e^{-\tau} H e^{\tau}\) to fold short-range correlation into the Hamiltonian, then solves with Selected-CI.

All methods consume the same Hamiltonian object — you can swap solvers without rebuilding integrals.

Selected-CI walkthrough — H₂O convergence scan

Selected-CI is the workhorse. It starts with the HF determinant, iteratively adds the most important singles and doubles (ranked by their perturbative contribution), and re-diagonalises. As you increase target_size, the energy converges smoothly toward the FCI limit.

Here’s the full convergence scan for H₂O / STO-3G:

import time
import vibeqc as vq
from vibeqc.solvers import (
    SelectedCIOptions,
    build_hamiltonian_mo,
    get_hf_orbital_provider,
    solve_selected_ci,
)

# H₂O at the Pitzer/Pulay equilibrium geometry
import numpy as np
R_oh = 1.8089
theta = np.radians(104.52 / 2)
mol = vq.Molecule([
    vq.Atom(8, [0.0, 0.0, 0.0]),
    vq.Atom(1, [0.0, +R_oh * np.sin(theta), -R_oh * np.cos(theta)]),
    vq.Atom(1, [0.0, -R_oh * np.sin(theta), -R_oh * np.cos(theta)]),
])
basis = vq.BasisSet(mol, "sto-3g")

# RHF reference
C = get_hf_orbital_provider(mol, basis)
H = build_hamiltonian_mo(mol, basis, C)

rhf_result = vq.run_rhf(mol, basis, vq.RHFOptions())
print(f"E(RHF)  = {rhf_result.energy:.8f} Ha")
print(f"E_nuc   = {mol.nuclear_repulsion():.8f} Ha")
print(f"nbasis  = {basis.nbasis},  n_el = {mol.n_electrons()}")
print()

# ── Convergence scan: increasing target_size ──
header = f" {'target':>7s}  {'ndet':>6s}  {'E_var':>16s}  {'ΔE':>12s}  {'E_pt2':>12s}  {'E_tot':>16s}  {'ms':>6s}"
print(header)
print("-" * len(header))

prev_e_var = float("inf")
for target in [1, 3, 5, 10, 20, 50, 100, 200]:
    t0 = time.perf_counter()
    ci_opts = SelectedCIOptions(
        target_size=target,
        max_iter=30,
        conv_tol_energy=1e-10,
        do_pt2_correction=True,
        pt2_threshold=1e-8,
        verbose=0,
    )
    result = solve_selected_ci(H, ci_opts)
    dt = time.perf_counter() - t0
    ndet = len(result.ci_labels) if result.ci_labels else 0
    e_var = result.energy - (result.pt2_correction or 0.0)
    delta = e_var - prev_e_var if prev_e_var != float("inf") else 0.0
    prev_e_var = e_var
    pt2 = result.pt2_correction or 0.0
    print(
        f" {target:>7d}  {ndet:>6d}  {e_var:>16.10f}  {delta:>12.2e}  "
        f"{pt2:>12.2e}  {result.energy:>16.10f}  {dt * 1000:>5.0f}"
    )

What you should see:

E(RHF)  = -74.96306453 Ha

 target    ndet            E_var           ΔE        E_pt2           E_tot     ms
-----------------------------------------------------------------------------
      1       1   -74.9630645300   0.00e+00   0.00e+00   -74.9630645300      0
      3       3   -74.9652419726  -2.18e-03  -4.69e-02   -75.0121489285      1
      5       5   -74.9667159170  -1.47e-03  -4.54e-02   -75.0121315831      1
     10      10   -74.9686632147  -1.95e-03  -4.35e-02   -75.0121341488      1
     20      20   -74.9712110099  -2.55e-03  -4.09e-02   -75.0121318632      2
     50      50   -74.9765719143  -5.36e-03  -3.55e-02   -75.0121088645      4
    100     100   -74.9841836813  -7.61e-03  -2.79e-02   -75.0121058434      9
    200     200   -74.9912803868  -7.10e-03  -2.08e-02   -75.0120885498     19

A few things to notice:

  • target_size=1 gives you back the RHF energy — the reference determinant alone.

  • By ~10 determinants, you’ve captured 95% of the correlation energy.

  • The PT2 correction shrinks as the variational space grows — it’s estimating the energy from determinants you haven’t included yet, so it naturally decays as you include more.

  • The total energy (E_var + E_pt2) is remarkably stable across the whole scan — usually within ~0.1 mHa after the first few determinants. This is the hallmark of a well-behaved CIPSI selection.

How it works under the hood

Each Selected-CI iteration does three things:

  1. Selection — from the current wavefunction \(|\Psi\rangle\), generate all singles and doubles. For each candidate \(|D\rangle\), estimate its first-order perturbative weight:

    \[w_D \approx \frac{|\langle D|H|\Psi\rangle|^2}{E_{\rm var} - \langle D|H|D\rangle}\]

  2. Pruning — keep the top-N candidates (controlled by target_size and pt2_threshold).

  3. Diagonalisation — build and diagonalise the Hamiltonian in the expanded determinant space. This is the variational step; it guarantees an upper bound to the exact-in-basis energy.

The optional PT2 correction uses Epstein-Nesbet denominators to estimate the residual correlation from all singles and doubles outside the variational space.

Full CI — exact-in-basis for small molecules

Full CI diagonalises the Hamiltonian over the complete determinant space. It’s the exact solution within the given one-electron basis — the gold standard that Selected-CI, DMRG, and v2RDM approximate.

For H₂O / STO-3G (7 spatial orbitals, 10 electrons), the full determinant space has \(C(7,5) \times C(7,5) = 21 \times 21 = 441\) Slater determinants — small enough for dense diagonalisation.

Build and diagonalise the full FCI matrix manually

import numpy as np
from scipy.linalg import eigh
from vibeqc.solvers import (
    build_hamiltonian_matrix,
    build_hamiltonian_mo,
    generate_determinants,
    get_hf_orbital_provider,
)
from vibeqc.solvers._determinant import excitation_rank
from vibeqc.solvers._slater_condon import (
    diagonal_matrix_element_unrestricted,
    single_excitation_matrix_element,
    double_excitation_matrix_element,
)

# Assumes `mol`, `basis`, `H` from the selected-CI example above
norb = H.norb
nalpha = nbeta = mol.n_electrons() // 2

# Generate all determinants (441 for H₂O/STO-3G)
all_dets = generate_determinants(norb, nalpha, nbeta)
ndet = len(all_dets)
print(f"Full FCI space: {ndet} determinants")

# Build the Hamiltonian matrix
H_full = np.zeros((ndet, ndet))

# Diagonal elements
for I, (aI, bI) in enumerate(all_dets):
    H_full[I, I] = diagonal_matrix_element_unrestricted(aI, bI, H.h1e, H.h2e)

# Off-diagonal: only connect determinants with excitation rank ≤ 2
for I in range(ndet):
    aI, bI = all_dets[I]
    for J in range(I + 1, ndet):
        aJ, bJ = all_dets[J]
        rank_a = excitation_rank(aI, aJ)
        rank_b = excitation_rank(bI, bJ)
        total = rank_a + rank_b
        if total > 2:
            continue

        val = 0.0
        if total == 1:
            if rank_a == 1:
                holes = sorted(set(aI) - set(aJ))
                parts = sorted(set(aJ) - set(aI))
                val = single_excitation_matrix_element(aI, holes[0], parts[0], H.h1e, H.h2e)
            else:
                holes = sorted(set(bI) - set(bJ))
                parts = sorted(set(bJ) - set(bI))
                val = single_excitation_matrix_element(bI, holes[0], parts[0], H.h1e, H.h2e)
        elif total == 2:
            if rank_a == 2:
                holes = sorted(set(aI) - set(aJ))
                parts = sorted(set(aJ) - set(aI))
                val = double_excitation_matrix_element(aI, holes[0], holes[1], parts[0], parts[1], H.h2e)
            elif rank_b == 2:
                holes = sorted(set(bI) - set(bJ))
                parts = sorted(set(bJ) - set(bI))
                val = double_excitation_matrix_element(bI, holes[0], holes[1], parts[0], parts[1], H.h2e)
            elif rank_a == 1 and rank_b == 1:
                from vibeqc.solvers._slater_condon import _phase_factor
                h_a = sorted(set(aI) - set(aJ)); p_a = sorted(set(aJ) - set(aI))
                h_b = sorted(set(bI) - set(bJ)); p_b = sorted(set(bJ) - set(bI))
                sign = _phase_factor(aI, [(h_a[0], p_a[0])]) * _phase_factor(bI, [(h_b[0], p_b[0])])
                val = sign * H.h2e[h_a[0], h_b[0], p_a[0], p_b[0]]

        H_full[I, J] = val
        H_full[J, I] = val

# Diagonalise
evals, evecs = eigh(H_full)
e_fci = evals[0] + H.nuclear_repulsion

print(f"E(FCI, full)  = {e_fci:.8f} Ha")
print(f"E(RHF)        = {rhf_result.energy:.8f} Ha")
print(f"E_corr        = {e_fci - rhf_result.energy:+.8f} Ha")

Literature benchmark — H₂O FCI

The vibe-qc FCI energy for H₂O / STO-3G is −75.01241 Ha, matching the literature reference of −75.01266 Ha to 0.12 mHa (Bauschlicher & Taylor, J. Chem. Phys. 85, 2779, 1986). The sub-milliHartree agreement validates the integral transformation, Slater–Condon rules, and Hamiltonian diagonalisation against an independent implementation — this is the gold-standard test for any post-HF code path.

The corresponding correlation energy is −0.04960 Ha (≈ 31.1 kcal/mol, or about 5% of the total electronic energy). Most of it comes from determinants where the α and β occupation patterns differ — the so-called “open-shell singlet” configurations are essential for quantitative dynamic correlation.

If your FCI energy differs from −75.01241 Ha by more than 1 μHa, check your geometry (the Pitzer/Pulay equilibrium is R(OH) = 1.8089 bohr, ∠HOH = 104.52°) and basis set (STO-3G).

DMRG — bond-dimension scheduling

The Density Matrix Renormalisation Group represents the wavefunction as a matrix product state (MPS) — a chain of tensors, one per orbital, with a controllable bond dimension \(M\) that sets the maximum entanglement captured.

The key idea is bond-dimension scheduling: start cheap (\(M=4\)), then progressively increase \(M\) to refine the wavefunction. Each sweep uses the previous result as a warm start.

from vibeqc.solvers import (
    DMRGOptions,
    build_hamiltonian_mo,
    get_hf_orbital_provider,
    solve_dmrg,
)

# Assumes `mol`, `basis` from above
C = get_hf_orbital_provider(mol, basis)
H = build_hamiltonian_mo(mol, basis, C)

opts = DMRGOptions(
    bond_dim_schedule=[4, 8, 16, 32, 64],
    n_sweeps=8,
    conv_tol_energy=1e-6,
    verbose=1,
)
result = solve_dmrg(H, opts)

print(f"E(DMRG)     = {result.energy:.10f} Ha")
print(f"bond_dim    = {result.bond_dim}")
print(f"trunc_err   = {result.truncation_error}")

For H₂ / STO-3G (2 orbitals, 2 electrons), DMRG recovers the FCI energy at \(M=4\). For larger active spaces (12–16 orbitals), you’ll need \(M\) in the hundreds to approach chemical accuracy.

The Python DMRG backend is a teaching tool. It uses NumPy-level tensor contractions without quantum-number (particle, spin) symmetries. For production DMRG on more than ~12 orbitals, use an external DMRG engine like block2 as the backend — the Hamiltonian object is designed to interoperate.

How DMRG relates to Selected-CI

  • Selected-CI is best when the important determinants are sparse (a few dominate). It’s an additive approach — you explicitly list determinants.

  • DMRG is best when the entanglement is spread across many orbitals (strong correlation, bond breaking). It’s a compressive approach — the MPS implicitly encodes an exponential number of determinants.

  • For systems with both sparse and entangled correlation (most real molecules at equilibrium), the two give similar accuracy at similar cost. Pick the one you’re more comfortable tuning.

v2RDM — variational 2-RDM under N-representability

Instead of working with the \(N\)-electron wavefunction, v2RDM works directly with the two-electron reduced density matrix (2-RDM):

\[E[{}^1D, {}^2D] = \sum_{ij} h_{ij} \, {}^1D_{ij} + \frac{1}{4} \sum_{ijkl} g_{ijkl} \, {}^2D_{ijkl} + E_{\rm nuc}\]

The 1-RDM is obtained by contraction: \({}^1D_{ij} = \frac{1}{N-1} \sum_k {}^2D_{ikjk}\).

The catch: not every 2-RDM comes from a valid \(N\)-electron wavefunction. The N-representability problem is solved by enforcing positivity constraints:

Constraint

Ensures

\({}^2D \succeq 0\)

Two-particle RDM is PSD

\({}^2Q \succeq 0\)

Two-hole RDM is PSD

\({}^2G \succeq 0\)

Particle-hole RDM is PSD

These are the P, Q, and G conditions. Together they give energies that agree with FCI to within ~0.1 mHa for small molecules.

from vibeqc.solvers import (
    V2RDMOptions,
    build_hamiltonian_mo,
    get_hf_orbital_provider,
    solve_v2rdm,
)

C = get_hf_orbital_provider(mol, basis)
H = build_hamiltonian_mo(mol, basis, C)

opts = V2RDMOptions(
    constraints="pqg",       # full P+Q+G
    mu=10.0,                 # initial penalty parameter
    mu_factor=1.2,           # penalty growth factor
    outer_max_iter=500,
    conv_tol_primal=1e-6,
    verbose=1,
)
result = solve_v2rdm(H, opts)

print(f"E(v2RDM)           = {result.energy:.8f} Ha")
print(f"constraint_residual = {result.constraint_residual:.2e}")

You can trade accuracy for speed by relaxing constraints:

  • constraints="p" — only \({}^2D \succeq 0\) (fast, approximate)

  • constraints="pq" — add the Q condition (most of the correlation)

  • constraints="pqg" — full P+Q+G (close to FCI accuracy)

v2RDM uses \(O(n^6)\) storage (the 2-RDM is a 4-index tensor). For more than ~10–12 orbitals, use an active space. A warning is emitted automatically.

Transcorrelated — similarity-transformed Hamiltonian

Transcorrelated methods take a different approach. Instead of expanding the wavefunction, they apply a similarity transformation to the Hamiltonian:

\[\tilde{H} = e^{-\tau} \, H \, e^{\tau}\]

where \(\tau\) is a correlation factor (a Jastrow factor, typically depending on inter-electronic distances \(r_{ij}\)). The transformed Hamiltonian \(\tilde{H}\) has:

  • Screened two-electron integrals — the electron repulsion is damped at short range, removing the Coulomb cusp that’s hard to capture with a determinant expansion.

  • Three-body terms — from commutators \([\tau, [T, \tau]]\). These are folded into lower-rank contributions via the normal- ordered two-body approximation (NO2B).

  • Non-Hermiticity — \(\tilde{H}\) is generally non-Hermitian, but you can symmetrise it for standard eigensolvers (at the cost of some theoretical rigour).

After the transformation, you solve \(\tilde{H}\) with any two-body solver — typically Selected-CI, since the transformed Hamiltonian’s determinant expansion converges much faster than the bare one.

run_job one-liner

import vibeqc as vq

mol = vq.Molecule([
    vq.Atom(1, [0.0, 0.0, 0.0]),
    vq.Atom(1, [0.0, 0.0, 1.4]),
])

vq.run_job(
    mol,
    basis="sto-3g",
    method="transcorrelated_ci",
    selected_ci_options=vq.SelectedCIOptions(target_size=50, verbose=1),
    transcorrelated_options=vq.TranscorrelatedOptions(
        form="simple_gaussian",
        gamma=0.5,
        no2b=True,
        symmetrize=True,
    ),
    output="h2_tc",
)

Manual API

from vibeqc.solvers import (
    TranscorrelatedOptions,
    build_hamiltonian_mo,
    build_transcorrelated_hamiltonian,
    get_hf_orbital_provider,
    solve_selected_ci,
    SelectedCIOptions,
)

C = get_hf_orbital_provider(mol, basis)
H_bare = build_hamiltonian_mo(mol, basis, C)

tc_opts = TranscorrelatedOptions(
    form="simple_gaussian",
    gamma=0.5,
    no2b=True,
    symmetrize=False,       # keep non-Hermitian for formal rigour
    report_diagnostics=True,
)
H_tc = build_transcorrelated_hamiltonian(H_bare, tc_opts)

# Solve the transcorrelated Hamiltonian with Selected-CI
ci_opts = SelectedCIOptions(target_size=100, verbose=1)
result = solve_selected_ci(H_tc, ci_opts)

print(f"E(TC-SelCI) = {result.energy:.8f} Ha")

The gamma parameter controls correlation strength — larger values fold more short-range correlation into the Hamiltonian, making the CI converged faster but increasing the three-body contribution (and the NO2B approximation error). Typical values are in the range 0.1–1.0.

Active spaces

For systems where a full-orbital treatment is too expensive (DMRG with many virtuals, or v2RDM on more than ~12 orbitals), you can restrict the calculation to an active space: a subset of orbitals and electrons that captures the dominant correlation effects.

Through run_job:

vq.run_job(
    mol,
    basis="cc-pVDZ",
    method="selected_ci",
    active_space=(6, 4),   # 6 active orbitals, 4 active electrons
    selected_ci_options=vq.SelectedCIOptions(target_size=500, verbose=1),
    output="h2o_active",
)

The tuple (n_active, n_elec) means:

  • n_active — number of spatial orbitals in the active space (taken from around the HOMO–LUMO gap).

  • n_elec — number of electrons in the active space.

In the low-level API, active-space selection is a separate step:

# After building the full MO Hamiltonian:
from vibeqc.solvers import canonical_orthogonalize
import numpy as np

# The Hamiltonian.active_space() method truncates to the
# top n_active orbitals with n_elec electrons.
# (Currently exposed via run_job; the low-level API is in
#  active development — see the examples/wavefunction/
#  directory for the current pattern.)

Choosing an active space

For most organic molecules near equilibrium:

  1. Count valence orbitals: for each first-row atom, include the 2s and 2p (4 orbitals per heavy atom). Add the 1s for hydrogens.

  2. Count valence electrons: all electrons beyond the core (1s² for C/N/O, none for H).

  3. Freeze the core: the 1s orbitals of heavy atoms are chemically inert and can be frozen at the HF level.

  4. Include a few low-lying virtuals for dynamic correlation.

For H₂O / cc-pVDZ: 4 occupied orbitals (O 2s, O 2p×3 + H 1s×2, but one is the bonding combination), 2 core (O 1s), leaving ~5–6 valence orbitals with 6–8 electrons as a natural active space.

Non-HF orbitals — canonical orthogonalisation

You don’t have to use HF orbitals. Any orthonormal set of orbitals spanning the AO basis is valid — the solvers only require that the Hamiltonian integrals are in an orthonormal basis.

The simplest alternative is canonical (Löwdin) orthogonalisation:

\[X = U s^{-1/2} U^T\]

where \(S = U s U^T\) is the eigen-decomposition of the AO overlap matrix. This produces orbitals closest (in a least-squares sense) to the original AOs while enforcing orthonormality.

from vibeqc.solvers import (
    Hamiltonian,
    SelectedCIOptions,
    build_hamiltonian_ao,
    canonical_orthogonalize,
    solve_selected_ci,
    transform_hamiltonian,
)
import vibeqc as vq
import numpy as np

mol = vq.Molecule([
    vq.Atom(8, [ 0.0,  0.00,  0.00]),
    vq.Atom(1, [ 0.0,  1.43, -0.98]),
    vq.Atom(1, [ 0.0, -1.43, -0.98]),
])
basis = vq.BasisSet(mol, "sto-3g")

# 1. Build AO Hamiltonian
H_ao = build_hamiltonian_ao(mol, basis)

# 2. Get AO overlap and orthogonalise
S = np.asarray(vq.compute_overlap(basis))
X = canonical_orthogonalize(S)

# 3. Transform to orthonormal basis (non-HF!)
H_orth = transform_hamiltonian(H_ao, X)

# 4. Solve
result = solve_selected_ci(H_orth, SelectedCIOptions(target_size=100, verbose=1))
print(f"E(Sel-CI, Löwdin) = {result.energy:.8f} Ha")

Canonical orbitals are:

  • Guaranteed to exist — no SCF convergence problems.

  • Fragment-consistent — the orbitals are local and resemble AOs (useful for embedding / fragment methods).

  • Less compact than HF — the determinant expansion is larger because the reference determinant is farther from the exact wavefunction. You’ll need a larger target_size to reach the same accuracy.

This is also how you’d use orbitals from any external program: grab the MO coefficients (as an (n_ao, n_mo) matrix), feed them to transform_hamiltonian, and solve. The solver doesn’t care where the orbitals came from — it just needs orthonormality.

Where to go next

  • Examples directory — runnable scripts for each method live in examples/wavefunction/, including H₂ dissociation curves, H₄ linear chain, OH radical, and CH₂ singlet-triplet splitting.

  • API reference — all solver classes and functions are documented in api/index under vibeqc.solvers.

  • SCF convergence — if your HF orbitals aren’t converging, check SCF convergence before feeding them to a post-HF solver.

  • Regression tests — compare your results against the benchmark values in tests/.

  • Roadmap — C++ DMRG with SU(2) symmetry and GPU-accelerated FCI are on the roadmap.