Skip to content
Tutorials
Dynamical Mean Field Theory (DMFT) Solvers
DMFT-02 Hybridization

DMFT-02 Hybridization

Hybridization Expansion CT-HYB

We start by running a continuous-time quantum Monte Carlo code: the hybridization-expansion algorithm CT-HYB. As an example we reproduce Fig. 11 in the DMFT review by Georges et al., Rev. Mod. Phys. 68, 13 (1996). The series of six curves shows how the system, a Hubbard model on the Bethe lattice with interaction U=3D/2U=3D/\sqrt{2} at half filling, enters an antiferromagnetic phase upon cooling. In Tutorials 03 and 07 we will reproduce the same results with the interaction-expansion continuous-time solver and with the discrete-time Hirsch-Fye Quantum Monte Carlo code, respectively. The input parameters are the same, apart from a few solver-related parameters.

Running the simulation

The CT-HYB simulation will run for roughly 1 hour in total if you want to reproduce all 6 curves in Fig. 11 mentioned above. The files for this tutorial can be found in the directory tutorials/dmft-02-hybridization.

All DMFT tutorials can be started using a python script. The python script generates parameter files, runs them, and plots the results. You can run the short script tutorial2.py, reproducing only 2 out of the 6 curves (runtime: roughly 20 minutes), or the long version tutorial2_long.py, reproducing all 6 curves in the figure (runtime: roughly 1 hour).

The python script tutorial2.py automatically prepares the input files for the 2 simulations, parm_beta_6.0 and parm_beta_12.0, and runs them (/path-to-alps-installation/bin/dmft parm_beta_x).

import pyalps
import numpy as np
import matplotlib.pyplot as plt
import pyalps.plot

#prepare the input parameters
parms=[]
for b in [6., 12.]:
    parms.append(
        {
             'ANTIFERROMAGNET'     : 1,
             'CONVERGED'           : 0.003,
             'FLAVORS'             : 2,
             'H'                   : 0,
             'H_INIT'              : 0.03*b/8.,
             'MAX_IT'              : 6,
             'MAX_TIME'            : 300,
             'MU'                  : 0,
             'N'                   : 250,
             'NMATSUBARA'          : 250,
             'N_MEAS'              : 10000,
             'OMEGA_LOOP'          : 1,
             'SEED'                : 0,
             'SITES'               : 1,
             'SOLVER'              : 'hybridization',
             'SC_WRITE_DELTA'      : 1,
             'SYMMETRIZATION'      : 0,
             'U'                   : 3,
             't'                   : 0.707106781186547,
             'SWEEPS'              : int(10000*b/16.),
             'THERMALIZATION'      : 1000,
             'BETA'                : b
        }
    )


#write the input file and run the simulation
for p in parms:
    input_file = pyalps.writeParameterFile('parm_beta_'+str(p['BETA']),p)
    res = pyalps.runDMFT(input_file)

Understanding the input parameters

The input file parm_beta_6.0 produced by the script above, with added comments on the parameters:

H_INIT = 0.0225   //  initial magnetic field in the direction of quantization axis, to produce the initial Weiss field
ANTIFERROMAGNET = 1   // allow antiferromagnetic ordering; in single site DMFT it is meaningfull only for bipartite lattices
SEED = 0   // Monte Carlo Random Number Seed 
CONVERGED = 0.003   // criterion for the convergency of the iterations
MAX_IT = 6   // upper limit on the number of iterations (the selfconsistency may be stopped before if criterion based on CONVERGED is reached)
SWEEPS = 3750   // Total number of sweeps to be computed (the solver may be stopped before reaching this limit on run-time limit set by MAX_TIME)
FLAVORS = 2   // flavors 0 and 1 correspond to spin up and down
SYMMETRIZATION = 0   // We are not enforcing a paramagnetic self consistency condition (symmetry in flavor 0 and 1)
NMATSUBARA = 250   // The cut-off for Matsubara frequencies 
H = 0   // Magnetic field in the direction of quantization axis
OMEGA_LOOP = 1   // the selfconsistency runs in Matsubara frequencies
SITES = 1   // number of sites of the impurity: for single site DMFT simulation it is 1
N = 250   // auxiliary discretization of the imaginary-time Green's function
BETA = 6.0   // Inverse temperature
U = 3   // Interaction strength
MAX_TIME = 300   // Upper time limit in seconds to run the impurity solver (per iteration)
SC_WRITE_DELTA = 1   // option for selfconsistency to write the hybridization function for the impurity solver
N_MEAS = 10000   // number of updates in between measurements
SOLVER = "hybridization"   // The Hybridization solver
THERMALIZATION = 1000   // Thermalization Sweeps 
MU = 0   // Chemical potential; for particle-hole symmetric models corresponds MU=0 to half-filled case
t = 0.707106781187   // hopping parameter; for the Bethe lattice considered here $W=2D=4t$

Note that there is no parameter specifying the band structure or lattice type. By default a Bethe lattice is assumed, but this can be changed (see DMFT-08 Lattices).

