[triqs/hilbert_space] Many-body states and Hilbert spaces

hilbert_space toolbox provides the user with a small set of abstractions, useful for doing calculations in the Hamiltonian picture for finite fermionic systems. Among these abstarctions are Hilbert spaces, their subspaces, states and operators acting on them. There is also an efficient implementation of the automatic partitioning algorithm that allows to split a Hilbert space into a number of subspaces, each invariant under action of a given operator.

Example of use

#include <bitset>

#include <triqs/operators/many_body_operator.hpp>

#include <triqs/hilbert_space/fundamental_operator_set.hpp>
#include <triqs/hilbert_space/hilbert_space.hpp>
#include <triqs/hilbert_space/state.hpp>
#include <triqs/hilbert_space/imperative_operator.hpp>
#include <triqs/hilbert_space/space_partition.hpp>

using namespace triqs::operators;
using namespace triqs::hilbert_space;
using namespace std;

int main() {

  // Model system: half-filled Hubbard dimer

  // Parameters
  double U  = 2.0;   // Coulomb repulsion
  double mu = U / 2; // Chemical potential
  double t  = 0.2;   // Hopping amplitude

  // Set of fundamental operators
  // Every fundamental operator carries two indices, one atom index and one spin projection
  fundamental_operator_set fops;
  fops.insert(1, "up");
  fops.insert(1, "down");
  fops.insert(2, "up");
  fops.insert(2, "down");

  // Full Hilbert space of the dimer
  hilbert_space hs(fops);

  // Hamiltonian
  auto H = U * n(1, "up") * n(1, "down") + U * n(2, "up") * n(2, "down");
  H += -mu * (n(1, "up") + n(1, "down") + n(2, "up") + n(2, "down"));

  for (auto s : {"up", "down"}) H += -t * (c_dag(1, s) * c(2, s) + c_dag(2, s) * c(1, s));

  // Imperative form of the Hamiltonian acting on the full Hilbert space
  using full_hs_operator_t = imperative_operator<hilbert_space, double, false>;
  full_hs_operator_t hamiltonian(H, fops);

  // Many-body state in the full Hilbert space
  using state_t = state<hilbert_space, double, true /* using sparse storage */>;
  state_t st(hs);

  // Partition the Hilbert space
  space_partition<state_t, full_hs_operator_t> partition(st, hamiltonian);

  // Preallocate a container for the invariant subspaces of H
  vector<sub_hilbert_space> inv_subspaces;
  for (int n = 0; n < partition.n_subspaces(); ++n) inv_subspaces.emplace_back(n);

  // Fill the subspaces with Fock spaces
  foreach (partition, [&](int s, int spn) { inv_subspaces[spn].add_fock_state(hs.get_fock_state(s)); })
    ;

  // Observables: N_up and N_dn
  auto N_up = n(1, "up") + n(2, "up"), N_down = n(1, "down") + n(2, "down");

  // Subspace connections generated by the observables
  auto N_up_conn   = partition.find_mappings(full_hs_operator_t(N_up, fops));
  auto N_down_conn = partition.find_mappings(full_hs_operator_t(N_down, fops));

  // Imperative versions of N_up and N_down acting between subspaces
  using sub_hs_operator_t = imperative_operator<sub_hilbert_space, double, true>;

  vector<int> N_up_hilbert_map(inv_subspaces.size(), -1);
  for (auto const &conn : N_up_conn) N_up_hilbert_map[conn.first] = conn.second;
  sub_hs_operator_t N_up_op(N_up, fops, N_up_hilbert_map, &inv_subspaces);

  vector<int> N_down_hilbert_map(inv_subspaces.size(), -1);
  for (auto const &conn : N_down_conn) N_down_hilbert_map[conn.first] = conn.second;
  sub_hs_operator_t N_down_op(N_down, fops, N_down_hilbert_map, &inv_subspaces);

  // Calculate and print expectation values of N_up and N_down
  for (int spn = 0; spn < inv_subspaces.size(); ++spn) {
    // Iterate over invariant subspaces
    auto const &subspace = inv_subspaces[spn];
    cout << "Invariant subspace " << spn << endl;
    cout << "====================" << endl;
    cout << "State    N_up   N_dn" << endl;
    cout << "====================" << endl;

    // Iterate over basis Fock states in the selected subspace
    for (int i = 0; i < subspace.size(); ++i) {

      // Many-body state defined on the subspace
      state<sub_hilbert_space, double, false> ket(subspace);
      ket(i) = 1;

      cout << "|" << bitset<4>(subspace.get_fock_state(i)) << ">   ";
      cout << ((N_up_hilbert_map[spn] != -1) ? dot_product(ket, N_up_op(ket)) : 0) << "      ";
      cout << ((N_down_hilbert_map[spn] != -1) ? dot_product(ket, N_down_op(ket)) : 0) << endl;
    }
  }

  return 0;
}