# Tutorial 02: 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 it et al.. The series of six curves shows how the system, a Hubbard model on the Bethe lattice with interaction $U=3D/\sqrt{2}$ at half filling, enters an antiferromagnetic phase upon cooling.

The CT-HYB simulation will run for about a minute per iteration. The parameter files for running this simulation can be found in the directory tutorials/dmft-02-hybridization.

All DMFT tutorials can either be started using vistrails, or directly by invoking a command on the command line. The vistrails scripts described in the following automatically generate parameter files, run them, and plot the results. Running one of the dmft parameter sets manually, e.g. by entering the directory 'beta_14_U3_tsqrt2' in tutorials/dmft-02-hybridization, and running the dmft code '/opt/alps/bin/dmft hybrid.param' leads to the same results.

The main parameters are:

SEED = 0; Monte Carlo Random Number Seed
THERMALIZATION = 1000;  Thermalization Sweeps
SWEEPS = 100000000; Total Sweeps to be computed
MAX_TIME = 60;  Maximum time to run the impurity solver
BETA = 12.;  Inverse temperature
SITES = 1;  This is a single site DMFT simulation, so Sites is 1
N = 1000;  Number of time bins for Green's function measurement
NMATSUBARA = 1000;  The number of Matsubara frequencies
U = 3;  Interaction energy
t = 1;  hopping parameter. For the Bethe lattice considered here $W=2D=4t$
MU = 0;  Chemical potential
H = 0;  Magnetic field
SYMMETRIZATION = 0;  We are not enforcing a paramagnetic self consistency condition
SOLVER = Hybridization;  The Hybridization solver


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 later). To start a simulation type:

dmft hybrid.param


The code will run for up to 10 self-consistency iterations. In the directory in which you run the program you will find Green's functions files G_tau_? as well the self energies (selfenergy_?) and Green's functions in frequency space G_omega_? in your output directory. G_tau in these examples has two entries: a spin-up and a spin-down column. The entry at \beta is the negative density; where it is different outside of error bars the system is in an antiferromagnetic phase. You can run the following lines in the python shell in order to plot the Green's functions for different $\beta$ and compare your result to Fig. 11 of Georges it et al.. In tutorial 07 we will reproduce the same results with a discrete-time Quantum Monte Carlo code: the Hirsch Fye code. The parameters are the same, apart from the command for the solver.

You can find the parameter files (called *tsqrt2) in the directory tutorials/dmft-02-hybridization in the examples. If you want to use Vistrails to run your DMFT-simulations, you can use dmft-02-hybridization.vt. Alternatively you can run the python script tutorial2a.py:

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

#prepare the input parameters
parms=[]
for b in [6.,8.,10.,12.,14.,16.]:
parms.append(
{
'ANTIFERROMAGNET'     : 1,
'CONVERGED'           : 0.03,
'F'                   : 10,
'FLAVORS'             : 2,
'H'                   : 0,
'H_INIT'              : 0.2,
'MAX_IT'              : 10,
'MAX_TIME'            : 60,
'MU'                  : 0,
'N'                   : 1000,
'NMATSUBARA'          : 1000,
'N_FLIP'              : 0,
'N_MEAS'              : 10000,
'N_MOVE'              : 0,
'N_ORDER'             : 50,
'N_SHIFT'             : 0,
'OMEGA_LOOP'          : 1,
'OVERLAP'             : 0,
'SEED'                : 0,
'SITES'               : 1,
'SOLVER'              : 'Hybridization',
'SYMMETRIZATION'      : 0,
'TOLERANCE'           : 0.01,
'U'                   : 3,
't'                   : 0.707106781186547,
'SWEEPS'              : 100000000,
'THERMALIZATION'      : 1000,
'BETA'                : b,
'G0TAU_INPUT'         : 'G0_tau_input_beta_'+str(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)



After running these simulations compare the output to the Hirsch Fye results of tutorials 07 or the DMFT review, or to the interaction expansion results of the next tutorial. To rerun a simulation, you can specify a starting solution by defining G0OMEGA_INPUT, e.g. copy G0omga_output to G0_omega_input, specify G0OMEGA_INPUT = G0_omega_input in the parameter file and rerun the code. As in the Tutorial DMFT-07 Hirsch-Fye you can observe the transition to the antiferromagnetic phase by plotting the Green's functions using the script:

flavors=parms[0]['FLAVORS']
listobs=[]
for f in range(0,flavors):
listobs.append('Green_'+str(f))

data = ll.ReadMeasurementFromFile(pyalps.getResultFiles(pattern='parm_beta_*h5'), respath='/simulation/results/G_tau', measurements=listobs, verbose=True)
for d in pyalps.flatten(data):
d.x = d.x*d.props["BETA"]/float(d.props["N"])
d.props['label'] = r'$\beta=$'+str(d.props['BETA'])
plt.figure()
plt.xlabel(r'$\tau$')
plt.ylabel(r'$G(\tau)$')
plt.title('Hubbard model on the Bethe lattice')
pyalps.pyplot.plot(data)
plt.legend()
plt.show()


You will notice that the results are relatively noisy. The reason for that is that the expansion order at such high temperatures is very small, which renders the measurement procedure inefficient. You can improve statistics by increasing the total run time (MAX_TIME) or by running it on more than one CPU. For running it with MPI, try mpirun -np procs dmft parameter_file or consult the man page of your mpi installation.

If you want to check the convergence of your DMFT self-consistency, you can plot the Green's functions of different iterations using tutorial2b.py:

ll=pyalps.load.Hdf5Loader()
for p in parms:
grouped = pyalps.groupSets(pyalps.flatten(data), ['iteration'])
nd=[]
for group in grouped:
r = pyalps.DataSet()
r.y = np.array(group[0].y)
r.x = np.array([e*group[0].props['BETA']/float(group[0].props['N']) for e in group[0].x])
r.props = group[0].props
r.props['label'] = r.props['iteration']
nd.append( r )
plt.figure()
plt.xlabel(r'$\tau$')
plt.ylabel(r'$G(\tau)$')
plt.title(r'\beta = %.4s' %nd[0].props['BETA'])
pyalps.pyplot.plot(nd)
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.