Operator Interface for Symbolic Many-Body Algebra¶

Motivation¶

The mbpt::Operator class represents an abstract Fock-space many-body operator, such as the Hamiltonian operator, that can be used within SeQuant expressions. We already saw earlier that the lower-level tensor representation of such operators involves specific index labels which makes tensor representation inconvenient when multiple instances of the same operator occur. For example, consider a general 1-body operator

\[\hat{o}_1 \equiv o^{p_1}_{p_2} a_{p_2}^{p_1}.\]

The abstract representation \(\hat{o}_1\) is convenient for efficient algebraic manipulations such as expansion and simplication:

\[\left(1 + \hat{o}_1\right)^2 = 1 + \hat{o}_1 + \hat{o}_1 + \hat{o}_1^2 = 1 + 2 \hat{o}_1 + \hat{o}_1^2\]

The tensor form is less convenient for such manipulations for 2 reasons:

lowering multiple instances of \(\hat{o}_1\) such as in its square \(\hat{o}_1^2 \equiv \hat{o}_1 \hat{o}_1\) to the tensor form is context-sensitive since one must avoid index reuse:

\[\hat{o}_1^2 = o^{p_1}_ {p_2} a_{p_2}^{p_1} o^{p_3}_{p_4} a_{p_4}^{p_3}\]
operating on tensor forms of many-body operators involves frequent use of canonicalization, such as to be able to simplify \(\hat{o}_1 + \hat{o}_1\) in the following tensor form

\[\hat{o}_1 + \hat{o}_1 = o^{p_1}_ {p_2} a_{p_2}^{p_1} + o^{p_3}_{p_4} a_{p_4}^{p_3}\]

to

\[2 \hat{o}_1 = 2 o^{p_1}_ {p_2} a_{p_2}^{p_1}\]

requires canonicalization of each term before attempting simplification. In the abstract form simplifying \(\hat{o}_1 + \hat{o}_1\) to \(2 \hat{o}_1\) is as simple as if \(\hat{o}_1\) were a scalar.

To be useful for more than just algebraic simplification the abstract operator form must include enough detail to determine what products of such operators do to quantum numbers (e.g., occupation numbers) of a Fock-space state. Each mbpt::Operator thus denotes its effect on the quantum numbers. This allows efficient evaluation of vacuum averages of products of mbpt::Operator objects by identifying only the non-zero contributions.

`mbpt::Operator` Anatomy¶

mbpt::Operator is constructed from 3 callables (e.g., lambdas):

void -> std::wstring_view: returns the operator label such as “t” or “λ”
void -> ExptPtr: returns the tensor form of the operator, and
QuantumNumberChange<>& -> void: encodes the action on the given input state by mutating its quantum numbers given as the argument to the lambda.

The quantum numbers of a state and their changes are represented by the same class, mbpt::QuantumNumberChange. For the Fock-space states the quantum numbers are the number of creators/annihilators. In general context the number of creators/annihilators must be tracked separately for multiple IndexSpaces (e.g., for each hole and particle states in the particle-hole context). To support general operators which can create a variable number of particles and holes the numbers of creators/annihilators must be encoded by intervals. For example, if space \(\{p\}\) is a union of hole (fully occupied in the vacuum) space \(\{i\}\) and particle (fully empty in the vacuum) space \(\{a\}\) operator \(\hat{o}_1\) includes 4 components,

\[o^{p_1}_{p_2} a_{p_2}^{p_1} = o^{i_1}_{i_2} a_{i_2}^{i_1} + o^{i_1}_{a_2} a_{a_2}^{i_1} + o^{a_1}_{i_2} a_{i_2}^{a_1} + o^{a_1}_{a_2} a_{a_2}^{a_1},\]

each of which has 1 (particle or hole) creator and 1 (particle or hole) annihilator. The net effect of this operator on a Fock-space state is to leave its number of hole/particle creators/annihilators unchanged or increased by 1. E.g., its effect on a state with \(n\) particle creators is produce a state with \([n,n+1]\) particle creators; similarly, action on a state with \([n,m]\) hole annihilators (\(n \leq m\)) will produce a state with \([n,m+1]\) hole annihilators.

