Documentation and tutorial update

doc/Ce-HI.rst

@@ -178,7 +178,7 @@ Then one needs to load projectors needed for calculations of corresponding proje
 Then after the solver initialization and setting up atomic levels we compute 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)
+   S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund)
 
 put it into SK class and then calculated the actual DOS::
 

doc/Ce-gamma.py

@@ -2,22 +2,25 @@ from pytriqs.applications.dft.sumk_lda import *
 from pytriqs.applications.dft.converters.wien2k_converter import *
 from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
 
-LDAFilename = 'Ce'
+lda_filename = 'Ce-gamma'
-Beta = 40
+beta = 40
-Uint = 6.00
+U_int = 6.00
-JHund = 0.70
+J_hund = 0.70
 Loops =  2                       # Number of DMFT sc-loops
 Mix = 0.7                        # Mixing factor in QMC
 DC_type = 0                      # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
 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
 
-HDFfilename = LDAFilename+'.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=LDAFilename)
+Converter = Wien2kConverter(filename=lda_filename)
 Converter.convert_dmft_input()
 
 #check if there are previous runs:
@@ -39,13 +42,13 @@ previous_runs    = mpi.bcast(previous_runs)
 previous_present = mpi.bcast(previous_present)
 
 # Init the SumK class
-SK=SumkLDA(hdf_file=LDAFilename+'.h5',use_lda_blocks=False)
+SK=SumkLDA(hdf_file=lda_filename+'.h5',use_lda_blocks=False)
 
 Norb = SK.corr_shells[0][3]
 l    = SK.corr_shells[0][2]
 
 # Init the Solver:
-S = Solver(beta = Beta, l = l)
+S = Solver(beta = beta, l = l)
 
 if (previous_present):
     # load previous data:
@@ -77,7 +80,7 @@ for Iteration_Number in range(1,Loops+1):
         dm = S.G.density()
 
         if ((Iteration_Number==1)and(previous_present==False)):
-	    SK.set_dc( dens_mat=dm, U_interact = Uint, J_hund = JHund, orb = 0, use_dc_formula = DC_type)
+	    SK.set_dc( dens_mat=dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_val=use_val)
 
         # set atomic levels:
         eal = SK.eff_atomic_levels()[0]
@@ -90,7 +93,7 @@ for Iteration_Number in range(1,Loops+1):
             del ar
 
         # solve it:
-        S.solve(U_int = Uint, J_hund = JHund, verbosity = 1)
+        S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1)
 
         if (mpi.is_master_node()):
             ar = HDFArchive(HDFfilename)
@@ -116,7 +119,7 @@ for Iteration_Number in range(1,Loops+1):
         
         # after the Solver has finished, set new double counting: 
         dm = S.G.density()
-        SK.set_dc( dm, U_interact = Uint, J_hund = JHund, orb = 0, use_dc_formula = DC_type )
+        SK.set_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type , use_val=use_val)
         # correlation energy calculations:
         correnerg = 0.5 * (S.G * S.Sigma).total_density()
         mpi.report("Corr. energy = %s"%correnerg)
@@ -151,7 +154,7 @@ for Iteration_Number in range(1,Loops+1):
 # find exact chemical potential
 if (SK.density_required):
     SK.chemical_potential = SK.find_mu( precision = 0.000001 )
-dN,d = SK.calc_density_correction(filename = LDAFilename+'.qdmft')
+dN,d = SK.calc_density_correction(filename = lda_filename+'.qdmft')
 
 mpi.report("Trace of Density Matrix: %s"%d)
 
@@ -163,7 +166,6 @@ if (mpi.is_master_node()):
     DCenerg = ar['DCenerg%s'%itn]
     del ar
     correnerg -= DCenerg[0]
-    f=open(LDAFilename+'.qdmft','a')
+    f=open(lda_filename+'.qdmft','a')
     f.write("%.16f\n"%correnerg)
     f.close()
-

doc/Ce-gamma_DOS.py

@@ -1,13 +1,13 @@
 from pytriqs.applications.dft.sumk_lda_tools import *
 from pytriqs.applications.dft.converters.wien2k_converter import *
-from pytriqs.applications.impurity_solvers.hubbard_I.solver import Solver
+from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
 
 # Creates the data directory, cd into it:
 #Prepare_Run_Directory(DirectoryName = "Ce-Gamma") 
-LDAFilename = 'Ce-gamma'
+lda_filename = 'Ce-gamma'
 Beta =  40
