Getting started¶
SeQuant is a general-purpose symbolic tensor algebra, but the primary use case is in quantum many-body physics and quantum chemistry. The following is a brief tutorial on using SeQuant for this purpose: applying Wick’s theorem to simple operator products, then giving indices semantic meaning by registering index spaces, and finally composing full many-body equations out of predefined operators.
Wick’s Theorem¶
To get started let’s use SeQuant to apply Wick’s theorem to a simple product of elementary (creation and annihilation) fermionic operators:
This is achieved by the following SeQuant program:
Index p1(L"p_1"), p2(L"p_2"), p3(L"p_3"), p4(L"p_4");
auto cp1 = fcrex(p1), cp2 = fcrex(p2);
auto ap3 = fannx(p3), ap4 = fannx(p4);
std::wcout << to_latex(ap3 * ap4 * cp1 * cp2) << " = "
<< to_latex(FWickTheorem{ap3 * ap4 * cp1 * cp2}
.full_contractions(false)
.compute())
<< std::endl;
Hint
All core user-facing SeQuant code lives in C++ namespace sequant, and the shown code assumes this namespace has been imported via
using namespace sequant.
Running this program should produce a LaTeX expression for this formula:
Wick’s theorem reduces this to a sum of terms in which the elementary operators have been reordered so that every creation operator sits to the left of every annihilation operator — the normal order, or normal-ordered form — plus terms in which some operator pairs have been contracted into overlaps instead. Normal-ordered operators are useful because their vacuum expectation value vanishes trivially, which is why SeQuant expresses Wick’s-theorem output in terms of them. The following tensor notation is used to denote elementary and composite normal-ordered (or, shortly, normal) operators:
\(s^p_q \equiv \langle q | p \rangle\) denotes inner products (“overlaps”) of 1-particle states. Wick’s theorem can of course be applied directly to products of normal composite operators, e.g,
auto nop1 = ex<FNOperator>(cre(p1, p2), ann(p3, p4));
// OR
// auto nop1 = ex<FNOperator>(cre({p1, p2}), ann({p3, p4}));
// OR
// auto nop1 = ex<FNOperator>(cre(std::vector{p1, p2}),
// ann(std::array{L"p3", L"p4"}));
// OR
// auto nop1 = ex<FNOperator>(cre(std::set{"p1", "p2"}),
// ann(std::vector{L"p3", L"p4"}));
auto nop2 = ex<FNOperator>(cre({L"p_5"}), ann({L"p_6", L"p_7"}));
std::wcout
<< to_latex(nop1 * nop2) << " = "
<< to_latex(FWickTheorem{nop1 * nop2}.full_contractions(false).compute())
<< std::endl;
produces
A number-nonconserving operator has an unequal number of upper (creation) and lower (annihilation) indices; \(␣\) marks the resulting empty “slots” that would otherwise misalign the tensor’s rows.
Same algebra can be performed for bosons:
auto nop3 = ex<BNOperator>(cre({p1, p2}), ann({p3, p4}));
auto nop4 = ex<BNOperator>(cre({L"p_5", L"p_6"}), ann({L"p_7"}));
std::wcout
<< to_latex(nop3 * nop4) << " = "
<< to_latex(BWickTheorem{nop3 * nop4}.full_contractions(false).compute())
<< std::endl;
where \(b\) denotes normal bosonic operators constructed analogously with the normal fermionic operators \(a\).
Extended Wick’s theorem¶
With Vacuum::MultiProduct the reference is a general (e.g. multiconfigurational) state that commutes with the number operator,
and operators, which need not conserve particle number, are normal-ordered relative to it (generalized normal order). The same
WickTheorem then evaluates products of such operators into one-body densities \(\gamma\), one-body hole densities
\(\eta\) and density cumulants \(\kappa_k\) over the active (partially occupied) orbitals; max_cumulant_rank bounds the
rank of the cumulants and eta_as_delta_minus_gamma spells \(\eta\) as \(\delta - \gamma\):
// a multireference vocabulary of index spaces: core (o, i), active (u),
// virtual (a, g); the reference is a general state in the active space
set_default_context(
Context({.index_basis_registry_shared_ptr = mbpt::make_mr_spaces(),
.vacuum = Vacuum::MultiProduct}));
// a product of two generalized-normal-ordered one-body operators
auto expr = ex<FNOperator>(cre({L"u_1"}), ann({L"u_2"})) *
ex<FNOperator>(cre({L"u_3"}), ann({L"u_4"}));
// its reference expectation value: a γη pair plus a 2-body cumulant κ
auto vev = FWickTheorem{expr}.compute();
std::wcout << to_latex(vev) << std::endl;
// its generalized-normal-ordered form, cumulants truncated at rank 2
auto gno = FWickTheorem{expr}
.full_contractions(false)
.max_cumulant_rank(2)
.compute();
std::wcout << to_latex(gno) << std::endl;
This prints the reference expectation value
followed by the generalized-normal-ordered form of the product, in which \(\tilde{a}\) denotes operators normal-ordered relative to the reference:
The labels γ, η and κ are reserved: build such tensors with the factories in SeQuant/core/density.hpp (or
parse them, e.g. γ{u_1;u_2}), which give them their fixed symmetries. Note that parsing reads the context in effect: a
parsed ã takes its vacuum, and a parsed density takes the bra-ket symmetry its field implies, so parse an expression under
the context it is meant for.
Register Index Spaces¶
Tensor expressions annotated by abstract indices are common. In some contexts all tensor modes refer to the same range or underlying vector space (as in all examples shown so far); then there is no need to distinguish modes of different types. But in some contexts indices carry important semantic meaning. For example, the energy expression in the coupled-cluster method,
contains tensors with 2 types of modes, denoted by \(i\) and \(a\), that represent single-particle (SP) states occupied and unoccupied in the reference state, respectively. To simplify symbolic manipulation of such expressions SeQuant allows to define a custom vocabulary of index spaces and to define their set-theoretic relationships. The following example illustrates the full space denoted by \(p\) partitioned into occupied \(i\) and unoccupied \(a\) base subspaces:
using namespace sequant;
IndexBasisRegistry isr;
// base spaces
isr.add(L"i", 0b01).add(L"a", 0b10);
// union of 2 base spaces
// can create manually, as isr.add(L"p", 0b11) , or explicitly ...
isr.add_union(L"p", {L"i", L"a"}); // union of i and a
// can access unions and intersections of base and composite spaces
SEQUANT_ASSERT(isr.unIon(L"i", L"a") == isr.retrieve_space(L"p"));
SEQUANT_ASSERT(isr.intersection(L"p", L"i") == isr.retrieve_space(L"i"));
// to use the vocabulary defined by isr use it to make a Context object and
// make it the default
set_default_context({.index_basis_registry = std::move(isr)});
// now can use space labels to construct Index objects representing said
// spaces
Index i1(L"i_1");
Index a1(L"a_1");
Index p1(L"p_1");
// set theoretic operations on spaces
SEQUANT_ASSERT(i1.space().type().includes(a1.space().type()) == false);
This and other vocabularies commonly used in quantum many-body context are supported out-of-the-box by SeQuant; their definitions are in SeQuant/domain/mbpt/convention.hpp. The previous example is equivalent to the following:
using namespace sequant;
using namespace sequant::mbpt;
// makes 2 base spaces, i and a, and their union
set_default_context(
{.index_basis_registry_shared_ptr = make_min_sr_spaces()});
// set theoretic operations on spaces
auto i1 = Index(L"i_1");
auto a1 = Index(L"a_1");
SEQUANT_ASSERT(i1.space().attr().intersection(a1.space().attr()).type() ==
IndexSpace::Type::null);
SEQUANT_ASSERT(i1.space().attr().intersection(a1.space().attr()).qns() ==
mbpt::Spin::any);
Bitset representation of index spaces allows to define set-theoretic operations naturally. Bitset-based representation is used not only for index space type attribute (IndexSpace::Type) but also for the quantum numbers attribute (IndexSpace::QuantumNumbers). The latter can be used to represent spin quantum numbers, particle types, etc.
The main difference of the last example with the original example is that the make_min_sr_spaces() factory changes the quantum numbers used by default (mbpt::Spin::any) to make spin algebraic manipulations (like tracing out spin degrees of freedom) easier. Users can create their own definitions to suit their needs, but the vast majority of users will not need to venture outside of the predefined vocabularies.
Notice that the set-theoretic operations are only partially automated. It is the user’s responsibility to define any and all unions and intersections of base spaces that they may encounter in their context. For this reason sequant::IndexBasisRegistry has its own add_unIon() and add_intersection() methods that perform error checking to ensure that only registered spaces are defined.
Named basis instances¶
An sequant::Index can run over a specific basis instance of its space (see Context and Configuration). A
basis instance can be registered under a label of its own with sequant::IndexBasisRegistry::add(); the entry
carries its own extent (the number of functions in the basis), metric (whether it is orthonormal) and field, read
off an index by sequant::IndexBasis::extent(), sequant::IndexBasis::metric() and
sequant::IndexBasis::field(). That label, like a space’s, must be one that an index can be parsed with
(sequant::io::serialization::v1::is_base_key()); registration throws otherwise. The space views
(sequant::IndexBasisRegistry::spaces()) see spaces only, and sequant::IndexBasisRegistry::retrieve_space() throws
on such a label, use sequant::IndexBasisRegistry::retrieve(). The set algebra and the occupancy queries take
the label of a named basis as that of its space; the setters that designate a space’s role
(sequant::IndexBasisRegistry::vacuum_occupied_space() and its kin) set a type, so they reject it as
sequant::IndexBasisRegistry::retrieve_space() does. Indices in a named basis are printed, serialized and deserialized
by that label, and constructed from it; the name is part of the basis’s identity, and code that keys an axis by
label tells them apart from the space’s own basis by sequant::IndexBasis::base_key(). An unnamed basis instance has no entry and no metadata of its own: it
is the space’s own basis up to an orthonormal rotation, so its extent is the dimension of the space, its metric is
unit and its field is the space’s, and its basis key is its space’s, so evaluation sizes, tiles and slices it as the
space. That is right for a rotation of the space, such as a localized basis; a basis of a different extent, such as a
truncated or an overcomplete set, or a non-orthonormal one must be registered under a name. A named basis can also be
made to follow another named basis (sequant::IndexBasisRegistry::follow()): its extent, metric and field are then
those of the basis it follows, set through that one’s label, as the PAO bases follow the AO basis they are projected
from.
using namespace sequant;
auto isr = mbpt::make_min_sr_spaces();
// a basis instance of a registered space can be registered under its own
// label: here a localized basis of the space i, as its basis instance 1
const IndexSpace i = isr->retrieve_space(L"i");
isr->add(L"ĩ", IndexBasis{i, 1});
// a Context owns its registry: register everything first, then hand it over
// (it is moved in)
set_default_context({.index_basis_registry_shared_ptr = std::move(isr)});
const auto& registry = *get_default_context().index_basis_registry();
Index loc1(registry.retrieve(L"ĩ"), 1);
// same space ...
SEQUANT_ASSERT(loc1.space() == i);
// ... a specific basis of it, which carries its name: the name is part of
// the basis, as a space's label is of the space, so the registry's entry
// is the basis and a bare {space, instance} pair is another, unnamed one
SEQUANT_ASSERT(loc1.basis() == registry.retrieve(L"ĩ"));
SEQUANT_ASSERT(loc1.basis() != IndexBasis(i, 1));
SEQUANT_ASSERT(same_instance(loc1.basis(), IndexBasis(i, 1)));
SEQUANT_ASSERT(registry.resolve(IndexBasis(i, 1)) == loc1.basis());
// a name is not a space
SEQUANT_ASSERT(registry.contains(L"ĩ") &&
registry.retrieve_space_ptr(L"ĩ") == nullptr);
// the space views see i once: the named entry is not a space
SEQUANT_ASSERT(
std::ranges::count_if(registry.spaces(), [](const IndexSpace& s) {
return s.base_key() == L"i";
}) == 1);
// indices in a named basis are printed, serialized and deserialized by that
// label, and constructed from it
SEQUANT_ASSERT(Index(L"ĩ_1") == loc1);
SEQUANT_ASSERT(loc1.full_label() == L"ĩ_1"); // not i_1<;1>
SEQUANT_ASSERT(loc1.basis().base_key() == L"ĩ");
SEQUANT_ASSERT(serialize(deserialize(L"t{ĩ_1;i_1}"), {.annot_symm = false}) ==
L"t{ĩ_1;i_1}");
Quasiparticles¶
Index spaces interact with another important piece of context: the choice of reference vacuum state. In most cases we are interested in using SeQuant to manipulate expressions involving operators in normal order relative to a vacuum state with a finite number of particles — whose excitations are called quasiparticles — rather than with respect to the genuine vacuum with zero particles. The choice of vacuum state as well as other related traits (whether the SP states are orthonormal, etc.) is defined by the implicit global context. The SeQuant programs until now used the genuine vacuum. The active context can be examined by calling get_default_context(), changed via set_default_context(), and reset to the default via reset_default_context() (see the User Guide’s Context and Configuration for the full picture, including the one-line mbpt::load() shortcut for standard quantum-chemistry conventions):
using namespace sequant;
// the default is to use genuine vacuum
SEQUANT_ASSERT(get_default_context().vacuum() == Vacuum::Physical);
// now set the context to a single product of SP states
set_default_context({.index_basis_registry = IndexBasisRegistry{},
.vacuum = Vacuum::SingleProduct});
SEQUANT_ASSERT(get_default_context().vacuum() == Vacuum::SingleProduct);
// reset the context back to the default
reset_default_context();
SEQUANT_ASSERT(get_default_context().vacuum() == Vacuum::Physical);
However, to deal with the single-product vacuum it is necessary to register at least one space and announce it as occupied in the vacuum state:
IndexBasisRegistry isr;
isr.add(L"y", 0b01, is_vacuum_occupied);
It is also necessary to specify the complete space (union of all base spaces) so that the space of unoccupied SP states can be determined:
isr.add(L"y", 0b01, is_vacuum_occupied)
.add(L"z", 0b10)
.add(L"p", 0b11, is_complete);
The Wick’s theorem code itself is independent of the choice of vacuum:
#include <SeQuant/core/context.hpp>
#include <SeQuant/core/io/shorthands.hpp>
#include <SeQuant/core/wick.hpp>
int main() {
using namespace sequant;
set_default_context(
{.index_basis_registry = IndexBasisRegistry{}
.add(L"y", 0b01, is_vacuum_occupied)
.add(L"z", 0b10)
.add(L"p", 0b11, is_complete),
.vacuum = Vacuum::SingleProduct});
auto cp1 = fcrex(L"p_1"), cp2 = fcrex(L"p_2");
auto ap3 = fannx(L"p_3"), ap4 = fannx(L"p_4");
std::wcout << to_latex(ap3 * cp1 * ap4 * cp2) << " = "
<< to_latex(FWickTheorem{ap3 * cp1 * ap4 * cp2}
.full_contractions(false)
.compute())
<< std::endl;
return 0;
}
produces
Note that:
the tilde in \(\tilde{a}\) denotes normal order with respect to single-product vacuum
Einstein summation convention is implied, i.e., indices that appear twice in a given product (once in superscript, once in subscript) are summed over.
Operators¶
Development of SeQuant is primarily motivated by the perturbative many-body methods, collectively referred to here as Many-Body Perturbation Theory (MBPT). Examples of such methods include the coupled-cluster (CC) method and GW. The typical use case is to compute canonical forms of products of operators. For example, consider the coupled-cluster doubles (CCD) method. Amplitudes \(t^{i_1 i_2}_ {a_1 a_2}\) of the cluster operator,
are determined by solving the CCD equations:
A pedestrian way to compose such an expression is to define a cluster operator object using SeQuant tensors and normal-ordered operators. See the User Guide’s The Expression Tree for the full vocabulary of expression nodes (tensors, products, sums, and more) used below.
auto t2 = ex<Constant>(rational(1, 4)) *
ex<Tensor>(L"t", bra{L"a_1", L"a_2"}, ket{L"i_1", L"i_2"},
Symmetry::Antisymm) *
ex<FNOperator>(cre{L"a_1", L"a_2"}, ann{L"i_1", L"i_2"});
The normal-ordered Hamiltonian is defined similarly as a sum of 1- and 2-body contributions:
auto H = ex<Tensor>(L"f", bra{L"p_1"}, ket{L"p_2"}, Symmetry::Nonsymm) *
ex<FNOperator>(cre{L"p_1"}, ann{L"p_2"}) +
ex<Constant>(rational(1, 4)) *
ex<Tensor>(L"g", bra{L"p_1", L"p_2"}, ket{L"p_3", L"p_4"},
Symmetry::Antisymm) *
ex<FNOperator>(cre{L"p_1", L"p_2"}, ann{L"p_3", L"p_4"});
Note that the compact definition of the Hamiltonian is due to the use of the union (\(p\)) of base occupied (\(i\)) and unoccupied (\(a\)) spaces. Many other symbolic algebras only support use of nonoverlapping base spaces, in terms of which Hamiltonian and other tensor expressions would have a much more verbose form. The commutator of the Hamiltonian and cluster operators is trivially composed:
auto commutator = [](auto op1, auto op2) {
return simplify(op1 * op2 - op2 * op1);
};
auto c_ht = commutator(H, t2);
Note the use of simplify to rewrite an expression in a simpler form; see the User Guide’s Canonicalization for why this is needed and how it works.
Unfortunately, we immediately run into the limitation of the “pedestrian” approach. Namely, the double commutator cannot be correctly obtained as
auto c_htt_v1 =
ex<Constant>(rational(1, 2)) * commutator(commutator(H, t2), t2);
due to the explicit use of specific dummy indices in the definition of t2. Using it more than once in a given product will produce an expression where each dummy index appears more than 2 times, breaking the Einstein summation convention.
The issue is actually not the reuse of the same of operator object, but the reuse of dummy indices. A straightforward, but brittle, solution is to ensure that each dummy index is only used once. E.g., to use \(\hat{t}_2\) more than once in an expression we must make several versions of it, each with a separate set of dummy indices:
auto t2_0 = ex<Constant>(rational(1, 4)) *
ex<Tensor>(L"t", bra{L"a_1", L"a_2"}, ket{L"i_1", L"i_2"},
Symmetry::Antisymm) *
ex<FNOperator>(cre{L"a_1", L"a_2"}, ann{L"i_1", L"i_2"});
auto t2_1 = ex<Constant>(rational(1, 4)) *
ex<Tensor>(L"t", bra{L"a_3", L"a_4"}, ket{L"i_3", L"i_4"},
Symmetry::Antisymm) *
ex<FNOperator>(cre{L"a_3", L"a_4"}, ann{L"i_3", L"i_4"});
auto c_htt =
ex<Constant>(rational(1, 4)) * commutator(commutator(H, t2_0), t2_1);
This is too error-prone for making complex expressions. A better way is to represent \(\hat{t}_2\) by an object that generates tensor form with unique dummy indices generated on the fly. In SeQuant such MBPT operators live in mbpt namespace. The entire CCD amplitude equation is evaluated as follows:
#include <SeQuant/core/context.hpp>
#include <SeQuant/core/expr.hpp>
#include <SeQuant/core/io/shorthands.hpp>
#include <SeQuant/core/op.hpp>
#include <SeQuant/domain/mbpt/context.hpp>
#include <SeQuant/domain/mbpt/convention.hpp>
#include <SeQuant/domain/mbpt/op.hpp>
#include <SeQuant/domain/mbpt/vac_av.hpp>
inline auto commutator(auto op1, auto op2) { return op1 * op2 - op2 * op1; }
int main() {
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()});
auto hbar =
H(2) + commutator(H(2), t(2)) +
ex<Constant>(rational(1, 2)) * commutator(commutator(H(2), t(2)), t(2));
auto ccd_eq = vac_av(P(2) * hbar);
std::wcout << "<" << to_latex(P(2) * hbar) << "> = " << to_latex(ccd_eq)
<< std::endl;
return 0;
}
The result is
The use of MBPT operators rather than their tensor-level forms not only solves problems with the reuse of dummy indices, but also allows to implement additional optimizations such as algebraic simplifications of complex operator expressions and avoiding evaluation of operator products whose vacuum expectation values are guaranteed to vanish. This allows very efficient derivation of complex equations, e.g. CC equations through CCSDTQ are derived in a fraction of a second on a laptop:
$ cmake --build build --target srcc
$ time build/srcc time ./srcc 4 t std so
CC equations [type=t,rank=1,spinfree=false,screen=true,use_topology=true,use_connectivity=true,canonical_only=true] computed in 0.0027805 seconds
R1(expS1) has 8 terms:
CC equations [type=t,rank=2,spinfree=false,screen=true,use_topology=true,use_connectivity=true,canonical_only=true] computed in 0.012890749999999999 seconds
R1(expS2) has 14 terms:
R2(expS2) has 31 terms:
CC equations [type=t,rank=3,spinfree=false,screen=true,use_topology=true,use_connectivity=true,canonical_only=true] computed in 0.039590500000000001 seconds
R1(expS3) has 15 terms:
R2(expS3) has 37 terms:
R3(expS3) has 47 terms:
CC equations [type=t,rank=4,spinfree=false,screen=true,use_topology=true,use_connectivity=true,canonical_only=true] computed in 0.107501417 seconds
R1(expS4) has 15 terms:
R2(expS4) has 38 terms:
R3(expS4) has 53 terms:
R4(expS4) has 74 terms:
./srcc 4 t std so 0.27s user 0.41s system 100% cpu 0.674 total