Bernoulli Expansion based Unitary Coupled-Cluster - #581
Merged
Conversation
ajay-mk
force-pushed
the
ajay/feat/bernoulli-v2
branch
from
July 27, 2026 12:08
2b73c55 to
e78e554
Compare
ajay-mk
force-pushed
the
ajay/feat/bernoulli-v2
branch
from
August 5, 2026 05:17
65112a0 to
c667779
Compare
ajay-mk
force-pushed
the
ajay/feat/cc-rdm
branch
from
August 6, 2026 18:41
5ce4bfd to
2ea26c0
Compare
Adds the bottom layer of the Bernoulli expansion of the unitary-CC similarity-transformed Hamiltonian: a Wick reduction that retains partial contractions, so a product of normal-ordered operators reduces to a sum of normal-ordered operators rather than collapsing to a scalar vacuum average, and the normal-ordered commutator built on it. WickTheorem::use_topology is disabled explicitly rather than left alone: it defaults to ON (wick.hpp), and its one-representative-times-multiplicity bookkeeping is only exercised by the fully-contracted path. On this partial-contraction path it rescales terms whose amplitude pairs are symmetric, which leaves vacuum averages correct while corrupting projections onto excited manifolds. wick_commutator reindexes B's summed indices to fresh temporaries before forming A*B, since A and B are independently constructed and may otherwise share labels, which would fuse two independent summations.
Adds the second layer: the split of an operator O into O_N, "the non-diagonal part containing all the excitation and de-excitation operators" (defined above Eq. (43) of 10.1063/1.5030344), and the rank-preserving remainder O_R = O - O_N. The Bernoulli expansion's inner commutators carry N/R subscripts, so every nesting level needs this classification. Classification needs definite index spaces, so expand_to_blocks first rewrites each general index of the residual NormalOperator as a sum over the base spaces it spans. Only the hole and particle spaces are expanded over: in the single-reference setting the remaining base spaces are empty, so restricting to those keeps the expansion 2-way per index instead of compounding across the nested commutators. That makes the routine single-reference only, which the header warns about. The rank cutoff mirrors pdaggerq (nt_bra > bernoulli_excitation_level -> R) rather than the paper's uncapped O_N, since that is the convention defining qUCCSD and the one the numbers are validated against; terms above the cutoff fall to R rather than being dropped.
Adds the top layer: hbar(N, rank, skip1) sums H̄⁰..H̄^rank of 10.1063/1.5030344 Eq. (45), each order transcribed from its equation with the published coefficients and per-level N/R subscripts. Bernoulli numbers B₁=-1/2, B₂=1/12, B₃=0, B₄=-1/720 (Eq. 40) enter as those coefficients; a subscript R/N on a commutator means "form the commutator, then keep only its R/N part before the next nesting", which is what the split from the previous commit provides. Two cancellations from the paper are relied on and noted in place: F enters H̄ only at first order (stated just below Eq. (50)), and the higher orders carry only R-subscripted inner commutators. Every term is a nested commutator whose prefix is shared with other terms, within a rank and across ranks, so nest() memoizes each prefix (keyed by the base operator plus the tags applied so far). The nine rank-4 terms have only 3 distinct level-1 and 6 distinct level-2 nodes. Reusing a memoized ExprPtr is safe because expression composition deep-copies its operands. Contributions accumulate through Sum::append rather than chained operator+, which deep-copies the whole accumulated Sum on every call and is quadratic in the term count at high rank.
Adds CC::Options::hbar_expansion (BCH by default, Bernoulli opt-in) and dispatches CC::hbar() to bernoulli::hbar() when it is selected. Two constructor assertions guard the combination: the Bernoulli expansion is defined for the unitary ansatz only, and it requires an explicit hbar_comm_rank, since CC::hbar() otherwise falls back to rank 4 and would silently select the most expensive and least exercised order. CC::energy() takes the plain reference expectation value under this expansion: the tensor-level H̄ is already fully expanded, so no operator connectivity remains to constrain. Its comm_rank argument defaults to the amplitude rank and is passed explicitly for the qUCCSD [2|3] split, where the energy is taken at H̄³ while the amplitudes stop at H̄².
Pins the derived equations at Bernoulli ranks 1-3 for the unitary ansatz: term counts for the energy and for the singles/doubles residuals, plus the guards on invalid configurations (Bernoulli with a non-unitary ansatz, and Bernoulli without an explicit hbar_comm_rank). The rank-3 numbers (46 energy, 32 singles, 38 doubles terms) are the ones cross-checked term-by-term against pdaggerq, so a change here means the derivation changed.
Corrections: - Eq. (45) is the assembly H̄ = Σ_k H̄^k; (46)-(50) are the rank-by-rank operators. The file header attributed (45)-(50) to the latter. - expand_to_blocks: the SR "o"/"g" base spaces are not empty, so the old justification for dropping them was wrong. They are droppable because the single-reference projection annihilates those terms. The header @warning said the same wrong thing. - R_part's result is not block-resolved; it stays in compact general-index form. The header claimed the opposite. - wick_reduce leaves at most one residual NormalOperator, not exactly one -- fully-contracted terms carry none, which find_nop already handled. - The memo shares level-1 nodes across ranks; it is not a prefix relation. - The use_topology rescaling set is {2, 1/2, 1/3, 8/3, 2/3}; the comment said 3 where it should have said 1/3.
… key
A character outside {A,N,R} in a nest() tag string silently read as 'A'
(no filter) and would have yielded the wrong H̄; the whole Eq. (46)-(50)
transcription lives in these strings, so assert on every tag.
Grow the memo key in place instead of deriving it from the memo iterator:
container::map is a flat_map, whose insertions invalidate iterators, so
the read-back was correct only by the accident of no insertion happening
in between.
Include <algorithm> for std::max and range/v3's primitives for
ranges::distance rather than relying on them arriving transitively.
CC::Options::screen and use_topology reach the derivation only through CC::ref_av(); the Bernoulli path calls op::tensor::ref_av() directly and so picks up that function's own defaults instead.
The srcc.cpp analogue for the unitary ansatz, covering both H̄ expansions. CC::t() yields the whole equation set in one derivation -- element 0 the energy, element R the residual -- and the term counts are pinned so a change in either expansion fails ctest. Registered variants run in seconds; the Bernoulli H̄⁴ pins are recorded but left out of ctest, since that configuration takes ~2 minutes against sub-second times for everything else in this directory.
Both accumulation sites in bernoulli.cpp built a Sum by append and left every duplicate for one final simplify. Sum::append flattens nested sums and adds up Constants but never merges like terms, so the nested commutators -- which overlap heavily by construction, the same fact that makes nest() memoize -- carried their duplicates all the way to the end. hbar() now accumulates into a HashingAccumulator, which keys summands by hash under proportional_to and merges them via Product::add_identical at insertion. The prefactor has to be folded into each summand rather than wrapped around the sum: appending Constant*Sum inserts the scaled sum as one opaque summand, since append's flatten splits a Sum but not a Product wrapping one, and nothing would collapse. Distributing it with expand() instead is shorter but materializes an intermediate Sum and gives back most of the gain (rank 4: 163.9 s vs 148.2 s). expand_to_blocks_reduced's outer loop now uses transform_sum_expr, which canonicalizes each mapped result before accumulating -- necessary here because the block assignments carry fresh temporary indices and so cannot hash-collide until canonical. It canonicalizes IN PLACE, hence expand_term now clones rather than returning one of its arguments in the two early-exit paths; that path is also parallel (std::execution::par_unseq), which is safe because Index::next_tmp_index is a static std::atomic. Derivation time, tests/integration/ucc 2 bernoulli <rank>, relwithdebinfo: rank 2 0.416 -> 0.364 s, rank 3 7.584 -> 7.216 s, rank 4 170.6 -> 148.2 s. Output is unchanged: serialize() of the rank-3 and rank-4 equations is byte-identical to the pre-change baseline (45 543 and 292 512 bytes), and repeat runs are byte-identical to each other despite the parallel index minting. Term counts alone would not have been sufficient evidence -- the use_topology bug rescaled terms while leaving counts and the VEV correct -- and to_latex() would not either, since it omits symmetry attributes.
Trim the development narrative out of bernoulli.{cpp,hpp} and keep the
reasons the code needs. The use_topology(false) comment keeps its cause
(the flag defaults to ON and silently rescales terms carrying a symmetric
amplitude pair on the partial-contraction path) and drops the stale
wick.hpp line reference. is_N_term keeps why rank > cutoff falls to R
rather than being dropped, and loses the paper quotes.
Also collapse hbar's five using-declarations into one and fix the range
notation in its error message.
CC::eom_r gains an optional block_ranks argument: a row-major K x K matrix
over the projection manifolds giving each block of the secular matrix its
own H̄ commutator truncation, instead of one uniform H̄ everywhere. The
manifolds are indexed by ASCENDING rank, so the qUCCSD ranks {2,1,1,0}
(10.1063/5.0062090 Table I, 10.1021/acs.jctc.5c01991 Table 1) serve EE, IP
and EA alike.
Each block is the sandwich <i|H̄|j> plus an explicit -E shift on the
diagonal at the block's own rank, not the commutator form <i|[H̄,r_j]|0>:
the commutator's extra -<i|r_j H̄^(k)|0> is manifold j's amplitude
residual, which vanishes only when k equals the rank the amplitudes were
converged against. Under the Bernoulli expansion each block's H̄ has its N
part removed for the same reason.
Empty block_ranks keeps the existing uniform path. That path commutes H̄
with an operator-level R, which the tensor-level Bernoulli H̄ cannot take
part in, so it now throws instead of aborting inside op.ipp. The block
shape and unitarity checks throw as well: SEQUANT_ASSERT compiles away
under SEQUANT_ASSERT_BEHAVIOR=IGNORE, and the shape check guards an
out-of-bounds read of block_ranks.
Pin the term counts of the {2,1,1,0} EE and IP sigma equations under the
Bernoulli expansion, and cover the three ways CC::eom_r rejects a
block_ranks argument: a non-square matrix, a non-unitary ansatz, and an
empty matrix under Bernoulli.
CMakeUserPresets.json is the documented per-developer companion to CMakePresets.json, and Notes is a symlink into a personal notes repo.
… guard Shortens the block-truncation derivation comments in cc.cpp/cc.hpp, adds missing paper section/equation references, and switches to std::ranges::reverse. The previous version of this commit also dropped the explicit Bernoulli-empty-block_ranks throw in CC::eom_r, reasoning that CC::hbar's own ctor-time hbar_comm_rank check made it redundant. That conflated two independent preconditions: hbar_comm_rank being set says nothing about block_ranks being non-empty. Without the guard, CC::eom_r() under Bernoulli falls through to the uniform path, which commutes the tensor-level Bernoulli H̄ with an operator-level R -- exactly the mixing op.ipp's commutes_with_atom() assert exists to catch. Restored the throw; caught by running the existing bernoulli_quccsd_eom unit test, which this repo's validation history (Notes/ucc/bernoulli-ucc/PR.md) never actually re-ran after this change was first made.
CMakeUserPresets.json and Notes are personal/local exclusions, not project-wide ignores, and don't belong in a Bernoulli-expansion feature PR. Use .git/info/exclude or a global excludesFile for these instead.
Trim prose that paraphrased or re-derived the paper's math in favor of citing the equation it corresponds to. Fix a few inaccuracies found by re-checking against 10.1063/1.5030344 and 10.1063/5.0062090 directly: R_part was mislabeled a 'rank-preserving remainder' (it isn't, for rank >= 2 operators), the DD EOM block was missing its bare-Fock content (f_ij, f_ab), and the 'higher orders carry only R-subscripted inner commutators' claim doesn't hold for the V_N-seeded terms. Also drops the invented 'Cancellation #2' label and renames the orphaned 'Cancellation #1' to 'the F-cancellation' throughout.
Ansatz, block_ranks shape and the Bernoulli/uniform-path mismatch are all caller errors on an internal precondition, matching how the rest of CC validates its configuration.
It simplified its argument in place, so the pointee the caller passed came back canonicalized; both in-library callers cloned first to avoid that. Clone inside instead.
The B_n listed at Eq. (40) are the paper's, i.e. B_n/n! in the textbook normalization; say so. F-cancellation needs Brillouin, not canonicality, and H̄⁰ does carry F. expand_to_blocks drops the other base spaces rather than finding them empty; the projection is what makes that harmless.
tʼ builds an operator-level similarity transform, which a tensor-level H̄ cannot enter; it reached the mismatch and aborted deep inside op algebra.
… assert tests Term counts are blind to the coefficients, so pin one projected equation in full. For a single manifold the blocked path must reproduce the uniform one, which shares no code with it. REQUIRE_THROWS_AS on a SEQUANT_ASSERT only holds in a THROW build.
The hole/particle candidate list falls back to all base spaces when empty; if that is empty too, the assignment loop indexed an empty vector and emitted a term carrying a garbage index space.
Cut the debugging narrative, the restated code, and the per-rank node counts; tighten the paragraphs that stay.
ajay-mk
force-pushed
the
ajay/feat/bernoulli-v2
branch
from
August 11, 2026 19:14
b7b4686 to
d022fa6
Compare
It read six accessors off the CC it was handed and needed N threaded in by hand because N is private, so it was already a member in all but name. Taking it as one drops the cc and N parameters. Also add the <algorithm> include std::ranges::reverse needs; it was only ever reaching it transitively.
The candidate-space loop kept a second vector to fall back on when the registry defined no hole/particle split. Every registry in convention.cpp sets is_hole and is_particle, so that branch was unreachable. Use the throwing accessors instead and drop the fallback; the assert now names what is actually required.
SEQUANT_ASSERT(EXPR, ...) takes the message itself, so the legacy `cond && "msg"` idiom leaves the string in the stringified condition.
hbar's nest() carried the tag validation twice and the tag-to-filter cascade twice, once for the V base and once per nesting level. One `part` lambda covers both, taking a flag for whether the input is already wick_reduce'd. Rename the loop variable to `cur`; it was called `op`, which shadows the mbpt::op namespace the surrounding function uses.
Checked against the cited papers:
- "F enters H̄ ONLY here" on H̄¹ contradicted the H̄⁰ = F + V line three
lines above. Only the F-commutators truncate at the first power of σ,
which is what 10.1063/1.5030344 states below Eq. (50).
- The {2,1,1,0} ranks were credited to Table 1 of 10.1021/acs.jctc.5c01991,
which lists which H̄ components enter each block and no ranks at all.
The ranks are the superscripts in its Fig. 1. Cite each for what it
carries, here and in cc.cpp and the unit test.
- eom_r no longer presents {2,1,1,0} as "qUCCSD". It reproduces the
published method only under the Bernoulli expansion, since BCH keeps
the N part the paper drops. Give the numbers as the ranks that paper
truncates qUCCSD at, and name its own UCCSD[2|2,1,0] notation.
- Drop the H_DD integral list, which was missing <ia||bj>, and point at
the paper instead of restating it.
Two conflicts, both in CC's assertions, both because master reworked the same lines this branch had added to. CC::CC: master rewrote the skip_singles assert to SEQUANT_ASSERT's two-argument form. Took master's wording and kept the Bernoulli precondition block after it. CC::tʼ: master dropped the redundant hbar_comm_rank assert, now covered by the constructor, and rewrote the pertbar one. Took master's version and kept the assert that rejects the Bernoulli expansion here. Master's op::N -> op::ã rename and its new CC::rdm do not touch the Bernoulli path; rdm builds its own transform via mbpt::lst rather than CC::hbar, so hbar_expansion_ does not reach it.
SEQUANT_ASSERT(EXPR, ...) takes the message itself, so the legacy `cond && "msg"` idiom folds the string into the stringified condition. Converts the twelve remaining sites in this file; conditions are unchanged. Master has been moving the same way, which is what both conflicts in the merge just before this were about. Roughly 90 `&&`-form sites remain elsewhere in SeQuant/ and tests/; converting those is a separate repo-wide sweep.
ajay-mk
marked this pull request as ready for review
August 17, 2026 23:43
Bernoulli required only unitary(), which admits oU; require U. CC::rdm ignored hbar_expansion_ and silently returned the BCH density, so reject it as tʼ does.
Removing it is a no-op only where the block rank equals hbar_comm_rank; below that rank the terms are off-shell, so the numbers change too.
Contributor
There was a problem hiding this comment.
Pull request overview
Adds a Bernoulli-number expansion path for the unitary-CC similarity-transformed Hamiltonian ( \bar H = e^{-\sigma} H e^{\sigma} ) (as an alternative to the existing BCH-based path), and extends EOM right-hand sigma construction to support per-block (\bar H) truncation (including Bernoulli’s tensor-level (\bar H) representation).
Changes:
- Introduces
mbpt::bernoulli::hbar()(and supporting Wick/N–R partition helpers) and wires it intombpt::CCviaCC::Options::hbar_expansion. - Adds a per-block-truncated EOM-CC right-hand path (
CC::eom_r(..., block_ranks)), used for Bernoulli and for explicit block-rank truncations. - Adds unit + integration tests that pin structural/term-count expectations across BCH vs Bernoulli.
Reviewed changes
Copilot reviewed 8 out of 9 changed files in this pull request and generated 1 comment.
Show a summary per file
| File | Description |
|---|---|
SeQuant/domain/mbpt/bernoulli.hpp |
Declares the Bernoulli (\bar H) API and supporting detail helpers (Wick reduction, block expansion, N/R split). |
SeQuant/domain/mbpt/bernoulli.cpp |
Implements Bernoulli (\bar H) assembly (orders 0–4) plus partial-contraction Wick machinery and N/R partitioning. |
SeQuant/domain/mbpt/models/cc.hpp |
Adds HbarExpansion, exposes it through CC::Options, and extends eom_r API docs for per-block truncation. |
SeQuant/domain/mbpt/models/cc.cpp |
Implements Bernoulli dispatch in hbar()/energy()/t(), blocks-based EOM path, and configuration validation. |
tests/unit/test_mbpt_cc.cpp |
Adds unit coverage for Bernoulli Wick behavior, block expansion, N/R split, (\bar H) structure, config validation, and EOM block truncation. |
tests/integration/ucc.cpp |
New integration executable pinning BCH vs Bernoulli term counts for UCC derivations. |
tests/integration/CMakeLists.txt |
Registers ucc.cpp integration runs (full and reduced sets). |
CMakeLists.txt |
Adds Bernoulli sources to the MBPT build. |
.gitignore |
Fixes the run entry formatting. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
evaleev
reviewed
Aug 24, 2026
evaleev
left a comment
evaleev
left a comment
Member
There was a problem hiding this comment.
Review pass over the whole PR (9 files, ~1.3k diff lines): the new bernoulli.{cpp,hpp}, the cc.{cpp,hpp} changes, and the supporting machinery they lean on.
Checked and found correct, for the record: the manifold enumeration in eom_r_blocked mirrors op::R's loop exactly and the result-index map min(bp,bh) is in range and collision-free for EE/IP/EA; the -E diagonal shift written as <i|r_i H̄|0> really does factor into E·<i|r_i|0> (the ket operators are all right-type and can only contract into the bra), and its absence off-diagonal is right; the nest memo keys are unambiguous and the memo/hbars iterator use is safe despite container::map being a flat_map; operator*/operator+/Product::append deep-copy, so memoized subexpressions survive later simplify calls; Product hashing excludes the scalar, so HashingAccumulator merges the add() wrappers correctly; and NormalOperator is an AbstractTensor, so get_used_indices/transform_expr reach its indices and TensorNetwork canonicalization preserves non-c-number ordering (the commutator does not collapse).
Findings are inline. Summary:
- high —
bernoulli.cpp: the N/R cutoff ignoresskip1, soskip_singles+ Bernoulli silently produces wrong UCCD equations. - medium —
cc.cpp: the diagonal shift uses each block's own truncation rank, so different diagonal blocks are referenced to different energy zeros. - low —
cc.cpp: on the BCH path, a uniformblock_ranksmatrix is not equivalent to an empty one, contrary to the header doc.
4-7. low — four preconditions expressed asSEQUANT_ASSERT, which compiles todo {} while(0)withoutSEQUANT_ASSERT_ENABLED; plus afind_nopgap onNormalOperatorSequenceand a dropped-proto-index issue under CSV.
Two claims were wrong against 10.1063/1.5030344: the expansion is Sec. II B, not III B, and F appears in H-bar^0 as well as H-bar^1. Others named fewer conditions than the code enforces. The rest is duplication, and three invented words give way to the block, E_gr and expansion that the papers use.
eom_r, eom_l and eom_r_blocked each descended (np, nh) with the same pair of break conditions. eom_manifolds returns the sequence lowest rank first, so eom_r_blocked no longer reverses it.
pins becomes a constexpr std::array, and the two-entry std::map becomes a function that also names the accepted values when the argument is not one of them.
is_N_term only checked an upper bound, so under skip_singles a rank-1 term was classified N and dropped even though sigma carries no singles amplitude to make it vanish.
Require an exact K x K row-major matrix and a unitary ansatz with release-safe Exceptions. Document the BCH-versus-Bernoulli dispatch and the scalar removed from each separately truncated diagonal block.
Reject CSV before partial Wick reduction can lose proto-index dependencies. Enforce the CC configuration in its constructor, and reject multiple residual normal operators even when assertions are disabled.
Clarify the single-reference projection limitation and document the indexed EOM result layout accurately.
Member
Author
|
This is now good to go. All review comments have been addressed:
All tests in SeQuant and MPQC still pass. |
evaleev
approved these changes
Aug 25, 2026
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.
Adds the Bernoulli expansion of the UCC similarity-transformed Hamiltonian$\bar{H} = e^{-\sigma} H e^{\sigma}$ , $\sigma = T - T^\dagger$ , as an alternative to the BCH expansion in $\bar H^0$ – $\bar H^4$ in Liu, Asthana, Cheng & Mukherjee, JCP 148, 244110 (2018)
mbpt::CC. Eq. (45) defines the sum, and Eqs. (46)–(50) defineCompanion PR: MPQC4 #792.
Changes
mbpt::CC, plus per-blockCC::Options::hbar_expansion: SelectsBCHorBernoulli; the former is the default. Bernoulli requiresAnsatz::Uandhbar_comm_rank.bernoulli::hbar(N, rank, skip1)assemblesdetail::wick_reduce— Wick-reduces while keeping partial contractions(needed so nesting produces operators, not scalars). Must run with
use_topology(false): its weight bookkeeping is only valid on the fully-contracted path, and it defaults to on.detail::N_part/detail::R_part— thebernoulli::hbar()returns a tensor-level expression (BCH'shbar()returns an operator-level one). Callers have to take care of this.CC::eom_r: per-block H̄ truncation (block_ranks):{2,1,1,0}. See Liu & Cheng, JCP 155, 174102 (2021), Sec. II C; Zhang & Liu, JCTC 22, 3341 (2026), Fig. 1, for the IP/EA analogues.hbar_comm_rank; below that rank, the values change too.Tests
tests/unit/test_mbpt_cc.cpp: Wick reduction, N/R split,tests/integration/ucc.cpp: BCH/Bernoulli term-count pins at ranks 2–4.qUCCSD Example
Important References: