CASSCF groundwork: revive the Python bindings, add the shared C++ CAS engine and both geometry adapters

[susilehtola]

Groundwork for CASSCF. Started as a Python-bindings repair and has grown into the geometry-agnostic C++ CAS integral engine plus both adapters. Draft because the CI solver it feeds depends on reci, which is not public yet — everything here is the orbital side, which does not.

The branch name is now too narrow for its contents; happy to split the C++ CAS commits (e44e155, 43a3445, a4bfd99) onto their own branch if you would rather review them apart from the Python work.

The bindings were dead code

python/bindings.cpp has not compiled since the Armadillo→Eigen migration — it includes <armadillo> and libhelfem/include/ArmaEigen.h, neither of which still exists. HELFEM_PYTHON defaults OFF, so nothing ever tried. Five accessors had also been renamed underneath it: get_basis→make_basis, get_nbf→nbf, get_lval→lval, get_mval→mval, get_rad_Nel→rad_Nel.

The port removes more than it adds (−46/+43 in bindings.cpp): the SCF surface returns helfem::Matrix natively now, so the Armadillo round-trip at the numpy boundary was pure overhead.

helfem.pyscf_driver was missing

python/test_radial_df_factors.py imports install_full_eri, helfem_scf and build_active_eri from helfem.pyscf_driver — a module that is not in the tree. Both DF tests failed at their import line. This writes it to that contract.

The module rests on one identity. HelFEM exposes no four-index tensor, but coulomb contracts the second pair, so feeding it a pair density and contracting the free pair yields MO integrals directly:

(tu|vw) = C_v^T . coulomb(C_t C_u^T + h.c.) . C_w / 2

One coulomb call per pair, no AO→MO transform. With C = I that extracts the full AO tensor (tests and small bases); restricted to an active space it is precisely the integral list a CASSCF needs — which is the point, and why no new integral code is required for it. pyscf imports sit inside the functions that use them, so the tests stay runnable on a machine without pyscf.

Tests

All three Python tests now pass, having previously been unable to import:

test	result
test_smoke	H 1s eigenvalue −0.5000000000 against analytic
test_radial_df_factors	DF factors vs reassembled AO ERI, max err 1.7e-21; FE sparsity 16.3% nonzero, 17.3× compression
test_radial_df_exhaustive	intra max err 2.2e-16, inter 2.8e-17

Registered with CTest under the label python, 4 s total: ctest -L python. Total registration goes 201 → 204; nothing else moved.

They are copied into the build tree rather than run where they sit. Python puts a script's own directory first on sys.path, so running them in place resolves import helfem to the source package — which holds no extension module — and fails however PYTHONPATH is set.

A latent CMake bug, fixed first

enable_testing() sat below the build subdirectories. CMake silently discards add_test() calls made in a directory reached before testing was enabled — no warning, no error, the directory's CTestTestfile.cmake simply is not generated. Any test registered from libhelfem/, libhelfemqc/, python/ or src/ would have vanished without trace. Nothing registered tests there yet, so nothing was lost; registering the Python tests is how it surfaced. Moved to its own commit since it is independent of the rest.

Measured

On He (Nbf=14, lmax=0, nelem=3, nnodes=6, Rmax=20):

check	result
J, K rebuilt from the extracted tensor vs coulomb()/exchange()	2e-16 / 5e-16
8-fold permutational symmetry	2.2e-16
pyscf RHF on these integrals vs HelFEM's own atomic binary	−2.859592542390 vs −2.8595925424
patched-get_jk vs explicit-tensor SCF route	9.3e-15
build_active_eri vs pyscf ao2mo	1.1e-16
CAS(2,1) == RHF	1.0e-14
ninact=1/nact=0 core energy == RHF	8.9e-16
energy re-contracted from the CI RDMs; Tr(D)	0.0e+00; 2.000000000000

The last row is the identity a CI solver's RDM contract has to satisfy, so it is confirmed on HelFEM's side before any solver is wired in.

Full ctest quick tier at OMP_NUM_THREADS=1: see the comment below.

The shared C++ CAS engine