A specification of the initial Weiss field (set by the variables G0OMEGA_INPUT or G0TAU_INPUT) is missing as well - the program will thus at initialization compute the non-interacting Green’s function. It will use the initial magnetic field H_INIT, which produces in this case a small difference between flavors (0 and 1, representing ↑\uparrow and ↓\downarrow) to start away from the paramagnetic solution. The reason for this is that in very short simulations (like this tutorial), starting from a paramagnetic Weiss field could mean that the random noise would not produce enough difference in the first few iterations to move the system away from the paramagnetic regime; a badly converged paramagnet would then appear as a solution. The dependence of H_INIT on BETA serves to optimize the run, lowering the number of iterations needed.

Iterating to self-consistency

The code will run for up to 6 self-consistency iterations. For a more precise simulation you can raise this number; the simulation will stop early once the convergence criterion specified by the parameter CONVERGED is reached. In the directory in which you run the program you will find Green’s function files G_tau_i, the self-energies (selfenergy_i), and Green’s functions in Matsubara representation (frequency space) G_omega_i. G_tau in these examples has two entries: a spin-up and a spin-down column. The entry at τ=β−\tau=\beta^- is the negative occupation (density); from this we can obtain the magnetization of the system.

Error bars may be estimated from successive iterations on a converged system.

To rerun a simulation, you can specify a starting solution by defining the input parameter G0OMEGA_INPUT: copy the desired G0omega_output to filename_X and specify 'G0OMEGA_INPUT':'filename_X' in the python script (or G0OMEGA_INPUT=filename_X in the input file directly), then rerun the code.

Plotting the antiferromagnetic transition

As in Fig. 11 of the DMFT review by Georges et al., you can observe the transition to the antiferromagnetic phase by plotting the Green’s functions in the imaginary-time representation (part of tutorial2.py and tutorial2_long.py):

listobs=['0', '1']   # we will plot both flavors 0 and 1

# load the imaginary-time Green's function in final iteration from all result files
data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='parm_beta_*h5'), respath='/simulation/results/G_tau', what=listobs)
for d in pyalps.flatten(data):
    d.x = d.x*d.props["BETA"]/float(d.props["N"])   # rescale horizontal axis
    d.props['label'] = r'$\beta=$'+str(d.props['BETA'])+'; flavor='+str(d.props['observable'][len(d.props['observable'])-1])

plt.figure()
plt.xlabel(r'$\tau$')
plt.ylabel(r'$G_{flavor}(\tau)$')
plt.title('DMFT-02: Neel transition for the Hubbard model on the Bethe lattice\n(using the Hybridization expansion impurity solver)')
pyalps.plot.plot(data)
plt.legend()
plt.show()

You will notice that the results are relatively noisy. This is because the expansion order at such high temperatures is very small, which makes the measurement procedure inefficient. You can improve statistics by increasing the total runtime (MAX_TIME) and/or the number of SWEEPS. Note that running the solver directly under MPI via the SOLVER parameter (e.g. SOLVER = "mpirun -np procs /path-to-ALPS-installation/bin/hybridization") currently does not work, due to an issue with the path prefix; consult your MPI installation’s documentation (e.g. the mpirun man page) for the correct way to launch parallel jobs on your system.

Checking convergence

If you want to check the convergence of your DMFT self-consistency, you can plot the Green’s functions of different iterations using tutorial2eval.py, whose code is shown here:

listobs=['0']   # we look at a single flavor (=0) 
res_files = pyalps.getResultFiles(pattern='parm_*.h5')  # we look for result files

## load all iterations of G_{flavor=0}(tau)
data = pyalps.loadDMFTIterations(res_files, observable="G_tau", measurements=listobs, verbose=False)

## create a figure for each BETA
grouped = pyalps.groupSets(pyalps.flatten(data), ['BETA'])
for sim in grouped:
    common_props = pyalps.dict_intersect([ d.props for d in sim ])
    
    ## rescale x-axis and set label
    for d in sim:
        d.x = d.x * d.props['BETA']/float(d.props['N'])
        d.props['label'] = 'it'+d.props['iteration']
    
    ## plot all iterations for this BETA
    plt.figure()
    plt.xlabel(r'$\tau$')
    plt.ylabel(r'$G_{flavor=0}(\tau)$')
    plt.title('Simulation at ' + r'$\beta = {beta}$'.format(beta= common_props['BETA']))
    pyalps.plot.plot(sim)
    plt.legend()

plt.show()

Note that the iteration-resolved results are loaded by a different function (pyalps.loadDMFTIterations) than the final results (pyalps.loadMeasurements), because the iteration-resolved data is stored using a different folder structure (/simulation/iteration/number/results/) than the ALPS default (/simulation/results/) used for storing the final results.

As mentioned above, the occupation nfn_f is equal to Gf(τ=β−)G_f(\tau=\beta^-), the last entry of the imaginary-time Green’s function for flavor ff. Code for printing the final occupancies and plotting them versus β\beta is part of tutorial2eval.py:

