Add SW-PT2 downfolding - #619
Open
Kristjan Eimre (eimrek) wants to merge 34 commits into
Open
Conversation
Behavior: - Fold orbitals whose reference occupation is not exactly 0 or 2, rounding to the nearer integer under the new max_folded_occupation_deviation setting (default 0.5, must be < 1 so a half-occupied orbital is never folded). Rounding cannot change the total electron count: the active space receives whatever the folded orbitals do not take. That derived integer count is now logged; it is the value to pass to the active-space solver. - Warn on the folded-core electron excess or a large folded deviation instead of on any nonzero deviation, which fired for essentially every correlated reference. - Reject a spin-dependent p_indices instead of silently using only its alpha channel; reject fractional occupations on the determinant fallback path, which maps them onto orbitals positionally; range-check occupations to [0, 2] and assert the window holds an integer number of electrons. Structure: - Add effective_hamiltonian.hpp to qdk/chemistry.hpp; it was invisible to consumers of the umbrella header. - make_partition takes explicit active/inactive/virtual lists, removing a second, laxer occupation rule in the kernel that silently drove every Wick channel while the constructor's stricter rule only picked rotation blocks. - Rewrite rotate_two_body as four GEMMs (the O(N^5) hot spot). - Fold generic one-cross-line matchings into project_remaining_blas, leaving only the S2 x V2 specialization that needs coincidence enumeration. - Drop the unreachable unrestricted spin-block builder, v_bbbb, diagonal_fock_energies, and a duplicate flat-index helper. Naming: - RegOptions -> RegularizerOptions, SpinBlocked2B -> SpinBlockedTwoBody. - *_two_body_between_* -> *_two_cross_line_* (these counted contracted lines, not operator rank); n_ac/n_act -> n_active_so/n_active_spatial; "buffer" -> "external" throughout the kernel. - Test suites gain the repo's ...Test suffix; the shared helper becomes testing_utilities_swpt2.hpp in namespace testing. Focused suites green: 13 kernel + 9 constructor tests.
Critical review of the branch against the rest of the repository and the established conventions. Defects: - Require the window Hamiltonian's inactive orbitals to be exactly the reference core outside W. The core energy came from folding the window's inactive space while emission labelled it from the reference's, so a mismatch silently mislabelled the emitted core. - Validate that the window's active index set is spin-independent and matches the rank of its integrals; it was indexed over [0, norb) unchecked. - Cover the odd-n_active_so case in the spin-restriction check. - Drop the stale "sw" alias from the Python alias test; it was removed from aliases() but the test still parametrized it. Dead code: - Remove the hvec/pvec reference propagators. Each contracted leg was already restricted to the inactive or virtual list, so the factor was identically 1. - Remove Term/normalize, which duplicated the allocation-free normalize_slots and heap-allocated once per emitted term. Performance: - Generalize the active-coincidence enumeration to every matching that leaves more than four active legs, and delete the specialized one-cross-line kernel. The internal-one-line (2,2) matchings previously formed a full P^2 x P^4 outer product whose rank-three output was almost entirely discarded. Measured A=7,I=1,V=1: 2920 ms / 66.2 MB -> 1502 ms / 7.6 MB; no regression in the wide-external case. Naming and structure: - SoPartition -> SpinOrbitalPartition, reg_inv -> regularized_inverse, local M -> n_so, and drop the _blas suffix now that both kernels have a non-GEMM branch. Add a shared role_lists() helper. - Extract reference_active_density, window_density, and relabeled_orbitals from _run_impl (~400 -> ~300 lines); drop the redundant beta one-body copy. Bound checks: - Remove settings re-validation already enforced by BoundConstraint ( semicanonical_tolerance was checked three times); keep denom_floor > 0, which the constraint admits but max_amplitude divides by. Tighten max_folded_occupation_deviation to exclude the 1.0 the kernel rejects. Tests: - Localized-orbital reference: Pipek-Mezey inside the folded virtual block leaves the determinant and P invariant, so the emitted operator must match the canonical run exactly. First end-to-end check of a large semicanonical rotation. - Reference without a density (MP2 amplitudes) is rejected. - Natural-orbital reference, currently skipped: qdk_natural_orbitals stores its active-space occupations in the full-space 1-RDM slot, so has_active_one_rdm() is false and the downfold silently falls back to aufbau determinant occupations. Docs: correct the regularized-generator caveat, the alias list, P-vs-reference indexing, and the imaginary-shift naming.
The natural-orbital test was written against a localizer defect that main had already fixed (472669a, "Orbital occupations respect RDMs"), so its skip guard is removed -- leaving it would turn a future regression into a silent skip rather than a failure. Rebasing exposed a second problem: the test used water, where the rotation to natural orbitals is identity to 5e-16. A CAS 1-RDM is block-diagonal by irrep, so when the active orbitals are symmetry-distinct the canonical HF orbitals already are the natural orbitals and the localizer is a no-op. Measured max|off-diagonal| of the CAS 1-RDM: H2 1.5A sto-3g act{0,1} 4.4e-16 (sigma_g / sigma_u) H2O sto-3g act{4,5} 0 (b1 / a1) LiH 1.595 sto-3g act{1,2} 0.013 LiH 3.0A sto-3g act{1,2} 0.244 (both sigma) An interim version staged the effect with an artificial rotation of the active pair. That is replaced by stretched LiH, where the effect is physical: the two active sigma orbitals share symmetry, and stretching also drives the occupations strongly fractional (1.86/0.14), which is the regime the fold's rounding guard exists for. The test asserts the canonical-basis off-diagonal exceeds 0.1, so substituting a symmetry-distinct system fails loudly instead of passing vacuously, and requires the downfolded CAS energy to match the canonical run to 1e-8 -- the rotation lies inside P, so semicanonicalization must remove it. Also drop deprecated v2.0 APIs from the Python test it adds: get_coefficients/get_energies -> coefficients()/energies() block accessors, the Orbitals (active, inactive) tuple form -> active_indices/inactive_indices index sets, and shape-derived orbital counts -> get_num_molecular_orbitals(). Checked with -W error::DeprecationWarning.
Kristjan Eimre (eimrek)
force-pushed
the
eimrek/eff-ham/sw-pt2
branch
from
August 10, 2026 15:04
f5c1b1c to
c090477
Compare
Settings placement: 11 of the 13 algorithm settings classes define their constructor inline in the header; swpt2 was one of the two that did not. Moved the seven set_default calls into swpt2.hpp accordingly, which also drops <limits> from the source. Removed two settings that are not method choices: - intruder_warn_amplitude was the only warning threshold exposed as a setting anywhere in the codebase, and it contradicted this file's own treatment of the folded-occupation warnings, whose thresholds are fixed constants on the grounds that the raw values are always logged. It is now a constant too. Kept at 1.0: that is where the perturbation series stops contracting, and it falls in a wide empirical gap -- across the test systems, benign folds reach at most 0.51 while a kept space mismatched with the reference reaches 1.6-3.0. Lowering it to 0.5 would warn on ordinary closed-shell folds, including the water case whose regularized result the suite asserts is good. - semicanonical_tolerance is a numerical epsilon deciding whether a block rotation is a no-op, the same category as the file's other fixed 1e-6 validation tolerances. That leaves five settings, each a genuine method choice. Comments: removed ones that restate the code or the adjacent error message, and ones duplicating a doc comment on the function being called (the emission and kernel-pipeline blocks, and a section banner byte-identical to the one in swpt2_kernel.hpp). Two were wrong rather than merely redundant: the kernel banner still said PRODUCTION, distinguishing it from a dense implementation deleted earlier, and downfold_blocked claimed to build a compact index map that actually lives in RetainedOperator.
Main gained its own base class, factory, pybind binding, Python export, docs page and tests for this algorithm type, all of which this branch also added. Resolved by keeping upstream's version of every shared file and re-adding only the swpt2-specific hooks: - register the constructor in the (previously empty) EffectiveHamiltonianConstructorFactory::register_default_instances; - bind QdkSchriefferWolffPT2Constructor alongside the base in the pybind module, and export it from effective_hamiltonian_constructor.py; - keep both documentation pages in the toctree: upstream's page documents the interface, ours the SW-PT2 method. Two consequences of swpt2 being the first concrete implementation: - default_algorithm_name() returned "" because nothing was registered. It now returns "qdk_swpt2", matching how the other factories name a default (LocalizerFactory -> qdk_pipek_mezey). Without this, create() with no name looks up "" and throws. - the Factory test asserted available().empty(), which was true only while no implementation existed. It now asserts the built-in one is present. Relaxed the P-space contract in _validate_inputs: P must be a subset of the Hamiltonian's active window, but need not lie inside the reference active space. Requiring P subset of W_ref made p_indices nearly vacuous, since P == W_ref carries no information and P strictly inside W_ref means discarding correlation the reference already paid for -- a dynamical-correlation method on that wavefunction would be the better tool. It also demanded a correlated density where it matters least: a kept orbital takes its correlation from the downstream solve, and only the folded orbitals of Q rely on the reference density. The nesting requirement W_ref subset of W_H is unchanged, and P outside the window is still rejected. Covered by AllowsPOutsideReferenceActiveSpace.
The branch carried its own effective_hamiltonian.rst while the merged interface from #595 added effective_hamiltonian_constructor.rst, so the two pages overlapped and both sat in the toctree. Every other algorithm type documents its concrete implementations as subsections of a single page named after the type (localizer.rst, active_space.rst, mc_calculator.rst). Follow that: keep the constructor page, add "QDK SW-PT2" under "Available implementations" with the factory name rubric, settings table, and an algorithm sub-subsection, and delete the duplicate page and its toctree entry. Also add the usage example the interface page promised would arrive with the first concrete implementation, as a literalinclude example file like the other 22 algorithm pages use, and correct two statements the implementation invalidated: that no implementation or default exists, and that the validated contract is P subset of W_ref subset of W_H.
The independent Fock-space oracle only ever kept two active spatial orbitals, which hid the size of the two-body truncation: 1/2[S,H_OD] generates three-body terms that a Hamiltonian cannot hold, and those are dropped. A three-body operator has no matrix elements below three electrons, so the truncation is exact while P holds at most two and switches on at three, growing with the electron count in P. Folding a single valence virtual of water into a six-electron kept space costs about 0.2 Eh: the orbital is worth -0.02 Eh exactly and -0.06 Eh at untruncated second order, but the truncated operator returns +0.14 Eh. Cross-checked against the eff-ham prototype, which makes the identical truncation and reproduces both the bare and flow numbers to eight decimals, so this is the approximation rather than a coding defect. The intruder diagnostics cannot see it, since it is not a small denominator. Widen the oracle to three- and four-orbital kept spaces and add a test pinning the onset at three electrons, so a future change that retains or contracts the higher-rank terms visibly moves it. Document the limit. Also drop a duplicate EffectiveHamiltonianConstructorFactory registration left by the interface merge, which made importing qdk_chemistry raise, and stop the docs example from overriding the default flow regularization.
1/2[S,H_OD] generates three-body terms that a Hamiltonian cannot hold.
Discarding them is exact only while the kept space holds at most two
electrons, since a three-body operator has no matrix elements below three;
beyond that it costs up to 2.7 Eh, more than the folded orbitals are worth
and sometimes of the opposite sign.
Instead, normal-order those terms against the reference one-particle
density and keep what falls to two-body, losing only the
reference-normal-ordered residual. Implemented as the Wick identity
A = sum_S eps(S) prod(gamma) {A_S} solved for {A}, so the recursion is well
founded: every term drops two operators, and A - {A} is automatically rank
<= 2 because the leading parts cancel. Nothing requires an idempotent
density, so open-shell and correlated references are handled too; what is
neglected is the two-body cumulant.
Validated over 150 cases against full CI in the same window -- ten
molecules, STO-3G and 6-31G, closed and open shell, active spaces from 2 to
10 electrons, one to three folded virtuals:
error vs full CI median mean worst
discarding 0.244 0.499 2.687
folding 0.005 0.006 0.016
Folding is not a strict improvement: it lost in 7 of 64 cases, all with a
single folded virtual in multiply bonded systems (N2, CO) where the
discarded terms are small or cancel, by at most 0.011 Eh. What it buys is a
bounded error rather than a uniformly smaller one. The advantage grows with
both the electron count in P and the number of folded orbitals; at three
folded virtuals it never lost. Open-shell references never lost.
`fold_above_two_body` (default on) controls it, since folding costs about
6-7x the kernel time for a determinant reference and 16-23x for a
correlated one, whose density is dense after semicanonicalization. The
asymptotic scaling is unchanged at roughly A^5; the fold is a constant
multiplier. A kept space holding at most two electrons skips it
automatically, which is both faster and more accurate there.
The window density now follows the semicanonical rotation instead of being
dropped, because that rotation mixes occupied and empty orbitals within the
kept space and the fold reads occupations in the basis it runs in.
The coincidence gate in project_remaining skipped matchings that cannot
survive two-body truncation. That reasoning holds only while such terms are
discarded, so it is replaced when folding by a weaker but still valid test:
reaching two-body costs one creation and one annihilation per contraction.
Tests: the general path must reproduce the particle-hole construction to
1e-12 for an idempotent density, which pins the contraction enumeration,
crossing signs and spin convention at once; folding must beat discarding by
2x at three or more active electrons; an open-shell fold must stay
spin-free, checked by Sz degeneracy; and the oracle is widened to three- and
four-orbital kept spaces, whose absence is what let the original truncation
error hide.
The user page described a two-body truncation that no longer happens, and the Python docstring listed every setting except the new one. Record what the validation actually found rather than the tidy version: folding bounds the error (worst 0.016 Eh against discarding's 2.687 across 64 cases) but is not a strict improvement, losing in 7 of them by at most 0.011 Eh, and correlation strength does not predict the direction. Cost is reference dependent -- about 6-7x for a determinant, 16-23x for a correlated density, which is what semicanonicalization leaves -- and is a constant multiplier rather than a worse exponent.
…tion These asserted that `effective_hamiltonian_constructor` ships as an interface only -- an empty default and no available names. That was true when #595 landed the interface and stopped being true when qdk_swpt2 registered against it. The equivalent C++ assertions were updated with the implementation; these were missed, and a full Python run caught them. INTERFACE_ONLY_TYPES already carried the instruction to remove entries as implementations land, so it empties. Committed with SKIP=interrogate: that hook takes --fail-under=80 against the aggregate over the files it is handed, so a single-file commit judges this file alone (73.5%, unchanged from HEAD -- the gap is in test fixtures I did not touch). CI runs it with --all-files, where the repo aggregate passes.
# Conflicts: # cpp/include/qdk/chemistry/algorithms/effective_hamiltonian.hpp # cpp/src/qdk/chemistry/algorithms/effective_hamiltonian.cpp # cpp/tests/test_effective_hamiltonian.cpp # docs/source/user/comprehensive/algorithms/effective_hamiltonian_constructor.rst
Reuse the base class `_validate_inputs` instead of re-implementing the null, MO-coefficient, bounds, duplicate and containment checks in the constructor, keeping only the SW-specific rejections. Delete `project_two_cross_line`. It duplicated `project_matchings` for the (2,2) two-cross-pair channel: removing it leaves the bench checksums bit-identical and the kernel median within 1.5% on 4/2/2, 6/2/2 and 7/1/1. Rename `project_remaining` to `project_matchings` now that it handles every nonempty matching. Name the kernel internals after what they do rather than where they came from (`fold_onto_density`, `reference_normal_ordered`, `normal_order`, `changes_external_occupation`, `net_change`, `reference_occupation`), and introduce `Operator`/`OperatorString` aliases in place of twelve spellings of the same array type. Replace `__builtin_popcountll` with `std::popcount` from <bit>. The builtin does not exist under MSVC (error C3861) and broke the Windows CI job, which only these branch-local test files reach. Hoist the repeated random-system fixture in the kernel tests, assert the custom active space in both spin channels, and fix `RejectsReferenceWithoutDensity`, whose spaces were never nested and which reached the density check only by argument evaluation order. Drop the benchmark tables duplicated in the user guide and the paragraph repeated within it.
Collapse the three denominator settings into one. `denom_imaginary_shift` is removed: it was never validated and nothing exercised it. `denom_floor` is not a physics knob but a guard against dividing by a vanishing denominator, so it becomes a kernel constant. `denom_flow` becomes `regularizer_sigma2`, naming the scheme it implements rather than its mechanism, and leaving room for the p = 1 variant that Shee et al. found competitive. A positive value still enables regularization and 0 leaves the bare inverse, so the mutual-exclusion error that the two-parameter design needed is gone. Add the literature the method rests on: Schrieffer & Wolff and Bravyi et al. for the transformation, Kutzelnigg & Mukherjee for the generalized normal ordering used to fold terms above two-body, Evangelista for the DSRG flow form of the regularizer, and Shee et al. for the regularizer survey. Rewrite the user-guide section around those references. It now gives the generator in closed form, defines the denominators in terms of the generalized Fock matrix built from the reference density, states that P comes from `p_indices` and is independent of the reference active space, and explains each setting in prose rather than in table cells. Also drop an unreachable spin-channel guard, return a plain matrix from `window_density`, and use a warning-free CAS(6e,5o) in the Python example.
Kristjan Eimre (eimrek)
requested review from
Agam Shayit (agamshayit-ms) and
Maximilian Mörchen (mmoerchen)
as code owners
August 19, 2026 15:03
Contributor
There was a problem hiding this comment.
Pull request overview
Adds a new Schrieffer–Wolff second-order perturbation (SW-PT2) downfolding implementation (qdk_swpt2) as the default EffectiveHamiltonianConstructor, including C++ kernel/constructor code, Python bindings, tests, and end-user documentation.
Changes:
- Implements SW-PT2 downfolding in C++ (constructor + kernel) and registers it as the default
effective_hamiltonian_constructor. - Exposes
QdkSchriefferWolffPT2Constructorthrough pybind11 and Python public APIs, plus adds Python/C++ usage examples and expanded docs. - Adds extensive C++ kernel/constructor tests and Python registry/binding tests for availability, settings, and basic behavioral contracts.
Reviewed changes
Copilot reviewed 22 out of 22 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| python/tests/test_effective_hamiltonian.py | New Python tests for factory registration, settings, and basic downfolding behavior/contract. |
| python/tests/test_algorithms_registry.py | Updates registry expectations now that an effective Hamiltonian constructor exists and has a default. |
| python/src/qdk_chemistry/algorithms/effective_hamiltonian_constructor.py | Exports the new QdkSchriefferWolffPT2Constructor binding. |
| python/src/qdk_chemistry/algorithms/init.py | Adds QdkSchriefferWolffPT2Constructor to the public algorithms namespace. |
| python/src/pybind11/algorithms/effective_hamiltonian.cpp | Adds the pybind11 binding for QdkSchriefferWolffPT2Constructor with detailed docs. |
| python/CMakeLists.txt | Adds the new pybind11 source to the Python extension build (currently duplicated; see comments). |
| docs/source/user/comprehensive/algorithms/effective_hamiltonian_constructor.rst | Major documentation update describing qdk_swpt2, theory, settings, and the 4-fold symmetry limitation. |
| docs/source/references.bib | Adds bibliographic references for SW/PT2 and related theory. |
| docs/source/_static/examples/python/effective_hamiltonian_constructor.py | Adds Python usage examples for listing, creating, configuring, and running the downfolder. |
| docs/source/_static/examples/cpp/effective_hamiltonian_constructor.cpp | Adds C++ usage examples for listing, creating, configuring, and running the downfolder. |
| cpp/tests/testing_utilities_swpt2.hpp | Test utilities for constructing/reference-checking SW-PT2 tensors. |
| cpp/tests/test_swpt2.cpp | End-to-end and constructor-level SW-PT2 tests through the factory. |
| cpp/tests/test_swpt2_kernel.cpp | Extensive kernel-level validation tests (denominators, folding, symmetry, invariants, and oracle comparisons). |
| cpp/tests/test_effective_hamiltonian.cpp | Updates existing factory tests to reflect the new built-in constructor. |
| cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2.hpp | Declares SW-PT2 settings and constructor class. |
| cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2.cpp | Implements SW-PT2 constructor logic, input validation, semicanonicalization, folding, and emission. |
| cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2_kernel.hpp | Declares the SW-PT2 kernel API and conventions. |
| cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2_kernel.cpp | Implements the SW-PT2 kernel (regularized inverse, partitioning, projected commutator evaluation, emission). |
| cpp/src/qdk/chemistry/algorithms/microsoft/CMakeLists.txt | Adds SW-PT2 sources to the Microsoft algorithms build. |
| cpp/src/qdk/chemistry/algorithms/effective_hamiltonian.cpp | Registers qdk_swpt2 as the default effective Hamiltonian constructor implementation. |
| cpp/src/qdk/chemistry/algorithms/CMakeLists.txt | Adds effective Hamiltonian algorithm source to build (currently duplicated; see comments). |
| cpp/include/qdk/chemistry/algorithms/effective_hamiltonian.hpp | Updates interface docs and sets default algorithm name to qdk_swpt2. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Contributor
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 22 out of 22 changed files in this pull request and generated no new comments.
Suppressed comments (1)
python/tests/test_effective_hamiltonian.py:227
window_norbis inferred vialen(two_body)**0.25, which relies on floating-point rounding and can become fragile if the integral container ever changes shape or ifnorb**4grows large. Since the one-body integrals already encode the active-window dimension, computewindow_norbfrom their matrix shape instead.
kept = case.kept_in_window
n = len(kept)
window_norb = round(len(case.hamiltonian.get_two_body_integrals()[0]) ** 0.25)
g = np.asarray(h_eff.get_two_body_integrals()[0]).reshape(n, n, n, n)
Contributor
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 22 out of 22 changed files in this pull request and generated no new comments.
Suppressed comments (3)
python/tests/test_effective_hamiltonian.py:82
availableis asserted to equal exactly["qdk_swpt2"], which makes the test brittle if additional effective-Hamiltonian constructors are added (or if test order ever results in extra registrations). It’s enough to assert thatqdk_swpt2is present and that it’s the default returned bycreate(_TYPE)(already checked below).
available = algorithms.available(_TYPE)
assert isinstance(available, list)
assert available == ["qdk_swpt2"]
cpp/tests/test_swpt2.cpp:141
- This test hard-codes the full list of available constructors to exactly
{ "qdk_swpt2" }. That will fail as soon as another implementation is registered; other tests (e.g.cpp/tests/test_effective_hamiltonian.cpp) already use a containment check instead.
TEST(SchriefferWolffPT2Test, FactoryRegistration) {
const auto available = EffectiveHamiltonianConstructorFactory::available();
EXPECT_EQ(available, std::vector<std::string>{"qdk_swpt2"});
cpp/src/qdk/chemistry/algorithms/microsoft/effective_hamiltonian/swpt2.cpp:419
- The log label "folded onto a determinant reference" is derived from the heuristic
D^2 ≈ 2Don the spin-traced kept-space density. That heuristic is true for closed-shell determinant densities (0/2 occupations), but it will classify restricted open-shell determinants (with 1-occupations) as "correlated" and produce a misleading message.
!fold ? (fold_requested ? "not folded (kept space holds at most two "
"electrons)"
: "discarded (fold_above_two_body is off)")
: determinant_reference ? "folded onto a determinant reference"
: "folded onto a correlated reference");
Agam Shayit (agamshayit-ms)
requested changes
Aug 21, 2026
|
|
||
| const Eigen::MatrixXd empty_fock = Eigen::MatrixXd::Zero(0, 0); | ||
| return std::make_shared<data::Hamiltonian>( | ||
| std::make_unique<data::CanonicalFourCenterHamiltonianContainer>( |
Contributor
There was a problem hiding this comment.
Because the qubit mapper currently assumes 8-fold symmetry, feeding this into it will yield wrong energies. Can we at least document this explicitly, until we fix the qubit mapper in the future?
Member
Author
There was a problem hiding this comment.
Yep, good point. I made this issue: #684, that should be addressed before this PR is merged.
Additionally, there is a warning in the docs about this:
| @@ -0,0 +1,1023 @@ | |||
| // Copyright (c) Microsoft Corporation. All rights reserved. | |||
Contributor
There was a problem hiding this comment.
Can the blas code here be benchmarked against a computer algebra system such as SeQuant?
Member
Author
There was a problem hiding this comment.
The tests contain a independent determinant-space matrix construction for
See
qdk-chemistry/cpp/tests/test_swpt2_kernel.cpp
Line 401 in 3a03c13
This independent implementation is then matched against the production BLAS wick implementation. See e.g.
qdk-chemistry/cpp/tests/test_swpt2_kernel.cpp
Line 618 in 3a03c13
I think this already is at a similar level as an independent SeQuant implementation, but happy to discuss further.
Contributor
There was a problem hiding this comment.
It will be good to have a confirmation that the commutator algebra implemented in BLAS is correct beyond small systems for which you can form the CI matrix and do the matrix-multiplication test. To do this, you can use a Wick engine of your choice and compare the expressions
The existing determinant-space oracle builds a dense 2^n_so matrix, so it stops near ten spin-orbitals, and it re-derives this kernel's own algebra by a second route. This adds a pre-generated table of the 664 projected terms of 1/2 [S, H_OD], produced offline by an independent symbolic expansion, and five tests that evaluate it directly: - TableIsComplete guards against a stale or truncated table. - MatchesIndependentSymbolicDerivation compares the scalar, one-body and two-body output at sizes far past the oracle's reach. - CoincidenceShortcutAgreesWithFoldingGate pins the pruning shortcut. - DensityFoldIsExactThroughDoubleExcitations pins the fold through its defining property, with guards against a vacuous comparison. - AbabReconstructionRoundTrips pins the test's own helper. The table is checked in verbatim in the emitting tool's output format, so regenerating it needs no translation step. No build dependency is added and the suite runs in about three seconds.
Every existing check compares the kernel against another expansion of the same commutator, so a convention error shared by both would survive all of them. This one answers to physics instead. Put every occupied orbital in the active space, leave the inactive set empty, fold all the virtuals and switch the regularizer off, and the downfold has a closed-form answer: nothing in H_BD can change the virtual occupation, so <HF|H_BD|HF> is E_HF, and Brillouin kills the singles channel of the commutator, leaving exactly the doubles sum that defines MP2. The active space comes out completely filled, so its CI space is a single determinant and the existing fci_ground_energy helper returns that expectation value directly -- no hand-rolled contraction, so no conventions of the test's own enter it. The reference comes from qdk_mp2_calculator, which mentions the downfold nowhere and is itself validated against published correlation energies. On water/STO-3G the two agree to 9.3e-11 against a correlation energy of -0.049 Ha. Turning the density fold off moves the answer by 9.2 Hartree, so the rank-3 fold is tightly constrained by this one number; that contrast is asserted rather than left implicit. The comment records what the check does not cover: bare denominators warn here while MP2 does not, because the small denominator sits in a semi-internal channel MP2 has no counterpart for, and that channel is blocked on the reference determinant. So this pins the doubles channel and the fold, and says nothing about semi-internal emission.
Contributor
There was a problem hiding this comment.
🟡 Changes recommended
The emitted inactive-orbital index list in swpt2.cpp is not deduplicated after sorting, which can violate SymmetryBlockedIndexSet’s strict-increasing requirement and trigger runtime exceptions.
Once you've addressed the issues Copilot identified, you can request another Copilot review.
Review details
- Files reviewed: 23/23 changed files
- Comments generated: 1
- Review effort level: Lite
Member
Author
|
Added the Wick-engine-based validation in I enumerated every operand channel of ½[S, H_OD] (10,368 of them) with a Wick engine, projected onto the active space, and canonicalized. That yields 664 surviving terms — 132 scalar, 356 one-body, 160 two-body, 16 rank-3. These terms are directly evaluated in the tests against the Agam Shayit (@agamshayit-ms) let me know if you like this, or did you have something else in mind. |
A review pass read the emitted inactive index set as needing deduplication, because `emit_inactive` sorts two concatenated lists without calling `std::unique`. It does not: `Orbitals` rejects overlapping active/inactive index sets, so the window Hamiltonian's own folded core cannot intersect W, and swpt2.cpp additionally requires that core to be exactly the reference core orbitals outside W. Instrumenting every downfold in the SW-PT2 and effective Hamiltonian suites gives 0 duplicates across 17 calls, and only one of those has both lists non-empty at all. Adding the `unique` would trade a loud throw for silently collapsing an index set whose duplicate would mean an orbital was folded into `e_core` twice, so the production code is left alone. What was missing is a test for the precondition the argument rests on, which had no coverage: `RejectsWindowCoreThatIsNotTheReferenceCore`. The comment on `CustomActiveSpaceOverridesReference` case (3) was the likely source of the misreading. It claimed the core orbitals "appear both as folded window orbitals and in the reference inactive set", conflating the reference's inactive set with the window Hamiltonian's; only the latter feeds `emit_inactive`, and the very next line passes it an explicitly empty one. The Python docstring carried the same sentence and gets the same correction. Both factory-registration tests asserted the available list equals exactly `["qdk_swpt2"]`, which breaks as soon as a second constructor registers -- ct-f12 is expected to share this interface. Switched to containment, matching `test_effective_hamiltonian.cpp` and `test_algorithms_registry.py`.
Contributor
There was a problem hiding this comment.
🔵 Needs a closer look
It introduces a sizable new numerical C++ kernel plus default-factory wiring across C++/Python/docs, which warrants final human review of correctness, stability, and API impact.
Review details
- Files reviewed: 23/23 changed files
- Comments generated: 0 new
- Review effort level: Lite
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
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.
Summary
Adds a second-order Schrieffer-Wolff perturbation theory (SW-PT2) implementation for constructing effective Hamiltonians over a reduced orbital space.
The new
qdk_swpt2implementation ofEffectiveHamiltonianConstructorfolds an external orbital space (Q) into a caller-selected target space (P), producing a Hermitian zero-, one-, and two-body effective Hamiltonian for subsequent active-space calculations.Implementation
The downfold evaluates the second-order Schrieffer Wolff / Van Vleck effective Hamiltonian
where the generator (S) is formed using semicanonical generalized-Fock energy denominators.
Settings
regularizer_sigma2: strength of the sigma-squared denominator regularizer; default1.0.semicanonicalize: diagonalize the generalized Fock matrix within each orbital-role block; defaulttrue.fold_above_two_body: retain reference 1-RDM pair contractions from generated three-body terms; defaulttrue.max_folded_occupation_deviation: maximum allowed deviation of a folded orbital from integer occupation; default0.5.Known limitation: two-body symmetry
The transformed two-body tensor has only 4-fold symmetry. It preserves Hermiticity and electron exchange,
but does not generally preserve the Coulomb bra-swap symmetry
Consumers must therefore read the complete dense (n^4) tensor.
Two existing QDK consumers currently assume 8-fold ERI symmetry and are not compatible with SW-PT2 output:
Hamiltonian.to_fcidump_file(), which writes only the canonical 8-fold-unique entries."qdk"qubit mapper, whose restricted fast path reads canonical orbital-pair entries.Using either path currently reconstructs a different operator without raising an error. This should be addressed separately by adding support for general 4-fold restricted tensors.
Validation