src/general/cas_integrals.{h,cpp} carries the whole algebra once, geometry-agnostically. A geometry supplies exactly three things — coulomb, exchange, hcore — through cas::JKProvider, and everything else (inactive Fock, active ERIs, generalized Fock, orbital gradient, Fock response, the κκ Hessian block) is built on top of the same pair-density identity the Python driver uses. So an adapter is ~60 lines and no new integral code exists in either geometry.

src/general/cas_integrals_test.cpp pins it against the Python/PySCF reference at deterministic core-guess orbitals, agreeing to ~1e-13 across E_inactive, h_eff, the ERI list, the generalized Fock, the gradient norm and the κκ Hessian.

One convention needed reconciling: HelFEM's exchange() returns the signed contribution added to a spin channel's Fock matrix, where PySCF wants a positive K subtracted from a total density — a sign and a factor of two, which is why cas::closed_shell_veff is coulomb(P) + exchange(P/2) rather than the textbook expression.

The diatomic adapter, and the flag it refuses

DiatomicCASEngine (src/diatomic/cas_engine.h) throws if the basis has absm_symmetric set. With that flag, exchange() builds only the m ≥ 0 half of K and mirrors the rest — exact for a density symmetric under m → −m, wrong otherwise. A CAS forms pair densities C_t C_uᵀ between individual active orbitals, and one between the +m and −m partners of a π shell carries Δm = 2.

src/diatomic/cas_diatomic_test.cpp builds exactly that density and measures the damage rather than asserting it:

check	result
adapter refuses an absm_symmetric basis	throws
E_inactive == 0 at ninact=0	exact
h_eff == CᵀhC at ninact=0	exact
8-fold permutational symmetry of (tu|vw)	1.1e-16
π⁺/π⁻ pair density: |K_general − K_mirrored|max	6.93e-04, equal to |K_general|max
the same comparison on an m-symmetric density	0.0

The mirror overwrites the blocks carrying the whole result, with nothing to flag it at runtime — so the refusal is necessary, not defensive.

Scope, because the obvious reading is too strong: this is a statement about exchange(), not about --symmetry. coulomb() never consults the flag, it defaults to false, and the AO basis is identical for every --symmetry value — only the orthonormalization blocking and that flag differ. The CAS integrals are symmetry-agnostic, and a --symmetry=3 reference is usable input provided the flag is cleared. What --symmetry=3 cannot carry is the active space, since an |m| > 0 block there is one set of radial coefficients standing for two spatial orbitals; that is a separate question from the integrals.

An unrelated diatomic fix, surfaced by that test (d8f6af8)

converge_block's quadrature-cap warning computed its relative change as (cur - prev) after prev had been overwritten with cur three lines above, so it printed relative change still 0.000e+00 unconditionally — useless since it was written.

Fixing that immediately exposed a second thing: the warning fires on diatomic-H2-hf-r, a CI case, reporting 1.191e-01. Real but harmless. radial_integral(-1,0) is ∫BᵢBⱼ/sinh(μ)dμ, which diverges logarithmically at μ=0 — the block magnitude grows with quadrature order (5.57 → 6.96 → 7.90) while every other block reaches 1e-15 at n=160, and the growth is exactly logarithmic (the 320→512 step is a factor 1.6, and 1.39·log(512/320)/log(2) = 0.94, the measured difference). The divergence is confined to the row and column of the one basis function non-zero at μ=0 — which is exactly the function Nbf() drops from every m ≠ 0 shell, the only consumer. It is discarded before the matrix is seen, which is why it is dropped.

So converge_block gained an nskip parameter excluding those entries from the convergence test (the block is still returned whole). H2 is unchanged at −1.1336295702, the HF limit to all printed digits, and the warning no longer fires.

This commit touches only basis.cpp and is independent of the CAS work — it cherry-picks onto master on its own.

🤖 Generated with Claude Code

[susilehtola]

Full ctest quick tier at OMP_NUM_THREADS=1, -j4:

100% tests passed, 0 tests failed out of 61
Total Test time (real) = 1223.02 sec

61 = 8 unit + 49 quick integration + integration.invariants + the 3 new python cases. The ***Not Run (Disabled) lines in the log are the slow tier, registered-but-disabled as designed — not failures.

