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

From ALPS
Jump to: navigation, search
(Tutorial 02: Hybridization Expansion CT-HYB)
Line 5: Line 5:
 
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.
 
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.
  
The CT-HYB simulation will run for about a minute per iteration. The parameter files for running this simulation can be found in the directory <tt>tutorials/dmft-02-hybridization</tt>.  
+
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 <tt>tutorials/dmft-02-hybridization</tt>.  
  
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 hybrid.parm' leads to the same results.
+
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:
 
The main parameters are:
Line 13: Line 13:
 
  SEED = 0; Monte Carlo Random Number Seed  
 
  SEED = 0; Monte Carlo Random Number Seed  
 
  THERMALIZATION = 1000;  Thermalization Sweeps  
 
  THERMALIZATION = 1000;  Thermalization Sweeps  
  SWEEPS = 100000000; Total Sweeps to be computed  
+
  SWEEPS = 100000000; Total Sweeps to be computed (in fact, the solver will stop on run-time limit)
  MAX_TIME = 60;  Maximum time to run the impurity solver  
+
  MAX_TIME = 5;  Maximum time to run the impurity solver  
  BETA = 12.;  Inverse temperature  
+
  BETA = 14.0;  Inverse temperature  
 
  SITES = 1;  This is a single site DMFT simulation, so Sites is 1  
 
  SITES = 1;  This is a single site DMFT simulation, so Sites is 1  
  N = 1000;  Number of time bins for Green's function measurement
+
  N = 500;  Number of time bins for Green's function measurement
  NMATSUBARA = 1000;  The number of Matsubara frequencies  
+
  NMATSUBARA = 500;  The number of Matsubara frequencies  
 
  U = 3;  Interaction energy  
 
  U = 3;  Interaction energy  
  t = 1;  hopping parameter. For the Bethe lattice considered here $W=2D=4t$
+
  t = 0.707106781187;  hopping parameter. For the Bethe lattice considered here $W=2D=4t$
 
  MU = 0;  Chemical potential  
 
  MU = 0;  Chemical potential  
 
  H = 0;  Magnetic field  
 
  H = 0;  Magnetic field  
 
  SYMMETRIZATION = 0;  We are not enforcing a paramagnetic self consistency condition  
 
  SYMMETRIZATION = 0;  We are not enforcing a paramagnetic self consistency condition  
 
  SOLVER = Hybridization;  The Hybridization solver
 
  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).  
 
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:
 
To start a simulation type:
  
  dmft hybrid.parm
+
  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 <tt>G_tau_i</tt> as well the self energies (<tt>selfenergy_i</tt>) and Green's functions in frequency space <tt>G_omega_i</tt> in your output directory. <tt>G_tau</tt> in these examples has two entries: a spin-up and a  spin-down column. The entry at <tt>\beta</tt> 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).
  
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 <tt>G_tau_?</tt> as well the self energies (<tt>selfenergy_?</tt>) and Green's functions in frequency space <tt>G_omega_?</tt> in your output directory. <tt>G_tau</tt> in these examples has two entries: a spin-up and a  spin-down column. The entry at <tt>\beta</tt> 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 <math>\beta</math> and compare your result to Fig. 11 of [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''it et al.''].
 