Finally, for a product of mbpt::Operator it is possible to determine its effect on the quantum numbers of the input state using simple rules based on Wick’s theorem. This enables SeQuant to screen out expressions that would yield zero, before applying Wick’s theorem without much computational overhead. This logic can also be used to check if two operators commute with each other, which is crucial for reducing the expression to its compact canonical form.

The use of interval arithmetic for quantum number changes allows to treat operators in their natural compact form as long as possible. E.g., \(\hat{o}_1\) can be kept as a single term \(o^{p_1}_{p_2} a_{p_2}^{p_1}\) rather than having to decompose it into individual contributions, each with definite effect on quantum numbers. This is especially useful when the number of index spaces increases.

Operator Registry¶

The mbpt::OpRegistry maintains a mapping from operator labels to their operator classes. In this section we describe how to use and extend the operator registry.

Each operator must be registered with one of three classes:

OpClass::Ex: Excitation operators (e.g., cluster excitation operator)
OpClass::Deex: De-excitation operators (e.g., cluster de-excitation operator)
OpClass::Gen: General operators (e.g., Hamiltonian-like operators)

Note

Some operator labels are reserved for SeQuant’s internal MBPT and core functionalities and cannot be registered. See SeQuant/core/reserved.hpp for more details.

SeQuant provides two predefined registries:

make_minimal_registry(): Contains only a minimal set of operators
make_legacy_registry(): A backward-compatible registry that includes all operators used in legacy MBPT module.

Custom operators can be registered to extend SeQuant’s vocabulary. Once registered, operators can be used with both OpMaker and mbpt::Operator:

// Create a custom operator registry with a new operator
auto custom_registry = std::make_shared<OpRegistry>();

custom_registry->add(L"f", OpClass::Gen)
    .add(L"g", OpClass::Gen)
    .add(L"t", OpClass::Ex)
    .add(L"λ", OpClass::Deex)
    .add(L"X", OpClass::Ex);  // Custom excitation operator

// Use the custom registry to set the default MBPT context
set_default_mbpt_context({.op_registry_ptr = custom_registry});

// Create custom operator using OpMaker
auto X2_tensor = OpMaker<Statistics::FermiDirac>(L"X", 2)();

// Create custom operator using mbpt::Operator
auto X2_op = ex<op_t>(
    []() -> std::wstring_view { return L"X"; },
    []() -> ExprPtr { return OpMaker<Statistics::FermiDirac>(L"X", 2)(); },
    [](qns_t& qns) { qns += excitation_type_qns(2); });

// Use in expressions
auto custom_expr = h(2) * X2_op * t(2);

A registered (de)excitation operator can also be granted a basis instance per leg space with grant_basis(), which OpMaker then gives its legs in that space (see Context and Configuration); basis_grant() looks a grant up. General operators cannot be granted one.

The registry is stored in the mbpt::Context. See the unit test cases for further manipulations with mbpt::Context and mbpt::OpRegistry.

Normalization conventions¶

sequant::mbpt::NormalizationConvention in the MBPT context controls prefactors of coefficient-summed operators. For an operator with \(c\) creators and \(a\) annihilators, let \(m=c!a!\) in a spin-orbital basis; in a spin-free basis, \(c=a\) and \(m=c!\).

Default gives coefficient-summed operators the prefactor \(1/m\), while projectors A, S, and P have no additional prefactor.
Symmetric gives (de)excitation operators (OpClass::Ex and OpClass::Deex), whose coefficients are the amplitudes being solved for, and projectors the prefactor \(1/\sqrt{m}\). General operators (OpClass::Gen), such as the Hamiltonian, perturbations, and F12 geminals, have supplied coefficients and retain \(1/m\).

The prefactor follows from the operator’s sequant::mbpt::OpClass, so register a custom operator whose coefficients are solved for as Ex or Deex. Projectors are recognized by the reserved labels reserved::antisymm_label() (Â) and reserved::symm_label() (Ŝ): only operators with these labels receive the projector prefactor, and any other operator is normalized by its class. The replacement operator ã has explicit free indices and no summed coefficient tensor, so it carries no factorial prefactor.

Normalization is selected when tensor forms are generated; abstract operator identity is independent of the convention. Keep a scoped convention active through lowering to tensor form. See Context and Configuration for a compiled example of a scoped convention.

Lie Similarity Transformation¶

ExprPtr lst(ExprPtr A, ExprPtr B, size_t r, LSTOptions options = {});

