LCPQ
/
dft_tools
mirror of https://github.com/triqs/dft_tools
3
0
Fork 0
Code Issues Releases Wiki Activity

Add documentation

This commit is contained in:
Michel Ferrero 2013-08-07 16:40:18 +02:00
parent cd650884d3
commit 640941c326
18 changed files with 1382 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists.txt
View File

@ -21,3 +21,10 @@ add_subdirectory(fortran/dmftproj)
add_subdirectory(fortran/vertex)
add_subdirectory(python)
add_subdirectory(test)
option(BUILD_DOC "Build documentation" OFF)
if(${BUILD_DOC})
if(NOT TRIQS_WITH_DOCUMENTATION)
message("Error: TRIQS library has not been compiled with its documentation")
endif()
add_subdirectory(doc)
endif(${BUILD_DOC})

28
doc/CMakeLists.txt Normal file
View File

@ -0,0 +1,28 @@
# generate the conf.py
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in ${CMAKE_CURRENT_BINARY_DIR}/conf.py)
# all rst files of the documentation
file(GLOB_RECURSE doc_sources *.rst)
# create documentation target
set(sphinx_top ${CMAKE_CURRENT_BINARY_DIR}/html/index.html)
add_custom_command(OUTPUT ${sphinx_top} DEPENDS ${doc_sources} py_copy
COMMAND ${CMAKE_BINARY_DIR}/build_pytriqs ${TRIQS_SPHINXBUILD_EXECUTABLE} -c . -b html ${CMAKE_CURRENT_SOURCE_DIR} html)
add_custom_target(doc_sphinx ALL DEPENDS ${sphinx_top})
# install
install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/ COMPONENT documentation DESTINATION share/doc/dft_tools
FILES_MATCHING
PATTERN "*.html"
PATTERN "*.png"
PATTERN "*.js"
PATTERN "_*"
PATTERN "*.png"
PATTERN "*.gif"
PATTERN "*.xsl"
PATTERN "*.css"
PATTERN "*.pdf"
PATTERN "*.py"
PATTERN "*.txt"
PATTERN "*.bib"
)

203
doc/Ce-HI.rst Normal file
View File

