Skip to content

Add manuscript examples to test suite #451

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

Open
ajay-mk wants to merge 18 commits into
base: master
Choose a base branch
from
Open

Add manuscript examples to test suite #451

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
98c36a0
test(examples): add initial set of examples from MSII
ajay-mk Dec 10, 2025
62d8050
test(examples): add CC perturbed amplitude derivations
ajay-mk Dec 10, 2025
90a835a
refactor: rename PT related functions
ajay-mk Dec 15, 2025
05c8f61
feat(examples): add `ω` terms to PT amplitude examples
ajay-mk Dec 15, 2025
d52e771
refactor(CC): rename PT methods according to bf96824
ajay-mk Dec 15, 2025
1ce3ce4
fix: more fixes for PT related ops
ajay-mk Dec 15, 2025
7066a7e
test(examples): introduce `H̅()` method for sim. transformed Hamilton…
ajay-mk Dec 15, 2025
77718cb
test(examples): add ground state amplitude derivations
ajay-mk Dec 15, 2025
d62cd15
test(examples): proper verification of EOM examples
ajay-mk Dec 16, 2025
21f927a
test(examples): reorganize + docs
ajay-mk Dec 16, 2025
826da69
refac(mbpt): fixed rank operators use lowercase identifiers. `_` suff…
ajay-mk Dec 17, 2025
cb4d391
fix(examples): fix 62d8050
ajay-mk Dec 17, 2025
880e5c0
fix(benchmark): update benchmark cases according to 826da69
ajay-mk Dec 17, 2025
fabc159
fix(python): Trying to fix pysequant issue
ajay-mk Dec 17, 2025
a3199b0
fix(ci): Ninja uses only 2 threads
ajay-mk Dec 18, 2025
af7a776
Merge remote-tracking branch 'origin/master' into ajay/examples/MS-II
ajay-mk Dec 19, 2025
284b5f9
refac(CC): use arithmetic overloads of `Constant` and `Variable` for …
ajay-mk Dec 22, 2025
0a3fe82
cleanup
ajay-mk Dec 22, 2025
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
33 changes: 16 additions & 17 deletions SeQuant/domain/mbpt/models/cc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,26 +139,26 @@ std::vector<ExprPtr> CC::λ(size_t commutator_rank) {
return result;
}

