Measuring the two-particle Green’s function for our AIM example

Let us revisit the example of a one-orbital Anderson impurity embedded in a flat (Wilson) conduction bath and look at all the changes to the Python script needed to calculate the two-particle Green’s function in particle-hole frequency convention using worm sampling:

from triqs.gf import *
from triqs.operators import *
from h5 import HDFArchive
import triqs.utility.mpi as mpi

from w2dyn_cthyb import Solver

# Parameters
D, V, U = 1.0, 0.2, 4.0
e_f, beta = -U/2.0, 50

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

# Initialize the non-interacting Green's function S.G0_iw
for name, g0 in S.G0_iw: g0 << inverse(iOmega_n - e_f - V**2 * Wilson(D))

# Run the solver. The results will be in S.G2_iw_ph, errors (for parallel runs)
# in S.G2_iw_ph_error
S.solve(h_int = U * n('up',0) * n('down',0),     # Local Hamiltonian
        n_cycles  = 1000000,                     # Number of QMC cycles
        length_cycle = 20,                       # Length of one cycle
        n_warmup_cycles = 100000,                # Warmup cycles
        measure_G2_iw_ph = True,                 # Measure G2_iw_ph
        measure_G2_n_bosonic = 30,               # Number of non-negative bosonic frequencies
        measure_G2_n_fermionic = 30,             # Number of non-negative fermionic frequencies
        # Do not measure one-particle quantities to avoid Z-space sampling
        measure_G_l = False,                     # Do not measure G_l
        measure_G_tau = False)                   # Do not measure G_tau

# Save the results in an HDF5 file (only on the master node)
if mpi.is_master_node():
    with HDFArchive("aim_2p_solution.h5",'w') as Results:
        Results["G2_iw_ph"] = S.G2_iw_ph

Running this script takes about 5 minutes and generates an HDF5 archive file called aim_2p_solution.h5. This file contains the two-particle Green’s function on two fermionic and one bosonic frequency.

In detail, necessary or useful changes to calculate the two-particle Green’s function in addition to or instead of the one-particle Green’s function are:

        length_cycle = 20,                       # Length of one cycle

The length of one cycle can be reduced considerably as we use worm sampling, so even in cases where the hybridization operators of the configuration have not changed much since the last measurement, the change of the position of the worm operators can lead to a significantly different contribution to the measurement.

        measure_G2_iw_ph = True,                 # Measure G2_iw_ph
        measure_G2_n_bosonic = 30,               # Number of non-negative bosonic frequencies
        measure_G2_n_fermionic = 30,             # Number of non-negative fermionic frequencies

Here, we set the parameter that enables measurement of the two-particle Green’s function in particle-hole convention, which we have to do explicitly, and the number of bosonic and fermionic Matsubara frequencies. The numbers of frequencies here are the same as the default values, which will lead to a result with one bosonic Matsubara axis with a total (including negative) of 59 frequencies and two fermionic Matsubara axes with a total of 60 frequencies each.

        measure_G_l = False,                     # Do not measure G_l
        measure_G_tau = False)                   # Do not measure G_tau

We additionally turn off the measurement of the two one-particle quantities that would usually be measured by default, so the part of the calculation where partition function configurations are sampled can be skipped entirely, as this interface can currently only be used to measure the two-particle Green’s function using worm sampling. If we wanted to measure the one-particle Green’s function in the same solver call (and with the same parameters, including length_cycle, which may not be ideal), we could also leave these parameters at their default values.

        Results["G2_iw_ph"] = S.G2_iw_ph

Finally, in order to write it to an HDF5 file, we access the two-particle Green’s function, which is stored in the solver object’s attribute S.G2_iw_ph as a Block2Gf after the call to solve(). Additionally, if the calculation was parallelized using MPI, we could also read estimates for the standard errors from S.G2_iw_ph_error.

It should be noted that the measurement of the two-particle Green’s function is currently only available for models with diagonal hybridization and a diagonal quadratic part of the local Hamiltonian. As the interface does not try to identify this case automatically, all blocks in the block structure must have size 1 or the solver will terminate when the calculation is attempted.