LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-11-07 06:33:48 +01:00
Code Issues Releases Wiki Activity

new version of single-shot dft+dmft

Browse Source
This commit is contained in:
aichhorn 2015-08-13 16:50:48 +02:00
parent a7cc27cab3
commit 6fff56fe8d
3 changed files with 390 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

131
doc/guide/dftdmft_singleshot.rst
View File

@ -5,11 +5,10 @@
Single-shot DFT+DMFT Single-shot DFT+DMFT
==================== ====================
.. warning::
TO BE UPDATED!
After having set up the hdf5 archive, we can now do our DFT+DMFT calculation. It consists of After having set up the hdf5 archive, we can now do our DFT+DMFT calculation. It consists of
initialisation steps, and the actual DMFT self consistency loop. initialization steps, and the actual DMFT self-consistency loop, as is
discussed below.
Initialisation of the calculation Initialisation of the calculation
--------------------------------- ---------------------------------
@ -25,15 +24,27 @@ to get the local quantities used in DMFT. It is initialized by::
Setting up the impurity solver Setting up the impurity solver
------------------------------ ------------------------------
The next step is to setup the impurity solver. The next step is to setup an impurity solver. There are different
For more details here, see the `CTHYB <http://ipht.cea.fr/triqs/1.2/applications/cthyb/>`_ documentation. solvers available within the TRIQS framework. Below, we will discuss
the example of the hybridisation
expansion :ref:`CTHYB solver <triqscthyb:welcome>`. Later on, we will
see also the example of the Hubbard-I solver. They all have in common,
that they are called by a uniform command::
S.solve(params)
where `params` are the solver parameters and depend on the actual
solver that is used. Before going into the details of the solver, let
us discuss in the next section how to perform the DMFT loop using
the methods of :program:`dft_tools`, assuming that we have set up a
working solver instance.
Doing the DMFT loop Doing the DMFT loop
------------------- -------------------
Having initialised the SumK class and the Solver, we can proceed with the DMFT Having initialized the SumK class and the Solver, we can proceed with the DMFT
loop itself. As explained in the tutorial, we have to set up the loop over DMFT loop itself. We have to set up the loop over DMFT
iterations and the self-consistency condition:: iterations and the self-consistency condition::
n_loops = 5 n_loops = 5
@ -47,13 +58,18 @@ iterations and the self-consistency condition::
S.solve(h_int=h_int, **p) # now solve the impurity problem S.solve(h_int=h_int, **p) # now solve the impurity problem
dm = S.G_iw.density() # Density matrix of 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=dc_type) # Set the double counting term 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 SK.save(['chemical_potential','dc_imp','dc_energ']) # Save data in the hdf5 archive
These basic steps are enough to set up the basic DMFT Loop. For a detailed These basic steps are enough to set up the basic DMFT Loop. For a detailed
description of the :class:`SumkDFT` routines, see the reference manual. After description of the :class:`SumkDFT` routines, see the reference
the self-consistency steps, the solution of the Anderson impurity problem is manual.
calculation by CTQMC. Different to model calculations, we have to do a few
After
the self-consistency steps (extracting a new :math:`G^0(i\omega)`),
the Anderson impurity problem is solved.
Different to model calculations, we have to do a few
more steps after this, because of the double-counting correction. We first more steps after this, because of the double-counting correction. We first
calculate the density of the impurity problem. Then, the routine `calc_dc` calculate the density of the impurity problem. Then, the routine `calc_dc`
takes as parameters this density matrix, the Coulomb interaction, Hund's rule takes as parameters this density matrix, the Coulomb interaction, Hund's rule
@ -80,11 +96,21 @@ For full charge-self consistent calculations, there are some more things
to consider, which we will see later on. to consider, which we will see later on.
A more advanced example A full DFT+DMFT calculation
----------------------- ---------------------------
We will discuss now how to set up a full working calculation,
including setting up the CTHYB solver, and specifying some more parameters
in order to make the calculation more efficient. Here, we
will see a more advanced example, which is also suited for parallel
execution. For the convenience of the user, we provide also two
working python scripts in this documentation. One for a calculation
using Kanamori definitions (:download:`dft_dmft_cthyb.py
<images_scripts/dft_dmft_cthyb.py>`) and one with a
rotational-invariant Slater interaction hamiltonian (:download:`dft_dmft_cthyb_slater.py
<images_scripts/dft_dmft_cthyb.py>`). The user has to adapt these
scripts to his own needs.
Normally, one wants to adjust some more parameters in order to make the calculation more efficient. Here, we
will see a more advanced example, which is also suited for parallel execution.
First, we load the necessary modules:: First, we load the necessary modules::
from pytriqs.applications.dft.sumk_dft import * from pytriqs.applications.dft.sumk_dft import *
@ -94,16 +120,17 @@ First, we load the necessary modules::
from pytriqs.operators.util import * from pytriqs.operators.util import *
from pytriqs.applications.impurity_solvers.cthyb import * from pytriqs.applications.impurity_solvers.cthyb import *
The last two lines load the modules for the construction of the CTHYB
solver.
Then we define some parameters:: Then we define some parameters::
dft_filename='srvo3' dft_filename='SrVO3'
U = 2.7 U = 4.0
J = 0.65 J = 0.65
beta = 40 beta = 40
loops = 10 # Number of DMFT sc-loops loops = 10 # Number of DMFT sc-loops
sigma_mix = 0.8 # Mixing factor of Sigma after solution of the AIM sigma_mix = 0.8 # Mixing factor of Sigma after solution of the AIM
delta_mix = 1.0 # Mixing factor of Delta as input for the AIM
dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF
use_blocks = True # use bloc structure from DFT input use_blocks = True # use bloc structure from DFT input
prec_mu = 0.0001 prec_mu = 0.0001
@ -114,7 +141,11 @@ Then we define some parameters::
p["n_warmup_cycles"] = 2000 p["n_warmup_cycles"] = 2000
p["n_cycles"] = 20000 p["n_cycles"] = 20000
Most of these parameters are self-explanatory. The first, `dft_filename`, gives the filename of the input files. Most of these parameters are self-explanatory. The first,
`dft_filename`, gives the filename of the input files. For more
details on the solver parameters, we refer the user to
the :ref:`CTHYB solver <triqscthyb:welcome>` documentation.
The next step, as described in the previous section, is to convert the input files:: The next step, as described in the previous section, is to convert the input files::
Converter = Wien2kConverter(filename=dft_filename, repacking=True) Converter = Wien2kConverter(filename=dft_filename, repacking=True)
@ -140,24 +171,57 @@ from scratch::
previous_runs = mpi.bcast(previous_runs) previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present) previous_present = mpi.bcast(previous_present)
You can see in this code snipet, that all results of this calculation
will be stored in a separate subgroup in the hdf file, called
`dmft_output`. Removing this subgroup allows you to reset your
calculation to the starting point easily.
Now we can use all this information to initialise the :class:`SumkDFT` class:: Now we can use all this information to initialise the :class:`SumkDFT` class::
SK = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks) SK = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
The next step is to initialise the :class:`Solver` class:: The next step is to initialise the :class:`Solver` class. It consist
of two steps
#. Calculating the multi-band interaction matrix, and setting up the
interaction hamiltonian
#. Setting up the solver class
The first step is done using methods of
the :ref:`TRIQS <triqslibs:welcome>` library::
n_orb = SK.corr_shells[0]['dim'] n_orb = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l'] l = SK.corr_shells[0]['l']
spin_names = ["up","down"] spin_names = ["up","down"]
orb_names = [i for i in range(n_orb)] orb_names = [i for i in range(n_orb)]
# Use GF structure determined by DFT blocks # Use GF structure determined by DFT blocks:
gf_struct = SK.gf_struct_solver[0] gf_struct = SK.gf_struct_solver[0]
# Construct U matrix for density-density calculations # Construct U matrix for density-density calculations:
Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J) Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J)
# Construct Hamiltonian and solver
h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat, H_dump="H.txt") We assumed here that we want to use an interaction matrix with
Kanamori definitions of :math:`U` and :math:`J`. For
other choices (Slater interaction matrix for instance), and other
parameters, we refer to the reference manual
of the :ref:`TRIQS <triqslibs:welcome>` library.
Next, we construct the hamiltonian and the solver::
h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
S = Solver(beta=beta, gf_struct=gf_struct) S = Solver(beta=beta, gf_struct=gf_struct)
As you see, we take only density-density interactions into
account. Other choices for the hamiltonian are
* h_int_kanamori
* h_int_slater
These two include full rotational invariant interactions. Again,
options can be found in the :ref:`TRIQS <triqslibs:welcome>` library
reference manual.
If there are previous runs stored in the hdf5 archive, we can now load the self energy If there are previous runs stored in the hdf5 archive, we can now load the self energy
of the last iteration:: of the last iteration::
@ -174,7 +238,7 @@ last saved chemical potential and double counting values are read in and set.
Now we can go to the definition of the self-consistency step. It consists again Now we can go to the definition of the self-consistency step. It consists again
of the basic steps discussed in the previous section, with some additional of the basic steps discussed in the previous section, with some additional
refinement:: refinements::
for iteration_number in range(1,loops+1): for iteration_number in range(1,loops+1):
if mpi.is_master_node(): print "Iteration = ", iteration_number if mpi.is_master_node(): print "Iteration = ", iteration_number
@ -192,24 +256,13 @@ refinement::
S.Sigma_iw << SK.dc_imp[0]['up'][0,0] S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
# Calculate new G0_iw to input into the solver: # Calculate new G0_iw to input into the solver:
if mpi.is_master_node(): S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
# We can do a mixing of Delta in order to stabilize the DMFT iterations: S.G0_iw << inverse(S.G0_iw)
S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
if (iteration_number>1 or previous_present):
mpi.report("Mixing input Delta with factor %s"%delta_mix)
Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['Delta_iw']
S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta
ar['Delta_iw'] = delta(S.G0_iw)
S.G0_iw << inverse(S.G0_iw)
del ar
S.G0_iw << mpi.bcast(S.G0_iw)
# Solve the impurity problem: # Solve the impurity problem:
S.solve(h_int=h_int, **p) S.solve(h_int=h_int, **p)
# Solved. Now do post-processing: # Solved. Now do post-solution stuff:
mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density()) mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density())
# Now mix Sigma and G with factor sigma_mix, if wanted: # Now mix Sigma and G with factor sigma_mix, if wanted:

