Repository navigation
Python bindings revived, with a PySCF-driven CASSCF for the atomic basis - #353
Open
susilehtola wants to merge 5 commits into
Open
susilehtola wants to merge 5 commits into
susilehtola wants to merge 5 commits into
Conversation
python/bindings.cpp had not compiled since the Armadillo-to-Eigen
migration: it included <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: the SCF surface is natively
helfem::Matrix now, so the arma round-trip at the numpy boundary was
pure overhead.
python/helfem/pyscf_driver.py was missing entirely -- two of the three
existing tests import install_full_eri, helfem_scf and build_active_eri
from it, so they could not even get past their import line. It is
written here 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 gives MO
integrals directly,
(tu|vw) = C_v^T . coulomb(C_t C_u^T + h.c.) . C_w / 2
at one coulomb() call per pair and no AO->MO transform. With C = I that
extracts the full AO tensor; restricted to an active space it is the
integral list a CASSCF needs. pyscf imports are kept inside the
functions that use them so the tests stay runnable without pyscf.
Measured on He (Nbf=14, lmax=0, nelem=3, nnodes=6, Rmax=20):
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 0.0e+00, Tr(D) = 2
The three tests are now registered with CTest under the label "python"
(4 s total). 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.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The CAS energy gradient with respect to an orbital rotation, at fixed CI
coefficients, is g = 2(F - F^T) with the generalized Fock matrix
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 costing one
coulomb() per active (t,u) pair with the 2-RDM slice sum_vw d_tuvw C_v C_w^T
as its argument. As with the active integrals, no AO->MO transform appears
anywhere: the 1- and 2-RDMs plus the existing contractions are the whole
input, which is what makes this portable to the diatomic basis unchanged.
Verified by central-differencing frozen_ci_energy, which rebuilds the CAS
energy through active_hamiltonian and so shares no algebra with the gradient
formula. That separation is deliberate: CLAUDE.md records that a derivative
check differentiating one expression against finite differences of itself
agrees perfectly with a consistently wrong one.
On Be, CAS(2,3) over a 1s core (Nbf=14):
inactive-active worst |fd - analytic| 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 is allowed to relax. Check 3 is that statement quantitatively --
at a converged CASCI the CI is variational, so the frozen-CI gradient already
equals the relaxed one.
Normalizing the probe direction matters and is worth recording: the
central-difference truncation error scales as ||K||^3 h^2, so an unnormalized
random direction over all 91 pairs disagrees by 6e-2 relative for entirely
uninteresting reasons. A step-size study showed the clean 4.00x-per-halving
signature of O(h^2), which is how that was settled rather than guessed at.
The new test needs pyscf and scipy, unlike the numpy-only DF tests, so
CMakeLists probes for them separately rather than failing on a numpy-only
machine. Python cases: 3 -> 4, total registration 204 -> 205.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…grals Solve the CI at fixed orbitals, build the orbital gradient from the resulting RDMs, step the orbitals, repeat. This is the reference implementation the second-order path will be diffed against, and the first point in the arc that produces a working CASSCF. Verification is the tightest available here: PySCF's mcscf.CASSCF running on HelFEM's own AO integrals against this optimizer on the same integrals. Two independent optimizers, identical input, nothing shared but the integrals -- so unlike a derivative check it cannot be fooled by a consistently wrong energy expression. On Be CAS(2,3) over a 1s core (Nbf=14): HelFEM -14.563849375547 (71 macro, 106 CI solves, |g| 2.3e-07) PySCF -14.563849375707 diff 1.6e-10 Re-referencing is the whole design and was learned the hard way. dE/dx equals the generalized-Fock gradient ONLY at x = 0; away from it the Frechet derivative of the matrix exponential enters. A 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: re-checking it at a random point well away from the reference (||K|| = 2.1) 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, so the energy falls monotonically. An additional exact redundancy shows up here, beyond the active-active rotations already excluded: the converged solution has an active natural orbital with occupation 0.00000000, so 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 degrees. The test therefore 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. The second-order solver will meet this as a genuine null space of the orbital Hessian -- and HelFEM's conditioning report, which treats exact zeros in the Hessian diagonal as a bug, will flag it as one. Python cases: 4 -> 5, total registration 205 -> 206. ctest -L python, 9.5 s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Orbital rotations and CI coefficients optimized together rather than alternating. All four Hessian blocks are analytic and built from the same coulomb()/exchange() pair-density primitive as the integrals and the gradient, so no AO->MO transform appears anywhere and the construction carries to the diatomic basis unchanged: H_kk d/dt of the generalized Fock under a one-index transformation H_ss 2 P_perp (H - E) P_perp -- a sigma-vector H_ks the generalized Fock evaluated on TRANSITION RDMs H_sk = H_ks^T On Be CAS(2,3) over a 1s core (Nbf=14), against the two-step reference: two-step -14.563849375547 71 macro, |g| 2.3e-07 coupled -14.563849376675 22 macro, |g| 2.0e-08 closing |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 had to be got right rather than assumed, each caught by a check that could fail. The transition case is NOT the state case with different numbers substituted. For two states of the same CAS, <I|E_ij|0> = delta_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. Differentiating the gradient walks exp(tK')exp(sK) while the Hessian is defined on exp(sK + tK'), so by BCH the two differ by (1/2) g . [K', K] and agree only where the gradient vanishes. Rather than dodge this by testing at a stationary point, the test predicts the discrepancy: the antisymmetric part of the raw derivative reproduces the commutator to five significant figures. A spurious 1/2 in the kappa-kappa assembly halved that block -- the scalar form contracts the full antisymmetric matrix and double-counts each pair, extracting pair entries directly does not. It did not look like a factor error: the optimizer merely stalled 8.8e-04 high. Splitting the Hessian check by segment identified it immediately as ratio exactly 2.000000 on the orbital block against 1.000000 on the CI block. Step control needs the null space PROJECTED OUT, not shifted. 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, and 10 of the 51 modes here have |eigenvalue| < 1e-6 carrying a gradient of 1.7e-09. A Levenberg shift 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, stalling at |g| ~ 1e-02 while looking merely slow. Python cases: 5 -> 6, total registration 206 -> 207. ctest -L python, 23 s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
build_active_eri passed coulomb() the SYMMETRISED pair density C_t C_u^T + C_u C_t^T and halved the result, storing ((tu|vw) + (ut|vw))/2. That is exact for real orbitals only. HelFEM's atomic basis uses complex spherical harmonics, so any orbital with m != 0 is complex, and for it the antisymmetric part of the pair density -- a purely imaginary density -- is discarded. The C++ CAS engine was fixed to keep that channel; this route cannot, because it stores the real 8-fold tensor PySCF's ao2mo and FCI take. The damage was not confined to CAS: install_full_eri calls build_active_eri with C = I, so for lmax > 0 the whole AO tensor handed to PySCF was symmetrised, and even its RHF exchange would have been wrong. Every test here ran at lmax = 0, where nothing can tell. So it now REFUSES such pairs with ComplexOrbitalError. The check is exact and costs nothing: with a real tensor coulomb(P^T) = coulomb(P)^T, so J = coulomb(C_t C_u^T) is symmetric iff the pair density is a real function. 0.5 (J + J^T) equals the old symmetrised value by linearity, so accepted numbers are unchanged to roundoff (5.6e-17 on Be lmax=1 -- not to the bit, since coulomb() is evaluated on a different argument). The other symmetrising sites (generalized Fock and its response) contract only active pairs, and both CASSCF drivers pass every accepted orbital set through build_active_eri first -- checked in the code, and stated as exactly that at each site rather than as a blanket guarantee. test_real_orbital_guard, on an lmax = 1 basis: a pure m = +1 orbital is refused; pure m = 0 orbitals in the SAME basis are accepted, so the guard judges the orbitals rather than the basis; and the accepted tensor matches the old construction to 1e-14 relative. Also restores enable_testing() before add_subdirectory(python). Master calls it after, so python/'s tests were silently dropped -- configure reported them registered and `ctest -L python` found none. #344 carried this fix; it was nearly lost in the split, because the check that it was "already on master" was run against a checkout of the #344 branch. Also: the bindings' includes move to <helfem/X.h>, master's convention since the header namespacing. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Revives HelFEM's Python bindings and adds a PySCF-driven CASSCF on top of them — two-step and fully coupled (κ, c) — for the atomic basis. Split out of #344, which also carried the C++ CAS engine (now #352).
The bindings were dead code
python/bindings.cpphadn't compiled since the Armadillo→Eigen migration: it included<armadillo>and a removedArmaEigen.h, and five accessors had been renamed underneath it.HELFEM_PYTHONdefaults OFF, so nothing noticed. Ported to Eigen (−46/+43), andhelfem.pyscf_driver— imported by two existing tests but never committed — is written to that contract.The Python tests were silently not running
Master calls
enable_testing()afteradd_subdirectory(python), and CMake silently dropsadd_test()calls made before testing is enabled. Configure reported the Python tests registered;ctest -L pythonfound none.enable_testing()now comes before any subdirectory that registers tests, still gated onHELFEM_BUILD_TESTS.CASSCF
Every integral comes from
coulomb()/exchange()on MO pair densities — no four-index AO tensor, no AO→MO transform:Two-step CASSCF re-solves the CI inside each orbital step. The coupled version puts orbital rotations and CI coefficients in one Newton step and converges quadratically: 71 → 22 macroiterations, ‖g‖ 1.4e-4 → 7.5e-8 → 2.0e-8. Two traps found along the way are documented where they apply: the CAS Hessian has two exact null spaces (active-active rotations, and any rotation touching a zero-occupation orbital), which must be projected out rather than Levenberg-shifted.
Real orbitals only — and now refused, not silently wrong
This route stores
(tu|vw)as a real tensor with the 8-fold symmetry PySCF'sao2moand FCI assume. That's exact for real orbitals. HelFEM's atomic basis uses complex spherical harmonics, so an m ≠ 0 orbital is complex, and for it(tu|vw) ≠ (ut|vw).build_active_eriused to symmetrise the pair density and carry on. The damage wasn't confined to CAS:install_full_ericalls it with C = I, so for lmax > 0 the whole AO tensor handed to PySCF was symmetrised and even its RHF exchange would have been wrong. Every test ran at lmax = 0, where nothing can tell.It now raises
ComplexOrbitalErrorinstead. The check costs nothing: with a real tensorcoulomb(Pᵀ) = coulomb(P)ᵀ, soJ = coulomb(C_tC_uᵀ)is symmetric iff the pair density is a real function, and½(J + Jᵀ)is the old value to roundoff. Both CASSCF drivers pass every accepted orbital set through it before reaching the other symmetrising sites, checked in the code and stated at each site. The C++ engine (#352) keeps the imaginary channel; this route can't, because PySCF can't take it.test_real_orbital_guardon an lmax = 1 basis: a pure m = +1 orbital is refused; pure m = 0 orbitals in the same basis are accepted, so the guard judges orbitals, not the basis; and the accepted tensor matches the old construction to 2.2e-17.Tests
test_smoketest_radial_df_factors,_exhaustivetest_cas_gradienttest_casscfmcscf.CASSCFon the same integralstest_casscf_coupledtest_real_orbital_guardctest -L python7/7,ctest -L unit11/11.🤖 Generated with Claude Code