Non-mean-field solvers (vibeqc.solvers)

vibe-qc v0.9.0 — Developer documentation

Overview

The vibeqc.solvers package provides wavefunction-based electronic-structure methods that do not conceptually rely on a Hartree–Fock reference state. All solvers share a common interface:

result = solver.solve(hamiltonian, options) -> SolverResult

where hamiltonian is a Hamiltonian object carrying one- and two-electron integrals in an orthonormal spatial-orbital basis. The orbital source (HF, canonical orthogonalisation, localised orbitals, …) is a separate concern handled by the orbital-provider layer.

Architecture

vibeqc.solvers
├── _common.py          # Hamiltonian, SolverOptions, SolverResult
├── _hamiltonian.py     # AO → MO integral construction
├── _determinant.py     # Determinant / CSF data structures
├── _slater_condon.py   # Slater–Condon matrix-element rules
├── _selected_ci.py     # CIPSI-style selected-CI solver
├── _dmrg.py            # Two-site DMRG / MPS solver
├── _v2rdm.py           # Variational 2-RDM via SDP
└── _transcorrelated.py # Similarity-transformed Hamiltonian

Layers

  1. Integrals / orbital basis / Hamiltonian construction (_hamiltonian.py): Builds Hamiltonian objects from vibe-qc’s native C++ integral bindings. Supports AO → MO transformation and canonical (Löwdin) orthogonalisation. The get_hf_orbital_provider helper is a convenience, not a requirement.

  2. Many-body objects (_determinant.py, _slater_condon.py): Determinant representations as tuples of occupied indices, excitation-rank logic, Slater–Condon matrix elements (diagonal, single, double), and full Hamiltonian-matrix construction.

  3. Solver backends (_selected_ci.py, _dmrg.py, _v2rdm.py): Each solver implements solve(Hamiltonian, Options) → SolverResult.

  4. Diagnostics and benchmarking: The SolverResult object carries method-specific diagnostics (PT2 corrections, truncation errors, constraint residuals, RDM objects) alongside energy traces.

  5. Public API / CLI exposure: vibeqc.solvers.__init__ exports all public classes and convenience functions. The package is importable as from vibeqc import solvers or from vibeqc.solvers import ….

run_job integration

The wavefunction solvers are reachable both as a standalone vibeqc.solvers import and through the high-level run_job entry point:

from vibeqc import Molecule, run_job

mol = Molecule.from_xyz("h2.xyz")
result = run_job(mol, basis="sto-3g", method="fci", output="output-h2-fci")

run_job(method=…) accepts "fci", "selected_ci", "dmrg", "v2rdm", and "transcorrelated_ci" alongside the mean-field "rhf" / "uhf" / "rks" / "uks". The solver-specific option objects are passed via the matching run_job kwargs (selected_ci_options=, dmrg_options=, v2rdm_options=, transcorrelated_options=); an active_space=(n_orbitals, n_electrons) kwarg applies a frozen-core / CAS-style truncation before the solver runs. run_job builds the Hartree-Fock orbital basis internally (UHF when multiplicity > 1, RHF otherwise) and hands the resulting Hamiltonian to the solver.

Where HF still lives

The HF reference is a convenience default, not a structural dependency of the solver layer:

  1. get_hf_orbital_provider in _hamiltonian.py — purely a convenience. Callers can pass any orthonormal orbital basis via transform_hamiltonian(ham_ao, C) or construct the Hamiltonian directly.

  2. The run_job entry point builds an HF orbital basis as the default orbital source before dispatching to a solver — but the solver itself only sees the abstract Hamiltonian. Supply a different orbital set by building the Hamiltonian yourself and calling the solve_* function directly.

  3. The v2RDM initial guess — _init_rdm2 builds the HF 2-RDM from aufbau occupation. This is a starting point, not a structural dependency; the RDM is immediately modified by the solver.

The solver layer itself (_selected_ci.py, _dmrg.py, _v2rdm.py, _transcorrelated.py) has zero HF-specific code. All methods operate on the abstract Hamiltonian interface.

Method summaries

Selected CI (_selected_ci.py)

