From 00a775a93d89200a5dcbdae17bd4a4b85ed4ed4c Mon Sep 17 00:00:00 2001 From: Manuel Zingl Date: Thu, 20 Aug 2015 15:46:14 +0200 Subject: [PATCH] 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!? --- doc/guide/analysis.rst | 89 +++++++++++------ doc/guide/transport.rst | 2 + python/build_sigma_from_txt.py | 178 ++++++++++++++++----------------- 3 files changed, 147 insertions(+), 122 deletions(-) diff --git a/doc/guide/analysis.rst b/doc/guide/analysis.rst index 24c46f97..67c0f2bf 100644 --- a/doc/guide/analysis.rst +++ b/doc/guide/analysis.rst @@ -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: - * :meth:`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:`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. 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:`spaghettis` for the momentum-resolved spectral function (i.e. ARPES) + * :meth:`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) .. warning:: 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. 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 -------------- @@ -40,18 +36,40 @@ class:: 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, -your self energy is already stored as a real frequency :class:`BlockGf` object +If required, we have to load and initialise the real frequency self energy. Most conveniently, +you have your self energy already stored as a real frequency :class:`BlockGf` object in a hdf5 file:: ar = HDFArchive('case.h5', 'a') SigmaReFreq = ar['dmft_output']['Sigma_w'] - SK.put_Sigma(Sigma_imp = [SigmaReFreq]) -Additionally, the chemical potential and the double counting -correction from the DMFT calculation are set, and the archive is closed again:: +You may also have your self energy stored in text files. For this case we provide the function +: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_dc(dc_imp,dc_energ) 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) 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 -density of states after the projection set `with_Sigma` and `with_dc` to `False`. +`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`. 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 --------------- @@ -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:: 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 `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) -------------------------------------------------------------- -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) - -This produces both the momentum-integrated (total density of states or DOS) and orbitally-resolved (partial/projected DOS) spectral functions. -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. +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 +`DOS_parproj_*` instead. 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 -experiments. We assume here that we already converted the output of the :program:`dmftproj` program with the -converter routines (see :ref:`conversion`). The spectral function is calculated by:: +experiments. First we have to execute `lapw1`, `lapw2 -almd` and :program:`dmftproj` with the `-band` +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 diff --git a/doc/guide/transport.rst b/doc/guide/transport.rst index d1f3025b..4d32b0af 100644 --- a/doc/guide/transport.rst +++ b/doc/guide/transport.rst @@ -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) +The converter :meth:`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:: ar = HDFArchive('case.h5', 'a') diff --git a/python/build_sigma_from_txt.py b/python/build_sigma_from_txt.py index 1f8ab975..387d810b 100644 --- a/python/build_sigma_from_txt.py +++ b/python/build_sigma_from_txt.py @@ -1,99 +1,99 @@ - def read_fortran_file (filename): - """ Returns a generator that yields all numbers in the Fortran file as float, one by one""" - import os.path - if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename - for line in open(filename,'r') : - for x in line.replace('D','E').split() : - yield string.atof(x) +def read_fortran_file (filename): + """ Returns a generator that yields all numbers in the Fortran file as float, one by one""" + import os.path + if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename + for line in open(filename,'r') : + for x in line.replace('D','E').split() : + 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): - """Uses Data from files to construct Sigma (or GF) on the real axis.""" +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.""" - 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: - bl = self.gf_struct_solver[orb].items()[0][0] # block name - ol = self.gf_struct_solver[orb].items()[0][1] # list of orbital indices - if (len(ol)==1): # if blocks are of size one - Fname = filename+'_'+bl+'.dat' - else: - Fname = filename+'_'+bl+'/'+str(ol[0])+'_'+str(ol[0])+'.dat' + # first get the mesh out of any one of the files: + bl = self.gf_struct_solver[orb].items()[0][0] # block name + ol = self.gf_struct_solver[orb].items()[0][1] # list of orbital indices + if (len(ol)==1): # if blocks are of size one + Fname = filename+'_'+bl+'.dat' + else: + Fname = filename+'_'+bl+'/'+str(ol[0])+'_'+str(ol[0])+'.dat' - R = read_fortran_file(Fname) - mesh = numpy.zeros([n_om],numpy.float_) - 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) + R = read_fortran_file(Fname) + mesh = numpy.zeros([n_om],numpy.float_) + try: 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()] - 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) + # pass SigmaME to other nodes + SigmaME = mpi.bcast(SigmaME) + mpi.barrier() - #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() + SigmaME.note='ReFreq' - - 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 + return SigmaME