[doc] corrections to tutorials

doc/Ce-HI.rst

@@ -95,23 +95,16 @@ where the solver is initialized with the value of `beta`, and the orbital quantu
 The Hubbard-I initialization `Solver` has also optional parameters one may use:
 
   * `n_msb`: the number of Matsubara frequencies used. The default is `n_msb=1025`.
-  * `use_spin_orbit`: if set 'True' the solver is run with spin-orbit coupling
+  * `use_spin_orbit`: if set 'True' the solver is run with spin-orbit coupling included. To perform actual DFT+DMFT calculations with spin-orbit one should also run   :program:`Wien2k` and :program:`dmftproj` in spin-polarized mode and with spin-orbit included. By default, `use_spin_orbit=False`.
-included. To perform actual DFT+DMFT calculations with spin-orbit one should
+  * `Nmoments`: the number of moments used to describe high-ferquency tails of the Hubbard-I Green's function and self-energy. By default `Nmoments = 5`
-also run   :program:`Wien2k` and :program:`dmftproj` in spin-polarized mode and
-with spin-orbit included. By default, `use_spin_orbit=False`.
 
-The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the
+The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the Hubbard U parameter `U_int` and Hund's rule coupling `J_hund`. Notice that the solver constructs the full 4-index `U`-matrix by default, and the `U_int` parameter is in fact the Slatter `F0` integral. Other optional parameters are:
-Hubbard U parameter `U_int` and Hund's rule coupling `J_hund`. Notice that the
-solver constructs the full 4-index `U`-matrix by default, and the `U_int` parameter
-is in fact the Slatter `F0` integral. Other optional parameters are:
 
-  * `T`: matrix that transforms the interaction matrix from complex spherical
+  * `T`: matrix that transforms the interaction matrix from complex spherical harmonics to a symmetry adapted basis. By default, the complex spherical harmonics basis is used and `T=None`.
-harmonics to a symmetry adapted basis. By default, the complex spherical harmonics
+  * `verbosity`: tunes output from the solver. If `verbosity=0` only basic information is printed, if `verbosity=1` the ground state atomic occupancy and its energy are printed, if `verbosity=2` additional information is printed for all occupancies that were diagonalized. By default, `verbosity=0`. 
-basis is used and `T=None`.
+
-  * `verbosity`: tunes output from the solver. If `verbosity=0` only basic
+  * `Iteration_Number`: the iteration number of the DMFT loop. Used only for printing. By default   `Iteration_Number=1`
-information is printed, if `verbosity=1` the ground state atomic occupancy and
+  * `Test_Convergence`: convergence criterion. Once the self-energy is converged below `Test_Convergence` the Hubbard-I solver is not called anymore. By default `Test_Convergence=0.0001`.
-its energy are printed, if `verbosity=2` additional information is printed for
-all occupancies that were diagonalized. By default, `verbosity=0`.
 
 We need also to introduce some changes in the DMFT loop with respect that used for CT-QMC calculations in :ref:`advanced`. 
 The hybridization function is neglected in the Hubbard-I approximation, and only non-interacting level 
@@ -160,22 +153,23 @@ use here the default convergence criterion in :program:`Wien2k` (convergence to
 After calculations are done we may check the value of correlation ('Hubbard') energy correction to the total energy::
     
    >grep HUBBARD Ce-gamma.scf|tail -n 1
-   HUBBARD ENERGY(included in SUM OF EIGENVALUES):           -0.220502
+   HUBBARD ENERGY(included in SUM OF EIGENVALUES):           -0.012866
 
-and the band ("kinetic") energy with DMFT correction::
+In the case of Ce, with the correlated shell occupancy close to 1 the Hubbard energy is close to 0, while the DC correction to energy is about J/4 in accordance with the fully-localized-limit formula, hence, giving the total correction :math:`\Delta E_{HUB}=E_{HUB}-E_{DC} \approx -J/4`, which is in our case is equal to -0.175 eV :math:`\approx`-0.013 Ry.
+
+The band ("kinetic") energy with DMFT correction is ::
 
    >grep DMFT Ce-gamma.scf |tail -n 1
-   KINETIC ENERGY with DMFT correction:                      -5.329087
+   KINETIC ENERGY with DMFT correction:                      -5.370632
 
-as well as the convergence in total energy::
+One may also check the convergence in total energy::
    
    >grep :ENE Ce-gamma.scf |tail -n 5
-   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.77119670
+   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.56318334
-   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.77050935
+   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.56342250
-   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.77040176
+   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.56271503
-   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.77020712
+   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.56285812
-   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.77037540
+   :ENE  : ********** TOTAL ENERGY IN Ry =       -17717.56287381
-
 
 Calculating DOS with Hubbard-I
 ------------------------------
@@ -197,7 +191,7 @@ Then one needs to load projectors needed for calculations of corresponding proje
 
    SK = SumkDFTTools(hdf_file=dft_filename+'.h5', use_dft_blocks=False)
 
-Then after the solver initialization and setting up atomic levels we compute atomic Green's function and self-energy on the real axis::
+Then after the solver initialization we load the previously calculated chemical potential and double-counting correction.  Having set up atomic levels we then compute the atomic Green's function and self-energy on the real axis::
 
    S.set_atomic_levels( eal = eal )
    S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund)