-Uint = 6.00
+U_int = 6.00
-JHund = 0.70
+J_hund = 0.70
 DC_type = 0                      # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
 load_previous = True              # load previous results
 useBlocs = False                 # use bloc structure from LDA input
@@ -17,13 +17,13 @@ ommax=6.0
 N_om=2001
 broadening = 0.02
 
-HDFfilename = LDAFilename+'.h5'
+HDFfilename = lda_filename+'.h5'
 
 # Convert DMFT input:
 # Can be commented after the first run
-Converter = Wien2kConverter(filename=LDAFilename,repacking=True)
+Converter = Wien2kConverter(filename=lda_filename,repacking=True)
 Converter.convert_dmft_input()
-Converter.convert_par_proj_input()
+Converter.convert_parproj_input()
 
 #check if there are previous runs:
 previous_runs = 0
@@ -48,7 +48,7 @@ previous_present = mpi.bcast(previous_present)
 # from a converted h5 archive.
 
 # Init the SumK class
-SK = SumkLDATools(hdf_file=LDAFilename+'.h5',use_lda_blocks=False)
+SK = SumkLDATools(hdf_file=lda_filename+'.h5',use_lda_blocks=False)
 
 
 if (mpi.is_master_node()):
@@ -58,13 +58,12 @@ N = SK.corr_shells[0][3]
 l = SK.corr_shells[0][2]
 
 # Init the Solver:
-S = Solver(Beta = Beta, Uint = Uint, JHund = JHund, l = l)
+S = Solver(beta = Beta, l = l)
-S.Nmoments=10
+S.Nmoments= 8
 
 # set atomic levels:
 eal = SK.eff_atomic_levels()[0]
 S.set_atomic_levels( eal = eal )
-S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om)
+S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om,U_int=U_int,J_hund=J_hund)
-S.Sigma.save('S.Sigma')
+SK.put_Sigma(Sigma_imp = [S.Sigma])
-SK.put_Sigma(Sigmaimp = [S.Sigma])
 SK.dos_partial(broadening=broadening)

doc/LDADMFTmain.rst

@@ -53,10 +53,12 @@ specific needs. They are:
 
 The solver method is called later by this statement::
 
-  S.solve(U_interact = U, J_hund = J)
+  S.solve(U_interact,J_hund,use_spinflip=False,use_matrix=True,
+                   l=2,T=None, dim_reps=None, irep=None, deg_orbs=[],n_cycles =10000,
+                   length_cycle=200,n_warmup_cycles=1000)
 
-The parameters for the Coulomb interaction `U_interact` and the Hunds coupling `J_hund` are necessary parameters.
+The parameters for the Coulomb interaction `U_interact` and the Hunds coupling `J_hund` are necessary parameters. The rest are optional parameters, for which default values are set. 
-The following parameters are optional, by highly recommended to be set:
+They denerally should be reset for a given problem. Their meaning is as follows:
 
   * `use_matrix`: If `True`, the interaction matrix is calculated from Slater integrals, which are calculated from `U_interact` and 
     `J_hund`. Otherwise, a Kanamori representation is used. Attention: We define the intraorbital interaction as 
@@ -65,22 +67,20 @@ The following parameters are optional, by highly recommended to be set:
   * `T`: A matrix that transforms the interaction matrix from spherical harmonics, to a symmetry adapted basis. Only effective, if 
     `use_matrix=True`.
   * `l`: Orbital quantum number. Again, only effective for Slater parametrisation.
-  * `deg_shells`: A list that gives the degeneracies of the orbitals. It is used to set up a global move of the CTQMC solver.
+  * `deg_orbs`: A list that gives the degeneracies of the orbitals. It is used to set up a global move of the CTQMC solver.
   * `use_spinflip`: If `True`, the full rotationally-invariant interaction is used. Otherwise, only density-density terms are
     kept in the local Hamiltonian.
   * `dim_reps`: If only a subset of the full d-shell is used a correlated orbtials, one can specify here the dimensions of all the subspaces
     of the d-shell, i.e. t2g and eg. Only effective for Slater parametrisation.
   * `irep`: The index in the list `dim_reps` of the subset that is used. Only effective for Slater parametrisation.
+  * `n_cycles`: Number of CTQMC cycles (a sequence of moves followed by a measurement) per core. The default value of 10000 is the minimum, and generally should be incresed
+  * `length_cycle`: Number of CTQMC moves per one cycle
+  * `n_warmup_cycles`: Number of initial CTQMC cycles before measurements start. Usually of order of 10000, sometimes needs to be increased significantly.
 
 Most of above parameters can be taken directly from the :class:`SumkLDA` class, without defining them by hand. We will see a specific example 
 at the end of this tutorial.
 
