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;
IndexSpaceRegistry 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
assert(isr.unIon(L"i", L"a") == isr.retrieve(L"p"));
assert(isr.intersection(L"p", L"i") == isr.retrieve(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_space_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
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_space_registry_shared_ptr = make_min_sr_spaces()});

// set theoretic operations on spaces
auto i1 = Index(L"i_1");
auto a1 = Index(L"a_1");
assert(i1.space().attr().intersection(a1.space().attr()).type() ==
       IndexSpace::Type::null);
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::IndexSpaceRegistry has its own unIon() and intersection() methods that perform error checking to ensure that only registered spaces are defined.

Quasiparticles

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, 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():

using namespace sequant;

// the default is to use genuine vacuum
assert(get_default_context().vacuum() == Vacuum::Physical);
// make default IndexSpaceRegistry
auto isr = std::make_shared<IndexSpaceRegistry>();
// now set the context to a single product of SP states
set_default_context({.index_space_registry_shared_ptr = isr,
                     .vacuum = Vacuum::SingleProduct});
assert(get_default_context().vacuum() == Vacuum::SingleProduct);
// reset the context back to the default
reset_default_context();
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:

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 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/wick.hpp>

int main() {
  using namespace sequant;

  set_default_context(
      {.index_space_registry = IndexSpaceRegistry{}
                                   .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.