Single-shot DFT+DMFT

After having set up the hdf5 archive, we can now proceed to our first DFT+DMFT calculation. It consists of initialization steps, and the actual DMFT self-consistency loop, With the code snippets below you can build your own script and target it to your needs. Little examples on mixing and on restarting from a previous calculation at the end of this page should also demonstrate how simple you can modify your own DMFT script. A full working calculation for SrVO3 is discussed in the next section.

Initialization of the calculation

Before doing the actual calculation, we have to initialize all needed objects. The first thing is the SumkDFT class. It contains all basic routines that are necessary to perform a summation in k-space to get the local quantities used in DMFT. It is initialized by:

from triqs_dft_tools.sumk_dft import *
SK = SumkDFT(hdf_file = filename + '.h5')

Setting up the impurity solver

The next step is to setup an impurity solver. There are different solvers available within the TRIQS framework. E.g. for SrVO3, we will use the hybridization expansion CTHYB solver. Later on, we will see also the example of the Hubbard-I solver. They all have in common, that they are called by an uniform command:

S.solve(params)

where params are the solver parameters and depend on the actual solver. Setting up the CTHYB solver for SrVO3 is discussed on the next page. Here, let us now perform the DMFT loop using the methods of DFTTools, assuming that we have already set up a working solver instance.

Doing the DMFT loop

Having initialized the Sumk class and the solver, we can proceed with the actual DMFT part of the calculation. We set up the loop over DMFT iterations and the self-consistency condition:

n_loops = 15
for iteration_number in range(n_loops) :        # start the DMFT loop
    SK.set_Sigma([ S.Sigma ])                   # Put self energy to the SumK class
    chemical_potential = SK.calc_mu()           # calculate the chemical potential for the given density
    S.G_iw << SK.extract_G_loc()[0]                     # extract the local Green function
    S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))    # finally get G0, the input for the solver

    S.solve(h_int=h_int, **p)                   # now solve the impurity problem

    dm = S.G_iw.density()                                           # Density matrix of the impurity problem
    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=1) # Set the double counting term
    SK.save(['chemical_potential','dc_imp','dc_energ'])             # Save data in the hdf5 archive

These steps are enough for a basic DMFT Loop. After the self-consistency steps, which lead to a new \(G^0(i\omega)\), the impurity solver is called. Different to model calculations, we have to do a few more steps after this, because of the double-counting correction. We first calculate the density of the impurity problem. Then, the routine calc_dc takes as parameters this density matrix, the Coulomb interaction, Hund’s rule coupling, and the type of double-counting that should be used. Possible values for use_dc_formula are:

• 0: Full-localised limit (FLL)

• 1: DC formula as given in K. Held, Adv. Phys. 56, 829 (2007).

• 2: Around-mean-field (AMF)

At the end of the calculation, we can save the Green function and self energy into a file:

from h5 import HDFArchive
import triqs.utility.mpi as mpi
if mpi.is_master_node():
    ar = HDFArchive("YourDFTDMFTcalculation.h5",'w')
    ar["G"] = S.G_iw
    ar["Sigma"] = S.Sigma_iw

These are the essential steps necessary for a one-shot DFT+DMFT calculation. For a detailed description of the SumkDFT routines, see the reference manual. To perform full charge self-consistent calculations, there are some more things to consider, which we will see later on.

Restarting a calculation

Often only a few DMFT iterations are performed first, and thus, it is desirable to carry out further iterations, e.g. to improve on the convergence. With a little modification at the initialization stage (before the DMFT loop) it is possible to detect if previous runs are present, or if the calculation should start from scratch:

previous_runs = 0
previous_present = False
if mpi.is_master_node():
    with HDFArchive(dft_filename+'.h5','a') as f:
        if 'dmft_output' in f:
            ar = f['dmft_output']
            if 'iterations' in ar:
                previous_present = True
                previous_runs = ar['iterations']
        else:
            f.create_group('dmft_output')

previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present)

You can see from this code snippet, that removing the subgroup dmft_results from the hdf file has the effect of reseting the calculation to the starting point. If there are previous runs stored in the hdf5 archive, we can now load the self energy, the chemical potential and double counting values of the last iteration:

if previous_present:
    if mpi.is_master_node():
        with HDFArchive(dft_filename+'.h5','r') as ar:
            S.Sigma_iw << ar['dmft_output']['Sigma_iw']

    S.Sigma_iw << mpi.bcast(S.Sigma_iw)
    chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
    SK.set_mu(chemical_potential)
    SK.set_dc(dc_imp,dc_energ)

The data is loaded only on the master node, and therefore we broadcast it to the slave nodes. Be careful when storing the iteration_number as we also have to add the previous iteration count:

ar['dmft_output']['iterations'] = iteration_number + previous_runs

Mixing

In some cases a mixing of two consecutive self energies (or alternatively two hybridization functions) can be necessary in order to ensure convergence:

mix = 0.8 # mixing factor
if (iteration_number>1 or previous_present):
    if mpi.is_master_node():
        with HDFArchive(dft_filename+'.h5','r') as ar:
            mpi.report("Mixing Sigma and G with factor %s"%mix)
            S.Sigma_iw << mix * S.Sigma_iw + (1.0-mix) * ar['dmft_output']['Sigma_iw']
            S.G_iw << mix * S.G_iw + (1.0-mix) * ar['dmft_output']['G_iw']
    S.G_iw << mpi.bcast(S.G_iw)
    S.Sigma_iw << mpi.bcast(S.Sigma_iw)

In this little piece of code, which should be placed after calling the solver, two consecutive self energies are linearly mixed with the factor mix. Of course, it is possible to implement more advanced mixing schemes (e.g. Broyden’s methods), however, in most cases simple linear mixing or even no mixing is sufficient for a reasonably fast convergence.