No C++ outside python/ changed and bindings.cpp links into no binary, so this was confirming rather than discovering; recorded for the record.

[susilehtola]

Added: Phase 2, the CASSCF orbital gradient (139f0e2)

Scope extended beyond the original "revive the bindings" framing — say the word and I'll split this into its own PR.

g = 2(F − Fᵀ) with the generalized Fock built from the active-space RDMs:

F[m][i] = 2 ( F^I[m][i] + F^A[m][i] )            i inactive
F[m][t] = sum_u D_tu F^I[m][u] + Q[m][t]         t active
Q[m][t] = sum_uvw d_tuvw (mu|vw)

Every term is a coulomb()/exchange() call on a pair density — Q costs one coulomb() per active (t,u) pair, contracting the 2-RDM slice sum_vw d_tuvw C_v C_w^T. Still no AO→MO transform anywhere, so this ports to the diatomic basis unchanged.

Verified by central-differencing frozen_ci_energy, which rebuilds the energy through active_hamiltonian and shares no algebra with the gradient formula. On Be, CAS(2,3) over a 1s core:

check	worst |fd − analytic|
inactive–active	1.9e-09
inactive–virtual	1.1e-08
active–virtual	1.5e-09
active–active	1.2e-09
unit-norm direction, Richardson to h→0	1.5e-10
CI-relaxed vs frozen-CI gradient	2.5e-11
CAS(4,14) vs FCI in the same space	0.0e+00

Redundant classes are checked too and must match: at frozen CI an active–active rotation does change the energy, and only becomes redundant once the CI relaxes. The CI-relaxed row is that statement quantitatively.

One thing worth recording: normalizing the probe direction matters. The central-difference truncation error scales as ‖K‖³h², so an unnormalized random direction over all 91 pairs disagreed by 6e-2 relative — for entirely uninteresting reasons. A step-size study gave the clean 4.00×-per-halving signature of O(h²), which is how that was settled rather than assumed.

The new test needs pyscf and scipy, unlike the numpy-only DF tests, so CMake probes for them separately. Python cases 3 → 4; total registration 204 → 205; ctest -L python 4.5 s, all passing.

[susilehtola]

Added: Phase 3, two-step CASSCF (95f77d9)

The plan's load-bearing verification. PySCF's mcscf.CASSCF run on HelFEM's own AO integrals, against a native optimizer on the same integrals — two independent optimizers, nothing shared but the input. On Be CAS(2,3) over a 1s core (Nbf=14):

HelFEM two-step CASSCF	−14.563849375547 (71 macro, 106 CI solves, |g| 2.3e-07)
PySCF mcscf.CASSCF	−14.563849375707
difference	1.6e-10

Re-referencing, learned the hard way

dE/dx equals the generalized-Fock gradient only at x = 0 — away from it the Fréchet derivative of the matrix exponential enters. My first attempt handed scipy L-BFGS the gradient from the rotated point under a frozen-reference parametrization, an inconsistent objective/gradient pair. It thrashed (‖g‖ oscillating between 1e-3 and 2e-2) and crawled to −14.56062 against a true minimum of −14.56385.

The gradient was not at fault, and I checked that rather than assuming it: re-verified at a random point well away from the reference (‖K‖ = 2.1), it reproduced finite differences to 3.7e-07 on values of order 10–50. Each macroiteration now works at x = 0 and folds the accepted step into C, with a backtracking line search on the true energy — monotone by construction.

A second exact redundancy

Worth flagging for the second-order path. Beyond the active-active rotations already excluded, this solution has an active natural orbital with occupation 0.00000000: all its 1- and 2-RDM elements vanish, its column of the generalized Fock is identically zero, and every rotation involving it is exactly flat.

That is why the two codes agree on the energy to 1.6e-10 while their occupied spaces differ in one direction by 3.9°. Not a convergence failure — they stopped at different points along a genuinely flat direction.

So the test asserts on the energy and on the core density (which agree to 3.7e-06) and reports the natural occupations, rather than demanding the orbitals match. A CASSCF code comparison that checks orbitals will produce false failures.

