Second-quantization operators

many_body_operator_generic is a templated class, which implements the algebra of fermionic operators. An object of this class represents a general fermionic operator and supports all standard algebraic operations (sums, products, multiplication by a scalar). It allows to write readable and clean C++ code involving various operators, such as Hamiltonians and observables of many-body systems.


The internal storage of the many_body_operator_generic object is not based on a matrix representation. Instead of that, the object stores a list of normally-ordered monomials in basis elements (creation and annihilation operators), accompanied by scalar coefficients. This approach allows to minimize the required storage space.

many_body_operator_generic is declared as

namespace triqs {
namespace operators {
    template<typename ScalarType> class many_body_operator_generic;

Template parameters

  • ScalarType determines the scalar type to construct the algebra.

ScalarType will be double, std::complex<double>, or triqs::utility::real_or_complex most of the time. triqs::utility::real_or_complex is a variant numeric type that can carry either a real or a complex value. Type of the result of an arithmetic expression involving real_or_complex objects is dynamically deduced at runtime.

One may, however, use an arbitrary user-defined type, as long as it meets two conditions:

  • Objects of the type form a field, i.e. they support operations +, -, *, = and compound assignments +=, -=, *=, /=.
  • There is a free function bool is_zero(ScalarType [const&] x), declared in namespace triqs::utility, which detects if an object of the type ScalarType is zero.

One more requirement must be met to make free function dagger() (Hermitian conjugate) available:

  • There is a free function ScalarType conj(ScalarType [const&] x), declared in namespace triqs::utility, which returns a complex conjugate of x.

For the sake of convenience, three type aliases are declared:

// Operator with real matrix elements
using many_body_operator_real = many_body_operator_generic<double>;
// Operator with complex matrix elements
using many_body_operator_complex = many_body_operator_generic<std::complex<double>>;
// Operator with polymorphic matrix elements
using many_body_operator = many_body_operator_generic<real_or_complex>;


many_body_operator_generic provides a minimal set of constructors:

// Default constructor; constructs a zero operator
// Construct a constant operator
many_body_operator_generic(ScalarType const& x);
// Copy-constructor
many_body_operator_generic(many_body_operator_generic const&);
// Move-constructor
many_body_operator_generic(many_body_operator_generic&&) = default;

Three factory functions can be used to construct nontrivial operators:

// Annihilation
template<typename ScalarType = real_or_complex, typename... IndexTypes> many_body_operator_generic<ScalarType> c(IndexTypes... ind);
// Creation
template<typename ScalarType = real_or_complex, typename... IndexTypes> many_body_operator_generic<ScalarType> c_dag(IndexTypes... ind);
// Number of particles
template<typename ScalarType = real_or_complex, typename... IndexTypes> many_body_operator_generic<ScalarType> n(IndexTypes... ind);

IndexTypes is an arbitrarily long sequence of index types, each being int or std::string.

Creation and annihilation operators obey the canonical anticommutation relation

\[\hat c^\dagger_{\mathrm{ind}_1} \hat c_{\mathrm{ind}_2} + \hat c_{\mathrm{ind}_2} \hat c^\dagger_{\mathrm{ind}_1} = \delta_{{\mathrm{ind}_1},{\mathrm{ind}_2}},\]

and the number of particle is defined as

\[\hat n_\mathrm{ind} = \hat c^\dagger_\mathrm{ind} \hat c_\mathrm{ind}.\]

There is no need to preregister valid values of ind before they are used to create an elementary operator. This means, that an algebra can be extended with new basis elements on-the-fly, after some expressions have been created.

Overloaded operations

many_body_operator_generic class defines a number of arithmetic operations with objects of the class and constants of type ScalarType. If A and B are objects of class many_body_operator_generic (instantiated with the same scalar and index types) and x is an instance of ScalarType, then the following expressions are valid:

// Addition
A + B; A + x; x + A;
A += B; A += x;

// Subtraction
A - B; A - x; x - A;
A -= B; A -= x; -A;

// Multiplication
A*B; x*A; A*x;
A *= B; A *= x;

// Division by scalar
A / x;
A /= x;

The result of any of the defined operations is guaranteed to preserve its normally ordered form.

many_body_operator_generic can be copy-constructed and assigned from another many_body_operator_generic instantiation with a compatible scalar type. For example, it is possible to copy-construct many_body_operator_complex from many_body_operator_real, but not vice versa.

An instance of many_body_operator_generic can be inserted into an output stream, provided ScalarType supports insertion into the stream.

many_body_operator_generic<double> x = c(0);
many_body_operator_generic<double> y = c_dag(1);

std::cout << (x + y)*(x - y) << std::endl; // prints "2*C^+(1)C(0)"

Member types

using scalar_t = ScalarType;

Accessor to the ScalarType.


bool is_zero() const;

Returns true if this operator is a precise zero.

triqs::hilbert_space::fundamental_operator_set make_fundamental_operator_set() const;

Returns fundamental_operator_set containing all indices met within this operator.

static many_body_operator_generic make_canonical(bool is_dag, indices_t indices);

Returns a canonical operator (creation, if is_dag = true, annihilation otherwise) with given indices.

Free functions

many_body_operator_generic<ScalarType> dagger(many_body_operator_generic<ScalarType> const& op);

Returns the Hermitian conjugate of op.

many_body_operator_generic<ScalarType> real(many_body_operator_generic<ScalarType> const& op);

Returns a copy of op with the imaginary parts of all monomial coefficients set to zero.

many_body_operator_generic<ScalarType> imag(many_body_operator_generic<ScalarType> const& op);

Returns a copy of op with the real parts of all monomial coefficients set to zero.

template<typename L>
many_body_operator_generic<ScalarType> transform(many_body_operator_generic<ScalarType> const& op, Lambda&& L);

Transforms op by applying a given functor L to each monomial. The functor must take two arguments convertible from monomial_t (see next paragraph) and ScalarType respectively, and return a new coefficient of the monomial.

Iteration over monomials

The aim of many_body_operator_generic is to have a class allowing to encode different operator expressions in C++ in the form closest to the mathematical notation. At the same time, one would like to explicitly extract the structure of a given operator (to calculate its matrix elements, for instance). For this purpose many_body_operator_generic exposes the following part of its interface:

