diff --git a/doc/guide/dftdmft_singleshot.rst b/doc/guide/dftdmft_singleshot.rst
index 51330c3f..02d68466 100644
--- a/doc/guide/dftdmft_singleshot.rst
+++ b/doc/guide/dftdmft_singleshot.rst
@@ -5,11 +5,10 @@
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
-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
---------------------------------
@@ -25,15 +24,27 @@ to get the local quantities used in DMFT. It is initialized by::
Setting up the impurity solver
------------------------------
-The next step is to setup the impurity solver.
-For more details here, see the `CTHYB `_ documentation.
+The next step is to setup an impurity solver. There are different
+solvers available within the TRIQS framework. Below, we will discuss
+the example of the hybridisation
+expansion :ref:`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 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
-------------------
-Having initialised 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
+Having initialized the SumK class and the Solver, we can proceed with the DMFT
+loop itself. We have to set up the loop over DMFT
iterations and the self-consistency condition::
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
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
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
-the self-consistency steps, the solution of the Anderson impurity problem is
-calculation by CTQMC. Different to model calculations, we have to do a few
+description of the :class:`SumkDFT` routines, see the reference
+manual.
+
+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
calculate the density of the impurity problem. Then, the routine `calc_dc`
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.
-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
+`) and one with a
+rotational-invariant Slater interaction hamiltonian (:download:`dft_dmft_cthyb_slater.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::
from pytriqs.applications.dft.sumk_dft import *
@@ -94,16 +120,17 @@ First, we load the necessary modules::
from pytriqs.operators.util 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::
- dft_filename='srvo3'
- U = 2.7
+ dft_filename='SrVO3'
+ U = 4.0
J = 0.65
beta = 40
loops = 10 # Number of DMFT sc-loops
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
use_blocks = True # use bloc structure from DFT input
prec_mu = 0.0001
@@ -114,7 +141,11 @@ Then we define some parameters::
p["n_warmup_cycles"] = 2000
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 ` documentation.
+
The next step, as described in the previous section, is to convert the input files::
Converter = Wien2kConverter(filename=dft_filename, repacking=True)
@@ -140,24 +171,57 @@ from scratch::
previous_runs = mpi.bcast(previous_runs)
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::
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 ` library::
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
+ # Use GF structure determined by DFT blocks:
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)
- # 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 ` 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)
+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 ` library
+reference manual.
+
+
If there are previous runs stored in the hdf5 archive, we can now load the self energy
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
of the basic steps discussed in the previous section, with some additional
-refinement::
+refinements::
for iteration_number in range(1,loops+1):
if mpi.is_master_node(): print "Iteration = ", iteration_number
@@ -192,24 +256,13 @@ refinement::
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)
+ S.G0_iw << S.Sigma_iw + inverse(S.G_iw)
+ S.G0_iw << inverse(S.G0_iw)
# Solve the impurity problem:
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())
# Now mix Sigma and G with factor sigma_mix, if wanted:
diff --git a/doc/guide/images_scripts/dft_dmft_cthyb.py b/doc/guide/images_scripts/dft_dmft_cthyb.py
new file mode 100644
index 00000000..808a31c9
--- /dev/null
+++ b/doc/guide/images_scripts/dft_dmft_cthyb.py
@@ -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
diff --git a/doc/guide/images_scripts/dft_dmft_cthyb_slater.py b/doc/guide/images_scripts/dft_dmft_cthyb_slater.py
new file mode 100644
index 00000000..8ed8ccc4
--- /dev/null
+++ b/doc/guide/images_scripts/dft_dmft_cthyb_slater.py
@@ -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