## load the final iteration of G_{flavor=0}(tau)
data_G_tau = pyalps.loadMeasurements(res_files, respath='/simulation/results/G_tau', what=listobs, verbose=False)  

print("Occupation in the last iteration at flavor=0")
for d in pyalps.flatten(data_G_tau):
    # obtain occupation using relation: <n_{flavor=0}> = -<G_{flavor=0}(tau=beta)>
    d.y = np.array([-d.y[-1]])
    print("n_0(beta =",d.props['BETA'],") =",d.y[0])
    d.x = np.array([0])
    d.props['observable'] = 'occupation'

occupation = pyalps.collectXY(data_G_tau, 'BETA', 'occupation')
for d in occupation:
    d.props['line']="scatter"

plt.figure()
pyalps.plot.plot(occupation)
plt.xlabel(r'$\beta$')
plt.ylabel(r'$n_{flavor=0}$')
plt.title('Occupation versus BETA')

plt.show()

Matsubara-frequency Green’s function and self-energy

Since our self-consistency runs in Matsubara frequencies (recall parameter OMEGA_LOOP=1), the convergence criterion is max∣Gfit(iωn)−Gfit+1(iωn)∣<\mathrm{max}|G_{f}^{it}(i\omega_n)-G_{f}^{it+1}(i\omega_n)|\lt CONVERGED. The imaginary part (the real part analogously) of the Matsubara-frequency Green’s function is plotted by

from math import pi

## load all iterations of G_{flavor=0}(i omega_n)
data = pyalps.loadDMFTIterations(pyalps.getResultFiles(pattern='parm_*.h5'), observable="G_omega", measurements=listobs, verbose=False)

## create a figure for each BETA
grouped = pyalps.groupSets(pyalps.flatten(data), ['BETA'])
for sim in grouped:
    common_props = pyalps.dict_intersect([ d.props for d in sim ])
    
    ## rescale x-axis and set label
    for d in sim:
        d.x = np.array([(2.*n+1)*pi/common_props['BETA'] for n in d.x])
        d.y = np.array(d.y.imag)
        d.props['label'] = "it"+d.props['iteration']
        d.props['line']="scatter"
        d.props['fillmarkers'] = False
    
    ## plot all iterations for this BETA
    plt.figure()
    plt.xlabel(r'$i\omega_n$')
    plt.ylabel(r'$Im\ G_{flavor=0}(i\omega_n)$')
    plt.title('Simulation at ' + r'$\beta = {beta}$'.format(beta=common_props['BETA']))
    pyalps.plot.plot(sim)
    plt.legend()

plt.show()

It is usually best to observe convergence in the self-energy, which is much more sensitive. Note that longer simulations are required to obtain smoother Green’s functions and self-energies; in this simulation the noise in the intermediate range of Matsubara frequencies is very strong. The self-energy is obtained via the Dyson equation as Σfit(iωn)=G0fit(iωn)−1−Gfit(iωn)−1\Sigma_f^{it}(i\omega_n)=G0_f^{it}(i\omega_n)^{-1}-G_f^{it}(i\omega_n)^{-1}, and its imaginary part is plotted by this fragment of tutorial2eval.py (real part analogously):

## load all iterations of G_{flavor=0}(i omega_n) and G0_{flavor=0}(i omega_n)
data_G = pyalps.loadDMFTIterations(pyalps.getResultFiles(pattern='parm_*.h5'), observable="G_omega", measurements=listobs, verbose=False)
data_G0 = pyalps.loadDMFTIterations(pyalps.getResultFiles(pattern='parm_*.h5'), observable="G0_omega", measurements=listobs, verbose=False)

## create a figure for each BETA
grouped_G = pyalps.groupSets(pyalps.flatten(data_G), ['BETA','observable'])
for sim in grouped_G:
    common_props = pyalps.dict_intersect([ d.props for d in sim ])
    
    ## compute selfenergy using the Dyson equation, rescale x-axis and set label
    for d_G in sim:
        # find corresponding dataset from data_G0
        d_G0 = [s for s in pyalps.flatten(data_G0) if s.props['iteration']==d_G.props['iteration'] and s.props['BETA']==common_props['BETA']][0]
        d_G.x = np.array([(2.*n+1)*pi/common_props['BETA'] for n in d_G.x])
        # Dyson equation
        Sigma = np.array([1./d_G0.y[w] - 1./d_G.y[w] for w in range(len(d_G.y))])
        d_G.y = np.array(Sigma.imag)
        d_G.props['label'] = "it"+d_G.props['iteration']
        d_G.props['line']="scatter"
        d_G.props['fillmarkers'] = False
    
    ## plot all iterations for this BETA
    plt.figure()
    plt.xlabel(r'$i\omega_n$')
    plt.ylabel(r'$Im\ \Sigma_{flavor=0}(i\omega_n)$')
    plt.title('Simulation at ' + r'$\beta = {beta}$'.format(beta=common_props['BETA']))
    pyalps.plot.plot(sim)
    plt.legend()

plt.show()
DMFT-01 IntroDMFT-03 Interaction