From 867418ce198eee7aeed445af259cc5e9ebc5f487 Mon Sep 17 00:00:00 2001 From: aichhorn Date: Thu, 9 Jan 2014 20:41:31 +0530 Subject: [PATCH] Updated DOC to version 1.0, bugfix in sumk_lda.py --- doc/Ce-HI.rst | 17 +++++++++-------- doc/Ce-gamma.py | 36 +++++++++++++++++------------------- doc/LDADMFTmain.rst | 20 ++++++++++++++------ doc/advanced.rst | 14 +++++++------- doc/analysis.rst | 2 +- doc/interface.rst | 2 +- python/sumk_lda.py | 24 ++++++++++++++---------- 7 files changed, 63 insertions(+), 52 deletions(-) diff --git a/doc/Ce-HI.rst b/doc/Ce-HI.rst index b9591423..7e1f187a 100644 --- a/doc/Ce-HI.rst +++ b/doc/Ce-HI.rst @@ -77,7 +77,7 @@ In order to run LDA+DMFT calculations within Hubbard-I we need the corresponding It is generally similar to the script for the case of DMFT calculations with the CT-QMC solver (see :ref:`advanced`), however there are also some differences. First, instead of *pytriqs.applications.dft.solver_multiband* we import Hubbard-I solver :: - from pytriqs.applications.impurity_solvers.hubbard_I.solver import Solver + from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver The Hubbard-I solver is very fast and we do not need to take into account the LDA blocking structure or use any approximation for the *U*-matrix :: @@ -87,18 +87,19 @@ The Hubbard-I solver is very fast and we do not need to take into account the LD We load and convert the :program:`dmftproj` output and initialize the *SumkLDA* class as described in :ref:`LDADMFTmain` and :ref:`advanced` and then set up the Hubbard-I solver :: - S = Solver(beta = beta, U_int = U_int, J_hund = J_hund, l = l) - S.Nmoments=10 + S = Solver(beta = beta, l = l) -where the solver is initialized with the value of `beta` as well as the `U` parameter (`U_int`) and Hund's rule coupling `J_hund`. Notice that `Solver_Hubbard-I` constructs the full 4-index `U`-matrix by default, and the `U` parameter is in fact the Slatter `F0` integral. -The last necessary parameter is the orbital quantum number `l` (equal to 3 in our case). -The next line gives the number of self-energy momenta used to compute contribution from the high-frequency tails. +where the solver is initialized with the value of `beta`, and the orbital quantum number `l` (equal to 3 in our case). -The Hubbard-I solver initialization `Solver` has also several optional parameters one may use: + +The Hubbard-I initialization `Solver` has also optional parameters one may use: * `n_msb`: is the number of Matsubara frequencies used (default is `n_msb=1025`) - * `T`: A matrix that transforms the interaction matrix from complex spherical harmonics to a symmetry adapted basis. By default complex spherical harmonics basis is used and `T=None` * `use_spin_orbit`: if set 'True' the solver is run with spin-orbit coupling included. To perform actual LDA+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` + +The `Solver.solve(U_int, J_hund)` statement has two necessary parameters, the `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` parameter is in fact the Slatter `F0` integral. Other optional parameters are: + + * `T`: A matrix that transforms the interaction matrix from complex spherical harmonics to a symmetry adapted basis. By default complex spherical harmonics basis is used and `T=None` * `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` We need also to introduce some changes in the DMFT loop with respect to the ones used for CT-QMC calculations in :ref:`advanced`. diff --git a/doc/Ce-gamma.py b/doc/Ce-gamma.py index bd40fcf9..8586b765 100644 --- a/doc/Ce-gamma.py +++ b/doc/Ce-gamma.py @@ -1,12 +1,12 @@ from pytriqs.applications.dft.sumk_lda 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 -LDAFilename = 'Ce-gamma' +LDAFilename = 'Ce' Beta = 40 Uint = 6.00 JHund = 0.70 -Loops = 3 # Number of DMFT sc-loops +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 @@ -17,7 +17,7 @@ HDFfilename = LDAFilename+'.h5' # Convert DMFT input: # Can be commented after the first run -Converter = Wien2kConverter(filename=LDAFilename,repacking=True) +Converter = Wien2kConverter(filename=LDAFilename) Converter.convert_dmft_input() #check if there are previous runs: @@ -45,8 +45,7 @@ Norb = SK.corr_shells[0][3] l = SK.corr_shells[0][2] # Init the Solver: -S = Solver(Beta = Beta, Uint = Uint, JHund = JHund, l = l, Verbosity=2) -S.Nmoments=10 +S = Solver(beta = Beta, l = l) if (previous_present): # load previous data: @@ -64,21 +63,21 @@ for Iteration_Number in range(1,Loops+1): itn = Iteration_Number + previous_runs # put Sigma into the SumK class: - SK.put_Sigma(Sigmaimp = [ S.Sigma ]) + SK.put_Sigma(Sigma_imp = [ S.Sigma ]) # Compute the SumK, possibly fixing mu by dichotomy - if SK.Density_Required and (Iteration_Number > 0): - Chemical_potential = SK.find_mu( precision = 0.000001 ) + if SK.density_required and (Iteration_Number > 0): + Chemical_potential = SK.find_mu( precision = 0.01 ) else: mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=Chemical_potential)) # Density: - S.G <<= SK.extract_Gloc()[0] + S.G <<= SK.extract_G_loc()[0] mpi.report("Total charge of Gloc : %.6f"%S.G.total_density()) dm = S.G.density() if ((Iteration_Number==1)and(previous_present==False)): - SK.SetDoubleCounting( dm, U_interact = Uint, J_Hund = JHund, orb = 0, useDCformula = DC_type) + SK.set_dc( dens_mat=dm, U_interact = Uint, J_hund = JHund, orb = 0, use_dc_formula = DC_type) # set atomic levels: eal = SK.eff_atomic_levels()[0] @@ -91,7 +90,7 @@ for Iteration_Number in range(1,Loops+1): del ar # solve it: - S.Solve() + S.solve(U_int = Uint, J_hund = JHund, verbosity = 1) if (mpi.is_master_node()): ar = HDFArchive(HDFfilename) @@ -99,7 +98,7 @@ for Iteration_Number in range(1,Loops+1): # Now mix Sigma and G: if ((itn>1)or(previous_present)): - if (mpi.is_master_node()): + if (mpi.is_master_node()and (Mix<1.0)): mpi.report("Mixing Sigma and G with factor %s"%Mix) if ('SigmaF' in ar): S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF'] @@ -117,13 +116,13 @@ for Iteration_Number in range(1,Loops+1): # after the Solver has finished, set new double counting: dm = S.G.density() - SK.SetDoubleCounting( dm, U_interact = Uint, J_Hund = JHund, orb = 0, useDCformula = DC_type ) + SK.set_dc( dm, U_interact = Uint, J_hund = JHund, 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) if (mpi.is_master_node()): ar['correnerg%s'%itn] = correnerg - ar['DCenerg%s'%itn] = SK.DCenerg + ar['DCenerg%s'%itn] = SK.dc_energ del ar @@ -150,9 +149,9 @@ 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_DensityCorrection(Filename = LDAFilename+'.qdmft') +if (SK.density_required): + SK.chemical_potential = SK.find_mu( precision = 0.000001 ) +dN,d = SK.calc_density_correction(filename = LDAFilename+'.qdmft') mpi.report("Trace of Density Matrix: %s"%d) @@ -168,4 +167,3 @@ if (mpi.is_master_node()): f.write("%.16f\n"%correnerg) f.close() - diff --git a/doc/LDADMFTmain.rst b/doc/LDADMFTmain.rst index 9bbd6d3d..f7c5c4f7 100644 --- a/doc/LDADMFTmain.rst +++ b/doc/LDADMFTmain.rst @@ -41,15 +41,23 @@ Setting up the Multi-Band Solver There is a module that helps setting up the multiband CTQMC solver. It is loaded and initialized by:: from pytriqs.applications.dft.solver_multiband import * - S = SolverMultiBand(Beta, U_interact, J_Hund, Norb) + S = SolverMultiBand(beta, n_orb, gf_struct = SK.gf_struct_solver[0], map=SK.map[0]) The necessary parameters are the inverse temperature `beta`, the Coulomb interaction `U_interact`, the Hund's rule coupling `J_hund`, and the number of orbitals `n_orb`. There are again several optional parameters that allow to modify the local Hamiltonian to specific needs. They are: * `gf_struct`: Contains the block structure of the local density matrix. Has to be given in the format as calculated by :class:`SumkLDA`. - * `map`: If `gf_Struct` is given as parameter, also `map` has to be given. This is the mapping from the block structure to a general + * `map`: If `gf_struct` is given as parameter, also `map` has to be given. This is the mapping from the block structure to a general up/down structure. + +The solver method is called later by this statement:: + + S.solve(U_interact = U, J_hund = J) + +The parameters for the Coulomb interaction `U_interact` and the Hunds coupling `J_hund` are necessary parameters. +The following parameters are optional, by highly recommended to be set: + * `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 `U_interact+2J_hund`, the interorbital interaction for opposite spins as `U_interact`, and interorbital for equal spins as @@ -69,8 +77,8 @@ at the end of this tutorial. After initialisation, several other CTQMC parameters can be set (see CTQMC doc). The most important are: - * `S.N_Cycles`: Number of QMC cycles per node. - * `S.N_Warmup_Cycles`: Number of iterations used for thermalisation + * `S.n_cycles`: Number of QMC cycles per node. + * `S.n_warmup_cycles`: Number of iterations used for thermalisation @@ -91,7 +99,7 @@ set up the loop over DMFT iterations and the self-consistency condition:: S.G <<= SK.extract_G_loc()[0] # extract the local Green function S.G0 <<= inverse(S.Sigma + inverse(S.G)) # finally get G0, the input for the Solver - S.Solve() # now solve the impurity problem + S.Solve(U_interact = U, J_hund = J) # now solve the impurity problem dm = S.G.density() # density matrix of the impurity problem SK.set_dc( dm, U_interact = U, J_hund = J, use_dc_formula = 0) # Set the double counting term @@ -112,7 +120,7 @@ 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("single_site_bethe.h5",'w') + R = HDFArchive("YourLDADMFTcalculation.h5",'w') R["G"] = S.G R["Sigma"] = S.Sigma diff --git a/doc/advanced.rst b/doc/advanced.rst index 8f366186..7d113a26 100644 --- a/doc/advanced.rst +++ b/doc/advanced.rst @@ -64,16 +64,14 @@ The next step is to initialise the Solver:: Norb = SK.corr_shells[0][3] l = SK.corr_shells[0][2] - S = SolverMultiBand(beta=beta,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) + 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:: - S.N_Cycles = qmc_cycles - S.Length_Cycle = length_cycle - S.N_Warmup_Cycles = warming_iterations + 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 +119,9 @@ previous section, with some additional refinement:: S.G0 = mpi.bcast(S.G0) # Solve the impurity problem: - S.Solve() + 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)) # solution done, do the post-processing: mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density()) diff --git a/doc/analysis.rst b/doc/analysis.rst index f03b96a3..18548e2a 100644 --- a/doc/analysis.rst +++ b/doc/analysis.rst @@ -43,7 +43,7 @@ density of states of the Wannier orbitals, you simply type:: SK.check_input_dos(om_min, om_max, n_om) -which produces plots between real frequencies `ommin` and `ommax`, using a mesh of `N_om` points. There +which produces plots between real frequencies `om_min` and `om_max`, using a mesh of `n_om` points. There is an optional parameter, `broadening`, which defines an additional Lorentzian broadening, and is set to `0.01` by default. diff --git a/doc/interface.rst b/doc/interface.rst index ae9eaad1..21759dfb 100644 --- a/doc/interface.rst +++ b/doc/interface.rst @@ -5,7 +5,7 @@ The interface The basic function of the interface to the Wien2k program package is to take the output of the program that constructs the projected local -orbitals (:program:`dmftproj`), and to store all the necessary information into +orbitals (:program:`dmftproj`, for documentation see :download:`TutorialDmftproj.pdf `), and to store all the necessary information into an hdf5 file. This latter file is then used to do the DMFT calculation. The reason for this structure is that this enables the user to have everything that is necessary to reproduce the calculation in one single hdf5 arxive. diff --git a/python/sumk_lda.py b/python/sumk_lda.py index 8401e478..77091304 100644 --- a/python/sumk_lda.py +++ b/python/sumk_lda.py @@ -654,13 +654,13 @@ class SumkLDA: #if (not hasattr(self,"dc_imp")): self.__init_dc() - - - dm = [ {} for i in xrange(self.n_corr_shells)] - for i in xrange(self.n_corr_shells): - l = self.corr_shells[i][3] #*(1+self.corr_shells[i][4]) - for j in xrange(len(self.gf_struct_corr[i])): - dm[i]['%s'%self.gf_struct_corr[i][j][0]] = numpy.zeros([l,l],numpy.float_) + + + #dm = [ {} for i in xrange(self.n_corr_shells)] + #for i in xrange(self.n_corr_shells): + # l = self.corr_shells[i][3] #*(1+self.corr_shells[i][4]) + # for j in xrange(len(self.gf_struct_corr[i])): + # dm[i]['%s'%self.gf_struct_corr[i][j][0]] = numpy.zeros([l,l],numpy.float_) for icrsh in xrange(self.n_corr_shells): @@ -676,10 +676,14 @@ class SumkLDA: self.dc_imp[icrsh]['%s'%self.gf_struct_corr[icrsh][j][0]] = numpy.identity(l,numpy.float_) blname = self.gf_struct_corr[icrsh][j][0] Ncr[blname] = 0.0 - for bl in self.map[iorb][blname]: - Ncr[blname] += dens_mat[bl].real.trace() - + for a,al in self.gf_struct_solver[iorb]: + #for bl in self.map[iorb][blname]: + bl = self.map_inv[iorb][a] + #print 'bl, valiue = ',bl,dens_mat[a].real.trace() + Ncr[bl] += dens_mat[a].real.trace() + + #print 'Ncr=',Ncr M = self.corr_shells[icrsh][3] Ncrtot = 0.0