Implements a CIPSI-style (Configuration Interaction using a Perturbative Selection made Iteratively) algorithm:

  1. Start from a reference determinant.

  2. Generate all singles and doubles from the current space.

  3. Estimate PT2 weights: ΔE ≈ |⟨D|H|Ψ₀⟩|² / (E₀ − ⟨D|H|D⟩).

  4. Add top-N candidates to the variational space.

  5. Diagonalise H in the expanded space.

  6. Repeat until energy change < conv_tol_energy or target_size reached.

Key options: target_size, pt2_threshold, do_pt2_correction, spin_restricted.

Scaling: O(N_det² + N_cand · n_occ² · n_vir²). For small active spaces (≤ 12 orbitals) this runs in seconds.

DMRG (_dmrg.py)

Implements a two-site DMRG sweep with MPS wavefunction:

  • MPO construction via complementary-operator approach.

  • Two-site variational optimisation (power iteration).

  • SVD truncation with discarded-weight reporting.

  • Bond-dimension schedule (e.g. [4, 8, 16, 32, 64]).

Limitations (minimal implementation):

  • Python-level tensor contractions only — adequate for ≤ 12 orbitals and bond dimension ≤ 128. Production use requires C++/GPU.

  • No quantum-number symmetry blocking.

  • The MPO is built as a placeholder; the two-site solver uses power iteration on the reduced wavefunction.

Key options: bond_dim_schedule, n_sweeps, truncation_tol, orbital_order.

Variational 2-RDM (_v2rdm.py)

Minimises E[²D] = Tr(h·¹D) + ¼ Tr(g·²D) under the PSD constraint ²D ≥ 0 and Tr(¹D) = N via an augmented-Lagrangian method:

  1. Compute energy gradient ∂E/∂²D.

  2. Add penalty gradient for Tr(¹D) − N.

  3. Take a gradient step.

  4. Enforce antisymmetry.

  5. Project onto the PSD cone (eigenvalue thresholding).

  6. Update Lagrange multiplier.

Known limitation: Only ²D ≥ 0 is enforced (the P in PQG). The ²Q ≥ 0 and ²G ≥ 0 constraints are not currently implemented, which means the energy is not strictly variational (can be above the exact ground state). A full SDP solver (e.g., SDPA, MOSEK) or the boundary-point method with all three PSD constraints would be needed for quantitative accuracy.

Key options: mu (penalty), mu_factor, conv_tol_primal.

Transcorrelated Hamiltonian (_transcorrelated.py)

Applies a similarity transformation H̃ = e^{−τ} H e^{τ} with a correlator τ:

  • Simple Gaussian: Damp two-electron integrals by (1 − γ·exp(−Δ)) based on orbital-index proximity.

  • Exponential Jastrow: Longer-tailed damping, (1 − γ·(1 − exp(−βΔ))/β).

  • Custom: Pass-through for user-supplied transformations.

  • NO2B: Normal-ordered two-body approximation — discards the 3-body terms that arise from double commutators.

  • Symmetrisation: Optional Hermitisation for standard eigensolvers.

Diagnostics: Reports norms of original/transformed 1e and 2e terms, discarded 3-body norm, and anti-Hermitian part magnitude.

Key options: form, gamma, no2b, symmetrize, report_diagnostics.

Usage example

from vibeqc import Atom, Molecule, BasisSet
from vibeqc.solvers import (
    build_hamiltonian_mo,
    get_hf_orbital_provider,
    SelectedCIOptions,
    solve_selected_ci,
)

# Build the system
mol = Molecule(
    [Atom(1, [0, 0, 0]), Atom(1, [0, 0, 1.4])],
)
basis = BasisSet(mol, "sto-3g")
C = get_hf_orbital_provider(mol, basis)  # any orbitals work
H = build_hamiltonian_mo(mol, basis, C)

# Run Selected-CI
opts = SelectedCIOptions(target_size=50, verbose=1)
result = solve_selected_ci(H, opts)
print(f"E = {result.energy:.10f} Ha  ({len(result.ci_labels)} determinants)")

Testing

.venv/bin/python -m pytest tests/test_solvers_*.py -v

36 tests covering: common interfaces, determinant generation, Slater–Condon rules, Hamiltonian construction, Selected-CI convergence, v2RDM constraint enforcement, and transcorrelated workflows.

References

See individual module docstrings for method-specific references. Key architectural references for each method are listed in the __init__.py module docstring.