ALPS 2.1 Tutorials:DMFT-07 Hirsch-Fye
Languages: |
English • 日本語 (ja) • 繁體中文 (zh-tw) • 简体中文 (zh) |
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 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-01-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-01-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 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)) ll=pyalps.load.Hdf5Loader() 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 - 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-01-hirschfye.vt.
Tutorial by Emanuel - Please don't hesitate to ask!