You can run the following lines in the python shell in order to plot the Green's functions for different <math>\beta</math> and compare your result to Fig. 11 of [http://dx.doi.org/10.1103/RevModPhys.68.13 Georges ''it et al.''].
In tutorial 07 we will reproduce the same results with a discrete-time Quantum Monte Carlo code: the Hirsch Fye code. The parameters are the same, apart from the command for the solver.
+
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 <tt>*tsqrt2</tt>) in the directory <tt>tutorials/dmft-02-hybridization</tt> in the examples. If you want to use Vistrails to run your DMFT-simulations, you can use [http://alps.comp-phys.org/static/tutorials2.1.0/dmft-02-hybridization/dmft-02-hybridization.vt dmft-02-hybridization.vt].
+
You can find the parameter files (called <tt>parm_beta_x</tt>) in the directory <tt>tutorials/dmft-02-hybridization</tt> in the subdirectories. If you want to use Vistrails to run your DMFT-simulations, you can use [http://alps.comp-phys.org/static/tutorials2.1.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.1.0/dmft-02-hybridization/tutorial2a.py tutorial2a.py]:
 
Alternatively you can run the python script [http://alps.comp-phys.org/static/tutorials2.1.0/dmft-02-hybridization/tutorial2a.py tutorial2a.py]:
  
Line 48: Line 55:
 
     parms.append(
 
     parms.append(
 
             {
 
             {
              'ANTIFERROMAGNET'    : 1,
+
              'ANTIFERROMAGNET'    : 1,
              'CONVERGED'          : 0.03,
+
              'CONVERGED'          : 0.005,
              'F'                  : 10,
+
              'FLAVORS'            : 2,
              'FLAVORS'            : 2,
+
              'H'                  : 0,
              'H'                  : 0,
+
              'H_INIT'              : 0.05,
              'H_INIT'              : 0.2,
+
              'MAX_IT'              : 15,
              'MAX_IT'              : 10,
+
              'MAX_TIME'            : 5,
              'MAX_TIME'            : 60,
+
              'MU'                  : 0,
              'MU'                  : 0,
+
              'N'                  : 500,
              'N'                  : 1000,
+
              'NMATSUBARA'          : 500,
              'NMATSUBARA'          : 1000,
+
              'N_MEAS'              : 10000,
              'N_FLIP'              : 0,
+
              'N_ORDER'            : 50,
              'N_MEAS'              : 10000,
+
              'OMEGA_LOOP'          : 1,
              'N_MOVE'              : 0,
+
              'SEED'                : 0,
              'N_ORDER'            : 50,
+
              'SITES'              : 1,
              'N_SHIFT'            : 0,
+
              'SOLVER'              : 'Hybridization',
              'OMEGA_LOOP'          : 1,
+
              'SYMMETRIZATION'      : 0,
              'OVERLAP'            : 0,
+
              'U'                  : 3,
              'SEED'                : 0,
+
              't'                  : 0.707106781186547,
              'SITES'              : 1,
+
              'SWEEPS'              : 100000000,
              'SOLVER'              : 'Hybridization',
+
              'THERMALIZATION'      : 1000,
              'SYMMETRIZATION'      : 0,
+
              'BETA'                : b,
              'TOLERANCE'          : 0.01,
+
              'CHECKPOINT'          : 'solverdump_beta_'+str(b)
              'U'                  : 3,
 
              't'                  : 0.707106781186547,
 
              'SWEEPS'              : 100000000,
 
              'THERMALIZATION'      : 1000,
 
              'BETA'                : b,
 
              'CHECKPOINT'          : 'solverdump_beta_'+str(b)+'.task1.out.h5',
 
              'G0TAU_INPUT'        : 'G0_tau_input_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
 
  #write the input file and run the simulation
 
  for p in parms:
 
  for p in parms:
Line 86: Line 94:
 
     res = pyalps.runDMFT(input_file)
 
     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 <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.  
+
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 <tt>G0OMEGA_INPUT</tt>, e.g. copy the desired <tt>G0omega_output</tt> to <tt>G0_omega_input</tt>, specify <tt>G0OMEGA_INPUT = G0_omega_input</tt> in the parameter file and rerun the code.  
As in the [http://alps.comp-phys.org/mediawiki/index.php/ALPS_2_Tutorials:DMFT-01_Hirsch-Fye Tutorial DMFT-07 Hirsch-Fye] you can observe the transition to the antiferromagnetic phase by plotting the Green's functions using the script:
+
As in the [[ALPS_2_Tutorials:DMFT-07_Hirsch-Fye DMFT-07 The Hirsch-Fye solver  ]] you can observe the transition to the antiferromagnetic phase by plotting the Green's functions using the script:
  
 
  flavors=parms[0]['FLAVORS']
 
  flavors=parms[0]['FLAVORS']
Line 98: Line 106:
 
  for d in pyalps.flatten(data):
 
  for d in pyalps.flatten(data):
 
     d.x = d.x*d.props["BETA"]/float(d.props["N"])
 
     d.x = d.x*d.props["BETA"]/float(d.props["N"])
     d.props['label'] = r'$\beta=$'+str(d.props['BETA'])
+
     d.props['label'] = r'$\beta=$'+str(d.props['BETA'])+'; flavor='+str(d.props['observable'][len(d.props['observable'])-1])
 +
 
  plt.figure()
 
  plt.figure()
 
  plt.xlabel(r'$\tau$')
 
  plt.xlabel(r'$\tau$')
  plt.ylabel(r'$G(\tau)$')
+
  plt.ylabel(r'$G_{flavor}(\tau)$')
  plt.title('Hubbard model on the Bethe lattice')
+
  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)
 
  pyalps.plot.plot(data)
 
  plt.legend()
 
  plt.legend()
Line 110: Line 119:
  
 
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.1.0/dmft-02-hybridization/tutorial2b.py tutorial2b.py]:
 
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.1.0/dmft-02-hybridization/tutorial2b.py tutorial2b.py]:
 +
 +
listobs=['Green_0']  # we look at convergence of a single flavor (=0)
  
 
  ll=pyalps.load.Hdf5Loader()
 
  ll=pyalps.load.Hdf5Loader()
Line 121: Line 132:
 
         r.x = np.array([e*group[0].props['BETA']/float(group[0].props['N']) for e in group[0].x])
 
         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 = group[0].props
         r.props['label'] = r.props['iteration']
+
         r.props['label'] = 'it'+r.props['iteration']
 
         nd.append( r )
 
         nd.append( r )
 
     plt.figure()
 
     plt.figure()
 
     plt.xlabel(r'$\tau$')
 
     plt.xlabel(r'$\tau$')
     plt.ylabel(r'$G(\tau)$')
+
     plt.ylabel(r'$G_{flavor=0}(\tau)$')
     plt.title(r'\beta = %.4s' %nd[0].props['BETA'])
+
     plt.title('DMFT-02: ' + r'$\beta = %.4s$' %nd[0].props['BETA'])
 
     pyalps.plot.plot(nd)
 
     pyalps.plot.plot(nd)
 
     plt.legend()
 
     plt.legend()
 +
 
  plt.show()
 
  plt.show()
  
 
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.
 
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.

Revision as of 18:51, 29 June 2012


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 tutorial2a.py:

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

#prepare the input parameters
parms=[]
for b in [6.,8.,10.,12.,14.,16.]:
    parms.append(
            {
             '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:

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'])+'; 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 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 tutorial2b.py:

listobs=['Green_0']   # we look at convergence of a single flavor (=0)
ll=pyalps.load.Hdf5Loader()
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'])
    nd=[]
    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.figure()
    plt.xlabel(r'$\tau$')
    plt.ylabel(r'$G_{flavor=0}(\tau)$')
    plt.title('DMFT-02: ' + r'$\beta = %.4s$' %nd[0].props['BETA'])
    pyalps.plot.plot(nd)
    plt.legend()

plt.show()

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.