LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-10-31 11:13:46 +01:00
Code Issues Releases Wiki Activity

Updated DOC to version 1.0, bugfix in sumk_lda.py

Browse Source
This commit is contained in:
aichhorn 2014-01-09 20:41:31 +05:30
parent b7dca200ec
commit 867418ce19
7 changed files with 63 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
doc/Ce-HI.rst
View File

@ -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`), 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 :: 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 :: 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 :: 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 = Solver(beta = beta, l = l)
S.Nmoments=10
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. where the solver is initialized with the value of `beta`, and the orbital quantum number `l` (equal to 3 in our case).
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.
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`) * `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` * `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` * `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`. We need also to introduce some changes in the DMFT loop with respect to the ones used for CT-QMC calculations in :ref:`advanced`.

36
doc/Ce-gamma.py
View File

@ -1,12 +1,12 @@
from pytriqs.applications.dft.sumk_lda import * from pytriqs.applications.dft.sumk_lda import *
from pytriqs.applications.dft.converters.wien2k_converter 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 Beta = 40
Uint = 6.00 Uint = 6.00
JHund = 0.70 JHund = 0.70
Loops = 3 # Number of DMFT sc-loops Loops = 2 # Number of DMFT sc-loops
Mix = 0.7 # Mixing factor in QMC Mix = 0.7 # Mixing factor in QMC
DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein 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 DC_Mix = 1.0 # 1.0 ... all from imp; 0.0 ... all from Gloc
@ -17,7 +17,7 @@ HDFfilename = LDAFilename+'.h5'
# Convert DMFT input: # Convert DMFT input:
# Can be commented after the first run # Can be commented after the first run
Converter = Wien2kConverter(filename=LDAFilename,repacking=True) Converter = Wien2kConverter(filename=LDAFilename)
Converter.convert_dmft_input() Converter.convert_dmft_input()
#check if there are previous runs: #check if there are previous runs:
@ -45,8 +45,7 @@ Norb = SK.corr_shells[0][3]
l = SK.corr_shells[0][2] l = SK.corr_shells[0][2]
# Init the Solver: # Init the Solver:
S = Solver(Beta = Beta, Uint = Uint, JHund = JHund, l = l, Verbosity=2) S = Solver(beta = Beta, l = l)
S.Nmoments=10
if (previous_present): if (previous_present):
# load previous data: # load previous data:
@ -64,21 +63,21 @@ for Iteration_Number in range(1,Loops+1):
itn = Iteration_Number + previous_runs itn = Iteration_Number + previous_runs
# put Sigma into the SumK class: # 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 # Compute the SumK, possibly fixing mu by dichotomy
if SK.Density_Required and (Iteration_Number > 0): if SK.density_required and (Iteration_Number > 0):
Chemical_potential = SK.find_mu( precision = 0.000001 ) Chemical_potential = SK.find_mu( precision = 0.01 )
else: 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: # 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()) mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
dm = S.G.density() dm = S.G.density()
if ((Iteration_Number==1)and(previous_present==False)): 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: # set atomic levels:
eal = SK.eff_atomic_levels()[0] eal = SK.eff_atomic_levels()[0]
@ -91,7 +90,7 @@ for Iteration_Number in range(1,Loops+1):
del ar del ar
# solve it: # solve it:
S.Solve() S.solve(U_int = Uint, J_hund = JHund, verbosity = 1)
if (mpi.is_master_node()): if (mpi.is_master_node()):
ar = HDFArchive(HDFfilename) ar = HDFArchive(HDFfilename)
@ -99,7 +98,7 @@ for Iteration_Number in range(1,Loops+1):
# Now mix Sigma and G: # Now mix Sigma and G:
if ((itn>1)or(previous_present)): 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) mpi.report("Mixing Sigma and G with factor %s"%Mix)
if ('SigmaF' in ar): if ('SigmaF' in ar):
S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF'] 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: # after the Solver has finished, set new double counting:
dm = S.G.density() 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: # correlation energy calculations:
correnerg = 0.5 * (S.G * S.Sigma).total_density() correnerg = 0.5 * (S.G * S.Sigma).total_density()
mpi.report("Corr. energy = %s"%correnerg) mpi.report("Corr. energy = %s"%correnerg)
if (mpi.is_master_node()): if (mpi.is_master_node()):
ar['correnerg%s'%itn] = correnerg ar['correnerg%s'%itn] = correnerg
ar['DCenerg%s'%itn] = SK.DCenerg ar['DCenerg%s'%itn] = SK.dc_energ
del ar del ar
@ -150,9 +149,9 @@ for Iteration_Number in range(1,Loops+1):
# find exact chemical potential # find exact chemical potential
if (SK.Density_Required): if (SK.density_required):
SK.Chemical_potential = SK.find_mu( precision = 0.000001 ) SK.chemical_potential = SK.find_mu( precision = 0.000001 )
dN,d = SK.calc_DensityCorrection(Filename = LDAFilename+'.qdmft') dN,d = SK.calc_density_correction(filename = LDAFilename+'.qdmft')
mpi.report("Trace of Density Matrix: %s"%d) mpi.report("Trace of Density Matrix: %s"%d)
@ -168,4 +167,3 @@ if (mpi.is_master_node()):
f.write("%.16f\n"%correnerg) f.write("%.16f\n"%correnerg)
f.close() f.close()

