ALPS 2 Tutorials:DMFT-02 Hybridization

Revision as of 18:51, 29 June 2012 by Jimriska (talk | contribs) (Tutorial 02: Hybridization Expansion CT-HYB)

Jump to: navigation, search

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 5 seconds (1 minute if you have revision of ALPS older than 6238) per iteration. The parameter files for running this simulation may 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 parm_beta_x' 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 (in fact, the solver will stop on run-time limit)
MAX_TIME = 5;  Maximum time to run the impurity solver 
BETA = 14.0;  Inverse temperature 
SITES = 1;  This is a single site DMFT simulation, so Sites is 1 
N = 500;  Number of time bins for Green's function measurement
NMATSUBARA = 500;  The number of Matsubara frequencies 
U = 3;  Interaction energy 
t = 0.707106781187;  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
FLAVORS = 2;  spin up and down
H_INIT = 0.05;   initial magnetic field to produce the initial Weiss field
ANTIFERROMAGNET = 1;  allow antiferromagnetic ordering
CONVERGED = 0.005;  criterion for the convergency of the iterations
MAX_IT = 15;  maximum number of iterations
OMEGA_LOOP = 1;   the selfconsistency runs in Matsubara frequencies

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

The code will run for up to 15 self-consistency iterations. 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 frequency space G_omega_i 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 (this statement does hold for converged solution and for selfconsistency in imaginary time - with OMEGA_LOOP = 0; in other case you have to observe fluctuations in successive iterations).

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 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 the same, apart from few solver-related parameters.

You can find the parameter files (called parm_beta_x) in the directory tutorials/dmft-02-hybridization in the subdirectories. 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

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

#prepare the input parameters
for b in [6.,8.,10.,12.,14.,16.]:
             'ANTIFERROMAGNET'     : 1,
             'CONVERGED'           : 0.005,
             'FLAVORS'             : 2,
             'H'                   : 0,
             'H_INIT'              : 0.05,
             'MAX_IT'              : 15,
             'MAX_TIME'            : 5,
             'MU'                  : 0,
             'N'                   : 500,
             'NMATSUBARA'          : 500,
             'N_MEAS'              : 10000,
             'N_ORDER'             : 50,
             'OMEGA_LOOP'          : 1,
             'SEED'                : 0,
             'SITES'               : 1,
             'SOLVER'              : 'Hybridization',
             'SYMMETRIZATION'      : 0,
             'U'                   : 3,
             't'                   : 0.707106781186547,
             'SWEEPS'              : 100000000,
             'THERMALIZATION'      : 1000,
             'BETA'                : b,
             'CHECKPOINT'          : 'solverdump_beta_'+str(b)

# NOTE: in revision of ALPS older than 6238, the MAX_TIME will effectively be 60 seconds.        
# For more precise calculations we propose to you to:
#   enhance the MAX_TIME (to 60), 
#   lower the CONVERGED (to 0.003), 
#   increase MAX_IT (to 20)
#   raise N and NMATSUBARA (to 1000)
# ( the runtime of the script with changed parameters will be roughly 2 hours )

#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 the input parameter G0OMEGA_INPUT, e.g. copy the desired G0omega_output to G0_omega_input, specify G0OMEGA_INPUT = G0_omega_input in the parameter file and rerun the code. As in the DMFT-07 The Hirsch-Fye solver you can observe the transition to the antiferromagnetic phase by plotting the Green's functions using the script:

for f in range(0,flavors):

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'])+'; flavor='+str(d.props['observable'][len(d.props['observable'])-1])

plt.title('DMFT-02: Neel transition for the Hubbard model on the Bethe lattice\n(using the Hybridization expansion impurity solver)')

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

listobs=['Green_0']   # we look at convergence of a single flavor (=0)
for p in parms:
    data = ll.ReadDMFTIterations(pyalps.getResultFiles(pattern='parm_beta_'+str(p['BETA'])+'.h5'), measurements=listobs, verbose=True)
    grouped = pyalps.groupSets(pyalps.flatten(data), ['iteration'])
    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'] = 'it'+r.props['iteration']
        nd.append( r )
    plt.title('DMFT-02: ' + r'$\beta = %.4s$' %nd[0].props['BETA'])

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.