feat: CC RDMs - #578
Conversation
Derives the rank-r RDM as the reference expectation value of a
similarity-transformed number operator (op::N). Traditional ansatz uses
<0|(1+Λ) e^{-T} N e^{T}|0> with the linked-cluster connection constraint
enforced; unitary ansatz drops Λ and uses σ = T − T⁺ with genuine
commutators, needing no constraint.
Tests cover the traditional (rank 1 and 2) and unitary (rank 1) cases.
For a 1-body operator, we only need to consider up to 2 nested commutators (in traditional coupled-cluster) and for higher operators we use up to 4. For UCC Ansatz the user sets the commutator truncation rank.
CC::rdm now takes an optional commutator truncation rank, mirroring
CC::energy. The traditional-ansatz default is the exact termination
point of e^{-T} N e^{T} in <0|(1+Λ) N̄|0>: T-T contractions vanish, so
each T must contract with one of the number operator's 2*rank legs, and
the T legs left over can only be closed by the single Λ of rank <= N,
giving min(2*rank, rank + N).
The unitary ansatz has no such termination and keeps using
hbar_comm_rank.
|
Weird error because I don't use Copilot for anything 🙂 |
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 6 out of 6 changed files in this pull request and generated no new comments.
Suppressed comments (4)
SeQuant/domain/mbpt/models/cc.cpp:215
- CC::rdm() uses hbar_comm_rank_.value() when no comm_rank is provided for the unitary ansatz. In builds where SEQUANT_ASSERT is disabled, constructing a unitary CC without hbar_comm_rank will throw std::bad_optional_access here; it’s better to fail with a clear runtime error (or otherwise handle the missing value) instead of relying on asserts.
const auto commutator_rank = comm_rank.value_or(
unitary() ? hbar_comm_rank_.value() : std::min(2 * rank, rank + N));
auto Nbar = mbpt::lst(N_op, T(N, skip_singles()), commutator_rank,
{.unitary = unitary()});
SeQuant/domain/mbpt/models/cc.cpp:203
- Minor typo in the comment: there’s an extra ';' before “(see op::N)”, which makes the sentence read oddly.
This issue also appears on line 212 of the same file.
// 1. rank-`rank` number operator N = {a†_{p1..pr} a_{q1..qr}};(see op::N).
SeQuant/domain/mbpt/models/cc.hpp:168
- Doc comment has an extra space before the hyphen in “@p rank -particle”, which will render strangely in generated docs.
/// @brief derives the @p rank -particle reduced density matrix (RDM) as the
/// reference expectation value of a similarity-transformed number operator.
SeQuant/domain/mbpt/op.hpp:1169
- This Doxygen brief currently expands to something like “rank-rank number operator” because it uses “rank-\p rank”. Rephrasing avoids the duplicated word and reads more clearly in generated docs.
/// @brief Makes a rank-\p rank number (replacement) operator
/// \f$ \{ a^{p_1 \dots p_r}_{q_1 \dots q_r} \} \f$ over the complete
/// space.
…ion rank Lets a caller reuse an existing engine's other options while evaluating at a different commutator rank, without mutating an engine other derivations share.
…comm_rank access op::N was reachable but unregistered, so printing/introspecting it threw. Register it like its sibling operators, and add CC::rdm's own guard on hbar_comm_rank_ to match the other unitary-ansatz call sites in this class.
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 7 out of 7 changed files in this pull request and generated no new comments.
Suppressed comments (3)
SeQuant/domain/mbpt/models/cc.cpp:67
- CC::with_hbar_comm_rank accepts rank==0, which can silently disable the BCH expansion (mbpt::lst returns the untransformed operator when commutator_rank==0) and later break assumptions like CC::λ requiring hbar_comm_rank>=1. Consider validating the input here to keep the API consistent with other CC entry points.
CC CC::with_hbar_comm_rank(size_t rank) const {
CC result(*this);
result.hbar_comm_rank_ = rank;
return result;
}
SeQuant/domain/mbpt/models/cc.cpp:223
- For unitary ansatz, CC::rdm uses hbar_comm_rank_ as the default commutator truncation, but the code doesn’t guard against hbar_comm_rank_ being 0 (which makes N̄ == N and the VEV trivially 0 for the normal-ordered number operator). Adding a sanity check when the default is used prevents surprising zero RDMs.
if (unitary())
SEQUANT_ASSERT(hbar_comm_rank_ &&
"hbar_comm_rank must be specified for unitary ansatz in "
"CC::rdm");
const auto commutator_rank = comm_rank.value_or(
unitary() ? hbar_comm_rank_.value() : std::min(2 * rank, rank + N));
SeQuant/domain/mbpt/models/cc.hpp:173
- The Doxygen text "@p rank -particle" is awkward/ambiguous (the parameter markup breaks the hyphenation). Rewording avoids a formatting glitch and reads more clearly.
/// @brief derives the @p rank -particle reduced density matrix (RDM) as the
/// reference expectation value of a similarity-transformed number operator.
|
This will need changes once #584 goes through. |
evaleev
left a comment
There was a problem hiding this comment.
Inline version of the main points from my earlier summary comment. Design is sound overall — op::N deliberately bypassing OpMaker so its indices survive as γ's free indices is the right call, and the connectivity machinery in vac_av.cpp lines up correctly with an op that lowers to a bare FNOperator. The items below are ordered roughly by how hard they'd be to change after MPQC binds against this.
| /// A general operator of rank \p K | ||
| ExprPtr θ(std::size_t K); | ||
|
|
||
| /// @brief Makes a rank-\p rank number (replacement) operator |
There was a problem hiding this comment.
Naming — the item hardest to change later.
N conventionally denotes the particle-number operator Σ_p a†_p a_p. This is a rank-r replacement operator (rank-1 is Ê^p_q), which is why the docs have to say "number (replacement) operator" everywhere. Two concrete costs:
CC::rdmhas to work around the collision withCC::N(the cluster rank) via anop::qualification plus an explanatory comment.OpRegistry::validate_opthrows if a label is already registered, so addingNto both shipped registries breaks any downstream code that registers its ownN— a very likely user label.
Would E, Ê, or γ̂ work? Renaming after MPQC ships against this is much more expensive than now.
There was a problem hiding this comment.
Clarifying what this actually is, since it bears on the name: it's a single replacer, not the RDM operator Σ γ^{i…}_{j…} a^{j…}_{i…}.
The tensor-form generator returns ex<FNOperator>(cre(cres), ann(anns)) and nothing else — no Tensor factor, no normalization, no summation — over fixed non-tmp ordinals, so the indices stay free. It is exactly {ã^{p_1…p_r}_{p_{r+1}…p_{2r}}}, one component.
The contrast with θ makes it sharp: tensor::θ(K) goes through OpMaker::make, which returns apply_normalization(θ_tensor * NormalOperator(...)) over make_tmp_index indices, i.e. 1/(r!)² Σ θ^{…}_{…} ã^{…}_{…}. op::N is θ with the amplitude tensor and the normalization stripped off — which is precisely why it has to bypass OpMaker.
So CC::rdm(r) returns the element γ^{p_1…p_r}{p{r+1}…p_{2r}} as an expression in free complete-space indices, not Σ γ a. That's the right primitive for the stated purpose (codegen of a γ tensor); the Σ γ a form would serve the opposite direction — evaluating expectation values against a density you already have. The construction is right; only the name is wrong.
On the name
E is out, as noted elsewhere: SeQuant already prints the spin-free replacer as ã, and E^p_q is just the textbook name for that same object, so E would be a synonym for something that exists.
Out for the same class of reason: γ̂/Γ̂ mislabel it as the density operator — the exact error N makes in the other direction.
Reserved labels are only {Â, Ŝ, P̂, δ, s}, so the field is otherwise open; the real constraints are avoiding index base keys (p, a, i, …) and existing registry entries (h, f, f̃, g, θ, t, λ, R, L, …).
Suggestion: make the label mirror the notation it is — registry label L"ã", function op::replacer(rank) (or op::ã). Seeing ã at the operator level and ã^{p_1}_{p_2} in the lowered form is the same object twice, which is the point. It's unused, unreserved, and not an index key.
If a Greek single letter is preferred for symmetry with θ, ν is free and semantically uncommitted. I'd avoid τ — CC literature already uses it for the T₂+T₁T₁ combination.
There was a problem hiding this comment.
Renamed to ã: registry label L"ã" in both shipped registries, and function op::ã in mbpt namespace. 31289a6
| } | ||
|
|
||
| ExprPtr N(std::size_t rank) { | ||
| if (rank < 1) throw Exception("mbpt::op::N: rank must be >= 1"); |
There was a problem hiding this comment.
Diverges from its siblings on precondition handling: θ, t, λ all do SEQUANT_ASSERT(K > 0) then the registry check. Here it's a throw first, registry assert second.
Also rank < 1 on a std::size_t is just rank == 0 — worth spelling that way.
There was a problem hiding this comment.
Done: SEQUANT_ASSERT(rank > 0, ...) then the registry check, matching θ/t/λ. rank < 1 is now rank > 0. 31289a6
| for (std::size_t i = 0; i < rank; ++i) | ||
| cres.emplace_back(space, i + 1); | ||
| for (std::size_t i = 0; i < rank; ++i) | ||
| anns.emplace_back(space, rank + i + 1); |
There was a problem hiding this comment.
The free-index convention is the thing callers actually need, and it's currently undocumented.
This emits creators p_1…p_r and annihilators p_{r+1}…p_{2r}, so CC::rdm(2) returns an expression whose ket indices are literally p_3, p_4. The docs (both here and in cc.hpp) write q_1…q_r, which corresponds to no label in the output — MPQC has to hard-code the real mapping. Please state the exact labels.
(Non-tmp ordinals are the right choice — they can't collide with the Index::make_tmp_index ordinals used by T/Λ. It's only the documentation that's missing.)
Nit: the two loops can be one.
There was a problem hiding this comment.
Stated on op::ã as a contract. Creators use ordinals 1 ... rank, annihilators rank+1 ... 2 * rank. MPQC already relies on it for closed_shell_spintrace call.
Two loops are now unified. Two checks on the indices, not one: to_latex(ã(2)->tensor_form()) in test_mbpt.cpp pins the operator's own ordinals, and at the RDM level the traditional Γ's free indices are pinned with get_unique_indices. c3d1db0
| /// \f$ \gamma^{p_1 \dots p_r}_{q_1 \dots q_r} = \langle 0| (1 + | ||
| /// \hat{\Lambda}) e^{-\hat{T}} \{ a^{p_1 \dots p_r}_{q_1 \dots q_r} \} | ||
| /// e^{\hat{T}} |0 \rangle \f$; unitary ansatz drops \f$ \hat{\Lambda} \f$ and | ||
| /// uses \f$ e^{-\hat{\sigma}} \dots e^{\hat{\sigma}} \f$ with \f$ | ||
| /// \hat{\sigma} = \hat{T} - \hat{T}^\dagger \f$. The number-operator indices | ||
| /// are in the complete space, so they survive as the free indices of the | ||
| /// result. Returns the *linked* density (the number | ||
| /// operator is required to connect to the cluster amplitudes). |
There was a problem hiding this comment.
Three doc issues in this block:
\gamma^{p_1 \dots p_r}_{q_1 \dots q_r}— there is noq_*in the output; the ket indices arep_{r+1}…p_{2r}(see my comment onop.cpp).- The reference contribution is not included.
ref_avreturns 0 for a bare normal-ordered operator, so this is the correlation part only and the caller addsδ. The@returnline hints at it ("correlation part"), but the leading formula doesn't, and that's the part people read. - "Returns the linked density (the number operator is required to connect to the cluster amplitudes)" is true only for the traditional branch — the unitary branch passes an empty
OpConnections{}. Consistent withCC::energy, but the sentence reads as unconditional.
Worth adding: the result is not manifestly antisymmetric under p_1↔p_2 / p_{r+1}↔p_{r+2} for rank ≥ 2, since N carries no antisymmetrizer.
There was a problem hiding this comment.
The formula reads \gamma^{p_1 \dots p_r}_{p_{r+1} \dots p_{2r}} over \tilde{a}, so it matches what is emitted. The missing reference contribution, the traditional-only scope of "linked", and the rank >= 2 antisymmetry caveat each got their own @note.
Kept "not manifestly antisymmetric": the claim is that the form carries no antisymmetrizer, not that the object lacks the symmetry.
Also in this block: @param comm_rank referred to \bar{N}, which no longer names anything, and the min(2*rank, rank + N) argument moved out of the commit message into cc.cpp.
| if (unitary()) | ||
| SEQUANT_ASSERT(hbar_comm_rank_ && | ||
| "hbar_comm_rank must be specified for unitary ansatz in " | ||
| "CC::rdm"); |
There was a problem hiding this comment.
This guard is dead code: the CC constructor already asserts hbar_comm_rank_ for the unitary ansatz, and with_hbar_comm_rank only ever sets it — so it can never fire.
Worse, in assert-disabled builds both asserts vanish and the actual guard becomes hbar_comm_rank_.value() on the next line → std::bad_optional_access, which conflicts with the repo's preference for sequant::Exception. Either drop this, or make it a real throw Exception(...) that survives release builds.
Minor: SEQUANT_ASSERT now takes a message argument, so SEQUANT_ASSERT(hbar_comm_rank_, "…") beats the cond && "msg" idiom — though the rest of this file uses the old form.
There was a problem hiding this comment.
Agree with splitting the assertion, but SEQUANT_ASSERT is used throughout in the CC class and if we want to switch to sequant::Exception we should switch everything.
There was a problem hiding this comment.
Dropped, per the first of your two options. 5743c9e
It was dead in three places rather than one — tʼ and eom_r carried the same copy — so the constructor is now the sole guard. tʼ keeps its pertbar_comm_rank check, which the constructor does not cover. Message-argument form used on the lines this touches.
| "hbar_comm_rank must be specified for unitary ansatz in " | ||
| "CC::rdm"); | ||
| const auto commutator_rank = comm_rank.value_or( | ||
| unitary() ? hbar_comm_rank_.value() : std::min(2 * rank, rank + N)); |
There was a problem hiding this comment.
The min(2*rank, rank + N) bound is the load-bearing correctness claim in this PR, and right now it's justified only in a commit message — nothing in CI checks it.
The argument checks out as I read it (each T must consume ≥1 of N's 2r legs ⇒ k ≤ 2r; and T's 2k quasi-creators must be absorbed by N's 2r + Λ's 2N legs ⇒ k ≤ r + N), but it deserves a test rather than a reader's re-derivation. See my comment on the test file.
| /// @brief returns a copy of this object whose H̄ commutator truncation rank | ||
| /// is @p rank, leaving this engine unchanged; all other options are carried | ||
| /// over | ||
| [[nodiscard]] CC with_hbar_comm_rank(size_t rank) const; |
There was a problem hiding this comment.
No validation that rank >= 1. CC::λ computes hbar(commutator_rank - 1) guarded only by a SEQUANT_ASSERT; with asserts disabled and hbar_comm_rank_ == 0 that underflows to SIZE_MAX, and mbpt::lst then loops for (k = 1; k <= SIZE_MAX; ++k) — a hang.
This is reachable today via Options{.hbar_comm_rank = 0}, so it's pre-existing, but the new setter adds a second door and is the natural place to close it. A SEQUANT_ASSERT(rank >= 1) here (and ideally in the Options validation) is cheap.
Separately: this has no caller in the PR — rdm takes comm_rank directly. Fine as a general utility, just noting it's unrelated scope.
There was a problem hiding this comment.
This is actually fine and hbar_comm_rank = 0 is a valid use case, and only is a problem for CC::λ. CC::λ has an assertion to catch this case. The convention in the class is to use SEQUANT_ASSERT over sequant::Exception.
with_hbar_comm_rank does not have a caller, but is used in MPQC for a non-RDM case. I can move this out of the PR is that is a blocker.
| auto g = CC{N}.rdm(1); | ||
| REQUIRE(g); | ||
| REQUIRE(size(g) == 10); | ||
| auto G = CC{N}.rdm(2); | ||
| REQUIRE(G); | ||
| REQUIRE(size(G) == 94); | ||
| // an explicit comm_rank truncates early: cutting Γ at the 2nd nested | ||
| // commutator drops the terms the default (4th, exact) picks up | ||
| auto G2 = CC{N}.rdm(2, 2); | ||
| REQUIRE(G2); | ||
| REQUIRE(size(G2) == 71); // 23 terms short of the exact 94 |
There was a problem hiding this comment.
These are term-count-only assertions, while the rest of this file uses REQUIRE_THAT(..., EquivalentTo(L"...")) for correctness (see the λ and energy sections). Counts drift silently with any canonicalization or simplify change and don't verify the equations at all — so correctness currently rests entirely on the MPQC cross-check, which this repo's CI can't see.
Highest-value additions, in order:
- Prove the default
comm_rankis exact. One assertion turns the commit-message argument into something CI enforces:Much stronger thanREQUIRE_THAT(CC{2}.rdm(2), EquivalentTo(CC{2}.rdm(2, 6))); // default is converged REQUIRE_THAT(CC{2}.rdm(1), EquivalentTo(CC{2}.rdm(1, 5)));
== 94, and it's the claim most likely to be wrong. - Spell out the 1-RDM with
EquivalentTo— 10 terms is small enough, and it pins the free-index convention (p_1/p_2) into the suite. - Regression test for the bug fixed in the last commit:
op::Nwas unregistered and threw on printing. Nothing here exercisesto_latex(N(1)). - Missing:
op::N(0)throw, UCC rank-2.
The test_mbpt.cpp addition is the right shape by contrast — it checks structure (ncreators, nannihilators, spaces) rather than counts.
There was a problem hiding this comment.
Both 1-RDMs are now spelled out with EquivalentTo; the to_latex regression, the rank-0 precondition, and the unitary rank-2 density are now tested.
For the 2-RDMs the expressions stay term counts (94, 71, 24) — 94 terms is not worth spelling out. The rank-2 path is not count-only though: TCC RDM default comm_rank compares rdm(2) against rdm(2, 6) as full expressions, and the traditional Γ has its free indices pinned with get_unique_indices (ket {p_1, p_2}, bra {p_3, p_4}, no aux). c3d1db0
evaleev
left a comment
There was a problem hiding this comment.
The physics and the construction are right — op::N bypassing OpMaker so its indices survive as γ's free indices is the correct primitive, the connectivity machinery lines up with an op that lowers to a bare FNOperator, and CI is green across all 16 checks. Details are in the inline comments; blocking items collected here.
Blocking
-
The name
N. It denotes a single replacer{ã^{p_1…p_r}_{p_{r+1}…p_{2r}}}, not the number operator and not the density operator — and it squats a very likely user label in both shipped registries, whereOpRegistry::validate_opturns a downstream re-registration into a throw. This is the piece that gets expensive once MPQC binds against it. See the thread onop.hppfor the reasoning and concrete alternatives. -
Document the free-index convention. The result's ket indices are literally
p_{r+1}…p_{2r}; the docs sayq_1…q_r, which matches no label in the output. Callers can't use this without knowing the real mapping. Same block should state that the reference δ contribution is not included and that the result isn't manifestly antisymmetric for rank ≥ 2. -
Test that the default
comm_rankis exact.min(2*rank, rank + N)is the load-bearing correctness claim and is currently justified only in a commit message.REQUIRE_THAT(CC{2}.rdm(2), EquivalentTo(CC{2}.rdm(2, 6)))converts it into something CI enforces — much stronger thansize(G) == 94, which drifts silently with any canonicalization change.
Should fix, not blocking
- The unitary guard in
CC::rdmis dead code (the ctor already asserts it), and in assert-disabled builds the real guard becomes.value()→std::bad_optional_access. with_hbar_comm_rankdoesn't validaterank >= 1; viaCC::λ'scommutator_rank - 1that underflows toSIZE_MAXand hangsmbpt::lstin a release build.op::Nthrows whereθ/t/λassert, with the registry check in the opposite order;rank < 1on an unsigned isrank == 0.- The "linked density" sentence is true only for the traditional branch.
Happy to re-review promptly — most of the list is small.
`N` denotes the particle-number operator by convention, but this is a
rank-r replacement operator -- a single replacer {ã^{p_1..p_r}_{p_{r+1}..p_{2r}}},
which is why the docs had to keep saying "number (replacement) operator".
Two concrete costs: `CC::rdm` had to qualify it as `op::N` to dodge the
`CC::N` member, and registering `N` in both shipped registries would make
`OpRegistry::validate_op` throw for any downstream code registering its own
`N` -- a very likely user label.
The registry label mirrors the notation the operator lowers to, so `ã`
appears at both the operator and the tensor level.
Also fold the two index loops into one and match the sibling operators
(`θ`, `t`, `λ`) by asserting the rank precondition instead of throwing.
Documentation. The formula wrote the ket indices as q_1..q_r, which match
no label in the output -- they are p_{r+1}..p_{2r}, and MPQC has to
hard-code that mapping. Point at op::ã, which now states the ordinals.
Three caveats a caller cannot infer were also missing: the reference
contribution is absent (ref_av returns 0 for a bare normal-ordered
operator, so this is the correlation part and the caller adds delta), the
result is not manifestly antisymmetric for rank >= 2, and "linked" holds
only for the traditional branch since the unitary one passes empty
connections. Record the argument for the min(2*rank, rank + N) default in
the code rather than leaving it in a commit message.
Assertions. The unitary hbar_comm_rank check was dead in three places --
the constructor already asserts it and nothing else clears it -- so drop
the copies in rdm, tʼ and eom_r and keep the constructor as the sole
guard. tʼ keeps its pertbar_comm_rank check, which the constructor does
not cover. Use the message argument form on the lines this touches.
The rdm coverage was term counts only, which drift silently with any canonicalization or simplify change and verify none of the equations. The load-bearing claim is that the default comm_rank, min(2*rank, rank + N), is already exact; raising it must add nothing. Assert that directly instead of leaving the argument in a commit message. Spell out both 1-RDMs with EquivalentTo. Ten and eight terms are small enough to read, and they pin the free-index convention -- p_1 is a ket index and p_2 a bra index, since ã's creators reach the expression through ket slots and its annihilators through bra slots -- which MPQC reconstructs by hand. The 2-RDMs stay counts -- 94 terms is not worth spelling out -- but the traditional one also checks its free indices directly, since that is the shape MPQC evaluates. Also cover what was untested: to_latex on the operator form, which threw before op::ã was registered; the rank 0 precondition; and the unitary rank-2 density.
CC::rdm was the one mbpt::lst() call site passing options inline rather
than lst_options(), which left use_connected_form at its false default
while every other site got the connected-product form. That contradicted
lst_options()'s own doc ("the LSTOptions this engine uses for every
mbpt::lst() call") and did redundant work: the traditional branch already
hands ref_av the {ã,t} connectivity that makes the cheaper form
equivalent, and the unitary branch gets use_connected_form == false from
lst_options() anyway, so nothing changes there.
Verified inert rather than assumed: the full to_latex serialization of
rdm(1), rdm(2), rdm(2,2), rdm(1,5), rdm(2,6) and both unitary densities
is byte-identical before and after.
c3d1db0 to
5ce4bfd
Compare
Resolves #572. Needed for for the real-time project.
The equations are numerically verified in MPQC (T and U ansatz, SO/CS, CCSD and CCSDT). MPQC PR: https://github.com/ValeevGroup/mpqc4/pull/789