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

From ALPS
Jump to: navigation, search
(Tutorial 02: Hybridization Expansion CT-HYB)
 
(4 intermediate revisions by 2 users not shown)
Line 3: Line 3:
  
 
=Tutorial 02: Hybridization Expansion CT-HYB =
 
=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 ''it 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.
+
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.
  
The CT-HYB simulation will run in total roughly 2 hours 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>.  
+
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>.  
  
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: 2 hours) or its short version reproducing only 2 out of the 6 curves (running roughly 20 minutes)
+
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)
  
 
  alpspython tutorial2.py
 
  alpspython tutorial2.py
  
(vispython for Mac instead of alpspython) or the longer version to reproduce all curves in the Fig.11 of [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''it et al.'']
+
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.'']
  
 
  alpspython tutorial2_long.py
 
  alpspython tutorial2_long.py
Line 29: Line 29:
 
             {
 
             {
 
               'ANTIFERROMAGNET'    : 1,
 
               'ANTIFERROMAGNET'    : 1,
               'CONVERGED'          : 0.005,
+
               'CONVERGED'          : 0.003,
 
               'FLAVORS'            : 2,
 
               'FLAVORS'            : 2,
 
               'H'                  : 0,
 
               'H'                  : 0,
               'H_INIT'              : 0.05,
+
               'H_INIT'              : 0.03*b/8.,
               'MAX_IT'              : 10,
+
               'MAX_IT'              : 6,
 
               'MAX_TIME'            : 300,
 
               'MAX_TIME'            : 300,
 
               'MU'                  : 0,
 
               'MU'                  : 0,
               'N'                  : 500,
+
               'N'                  : 250,
               'NMATSUBARA'          : 500,
+
               'NMATSUBARA'          : 250,
 
               'N_MEAS'              : 10000,
 
               'N_MEAS'              : 10000,
 
               'OMEGA_LOOP'          : 1,
 
               'OMEGA_LOOP'          : 1,
Line 47: Line 47:
 
               'U'                  : 3,
 
               'U'                  : 3,
 
               't'                  : 0.707106781186547,
 
               't'                  : 0.707106781186547,
               'SWEEPS'              : 2500,
+
               'SWEEPS'              : int(10000*b/16.),
 
               'THERMALIZATION'      : 1000,
 
               'THERMALIZATION'      : 1000,
 
               'BETA'                : b
 
               'BETA'                : b
Line 61: Line 61:
 
The input file <tt>parm_beta_6.0</tt> produced by the script above with added comments on the parameters:
 
The input file <tt>parm_beta_6.0</tt> produced by the script above with added comments on the parameters:
  
  H_INIT = 0.05   //  initial magnetic field in the direction of quantization axis, to produce the initial Weiss field
+
  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
 
  ANTIFERROMAGNET = 1  // allow antiferromagnetic ordering; in single site DMFT it is meaningfull only for bipartite lattices
 
  SEED = 0  // Monte Carlo Random Number Seed  
 
  SEED = 0  // Monte Carlo Random Number Seed  
  CONVERGED = 0.005   // criterion for the convergency of the iterations
+
  CONVERGED = 0.003   // criterion for the convergency of the iterations
  MAX_IT = 10   // upper limit on the number of iterations (the selfconsistency may be stopped before if criterion based on CONVERGED is reached)
+
  MAX_IT = 6   // upper limit on the number of iterations (the selfconsistency may be stopped before if criterion based on CONVERGED is reached)
  SWEEPS = 2500   // Total number of sweeps to be computed (the solver may be stopped before reaching this limit on run-time limit set by MAX_TIME)
+
  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
 
  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)
 
  SYMMETRIZATION = 0  // We are not enforcing a paramagnetic self consistency condition (symmetry in flavor 0 and 1)
  NMATSUBARA = 500   // The cut-off for Matsubara frequencies  
+
  NMATSUBARA = 250   // The cut-off for Matsubara frequencies  
 
  H = 0  // Magnetic field in the direction of quantization axis
 
  H = 0  // Magnetic field in the direction of quantization axis
 
  OMEGA_LOOP = 1  // the selfconsistency runs in Matsubara frequencies
 
  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
 
  SITES = 1  // number of sites of the impurity: for single site DMFT simulation it is 1
  N = 500   // auxiliary discretization of the imaginary-time Green's function
+
  N = 250   // auxiliary discretization of the imaginary-time Green's function
 
  BETA = 6.0  // Inverse temperature
 
  BETA = 6.0  // Inverse temperature
 
  U = 3  // Interaction strength
 
  U = 3  // Interaction strength
Line 86: Line 86:
 
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 ]]).  
 
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 ]]).  
  
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.
+
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.
  
The code will run for up to 10 (15 in long version <tt>tutorial2_long.py</tt>) self-consistency iterations. 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.  
+
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.  
  
 
Error bars may be estimated via successive iterations on a converged system.
 
Error bars may be estimated via successive iterations on a converged system.
Line 112: Line 112:
 
  plt.show()
 
  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 (<tt>MAX_TIME</tt>). 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.
+
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:
 
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:
Line 144: Line 144:
 
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.
 
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 occupancies and plotting then is a part of [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2eval.py tutorial2eval.py],
+
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)
 
  ## load the final iteration of G_{flavor=0}(tau)
Line 200: Line 200:
 
  plt.show()
 
  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. 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 plotted by this fragment of [http://alps.comp-phys.org/static/tutorials2.2.0/dmft-02-hybridization/tutorial2eval.py tutorial2eval.py]:
+
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)
 
  ## load all iterations of G_{flavor=0}(i omega_n) and G0_{flavor=0}(i omega_n)

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