[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.

• [atom_diag] Exact diagonalization of an atomic problem
• [atom_diag_functions] Utility functions for atom_diag

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)

# Add chemical potential
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
ad = AtomDiag(H, fops, [N_up,N_dn])
print(ad.n_subspaces) # Number of invariant subspaces, 4 * 4 = 16

# Now split using the total number of particles, N = N_up + N_dn
ad = AtomDiag(H, fops, [N_up+N_dn])
print(ad.n_subspaces) # 7

# Split the Hilbert space automatically
ad = AtomDiag(H, fops)
print(ad.n_subspaces) # 28

# Partition function for inverse temperature \beta=3
beta = 3
print(partition_function(ad, beta))

# Equilibrium density matrix
dm = atomic_density_matrix(ad, beta)

# Expectation values of orbital double occupancies
print(trace_rho_op(dm, n('up',0) * n('dn',0), ad))
print(trace_rho_op(dm, n('up',1) * n('dn',1), ad))
print(trace_rho_op(dm, n('up',2) * n('dn',2), ad))

# 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:
    ar['ad'] = ad

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)
  auto ad = triqs::atom_diag::atom_diag<false>(h, fops);

  // 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;
}