Skip to content

feat(mbpt): square-root-normalized δr, δl projector operators #555

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
ajay-mk wants to merge 9 commits into
base: master
Choose a base branch
from
Draft

feat(mbpt): square-root-normalized δr, δl projector operators #555

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
262b45c
feat(core): handle Power in get_unique_indices
ajay-mk Jun 13, 2026
ef66412
feat(mbpt): add SquareRoot normalization support to OpMaker
ajay-mk Jun 13, 2026
290e4d1
feat(mbpt): make Normalization a parameter of operator builders
ajay-mk Jun 13, 2026
ddbe97d
feat(mbpt): add δr/δl projector ops
ajay-mk Jun 13, 2026
f139b34
refactor(mbpt/CC): use δl/δr in EOM derivations
ajay-mk Jun 13, 2026
a33e1de
docs: improve docs on `mbpt::Normalization` [skip ci]
ajay-mk Jun 16, 2026
d75c82e
fix: assert Constant or Variable base for Power in get_unique_indices
ajay-mk Jun 16, 2026
2edde54
docs: improve docs on `mbpt::Normalization`
ajay-mk Jun 16, 2026
7f92a62
fix: add `[[maybe_unused]]` in `get_unique_indices(Power)`
ajay-mk Jun 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions SeQuant/core/utility/indices.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,12 @@ IndexGroups<Container> get_unique_indices(const Variable&) {
return {};
}

template <typename Container = std::vector<Index>>
IndexGroups<Container> get_unique_indices([[maybe_unused]] const Power& power) {
SEQUANT_ASSERT(power.base()->is<Constant>() || power.base()->is<Variable>());
return {};
Comment thread
ajay-mk marked this conversation as resolved.
}

/// @returns Lists of non-contracted indices arising when contracting the two
/// given tensors in the order bra, ket, auxiliary
template <typename Container = container::svector<Index>>
Expand Down Expand Up @@ -249,6 +255,8 @@ IndexGroups<Container> get_unique_indices(const Expr& expr) {
return get_unique_indices<Container>(expr.as<Sum>());
} else if (expr.is<Product>()) {
return get_unique_indices<Container>(expr.as<Product>());
} else if (expr.is<Power>()) {
return get_unique_indices<Container>(expr.as<Power>());
} else {
throw Exception(
"Encountered unsupported expression type in get_unique_indices");
Expand Down
20 changes: 13 additions & 7 deletions SeQuant/domain/mbpt/models/cc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,6 +322,11 @@ std::vector<ExprPtr> CC::λʼ(size_t rank, size_t order,
return result;
}

namespace {
// EOM eigenvector operators R and L use SquareRoot normalization s
constexpr Normalization eom_norm = Normalization::SquareRoot;
} // namespace

std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
SEQUANT_ASSERT((np > 0 || nh > 0) && "Unsupported excitation order");
if (np != nh)
Expand All @@ -341,9 +346,9 @@ std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
// for unitary ansatz, we need to compute the commutator [hbar, R], otherwise
// just hbar * R is sufficient because ref_av uses connectivity
if (this->unitary()) {
hbar_R = commutator(hbar, R(np, nh));
hbar_R = commutator(hbar, R(np, nh, eom_norm));
} else {
hbar_R = hbar * R(np, nh);
hbar_R = hbar * R(np, nh, eom_norm);
}

// connectivity: empty for unitary ansatz, build otherwise
Expand All @@ -361,9 +366,10 @@ std::vector<ExprPtr> CC::eom_r(nₚ np, nₕ nh) const {
std::int64_t rp = np, rh = nh;
while (rp >= 0 && rh >= 0) {
if (rp == 0 && rh == 0) break;
// project with <rp, rh| (i.e., multiply P(rp, rh)) and compute VEV
// project with <rp, rh| (i.e., multiply δl(rp, rh)) and compute VEV
// use δl for consistent normalization
result.at(min(rp, rh)) =
this->ref_av(P(nₚ(rp), nₕ(rh)) * hbar_R, op_connect);
this->ref_av(δl(nₚ(rp), nₕ(rh)) * hbar_R, op_connect);
if (rp == 0 || rh == 0) break;
rp--;
rh--;
Expand All @@ -386,7 +392,7 @@ std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) const {
const auto hbar = this->hbar();

// L * hbar
const auto L_hbar = L(np, nh) * hbar;
const auto L_hbar = L(np, nh, eom_norm) * hbar;

// connectivity:
// default connections + connect H with projectors
Expand All @@ -405,9 +411,9 @@ std::vector<ExprPtr> CC::eom_l(nₚ np, nₕ nh) const {
std::int64_t rp = np, rh = nh;
while (rp >= 0 && rh >= 0) {
if (rp == 0 && rh == 0) break;
// right project with |rp,rh> (i.e., multiply P(-rp, -rh)) and compute VEV
// right project with |rp,rh> (i.e., multiply δr(rp, rh)) and compute VEV
result.at(min(rp, rh)) =
this->ref_av(L_hbar * P(nₚ(-rp), nₕ(-rh)), op_connect);
this->ref_av(L_hbar * δr(nₚ(rp), nₕ(rh)), op_connect);
if (rp == 0 || rh == 0) break;
rp--;
rh--;
Expand Down
Loading
Loading