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,

\[E_\mathrm{CC} = F^{a_1}_ {i_1} t^{i_1}_ {a_1} + \frac{1}{4} \bar{g}^{a_1 a_2}_ {i_1 i_2} (t^{i_1 i_2}_ {a_1 a_2} + 2 t^{i_1}_ {a_1} t^{i_2}_ {a_2})\]

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

\[{{\tilde{a}_ {{p_3}}}{\tilde{a}^{{p_1}}}{\tilde{a}_ {{p_4}}}{\tilde{a}^{{p_2}}}} = \bigl({\tilde{a}^{{p_1}{p_2}}_ {{p_3}{p_4}}} - {{s^{{p_1}}_ {{z_1}}}{s^{{z_1}}_ {{p_3}}}{\tilde{a}^{{p_2}}_ {{p_4}}}} - {{s^{{p_2}}_ {{z_1}}}{s^{{z_1}}_ {{p_4}}}{\tilde{a}^{{p_1}}_ {{p_3}}}} - {{s^{{p_1}}_ {{y_1}}}{s^{{y_1}}_ {{p_4}}}{\tilde{a}^{{p_2}}_ {{p_3}}}} + {{s^{{p_2}}_ {{z_1}}}{s^{{z_1}}_ {{p_3}}}{\tilde{a}^{{p_1}}_ {{p_4}}}} + {{s^{{p_1}}_ {{z_1}}}{s^{{p_2}}_ {{z_2}}}{s^{{z_1}}_ {{p_3}}}{s^{{z_2}}_ {{p_4}}}} + {{s^{{p_1}}_ {{y_1}}}{s^{{p_2}}_ {{z_1}}}{s^{{z_1}}_ {{p_3}}}{s^{{y_1}}_ {{p_4}}}}\bigr)\]

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.