mirror of
https://github.com/triqs/dft_tools
synced 2024-08-29 23:33:43 +02:00
analysis.rst done. Minor change in transport.rst
I also described how one can read a self energy form a data file. However, this needs to be tested and also included in the reference manual. Maybe the function should move back into sumk_dft_tools!?
This commit is contained in:
parent
6a50b57c59
commit
00a775a93d
@ -7,13 +7,13 @@ This section explains how to use some tools of the package in order to analyse t
|
|||||||
|
|
||||||
There are two practical tools for which a self energy on the real axis is not needed, namely:
|
There are two practical tools for which a self energy on the real axis is not needed, namely:
|
||||||
|
|
||||||
* :meth:`dos_wannier_basis` for the density of states of the Wannier orbitals and
|
* :meth:`dos_wannier_basis <pytriqs.applications.dft.sumk_dft_tools.dos_wannier_basis>` for the density of states of the Wannier orbitals and
|
||||||
* :meth:`partial_charges` for the partial charges according to the Wien2k definition.
|
* :meth:`partial_charges <pytriqs.applications.dft.sumk_dft_tools.partial_charges>` for the partial charges according to the Wien2k definition.
|
||||||
|
|
||||||
However, a real frequency self energy has to be provided by the user to use the methods:
|
However, a real frequency self energy has to be provided by the user to use the methods:
|
||||||
|
|
||||||
* :meth:`dos_parproj_basis` for the momentum-integrated spectral function including self energy effects and
|
* :meth:`dos_parproj_basis <pytriqs.applications.dft.sumk_dft_tools.dos_parproj_basis>` for the momentum-integrated spectral function including self energy effects and
|
||||||
* :meth:`spaghettis` for the momentum-resolved spectral function (i.e. ARPES)
|
* :meth:`spaghettis <pytriqs.applications.dft.sumk_dft_tools.spaghettis>` for the momentum-resolved spectral function (i.e. ARPES)
|
||||||
|
|
||||||
.. warning::
|
.. warning::
|
||||||
This package does NOT provide an explicit method to do an **analytic continuation** of the
|
This package does NOT provide an explicit method to do an **analytic continuation** of the
|
||||||
@ -21,10 +21,6 @@ However, a real frequency self energy has to be provided by the user to use the
|
|||||||
There are methods included e.g. in the :program:`ALPS` package, which can be used for these purposes.
|
There are methods included e.g. in the :program:`ALPS` package, which can be used for these purposes.
|
||||||
Keep in mind that all these methods have to be used very carefully!
|
Keep in mind that all these methods have to be used very carefully!
|
||||||
|
|
||||||
.. note::
|
|
||||||
Add the doc for loading the self energy from a data file. We have to provide this option, because
|
|
||||||
in general the user won't has it stored in h5 file!!
|
|
||||||
|
|
||||||
Initialisation
|
Initialisation
|
||||||
--------------
|
--------------
|
||||||
|
|
||||||
@ -40,18 +36,40 @@ class::
|
|||||||
|
|
||||||
Note that all routines available in :class:`SumkDFT` are also available here.
|
Note that all routines available in :class:`SumkDFT` are also available here.
|
||||||
|
|
||||||
If required, a self energy is load and initialise in the next step. Most conveniently,
|
If required, we have to load and initialise the real frequency self energy. Most conveniently,
|
||||||
your self energy is already stored as a real frequency :class:`BlockGf` object
|
you have your self energy already stored as a real frequency :class:`BlockGf` object
|
||||||
in a hdf5 file::
|
in a hdf5 file::
|
||||||
|
|
||||||
ar = HDFArchive('case.h5', 'a')
|
ar = HDFArchive('case.h5', 'a')
|
||||||
SigmaReFreq = ar['dmft_output']['Sigma_w']
|
SigmaReFreq = ar['dmft_output']['Sigma_w']
|
||||||
SK.put_Sigma(Sigma_imp = [SigmaReFreq])
|
|
||||||
|
|
||||||
Additionally, the chemical potential and the double counting
|
You may also have your self energy stored in text files. For this case we provide the function
|
||||||
correction from the DMFT calculation are set, and the archive is closed again::
|
:meth:`constr_Sigma_real_axis`, which loads the data and puts it into a real frequency :class:`BlockGf` object::
|
||||||
|
|
||||||
|
from pytriqs.applications.dft.build_sigma_from_txt import *
|
||||||
|
SigmaReFreq = constr_Sigma_real_axis(SK, filename, hdf=False, hdf_dataset='SigmaReFreq',n_om=0, orb=0)
|
||||||
|
|
||||||
|
where:
|
||||||
|
|
||||||
|
* `filename`: the name of the hdf5 archive file or the `fname` pattern in text files names as described above,
|
||||||
|
* `hdf`: if `True`, the real axis self energy will be read from the hdf5 file, otherwise from the text files,
|
||||||
|
* `hdf_dataset`: the name of dataset where the self energy is stored in the hdf5 file,
|
||||||
|
* `orb`: index of an inequivalent shell,
|
||||||
|
* `n_om`: the number of points in the real-axis mesh (used only if `hdf=False`).
|
||||||
|
|
||||||
|
|
||||||
|
It is important that some rules concerning the structure of the data is followed:
|
||||||
|
* Each data file should contain the three columns: real frequency, real part and imaginary part of the self-energy in exactly this order.
|
||||||
|
* 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, where `(i)` and `(j)` are the zero based orbital indices, in the `filename_(block)` directory.
|
||||||
|
|
||||||
chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
|
Finally, we put the self energy into the `SK` object::
|
||||||
|
|
||||||
|
SK.put_Sigma(Sigma_imp = [SigmaReFreq])
|
||||||
|
|
||||||
|
and additionally set the chemical potential and the double counting correction from the DMFT calculation::
|
||||||
|
|
||||||
|
chemical_potential, dc_imp, dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
|
||||||
SK.set_mu(chemical_potential)
|
SK.set_mu(chemical_potential)
|
||||||
SK.set_dc(dc_imp,dc_energ)
|
SK.set_dc(dc_imp,dc_energ)
|
||||||
del ar
|
del ar
|
||||||
@ -65,8 +83,18 @@ For plotting the density of states of the Wannier orbitals, you type::
|
|||||||
SK.dos_wannier_basis(broadening=0.03, mesh=[om_min, om_max, n_om], with_Sigma=False, with_dc=False, save_to_file=True)
|
SK.dos_wannier_basis(broadening=0.03, mesh=[om_min, om_max, n_om], with_Sigma=False, with_dc=False, save_to_file=True)
|
||||||
|
|
||||||
which produces plots between the real frequencies `om_min` and `om_max`, using a mesh of `n_om` points. The parameter
|
which produces plots between the real frequencies `om_min` and `om_max`, using a mesh of `n_om` points. The parameter
|
||||||
`broadening` defines an additional Lorentzian broadening, and has the default value of `0.01` eV. To check the Wannier
|
`broadening` defines an additional Lorentzian broadening, and has the default value of `0.01 eV`. To check the Wannier
|
||||||
density of states after the projection set `with_Sigma` and `with_dc` to `False`.
|
density of states after the projection set `with_Sigma` and `with_dc` to `False`. If `save_to_file` is set to `True`
|
||||||
|
the output is printed into the files
|
||||||
|
|
||||||
|
* `DOS_wannier_(sp).dat`: The total DOS, where `(sp)` stands for `up`, `down`, or combined `ud`. The latter case
|
||||||
|
is relevant for calculations including spin-orbit interaction.
|
||||||
|
* `DOS_wannier_(sp)_proj(i).dat`: The DOS projected to an orbital with index `(i)`. The index `(i)` refers to
|
||||||
|
the indices given in ``SK.shells``.
|
||||||
|
* `DOS_wannier_(sp)_proj(i)_(m)_(n).dat`: As above, but printed as orbitally-resolved matrix in indices
|
||||||
|
`(m)` and `(n)`. For `d` orbitals, it gives the DOS seperately for, e.g., :math:`d_{xy}`, :math:`d_{x^2-y^2}`, and so on,
|
||||||
|
|
||||||
|
otherwise, the ouptut is returend by the function for further use in python.
|
||||||
|
|
||||||
Partial charges
|
Partial charges
|
||||||
---------------
|
---------------
|
||||||
@ -75,7 +103,7 @@ Since we can calculate the partial charges directly from the Matsubara Green's f
|
|||||||
real frequency self energy for this purpose. The calculation is done by::
|
real frequency self energy for this purpose. The calculation is done by::
|
||||||
|
|
||||||
SK.put_Sigma(Sigma_imp = SigmaImFreq)
|
SK.put_Sigma(Sigma_imp = SigmaImFreq)
|
||||||
dm = SK.partial_charges(beta=40.0 with_Sigma=True, with_dc=True)
|
dm = SK.partial_charges(beta=40.0, with_Sigma=True, with_dc=True)
|
||||||
|
|
||||||
which calculates the partial charges using the self energy, double counting, and chemical potential as set in the
|
which calculates the partial charges using the self energy, double counting, and chemical potential as set in the
|
||||||
`SK` object. On return, `dm` is a list, where the list items correspond to the density matrices of all shells
|
`SK` object. On return, `dm` is a list, where the list items correspond to the density matrices of all shells
|
||||||
@ -85,29 +113,24 @@ in the hdf5 archive. For the detailed structure of `dm`, see the reference manua
|
|||||||
Correlated spectral function (with real frequency self energy)
|
Correlated spectral function (with real frequency self energy)
|
||||||
--------------------------------------------------------------
|
--------------------------------------------------------------
|
||||||
|
|
||||||
With this self energy, we can now execute::
|
To produce both the momentum-integrated (total density of states or DOS) and orbitally-resolved (partial/projected DOS) spectral functions
|
||||||
|
we can execute::
|
||||||
|
|
||||||
|
SK.dos_parproj_basis(broadening=0.0, with_Sigma=True, with_dc=True, save_to_file=True)
|
||||||
|
|
||||||
SK.dos_parproj_basis(broadening=broadening)
|
The variable `broadening` is an additional Lorentzian broadening (default: `0.01 eV`) applied to the resulting spectra.
|
||||||
|
The output is written in the same way as described above for Wannier density of states, but with file names
|
||||||
This produces both the momentum-integrated (total density of states or DOS) and orbitally-resolved (partial/projected DOS) spectral functions.
|
`DOS_parproj_*` instead.
|
||||||
The variable `broadening` is an additional Lorentzian broadening applied to the resulting spectra.
|
|
||||||
The output is printed into the files
|
|
||||||
|
|
||||||
* `DOScorr(sp).dat`: The total DOS, where `(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`: As above, but printed as orbitally-resolved matrix in indices
|
|
||||||
`(m)` and `(n)`. For `d` orbitals, it gives the DOS seperately for, e.g., :math:`d_{xy}`, :math:`d_{x^2-y^2}`, and so on.
|
|
||||||
|
|
||||||
Momentum resolved spectral function (with real frequency self energy)
|
Momentum resolved spectral function (with real frequency self energy)
|
||||||
---------------------------------------------------------------------
|
---------------------------------------------------------------------
|
||||||
|
|
||||||
Another quantity of interest is the momentum-resolved spectral function, which can directly be compared to ARPES
|
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
|
experiments. First we have to execute `lapw1`, `lapw2 -almd` and :program:`dmftproj` with the `-band`
|
||||||
converter routines (see :ref:`conversion`). The spectral function is calculated by::
|
option and use the :meth:`convert_bands_input()` routine to convert the required files. For a detailed description
|
||||||
|
see :ref:`conversion`. The spectral function is then calculated by::
|
||||||
|
|
||||||
SK.spaghettis(broadening)
|
SK.spaghettis(broadening=0.01,plot_shift=0.0,plot_range=None,ishell=None,save_to_file='Akw_')
|
||||||
|
|
||||||
Optional parameters are
|
Optional parameters are
|
||||||
|
|
||||||
|
@ -81,6 +81,8 @@ First we have to read the Wien2k files and store the relevant information in the
|
|||||||
|
|
||||||
SK = SumkDFTTools(hdf_file='case.h5', use_dft_blocks=True)
|
SK = SumkDFTTools(hdf_file='case.h5', use_dft_blocks=True)
|
||||||
|
|
||||||
|
The converter :meth:`convert_transport_input <pytriqs.applications.dft.converters.wien2k_converter.Wien2kConverter.convert_transport_input>`
|
||||||
|
reads the required data of the Wien2k output and stores it in the `dft_transp_input` subgroup of your hdf file.
|
||||||
Additionally we need to read and set the self energy, the chemical potential and the double counting::
|
Additionally we need to read and set the self energy, the chemical potential and the double counting::
|
||||||
|
|
||||||
ar = HDFArchive('case.h5', 'a')
|
ar = HDFArchive('case.h5', 'a')
|
||||||
|
@ -1,99 +1,99 @@
|
|||||||
def read_fortran_file (filename):
|
def read_fortran_file (filename):
|
||||||
""" Returns a generator that yields all numbers in the Fortran file as float, one by one"""
|
""" Returns a generator that yields all numbers in the Fortran file as float, one by one"""
|
||||||
import os.path
|
import os.path
|
||||||
if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename
|
if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename
|
||||||
for line in open(filename,'r') :
|
for line in open(filename,'r') :
|
||||||
for x in line.replace('D','E').split() :
|
for x in line.replace('D','E').split() :
|
||||||
yield string.atof(x)
|
yield string.atof(x)
|
||||||
|
|
||||||
def constr_Sigma_real_axis(self, filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0,orb=0, tol_mesh=1e-6):
|
def constr_Sigma_real_axis(self, filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0,orb=0, tol_mesh=1e-6):
|
||||||
"""Uses Data from files to construct Sigma (or GF) on the real axis."""
|
"""Uses Data from files to construct Sigma (or GF) on the real axis."""
|
||||||
|
|
||||||
if not hdf: # then read sigma from text files
|
if not hdf: # then read sigma from text files
|
||||||
|
|
||||||
# first get the mesh out of any one of the files:
|
# first get the mesh out of any one of the files:
|
||||||
bl = self.gf_struct_solver[orb].items()[0][0] # block name
|
bl = self.gf_struct_solver[orb].items()[0][0] # block name
|
||||||
ol = self.gf_struct_solver[orb].items()[0][1] # list of orbital indices
|
ol = self.gf_struct_solver[orb].items()[0][1] # list of orbital indices
|
||||||
if (len(ol)==1): # if blocks are of size one
|
if (len(ol)==1): # if blocks are of size one
|
||||||
Fname = filename+'_'+bl+'.dat'
|
Fname = filename+'_'+bl+'.dat'
|
||||||
else:
|
else:
|
||||||
Fname = filename+'_'+bl+'/'+str(ol[0])+'_'+str(ol[0])+'.dat'
|
Fname = filename+'_'+bl+'/'+str(ol[0])+'_'+str(ol[0])+'.dat'
|
||||||
|
|
||||||
R = read_fortran_file(Fname)
|
R = read_fortran_file(Fname)
|
||||||
mesh = numpy.zeros([n_om],numpy.float_)
|
mesh = numpy.zeros([n_om],numpy.float_)
|
||||||
try:
|
try:
|
||||||
for i in xrange(n_om):
|
|
||||||
mesh[i] = R.next()
|
|
||||||
sk = R.next()
|
|
||||||
sk = R.next()
|
|
||||||
|
|
||||||
except StopIteration : # a more explicit error if the file is corrupted.
|
|
||||||
raise "SumkDFT.read_Sigma_ME : reading mesh failed!"
|
|
||||||
R.close()
|
|
||||||
|
|
||||||
# check whether the mesh is uniform
|
|
||||||
bin = (mesh[n_om-1]-mesh[0])/(n_om-1)
|
|
||||||
for i in xrange(n_om):
|
for i in xrange(n_om):
|
||||||
assert abs(i*bin+mesh[0]-mesh[i]) < tol_mesh, 'constr_Sigma_ME: real-axis mesh is non-uniform!'
|
mesh[i] = R.next()
|
||||||
|
sk = R.next()
|
||||||
|
sk = R.next()
|
||||||
|
|
||||||
# construct Sigma
|
except StopIteration : # a more explicit error if the file is corrupted.
|
||||||
|
raise "SumkDFT.read_Sigma_ME : reading mesh failed!"
|
||||||
|
R.close()
|
||||||
|
|
||||||
|
# check whether the mesh is uniform
|
||||||
|
bin = (mesh[n_om-1]-mesh[0])/(n_om-1)
|
||||||
|
for i in xrange(n_om):
|
||||||
|
assert abs(i*bin+mesh[0]-mesh[i]) < tol_mesh, 'constr_Sigma_ME: real-axis mesh is non-uniform!'
|
||||||
|
|
||||||
|
# construct Sigma
|
||||||
|
a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
|
||||||
|
glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[n_om-1]),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
|
||||||
|
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
|
||||||
|
|
||||||
|
#read Sigma
|
||||||
|
for i,g in SigmaME:
|
||||||
|
mesh=[w for w in g.mesh]
|
||||||
|
for iL in g.indices:
|
||||||
|
for iR in g.indices:
|
||||||
|
if (len(g.indices) == 1):
|
||||||
|
Fname = filename+'_%s'%(i)+'.dat'
|
||||||
|
else:
|
||||||
|
Fname = 'SigmaME_'+'%s'%(i)+'_%s'%(iL)+'_%s'%(iR)+'.dat'
|
||||||
|
R = read_fortran_file(Fname)
|
||||||
|
try:
|
||||||
|
for iom in xrange(n_om):
|
||||||
|
sk = R.next()
|
||||||
|
rsig = R.next()
|
||||||
|
isig = R.next()
|
||||||
|
g.data[iom,iL,iR]=rsig+1j*isig
|
||||||
|
except StopIteration : # a more explicit error if the file is corrupted.
|
||||||
|
raise "SumkDFT.read_Sigma_ME : reading Sigma from file failed!"
|
||||||
|
R.close()
|
||||||
|
|
||||||
|
|
||||||
|
else: # read sigma from hdf
|
||||||
|
|
||||||
|
omega_min=0.0
|
||||||
|
omega_max=0.0
|
||||||
|
n_om=0
|
||||||
|
if (mpi.is_master_node()):
|
||||||
|
ar = HDFArchive(filename)
|
||||||
|
SigmaME = ar[hdf_dataset]
|
||||||
|
del ar
|
||||||
|
|
||||||
|
#OTHER SOLUTION FIXME
|
||||||
|
#else:
|
||||||
|
# SigmaME=0
|
||||||
|
#SigmaME = mpi.broadcast..
|
||||||
|
|
||||||
|
# we need some parameters to construct Sigma on other nodes
|
||||||
|
omega_min=SigmaME.mesh.omega_min
|
||||||
|
omega_max=SigmaME.mesh.omega_max
|
||||||
|
n_om=len(SigmaME.mesh)
|
||||||
|
omega_min=mpi.bcast(omega_min)
|
||||||
|
omega_max=mpi.bcast(omega_max)
|
||||||
|
n_om=mpi.bcast(n_om)
|
||||||
|
mpi.barrier()
|
||||||
|
# construct Sigma on other nodes
|
||||||
|
if (not mpi.is_master_node()):
|
||||||
a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
|
a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
|
||||||
glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[n_om-1]),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
|
glist = lambda : [ GfReFreq(indices = al, window=(omega_min,omega_max),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
|
||||||
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
|
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
|
||||||
|
# pass SigmaME to other nodes
|
||||||
|
SigmaME = mpi.bcast(SigmaME)
|
||||||
|
mpi.barrier()
|
||||||
|
|
||||||
#read Sigma
|
SigmaME.note='ReFreq'
|
||||||
for i,g in SigmaME:
|
|
||||||
mesh=[w for w in g.mesh]
|
|
||||||
for iL in g.indices:
|
|
||||||
for iR in g.indices:
|
|
||||||
if (len(g.indices) == 1):
|
|
||||||
Fname = filename+'_%s'%(i)+'.dat'
|
|
||||||
else:
|
|
||||||
Fname = 'SigmaME_'+'%s'%(i)+'_%s'%(iL)+'_%s'%(iR)+'.dat'
|
|
||||||
R = read_fortran_file(Fname)
|
|
||||||
try:
|
|
||||||
for iom in xrange(n_om):
|
|
||||||
sk = R.next()
|
|
||||||
rsig = R.next()
|
|
||||||
isig = R.next()
|
|
||||||
g.data[iom,iL,iR]=rsig+1j*isig
|
|
||||||
except StopIteration : # a more explicit error if the file is corrupted.
|
|
||||||
raise "SumkDFT.read_Sigma_ME : reading Sigma from file failed!"
|
|
||||||
R.close()
|
|
||||||
|
|
||||||
|
return SigmaME
|
||||||
else: # read sigma from hdf
|
|
||||||
|
|
||||||
omega_min=0.0
|
|
||||||
omega_max=0.0
|
|
||||||
n_om=0
|
|
||||||
if (mpi.is_master_node()):
|
|
||||||
ar = HDFArchive(filename)
|
|
||||||
SigmaME = ar[hdf_dataset]
|
|
||||||
del ar
|
|
||||||
|
|
||||||
#OTHER SOLUTION FIXME
|
|
||||||
#else:
|
|
||||||
# SigmaME=0
|
|
||||||
#SigmaME = mpi.broadcast..
|
|
||||||
|
|
||||||
# we need some parameters to construct Sigma on other nodes
|
|
||||||
omega_min=SigmaME.mesh.omega_min
|
|
||||||
omega_max=SigmaME.mesh.omega_max
|
|
||||||
n_om=len(SigmaME.mesh)
|
|
||||||
omega_min=mpi.bcast(omega_min)
|
|
||||||
omega_max=mpi.bcast(omega_max)
|
|
||||||
n_om=mpi.bcast(n_om)
|
|
||||||
mpi.barrier()
|
|
||||||
# construct Sigma on other nodes
|
|
||||||
if (not mpi.is_master_node()):
|
|
||||||
a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
|
|
||||||
glist = lambda : [ GfReFreq(indices = al, window=(omega_min,omega_max),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
|
|
||||||
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
|
|
||||||
# pass SigmaME to other nodes
|
|
||||||
SigmaME = mpi.bcast(SigmaME)
|
|
||||||
mpi.barrier()
|
|
||||||
|
|
||||||
SigmaME.note='ReFreq'
|
|
||||||
|
|
||||||
return SigmaME
|
|
||||||
|
Loading…
Reference in New Issue
Block a user