Returns Lie similarity transformation, \(\bar{A} = e^{-\hat{B}} \hat{A} e^{\hat{B}}\), as a series of nested commutators, \([\hat{A},\hat{B}]\), \([[\hat{A},\hat{B}],\hat{B}]\), etc., up to order r.

The method lst is declared in <SeQuant/domain/mbpt/utils.hpp>. The CC class uses it to construct the similarity-transformed Hamiltonian, but it is a standalone utility and can be used on any expression built from mbpt::Operator.

The behavior of lst is controlled by LSTOptions:

struct LSTOptions {
  SEQUANT_DESIGNATED_INIT_ONLY;  // see below
  bool unitary = false;
  bool use_connected_form = false;
  bool skip_clone = false;
};

LSTOptions must be written with C++20 designated initializers – LSTOptions{.unitary = true}, or {} for all defaults. Positional initialization, LSTOptions{true, false}, is rejected at compile time by the SEQUANT_DESIGNATED_INIT_ONLY member. This is deliberate: fields have been reordered and repurposed before (use_connected_form took the position of the older use_commutators and inverted its meaning), and a positional call site would have kept compiling with its meaning silently changed. Naming the fields makes any such change a compile error instead.

unitary: if true, \(\hat{B}\) is replaced by \(\hat{B} - \hat{B}^\dagger\), i.e. the method returns \(e^{-(\hat{B} - \hat{B}^\dagger)} \hat{A} e^{\hat{B} - \hat{B}^\dagger}\). Use this for unitary ansätze.
use_connected_form: selects how each commutator is written, see below.
skip_clone: if true, A is transformed in place instead of being cloned first. Set it only if nothing else holds a reference to A.

Commutators or connected products¶

Each commutator can be written in two ways, and use_connected_form selects between them:

false (the default) writes it explicitly, \([\hat{A},\hat{B}] = \hat{A}\hat{B} - \hat{B}\hat{A}\).
true writes it as a connected product, \((\hat{A}\hat{B})_c\). This gives fewer terms, but only reproduces the commutator once the same operators are connected downstream when taking the expectation value, using OpConnections.

The connectivity discussion below assumes the reference is the Wick vacuum. Otherwise use explicit commutators: the reference average of \(\hat{B}\hat{A}\) need not vanish, so the connected product need not reproduce the commutator under any connectivity. Unless the vacuum is Vacuum::MultiProduct, ref_av then also rejects connectivity lists (see Vacuum averaging and final expression), even when assertions are disabled.

Both forms give the same equations, given the right connectivity. They differ in where the disconnected terms are removed: the commutator removes them algebraically, the connected product relies on the connectivity you supply to vac_av/ref_av.

Both the operator-level and tensor-level vac_av/ref_av overloads default to empty connectivity constraints. Omitting the options argument is equivalent to passing {}. To use the connected-product form, supply the required connections explicitly through EVOptions::connect. default_op_connections() provides conventional CC connections between the Hamiltonian (h, f, f̃, g) and the cluster operator t. Pairing use_connected_form = true with empty connectivity keeps disconnected terms that the commutator would have cancelled, so the result is wrong.

The same trade-off shows up in CC::hbar(), which returns the connected form for a non-unitary ansatz when the reference is the Wick vacuum; see Using \bar{H} outside the CC class.

Examples¶

Let’s look at some examples using the following SeQuant context:

using namespace sequant;
using namespace sequant::mbpt;
set_default_context({.index_basis_registry_shared_ptr = make_min_sr_spaces(),
                     .vacuum = Vacuum::SingleProduct});
set_default_mbpt_context({.op_registry_ptr = make_legacy_registry()});

Constructing an `Operator`¶

using sequant::mbpt::op_t;
using sequant::mbpt::qns_t;

// constructor takes label, tensor form and QuantumNumberChange
op_t fock_op([]() -> std::wstring_view { return L"f"; },
             []() -> ExprPtr {
               return ex<Tensor>(L"f", bra{L"p_1"}, ket{L"p_2"}) *
                      ex<FNOperator>(cre({L"p_1"}), ann({L"p_2"}));
             },
             [](qns_t& qns) { qns += sequant::mbpt::general_type_qns(1); });

The components of the Operator can be accessed as follows:

// get the label
auto label = fock_op.label();