doc/Ce-gamma.indmftpr

@@ -6,4 +6,3 @@ complex
 0 0 0 0          ! l included for each sort
 0
 -.40 0.40         ! Energy window relative to E_f 
-

doc/Ce-gamma.py

@@ -2,25 +2,23 @@ from pytriqs.applications.dft.sumk_dft import *
 from pytriqs.applications.dft.converters.wien2k_converter import *
 from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
 
-dft_filename = 'Ce-gamma'
+lda_filename = 'Ce-gamma'
 beta = 40
 U_int = 6.00
 J_hund = 0.70
-Loops =  2                       # Number of DMFT sc-loops
+Loops =  5                       # Number of DMFT sc-loops
 Mix = 0.7                        # Mixing factor in QMC
-                                 # 1.0 ... all from imp; 0.0 ... all from Gloc
 DC_type = 0                      # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
-useBlocs = False                 # use bloc structure from DFT input
+DC_Mix = 1.0                     # 1.0 ... all from imp; 0.0 ... all from Gloc   
+useBlocs = False                 # use bloc structure from LDA input
 useMatrix = True                 # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
-Natomic = 1
+chemical_potential_init=0.0      # initial chemical potential
 
-HDFfilename = dft_filename+'.h5'
+HDFfilename = lda_filename+'.h5'
-
-use_val= U_int * (Natomic - 0.5) - J_hund * (Natomic * 0.5 - 0.5)
 
 # Convert DMFT input:
 # Can be commented after the first run
-Converter = Wien2kConverter(filename=dft_filename)
+Converter = Wien2kConverter(filename=lda_filename)
 Converter.convert_dft_input()
 
 #check if there are previous runs:
@@ -42,23 +40,30 @@ previous_runs    = mpi.bcast(previous_runs)
 previous_present = mpi.bcast(previous_present)
 
 # Init the SumK class
-SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=False)
+SK=SumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
 
 Norb = SK.corr_shells[0]['dim']
 l    = SK.corr_shells[0]['l']
 
-# Init the Solver:
+# Init the Hubbard-I solver:
 S = Solver(beta = beta, l = l)
 
+chemical_potential=chemical_potential_init
+# load previous data: old self-energy, chemical potential, DC correction
 if (previous_present):
-    # load previous data:
     mpi.report("Using stored data for initialisation")
     if (mpi.is_master_node()):
         ar = HDFArchive(HDFfilename,'a')
-        S.Sigma << ar['SigmaImFreq']
+        S.Sigma <<= ar['SigmaF']
         del ar
+        things_to_load=['chemical_potential','dc_imp']
+        old_data=SK.load(things_to_load)
+        chemical_potential=old_data[0]
+        SK.dc_imp=old_data[1]
     S.Sigma = mpi.bcast(S.Sigma)
-    SK.load()
+    chemical_potential=mpi.bcast(chemical_potential)
+    SK.dc_imp=mpi.bcast(SK.dc_imp)
+
 
 # DMFT loop:
 for Iteration_Number in range(1,Loops+1):