@ -0,0 +1,203 @@
.. index:: tutorial on Ce within Hub.-I approximation
.. _LDADMFTtutorial:
LDA+DMFT tutorial: Ce with Hubbard-I approximation
==================================================
In this tutorial we will perform LDA+DMFT :program:`Wien2k` calculations of the high-temperature :math:`\gamma`-phase of Ce employing the
Hubbard-I approximation for its localized *4f* shell.
First we create the Wien2k :file:`Ce-gamma.struct` file as described in `Wien2k manual <http://www.wien2k.at/reg_user/textbooks/usersguide.pdf>`_
for the :math:`\gamma`-Ce fcc structure with lattice parameter of 9.75 a.u.
.. literalinclude:: Ce-gamma.struct
We initalize non-magnetic :program:`Wien2k` calculations using the :program:`init` script as described in the same manual.
For this example we specify 3000 *k*-points in the full Brillouin zone
and LDA exchange-correlation potential (*vxc=5*), other parameters are defaults.
The Ce *4f* electrons are treated as valence states.
Hence, the initialization script is executed as follows ::
init -b -vxc 5 -numk 3000
and then LDA calculations of non-magnetic :math:`\gamma`-Ce are performed by launching the :program:`Wien2k` :program:`run` script.
These self-consistent LDA calculations will typically take a couple of minutes.
Then we create :file:`Ce-gamma.indmftpr` file specifying parameters for construction of Wannier orbitals representing *4f* states:
.. literalinclude:: Ce-gamma.indmftpr
First three lines give the number of inequivalent sites, their multiplicity (to be in accordance with the *struct* file) and the maximum orbital quantum number :math:`l_{max}`
The following four lines describe the treatment of Ce *spdf* orbitals by the :program:`dmftproj` program ::
complex
1 1 1 2 ! l included for each sort
0 0 0 0 ! l included for each sort
0
where `complex` is the choice for the angular basis to be used (spherical complex harmonics), in the next line we specify, for each orbital
quantum number, whether it is treated as correlated ('2') and, hence, the corresponding Wannier orbitals will be generated or uncorrelated ('1').
In the latter case the :program:`dmftproj` program will generate projectors to be used in calculations of corresponding partial densities of states (see below).
In the present case we choose the fourth (i. e. *f*) orbitals as correlated.
The next line specify the number of irreducible representations into which a given correlated shell should be split (or
'0' if no splitting is desired, as in the present case). The fourth line specifies whether the spin-orbit interaction should be switched on ('1') or off ('0', as in the present case).
Finally, the last line iof the file ::
-.40 0.40 ! Energy window relative to E_f
specify the energy window for Wannier functions' construction. For a more complete description of :program:`dmftproj` options see its manual.
To prepaire input data for :program:`dmftproj` we execute :program:`lapw2` with the `-almd` option ::
lapw2 -almd
Then :program:`dmftproj` is executed in its default mode (i.e. without spin-polarization or spin-orbit included) ::
dmftproj
This program produces the following files:
* :file:`Ce-gamma.ctqmcout` and :file:`Ce-gamma.symqmc` containing projector operators and symmetry operations for orthonormalized Wannier orbitals, respectively.
* :file:`Ce-gamma.parproj` and :file:`Ce-gamma.sympar` containing projector operators and symmetry operations for uncorrelated states, respectively. These files are needed for projected density-of-states or spectral-function calculations.
* :file:`Ce-gamma.oubwin` needed for the charge desity recalculation in the case of fully self-consistent LDA+DMFT run (see below).
Now we have all necessary input from :program:`Wien2k` for running DMFT calculations.
.. index:: Hubbard-I in TRIQS
.. _HubITRIQS:
Hubbard-I calculations in TRIQS
-------------------------------
In order to run LDA+DMFT calculations within Hubbard-I we need the corresponding python script, :ref:`Ce-gamma-script`.
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
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 ::
use_blocks = False # use bloc structure from LDA input
use_matrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
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
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.
The Hubbard-I solver initialization `Solver` has also several 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`
* `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`.
The hybridization function is neglected in the Hubbard-I approximation, and only non-interacting level
positions (:math:`\hat{\epsilon}=-\mu+\langle H^{ff} \rangle - \Sigma_{DC}`) are required.
Hence, instead of computing `S.G0` as in :ref:`advanced` we set the level positions ::
# set atomic levels:
eal = SK.eff_atomic_levels()[0]
S.set_atomic_levels( eal = eal )
The part after the solution of the impurity problem remains essentially the same: we mix the self-energy and local
Green's function and then save them in the hdf5 file .
Then the double counting is recalculated and the correlation energy is computed with the Migdal formula and stored in hdf5.
Finally, we compute the modified charge density and save it as well as correlational correction to the total energy in
:file:`Ce-gamma.qdmft` file, which is then read by :program:`lapw2` in the case of self-consistent LDA+DMFT calculations.
.. index:: running LDA+DMFT calculations
Running LDA+DMFT calculations
-----------------------------
After having prepaired the script one may run one-shot DMFT calculations by executing :ref:`Ce-gamma-script` with :program:`pytriqs` in one-processor ::
pytriqs Ce-gamma.py
or parallel mode ::
mpirun pytriqs Ce-gamma.py
where :program:`mpirun` launches these calculations in parallel mode and enables MPI. The exact form of this command will, of course, depend on mpi-launcher installed in your system.
Instead of doing one-shot run one may also perform fully self-consistent LDA+DMFT calculations, as we will do in this tutorial. We launch these calculations as follows ::
run_para -qdmft
where `-qdmft` flag turns on LDA+DMFT calculations with :program:`Wien2k`. We use here the default convergence criterion in :program:`Wien2k` (convergence to 0.1 mRy in energy).
After calculations are done we may check the value of correlational ('Hubbard') energy correction to the total energy::
>grep HUBBARD Ce-gamma.scf|tail -n 1
HUBBARD ENERGY(included in SUM OF EIGENVALUES): -0.012875
and the band("kinetic") energy with DMFT correction::
>grep DMFT Ce-gamma.scf |tail -n 1
KINETIC ENERGY with DMFT correction: -3.714346
as well as the convergence in total energy::
>grep :ENE Ce-gamma.scf |tail -n 5
:ENE : ********** TOTAL ENERGY IN Ry = -17717.554865
:ENE : ********** TOTAL ENERGY IN Ry = -17717.554577
:ENE : ********** TOTAL ENERGY IN Ry = -17717.554420
:ENE : ********** TOTAL ENERGY IN Ry = -17717.554371
:ENE : ********** TOTAL ENERGY IN Ry = -17717.554419
Calculating DOS with Hubbard-I
------------------------------
Within Hubbard-I one may also easily obtain the spectral function ("band structure") and integrated spectral function ("density of states, DOS").
In difference with the CT-QMC approach one does not need to provide the real-frequency self-energy (see :ref:`analysis`) it can be calculated directly by the Hubbard-I solver.
The corresponding script :ref:`Ce-gamma_DOS-script` contains several new parameters ::
ommin=-4.0 # bottom of the energy range for DOS calculations
ommax=6.0 # top of the energy range for DOS calculations
N_om=2001 # number of points on the real-energy axis mesh
broadening = 0.02 # broadening (the imaginary shift of the real-energy mesh)
Then one needs to load projectors needed for calculations of corresponding projected densities of states, as well as corresponding symmetries. To get access to analysing tools we initialize the `SumkLDATools` class ::
SK = SumkLDATools(hdf_file=lda_filename+'.h5', use_lda_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::
S.set_atomic_levels( eal = eal )
S.GF_realomega(ommin=ommin, ommax = ommax, N_om=N_om)
put it into SK class and then calculated the actual DOS::
SK.dos_partial(broadening=broadening)
We may first increase the number of **k**-points in BZ to 10000 by executing :program:`Wien2k` program :program:`kgen` ::
x kgen
and then by executing :ref:`Ce-gamma_DOS-script` with :program:`pytriqs`::
pytriqs Ce-gamma_DOS.py
In result we get the total DOS for spins `up` and `down` (identical in our paramagnetic case) in :file:`DOScorrup.dat` and :file:`DOScorrdown.dat` files, respectively, as well as projected DOSs written in the corresponding files as described in :ref:`analysis`.
In our case, for example, the files :file:`DOScorrup.dat` and :file:`DOScorrup_proj3.dat` contain the total DOS for spin *up* and the corresponding projected DOS for Ce *4f* orbital, respectively. They are plotted below.
.. image:: Ce_DOS.png
:width: 700
:align: center
As one may clearly see, the Ce *4f* band is split by the local Coulomb interaction into the filled lower Hubbard band and empty upper Hubbard band (the latter is additionally split into several peaks due to the Hund's rule coupling and multiplet effects).

9
doc/Ce-gamma.indmftpr Normal file
View File

@ -0,0 +1,9 @@
1 ! Nsort
1 ! Mult(Nsort)
3 ! lmax
complex
1 1 1 2 ! l included for each sort
0 0 0 0 ! l included for each sort
0
-.40 0.40 ! Energy window relative to E_f

171
doc/Ce-gamma.py Normal file
View File

@ -0,0 +1,171 @@
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
LDAFilename = 'Ce-gamma'
Beta = 40
Uint = 6.00
JHund = 0.70
Loops = 3 # 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)
HDFfilename = LDAFilename+'.h5'
# Convert DMFT input:
# Can be commented after the first run
Converter = Wien2kConverter(filename=LDAFilename,repacking=True)
Converter.convert_dmft_input()
#check if there are previous runs:
previous_runs = 0
previous_present = False
if mpi.is_master_node():
ar = HDFArchive(HDFfilename,'a')
if 'iterations' in ar:
previous_present = True
previous_runs = ar['iterations']
else:
previous_runs = 0
previous_present = False
del ar
mpi.barrier()
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)
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
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['SigmaF']
del ar
S.Sigma = mpi.bcast(S.Sigma)
SK.load()
# DMFT loop:
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 ])
# Compute the SumK, possibly fixing mu by dichotomy
if SK.Density_Required and (Iteration_Number > 0):
Chemical_potential = SK.find_mu( precision = 0.000001 )
else:
mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=Chemical_potential))
# Density:
S.G <<= SK.extract_Gloc()[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)
# set atomic levels:
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()
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()):
mpi.report("Mixing Sigma and G with factor %s"%Mix)
if ('SigmaF' in ar):
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 = mpi.bcast(S.G)
S.Sigma = mpi.bcast(S.Sigma)
if (mpi.is_master_node()):
ar['SigmaF'] = S.Sigma
ar['GF'] = S.G
# 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 )
# 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
del ar
#Save stuff:
SK.save()
if (mpi.is_master_node()):
print 'DC after solver: ',SK.dc_imp[SK.invshellmap[0]]
# do some analysis:
mpi.report("Orbital densities of impurity Green function:")
dm1 = S.G.density()
for s in dm1:
mpi.report("Block %s: "%s)
for ii in range(len(dm1[s])):
str = ''
for jj in range(len(dm1[s])):
if (dm1[s][ii,jj].real>0):
str += " %.4f"%(dm1[s][ii,jj].real)
else:
str += " %.4f"%(dm1[s][ii,jj].real)
mpi.report(str)
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
# 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')
mpi.report("Trace of Density Matrix: %s"%d)
#correlation energy:
if (mpi.is_master_node()):
ar = HDFArchive(HDFfilename)
itn = ar['iterations']
correnerg = ar['correnerg%s'%itn]
DCenerg = ar['DCenerg%s'%itn]
del ar
correnerg -= DCenerg[0]
f=open(LDAFilename+'.qdmft','a')
f.write("%.16f\n"%correnerg)
f.close()

11
doc/Ce-gamma.struct Normal file
View File

@ -0,0 +1,11 @@
Title
F LATTICE,NONEQUIV.ATOMS: 1225_Fm-3m
MODE OF CALC=RELA unit=bohr
9.750000 9.750000 9.750000 90.000000 90.000000 90.000000
ATOM 1: X=0.00000000 Y=0.00000000 Z=0.00000000
MULT= 1 ISPLIT= 2
Ce NPT= 781 R0=0.00010000 RMT= 2.0000 Z: 58.0
LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000
0.0000000 1.0000000 0.0000000
0.0000000 0.0000000 1.0000000
0 NUMBER OF SYMMETRY OPERATIONS

70
doc/Ce-gamma_DOS.py Normal file
View File

@ -0,0 +1,70 @@
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
# Creates the data directory, cd into it:
#Prepare_Run_Directory(DirectoryName = "Ce-Gamma")
LDAFilename = 'Ce-gamma'
Beta = 40
Uint = 6.00
JHund = 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
useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
ommin=-4.0
ommax=6.0
N_om=2001
broadening = 0.02
HDFfilename = LDAFilename+'.h5'
# Convert DMFT input:
# Can be commented after the first run
Converter = Wien2kConverter(filename=LDAFilename,repacking=True)
Converter.convert_dmft_input()
Converter.convert_par_proj_input()
#check if there are previous runs:
previous_runs = 0
previous_present = False
if mpi.is_master_node():
ar = HDFArchive(HDFfilename,'a')
if 'iterations' in ar:
previous_present = True
previous_runs = ar['iterations']
else:
previous_runs = 0
previous_present = False
del ar
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 = SumkLDATools(hdf_file=LDAFilename+'.h5',use_lda_blocks=False)
if (mpi.is_master_node()):
print 'DC after reading SK: ',SK.dc_imp[SK.invshellmap[0]]
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.Nmoments=10
# 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.Sigma.save('S.Sigma')
SK.put_Sigma(Sigmaimp = [S.Sigma])
SK.dos_partial(broadening=broadening)

7
doc/Ce-gamma_DOS_script.rst Normal file
View File

@ -0,0 +1,7 @@
.. _Ce-gamma_DOS-script:
Ce-gamma_DOS.py
---------------
.. literalinclude:: Ce-gamma_DOS.py
:language: python

7
doc/Ce-gamma_script.rst Normal file
View File

@ -0,0 +1,7 @@
.. _Ce-gamma-script:
Ce-gamma.py
-----------
.. literalinclude:: Ce-gamma.py
:language: python

BIN
doc/Ce_DOS.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 17 KiB

122
doc/LDADMFTmain.rst Normal file
View File

@ -0,0 +1,122 @@
.. index:: LDA+DMFT calculation
.. _LDADMFTmain:
The LDA+DMFT calculation
========================
After having set up the hdf5 arxive, we can now do our LDA+DMFT calculation. It consists of
initialisation steps, and the actual DMFT self consistency loop.
.. index:: initialisation of LDA+DMFT
Initialisation of the calculation
---------------------------------
Before doing the calculation, we have to intialize all the objects that we will need. The first thing is the
:class:`SumkLDA` class. It contains all basic routines that are necessary to perform a summation in k-space
to get the local quantities used in DMFT. It is initialized by::
from pytriqs.applications.dft.sumk_lda import *
SK = SumkLDA(hdf_file = filename)
The only necessary parameter is the filename of the hdf5 archive. In addition, there are some optional parameters:
* `mu`: The chemical potential at initialization. This value is only used, if there is no other value found in the hdf5 arxive. Standard is 0.0
* `h_field`: External magnetic field, standard is 0.0
* `use_lda_blocks`: If true, the structure of the density matrix is analysed at initialisation, and non-zero matrix elements
are identified. The DMFT calculation is then restricted to these matrix elements, yielding a more efficient solution of the
local interaction problem. Also degeneracies in orbital and spin space are recognised, and stored for later use. Standard value is `False`.
* `lda_data`, `symm_corr_data`, `par_proj_data`, `symm_par_data`, `bands_data`: These string variables define the subgroups in the hdf5 arxive,
where the corresponding information is stored. The standard values are consistent with the standard values in :ref:`interfacetowien`.
At initialisation, the necessary data is read from the hdf5 file. If we restart a calculation from a previous one, also the information on
the degenerate shells, the block structure of the density matrix, the chemical potential, and double counting correction are read.
.. index:: Multiband solver
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)
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
up/down structure.
* `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
`U_interact-J_hund`!
* `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.
* `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.
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:
* `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
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 iterations and the self-consistency condition::
n_loops = 5
for iteration_number in range(n_loops) : # start the DMFT loop
SK.put_Sigma(Sigma_imp = [ S.Sigma ]) # Put self energy to the SumK class
chemical_potential = SK.find_mu() # find the chemical potential for the given density
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
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.save() # save everything to the hdf5 arxive
These basic steps are enough to set up the basic DMFT Loop. For a detailed description of the :class:`SumkLDA` 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 more steps after this, because of the double-counting correction. We first
calculate the density of the impurity problem. Then, the routine `set_dc` takes as parameters this density matrix, the
Coulomb interaction, Hund's rule coupling, and the type of double-counting that should be used. Possible values for `use_dc_formula` are:
* `0`: Full-localised limit
* `1`: DC formula as given in K. Held, Adv. Phys. 56, 829 (2007).
* `2`: Around-mean-field
At the end of the calculation, we can save the Greens function and self energy into a file::
from pytriqs.archive import HDFArchive
import pytriqs.utility.mpi as mpi
if mpi.is_master_node():
R = HDFArchive("single_site_bethe.h5",'w')
R["G"] = S.G
R["Sigma"] = S.Sigma
This is it!
These are the essential steps to do a one-shot LDA+DMFT calculation. For full charge-self consistent calculations, there are some more things
to consider, which we will see later on.

BIN
doc/TutorialDmftproj.pdf Normal file
View File

Binary file not shown.

158
doc/advanced.rst Normal file
View File

@ -0,0 +1,158 @@
.. _advanced:
A more advanced example
=======================
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_lda import *
from pytriqs.applications.dft.converters.wien2k_converter import *
from pytriqs.applications.dft.solver_multiband import *
from pytriqs.gf.local import *
from pytriqs.archive import *
Then we define some parameters::
lda_filename='srvo3'
U = 2.7
J = 0.65
beta = 40
loops = 10 # Number of DMFT sc-loops
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 LDA input
use_matrix = False # True: Slater parameters, False: Kanamori parameters U+2J, U, U-J
use_spinflip = False # use the full rotational invariant interaction?
prec_mu = 0.0001
qmc_cycles = 20000
length_cycle = 200
warming_iterations = 2000
Most of these parameters are self-explaining. The first, `lda_filename`, gives the filename of the input files.
The next step, as described in the previous section, is to convert the input files::
Converter = Wien2kConverter(filename=lda_filename, repacking=True)
Converter.convert_dmft_input()
mpi.barrier()
The command ``mpi.barrier()`` ensures that all nodes wait until the conversion of the input is finished on the master
node. After the conversion, we can check in the hdf5 archive, if previous runs are present, or if we have to start
from scratch::
previous_runs = 0
previous_present = False
if mpi.is_master_node():
ar = HDFArchive(lda_filename+'.h5','a')
if 'iterations' in ar:
previous_present = True
previous_runs = ar['iterations']
del ar
previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present)
# if previous runs are present, no need for recalculating the bloc structure:
calc_blocs = use_blocks and (not previous_present)
Now we can use all this information to initialise the :class:`SumkLDA` class::
SK=SumkLDA(hdf_file=lda_filename+'.h5',use_lda_blocks=calc_blocs)
If there was a previous run, we know already about the block structure, and therefore `UseLDABlocs` is set to `False`.
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)
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
If there are previous runs stored in the hdf5 archive, we can now load the self energy
of the last iteration::
if (previous_present):
if (mpi.is_master_node()):
ar = HDFArchive(lda_filename+'.h5','a')
S.Sigma <<= ar['SigmaF']
del ar
S.Sigma = mpi.bcast(S.Sigma)
The last command is the broadcasting of the self energy from the master node to the slave nodes.
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::
for iteration_number in range(1,loops+1) :
SK.symm_deg_gf(S.Sigma,orb=0) # symmetrise Sigma
SK.put_Sigma(Sigma_imp = [ S.Sigma ]) # put Sigma into the SumK class:
chemical_potential = SK.find_mu( precision = prec_mu ) # find the chemical potential
S.G <<= SK.extract_G_loc()[0] # calculation of the local Green function
mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
if ((iteration_number==1)and(previous_present==False)):
# Init the DC term and the real part of Sigma, if no previous run was found:
dm = S.G.density()
SK.set_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
S.Sigma <<= SK.dc_imp[0]['up'][0,0]
# now calculate new G0:
if (mpi.is_master_node()):
# We can do a mixing of Delta in order to stabilize the DMFT iterations:
S.G0 <<= S.Sigma + inverse(S.G)
ar = HDFArchive(lda_filename+'.h5','a')
if ((iteration_number>1) or (previous_present)):
mpi.report("Mixing input Delta with factor %s"%Delta_mix)
Delta = (Delta_mix * S.G0.delta()) + (1.0-Delta_mix) * ar['DeltaF']
S.G0 <<= S.G0 + S.G0.delta() - Delta
ar['DeltaF'] = S.G0.delta()
S.G0 <<= inverse(S.G0)
del ar
S.G0 = mpi.bcast(S.G0)
# Solve the impurity problem:
S.Solve()
# solution done, do the post-processing:
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
# Now mix Sigma and G with factor Mix, if wanted:
if ((iteratio_number>1) or (previous_present)):
if (mpi.is_master_node()):
ar = HDFArchive(lda_filename+'.h5','a')
mpi.report("Mixing Sigma and G with factor %s"%mix)
S.Sigma <<= mix * S.Sigma + (1.0-mix) * ar['SigmaF']
S.G <<= mix * S.G + (1.0-mix) * ar['GF']
del ar
S.G = mpi.bcast(S.G)
S.Sigma = mpi.bcast(S.Sigma)
# Write the final Sigma and G to the hdf5 archive:
if (mpi.is_master_node()):
ar = HDFArchive(lda_filename+'.h5','a')
ar['iterations'] = previous_runs + iteration_number
ar['SigmaF'] = S.Sigma
ar['GF'] = S.G
del ar
# Now set new double counting:
dm = S.G.density()
SK.set_dc( dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
#Save stuff:
SK.save()
This is all we need for the LDA+DMFT calculation. At the end, all results are stored in the hdf5 output file.

117
doc/analysis.rst Normal file
View File

@ -0,0 +1,117 @@
.. _analysis:
Analysing tools
===============
This section explains how to use some tools of the package in order to analyse the data.
.. warning::
The package does NOT provide an explicit method to do an analytic continuation of the
self energies and Green functions from Matsubara frequencies to the real frequancy axis!
There are methods included e.g. in the ALPS package, which can be used for these purposes. But
be careful: All these methods have to be used very carefully!!
The analysing tools can be found in an extension of the :class:`SumkLDA` class, they are
loaded by::
from pytriqs.applications.dft.sumk_lda_tools import *
This import the module ``SumkLDATools``. There are two practical tools, for which you don't
need a self energy on the real axis:
* The density of states of the Wannier orbitals.
* Partial charges according to the Wien2k definition.
Other routines need the self energy on the real frequency axis. If you managed to get them, you can
calculate
* the momentum-integrated spectral function including self-energy effects.
* the momentum-resolved spectral function (i.e. ARPES)
The initialisation of the class is completely equivalent to the initialisation of the :class:`SumkLDA`
class::
SK = SumkLDATools(hdf_file = filename)
By the way, all routines available in :class:`SumkLDA` are also available here.
Routines without real-frequency self energy
-------------------------------------------
For plotting the
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
is an optional parameter, `broadening`, which defines an additional Lorentzian broadening, and is set to `0.01`
by default.
Since we can calculate the partial charges directly from the Matsubara Green's functions, we also don't need a
real frequency self energy for this purpose. The calculation is done by::
ar = HDFArchive(SK.hdf_file)
SK.put_Sigma([ ar['SigmaF'] ])
del ar
dm = SK.partial_charges()
which calculates the partial charges using the input that is stored in the hdf5 file (self energy, double counting,
chemical potential). Here we assumed that the final self energy is stored as `SigmaF` in the archive.
On return, dm is a list, where the list items correspond to the density matrices of all shells
defined in the list ``SK.shells``. This list is constructed by the Wien2k converter routines and stored automatically
in the hdf5 archive. For the detailed structure of `dm`, see the reference manual.
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
keep in mind. The real frequency self energy has to carry the note `ReFreq`::
SigmaReFreq.note = 'ReFreq'
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
counting correction was already read in the initialisation process.
With this self energy, we can do now::
SK.dos_partial()
This produces the momentum-integrated spectral functions (density of states, DOS), also orbitally resolved.
The output is printed into the files
* `DOScorr(sp).dat`: The total DOS. `(sp)` stands for `up`, `down`, or combined `ud`. The latter case
is relevant for calculations including spin-orbit interaction.
* `DOScorr(sp)_proj(i).dat`: The DOS projected to an orbital with index `(i)`. The index `(i)` refers to
the indices given in ``SK.shells``.
* `DOScorr(sp)_proj(i)_(m)_(n).dat`: Sames as above, but printed as orbitally resolved matrix in indices
`(m)` and `(n)`. For `d` orbitals, it gives separately the DOS
for, e.g., :math:`d_{xy}`, :math:`d_{x^2-y^2}`, and so on.
Another quantity of interest is the momentum-resolved spectral function, which can directly be compared to ARPES
experiments. We assume here that we already converted the output of the :program:`dmftproj` program with the
converter routines, see :ref:`interfacetowien`. The spectral function is calculated by::
SK.spaghettis(broadening)
The variable `broadening`1 is an additional Lorentzian broadening that is added to the resulting spectra. The output is
written as the 3-column files ``Akw(sp).dat``, where `(sp)` has the same meaning as above. The output format is
`k`, :math:`\omega`, `value`. Optional parameters are
* `shift`: An additional shift, added as `(ik-1)*shift`, where `ik` is the index of the `k` point. Useful for plotting purposes,
standard value is 0.0.
* `plotrange`: A python list with two entries, first being :math:`\omega_{min}`, the second :math:`\omega_{max}`, setting the plot
range for the output. Standard value is `None`, in this case the momentum range as given in the self energy is plotted.
* `ishell`: If this is not `None` (standard value), but an integer, the spectral function projected to the orbital with index `ishell`
is plotted to the files. Attention: The spectra are not rotated to the local coordinate system as used in the :program:`Wien2k`
program (For experts).

242
doc/conf.py.in Normal file
View File

@ -0,0 +1,242 @@
# -*- coding: utf-8 -*-
#
# TRIQS documentation build configuration file, created by
# sphinx-quickstart on Tue Jul 23 15:29:36 2013.
#
# This file is execfile()d with the current directory set to its containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
import sys, os
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#sys.path.insert(0, os.path.abspath('.'))
# -- General configuration -----------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.pngmath', 'sphinx.ext.mathjax', 'matplotlib.sphinxext.plot_directive']
# Add any paths that contain templates here, relative to this directory.
#templates_path = ['_templates']
# The suffix of source filenames.
source_suffix = '.rst'
# The encoding of source files.
#source_encoding = 'utf-8-sig'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = u'TRIQS interface to Wien2k'
copyright = u'2013, M. Aichhorn, L. Pourovskii, V. Vildosola, C. Martins'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = '1.0'
# The full version, including alpha/beta/rc tags.
release = '1.0.0'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#language = None
# There are two options for replacing |today|: either, you set today to some
# non-false value, then it is used:
#today = ''
# Else, today_fmt is used as the format for a strftime call.
#today_fmt = '%B %d, %Y'
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = []
# The reST default role (used for this markup: `text`) to use for all documents.
#default_role = None
# If true, '()' will be appended to :func: etc. cross-reference text.
#add_function_parentheses = True
# If true, the current module name will be prepended to all description
# unit titles (such as .. function::).
#add_module_names = True
# If true, sectionauthor and moduleauthor directives will be shown in the
# output. They are ignored by default.
#show_authors = False
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
# -- Options for HTML output ---------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
html_theme = 'agogo'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#html_theme_options = {}
# Add any paths that contain custom themes here, relative to this directory.
html_theme_path = ['@TRIQS_THEMES_PATH@']
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".
#html_title = None
# A shorter title for the navigation bar. Default is the same as html_title.
#html_short_title = None
# The name of an image file (relative to this directory) to place at the top
# of the sidebar.
#html_logo = None
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
#html_favicon = None
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
#html_static_path = ['_static']
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
#html_last_updated_fmt = '%b %d, %Y'
# If true, SmartyPants will be used to convert quotes and dashes to
# typographically correct entities.
#html_use_smartypants = True
# Custom sidebar templates, maps document names to template names.
#html_sidebars = {}
# Additional templates that should be rendered to pages, maps page names to
# template names.
#html_additional_pages = {}
# If false, no module index is generated.
#html_domain_indices = True
# If false, no index is generated.
#html_use_index = True
# If true, the index is split into individual pages for each letter.
#html_split_index = False
# If true, links to the reST sources are added to the pages.
#html_show_sourcelink = True
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
#html_show_sphinx = True
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
#html_show_copyright = True
# If true, an OpenSearch description file will be output, and all pages will
# contain a <link> tag referring to it. The value of this option must be the
# base URL from which the finished HTML is served.
#html_use_opensearch = ''
# This is the file name suffix for HTML files (e.g. ".xhtml").
#html_file_suffix = None
# Output file base name for HTML help builder.
htmlhelp_basename = 'TRIQSdoc'
# -- Options for LaTeX output --------------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title, author, documentclass [howto/manual]).
latex_documents = [
('index', 'TRIQS.tex', u'TRIQS Documentation',
u'MF OP', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
# the title page.
#latex_logo = None
# For "manual" documents, if this is true, then toplevel headings are parts,
# not chapters.
#latex_use_parts = False
# If true, show page references after internal links.
#latex_show_pagerefs = False
# If true, show URL addresses after external links.
#latex_show_urls = False
# Documents to append as an appendix to all manuals.
#latex_appendices = []
# If false, no module index is generated.
#latex_domain_indices = True
# -- Options for manual page output --------------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
('index', 'triqs', u'TRIQS Documentation',
[u'MF OP'], 1)
]
# If true, show URL addresses after external links.
#man_show_urls = False
# -- Options for Texinfo output ------------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
('index', 'TRIQS', u'TRIQS Documentation',
u'MF OP', 'TRIQS', 'One line description of project.',
'Miscellaneous'),
]
# Documents to append as an appendix to all manuals.
#texinfo_appendices = []
# If false, no module index is generated.
#texinfo_domain_indices = True
# How to display URL addresses: 'footnote', 'no', or 'inline'.
#texinfo_show_urls = 'footnote'

41
doc/index.rst Normal file
View File

@ -0,0 +1,41 @@
.. index:: Wien2k
.. module:: pytriqs.applications.dft
.. _wien2k:
Wien2TRIQS
========================================================
`Connecting TRIQS to the Wien2k package`
In this section we will learn how we can connect our DMFT calculations with
realistic band structure calculations.
A priori TRIQS can be connected to various realistic band structure codes.
In this release, we provide the Wien2TRIQS extension module which contains an
interface to the `Wien2k package <http://www.wien2k.at>`_.
.. toctree::
:maxdepth: 1
interface
LDADMFTmain
advanced
analysis
selfcons
Ce-HI
In addition to the python-related modules, TRIQS also
provides the Wien2k add-on :program:`dmftproj`. It takes the
information about the wave functions calculated by the `Wien2k package
<http://www.wien2k.at>`_, and constructs projected Wannier functions
that are used as localised orbitals for the DMFT calculation.
The program :program:`dmftproj` is written in the flavor of the
`Wien2k package <http://www.wien2k.at>`_ without python
support. A detailed description of the usage and options of
:program:`dmftproj`
can be found in :download:`this extensive tutorial <TutorialDmftproj.pdf>`. In
addition, it contains also a description of the Wien2k scripts that
are necessary to do the full charge self-consistent calculations.