  • using indices_t = std::vector<triqs::utility::variant_int_string>;

    A vector of indices. Each index is a variant type with two options: int or std::string.

  • struct canonical_ops_t

    This structure represents an elementary operator (basis element of the algebra).

    struct canonical_ops_t {
        bool dagger;       // true = creation, false = annihilation
        indices_t indices; // values of indices
  • using monomial_t = std::vector<canonical_ops_t>;

    A normally ordered sequence of elementary operators (monomial).

  • using const_iterator = ...;

    A bidirectional constant iterator to the list of monomials. It can be dereferenced into a special proxy object, which contains two data members: scalar_t coef and monomial_t const& monomial.

  • begin()/cbegin()

    Returns const_iterator pointing at the first monomial.

  • end()/cend()

    Returns const_iterator pointing past the end.

Here is an example of use:

using Op = many_body_operator;
Op H = -0.5*(n(0) + n(1)) + n(0)*n(1);

for(Op::const_iterator it = H.begin(); it != H.end(); ++it){
    double coef = it->coef;
    auto monomial = it->monomial;

    std::cout << "Coefficient: " << coef << std::endl;
    std::cout << "Monomial: " << std::endl;
    for(auto const& o : monomial){
        std::cout << "dagger: " << o.dagger << " index: " << o.indices[0] << " "; // only 1 index per elementary operator
    std::cout << std::endl;

The output should be

Coefficient: -0.5
dagger: 1 index: 0 dagger: 0 index: 0
Coefficient: -0.5
dagger: 1 index: 1 dagger: 0 index: 1
Coefficient: 1
dagger: 1 index: 0 dagger: 1 index: 1 dagger: 0 index: 1 dagger: 0 index: 0

Serialization & HDF5

Objects of many_body_operator_generic are ready to be serialized/deserialized with Boost.Serialization. This also allows to transparently send/receive them through Boost.MPI calls.

Writing to/reading from HDF5 is supported for the polymorphic version of the operators (many_body_operator) provided they contain no terms beyond quartic.


Python wrapper for many_body_operator class is called Operator. It is found in module pytriqs.operators.operators :

from pytriqs.operators.operators import Operator, c, c_dag, n

It corresponds to a specialized version of many_body_operator_generic: real_or_complex as the scalar type and two indices. All arithmetic operations implemented in C++ are also available in Python as well as special methods __repr__() and __str__().

from pytriqs.operators.operators import *
from itertools import product

C_list = [c(1,0),c(2,0)]
Cd_list = [c_dag(1,0), c_dag(2,0)]

print "Anticommutators:"
for Cd,C in product(Cd_list,C_list):
    print "{", Cd, ",", C, "} =", Cd*C + C*Cd

print "Commutators:"
for Cd,C in product(Cd_list,C_list):
    print "[", Cd, ",", C, "] =", Cd*C - C*Cd

x = c('A',0)
y = c_dag('B',0)
print "x =", x
print "y =", y

print "Algebra:"

print "-x =", -x
print "x + 2.0 =", x + 2.0
print "2.0 + x =", 2.0 + x
print "x - 2.0 =", x - 2.0
print "2.0 - x =", 2.0 - x
print "3.0*y =", 3.0*y
print "y*3.0 =", y*3.0
print "x + y =", x + y
print "x - y =", x - y
print "(x + y)*(x - y) =", (x + y)*(x - y)

print "x*x is zero:", (x*x).is_zero()
print "dagger(x) = ", dagger(x)

{ 1*c_dag(1,0) , 1*c(1,0) } = 1
{ 1*c_dag(1,0) , 1*c(2,0) } = 0
{ 1*c_dag(2,0) , 1*c(1,0) } = 0
{ 1*c_dag(2,0) , 1*c(2,0) } = 1
[ 1*c_dag(1,0) , 1*c(1,0) ] = -1 + 2*c_dag(1,0)*c(1,0)
[ 1*c_dag(1,0) , 1*c(2,0) ] = 2*c_dag(1,0)*c(2,0)
[ 1*c_dag(2,0) , 1*c(1,0) ] = 2*c_dag(2,0)*c(1,0)
[ 1*c_dag(2,0) , 1*c(2,0) ] = -1 + 2*c_dag(2,0)*c(2,0)
x = 1*c('A',0)
y = 1*c_dag('B',0)
-x = -1*c('A',0)
x + 2.0 = 2 + 1*c('A',0)
2.0 + x = 2 + 1*c('A',0)
x - 2.0 = -2 + 1*c('A',0)
2.0 - x = 2 + -1*c('A',0)
3.0*y = 3*c_dag('B',0)
y*3.0 = 3*c_dag('B',0)
x + y = 1*c_dag('B',0) + 1*c('A',0)
x - y = -1*c_dag('B',0) + 1*c('A',0)
(x + y)*(x - y) = 2*c_dag('B',0)*c('A',0)
x*x is zero: True
dagger(x) =  1*c_dag('A',0)