# [triqs/atom_diag] Lightweight exact diagonalization solver and tools¶

atom_diag is a simple class allowing to find energy levels and eigenstates of finite systems of fermions (atomic-like problems).

Diagonalization is preceeded by reduction of the Hamiltonian to a block-diagonal form or, in other words, by finding its invariant subspaces in the full Hilbert space of the problem. This step is done either using user-supplied quantum number operators, or in an automatic fashion.

A diagonalization object can then be passed to a number of utility functions to compute equilibrium density matrix, thermal expectation values of observables and Matsubara Green’s functions at a given temperature.

atom_diag comes with C++ and Python interfaces. It may prove useful as a component to build quantum impurity solvers (hybridization expansion CT-QMC, Hubbard-I, NCA, …), which require solution of an auxiliary atomic problem. However, it is no substitute for a large scale exact diagonalization solver, since it can only treat problems of a moderate size.

## Example of use: Python¶

from triqs.gf import *
from triqs.operators import *
from triqs.operators.util.hamiltonians import h_int_kanamori
from h5 import HDFArchive
from triqs.atom_diag import *
import numpy as np
from itertools import product

# Definition of a 3-orbital atom
spin_names = ('up','dn')
orb_names = [0,1,2]
# Set of fundamental operators
fops = [(sn,on) for sn, on in product(spin_names,orb_names)]

# Numbers of particles with spin up/down
N_up = n('up',0) + n('up',1) + n('up',2)
N_dn = n('dn',0) + n('dn',1) + n('dn',2)

# Construct Hubbard-Kanamori Hamiltonian
U = 3.0 * np.ones((3,3))
Uprime = 2.0 * np.ones((3,3))
J_hund = 0.5

H = h_int_kanamori(spin_names, orb_names, U, Uprime, J_hund, True)

H += -4.0 * (N_up + N_dn)

# Add hopping terms between orbitals 0 and 1 with some complex amplitude
H += 0.1j * (c_dag('up',0) * c('up',1) - c_dag('up',1) * c('up',0))
H += 0.1j * (c_dag('dn',0) * c('dn',1) - c_dag('dn',1) * c('dn',0))

# Split H into blocks and diagonalize it using N_up and N_dn quantum numbers
print(ad.n_subspaces) # Number of invariant subspaces, 4 * 4 = 16

# Now split using the total number of particles, N = N_up + N_dn

# Split the Hilbert space automatically

# Partition function for inverse temperature \beta=3
beta = 3

# Equilibrium density matrix

# Expectation values of orbital double occupancies

# Atomic Green's functions
gf_struct = [['dn',orb_names],['up',orb_names]]
G_w = atomic_g_w(ad, beta, gf_struct, (-2, 2), 400, 0.01)
G_tau = atomic_g_tau(ad, beta, gf_struct, 400)
G_iw = atomic_g_iw(ad, beta, gf_struct, 100)
G_l = atomic_g_l(ad, beta, gf_struct, 20)

# Finally, we save our AtomDiag object for later use
with HDFArchive('atom_diag_example.h5') as ar:


## Example of use: C++¶

#include <triqs/atom_diag/atom_diag.hpp>
#include <triqs/atom_diag/functions.hpp>
#include <triqs/atom_diag/gf.hpp>

using namespace triqs::gfs;
using namespace triqs::operators;
using namespace triqs::atom_diag;

int main(int argc, const char *argv[]) {

// Hamiltonian: single band Hubbard atom + spin flips + anomalous terms
auto h = 2.0 * (n("up", 0) - 0.5) * (n("dn", 0) - 0.5);
h += 0.3 * (c_dag("up", 0) * c("dn", 0) + c_dag("dn", 0) * c("up", 0));
h += 0.1 * (c_dag("up", 0) * c_dag("dn", 0) - c("up", 0) * c("dn", 0));

// Set of fundamental operators C/C^+
triqs::hilbert_space::fundamental_operator_set fops;
fops.insert("up", 0);
fops.insert("dn", 0);

// Diagonalize the problem ('false' = no complex-valued matrix elements in h)

// Inverse temperature \beta = 10 and structure of Green's functions
double beta           = 10;
gf_struct_t gf_struct = {{"dn", {0}}, {"up", {0}}};

// Direct computation of atomic Green's functions in
// diferent representations

int n_tau        = 200;
int n_iw         = 100;
unsigned int n_l = 20;
int n_w          = 200;

// Imaginary time
auto G_tau = atomic_g_tau(ad, beta, gf_struct, n_tau);
// Matsubara frequencies
auto G_iw = atomic_g_iw(ad, beta, gf_struct, n_iw);
// Legendre polynomial basis
auto G_l = atomic_g_l(ad, beta, gf_struct, n_l);
// Real frequencies
auto G_w = atomic_g_w(ad, beta, gf_struct, {-2.0, 2.0} /* energy window */, n_w, 0.01 /* broadening */);

// Indirect computation of atomic Green's functions via
// Lehmann representation

auto lehmann = atomic_g_lehmann(ad, beta, gf_struct);

// Lehmann representation object can be reused but it requires
// extra memory to store
auto G_tau_ind = atomic_g_tau<false>(lehmann, gf_struct, {beta, Fermion, n_tau});
auto G_iw_ind  = atomic_g_iw<false>(lehmann, gf_struct, {beta, Fermion, n_iw});
auto G_l_ind   = atomic_g_l<false>(lehmann, gf_struct, {beta, Fermion, n_l});
auto G_w_ind   = atomic_g_w<false>(lehmann, gf_struct, {-2.0, 2.0, n_w}, 0.01);

return 0;
}