2016-05-09 10:19:56 +02:00
##########################################################################
2013-07-23 19:49:42 +02:00
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
2016-05-09 10:19:56 +02:00
##########################################################################
2022-03-01 23:51:39 +01:00
"""
Extension to the SumkDFT class with some analyiss tools
"""
2015-07-28 17:21:13 +02:00
import sys
2013-07-23 19:49:42 +02:00
from types import *
import numpy
2020-05-27 17:30:24 +02:00
from triqs . gf import *
import triqs . utility . mpi as mpi
2020-04-08 21:35:59 +02:00
from . symmetry import *
from . sumk_dft import SumkDFT
2015-07-28 17:21:13 +02:00
from scipy . integrate import *
from scipy . interpolate import *
2013-07-23 19:49:42 +02:00
2018-02-13 21:42:09 +01:00
if not hasattr ( numpy , ' full ' ) :
# polyfill full for older numpy:
numpy . full = lambda a , f : numpy . zeros ( a ) + f
2016-05-09 10:19:56 +02:00
2014-11-18 11:30:26 +01:00
class SumkDFTTools ( SumkDFT ) :
2015-03-12 00:01:12 +01:00
"""
Extends the SumkDFT class with some tools for analysing the data .
"""
2013-07-23 19:49:42 +02:00
2021-08-27 23:46:17 +02:00
def __init__ ( self , hdf_file , h_field = 0.0 , mesh = None , beta = 40 , n_iw = 1025 , use_dft_blocks = False , dft_data = ' dft_input ' , symmcorr_data = ' dft_symmcorr_input ' ,
2016-05-09 10:19:56 +02:00
parproj_data = ' dft_parproj_input ' , symmpar_data = ' dft_symmpar_input ' , bands_data = ' dft_bands_input ' ,
2023-04-15 00:43:23 +02:00
transp_data = ' dft_transp_input ' , misc_data = ' dft_misc_input ' , cont_data = ' dft_contours_input ' ) :
2015-03-12 00:01:12 +01:00
"""
Initialisation of the class . Parameters are exactly as for SumKDFT .
"""
2014-09-22 19:24:33 +02:00
2021-08-27 23:31:07 +02:00
SumkDFT . __init__ ( self , hdf_file = hdf_file , h_field = h_field , mesh = mesh , beta = beta , n_iw = n_iw ,
use_dft_blocks = use_dft_blocks , dft_data = dft_data , symmcorr_data = symmcorr_data ,
parproj_data = parproj_data , symmpar_data = symmpar_data , bands_data = bands_data ,
2023-04-15 00:43:23 +02:00
transp_data = transp_data , misc_data = misc_data , cont_data = cont_data )
2013-07-23 19:49:42 +02:00
2023-04-15 00:43:23 +02:00
def density_of_states ( self , mu = None , broadening = None , mesh = None , with_Sigma = True , with_dc = True , proj_type = None , dosocc = False , save_to_file = True ) :
2015-03-12 00:01:12 +01:00
"""
2023-06-07 15:53:19 +02:00
Calculates the density of states and the projected density of states .
The basis of the projected density of states is specified by proj_type .
2015-03-12 00:01:12 +01:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
mu : double , optional
Chemical potential , overrides the one stored in the hdf5 archive .
2023-06-19 15:34:47 +02:00
By default , this is automatically set to the chemical potential within the SK object .
2023-04-15 00:43:23 +02:00
broadening : double , optional
2023-06-19 15:34:47 +02:00
Lorentzian broadening of the spectra to avoid any numerical artifacts .
If not given , standard value of lattice_gf ( 0.001 eV ) is used .
2023-04-15 00:43:23 +02:00
mesh : real frequency MeshType , optional
Omega mesh for the real - frequency Green ' s function.
Given as parameter to lattice_gf .
with_Sigma : boolean , optional
If True , the self energy is used for the calculation .
If false , the DOS is calculated without self energy .
2023-06-19 15:34:47 +02:00
Both with_Sigma and with_dc equal to True is needed for DFT + DMFT A ( w ) calculated .
Both with_Sigma and with_dc equal to false is needed for DFT A ( w ) calculated .
2023-04-15 00:43:23 +02:00
with_dc : boolean , optional
If True the double counting correction is used .
proj_type : string , optional
2023-06-19 15:34:47 +02:00
The type of projection used for the orbital - projected DOS .
These projected spectral functions will be determined alongside the total spectral function .
By default , no projected DOS type will be calculated ( the corresponding projected arrays will be empty ) .
The following options are :
' None ' - Only total DOS calculated
2023-04-15 00:43:23 +02:00
' wann ' - Wannier DOS calculated from the Wannier projectors
' vasp ' - Vasp orbital - projected DOS only from Vasp inputs
' wien2k ' - Wien2k orbital - projected DOS from the wien2k theta projectors
dosocc : boolean , optional
If True , the occupied DOS , DOSproj and DOSproj_orb will be returned .
The prerequisite of this option is to have calculated the band - resolved
2023-06-19 15:34:47 +02:00
density matrices generated by the occupations ( ) routine .
2015-03-12 00:01:12 +01:00
save_to_file : boolean , optional
If True , text files with the calculated data will be created .
2023-06-09 17:47:36 +02:00
2015-03-12 00:01:12 +01:00
Returns
- - - - - - -
2023-04-15 00:43:23 +02:00
DOS : Dict of numpy arrays
2023-06-19 15:34:47 +02:00
Contains the full density of states with the form of DOS [ spn ] [ n_om ] where " spn " speficies the spin type of the calculation ( " up " , " down " , or combined " ud " which relates to calculations with spin - orbit coupling ) and " n_om " is the number of real frequencies as specified by the real frequency MeshType used in the calculation . This array gives the total density of states .
2023-04-15 00:43:23 +02:00
DOSproj : Dict of numpy arrays
2023-06-19 15:34:47 +02:00
DOS projected to atom ( shell ) with the form of DOSproj [ n_shells ] [ spn ] [ n_om ] where " n_shells " is the total number of correlated or uncorrelated shells ( depending on the input " proj_type " ) . This array gives the trace of the orbital - projected density of states . Empty if proj_type = None
DOSproj_orb : Dict of numpy arrays
Orbital - projected DOS projected to atom ( shell ) and resolved into orbital contributions with the form of DOSproj_orb [ n_shells ] [ spn ] [ n_om , dim , dim ] where " dim " specifies the orbital dimension of the correlated / uncorrelated shell ( depending on the input " proj_type " ) .
2023-04-15 00:43:23 +02:00
Empty if proj_type = None
2023-06-19 15:34:47 +02:00
The output files ( if ` save_to_file = True ` ) have two ( three in the orbital - resolved case ) columns representing the frequency and real part of the DOS ( and imaginary part of the DOS ) in that order .
The output files are as follows :
DOS_ ( spn ) . dat : The total DOS .
DOS_ ( proj_type ) _ ( spn ) _proj ( i ) . dat : The DOS projected to an orbital with index i which refers to the index given in SK . shells ( or SK . corr_shells for proj_type = " wann " ) .
DOS_ ( proj_type ) _ ( sp ) _proj ( i ) _ ( m ) _ ( n ) . dat : As above , but printed as orbitally - resolved matrix in indices " m " and " n " . For example , for " d " orbitals , it gives the DOS separately for each orbital ( e . g . , " d_(xy) " , " d_(x^2-y^2), and so on).
2015-03-12 00:01:12 +01:00
"""
2023-06-19 15:34:47 +02:00
# Note the proj_type = 'elk' (- Elk orbital-projected DOS only from Elk inputs) is not included for now.
# Brief description to why can be found in the comment above the currently commented out dft_band_characters() routine
# in converters/elk.py.
# code left here just in case it will be reused.
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
# assert proj_type in ('wann', 'vasp','wien2k','elk'), "'proj_type' must be either 'wann', 'vasp', 'wien2k', or 'elk'"
assert proj_type in ( ' wann ' , ' vasp ' , ' wien2k ' ,
) , " ' proj_type ' must be either ' wann ' , ' vasp ' , ' wien2k ' "
if ( proj_type != ' wann ' ) :
assert proj_type == self . dft_code , " proj_type must be from the corresponding dft inputs. "
2023-04-15 00:43:23 +02:00
if ( with_Sigma ) :
2023-06-07 15:53:19 +02:00
assert isinstance (
self . Sigma_imp [ 0 ] . mesh , MeshReFreq ) , " SumkDFT.mesh must be real if with_Sigma is True "
mesh = self . Sigma_imp [ 0 ] . mesh
2023-04-15 00:43:23 +02:00
elif mesh is not None :
assert isinstance ( mesh , MeshReFreq ) , " mesh must be of form MeshReFreq "
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2023-04-15 00:43:23 +02:00
elif self . mesh is not None :
assert isinstance ( self . mesh , MeshReFreq ) , " self.mesh must be of form MeshReFreq "
2023-06-07 15:53:19 +02:00
mesh = self . mesh
2023-04-15 00:43:23 +02:00
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2016-05-09 10:19:56 +02:00
else :
2023-04-15 00:43:23 +02:00
assert 0 , " ReFreqMesh input required for calculations without real frequency self-energy "
2023-06-07 15:53:19 +02:00
mesh_val = numpy . linspace ( mesh . w_min , mesh . w_max , len ( mesh ) )
2023-04-15 00:43:23 +02:00
n_om = len ( mesh )
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
2023-06-07 15:53:19 +02:00
# Read in occupations from HDF5 file if required
if ( dosocc ) :
2023-04-15 00:43:23 +02:00
mpi . report ( ' Reading occupations generated by self.occupations(). ' )
thingstoread = [ ' occik ' ]
subgroup_present , values_not_read = self . read_input_from_hdf (
subgrp = self . misc_data , things_to_read = thingstoread )
if len ( values_not_read ) > 0 and mpi . is_master_node :
2023-06-07 15:53:19 +02:00
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2013-07-23 19:49:42 +02:00
2023-06-07 15:53:19 +02:00
# initialise projected DOS type if required
2023-04-15 00:43:23 +02:00
spn = self . spin_block_names [ self . SO ]
2023-06-07 15:53:19 +02:00
n_shells = 1
2023-04-15 00:43:23 +02:00
if ( proj_type == ' wann ' ) :
n_shells = self . n_corr_shells
gf_struct = self . gf_struct_sumk . copy ( )
dims = [ self . corr_shells [ ish ] [ ' dim ' ] for ish in range ( n_shells ) ]
shells_type = ' corr '
elif ( proj_type == ' vasp ' ) :
2023-06-07 15:53:19 +02:00
n_shells = 1
2023-04-15 00:43:23 +02:00
gf_struct = [ [ ( sp , list ( range ( self . proj_mat_csc . shape [ 2 ] ) ) ) for sp in spn ] ]
dims = [ self . proj_mat_csc . shape [ 2 ] ]
shells_type = ' csc '
elif ( proj_type == ' wien2k ' ) :
self . load_parproj ( )
n_shells = self . n_shells
gf_struct = [ [ ( sp , self . shells [ ish ] [ ' dim ' ] ) for sp in spn ]
2023-06-07 15:53:19 +02:00
for ish in range ( n_shells ) ]
2023-04-15 00:43:23 +02:00
dims = [ self . shells [ ish ] [ ' dim ' ] for ish in range ( n_shells ) ]
shells_type = ' all '
2023-06-07 15:53:19 +02:00
# #commented out for now - unsure this produces DFT+DMFT PDOS
2023-04-15 00:43:23 +02:00
# elif (proj_type == 'elk'):
# n_shells = self.n_shells
# dims = [self.shells[ish]['dim'] for ish in range(n_shells)]
# gf_struct = [[(sp, self.shells[ish]['dim']) for sp in spn]
# for ish in range(n_shells)]
# things_to_read = ['band_dens_muffin']
# subgroup_present, values_not_read = self.read_input_from_hdf(
# subgrp=self.bc_data, things_to_read=things_to_read)
# if len(values_not_read) > 0 and mpi.is_master_node:
# raise ValueError(
# 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
2023-06-07 15:53:19 +02:00
# set-up output arrays
DOS = { sp : numpy . zeros ( [ n_om ] , float ) for sp in spn }
2023-04-15 00:43:23 +02:00
DOSproj = [ { } for ish in range ( n_shells ) ]
DOSproj_orb = [ { } for ish in range ( n_shells ) ]
2023-06-07 15:53:19 +02:00
# set-up Green's function object
2023-04-15 00:43:23 +02:00
if ( proj_type != None ) :
G_loc = [ ]
for ish in range ( n_shells ) :
glist = [ GfReFreq ( target_shape = ( block_dim , block_dim ) , mesh = mesh )
2023-06-07 15:53:19 +02:00
for block , block_dim in gf_struct [ ish ] ]
2023-04-15 00:43:23 +02:00
G_loc . append (
BlockGf ( name_list = spn , block_list = glist , make_copies = False ) )
G_loc [ ish ] . zero ( )
dim = dims [ ish ]
for sp in spn :
DOSproj [ ish ] [ sp ] = numpy . zeros ( [ n_om ] , float )
DOSproj_orb [ ish ] [ sp ] = numpy . zeros (
[ n_om , dim , dim ] , complex )
2023-06-07 15:53:19 +02:00
# calculate the DOS
2020-04-08 21:35:59 +02:00
ikarray = numpy . array ( list ( range ( self . n_k ) ) )
2015-03-13 11:10:38 +01:00
for ik in mpi . slice_array ( ikarray ) :
2016-05-09 10:19:56 +02:00
G_latt_w = self . lattice_gf (
2021-08-27 23:31:07 +02:00
ik = ik , mu = mu , broadening = broadening , mesh = mesh , with_Sigma = with_Sigma , with_dc = with_dc )
2014-12-07 00:29:39 +01:00
G_latt_w * = self . bz_weights [ ik ]
2023-06-07 15:53:19 +02:00
# output occupied DOS if nk inputted
if ( dosocc ) :
2023-04-15 00:43:23 +02:00
for bname , gf in G_latt_w :
2023-06-07 15:53:19 +02:00
G_latt_w [ bname ] . data [ : , : , : ] * = self . occik [ bname ] [ ik ]
# DOS
2023-04-15 00:43:23 +02:00
for bname , gf in G_latt_w :
DOS [ bname ] - = gf . data . imag . trace ( axis1 = 1 , axis2 = 2 ) / numpy . pi
2015-03-13 11:10:38 +01:00
# Projected DOS:
2023-04-15 00:43:23 +02:00
if ( proj_type != None ) :
for ish in range ( n_shells ) :
tmp = G_loc [ ish ] . copy ( )
tmp . zero ( )
tmp << self . proj_type_G_loc ( G_latt_w , tmp , ik , ish , proj_type )
G_loc [ ish ] + = tmp
mpi . barrier ( )
2013-07-23 19:49:42 +02:00
2015-03-13 11:10:38 +01:00
# Collect data from mpi:
for bname in DOS :
2023-06-02 16:51:48 +02:00
DOS [ bname ] = mpi . all_reduce ( DOS [ bname ] )
2023-04-15 00:43:23 +02:00
# Collect data from mpi and put in projected arrays
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
2023-04-15 00:43:23 +02:00
for ish in range ( n_shells ) :
2023-06-05 15:53:08 +02:00
G_loc [ ish ] << mpi . all_reduce ( G_loc [ ish ] )
2015-03-13 11:10:38 +01:00
# Symmetrize and rotate to local coord. system if needed:
2023-06-07 15:53:19 +02:00
if ( ( proj_type != ' vasp ' ) and ( proj_type != ' elk ' ) ) :
2023-04-15 00:43:23 +02:00
if self . symm_op != 0 :
2023-06-07 15:53:19 +02:00
if proj_type == ' wann ' :
2023-04-15 00:43:23 +02:00
G_loc = self . symmcorr . symmetrize ( G_loc )
else :
G_loc = self . symmpar . symmetrize ( G_loc )
if self . use_rotations :
for ish in range ( n_shells ) :
for bname , gf in G_loc [ ish ] :
G_loc [ ish ] [ bname ] << self . rotloc (
2023-06-07 15:53:19 +02:00
ish , gf , direction = ' toLocal ' , shells = shells_type )
2015-03-13 11:10:38 +01:00
# G_loc can now also be used to look at orbitally-resolved quantities
2023-04-15 00:43:23 +02:00
for ish in range ( n_shells ) :
for bname , gf in G_loc [ ish ] : # loop over spins
2021-10-13 20:09:48 +02:00
DOSproj [ ish ] [ bname ] = - gf . data . imag . trace ( axis1 = 1 , axis2 = 2 ) / numpy . pi
2016-05-09 10:19:56 +02:00
DOSproj_orb [ ish ] [ bname ] [
2023-06-07 15:53:19 +02:00
: , : , : ] + = ( 1.0 j * ( gf - gf . conjugate ( ) . transpose ( ) ) / 2.0 / numpy . pi ) . data [ : , : , : ]
2013-07-23 19:49:42 +02:00
2015-03-13 11:10:38 +01:00
# Write to files
if save_to_file and mpi . is_master_node ( ) :
2023-04-15 00:43:23 +02:00
for sp in spn :
f = open ( ' DOS_ %s .dat ' % sp , ' w ' )
2016-05-09 10:19:56 +02:00
for iom in range ( n_om ) :
2023-04-15 00:43:23 +02:00
f . write ( " %s %s \n " % ( mesh_val [ iom ] , DOS [ sp ] [ iom ] ) )
2014-09-22 19:24:33 +02:00
f . close ( )
2015-03-13 11:10:38 +01:00
# Partial
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
2023-04-15 00:43:23 +02:00
for ish in range ( n_shells ) :
2023-06-07 15:53:19 +02:00
f = open ( ' DOS_ ' + proj_type + ' _ %s _proj %s .dat ' % ( sp , ish ) , ' w ' )
2016-05-09 10:19:56 +02:00
for iom in range ( n_om ) :
f . write ( " %s %s \n " %
2023-04-15 00:43:23 +02:00
( mesh_val [ iom ] , DOSproj [ ish ] [ sp ] [ iom ] ) )
2014-09-22 19:24:33 +02:00
f . close ( )
2015-03-13 11:10:38 +01:00
# Orbitally-resolved
2023-04-15 00:43:23 +02:00
for i in range ( dims [ ish ] ) :
for j in range ( dims [ ish ] ) :
2023-06-07 15:53:19 +02:00
# For Elk with parproj - skip off-diagonal elements
# if(proj_type=='elk') and (i!=j): continue
2023-04-15 00:43:23 +02:00
f = open ( ' DOS_ ' + proj_type + ' _ ' + sp + ' _proj ' + str ( ish ) +
2016-05-09 10:19:56 +02:00
' _ ' + str ( i ) + ' _ ' + str ( j ) + ' .dat ' , ' w ' )
for iom in range ( n_om ) :
2016-09-22 18:32:45 +02:00
f . write ( " %s %s %s \n " % (
2023-06-07 15:53:19 +02:00
mesh_val [ iom ] , DOSproj_orb [ ish ] [ sp ] [ iom , i , j ] . real , DOSproj_orb [ ish ] [ sp ] [ iom , i , j ] . imag ) )
2013-07-23 19:49:42 +02:00
f . close ( )
2015-03-13 11:10:38 +01:00
return DOS , DOSproj , DOSproj_orb
2020-10-09 14:35:28 +02:00
2023-04-15 00:43:23 +02:00
def proj_type_G_loc ( self , G_latt , G_inp , ik , ish , proj_type = None ) :
2019-11-21 21:34:37 +01:00
"""
2023-04-15 00:43:23 +02:00
Internal routine which calculates the project Green ' s function subject to the
proj_type input .
2019-11-21 21:34:37 +01:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
G_latt : Gf
block of lattice Green ' s functions to be projected/downfolded
G_inp : Gf
block of local Green ' s functions used as a template for G_proj
ik : integer
integer specifing k - point index .
ish : integer
integer specifing shell index .
proj_type : string , optional
Output the orbital - projected DOS type from the following options :
' wann ' - Wannier DOS calculated from the Wannier projectors
' vasp ' - Vasp orbital - projected DOS only from Vasp inputs
' wien2k ' - Wien2k orbital - projected DOS from the wien2k theta projectors
2023-06-19 15:34:47 +02:00
2019-11-21 21:34:37 +01:00
Returns
- - - - - - -
2023-04-15 00:43:23 +02:00
G_proj : Gf
projected / downfolded lattice Green ' s function
Contains the band - resolved density matrices per k - point .
2019-11-21 21:34:37 +01:00
"""
2023-06-19 15:34:47 +02:00
# Note the proj_type = 'elk' (- Elk orbital-projected DOS only from Elk inputs) is not included for now.
# Brief description to why can be found in the comment above the currently commented out dft_band_characters() routine
# in converters/elk.py.
# code left here just in case it will be reused.
2023-04-15 00:43:23 +02:00
G_proj = G_inp . copy ( )
if ( proj_type == ' wann ' ) :
for bname , gf in G_proj :
G_proj [ bname ] << self . downfold ( ik , ish , bname ,
2023-06-07 15:53:19 +02:00
G_latt [ bname ] , gf ) # downfolding G
2023-04-15 00:43:23 +02:00
elif ( proj_type == ' vasp ' ) :
for bname , gf in G_latt :
G_proj [ bname ] << self . downfold ( ik , ish , bname , gf , G_proj [ bname ] , shells = ' csc ' )
elif ( proj_type == ' wien2k ' ) :
tmp = G_proj . copy ( )
for ir in range ( self . n_parproj [ ish ] ) :
tmp . zero ( )
for bname , gf in tmp :
tmp [ bname ] << self . downfold ( ik , ish , bname ,
2023-06-07 15:53:19 +02:00
G_latt [ bname ] , gf , shells = ' all ' , ir = ir )
2023-04-15 00:43:23 +02:00
G_proj + = tmp
# elif (proj_type == 'elk'):
# dim = self.shells[ish]['dim']
# ntoi = self.spin_names_to_ind[self.SO]
# for bname, gf in G_latt:
# n_om = len(gf.data[:,0,0])
# isp=ntoi[bname]
# nst=self.n_orbitals[ik,isp]
# #matrix multiply band resolved muffin density with
# #diagonal of band resolved spectral function and fill diagonal of
# #DOSproj_orb orbital dimensions with the result for each frequency
# bdm=self.band_dens_muffin[ik,isp,ish,0:dim,0:nst]
# tmp=[numpy.matmul(bdm, gf.data[iom,:,:].diagonal())
# for iom in range(n_om)]
# tmp=numpy.asarray(tmp)
# tmp2 = numpy.zeros([n_om,dim,dim], dtype=complex)
# if(dim==1):
# tmp2[:,0,0]=tmp[:,0]
# else:
# [numpy.fill_diagonal(tmp2[iom,:,:],tmp[iom,:])
# for iom in range(n_om)]
# G_proj[bname].data[:,:,:] = tmp2[:,:,:]
return G_proj
2023-06-07 15:53:19 +02:00
def load_parproj ( self , data_type = None ) :
2015-03-12 00:01:12 +01:00
"""
2023-04-15 00:43:23 +02:00
Internal routine which loads the n_parproj , proj_mat_all , rot_mat_all and
rot_mat_all_time_inv from parproj data from . h5 file .
2015-03-12 00:01:12 +01:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
data_type : string , optional
which data type desired to be read in .
' band ' - reads data converted by bands_convert ( )
None - reads data converted by parproj_convert ( )
2015-03-12 00:01:12 +01:00
"""
2023-06-07 15:53:19 +02:00
# read in the projectors
2023-04-15 00:43:23 +02:00
things_to_read = [ ' n_parproj ' , ' proj_mat_all ' ]
if data_type == ' band ' :
subgroup_present , values_not_read = self . read_input_from_hdf (
2023-06-07 15:53:19 +02:00
subgrp = self . bands_data , things_to_read = things_to_read )
2023-04-15 00:43:23 +02:00
else :
subgroup_present , values_not_read = self . read_input_from_hdf (
2023-06-07 15:53:19 +02:00
subgrp = self . parproj_data , things_to_read = things_to_read )
2023-04-15 00:43:23 +02:00
if self . symm_op :
self . symmpar = Symmetry ( self . hdf_file , subgroup = self . symmpar_data )
if len ( values_not_read ) > 0 and mpi . is_master_node :
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2023-06-07 15:53:19 +02:00
# read general data
2023-04-15 00:43:23 +02:00
things_to_read = [ ' rot_mat_all ' , ' rot_mat_all_time_inv ' ]
2023-03-23 22:40:13 +01:00
subgroup_present , values_not_read = self . read_input_from_hdf (
2016-05-09 10:19:56 +02:00
subgrp = self . parproj_data , things_to_read = things_to_read )
2023-03-23 22:40:13 +01:00
if len ( values_not_read ) > 0 and mpi . is_master_node :
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2013-07-23 19:49:42 +02:00
2023-04-15 00:43:23 +02:00
def occupations ( self , mu = None , with_Sigma = True , with_dc = True , save_occ = True ) :
2020-10-09 14:35:28 +02:00
"""
2023-04-15 00:43:23 +02:00
Calculates the band resolved density matrices ( occupations ) from the Matsubara
frequency self - energy .
2023-06-07 15:53:19 +02:00
2020-10-09 14:35:28 +02:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
mu : double , optional
Chemical potential , overrides the one stored in the hdf5 archive .
with_Sigma : boolean , optional
If True , the self energy is used for the calculation .
If false , the DOS is calculated without self energy .
with_dc : boolean , optional
If True the double counting correction is used .
save_occ : boolean , optional
If True , saves the band resolved density matrix in misc_data .
2020-10-09 14:35:28 +02:00
save_to_file : boolean , optional
2023-06-07 15:53:19 +02:00
If True , text files with the calculated data will be created . \
2020-10-09 14:35:28 +02:00
Returns
- - - - - - -
2023-04-15 00:43:23 +02:00
occik : Dict of numpy arrays
Contains the band - resolved density matrices per k - point .
2020-10-09 14:35:28 +02:00
"""
2023-04-15 00:43:23 +02:00
if with_Sigma :
mesh = self . Sigma_imp [ 0 ] . mesh
2020-10-09 14:35:28 +02:00
else :
2023-04-15 00:43:23 +02:00
mesh = self . mesh
2023-06-07 15:53:19 +02:00
assert isinstance (
mesh , MeshImFreq ) , " SumkDFT.mesh must be real if with_Sigma is True or mesh is not given "
2023-04-15 00:43:23 +02:00
2020-10-09 14:35:28 +02:00
if mu is None :
mu = self . chemical_potential
2023-04-15 00:43:23 +02:00
ntoi = self . spin_names_to_ind [ self . SO ]
2020-10-09 14:35:28 +02:00
spn = self . spin_block_names [ self . SO ]
2023-04-15 00:43:23 +02:00
occik = { }
for sp in spn :
2023-06-07 15:53:19 +02:00
# same format as gf.data ndarray
occik [ sp ] = [ numpy . zeros ( [ 1 , self . n_orbitals [ ik , ntoi [ sp ] ] ,
self . n_orbitals [ ik , ntoi [ sp ] ] ] , numpy . double ) for ik in range ( self . n_k ) ]
# calculate the occupations
2020-10-09 14:35:28 +02:00
ikarray = numpy . array ( range ( self . n_k ) )
2023-04-15 00:43:23 +02:00
for ik in range ( self . n_k ) :
G_latt = self . lattice_gf (
ik = ik , mu = mu , with_Sigma = with_Sigma , with_dc = with_dc )
for bname , gf in G_latt :
2023-06-07 15:53:19 +02:00
occik [ bname ] [ ik ] [ 0 , : , : ] = gf . density ( ) . real
2020-10-09 14:35:28 +02:00
# Collect data from mpi:
2023-04-15 00:43:23 +02:00
for sp in spn :
2023-06-05 15:53:08 +02:00
occik [ sp ] = mpi . all_reduce ( occik [ sp ] )
2023-04-15 00:43:23 +02:00
mpi . barrier ( )
2023-06-07 15:53:19 +02:00
# save to HDF5 file (if specified)
2023-04-15 00:43:23 +02:00
if save_occ and mpi . is_master_node ( ) :
things_to_save_misc = [ ' occik ' ]
# Save it to the HDF:
ar = HDFArchive ( self . hdf_file , ' a ' )
if not ( self . misc_data in ar ) :
ar . create_group ( self . misc_data )
for it in things_to_save_misc :
ar [ self . misc_data ] [ it ] = locals ( ) [ it ]
del ar
return occik
2020-10-09 14:35:28 +02:00
2023-04-15 00:43:23 +02:00
def spectral_contours ( self , mu = None , broadening = None , mesh = None , plot_range = None , FS = True , with_Sigma = True , with_dc = True , proj_type = None , save_to_file = True ) :
2020-10-09 14:35:28 +02:00
"""
2023-04-15 00:43:23 +02:00
Calculates the correlated spectral function at the Fermi level ( relating to the Fermi
surface ) or at specific frequencies .
2020-10-09 14:35:28 +02:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
mu : double , optional
Chemical potential , overrides the one stored in the hdf5 archive .
2023-06-19 15:34:47 +02:00
By default , this is automatically set to the chemical potential within the SK object .
2023-04-15 00:43:23 +02:00
broadening : double , optional
2023-06-19 15:34:47 +02:00
Lorentzian broadening of the spectra to avoid any numerical artifacts .
If not given , standard value of lattice_gf ( 0.001 eV ) is used .
2023-04-15 00:43:23 +02:00
mesh : real frequency MeshType , optional
Omega mesh for the real - frequency Green ' s function.
Given as parameter to lattice_gf .
2023-06-19 15:34:47 +02:00
plot_shift : double , optional
Offset [ = ( ik - 1 ) * plot_shift , where ik is the index of the k - point ] for each A ( k , w ) for stacked plotting of spectra .
2023-04-15 00:43:23 +02:00
plot_range : list of double , optional
2023-06-19 15:34:47 +02:00
Sets the energy window for plotting to ( plot_range [ 0 ] , plot_range [ 1 ] ) .
If not provided , the min and max values of the energy mesh is used .
2023-04-15 00:43:23 +02:00
FS : boolean
Flag for calculating the spectral function at the Fermi level ( omega ~ 0 )
If False , the spectral function will be generated for each frequency within
2023-06-19 15:34:47 +02:00
plot_range .
2023-04-15 00:43:23 +02:00
with_Sigma : boolean , optional
If True , the self energy is used for the calculation .
If false , the DOS is calculated without self energy .
2023-06-19 15:34:47 +02:00
Both with_Sigma and with_dc equal to True is needed for DFT + DMFT A ( k , w ) calculated .
Both with_Sigma and with_dc equal to false is needed for DFT A ( k , w ) calculated .
2023-04-15 00:43:23 +02:00
with_dc : boolean , optional
If True the double counting correction is used .
proj_type : string , optional
2023-06-19 15:34:47 +02:00
The type of projection used for the orbital - projected DOS .
These projected spectral functions will be determined alongside the total spectral function .
By default , no projected DOS type will be calculated ( the corresponding projected arrays will be empty ) .
The following options are :
' None ' - Only total DOS calculated
' wann ' - Wannier DOS calculated from the Wannier projectors
2020-10-09 14:35:28 +02:00
save_to_file : boolean , optional
If True , text files with the calculated data will be created .
2023-06-19 15:34:47 +02:00
2020-10-09 14:35:28 +02:00
Returns
- - - - - - -
2023-04-15 00:43:23 +02:00
Akw : Dict of numpy arrays
2023-06-19 15:34:47 +02:00
( Correlated ) k - resolved spectral function .
This dictionary has the form of ` Akw [ spn ] [ n_k , n_om ] ` where spn , n_k and n_om are the spin , number of k - points , and number of frequencies used in the calculation .
pAkw : Dict of numpy arrays
( Correlated ) k - resolved spectral function projected to atoms ( i . e . , the Trace of the orbital - projected A ( k , w ) ) .
This dictionary has the form of pAkw [ n_shells ] [ spn ] [ n_k , n_om ] where n_shells is the total number of correlated or uncorrelated shells .
Empty if proj_type = None
pAkw_orb : Dict of numpy arrays
( Correlated ) k - resolved spectral function projected to atoms and
resolved into orbital contributions .
This dictionary has the form of pAkw [ n_shells ] [ spn ] [ n_k , n_om , dim , dim ] where dim specifies the orbital dimension of the correlated / uncorrelated shell .
2023-04-15 00:43:23 +02:00
Empty if proj_type = None
2023-06-19 15:34:47 +02:00
The output files have three columns representing the k - point index , frequency and A ( k , w ) in that order . The output files are as follows :
Akw_ ( sp ) . dat : The total A ( k , w ) .
Akw_ ( proj_type ) _ ( spn ) _proj ( i ) . dat : The A ( k , w ) projected to shell with index ( i ) .
Akw_ ( proj_type ) _ ( spn ) _proj ( i ) _ ( m ) _ ( n ) . dat ` : As above , but for each ( m ) and ( n ) orbital contribution .
The files are prepended with either of the following :
For ` FS ` set to True : The output files ' name include _FS_ and these files contain four columns which are the cartesian reciprocal coordinates (kx, ky, kz) and Akw.
For ` FS ` set to False : The output files ' name include _omega_(iom) (with `iom` being the frequency mesh index). These files also contain four columns as described above along with a comment at the top of the file which gives the frequency value at which the spectral function was evaluated.
2020-10-09 14:35:28 +02:00
"""
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
2023-04-15 00:43:23 +02:00
assert proj_type in ( ' wann ' ) , " ' proj_type ' must be ' wann ' if not None "
2023-06-07 15:53:19 +02:00
# read in the energy contour energies and projectors
things_to_read = [ ' n_k ' , ' bmat ' , ' BZ_n_k ' , ' BZ_iknr ' , ' BZ_vkl ' ,
2020-10-09 14:35:28 +02:00
' n_orbitals ' , ' proj_mat ' , ' hopping ' ]
2023-03-23 22:40:13 +01:00
subgroup_present , values_not_read = self . read_input_from_hdf (
2023-06-07 15:53:19 +02:00
subgrp = self . cont_data , things_to_read = things_to_read )
2023-03-23 22:40:13 +01:00
if len ( values_not_read ) > 0 and mpi . is_master_node :
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2020-10-09 14:35:28 +02:00
if mu is None :
mu = self . chemical_potential
2023-04-15 00:43:23 +02:00
if ( with_Sigma ) :
2023-06-07 15:53:19 +02:00
assert isinstance (
self . Sigma_imp [ 0 ] . mesh , MeshReFreq ) , " SumkDFT.mesh must be real if with_Sigma is True "
mesh = self . Sigma_imp [ 0 ] . mesh
2023-04-15 00:43:23 +02:00
elif mesh is not None :
assert isinstance ( mesh , MeshReFreq ) , " mesh must be of form MeshReFreq "
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2023-04-15 00:43:23 +02:00
elif self . mesh is not None :
assert isinstance ( self . mesh , MeshReFreq ) , " self.mesh must be of form MeshReFreq "
2023-06-07 15:53:19 +02:00
mesh = self . mesh
2023-04-15 00:43:23 +02:00
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2023-04-15 00:43:23 +02:00
else :
assert 0 , " ReFreqMesh input required for calculations without real frequency self-energy "
2023-06-07 15:53:19 +02:00
mesh_val = numpy . linspace ( mesh . w_min , mesh . w_max , len ( mesh ) )
2023-04-15 00:43:23 +02:00
n_om = len ( mesh )
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
2023-06-07 15:53:19 +02:00
# for Fermi Surface calculations
2023-04-15 00:43:23 +02:00
if FS :
dw = abs ( mesh_val [ 1 ] - mesh_val [ 0 ] )
2023-06-07 15:53:19 +02:00
# ensure that a few frequencies around the Fermi level are included
2023-04-15 00:43:23 +02:00
plot_range = [ - 2 * dw , 2 * dw ]
mpi . report ( ' Generated A(k,w) will be evaluted at closest frequency to 0.0 in given mesh ' )
if plot_range is None :
2023-06-07 15:53:19 +02:00
n_om = len ( mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ] )
mesh_val2 = mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ]
2023-04-15 00:43:23 +02:00
else :
om_minplot = plot_range [ 0 ]
om_maxplot = plot_range [ 1 ]
2023-06-07 15:53:19 +02:00
n_om = len ( mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ] )
mesh_val2 = mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ]
# \omega ~= 0.0 index for FS file
2023-04-15 00:43:23 +02:00
abs_mesh_val = [ abs ( i ) for i in mesh_val2 ]
2023-06-07 15:53:19 +02:00
jw = [ i for i in range ( len ( abs_mesh_val ) ) if abs_mesh_val [ i ]
== numpy . min ( abs_mesh_val [ : ] ) ]
# calculate the spectral functions for the irreducible set of k-points
[ Akw , pAkw , pAkw_orb ] = self . gen_Akw ( mu = mu , broadening = broadening , mesh = mesh ,
plot_shift = 0.0 , plot_range = plot_range ,
shell_list = None , with_Sigma = with_Sigma , with_dc = with_dc ,
proj_type = proj_type )
2020-10-09 14:35:28 +02:00
if save_to_file and mpi . is_master_node ( ) :
2023-04-15 00:43:23 +02:00
spn = self . spin_block_names [ self . SO ]
vkc = numpy . zeros ( 3 , float )
2023-06-07 15:53:19 +02:00
mesh_val2 = mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ]
2023-04-15 00:43:23 +02:00
if FS :
n_om = 1
else :
n_om = len ( mesh_val2 )
for sp in spn :
2023-06-07 15:53:19 +02:00
# Open file for storage:
for iom in range ( n_om ) :
2023-04-15 00:43:23 +02:00
if FS :
f = open ( ' Akw_FS_ ' + sp + ' .dat ' , ' w ' )
2023-06-07 15:53:19 +02:00
jom = jw [ 0 ]
2023-04-15 00:43:23 +02:00
else :
f = open ( ' Akw_omega_ %s _ %s .dat ' % ( iom , sp ) , ' w ' )
2023-06-07 15:53:19 +02:00
jom = iom
f . write ( " #Spectral function evaluated at frequency = %s \n " % mesh_val2 [ jom ] )
2023-04-15 00:43:23 +02:00
for ik in range ( self . BZ_n_k ) :
2023-06-07 15:53:19 +02:00
jk = self . BZ_iknr [ ik ]
vkc [ : ] = numpy . matmul ( self . bmat , self . BZ_vkl [ ik , : ] )
2023-04-15 00:43:23 +02:00
f . write ( " %s %s %s %s \n " % ( vkc [ 0 ] , vkc [ 1 ] , vkc [ 2 ] , Akw [ sp ] [ jk , jom ] ) )
f . close ( )
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
n_shells = len ( pAkw [ : ] )
for iom in range ( n_om ) :
2023-04-15 00:43:23 +02:00
for sp in spn :
for ish in range ( n_shells ) :
if FS :
strng = ' Akw_FS ' + ' _ ' + proj_type + ' _ ' + sp + ' _proj ' + str ( ish )
2023-06-07 15:53:19 +02:00
jom = jw [ 0 ]
2023-04-15 00:43:23 +02:00
else :
2023-06-07 15:53:19 +02:00
strng = ' Akw_omega_ ' + str ( iom ) + ' _ ' + proj_type + \
' _ ' + sp + ' _proj ' + str ( ish )
jom = iom
2023-04-15 00:43:23 +02:00
f = open ( strng + ' .dat ' , ' w ' )
2023-06-07 15:53:19 +02:00
f . write ( " #Spectral function evaluated at frequency = %s \n " % mesh_val2 [ jom ] )
2023-04-15 00:43:23 +02:00
for ik in range ( self . BZ_n_k ) :
2023-06-07 15:53:19 +02:00
jk = self . BZ_iknr [ ik ]
vkc [ : ] = numpy . matmul ( self . bmat , self . BZ_vkl [ ik , : ] )
2023-04-15 00:43:23 +02:00
f . write ( " %s %s %s %s \n " % ( vkc [ 0 ] , vkc [ 1 ] , vkc [ 2 ] ,
2023-06-07 15:53:19 +02:00
pAkw [ ish ] [ sp ] [ jk , jom ] ) )
2023-04-15 00:43:23 +02:00
f . close ( )
2023-06-07 15:53:19 +02:00
dim = len ( pAkw_orb [ ish ] [ sp ] [ 0 , 0 , 0 , : ] )
2023-04-15 00:43:23 +02:00
for i in range ( dim ) :
for j in range ( dim ) :
strng2 = strng + ' _ ' + str ( i ) + ' _ ' + str ( j )
# Open file for storage:
f = open ( strng2 + ' .dat ' , ' w ' )
for ik in range ( self . BZ_n_k ) :
2023-06-07 15:53:19 +02:00
jk = self . BZ_iknr [ ik ]
vkc [ : ] = numpy . matmul ( self . bmat , self . BZ_vkl [ ik , : ] )
2023-04-15 00:43:23 +02:00
f . write ( " %s %s %s %s \n " % ( vkc [ 0 ] , vkc [ 1 ] , vkc [ 2 ] ,
2023-06-07 15:53:19 +02:00
pAkw_orb [ ish ] [ sp ] [ jk , jom , i , j ] ) )
2023-04-15 00:43:23 +02:00
f . close ( )
2020-10-09 14:35:28 +02:00
2023-04-15 00:43:23 +02:00
return Akw , pAkw , pAkw_orb
2020-10-09 14:35:28 +02:00
2023-04-15 00:43:23 +02:00
def spaghettis ( self , mu = None , broadening = None , mesh = None , plot_shift = 0.0 , plot_range = None , shell_list = None , with_Sigma = True , with_dc = True , proj_type = None , save_to_file = True ) :
2015-03-12 00:01:12 +01:00
"""
2023-04-15 00:43:23 +02:00
Calculates the k - resolved spectral function A ( k , w ) ( band structure )
2016-05-09 10:19:56 +02:00
2015-03-12 00:01:12 +01:00
Parameters
- - - - - - - - - -
2023-04-15 00:43:23 +02:00
mu : double , optional
Chemical potential , overrides the one stored in the hdf5 archive .
2023-06-19 15:34:47 +02:00
By default , this is automatically set to the chemical potential within the SK object .
2023-04-15 00:43:23 +02:00
broadening : double , optional
2023-06-19 15:34:47 +02:00
Lorentzian broadening of the spectra to avoid any numerical artifacts .
If not given , standard value of lattice_gf ( 0.001 eV ) is used .
2023-04-15 00:43:23 +02:00
mesh : real frequency MeshType , optional
2023-06-19 15:34:47 +02:00
Omega mesh for the real - frequency Green ' s function.
2023-04-15 00:43:23 +02:00
Given as parameter to lattice_gf .
plot_shift : double , optional
2023-06-19 15:34:47 +02:00
Offset [ = ( ik - 1 ) * plot_shift , where ik is the index of the k - point ] for each A ( k , w ) for stacked plotting of spectra .
2023-04-15 00:43:23 +02:00
plot_range : list of double , optional
2023-06-07 15:53:19 +02:00
Sets the energy window for plotting to ( plot_range [ 0 ] , plot_range [ 1 ] ) .
If not provided , the min and max values of the energy mesh is used .
2023-04-15 00:43:23 +02:00
shell_list : list of integers , optional
2023-06-07 15:53:19 +02:00
Contains the indices of the shells of which the projected spectral function
is calculated for .
If shell_list = None and proj_type is not None , then the projected spectral
2023-04-15 00:43:23 +02:00
function is calculated for all shells .
2023-06-19 15:34:47 +02:00
Note for experts : The spectra from Wien2k inputs are not rotated to the local coordinate system used in Wien2k .
2023-04-15 00:43:23 +02:00
with_Sigma : boolean , optional
2023-06-19 15:34:47 +02:00
If True , the self energy is used for the calculation .
2023-04-15 00:43:23 +02:00
If false , the DOS is calculated without self energy .
2023-06-19 15:34:47 +02:00
Both with_Sigma and with_dc equal to True is needed for DFT + DMFT A ( k , w ) calculated .
Both with_Sigma and with_dc equal to false is needed for DFT A ( k , w ) calculated .
2023-04-15 00:43:23 +02:00
with_dc : boolean , optional
If True the double counting correction is used .
proj_type : string , optional
2023-06-19 15:34:47 +02:00
The type of projection used for the orbital - projected DOS .
These projected spectral functions will be determined alongside the total spectral function .
By default , no projected DOS type will be calculated ( the corresponding projected arrays will be empty ) .
The following options are :
' None ' - Only total DOS calculated
' wann ' - Wannier DOS calculated from the Wannier projectors
' wien2k ' - Wien2k orbital - projected DOS from the wien2k theta projectors
2023-04-15 00:43:23 +02:00
save_to_file : boolean , optional
If True , text files with the calculated data will be created .
2023-06-19 15:34:47 +02:00
2015-03-12 00:01:12 +01:00
Returns
- - - - - - -
2023-04-15 00:43:23 +02:00
Akw : Dict of numpy arrays
2023-06-19 15:34:47 +02:00
( Correlated ) k - resolved spectral function .
This dictionary has the form of ` Akw [ spn ] [ n_k , n_om ] ` where spn , n_k and n_om are the spin , number of k - points , and number of frequencies used in the calculation .
pAkw : Dict of numpy arrays
( Correlated ) k - resolved spectral function projected to atoms ( i . e . , the Trace of the orbital - projected A ( k , w ) ) .
This dictionary has the form of pAkw [ n_shells ] [ spn ] [ n_k , n_om ] where n_shells is the total number of correlated or uncorrelated shells .
2023-04-15 00:43:23 +02:00
Empty if proj_type = None
2023-06-19 15:34:47 +02:00
pAkw_orb : Dict of numpy arrays
2023-06-07 15:53:19 +02:00
( Correlated ) k - resolved spectral function projected to atoms and
2023-06-19 15:34:47 +02:00
resolved into orbital contributions .
This dictionary has the form of pAkw [ n_shells ] [ spn ] [ n_k , n_om , dim , dim ] where dim specifies the orbital dimension of the correlated / uncorrelated shell .
Empty if proj_type = None
The output files have three columns representing the k - point index , frequency and A ( k , w ) in that order . The output files are as follows :
Akw_ ( sp ) . dat : The total A ( k , w ) .
Akw_ ( proj_type ) _ ( spn ) _proj ( i ) . dat : The A ( k , w ) projected to shell with index ( i ) .
Akw_ ( proj_type ) _ ( spn ) _proj ( i ) _ ( m ) _ ( n ) . dat ` : As above , but for each ( m ) and ( n ) orbital contribution .
2015-03-12 00:01:12 +01:00
"""
2023-06-07 15:53:19 +02:00
# initialisation
if ( proj_type != None ) :
2023-04-15 00:43:23 +02:00
assert proj_type in ( ' wann ' , ' wien2k ' ) , " ' proj_type ' must be either ' wann ' , ' wien2k ' "
2023-06-07 15:53:19 +02:00
if ( proj_type != ' wann ' ) :
assert proj_type == self . dft_code , " proj_type must be from the corresponding dft inputs. "
2023-04-15 00:43:23 +02:00
things_to_read = [ ' n_k ' , ' n_orbitals ' , ' proj_mat ' , ' hopping ' ]
2023-03-23 22:40:13 +01:00
subgroup_present , values_not_read = self . read_input_from_hdf (
2016-05-09 10:19:56 +02:00
subgrp = self . bands_data , things_to_read = things_to_read )
2023-03-23 22:40:13 +01:00
if len ( values_not_read ) > 0 and mpi . is_master_node :
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2023-06-07 15:53:19 +02:00
if ( proj_type == ' wien2k ' ) :
2023-04-15 00:43:23 +02:00
self . load_parproj ( data_type = ' band ' )
2016-05-09 10:19:56 +02:00
if mu is None :
mu = self . chemical_potential
2023-04-15 00:43:23 +02:00
if ( with_Sigma ) :
2023-06-07 15:53:19 +02:00
assert isinstance (
self . Sigma_imp [ 0 ] . mesh , MeshReFreq ) , " SumkDFT.mesh must be real if with_Sigma is True "
mesh = self . Sigma_imp [ 0 ] . mesh
2023-04-15 00:43:23 +02:00
elif mesh is not None :
assert isinstance ( mesh , MeshReFreq ) , " mesh must be of form MeshReFreq "
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2023-04-15 00:43:23 +02:00
elif self . mesh is not None :
assert isinstance ( self . mesh , MeshReFreq ) , " self.mesh must be of form MeshReFreq "
2023-06-07 15:53:19 +02:00
mesh = self . mesh
2023-04-15 00:43:23 +02:00
if broadening is None :
2023-06-07 15:53:19 +02:00
broadening = 0.001
2023-04-15 00:43:23 +02:00
else :
assert 0 , " ReFreqMesh input required for calculations without real frequency self-energy "
2023-06-07 15:53:19 +02:00
mesh_val = numpy . linspace ( mesh . w_min , mesh . w_max , len ( mesh ) )
2021-08-27 23:31:07 +02:00
n_om = len ( mesh )
2023-04-15 00:43:23 +02:00
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
2013-07-23 19:49:42 +02:00
if plot_range is None :
2023-04-15 00:43:23 +02:00
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
2013-07-23 19:49:42 +02:00
else :
om_minplot = plot_range [ 0 ]
om_maxplot = plot_range [ 1 ]
2023-06-07 15:53:19 +02:00
n_om = len ( mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ] )
2014-09-22 19:24:33 +02:00
2023-06-07 15:53:19 +02:00
[ Akw , pAkw , pAkw_orb ] = self . gen_Akw ( mu = mu , broadening = broadening , mesh = mesh ,
plot_shift = plot_shift , plot_range = plot_range ,
shell_list = shell_list , with_Sigma = with_Sigma , with_dc = with_dc ,
proj_type = proj_type )
2015-11-30 16:22:53 +01:00
2015-03-13 11:10:38 +01:00
if save_to_file and mpi . is_master_node ( ) :
2023-06-07 15:53:19 +02:00
mesh_val2 = mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ]
spn = self . spin_block_names [ self . SO ]
for sp in spn :
2023-04-15 00:43:23 +02:00
# Open file for storage:
f = open ( ' Akw_ ' + sp + ' .dat ' , ' w ' )
for ik in range ( self . n_k ) :
for iom in range ( n_om ) :
2023-06-07 15:53:19 +02:00
f . write ( ' %s %s %s \n ' % ( ik , mesh_val2 [ iom ] , Akw [ sp ] [ ik , iom ] ) )
2023-04-15 00:43:23 +02:00
f . write ( ' \n ' )
f . close ( )
2023-06-07 15:53:19 +02:00
if ( proj_type != None ) :
2023-04-15 00:43:23 +02:00
n_shells = len ( pAkw [ : ] )
2023-06-07 15:53:19 +02:00
if shell_list == None :
shell_list = [ ish for ish in range ( n_shells ) ]
2023-04-15 00:43:23 +02:00
for sp in spn :
for ish in range ( n_shells ) :
2023-06-07 15:53:19 +02:00
jsh = shell_list [ ish ]
2023-04-15 00:43:23 +02:00
f = open ( ' Akw_ ' + proj_type + ' _ ' +
2023-06-07 15:53:19 +02:00
sp + ' _proj ' + str ( jsh ) + ' .dat ' , ' w ' )
2013-07-23 19:49:42 +02:00
for ik in range ( self . n_k ) :
2023-04-15 00:43:23 +02:00
for iom in range ( n_om ) :
f . write ( ' %s %s %s \n ' % (
2023-06-07 15:53:19 +02:00
ik , mesh_val2 [ iom ] , pAkw [ ish ] [ sp ] [ ik , iom ] ) )
2023-04-15 00:43:23 +02:00
f . write ( ' \n ' )
2013-07-23 19:49:42 +02:00
f . close ( )
2023-06-07 15:53:19 +02:00
# get orbital dimension from the length of dimension of the array
dim = len ( pAkw_orb [ ish ] [ sp ] [ 0 , 0 , 0 , : ] )
2023-04-15 00:43:23 +02:00
for i in range ( dim ) :
for j in range ( dim ) :
2016-05-09 10:19:56 +02:00
# Open file for storage:
2023-04-15 00:43:23 +02:00
f = open ( ' Akw_ ' + proj_type + ' _ ' + sp + ' _proj ' + str ( jsh )
2023-06-07 15:53:19 +02:00
+ ' _ ' + str ( i ) + ' _ ' + str ( j ) + ' .dat ' , ' w ' )
2013-07-23 19:49:42 +02:00
for ik in range ( self . n_k ) :
2014-09-22 19:24:33 +02:00
for iom in range ( n_om ) :
2023-04-15 00:43:23 +02:00
f . write ( ' %s %s %s \n ' % (
2023-06-07 15:53:19 +02:00
ik , mesh_val2 [ iom ] , pAkw_orb [ ish ] [ sp ] [ ik , iom , i , j ] ) )
2013-07-23 19:49:42 +02:00
f . write ( ' \n ' )
f . close ( )
2023-04-15 00:43:23 +02:00
return Akw , pAkw , pAkw_orb
def gen_Akw ( self , mu , broadening , mesh , plot_shift , plot_range , shell_list , with_Sigma , with_dc , proj_type ) :
"""
Internal routine used by spaghettis and spectral_contours to Calculate the k - resolved spectral
function A ( k , w ) . For advanced users only .
Parameters
- - - - - - - - - -
mu : double
Chemical potential , overrides the one stored in the hdf5 archive .
broadening : double
Lorentzian broadening of the spectra .
mesh : real frequency MeshType , optional
Omega mesh for the real - frequency Green ' s function.
Given as parameter to lattice_gf .
plot_shift : double
Offset for each A ( k , w ) for stacked plotting of spectra .
plot_range : list of double
Sets the energy window for plotting to ( plot_range [ 0 ] , plot_range [ 1 ] ) .
shell_list : list of integers , optional
Contains the indices of the shells of which the projected spectral function
is calculated for .
If shell_list = None and proj_type is not None , then the projected spectral
function is calculated for all shells .
with_Sigma : boolean
If True , the self energy is used for the calculation .
If false , the DOS is calculated without self energy .
with_dc : boolean
If True the double counting correction is used .
proj_type : string
Output the orbital - projected A ( k , w ) type from the following :
' wann ' - Wannier A ( k , w ) calculated from the Wannier projectors
' wien2k ' - Wien2k orbital - projected A ( k , w ) from the wien2k theta projectors
2023-06-09 17:47:36 +02:00
2023-04-15 00:43:23 +02:00
Returns
- - - - - - -
Akw : Dict of numpy arrays
( Correlated ) k - resolved spectral function
pAkw : Dict of numpy arrays
( Correlated ) k - resolved spectral function projected to atoms .
Empty if proj_type = None
pAkw_orb : Dict of numpy arrays
( Correlated ) k - resolved spectral function projected to atoms and
resolved into orbital contributions . Empty if proj_type = None
"""
2023-06-05 15:38:57 +02:00
mesh_val = numpy . linspace ( mesh . w_min , mesh . w_max , len ( mesh ) )
2023-04-15 00:43:23 +02:00
n_om = len ( mesh )
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
if plot_range is None :
om_minplot = mesh_val [ 0 ] - 0.001
om_maxplot = mesh_val [ - 1 ] + 0.001
else :
om_minplot = plot_range [ 0 ]
om_maxplot = plot_range [ 1 ]
n_om = len ( mesh_val [ ( mesh_val > om_minplot ) & ( mesh_val < om_maxplot ) ] )
#set-up spectral functions
spn = self . spin_block_names [ self . SO ]
Akw = { sp : numpy . zeros ( [ self . n_k , n_om ] , float ) for sp in spn }
pAkw = [ ]
pAkw_orb = [ ]
#set-up projected A(k,w) and parameters if required
if ( proj_type ) :
if ( proj_type == ' wann ' ) :
n_shells = self . n_corr_shells
gf_struct = self . gf_struct_sumk . copy ( )
dims = [ self . corr_shells [ ish ] [ ' dim ' ] for ish in range ( n_shells ) ]
shells_type = ' corr '
elif ( proj_type == ' wien2k ' ) :
n_shells = self . n_shells
gf_struct = [ [ ( sp , self . shells [ ish ] [ ' dim ' ] ) for sp in spn ]
for ish in range ( n_shells ) ]
dims = [ self . shells [ ish ] [ ' dim ' ] for ish in range ( n_shells ) ]
shells_type = ' all '
#only outputting user specified parproj shells
if shell_list != None :
for ish in shell_list :
if ( ish > n_shells ) or ( ish < 0 ) :
raise IOError ( " indices in shell_list input do not correspond \
to existing self . shells indices " )
n_shells = len ( shell_list )
mpi . report ( " calculating spectral functions for following user specified shell_list: " )
[ mpi . report ( ' %s : %s ' % ( ish , self . shells [ ish ] ) ) for ish in shell_list ]
else :
shell_list = [ ish for ish in range ( n_shells ) ]
#projected Akw via projectors
pAkw = [ { } for ish in range ( n_shells ) ]
pAkw_orb = [ { } for ish in range ( n_shells ) ]
#set-up Green's function object
G_loc = [ ]
for ish in range ( n_shells ) :
jsh = shell_list [ ish ]
dim = dims [ ish ]
for sp in spn :
pAkw [ ish ] [ sp ] = numpy . zeros ( [ self . n_k , n_om ] , float )
pAkw_orb [ ish ] [ sp ] = numpy . zeros ( [ self . n_k , n_om , dim , dim ] , float )
glist = [ GfReFreq ( target_shape = ( block_dim , block_dim ) , mesh = mesh )
for block , block_dim in gf_struct [ ish ] ]
G_loc . append (
BlockGf ( name_list = spn , block_list = glist , make_copies = False ) )
G_loc [ ish ] . zero ( )
#calculate the spectral function
ikarray = numpy . array ( list ( range ( self . n_k ) ) )
for ik in mpi . slice_array ( ikarray ) :
G_latt_w = self . lattice_gf ( ik = ik , mu = mu , broadening = broadening , mesh = mesh , with_Sigma = with_Sigma , with_dc = with_dc )
# Non-projected A(k,w)
for bname , gf in G_latt_w :
Akw [ bname ] [ ik ] = - gf . data [ numpy . where ( ( mesh_val > om_minplot ) &
( mesh_val < om_maxplot ) ) ] . imag . trace ( axis1 = 1 , axis2 = 2 ) / numpy . pi
# shift Akw for plotting stacked k-resolved eps(k) curves
Akw [ bname ] [ ik ] + = ik * plot_shift
#project spectral functions
if ( proj_type != None ) :
# Projected A(k,w):
for ish in range ( n_shells ) :
G_loc [ ish ] . zero ( )
tmp = G_loc [ ish ] . copy ( )
tmp . zero ( )
tmp << self . proj_type_G_loc ( G_latt_w , tmp , ik , ish , proj_type )
G_loc [ ish ] + = tmp
# Rotate to local frame
if ( self . use_rotations ) :
for ish in range ( n_shells ) :
jsh = shell_list [ ish ]
for bname , gf in G_loc [ ish ] :
G_loc [ ish ] [ bname ] << self . rotloc (
jsh , gf , direction = ' toLocal ' , shells = shells_type )
for ish in range ( n_shells ) :
for bname , gf in G_loc [ ish ] : # loop over spins
pAkw_orb [ ish ] [ bname ] [ ik , : , : , : ] = - gf . data [ numpy . where ( ( mesh_val > om_minplot ) &
( mesh_val < om_maxplot ) ) , : , : ] . imag / numpy . pi
# shift pAkw_orb for plotting stacked k-resolved eps(k) curves
pAkw_orb [ ish ] [ sp ] [ ik ] + = ik * plot_shift
# Collect data from mpi
mpi . barrier ( )
for sp in spn :
2023-06-05 15:53:08 +02:00
Akw [ sp ] = mpi . all_reduce ( Akw [ sp ] )
2023-04-15 00:43:23 +02:00
if ( proj_type ) :
for ish in range ( n_shells ) :
2023-06-05 15:53:08 +02:00
pAkw_orb [ ish ] [ sp ] = mpi . all_reduce ( pAkw_orb [ ish ] [ sp ] )
2023-04-15 00:43:23 +02:00
pAkw [ ish ] [ sp ] = pAkw_orb [ ish ] [ sp ] . trace ( axis1 = 2 , axis2 = 3 )
mpi . barrier ( )
2013-07-23 19:49:42 +02:00
2023-04-15 00:43:23 +02:00
return Akw , pAkw , pAkw_orb
2021-08-27 23:31:07 +02:00
def partial_charges ( self , mu = None , with_Sigma = True , with_dc = True ) :
2015-03-12 00:01:12 +01:00
"""
Calculates the orbitally - resolved density matrix for all the orbitals considered in the input , consistent with
the definition of Wien2k . Hence , ( possibly non - orthonormal ) projectors have to be provided in the partial projectors subgroup of
the hdf5 archive .
2016-05-09 10:19:56 +02:00
2015-03-12 00:01:12 +01:00
Parameters
- - - - - - - - - -
2016-05-09 10:19:56 +02:00
2015-03-12 00:01:12 +01:00
with_Sigma : boolean , optional
If True , the self energy is used for the calculation . If false , partial charges are calculated without self - energy correction .
mu : double , optional
Chemical potential , overrides the one stored in the hdf5 archive .
with_dc : boolean , optional
If True the double counting correction is used .
2016-05-09 10:19:56 +02:00
2015-03-12 00:01:12 +01:00
Returns
- - - - - - -
dens_mat : list of numpy array
A list of density matrices projected to all shells provided in the input .
"""
2023-01-04 23:16:57 +01:00
assert self . dft_code in ( ' wien2k ' ) , " This routine has only been implemented for wien2k inputs "
2016-05-09 10:19:56 +02:00
things_to_read = [ ' dens_mat_below ' , ' n_parproj ' ,
' proj_mat_all ' , ' rot_mat_all ' , ' rot_mat_all_time_inv ' ]
2023-03-23 22:40:13 +01:00
subgroup_present , values_not_read = self . read_input_from_hdf (
2016-05-09 10:19:56 +02:00
subgrp = self . parproj_data , things_to_read = things_to_read )
2023-03-23 22:40:13 +01:00
if len ( values_not_read ) > 0 and mpi . is_master_node :
raise ValueError (
' ERROR: One or more necessary SumK input properties have not been found in the given h5 archive: ' , self . values_not_read )
2016-05-09 10:19:56 +02:00
if self . symm_op :
self . symmpar = Symmetry ( self . hdf_file , subgroup = self . symmpar_data )
2014-09-22 19:24:33 +02:00
2014-11-17 10:48:04 +01:00
spn = self . spin_block_names [ self . SO ]
2014-11-14 09:43:28 +01:00
ntoi = self . spin_names_to_ind [ self . SO ]
2015-03-13 11:10:38 +01:00
# Density matrix in the window
2023-01-23 21:40:57 +01:00
self . dens_mat_window = [ [ numpy . zeros ( [ self . shells [ ish ] [ ' dim ' ] , self . shells [ ish ] [ ' dim ' ] ] , complex )
2016-05-09 10:19:56 +02:00
for ish in range ( self . n_shells ) ]
for isp in range ( len ( spn ) ) ]
2015-03-13 11:10:38 +01:00
# Set up G_loc
2022-03-08 18:49:12 +01:00
gf_struct_parproj = [ [ ( sp , self . shells [ ish ] [ ' dim ' ] ) for sp in spn ]
2016-05-09 10:19:56 +02:00
for ish in range ( self . n_shells ) ]
2021-08-27 23:31:07 +02:00
G_loc = [ BlockGf ( name_block_generator = [ ( block , GfImFreq ( target_shape = ( block_dim , block_dim ) , mesh = self . mesh ) )
for block , block_dim in gf_struct_parproj [ ish ] ] , make_copies = False )
for ish in range ( self . n_shells ) ]
2016-05-09 10:19:56 +02:00
for ish in range ( self . n_shells ) :
G_loc [ ish ] . zero ( )
2013-07-23 19:49:42 +02:00
2020-04-08 21:35:59 +02:00
ikarray = numpy . array ( list ( range ( self . n_k ) ) )
2013-07-23 19:49:42 +02:00
for ik in mpi . slice_array ( ikarray ) :
2021-08-27 23:31:07 +02:00
G_latt_iw = self . lattice_gf ( ik = ik , mu = mu , with_Sigma = with_Sigma , with_dc = with_dc )
2015-03-13 11:10:38 +01:00
G_latt_iw * = self . bz_weights [ ik ]
2014-11-15 20:04:54 +01:00
for ish in range ( self . n_shells ) :
2015-03-13 11:10:38 +01:00
tmp = G_loc [ ish ] . copy ( )
2014-11-15 20:04:54 +01:00
for ir in range ( self . n_parproj [ ish ] ) :
2016-05-09 10:19:56 +02:00
for bname , gf in tmp :
tmp [ bname ] << self . downfold ( ik , ish , bname , G_latt_iw [
bname ] , gf , shells = ' all ' , ir = ir )
2015-03-13 11:10:38 +01:00
G_loc [ ish ] + = tmp
2014-09-22 19:24:33 +02:00
2015-03-13 11:10:38 +01:00
# Collect data from mpi:
2014-11-15 20:04:54 +01:00
for ish in range ( self . n_shells ) :
2023-06-02 16:51:48 +02:00
G_loc [ ish ] << mpi . all_reduce ( G_loc [ ish ] )
2013-07-23 19:49:42 +02:00
mpi . barrier ( )
2015-03-13 11:10:38 +01:00
# Symmetrize and rotate to local coord. system if needed:
2016-05-09 10:19:56 +02:00
if self . symm_op != 0 :
G_loc = self . symmpar . symmetrize ( G_loc )
2015-03-13 11:10:38 +01:00
if self . use_rotations :
for ish in range ( self . n_shells ) :
2016-05-09 10:19:56 +02:00
for bname , gf in G_loc [ ish ] :
G_loc [ ish ] [ bname ] << self . rotloc (
ish , gf , direction = ' toLocal ' , shells = ' all ' )
2014-09-22 19:24:33 +02:00
2014-11-15 20:04:54 +01:00
for ish in range ( self . n_shells ) :
2013-07-23 19:49:42 +02:00
isp = 0
2016-05-09 10:19:56 +02:00
for bname , gf in G_loc [ ish ] :
2015-03-13 11:10:38 +01:00
self . dens_mat_window [ isp ] [ ish ] = G_loc [ ish ] . density ( ) [ bname ]
isp + = 1
2014-09-22 19:24:33 +02:00
2015-03-13 11:10:38 +01:00
# Add density matrices to get the total:
2016-05-09 10:19:56 +02:00
dens_mat = [ [ self . dens_mat_below [ ntoi [ spn [ isp ] ] ] [ ish ] + self . dens_mat_window [ isp ] [ ish ]
for ish in range ( self . n_shells ) ]
for isp in range ( len ( spn ) ) ]
2013-07-23 19:49:42 +02:00
return dens_mat
2014-11-20 13:17:46 +01:00
2015-02-10 16:35:12 +01:00
def print_hamiltonian ( self ) :
2015-03-12 00:01:12 +01:00
"""
Prints the Kohn - Sham Hamiltonian to the text files hamup . dat and hamdn . dat ( no spin orbit - coupling ) , or to ham . dat ( with spin - orbit coupling ) .
"""
2015-02-10 16:35:12 +01:00
if self . SP == 1 and self . SO == 0 :
2016-05-09 10:19:56 +02:00
f1 = open ( ' hamup.dat ' , ' w ' )
f2 = open ( ' hamdn.dat ' , ' w ' )
2015-02-10 16:35:12 +01:00
for ik in range ( self . n_k ) :
2016-05-09 10:19:56 +02:00
for i in range ( self . n_orbitals [ ik , 0 ] ) :
f1 . write ( ' %s %s \n ' %
( ik , self . hopping [ ik , 0 , i , i ] . real ) )
for i in range ( self . n_orbitals [ ik , 1 ] ) :
f2 . write ( ' %s %s \n ' %
( ik , self . hopping [ ik , 1 , i , i ] . real ) )
2015-02-10 16:35:12 +01:00
f1 . write ( ' \n ' )
f2 . write ( ' \n ' )
f1 . close ( )
f2 . close ( )
else :
2016-05-09 10:19:56 +02:00
f = open ( ' ham.dat ' , ' w ' )
2015-02-10 16:35:12 +01:00
for ik in range ( self . n_k ) :
2016-05-09 10:19:56 +02:00
for i in range ( self . n_orbitals [ ik , 0 ] ) :
f . write ( ' %s %s \n ' %
( ik , self . hopping [ ik , 0 , i , i ] . real ) )
2015-02-10 16:35:12 +01:00
f . write ( ' \n ' )
f . close ( )