Measuring static observables

In addition to the interacting Green’s functions one is often interested in such observables, as orbital occupation numbers and double occupancy on the impurity. Expectation values of these time-independent operators are readily expressed in terms of the reduced density matrix of the system,

\[\hat\rho_\mathrm{imp} = \mathrm{Tr}_\mathrm{bath}[e^{-\beta\hat H}/Z].\]

Here, \(e^{-\beta\hat H}/Z\) is the density matrix of the full system (impurity + bath) described by the Hamiltonian \(\hat H\).

Let us consider a single-orbital Anderson impurity in a weak external magnetic field [script].

At first, we import all necessary modules, define some input parameters and construct a Solver object in the usual way.

from pytriqs.gf import *
from pytriqs.operators import *
from pytriqs.applications.impurity_solvers.cthyb import Solver

# Parameters
D = 1.0         # Half-bandwidth of the bath
V = 0.2         # Hybridisation amplitude
U = 4.0         # Coulomb interaction
e_f = -U/2      # Local energy level
h = 0.01        # External field
beta = 50       # Inverse temperature

# Construct the impurity solver with the inverse temperature
# and the structure of the Green's functions
S = Solver(beta = beta, gf_struct = [ ('up',[0]), ('down',[0]) ])

# Initialize the non-interacting Green's function S.G0_iw
# External magnetic field introduces Zeeman energy splitting between
# different spin components
S.G0_iw['up']   << inverse(iOmega_n - e_f + h/2 - V**2 * Wilson(D))
S.G0_iw['down'] << inverse(iOmega_n - e_f - h/2 - V**2 * Wilson(D))

We instruct the solve function to accumulate the density matrix by passing measure_density_matrix = True and use_norm_as_weight = True. The latter parameter tells the solver to employ a reweighting scheme – use a spherical norm instead of the trace to calculate the atomic weight. The reweighting allows to take into account important contributions to \(\hat\rho_\mathrm{imp}\), which otherwise would be missed.

# Run the solver
S.solve(h_int = U * n('up',0) * n('down',0),     # Local Hamiltonian
        n_cycles  = 500000,                      # Number of QMC cycles
        length_cycle = 200,                      # Length of one cycle
        n_warmup_cycles = 10000,                 # Warmup cycles
        measure_density_matrix = True,           # Measure the reduced density matrix
        use_norm_as_weight = True)               # Required to measure the density matrix

Results of density matrix accumulation are accessible via density_matrix attribute. We also need information about the structure of the local Hilbert space to compute expectation values of static observables. This information is stored as a special object in h_loc_diagonalization.

# Extract accumulated density matrix
rho = S.density_matrix

# Object containing eigensystem of the local Hamiltonian
h_loc_diag = S.h_loc_diagonalization

According to a well known formula, an expectation value of an operator \(\hat O\) acting on the impurity degrees of freedom is given by

\[\langle\hat O\rangle = \mathrm{Tr}_\mathrm{at}[\hat\rho_\mathrm{imp} \hat O].\]

There is a function called trace_rho_op(), which does this trace and returns the expectation value of \(\hat O\).

from pytriqs.atom_diag import trace_rho_op

# Evaluate occupations
print "<N_up> =", trace_rho_op(rho, n('up',0), h_loc_diag)
print "<N_down> = ", trace_rho_op(rho, n('down',0), h_loc_diag)

# Evaluate double occupancy
print "<N_up*N_down> =", trace_rho_op(rho, n('up',0)*n('down',0), h_loc_diag)

Typical output of the script may look like:

<N_up> = 0.620679103675
<N_down> =  0.379442861526
<N_up*N_down> = 0.00369097179335