std::vector<ExprPtr> CC::t_pt(size_t rank, size_t order,
std::optional<size_t> nbatch) {
std::vector<ExprPtr> CC::tʼ(size_t rank, size_t order,
std::optional<size_t> nbatch) {
SEQUANT_ASSERT(order == 1 &&
"sequant::mbpt::CC::t_pt(): only first-order perturbation is "
"sequant::mbpt::CC::tʼ(): only first-order perturbation is "
"supported now");
SEQUANT_ASSERT(rank == 1 &&
"sequant::mbpt::CC::t_pt(): only one-body perturbation "
"sequant::mbpt::CC::tʼ(): only one-body perturbation "
"operator is supported now");
SEQUANT_ASSERT(ansatz_ == Ansatz::T && "unitary ansatz is not yet supported");

// construct h1_bar
// truncate h1_bar at rank 2 for one-body perturbation
// operator and at rank 4 for two-body perturbation operator
const auto h1_truncate_at = rank == 1 ? 2 : 4;
const auto h1_bar = mbpt::lst(H_pt(rank, {.order = order, .nbatch = nbatch}),
const auto h1_bar = mbpt::lst(Hʼ(rank, {.order = order, .nbatch = nbatch}),
T(N), h1_truncate_at);

// construct [hbar, T(1)]
const auto hbar_pert =
mbpt::lst(H(), T(N), 3) * T_pt(N, {.order = order, .nbatch = nbatch});
mbpt::lst(H(), T(N), 3) * Tʼ(N, {.order = order, .nbatch = nbatch});

// [Eq. 34, WIREs Comput Mol Sci. 2019; 9:e1406]
const auto expr = simplify(h1_bar + hbar_pert);
Expand All @@ -176,20 +176,20 @@ std::vector<ExprPtr> CC::t_pt(size_t rank, size_t order,
std::vector<ExprPtr> result(N + 1);
for (auto p = N; p >= 1; --p) {
const auto freq_term = ex<Variable>(L"ω") * P(nₚ(p)) *
T_pt_(p, {.order = order, .nbatch = nbatch});
Tʼ_(p, {.order = order, .nbatch = nbatch});
result.at(p) =
this->ref_av(P(nₚ(p)) * expr, op_connect) - this->ref_av(freq_term);
}
return result;
}

std::vector<ExprPtr> CC::λ_pt(size_t rank, size_t order,
std::optional<size_t> nbatch) {
std::vector<ExprPtr> CC::λʼ(size_t rank, size_t order,
std::optional<size_t> nbatch) {
SEQUANT_ASSERT(order == 1 &&
"sequant::mbpt::CC::λ_pt(): only first-order perturbation is "
"sequant::mbpt::CC::λʼ(): only first-order perturbation is "
"supported now");
SEQUANT_ASSERT(rank == 1 &&
"sequant::mbpt::CC::λ_pt(): only one-body perturbation "
"sequant::mbpt::CC::λʼ(): only one-body perturbation "
"operator is supported now");
SEQUANT_ASSERT(ansatz_ == Ansatz::T && "unitary ansatz is not yet supported");

Expand All @@ -200,18 +200,17 @@ std::vector<ExprPtr> CC::λ_pt(size_t rank, size_t order,
// truncate h1_bar at rank 2 for one-body perturbation
// operator and at rank 4 for two-body perturbation operator
const auto h1_truncate_at = rank == 1 ? 2 : 4;
const auto h1_bar = mbpt::lst(H_pt(rank, {.order = order, .nbatch = nbatch}),
const auto h1_bar = mbpt::lst(Hʼ(rank, {.order = order, .nbatch = nbatch}),
T(N), h1_truncate_at);

// construct [hbar, T(1)]
const auto hbar_pert =
mbpt::lst(H(), T(N), 3) * T_pt(N, {.order = order, .nbatch = nbatch});
mbpt::lst(H(), T(N), 3) * Tʼ(N, {.order = order, .nbatch = nbatch});

// [Eq. 35, WIREs Comput Mol Sci. 2019; 9:e1406]
const auto One = ex<Constant>(1);
const auto expr =
simplify((One + Λ(N)) * (h1_bar + hbar_pert) +
Λ_pt(N, {.order = order, .nbatch = nbatch}) * hbar);
const auto expr = simplify((One + Λ(N)) * (h1_bar + hbar_pert) +
Λʼ(N, {.order = order, .nbatch = nbatch}) * hbar);

// connectivity:
// t and t1 with {h,f,g}
Expand All @@ -236,7 +235,7 @@ std::vector<ExprPtr> CC::λ_pt(size_t rank, size_t order,
std::vector<ExprPtr> result(N + 1);
for (auto p = N; p >= 1; --p) {
const auto freq_term = ex<Variable>(L"ω") *
Λ_pt_(p, {.order = order, .nbatch = nbatch}) *
Λʼ_(p, {.order = order, .nbatch = nbatch}) *
P(nₚ(-p));
result.at(p) =
this->ref_av(expr * P(nₚ(-p)), op_connect) + this->ref_av(freq_term);
Expand Down
4 changes: 2 additions & 2 deletions SeQuant/domain/mbpt/models/cc.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ class CC {
/// @pre `rank==1 && order==1`, only first order perturbation and one-body perturbation operator is supported now
/// @return std::vector of perturbed t amplitude equations
// clang-format on
[[nodiscard]] std::vector<ExprPtr> t_pt(
[[nodiscard]] std::vector<ExprPtr> tʼ(
size_t rank = 1, size_t order = 1,
std::optional<size_t> nbatch = std::nullopt);

Expand All @@ -100,7 +100,7 @@ class CC {
/// @pre `rank==1 && order==1`, only first order perturbation and one-body perturbation operator is supported now
/// @return std::vector of perturbed λ amplitude equations
// clang-format on
[[nodiscard]] std::vector<ExprPtr> λ_pt(
[[nodiscard]] std::vector<ExprPtr> λʼ(
size_t rank = 1, size_t order = 1,
std::optional<size_t> nbatch = std::nullopt);

Expand Down
64 changes: 32 additions & 32 deletions SeQuant/domain/mbpt/op.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -811,58 +811,58 @@ ExprPtr S(std::int64_t K) {
OpType::S, cre(creators), ann(annihilators))(dep, {Symmetry::Nonsymm});
}

ExprPtr H_pt(std::size_t R, const OpParams& params) {
ExprPtr Hʼ(std::size_t R, const OpParams& params) {
params.validate();
SEQUANT_ASSERT(params.order == 1 &&
"sequant::mbpt::tensor::H_pt: only supports first "
"sequant::mbpt::tensor::Hʼ: only supports first "
"order perturbation");
SEQUANT_ASSERT(R > 0 && "Operator rank must be > 0");
return OpMaker<Statistics::FermiDirac>(OpType::h_1, ncre(R), nann(R),
params)();
}

ExprPtr T_pt_(std::size_t K, const OpParams& params) {
ExprPtr Tʼ_(std::size_t K, const OpParams& params) {
params.validate();
SEQUANT_ASSERT(params.order == 1 &&
"sequant::mbpt::tensor::T_pt_: only supports first "
"sequant::mbpt::tensor::Tʼ_: only supports first "
"order perturbation");
SEQUANT_ASSERT(K > 0 && "Operator rank must be > 0");
return OpMaker<Statistics::FermiDirac>(OpType::t_1, ncre(K), nann(K),
params)();
}

ExprPtr T_pt(std::size_t K, const OpParams& params) {
ExprPtr Tʼ(std::size_t K, const OpParams& params) {
params.validate();
if (params.skip1) SEQUANT_ASSERT(K > 1);
ExprPtr result;
for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) {
result += tensor::T_pt_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
result += tensor::Tʼ_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
}
return result;
}

ExprPtr Λ_pt_(std::size_t K, const OpParams& params) {
ExprPtr Λʼ_(std::size_t K, const OpParams& params) {
params.validate();
SEQUANT_ASSERT(params.order == 1 &&
"sequant::mbpt::tensor::Λ_pt_: only supports first "
"sequant::mbpt::tensor::Λʼ: only supports first "
"order perturbation");
SEQUANT_ASSERT(K > 0 && "Operator rank must be > 0");
return OpMaker<Statistics::FermiDirac>(OpType::λ_1, ncre(K), nann(K),
params)();
}

ExprPtr Λ_pt(std::size_t K, const OpParams& params) {
ExprPtr Λʼ(std::size_t K, const OpParams& params) {
params.validate();
if (params.skip1) SEQUANT_ASSERT(K > 1);
ExprPtr result;
for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) {
result += tensor::Λ_pt_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
result += tensor::Λʼ_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
}
return result;
}
Expand Down Expand Up @@ -1053,62 +1053,62 @@ ExprPtr make_op_from_params(std::function<std::wstring_view()> label_gen,
}
} // anonymous namespace

ExprPtr H_pt(std::size_t R, const OpParams& params) {
ExprPtr Hʼ(std::size_t R, const OpParams& params) {
SEQUANT_ASSERT(R > 0);
SEQUANT_ASSERT(params.order == 1 &&
"only first order perturbation is supported now");

return make_op_from_params(
[]() -> std::wstring_view { return optype2label.at(OpType::h_1); },
[R, params]() -> ExprPtr { return tensor::H_pt(R, params); },
[R, params]() -> ExprPtr { return tensor::Hʼ(R, params); },
[R](qnc_t& qns) { qns = combine(general_type_qns(R), qns); }, params);
}

ExprPtr T_pt_(std::size_t K, const OpParams& params) {
ExprPtr Tʼ_(std::size_t K, const OpParams& params) {
SEQUANT_ASSERT(K > 0);
SEQUANT_ASSERT(params.order == 1 &&
"only first order perturbation is supported now");

return make_op_from_params(
[]() -> std::wstring_view { return optype2label.at(OpType::t_1); },
[K, params]() -> ExprPtr { return tensor::T_pt_(K, params); },
[K, params]() -> ExprPtr { return tensor::Tʼ_(K, params); },
[K](qnc_t& qns) { qns = combine(excitation_type_qns(K), qns); }, params);
}

ExprPtr T_pt(std::size_t K, const OpParams& params) {
ExprPtr Tʼ(std::size_t K, const OpParams& params) {
params.validate();
SEQUANT_ASSERT(K > (params.skip1 ? 1 : 0));
ExprPtr result;
for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) {
result += T_pt_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
result += Tʼ_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
}
return result;
}

ExprPtr Λ_pt_(std::size_t K, const OpParams& params) {
ExprPtr Λʼ_(std::size_t K, const OpParams& params) {
SEQUANT_ASSERT(K > 0);
SEQUANT_ASSERT(params.order == 1 &&
"only first order perturbation is supported now");

return make_op_from_params(
[]() -> std::wstring_view { return optype2label.at(OpType::λ_1); },
[K, params]() -> ExprPtr { return tensor::Λ_pt_(K, params); },
[K, params]() -> ExprPtr { return tensor::Λʼ_(K, params); },
[K](qnc_t& qns) { qns = combine(deexcitation_type_qns(K), qns); },
params);
}

ExprPtr Λ_pt(std::size_t K, const OpParams& params) {
ExprPtr Λʼ(std::size_t K, const OpParams& params) {
params.validate();
SEQUANT_ASSERT(K > (params.skip1 ? 1 : 0));
ExprPtr result;
for (auto k = (params.skip1 ? 2ul : 1ul); k <= K; ++k) {
result += Λ_pt_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
result += Λʼ_(k, {.order = params.order,
.nbatch = params.nbatch,
.batch_ordinals = params.batch_ordinals,
.skip1 = false});
}
return result;
}
Expand Down
28 changes: 14 additions & 14 deletions SeQuant/domain/mbpt/op.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1068,39 +1068,39 @@ ExprPtr S(std::int64_t K);
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr H_pt(std::size_t R, const OpParams& params = {.order = 1});
ExprPtr Hʼ(std::size_t R, const OpParams& params = {.order = 1});

/// @brief Makes perturbed particle-conserving excitation operator
/// @param K rank of the excitation operator
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr T_pt_(std::size_t K, const OpParams& params = {.order = 1});
ExprPtr Tʼ_(std::size_t K, const OpParams& params = {.order = 1});

/// @brief Makes sum of perturbed particle-conserving excitation operators
/// @param K rank up to which the sum is to be formed
/// @param params OpParams for operator construction. Default: order=1,
/// skip1=false
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr T_pt(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});
ExprPtr Tʼ(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});

/// @brief Makes perturbed particle-conserving deexcitation operator
/// @param K rank of the deexcitation operator
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr Λ_pt_(std::size_t K, const OpParams& params = {.order = 1});
ExprPtr Λʼ_(std::size_t K, const OpParams& params = {.order = 1});

/// @brief Makes sum of perturbed particle-conserving deexcitation operators
/// @param K rank up to which the sum is to be formed
/// @param params OpParams for operator construction. Default: order=1,
/// skip1=false
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr Λ_pt(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});
ExprPtr Λʼ(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});
} // namespace tensor
} // namespace op

Expand Down Expand Up @@ -1234,39 +1234,39 @@ ExprPtr S(std::int64_t K);
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr H_pt(std::size_t R, const OpParams& params = {.order = 1});
ExprPtr Hʼ(std::size_t R, const OpParams& params = {.order = 1});

/// @brief Makes perturbed particle-conserving excitation operator from OpParams
/// @param K rank of the excitation operator
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr T_pt_(std::size_t K, const OpParams& params = {.order = 1});
ExprPtr Tʼ_(std::size_t K, const OpParams& params = {.order = 1});

/// @brief Makes sum of perturbed particle-conserving excitation operators
/// @param K rank up to which the sum is to be formed
/// @param params OpParams for operator construction. Default: order=1,
/// skip1=false
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr T_pt(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});
ExprPtr Tʼ(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});

/// @brief Makes perturbed particle-conserving deexcitation operator
/// @param K rank of the deexcitation operator
/// @param params OpParams for operator construction. Default: order=1
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr Λ_pt_(std::size_t K, const OpParams& params = {.order = 1});
ExprPtr Λʼ_(std::size_t K, const OpParams& params = {.order = 1});

/// @brief Makes sum of perturbed particle-conserving deexcitation operators
/// @param K rank up to which the sum is to be formed
/// @param params OpParams for operator construction. Default: order=1,
/// skip1=false
/// @pre `params.order==1`, only first order perturbation is supported now
/// @pre If batching is used, ISR must contain batching space
ExprPtr Λ_pt(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});
ExprPtr Λʼ(std::size_t K,
const OpParams& params = {.order = 1, .skip1 = false});

/// @brief computes the quantum number change effected by a given Operator or
/// Operator Product when applied to the vacuum state
Expand Down
4 changes: 2 additions & 2 deletions doc/examples/user/cc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,12 +62,12 @@ int main() {

// start-snippet-3
// First-order perturbed amplitude equations
auto t_pt = CC{2}.t_pt(1, 1);
auto t_pt = CC{2}.tʼ(1, 1);
std::wcout << "T1 perturbed: " << to_latex(t_pt[1]) << "\n"
<< "T2 perturbed: " << to_latex(t_pt[2]) << "\n";

// First-order perturbed Lambda amplitude equations
auto l_pt = CC{2}.λ_pt(1, 1);
auto l_pt = CC{2}.λʼ(1, 1);
std::wcout << "λ1 perturbed: " << to_latex(l_pt[1]) << "\n"
<< "λ2 perturbed: " << to_latex(l_pt[2]) << "\n";
// end-snippet-3
Expand Down
4 changes: 2 additions & 2 deletions doc/user/guide/cc.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,8 @@ Coupled-Cluster Response

.. code-block:: cpp

std::vector<ExprPtr> t_pt(size_t order = 1, size_t rank = 1);
std::vector<ExprPtr> λ_pt(size_t order = 1, size_t rank = 1);
std::vector<ExprPtr> tʼ(size_t rank = 1, size_t order = 1);
std::vector<ExprPtr> λʼ(size_t rank = 1, size_t order = 1);

Derives perturbed amplitude equations for response theory calculations.

Expand Down
Loading