Bernoulli Expansion based Unitary Coupled-Cluster - #581
Draft
ajay-mk wants to merge 16 commits into
Draft
Conversation
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.
ajay-mk
force-pushed
the
ajay/feat/bernoulli-v2
branch
from
July 27, 2026 12:08
2b73c55 to
e78e554
Compare
… 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.
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.
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$ , Eqs. (45)–(50) of Ref. 1.
mbpt::CC. RanksGround state energies are validated, excited states pending.
More details to be added.
References
self-consistent polarization propagator theory: A third-order formulation and
pilot applications, J. Chem. Phys. 148, 244110 (2018).
doi:10.1063/1.5030344 — Sec. III B, the
reference this module implements.