20
doc/LDADMFTmain.rst
View File

@ -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:: 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 * 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`, 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 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: 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`. * `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. 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 * `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 `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 `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: 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_cycles`: Number of QMC cycles per node.
* `S.N_Warmup_Cycles`: Number of iterations used for thermalisation * `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.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.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 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 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 from pytriqs.archive import HDFArchive
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
if mpi.is_master_node(): if mpi.is_master_node():
R = HDFArchive("single_site_bethe.h5",'w') R = HDFArchive("YourLDADMFTcalculation.h5",'w')
R["G"] = S.G R["G"] = S.G
R["Sigma"] = S.Sigma R["Sigma"] = S.Sigma

14
doc/advanced.rst
View File

@ -64,16 +64,14 @@ The next step is to initialise the Solver::
Norb = SK.corr_shells[0][3] Norb = SK.corr_shells[0][3]
l = SK.corr_shells[0][2] l = SK.corr_shells[0][2]
S = SolverMultiBand(beta=beta,U_interact=U,J_hund=J,n_orb=Norb,use_matrix=use_matrix, S = SolverMultiBand(beta=beta,n_orb=Norb,gf_struct=SK.gf_struct_solver[0],map=SK.map[0])
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)
As we can see, many options of the solver are set by properties of the :class:`SumkLDA` class, so we don't have 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. We now set the basic parameters of the QMC solver::
S.N_Cycles = qmc_cycles S.n_cycles = qmc_cycles
S.Length_Cycle = length_cycle S.length_cycle = length_cycle
S.N_Warmup_Cycles = warming_iterations S.n_warmup_cycles = warming_iterations
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::
@ -121,7 +119,9 @@ previous section, with some additional refinement::
S.G0 = mpi.bcast(S.G0) S.G0 = mpi.bcast(S.G0)
# Solve the impurity problem: # 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: # solution done, do the post-processing:
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density()) mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())

2
doc/analysis.rst
View File

@ -43,7 +43,7 @@ density of states of the Wannier orbitals, you simply type::
SK.check_input_dos(om_min, om_max, n_om) 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` is an optional parameter, `broadening`, which defines an additional Lorentzian broadening, and is set to `0.01`
by default. by default.

2
doc/interface.rst
View File

@ -5,7 +5,7 @@ The interface
The basic function of the interface to the Wien2k program package is The basic function of the interface to the Wien2k program package is
to take the output of the program that constructs the projected local 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 <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 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 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. that is necessary to reproduce the calculation in one single hdf5 arxive.

24
python/sumk_lda.py
View File

@ -654,13 +654,13 @@ class SumkLDA:
#if (not hasattr(self,"dc_imp")): self.__init_dc() #if (not hasattr(self,"dc_imp")): self.__init_dc()
dm = [ {} for i in xrange(self.n_corr_shells)] #dm = [ {} for i in xrange(self.n_corr_shells)]
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]) # l = self.corr_shells[i][3] #*(1+self.corr_shells[i][4])
for j in xrange(len(self.gf_struct_corr[i])): # 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[i]['%s'%self.gf_struct_corr[i][j][0]] = numpy.zeros([l,l],numpy.float_)
for icrsh in xrange(self.n_corr_shells): 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_) 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] blname = self.gf_struct_corr[icrsh][j][0]
Ncr[blname] = 0.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] M = self.corr_shells[icrsh][3]
Ncrtot = 0.0 Ncrtot = 0.0