Difference between revisions of "ALPS 2 Tutorials:DMFT-02 Hybridization"

From ALPS
Jump to: navigation, search
(New page: =Tutorial 02: Hybridization Expansion CT-HYB / Tutorial 03: Interaction Expansion CT-INT= We will now reproduce the same result with a continuous-time Quantum Monte Carlo code: the Hybrid...)
 
(Tutorial 02: Hybridization Expansion CT-HYB)
 
(41 intermediate revisions by 7 users not shown)
Line 1: Line 1:
 +
{{Languages|ALPS_2_Tutorials:DMFT-02_Hybridization}}
  
=Tutorial 02: Hybridization Expansion CT-HYB / Tutorial 03: Interaction Expansion CT-INT=
 
We will now reproduce the same result with a continuous-time Quantum Monte Carlo code: the Hybridization or CT-HYB code  by [http://dx.doi.org/10.1103/PhysRevB.74.155107 Werner ''et al.'']. The parameters are the same, apart from the command for the solver:
 
  
SOLVER=Hybridization
+
=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 [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''et al.'']. The series of six curves shows how the system, a Hubbard model on the Bethe lattice with interaction <math>U=3D/\sqrt{2}</math> 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 Quantum Monte Carlo Hirsch-Fye code, respectively. The input parameters are there same, apart from few solver-related parameters.
  
You can find the parameter files (called <tt>*tsqrt2</tt>) in the directory <tt>tutorials/dmft-02-hybridization</tt> in the examples.
+
The CT-HYB simulation will run in total roughly 1 hour if you want to reproduce all 6 curves in the Fig. 11 mentioned above. The files for this tutorial may be found in the directory <tt>tutorials/dmft-02-hybridization</tt>.  
  
After running these simulations compare the output to the Hirsch Fye results. To rerun a simulation, you can specify a starting solution by defining <tt>G0OMEGA_INPUT</tt>, e.g. copy <tt>G0omga_output</tt> to <tt>G0_omega_input</tt>, specify <tt>G0OMEGA_INPUT = G0_omega_input</tt> in the parameter file and rerun the code.
+
All DMFT tutorials can either be started using vistrails or a python script. Both the vistrails scripts and the python scripts generate parameter files, run them, and plot the results. If you want to use Vistrails to run your DMFT-simulations, you can use [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/dmft-02-hybridization.vt dmft-02-hybridization.vt]. Alternatively you can run the python script [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2_long.py tutorial2_long.py] (runtime: roughly 1 hour) or its short version reproducing only 2 out of the 6 curves (running roughly 20 minutes)
  
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 (<tt>MAX_TIME</tt>) or by running it on more than one CPU. For running it with MPI, try <tt>mpirun -np procs dmft parameter_file</tt> or consult the man page of your mpi installation.
+
alpspython tutorial2.py
  
It is also instructive to run these calculations with a CT-INT code. This code performs an [http://link.aps.org/doi/10.1103/PhysRevB.72.035122 expansion in the interaction] (instead of the hybridization).
+
or the longer version to reproduce all curves in the Fig.11 of [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''et al.'']
The corresponding parameter files are very similar, you can find them in the directory <tt>tutorials/dmft-03-interaction</tt>.
 
  
=Tutorial 04: Mott Transition=
+
alpspython tutorial2_long.py
Mott transitions are metal insulator transitions that occur in many materials, e.g. transition metal compounds, as a function of pressure or doping. The [http://dx.doi.org/10.1103/RevModPhys.70.1039 review by Imada ''et al.''] gives an excellent introduction to the subject and mentions V2O3 and the organics as typical examples.
 
  
MIT are easily investigated by DMFT as the relevant physics is essentially local (or k-independent): At half filling the MIT can by modeled by a self energy with a pole at $\omega=0$ which splits the noninteracting band into an upper and a lower Hubbard band. In this context it is instructive to suppress antiferromagnetic long range order and enforce a paramagnetic solution in the DMFT simulation, to mimic the paramagnetic insulating phase. For this the up and down spin of the Green's functions are symmetrized (parameter <tt>SYMMETRIZATION = 1;</tt>).  
+
The python script <tt>tutorial2.py</tt> automatically prepares the input files for the 2 simulations, <tt>parm_beta_6.0</tt> and <tt>parm_beta_12.0</tt>, and runs them (/path-to-alps-installation/bin/dmft parm_beta_x).
  
We investigate the ''Mott'' transition in single-site DMFT, as a function of interaction at fixed temperature <tt>\beta t=20</tt> (see e.g. Fig. 2 [http://dx.doi.org/10.1103/PhysRevB.76.235123 in this paper]).
+
import pyalps
Starting from a non-interacting solution we see in the imaginary time Green's function that the solution is metallic for <tt>U/t \leq 4.5</tt>, and insulating for <tt>U/t\geq 5</tt>. A coexistence region could be found by starting from an insulating (or atomic) solution and trying to convert it for smaller $U$.
+
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)
  
Imaginary time Green's functions are not easy to interpret, and therefore many authors [http://dx.doi.org/10.1016/0370-1573(95)00074-7 employ analytic continuation methods]. There are however two clear features: the value at <tt>\beta</tt> corresponds to <tt>-n</tt>, the negative value of the density (per spin). And the value at <tt>\beta/2</tt> goes for decreasing temperature to <tt>A(\omega=0)</tt>, the spectral function at the Fermi energy: <tt>\beta G(\beta/2 \rightarrow A(0)</tt>. From a temperature dependence of the imaginary time Green's function we can therefore immediately see if the system is metallic or insulating.
+
The input file <tt>parm_beta_6.0</tt> produced by the script above with added comments on the parameters:
  
The <tt>beta20_U*</tt> parameter files in the directory <tt>tutorials/dmft-04-mott</tt> should show you how to go from a metallic (at small $U$) to an insulating (at large $U$) solution, at fixed $\beta$. The largest value of $U$ is deep within the insulating phase.
+
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$
  
=Tutorial 05: Orbitally Selective Mott Transition=
+
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 [[ALPS_2_Tutorials:DMFT-08_Lattices | DMFT-08 Setting a particular lattice ]]).  
An interesting phenomenon in multi-orbital models is the orbitally selective Mott transition, first examined [http://www.springerlink.com/content/l9fktwe5l2hdtq3w by Anisimov ''et al'']. A variant of this, a ''momentum-selective'' Mott transition, has recently been discussed in [http://link.aps.org/abstract/PRB/v80/e045120 cluster calculations] as a cluster representation of pseudogap physics.
 
  
In an orbitally selective Mott transition some of the orbitals involved become Mott insulating as a function of doping or interactions, while others stay metallic.
+
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 <math>\uparrow,\ \downarrow</math>) to start away from the paramagnetic solution - the reason for that is that in very short simulations (like this tutorial) starting from paramagnetic Weiss field it could happen that the random noise would not produce enough difference in first few iterations to get the system away from paramagnet. A badly converged paramagnet would then appear as a solution. The dependence of H_INIT on BETA serves for optimization of the run, lowering the number of needed iterations.
  
As a minimal model we consider two bands: a wide band  and a narrow band. In addition to the intra-orbital Coulomb repulsion $U$ we consider interactions $U'$, and $J$, with $U' = U-2J$. We limit ourselves to Ising-like interactions - a simplification that is often problematic for real compounds.
+
The code will run for up to 6 self-consistency iterations. For a precise simulation one shall raise this number and the simulation shall stop on the convergency criterion specified by the parameter CONVERGED. In the directory in which you run the program you will find Green's functions files <tt>G_tau_i</tt> as well the self energies (<tt>selfenergy_i</tt>) and Green's functions in Matsubara representation (frequency space) <tt>G_omega_i</tt>. <tt>G_tau</tt> in these examples has two entries: a spin-up and a  spin-down column. The entry at <tt>\tau=\beta^-</tt> is the negative occupation (density); by that we may get magnetization of the system.  
  
We choose here a case with two bandwidth $t1=0.5$ and $t2=1$ and density-density like interactions of $U'=U/2$, $J=U/4$, and $U$ between $1.6$ and $2.8$, where the first case shows a Fermi liquid-like behavior in both orbitals, the $U=2.2$ is orbitally selective, and $U=2.8$ is insulating in both orbitals.
+
Error bars may be estimated via successive iterations on a converged system.
  
The parameter files can be found in <tt>tutorials/dmft-05-osmt</tt> in the directories <tt>beta30_U1.8_2orbital</tt>, <tt>beta30_U2.2_2orbital</tt> and <tt>beta30_U2.8_2orbital</tt>. A paper using the same sample parameters can be found [http://dx.doi.org/10.1103/PhysRevB.72.081103 here].
+
To rerun a simulation, you can specify a starting solution by defining the input parameter G0OMEGA_INPUT, e.g. copy the desired <tt>G0omega_output</tt> to <tt>filename_X</tt> and specify input parameter 'G0OMEGA_INPUT':'filename_X' in the python script (or G0OMEGA_INPUT=filename_X in the input file directly) and rerun the code.  
  
Tutorial by [[User:Gullc|Emanuel]] - Please don't hesitate to ask!
+
As in the Fig. 11 in the DMFT review [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''it et al.''] you can observe the transition to the antiferromagnetic phase by plotting the Green's functions in imaginary-time represention (part of <tt>tutorial2.py</tt> and <tt>tutorial2_long.py</tt>):
 +
 
 +
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. 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 runtime (<tt>MAX_TIME</tt>) and/or the number of SWEEPS. The solver may run on more than one CPU using MPI, try <tt>SOLVER = "mpirun -np procs /path-to-ALPS-installation/bin/hybridization"</tt> (currently not working due to issue with path prefix) 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 [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2eval.py tutorial2eval.py], the corresponding 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 = %.4s$' % common_props['BETA'])
 +
    pyalps.plot.plot(sim)
 +
    plt.legend()
 +
 +
plt.show()
 +
 
 +
It has to be noted, that the iteration-resolved results are loaded by a different function (<tt>pyalps.loadDMFTIterations</tt>) as the final results (<tt>pyalps.loadMeasurements</tt>), because the iteration-resolved data is stored using a different folder structure (<tt>/simulation/iteration/number/results/</tt>) than the ALPS default (<tt>/simulation/results/</tt>) used for storage of the final results.
 +
 
 +
As already mentioned above, the occupation <math>n_f</math> equals to <math>G_f(\tau=\beta^-)</math>, which is the last entry of the imaginary time Green's function for flavor <math>f</math>. Code for printing the final occupancies and plotting them vs <math>\beta</math> is a part of [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2eval.py 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) 
 +
 +
beta, occupation = [], []
 +
print "Occupation in the last iteration at flavor=0"
 +
for i in range(len(data_G_tau)):
 +
    # obtain occupation using relation: <n_{flavor=0}> = -<G_{flavor=0}(tau=beta)>
 +
    beta_,n_ = float(data_G_tau[i][0].props['BETA']),-data_G_tau[i][0].y[-1]
 +
    print "n_0(beta =",beta_,") =",n_
 +
    beta.append(beta_)
 +
    occupation.append(n_)
 +
 +
d = pyalps.DataSet()
 +
d.y = occupation
 +
d.x = beta
 +
d.props['line']="scatter"
 +
plt.figure()
 +
plt.xlabel(r'$\beta$')
 +
plt.ylabel(r'$n_{flavor=0}$')
 +
plt.title('Occupation versus BETA')
 +
pyalps.plot.plot(d)
 +
 +
plt.show()
 +
 
 +
As our selfconsistency is in Matsubara frequencies (recall parameter OMEGA_LOOP=1), the criterion for convergency is <math>\mathrm{max}|G_{f}^{it}(i\omega_n)-G_{f}^{it+1}(i\omega_n)|<CONVERGED</math>. The imaginary part (real part analogously) 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 = %.4s$' % common_props['BETA'])
 +
    pyalps.plot.plot(sim)
 +
    plt.legend()
 +
 +
plt.show()
 +
 
 +
It is usually best to observe convergence in the selfenergy, which is much more sensitive. Note that longer simulations are required to obtain smoother Green's functions and selfenergies; in this simulation the noise in the intermediate range of Matsubara frequencies is very strong. The selfenergy is obtained via Dyson's equation as <math>\Sigma_f^{it}(i\omega_n)=G0_f^{it}(i\omega_n)^{-1}-G_f^{it}(i\omega_n)^{-1}</math> and  its imaginary part is plotted by this fragment of [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2eval.py 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 = %.4s$' % common_props['BETA'])
 +
    pyalps.plot.plot(sim)
 +
    plt.legend()
 +
 +
plt.show()

Latest revision as of 23:11, 30 May 2014


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 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. In tutorials 03 and 07 we will reproduce the same results with the interaction expansion continuous-time solver and with the discrete-time Quantum Monte Carlo Hirsch-Fye code, respectively. The input parameters are there same, apart from few solver-related parameters.

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

All DMFT tutorials can either be started using vistrails or a python script. Both the vistrails scripts and the python scripts generate parameter files, run them, and plot the results. 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 tutorial2_long.py (runtime: roughly 1 hour) or its short version reproducing only 2 out of the 6 curves (running roughly 20 minutes)

alpspython tutorial2.py

or the longer version to reproduce all curves in the Fig.11 of Georges et al.

alpspython tutorial2_long.py

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)

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 Setting a particular lattice ).

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,\ \downarrow) to start away from the paramagnetic solution - the reason for that is that in very short simulations (like this tutorial) starting from paramagnetic Weiss field it could happen that the random noise would not produce enough difference in first few iterations to get the system away from paramagnet. A badly converged paramagnet would then appear as a solution. The dependence of H_INIT on BETA serves for optimization of the run, lowering the number of needed iterations.

The code will run for up to 6 self-consistency iterations. For a precise simulation one shall raise this number and the simulation shall stop on the convergency criterion specified by the parameter CONVERGED. In the directory in which you run the program you will find Green's functions files G_tau_i as well 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); by that we may get magnetization of the system.

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

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

As in the Fig. 11 in the DMFT review Georges it et al. you can observe the transition to the antiferromagnetic phase by plotting the Green's functions in imaginary-time represention (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. 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 runtime (MAX_TIME) and/or the number of SWEEPS. The solver may run on more than one CPU using MPI, try SOLVER = "mpirun -np procs /path-to-ALPS-installation/bin/hybridization" (currently not working due to issue with path prefix) 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 tutorial2eval.py, the corresponding 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 = %.4s$' % common_props['BETA'])
   pyalps.plot.plot(sim)
   plt.legend()

plt.show()

It has to be noted, that the iteration-resolved results are loaded by a different function (pyalps.loadDMFTIterations) as 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 storage of the final results.

As already mentioned above, the occupation n_f equals to G_f(\tau=\beta^-), which is the last entry of the imaginary time Green's function for flavor f. Code for printing the final occupancies and plotting them vs \beta is a 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)  

beta, occupation = [], []
print "Occupation in the last iteration at flavor=0"
for i in range(len(data_G_tau)):
   # obtain occupation using relation: <n_{flavor=0}> = -<G_{flavor=0}(tau=beta)>
   beta_,n_ = float(data_G_tau[i][0].props['BETA']),-data_G_tau[i][0].y[-1]
   print "n_0(beta =",beta_,") =",n_
   beta.append(beta_)
   occupation.append(n_)

d = pyalps.DataSet()
d.y = occupation
d.x = beta
d.props['line']="scatter"
plt.figure()
plt.xlabel(r'$\beta$')
plt.ylabel(r'$n_{flavor=0}$')
plt.title('Occupation versus BETA')
pyalps.plot.plot(d)

plt.show()

As our selfconsistency is in Matsubara frequencies (recall parameter OMEGA_LOOP=1), the criterion for convergency is \mathrm{max}|G_{f}^{it}(i\omega_n)-G_{f}^{it+1}(i\omega_n)|<CONVERGED. The imaginary part (real part analogously) 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 = %.4s$' % common_props['BETA'])
   pyalps.plot.plot(sim)
   plt.legend()

plt.show()

It is usually best to observe convergence in the selfenergy, which is much more sensitive. Note that longer simulations are required to obtain smoother Green's functions and selfenergies; in this simulation the noise in the intermediate range of Matsubara frequencies is very strong. The selfenergy is obtained via Dyson's equation as \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 = %.4s$' % common_props['BETA'])
   pyalps.plot.plot(sim)
   plt.legend()

plt.show()