100
doc/interface.rst Normal file
View File

@ -0,0 +1,100 @@
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
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.
.. index:: Interface to Wien2k
.. _interfacetowien:
The interface to Wien2k
-----------------------
As explained above, this interface produces an hdf5 arxive out of the files that
were written by the band structure package :program:`Wien2k/dmftproj`.
For this purpose we
use the python module :class:`Wien2kConverter`. It is initialised as::
from pytriqs.applications.dft.converters.wien2k_converter import *
Converter = Wien2kConverter(filename = material_of_interest)
The only necessary parameter to this construction is the parameter `filename`.
It has to be the root of the files produces by dmftproj. For example, if you did a
calculation for TiO, the :program:`Wien2k` naming convention is that all files are called
:file:`TiO.*`, so you would give `filename = "TiO"`. The constructor opens
an hdf5 arxive, named :file:`material_of_interest.h5`, where all the data is stored.
There are three optional parameters to the Constructor:
* `lda_subgrp`: We store all data in sub groups of the hdf5 arxive. For the main data
that is needed for the DMFT loop, we use the sub group specified by this optional parameter.
If it is not given, the standard value `SumK_LDA` is used as sub group name.
* `symm_subgrp`: In this sub group we store all the data for applying the symmetry
operations in the DMFT loop. Standard value is `SymmCorr`.
* `repacking`: If true, and the hdf5 file already exists, the system command :program:`h5repack`
is invoked. This command ensures a minimal file size of the hdf5
file. Standard value is `False`. If you want to use this, be sure
that :program:`h5repack` is in your path variable!
After initialising the interface module, we can now convert the input text files into the
hdf5 arxive by::
Converter.convert_dmft_input()
This reads all the data, and stores it in the sub group `lda_subgrp`, as discussed above.
In this step, the files :file:`material_of_interest.ctqmcout` and :file:`material_of_interest.symqmc`
have to be present in the working directory.
After this step, all the necessary information for the DMFT loop is stored in the hdf5 arxive, where
the string variable `Converter.hdf_file` gives the file name of the arxive.
You can now proceed with :ref:`LDADMFTmain`.
Data for post-processing
------------------------
In order to calculate some properties using the DMFT self energy, several other routines are
used in order to convert the necessary input from :program:`Wien2k/dmftproj`. For instance, for
calculating the partial density of states or partial charges consistent with the definition
of :program:`Wien2k`, you have to use::
Converter.convert_parproj_input()
This reads the files :file:`material_of_interest.parproj` and :file:`material_of_interest.sympar`.
Again, there are two optional parameters
* `par_proj_subgrp`: The sub group, where the data for the partial projectors is stored. Standard
is `SumK_LDA_ParProj`.
* `symm_par_subgrp`: Sub group for the symmetry operations, standard value is `SymmPar`.
Another routine of the class allows to read the input for plotting the momentum-resolved
spectral function. It is done by::
Converter.convert_bands_input()
The optional parameter, which tells the routine where to store the data is here `bands_subgrp`,
and its standard value is `SumK_LDA_Bands`.
After having converted this input, you can further proceed with the :ref:`analysis`.
MPI issues
----------
The interface package is written such that all the operations are done only on the master node.
The broadcasting to the nodes has to be done by hand. The :class:`SumkLDA`, described in the
following section, takes care of this automatically.
Interfaces to other packages
----------------------------
Because of the modular structure, it is straight forward to extend the TRIQS package
in order to work with other band-structure codes. The only necessary requirement is that
the interface module produces an hdf5 arxive, that stores all the data in the specified
form. For the details of what data is stored in detail, see the reference manual.