@@ -69,84 +74,77 @@ for Iteration_Number in range(1,Loops+1):
         SK.put_Sigma(Sigma_imp = [ S.Sigma ])
 
         # Compute the SumK, possibly fixing mu by dichotomy
-        if SK.density_required and (Iteration_Number > 0):
+        if SK.density_required and (Iteration_Number > 1):
-            Chemical_potential = SK.calc_mu( precision = 0.01 )
+            chemical_potential = SK.calc_mu( precision = 0.000001 )
         else:
-            mpi.report("No adjustment of chemical potential\nTotal density  = %.3f"%SK.total_density(mu=Chemical_potential))
+            mpi.report("No adjustment of chemical potential\nTotal density  = %.3f"%SK.total_density(mu=chemical_potential))
 
         # Density:
-        S.G << SK.extract_G_loc()[0]
+        S.G <<= SK.extract_G_loc()[0]
         mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
-        dm = S.G.density()
 
+        # calculated DC at the first run to have reasonable initial non-interacting atomic level positions
         if ((Iteration_Number==1)and(previous_present==False)):
-	    SK.calc_dc( dens_mat=dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_val=use_val)
+            dc_value_init=U_int/2.0
+            dm=S.G.density()
+	    SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_dc_value=dc_value_init)
 
-        # set atomic levels:
+        # calculate non-interacting atomic level positions:
         eal = SK.eff_atomic_levels()[0]
         S.set_atomic_levels( eal = eal )
 
-        # update hdf5
-        if (mpi.is_master_node()):
-            ar = HDFArchive(HDFfilename,'a')
-            ar['Chemical_Potential%s'%itn] = Chemical_potential
-            del ar
-
         # solve it:
         S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1)
 
-        if (mpi.is_master_node()):
-            ar = HDFArchive(HDFfilename)
-            ar['iterations'] = itn
- 
         # Now mix Sigma and G:
         if ((itn>1)or(previous_present)):
             if (mpi.is_master_node()and (Mix<1.0)):
+                ar = HDFArchive(HDFfilename,'r')
                 mpi.report("Mixing Sigma and G with factor %s"%Mix)
-                if ('SigmaImFreq' in ar):
+                if ('SigmaF' in ar):
-                    S.Sigma << Mix * S.Sigma + (1.0-Mix) * ar['SigmaImFreq']
+                    S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF']
                 if ('GF' in ar):
-                    S.G << Mix * S.G + (1.0-Mix) * ar['GF']
+                    S.G <<= Mix * S.G + (1.0-Mix) * ar['GF']
-
+                del ar
             S.G = mpi.bcast(S.G)
             S.Sigma = mpi.bcast(S.Sigma)
-
-
-        
-        if (mpi.is_master_node()):
-            ar['SigmaImFreq'] = S.Sigma
-            ar['GF'] = S.G
         
         # after the Solver has finished, set new double counting: 
         dm = S.G.density()
-        SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type , use_val=use_val)
+        SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type )
+
         # correlation energy calculations:
         correnerg = 0.5 * (S.G * S.Sigma).total_density()
         mpi.report("Corr. energy = %s"%correnerg)
+
+        # store the impurity self-energy, GF as well as correlation energy in h5
         if (mpi.is_master_node()):
+            ar = HDFArchive(HDFfilename,'a')
+            ar['iterations'] = itn
+            ar['chemical_cotential%s'%itn] = chemical_potential
+            ar['SigmaF'] = S.Sigma
+            ar['GF'] = S.G
             ar['correnerg%s'%itn] = correnerg
             ar['DCenerg%s'%itn] = SK.dc_energ
             del ar
 
-
+        #Save essential SumkDFT data:
-        #Save stuff:
+        things_to_save=['chemical_potential','dc_energ','dc_imp']
-        SK.save(['chemical_potential','dc_imp','dc_energ'])
+        SK.save(things_to_save)
         if (mpi.is_master_node()):
-            print 'DC after solver: ',SK.dc_imp[SK.invshellmap[0]]
+            print 'DC after solver: ',SK.dc_imp[0]
 
-
+        # print out occupancy matrix of Ce 4f
-        # do some analysis:
         mpi.report("Orbital densities of impurity Green function:")
-        dm1 = S.G.density()
+        for s in dm:
-        for s in dm1:
             mpi.report("Block %s: "%s)