149
doc/guide/images_scripts/dft_dmft_cthyb.py Normal file
View File

@ -0,0 +1,149 @@
import pytriqs.utility.mpi as mpi
from pytriqs.operators.util import *
from pytriqs.archive import HDFArchive
from pytriqs.applications.impurity_solvers.cthyb import *
from pytriqs.gf.local import *
from pytriqs.applications.dft.sumk_dft import *
from pytriqs.applications.dft.converters.wien2k_converter import *
dft_filename='SrVO3'
U = U.0
J = 0.65
beta = 40
loops = 10 # Number of DMFT sc-loops
sigma_mix = 1.0 # Mixing factor of Sigma after solution of the AIM
delta_mix = 1.0 # Mixing factor of Delta as input for the AIM
dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF
use_blocks = True # use bloc structure from DFT input
prec_mu = 0.0001
h_field = 0.0
# Solver parameters
p = {}
p["max_time"] = -1
p["length_cycle"] = 50
p["n_warmup_cycles"] = 50
p["n_cycles"] = 5000
Converter = Wien2kConverter(filename=dft_filename, repacking=True)
Converter.convert_dft_input()
mpi.barrier()
previous_runs = 0
previous_present = False
if mpi.is_master_node():
f = HDFArchive(dft_filename+'.h5','a')
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')
del f
previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present)
SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field)
n_orb = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l']
spin_names = ["up","down"]
orb_names = [i for i in range(n_orb)]
# Use GF structure determined by DFT blocks
gf_struct = SK.gf_struct_solver[0]
# Construct U matrix for density-density calculations
Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J)
# Construct density-density Hamiltonian and solver
h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat, H_dump="H.txt")
S = Solver(beta=beta, gf_struct=gf_struct)
if previous_present:
chemical_potential = 0
dc_imp = 0
dc_energ = 0
if mpi.is_master_node():
S.Sigma_iw << HDFArchive(dft_filename+'.h5','a')['dmft_output']['Sigma_iw']
chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
chemical_potential = mpi.bcast(chemical_potential)
dc_imp = mpi.bcast(dc_imp)
dc_energ = mpi.bcast(dc_energ)
SK.set_mu(chemical_potential)
SK.set_dc(dc_imp,dc_energ)
for iteration_number in range(1,loops+1):
if mpi.is_master_node(): print "Iteration = ", iteration_number
SK.symm_deg_gf(S.Sigma_iw,orb=0) # symmetrise Sigma
SK.put_Sigma(Sigma_imp = [ S.Sigma_iw ]) # put Sigma into the SumK class
chemical_potential = SK.calc_mu( precision = prec_mu ) # find the chemical potential for given density
S.G_iw << SK.extract_G_loc()[0] # calc the local Green function
mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density())
# Init the DC term and the real part of Sigma, if no previous runs found:
if (iteration_number==1 and previous_present==False):
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
# Calculate new G0_iw to input into the solver:
if mpi.is_master_node():
# We can do a mixing of Delta in order to stabilize the DMFT iterations:
S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
if (iteration_number>1 or previous_present):
mpi.report("Mixing input Delta with factor %s"%delta_mix)
Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['Delta_iw']
S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta
ar['Delta_iw'] = delta(S.G0_iw)
S.G0_iw << inverse(S.G0_iw)
del ar
S.G0_iw << mpi.bcast(S.G0_iw)
# Solve the impurity problem:
S.solve(h_int=h_int, **p)
# Solved. Now do post-processing:
mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density())
# Now mix Sigma and G with factor sigma_mix, if wanted:
if (iteration_number>1 or previous_present):
if mpi.is_master_node():
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
mpi.report("Mixing Sigma and G with factor %s"%sigma_mix)
S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw']
S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw']
del ar
S.G_iw << mpi.bcast(S.G_iw)
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
# Write the final Sigma and G to the hdf5 archive:
if mpi.is_master_node():
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
if previous_runs: iteration_number += previous_runs
ar['iterations'] = iteration_number
ar['G_tau'] = S.G_tau
ar['G_iw'] = S.G_iw
ar['Sigma_iw'] = S.Sigma_iw
ar['G0-%s'%(iteration_number)] = S.G0_iw
ar['G-%s'%(iteration_number)] = S.G_iw
ar['Sigma-%s'%(iteration_number)] = S.Sigma_iw
del ar
# Set the new double counting:
dm = S.G_iw.density() # compute the density matrix of the impurity problem
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
# Save stuff into the dft_output group of hdf5 archive in case of rerun:
SK.save(['chemical_potential','dc_imp','dc_energ'])
if mpi.is_master_node():
ar = HDFArchive("dftdmft.h5",'w')
ar["G_tau"] = S.G_tau
ar["G_iw"] = S.G_iw
ar["Sigma_iw"] = S.Sigma_iw