89
doc/selfcons.rst Normal file
View File

@ -0,0 +1,89 @@
.. index:: full charge self consistency
Full charge self consistency
============================
.. warning::
Before using this tool, you should be familiar with the band-structure package :program:`Wien2k`, since
the calculation is controlled by the :program:`Wien2k` scripts! See also the :download:`dmftproj tutorial<TutorialDmftproj.pdf>`.
In order to do charge self-consistent calculations, we have to tell the band structure program about the
changes in the charge density due to correlation effects. In the following, we discuss how to use the
:program:`TRIQS` tools in combination with the :program:`Wien2k` program, although an extension to other
codes is also possible.
We can use the DMFT script as introduced in sections :ref:`LDADMFTmain` and :ref:`advanced`, with a few simple
modifications. First, in order to be compatible with the :program:`Wien2k` standards, the DMFT script has to be
named ``case.py``, where `case` is the name of the :program:`Wien2k` calculation, see the section
:ref:`interfacetowien` for details. Then we set the variable
`lda_filename` dynamically::
import os
lda_filename = os.getcwd().rpartition('/')[2]
This sets the `lda_filename` to the name of the current directory. The reminder of the scripts is completely the
same as in one-shot calculations. Only at the very end we have to calculate the modified charge density,
and store it in a format such that :program:`Wien2k` can read it. Therefore, after the DMFT loop that we saw in the
previous section, we symmetrise the self energy, and recalculate the impurity Green function::
SK.symm_deg_gf(S.Sigma,orb=0)
S.G <<= inverse(S.G0) - S.Sigma
S.G.invert()
These steps are not necessary, but can help to reduce fluctuation of the total energy.
Now we calculate the modified charge density::
# find exact chemical potential
SK.put_Sigma(Sigma_imp = [ S.Sigma ])
chemical_potential = SK.find_mu( precision = 0.000001 )
dN,d = SK.calc_density_correction(filename = lda_filename+'.qdmft')
SK.save()
First we find the chemical potential with high precision, and after that the routine
``SK.calc_density_correction(filename)`` calculates the density matrix including correlation effects. The result
is stored in the file `Filename`, which is later read by the :program:`Wien2k` program. The last statement saves
the chemical potential into the hdf5 archive.
We need also the correlation energy, which we evaluate by the Migdal formula::
correnerg = 0.5 * (S.G * S.Sigma).total_density()
From this value, we have to substract the double counting energy::
correnerg -= SK.dc_energ[0]
and save this value into the file::
if (mpi.is_master_node()):
f=open(lda_filename+'.qdmft','a')
f.write("%.16f\n"%correnerg)
f.close()
The above steps are valid for a calculation with only one correlated atom in the unit cell, the most likely case
where you will apply this method. That is the reason why we give the index `0` in the list `SK.dc_energ`.
If you have more than one correlated atom in the unit cell, but all of them
are equivalent atoms, you have to multiply the `correnerg` by their multiplicity, before writing it to the file.
The multiplicity is easily found in the main input file of the :program:`Wien2k` package, i.e. `case.struct`. In case of
non-equivalent atoms, the correlation energy has to be calculated for all of them separately (FOR EXPERTS ONLY).
As mentioned above, the calculation is controlled by the :program:`Wien2k` scripts and not by :program:`python`
routines. Therefore, you start your calculation for instance by::
me@home $ run -qdmft -i 10
The flag `-qdmft` tells the script, that the density matrix including correlation effects is read from the `case.qdmft`
file, and 10 self-consitency iterations are done. If you run the code on a parallel machine, you can specify the number of
nodes that are used::
me@home $ run -qdmft -np 64 -i 10
with the `-np` flag. In that case, you have to give the proper `MPI` execution statement, e.g. `mpiexec`, in the `run_lapw` script,
see the corresponding :program:`Wien2k` documentation. In many cases it is advisable to start from a converged one-shot
calculation.
For practical purposes, you keep the number of DMFT loops within one DFT cycle low, or even to `loops=1`. If you encouter
unstable convergence, you have to adjust the parameters such as
`loops`, `mix`, or `Delta_mix` to improve the convergence.
In the next section, :ref:`LDADMFTtutorial`, we will see in a detailed
example, how such a self consistent calculation is performed.