# [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;
}