// apply operator to another state
auto vac = qns_t{};  // vacuum state for example
auto modified = fock_op(vac);

// access the tensor form
auto fock_tensor = fock_op.tensor_form();

Expression construction and screening¶

using namespace sequant::mbpt::op;

// Construct a double excitation operator (T2)
auto t2 = t(2);

// Construct a two-body Hamiltonian
auto h2 = h(2);

// Product: h2 * t2 * t2
auto expr1 = h2 * t2 * t2;

// Screening: can an expression raise the vacuum to double excitation?
bool raises_to_double = raises_vacuum_to_rank(expr1, 2);  // true

// Screening: can an expression raise the vacuum up to double excitation?
bool raises_up_to_double = raises_vacuum_up_to_rank(expr1, 2);  // true

// Screening: can an expression lower a double excitation to the vacuum?
bool lowers_to_vacuum = lowers_rank_to_vacuum(λ(2) * h(1), 2);  // true

Vacuum averaging and final expression¶

The sequant::mbpt::op::vac_av function can be used to compute the vacuum average of an operator level expression. If reference state differs from the Wick vacuum sequant::mbpt::op::ref_av function should be used instead to compute the reference average. In this case partial contractions leave residual operators that become RDMs after Wick’s theorem has run, and both overloads of ref_av require empty connect and do_not_connect lists: an RDM of rank above one contains both connected and disconnected contributions, so whether the operators feeding it are connected is not a property of a term. An RDM index is external if it belongs to the context’s named indices (CanonicalizeOptions::named_indices) or, when no named indices are set, if it occurs in a single bra or ket slot of its term. An index in a single bra or ket slot cannot be summed over, so when named indices are set they must include every such index of a term that contains an RDM; ref_av fails otherwise. A dummy RDM index is summed over the active space and is relabeled into it; an external one keeps its label, and if its space extends beyond the active space the RDM takes an active-space index tied to the external one by a Kronecker delta. The active-space index carries no basis instance or protoindices of the index it replaces. An RDM index whose space does not overlap the active space, one lying entirely in the vacuum-occupied space or entirely outside the reference-occupied space, is left as is, although the normal-ordered expectation value it stands for vanishes there. vac_av always computes full contractions and honors these options, but it computes the Wick vacuum average, which may differ from the reference average. Under Vacuum::MultiProduct (see below) ref_av takes full contractions and honors the lists too.

Two operators are connected in a term if a factor that Wick’s theorem produces carries indices of both: a contraction or, under Vacuum::MultiProduct, also a density or a cumulant. Operator-level connection pairs use labels and apply to every matching pair with the first operator to the left of the second. Tensor-level pairs use zero-based normal-operator positions. Both conventions apply to connect and do_not_connect in EVOptions.

Both fermions genuine (Vacuum::Physical), single-determinant (Vacuum::SingleProduct), and multiconfigurational (Vacuum::MultiProduct) states can be used as the vacuum state; for bosons only the genuine vacuum is currently supported. When the vacuum is Vacuum::MultiProduct, ref_av (and vac_av) apply the extended Wick theorem instead, as sequant::WickTheorem does under that vacuum; the operators are then normal-ordered relative to the reference itself (generalized normal order), and the result is expressed in densities and density cumulants.

using namespace sequant::mbpt;

auto expr = op::H(2) * op::T(2) * op::T(2);
// connect the Hamiltonian with every t: H(2) is written in terms of f or h,
// depending on the vacuum, and default_op_connections() covers both
const EVOptions<std::wstring> opts{.connect = default_op_connections()};
auto result = op::vac_av(op::P(2) * expr, opts);
// vac_av is equivalent to ref_av for single-determinant reference:
// auto result = op::ref_av(op::P(2) * expr, opts);

std::wcout << "Result: " << result << "\n";

Note

op::P is a projection operator which can be used to construct an excited bra or ket manifold based on the input.

Operator Interface for Symbolic Many-Body Algebra¶

Motivation¶

mbpt::Operator Anatomy¶

Operator Registry¶

Normalization conventions¶

Lie Similarity Transformation¶

Commutators or connected products¶

Examples¶

Constructing an Operator¶

Expression construction and screening¶

Vacuum averaging and final expression¶

`mbpt::Operator` Anatomy¶

Constructing an `Operator`¶