The consequence for Phase 4: this is a real null space of the orbital Hessian that has to be projected out, and HelFEM's conditioning report — which treats exact zeros in the Hessian diagonal as a bug — will flag it. Correct for an SCF; expected for a CAS.

Python cases 4 → 5, total registration 205 → 206, ctest -L python 9.5 s all passing.

[susilehtola]

Added: Phase 4, fully coupled (κ, c) CASSCF (df29ff3)

All four Hessian blocks analytic, built from the same coulomb()/exchange() pair-density primitive as the integrals and the gradient — so still no AO→MO transform anywhere, and the construction carries to the diatomic basis unchanged.

block	how
H_κκ	d/dt of the generalized Fock under a one-index transformation
H_ss	2 P⊥(H−E)P⊥ — a σ-vector
H_κs	the generalized Fock evaluated on transition RDMs
H_sκ	= H_κsᵀ

On Be CAS(2,3) over a 1s core (Nbf=14):

	E	macro	final ‖g‖
two-step	−14.563849375547	71	2.3e-07
coupled	−14.563849376675	22	2.0e-08

Closing iterations go ‖g‖ = 1.4e-04 → 7.5e-08 → 2.0e-08 — the quadratic behaviour a correct second-order method should show. The two agree to 1.1e-09.

Three things that had to be got right rather than assumed

The transition case is not the state case with different numbers substituted. For two states of the same CAS, ⟨I|Ê_ij|0⟩ = δ_ij⟨I|0⟩ = 0 over the inactive block, so the inactive-Fock contribution to the inactive columns drops out while the active columns keep it. generalized_fock gains core_occ: 2 for a state density, 0 for a transition density.

