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 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. The typical use case is to compute canonical forms of products of operators. For example, consider the coupled-cluster doubles (CCD) method. Amplitudes $𝑡 𝑖 1 𝑖 2 𝑎 1 𝑎 2$ of the cluster operator,

ˆ 𝑡 \equiv ˆ 𝑡 2 = 14 𝑡 𝑖 1 𝑖 2 𝑎 1 𝑎 2 𝑎 𝑎 1 𝑎 2 𝑖 1 𝑖 2,

are determined by solving the CCD equations:

\forall 𝑖 1, 𝑖 2, 𝑎 1, 𝑎 2 : 0 = ⟨ 0 | 𝑎 𝑖 1 𝑖 2 𝑎 1 𝑎 2 e x p (- ˆ 𝑡 2) ˆ 𝐻 e x p (ˆ 𝑡 2) | 0 ⟩ = ⟨ 0 | 𝑎 𝑖 1 𝑖 2 𝑎 1 𝑎 2 (ˆ 𝐻 + [ˆ 𝐻, ˆ 𝑡 2] + 12 [[ˆ 𝐻, ˆ 𝑡 2], ˆ 𝑡 2]) | 0 ⟩ .

A pedestrian way to compose such expression is to define a cluster operator object using SeQuant tensors and normal-ordered operators:

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 ( $𝑝$ ) of base occupied ( $𝑖$ ) and unoccupied ( $𝑎$ ) 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. Commutator of the Hamiltonian and cluster operator 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. Its role will be emphasized later. 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 expresson 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 $ˆ 𝑡 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 $ˆ 𝑡 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/op.hpp>
#include <SeQuant/core/tensor.hpp>
#include <SeQuant/domain/mbpt/convention.hpp>
#include <SeQuant/domain/mbpt/op.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(Context(make_min_sr_spaces(), Vacuum::SingleProduct));

  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

(14 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑖 1 𝑖 2 𝑎 1 𝑎 2 - 12 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 𝑓 𝑎 3 𝑎 1 ¯ 𝑡 𝑖 1 𝑖 2 𝑎 2 𝑎 3 - 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑖 1 𝑎 3 𝑖 3 𝑎 1 ¯ 𝑡 𝑖 2 𝑖 3 𝑎 2 𝑎 3 + 18 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑎 1 𝑎 2 ¯ 𝑡 𝑖 1 𝑖 2 𝑎 3 𝑎 4 + 18 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑖 1 𝑖 2 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 3 𝑖 4 𝑎 1 𝑎 2 + 12 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 𝑓 𝑖 1 𝑖 3 ¯ 𝑡 𝑖 2 𝑖 3 𝑎 1 𝑎 2 + 116 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 2 𝑎 3 𝑎 4 ¯ 𝑡 𝑖 3 𝑖 4 𝑎 1 𝑎 2 - 14 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 3 𝑎 3 𝑎 4 ¯ 𝑡 𝑖 2 𝑖 4 𝑎 1 𝑎 2 - 14 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 2 𝑎 1 𝑎 3 ¯ 𝑡 𝑖 3 𝑖 4 𝑎 2 𝑎 4 - 14 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 3 𝑎 3 𝑎 4 ¯ 𝑡 𝑖 2 𝑖 4 𝑎 1 𝑎 2 - 14 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 2 𝑎 1 𝑎 3 ¯ 𝑡 𝑖 3 𝑖 4 𝑎 2 𝑎 4 + 12 𝐴 𝑎 1 𝑎 2 𝑖 1 𝑖 2 ¯ 𝑔 𝑎 3 𝑎 4 𝑖 3 𝑖 4 ¯ 𝑡 𝑖 1 𝑖 3 𝑎 1 𝑎 3 ¯ 𝑡 𝑖 2 𝑖 4 𝑎 2 𝑎 4)

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