149
doc/guide/images_scripts/dft_dmft_cthyb_slater.py Normal file
View File

@ -0,0 +1,149 @@
import pytriqs.utility.mpi as mpi
from pytriqs.operators.util import *
from pytriqs.archive import HDFArchive
from pytriqs.applications.impurity_solvers.cthyb import *
from pytriqs.gf.local import *
from pytriqs.applications.dft.sumk_dft import *
from pytriqs.applications.dft.converters.wien2k_converter import *
dft_filename='Gd_fcc'
U = 9.6
J = 0.8
beta = 40
loops = 10 # Number of DMFT sc-loops
sigma_mix = 1.0 # Mixing factor of Sigma after solution of the AIM
delta_mix = 1.0 # Mixing factor of Delta as input for the AIM
dc_type = 0 # DC type: 0 FLL, 1 Held, 2 AMF
use_blocks = True # use bloc structure from DFT input
prec_mu = 0.0001
h_field = 0.0
# Solver parameters
p = {}
p["max_time"] = -1
p["length_cycle"] = 50
p["n_warmup_cycles"] = 50
p["n_cycles"] = 5000
Converter = Wien2kConverter(filename=dft_filename, repacking=True)
Converter.convert_dft_input()
mpi.barrier()
previous_runs = 0
previous_present = False
if mpi.is_master_node():
f = HDFArchive(dft_filename+'.h5','a')
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')
del f
previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present)
SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field)
n_orb = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l']
spin_names = ["up","down"]
orb_names = [i for i in range(n_orb)]
# Use GF structure determined by DFT blocks
gf_struct = SK.gf_struct_solver[0]
# Construct Slater U matrix
Umat = U_matrix(n_orb=n_orb, U_int=U, J_hund=J, basis='cubic',)
# Construct Hamiltonian and solver
h_int = h_int_slater(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U_matrix=Umat)
S = Solver(beta=beta, gf_struct=gf_struct)
if previous_present:
chemical_potential = 0
dc_imp = 0
dc_energ = 0
if mpi.is_master_node():
S.Sigma_iw << HDFArchive(dft_filename+'.h5','a')['dmft_output']['Sigma_iw']
chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
chemical_potential = mpi.bcast(chemical_potential)
dc_imp = mpi.bcast(dc_imp)
dc_energ = mpi.bcast(dc_energ)
SK.set_mu(chemical_potential)
SK.set_dc(dc_imp,dc_energ)
for iteration_number in range(1,loops+1):
if mpi.is_master_node(): print "Iteration = ", iteration_number
SK.symm_deg_gf(S.Sigma_iw,orb=0) # symmetrise Sigma
SK.put_Sigma(Sigma_imp = [ S.Sigma_iw ]) # put Sigma into the SumK class
chemical_potential = SK.calc_mu( precision = prec_mu ) # find the chemical potential for given density
S.G_iw << SK.extract_G_loc()[0] # calc the local Green function
mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density())
# Init the DC term and the real part of Sigma, if no previous runs found:
if (iteration_number==1 and previous_present==False):
dm = S.G_iw.density()
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
# Calculate new G0_iw to input into the solver:
if mpi.is_master_node():
# We can do a mixing of Delta in order to stabilize the DMFT iterations:
S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
if (iteration_number>1 or previous_present):
mpi.report("Mixing input Delta with factor %s"%delta_mix)
Delta = (delta_mix * delta(S.G0_iw)) + (1.0-delta_mix) * ar['Delta_iw']
S.G0_iw << S.G0_iw + delta(S.G0_iw) - Delta
ar['Delta_iw'] = delta(S.G0_iw)
S.G0_iw << inverse(S.G0_iw)
del ar
S.G0_iw << mpi.bcast(S.G0_iw)
# Solve the impurity problem:
S.solve(h_int=h_int, **p)
# Solved. Now do post-processing:
mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density())
# Now mix Sigma and G with factor sigma_mix, if wanted:
if (iteration_number>1 or previous_present):
if mpi.is_master_node():
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
mpi.report("Mixing Sigma and G with factor %s"%sigma_mix)
S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['Sigma_iw']
S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['G_iw']
del ar
S.G_iw << mpi.bcast(S.G_iw)
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
# Write the final Sigma and G to the hdf5 archive:
if mpi.is_master_node():
ar = HDFArchive(dft_filename+'.h5','a')['dmft_output']
if previous_runs: iteration_number += previous_runs
ar['iterations'] = iteration_number
ar['G_tau'] = S.G_tau
ar['G_iw'] = S.G_iw
ar['Sigma_iw'] = S.Sigma_iw
ar['G0-%s'%(iteration_number)] = S.G0_iw
ar['G-%s'%(iteration_number)] = S.G_iw
ar['Sigma-%s'%(iteration_number)] = S.Sigma_iw
del ar
# Set the new double counting:
dm = S.G_iw.density() # compute the density matrix of the impurity problem
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
# Save stuff into the dft_output group of hdf5 archive in case of rerun:
SK.save(['chemical_potential','dc_imp','dc_energ'])
if mpi.is_master_node():
ar = HDFArchive("dftdmft.h5",'w')
ar["G_tau"] = S.G_tau
ar["G_iw"] = S.G_iw
ar["Sigma_iw"] = S.Sigma_iw