BCH ordering. Differentiating the gradient walks exp(tK')exp(sK) while the Hessian is defined on exp(sK+tK'), so the two differ by ½ g·[K',K] and agree only where the gradient vanishes. Rather than dodge this by testing at a stationary point — where it cannot appear — the test predicts the discrepancy, and the antisymmetric part of the raw derivative reproduces the commutator to five significant figures.

A spurious ½ halved the κκ block. The scalar form contracts the full antisymmetric matrix and double-counts each pair; extracting pair entries directly does not. It did not present as a factor error — the optimizer merely stalled 8.8e-04 high, which looks like a hundred other things. Splitting the Hessian check by segment identified it immediately: ratio exactly 2.000000 on the orbital block against 1.000000 on the CI block.

The null space must be projected out, not shifted

Worth recording, because it will recur in the C++ solver. A CAS Hessian is genuinely singular beyond the active-active rotations already excluded: an active orbital with zero natural occupation makes every rotation touching it exactly flat. Here 10 of 51 modes have |eigenvalue| < 1e-6, carrying a gradient of 1.7e-09 — i.e. nothing.

A Levenberg shift does not handle that. It leaves those modes with |g/w| up to 8e-02 from dividing ~0 by ~0; the trust radius then clips the whole step and throttles the informative directions. The optimizer stalls at ‖g‖ ~ 1e-02 while looking merely slow. Dropping the modes costs nothing, since the energy does not depend on them.

Python cases 5 → 6, total registration 206 → 207, ctest -L python 23 s all passing.

[susilehtola]

Added: the C++ CAS integral engine (e44e155)

The first C++ piece, and the part that needs no CI solver — so it proceeds while the reci interface request is outstanding.

src/general/cas_integrals.{h,cpp} holds the active-space Hamiltonian, the generalized Fock matrix and the orbital gradient. src/atomic/cas_engine.h is the entire atomic side: 60 lines.

That split is the payoff from the Python reference. None of this needs an AO→MO transform, because coulomb() contracts the second index pair — pair density in, contract the free pair, MO integrals out. So a geometry only has to supply J, K and hcore (cas::JKProvider, three virtuals) and the algebra lives once. The diatomic adapter will be the same size.

Validation

Against python/helfem/pyscf_driver.py, which is itself checked against PySCF (active integrals to 1e-16 vs ao2mo; CASSCF to 1.6e-10 vs mcscf.CASSCF on the same integrals). The two share no code, only the underlying J/K builds, so agreement is a real check rather than a tautology.

Orbitals come from the core guess — generalized eigenvectors of (hcore, S), deterministic in both languages, with an explicit phase convention since an eigenvector is only defined up to sign. Everything matches to ~1e-13, the print precision of the reference:

E_inactive	−13.487633038924
trace(h_eff)	−0.498155424803
sum(eri)	2.508862565558
sum(F)	−10.092279000769
norm(g)	5.799542776404

The exchange convention

The thing most likely to break, and worth stating. HelFEM's exchange() returns the signed contribution added to a spin channel's Fock matrix (src/atomic/main.cpp:383, Exx = 0.5 tr(P_spin K_spin)), while PySCF wants a positive K subtracted from a total density. The two differ by a sign and by a factor of two. cas::closed_shell_veff is the single place reconciling them, so nothing else has to think about it — and E_inactive above is what proves it right.

core_occ

generalized_fock takes it because a transition density is not a state density with different numbers substituted: for two states of the same CAS, ⟨I|Ê_ij|0⟩ = δ_ij⟨I|0⟩ = 0 over the inactive block, so the inactive-Fock term drops from the inactive columns while the active columns keep it. The test asserts exactly that — core_occ=0 must leave the active columns alone and must change the inactive ones.

Unit tests 8 → 9, total registration 207 → 208; cas_integrals_test runs in 0.04 s. Full unit + python labels: 15/15 passing.

[susilehtola]

Added: the C++ κκ Hessian block (43a3445)

fock_response_kappa is d/dt of the generalized Fock under a one-index transformation (dC = C dκ). Every piece of the Fock is a coulomb/exchange call on a density built from C, so its derivative is the same calls on the transformed densities plus contraction terms — still no new integral machinery.

Two ordering facts, both silent if wrong

hess_kappa_kappa_raw is not the Hessian. Differentiating the gradient walks exp(tK')exp(sK) while the Hessian is defined on exp(sK+tK'); by BCH the two differ by ½ g·[K',K] and agree only where the gradient vanishes. hess_kappa_kappa symmetrises, removing it exactly. The test doesn't dodge this by working at a stationary point where it can't appear — it predicts the discrepancy, and the raw derivative's antisymmetric part reproduces the commutator to 2.1e-15.

The contraction carries two separate halvings: one for the symmetrisation, one because contracting full antisymmetric matrices counts every rotation pair twice. Dropping the second halves the block — which doesn't present as a factor error. In the Python reference it surfaced only as an optimizer stalling 8.8e-04 high.

The test now finite-differences itself, with no CI solver

frozen_ci_energy is E_inactive + Σ h_eff·D + ½ Σ eri·d, so a fixed RDM pair suffices. It's rebuilt through active_hamiltonian and shares no algebra with the gradient or the Hessian — which is what makes differencing it a check rather than a re-expansion.

check
vs Python reference (sum(dF), norm(hess_kk_raw), elements)	~5e-13
K1·gradient, Richardson	9.1e-11
K1·H_κκ·K2, four-point	7.7e-08
BCH asymmetry vs ½g·[K',K]	2.1e-15

One note on method. A single central difference of the gradient is truncation-limited at 3.0e-07 for h=2e-4 and says nothing about the formula — my first tolerance flagged it as a failure. Rather than loosen the tolerance I checked the signature: the error falls by exactly 4.00× per halving, the O(h²) fingerprint. The test now asserts the Richardson-extrapolated value and that fourfold fall, so a real error can't hide behind a loose bound later.

unit label: 9/9 passing.

[susilehtola]

Split into two independent PRs, as discussed:

• C++ CAS integral engine with atomic and diatomic adapters #352 — the C++ CAS integral engine, its atomic and diatomic adapters, and the κκ Hessian block, including the complex-orbital fix (active_eri, generalized_fock and fock_response_kappa no longer symmetrise transition densities).
• Python bindings revived, with a PySCF-driven CASSCF for the atomic basis #353 — the revived Python bindings and the PySCF-driven CASSCF. It now refuses complex (m ≠ 0) orbitals with ComplexOrbitalError instead of silently symmetrising them, since the real 8-fold tensor PySCF takes cannot represent them.

No commit touched both halves, so the split is clean. d8f6af8 is dropped from both: #349 already put the quadrature-warning fix on master.