# Tutorial 07: Hirsch Fye Code

We start by running a discrete time Monte Carlo code: the Hirsch Fye code. 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 Hirsch Fye algorithm is described in here, and this review also provides an open source implementation for the codes. More information can also be found in Blümer's PhD. While many improvements have been developed (see e.g. Alvarez08 or Nukala09), the algorithm has been replaced by continuous-time algorithms.

The Hirsch Fye simulation will run for about a minute per iteration. The parameter files for running this simulation can be found in the directory tutorials/dmft-07-hirschfye.

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-07-hirschfye, adjusting the 'solver' path 'SOLVER = /opt/alps/bin/hirschfye' in hirschfye.param to point to the hirsch fye executable and running the dmft code '/opt/alps/bin/dmft hirschfye.param' leads to the same results.

The main parameters are:

SEED = 0; //Monte Carlo Random Number Seed
THERMALIZATION = 10000;  Thermalization Sweeps
SWEEPS = 1000000; Total Sweeps to be computed
MAX_TIME = 60;  Maximum time to run the simulation
BETA = 12.;  Inverse temperature
SITES = 1;  This is a single site DMFT simulation, so Sites is 1
N = 16;  Number of time slices (you will see that this parameter is rather small)
NMATSUBARA = 500;  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 = /opt/alps/bin/hirschfye;  The path to the external Hirsch Fye solver


To start a simulation type:

dmft hirschfye.param


or run the python script tutorial1a.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,
'FLAVORS'             : 2,
'H'                   : 0,
'MAX_IT'              : 10,
'MAX_TIME'            : 60,
'MU'                  : 0,
'N'                   : 16,
'NMATSUBARA'          : 500,
'OMEGA_LOOP'          : 1,
'SEED'                : 0,
'SITES'               : 1,
'SOLVER'              : '/opt/alps/bin/hirschfye',
'SYMMETRIZATION'      : 0,
'TOLERANCE'           : 0.01,
'U'                   : 3,
't'                   : 0.707106781186547,
'SWEEPS'              : 1000000,
'THERMALIZATION'      : 10000,
'BETA'                : b,
'G0OMEGA_INPUT'       : 'G0_omega_input_beta_'+str(b),
'BASENAME'            : 'hirschfye.param'
}
)

#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)


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

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()


As a discrete time method, HF suffers from $\Delta\tau$ - errors. Pick a set of parameters and run it for sucessively larger N! Also: you're running the DMFT simulation using an (almost) converged input bath function. By deleting the file G0_omega_input you can restart the calculation from the free solution and observe convergence.

If you want to use Vistrails to run your DMFT-simulations, you can use dmft-07-hirschfye.vt.