-After initialisation, several other CTQMC parameters can be set (see CTQMC doc). The most important are:
+After initialisation, several other CTQMC parameters can be set (see CTQMC doc). 
-
-  * `S.n_cycles`: Number of QMC cycles per node.
-  * `S.n_warmup_cycles`: Number of iterations used for thermalisation
-
-
 
 
 .. index:: LDA+DMFT loop, one-shot calculation
@@ -120,9 +120,9 @@ At the end of the calculation, we can save the Greens function and self energy i
   from pytriqs.archive import HDFArchive
   import pytriqs.utility.mpi as mpi
   if mpi.is_master_node():
-      R = HDFArchive("YourLDADMFTcalculation.h5",'w')
+      ar = HDFArchive("YourLDADMFTcalculation.h5",'w')
-      R["G"] = S.G
+      ar["G"] = S.G
-      R["Sigma"] = S.Sigma
+      ar["Sigma"] = S.Sigma
 
 This is it! 
 

doc/advanced.rst

@@ -67,11 +67,7 @@ The next step is to initialise the Solver::
   S = SolverMultiBand(beta=beta,n_orb=Norb,gf_struct=SK.gf_struct_solver[0],map=SK.map[0])
 
 As we can see, many options of the solver are set by properties of the :class:`SumkLDA` class, so we don't have
-to set them manually. We now set the basic parameters of the QMC solver::
+to set them manually. 
-
-  S.n_cycles  = qmc_cycles
-  S.length_cycle = length_cycle
-  S.n_warmup_cycles = warming_iterations
 
 If there are previous runs stored in the hdf5 archive, we can now load the self energy
 of the last iteration::
@@ -121,7 +117,8 @@ previous section, with some additional refinement::
         # Solve the impurity problem:
         S.Solve(U_interact=U,J_hund=J,n_orb=Norb,use_matrix=use_matrix, 
                 T=SK.T[0], gf_struct=SK.gf_struct_solver[0],map=SK.map[0], 
-                l=l, deg_orbs=SK.deg_shells[0], use_spinflip=use_spinflip))
+                l=l, deg_orbs=SK.deg_shells[0], use_spinflip=use_spinflip,
+                n_cycles =qmc_cycles,length_cycle=length_cycle,n_warmup_cycles=warming_iterations)
 
         # solution done, do the post-processing:
         mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())

doc/analysis.rst

@@ -66,19 +66,28 @@ Routines with real-frequency self energy
 ----------------------------------------
 
 In order to plot data including correlation effects on the real axis, one has to provide the real frequency self energy. 
-Most conveniently, it is stored as a real frequency :class:`BlockGf` object in the hdf5 file. There is one important thing to
+Most conveniently, it is stored as a real frequency :class:`BlockGf` object in the hdf5 file::
-keep in mind. The real frequency self energy has to carry the note `ReFreq`::
 
-  SigmaReFreq.note = 'ReFreq'
+  ar = HDFArchive(filename+'.h5','a')
-
+  ar['SigmaReFreq'] = Sigma_real
-This tells the SumkLDA routines, that it is indeed a real frequency Greens function. Supposed you have your self energy now
-in the archive, you can type::
-
-  ar=HDFArchive(SK.hdf_file)
-  SK.put_Sigma([ ar['SigmaReFreq'] ])
   del ar
 
-This loads the self energy and puts it into the :class:`SumkLDA` class. The chemical potential as well as the double
+You may also store it in text files. If all blocks of your self energy are of dimension 1x1  you store them in `filename_(block)0.dat` files. Here `(block)` is a block name (`up`, `down`, or combined `ud`). In the case when you have matrix blocks, you store them in `(i)_(j).dat` files in the `filename_(block)` directory
+
+
+This self energy is loaded and put into the :class:`SumkLDA` class by the function:: 
+
+  SK.constr_Sigma_real_axis(filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0)
+
+where:
+ 
+  * `filename` is the file name of the hdf5 archive file or the `fname` pattern in text files names as described above.  
+  * `hdf=True` the real-axis self energy will be read from the hdf5 file, `hdf=False`: from the text files
+  * `hdf_dataset` the name of dataset where the self energy is stored in the hdf5 file
+  * `n_om` number of points in the real-axis mesh (used only if `hdf=False`)
+  
+
+The chemical potential as well as the double
 counting correction was already read in the initialisation process.
 
 With this self energy, we can do now::