-            for ii in range(len(dm1[s])):
+            for ii in range(len(dm[s])):
                 str = ''
-                for jj in range(len(dm1[s])):
+                for jj in range(len(dm[s])):
-                    if (dm1[s][ii,jj].real>0):
+                    if (dm[s][ii,jj].real>0):
-                        str += "   %.4f"%(dm1[s][ii,jj].real)
+                        str += "   %.4f"%(dm[s][ii,jj].real)
                     else:
-                        str += "  %.4f"%(dm1[s][ii,jj].real)
+                        str += "  %.4f"%(dm[s][ii,jj].real)
                 mpi.report(str)
         mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
 
@@ -154,11 +152,13 @@ for Iteration_Number in range(1,Loops+1):
 # find exact chemical potential
 if (SK.density_required):
     SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
-dN,d = SK.calc_density_correction(filename = dft_filename+'.qdmft')
+
+# calculate and save occupancy matrix in the Bloch basis for Wien2k charge denity recalculation
+dN,d = SK.calc_density_correction(filename = lda_filename+'.qdmft')
 
 mpi.report("Trace of Density Matrix: %s"%d)
 
-#correlation energy:
+# store correlation energy contribution to be read by Wien2ki and then included to DFT+DMFT total energy
 if (mpi.is_master_node()):
     ar = HDFArchive(HDFfilename)
     itn = ar['iterations'] 
@@ -166,6 +166,6 @@ if (mpi.is_master_node()):
     DCenerg = ar['DCenerg%s'%itn]
     del ar
     correnerg -= DCenerg[0]
-    f=open(dft_filename+'.qdmft','a')
+    f=open(lda_filename+'.qdmft','a')
     f.write("%.16f\n"%correnerg)
     f.close()

doc/Ce-gamma_DOS.py

@@ -20,7 +20,6 @@ broadening = 0.02
 HDFfilename = dft_filename+'.h5'
 
 # Convert DMFT input:
-# Can be commented after the first run
 Converter = Wien2kConverter(filename=dft_filename,repacking=True)
 Converter.convert_dft_input()
 Converter.convert_parproj_input()
@@ -30,11 +29,11 @@ previous_runs = 0
 previous_present = False
 
 if mpi.is_master_node():
-    ar = HDFArchive(HDFfilename,'a')
+    ar = HDFArchive(HDFfilename)
     if 'iterations' in ar:
         previous_present = True
         previous_runs = ar['iterations']
-    else: 
+    else:
         previous_runs = 0
         previous_present = False
     del ar
@@ -43,27 +42,37 @@ mpi.barrier()
 previous_runs    = mpi.bcast(previous_runs)
 previous_present = mpi.bcast(previous_present)
 
-# if previous runs are present, no need for recalculating the bloc structure
-# It has to be commented, if you run this script for the first time, starting
-# from a converted h5 archive.
-
 # Init the SumK class
 SK = SumkDFTTools(hdf_file=dft_filename+'.h5',use_dft_blocks=False)
 
+# load old chemical potential and DC
+chemical_potential=0.0
+if mpi.is_master_node():
+    ar = HDFArchive(HDFfilename)
+    things_to_load=['chemical_potential','dc_imp']
+    old_data=SK.load(things_to_load)
+    chemical_potential=old_data[0]
+    SK.dc_imp=old_data[1]
+SK.chemical_potential=mpi.bcast(chemical_potential)
+SK.dc_imp=mpi.bcast(SK.dc_imp)
 
 if (mpi.is_master_node()):
-    print 'DC after reading SK: ',SK.dc_imp[SK.invshellmap[0]]
+    print 'DC after reading SK: ',SK.dc_imp[0]
 
 N = SK.corr_shells[0]['dim']
 l = SK.corr_shells[0]['l']
 
 # Init the Solver:
 S = Solver(beta = Beta, l = l)
-S.Nmoments= 8
 
 # set atomic levels:
 eal = SK.eff_atomic_levels()[0]
 S.set_atomic_levels( eal = eal )
+
+# Run the solver to get GF and self-energy on the real axis
 S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund)
 SK.put_Sigma(Sigma_imp = [S.Sigma])
+
+# compute DOS
 SK.dos_partial(broadening=broadening)
+