LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-08-16 01:08:29 +02:00
Code Issues Releases Wiki Activity

[style] format and doc strings

Browse Source
This commit is contained in:
Alexander Hampel 2023-06-07 09:53:19 -04:00
parent 6df646ea0b
commit 5ae4949313
1 changed files with 155 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

234
python/triqs_dft_tools/sumk_dft_tools.py
View File

@ -53,12 +53,11 @@ class SumkDFTTools(SumkDFT):
parproj_data=parproj_data, symmpar_data=symmpar_data, bands_data=bands_data, parproj_data=parproj_data, symmpar_data=symmpar_data, bands_data=bands_data,
transp_data=transp_data, misc_data=misc_data, cont_data=cont_data) transp_data=transp_data, misc_data=misc_data, cont_data=cont_data)
# Uses .data of only GfReFreq objects.
# Uses .data of only GfReFreq objects.
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): 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):
""" """
Calculates the density of states. The basis of the projected density of states is Calculates the density of states and the projected density of states.
specified by proj_type. The basis of the projected density of states is specified by proj_type.
Parameters Parameters
---------- ----------
mu : double, optional mu : double, optional
@ -96,33 +95,35 @@ class SumkDFTTools(SumkDFT):
DOS projected to atoms and resolved into orbital contributions. DOS projected to atoms and resolved into orbital contributions.
Empty if proj_type = None Empty if proj_type = None
""" """
if(proj_type!=None): 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','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'" assert proj_type in ('wann', 'vasp', 'wien2k',
if(proj_type!='wann'): ), "'proj_type' must be either 'wann', 'vasp', 'wien2k'"
assert proj_type==self.dft_code, "proj_type must be from the corresponding dft inputs." if (proj_type != 'wann'):
assert proj_type == self.dft_code, "proj_type must be from the corresponding dft inputs."
if (with_Sigma): if (with_Sigma):
assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True" assert isinstance(
mesh=self.Sigma_imp[0].mesh self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
mesh = self.Sigma_imp[0].mesh
elif mesh is not None: elif mesh is not None:
assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq" assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
elif self.mesh is not None: elif self.mesh is not None:
assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq" assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
mesh=self.mesh mesh = self.mesh
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
else: else:
assert 0, "ReFreqMesh input required for calculations without real frequency self-energy" assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh)) mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
n_om = len(mesh) n_om = len(mesh)
om_minplot = mesh_val[0] - 0.001 om_minplot = mesh_val[0] - 0.001
om_maxplot = mesh_val[-1] + 0.001 om_maxplot = mesh_val[-1] + 0.001
#Read in occupations from HDF5 file if required # Read in occupations from HDF5 file if required
if(dosocc): if (dosocc):
mpi.report('Reading occupations generated by self.occupations().') mpi.report('Reading occupations generated by self.occupations().')
thingstoread = ['occik'] thingstoread = ['occik']
subgroup_present, values_not_read = self.read_input_from_hdf( subgroup_present, values_not_read = self.read_input_from_hdf(
@ -131,16 +132,16 @@ class SumkDFTTools(SumkDFT):
raise ValueError( raise ValueError(
'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read) 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
#initialise projected DOS type if required # initialise projected DOS type if required
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
n_shells=1 n_shells = 1
if (proj_type == 'wann'): if (proj_type == 'wann'):
n_shells = self.n_corr_shells n_shells = self.n_corr_shells
gf_struct = self.gf_struct_sumk.copy() gf_struct = self.gf_struct_sumk.copy()
dims = [self.corr_shells[ish]['dim'] for ish in range(n_shells)] dims = [self.corr_shells[ish]['dim'] for ish in range(n_shells)]
shells_type = 'corr' shells_type = 'corr'
elif (proj_type == 'vasp'): elif (proj_type == 'vasp'):
n_shells=1 n_shells = 1
gf_struct = [[(sp, list(range(self.proj_mat_csc.shape[2]))) for sp in spn]] gf_struct = [[(sp, list(range(self.proj_mat_csc.shape[2]))) for sp in spn]]
dims = [self.proj_mat_csc.shape[2]] dims = [self.proj_mat_csc.shape[2]]
shells_type = 'csc' shells_type = 'csc'
@ -164,11 +165,11 @@ class SumkDFTTools(SumkDFT):
# raise ValueError( # raise ValueError(
# 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read) # 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
#set-up output arrays # set-up output arrays
DOS = {sp: numpy.zeros([n_om],float) for sp in spn} DOS = {sp: numpy.zeros([n_om], float) for sp in spn}
DOSproj = [{} for ish in range(n_shells)] DOSproj = [{} for ish in range(n_shells)]
DOSproj_orb = [{} for ish in range(n_shells)] DOSproj_orb = [{} for ish in range(n_shells)]
#set-up Green's function object # set-up Green's function object
if (proj_type != None): if (proj_type != None):
G_loc = [] G_loc = []
for ish in range(n_shells): for ish in range(n_shells):
@ -183,16 +184,16 @@ class SumkDFTTools(SumkDFT):
DOSproj_orb[ish][sp] = numpy.zeros( DOSproj_orb[ish][sp] = numpy.zeros(
[n_om, dim, dim], complex) [n_om, dim, dim], complex)
#calculate the DOS # calculate the DOS
ikarray = numpy.array(list(range(self.n_k))) ikarray = numpy.array(list(range(self.n_k)))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
G_latt_w = self.lattice_gf( G_latt_w = self.lattice_gf(
ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc) ik=ik, mu=mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc)
G_latt_w *= self.bz_weights[ik] G_latt_w *= self.bz_weights[ik]
#output occupied DOS if nk inputted # output occupied DOS if nk inputted
if(dosocc): if (dosocc):
for bname, gf in G_latt_w: for bname, gf in G_latt_w:
G_latt_w[bname].data[:,:,:] *= self.occik[bname][ik] G_latt_w[bname].data[:, :, :] *= self.occik[bname][ik]
# DOS # DOS
for bname, gf in G_latt_w: for bname, gf in G_latt_w:
DOS[bname] -= gf.data.imag.trace(axis1=1, axis2=2)/numpy.pi DOS[bname] -= gf.data.imag.trace(axis1=1, axis2=2)/numpy.pi
@ -209,13 +210,13 @@ class SumkDFTTools(SumkDFT):
for bname in DOS: for bname in DOS:
DOS[bname] = mpi.all_reduce(DOS[bname]) DOS[bname] = mpi.all_reduce(DOS[bname])
# Collect data from mpi and put in projected arrays # Collect data from mpi and put in projected arrays
if(proj_type != None): if (proj_type != None):
for ish in range(n_shells): for ish in range(n_shells):
G_loc[ish] << mpi.all_reduce(G_loc[ish]) G_loc[ish] << mpi.all_reduce(G_loc[ish])
# Symmetrize and rotate to local coord. system if needed: # Symmetrize and rotate to local coord. system if needed:
if((proj_type!='vasp') and (proj_type!='elk')): if ((proj_type != 'vasp') and (proj_type != 'elk')):
if self.symm_op != 0: if self.symm_op != 0:
if proj_type=='wann': if proj_type == 'wann':
G_loc = self.symmcorr.symmetrize(G_loc) G_loc = self.symmcorr.symmetrize(G_loc)
else: else:
G_loc = self.symmpar.symmetrize(G_loc) G_loc = self.symmpar.symmetrize(G_loc)
@ -223,13 +224,13 @@ class SumkDFTTools(SumkDFT):
for ish in range(n_shells): for ish in range(n_shells):
for bname, gf in G_loc[ish]: for bname, gf in G_loc[ish]:
G_loc[ish][bname] << self.rotloc( G_loc[ish][bname] << self.rotloc(
ish, gf, direction='toLocal',shells=shells_type) ish, gf, direction='toLocal', shells=shells_type)
# G_loc can now also be used to look at orbitally-resolved quantities # G_loc can now also be used to look at orbitally-resolved quantities
for ish in range(n_shells): for ish in range(n_shells):
for bname, gf in G_loc[ish]: # loop over spins for bname, gf in G_loc[ish]: # loop over spins
DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
DOSproj_orb[ish][bname][ DOSproj_orb[ish][bname][
:, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:] :, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:, :, :]
# Write to files # Write to files
if save_to_file and mpi.is_master_node(): if save_to_file and mpi.is_master_node():
@ -239,7 +240,7 @@ class SumkDFTTools(SumkDFT):
f.write("%s %s\n" % (mesh_val[iom], DOS[sp][iom])) f.write("%s %s\n" % (mesh_val[iom], DOS[sp][iom]))
f.close() f.close()
# Partial # Partial
if(proj_type!=None): if (proj_type != None):
for ish in range(n_shells): for ish in range(n_shells):
f = open('DOS_' + proj_type + '_%s_proj%s.dat' % (sp, ish), 'w') f = open('DOS_' + proj_type + '_%s_proj%s.dat' % (sp, ish), 'w')
for iom in range(n_om): for iom in range(n_om):
@ -249,18 +250,17 @@ class SumkDFTTools(SumkDFT):
# Orbitally-resolved # Orbitally-resolved
for i in range(dims[ish]): for i in range(dims[ish]):
for j in range(dims[ish]): for j in range(dims[ish]):
#For Elk with parproj - skip off-diagonal elements # For Elk with parproj - skip off-diagonal elements
#if(proj_type=='elk') and (i!=j): continue # if(proj_type=='elk') and (i!=j): continue
f = open('DOS_' + proj_type + '_' + sp + '_proj' + str(ish) + f = open('DOS_' + proj_type + '_' + sp + '_proj' + str(ish) +
'_' + str(i) + '_' + str(j) + '.dat', 'w') '_' + str(i) + '_' + str(j) + '.dat', 'w')
for iom in range(n_om): for iom in range(n_om):
f.write("%s %s %s\n" % ( f.write("%s %s %s\n" % (
mesh_val[iom], DOSproj_orb[ish][sp][iom, i, j].real,DOSproj_orb[ish][sp][iom, i, j].imag)) mesh_val[iom], DOSproj_orb[ish][sp][iom, i, j].real, DOSproj_orb[ish][sp][iom, i, j].imag))
f.close() f.close()
return DOS, DOSproj, DOSproj_orb return DOS, DOSproj, DOSproj_orb
def proj_type_G_loc(self, G_latt, G_inp, ik, ish, proj_type=None): def proj_type_G_loc(self, G_latt, G_inp, ik, ish, proj_type=None):
""" """
Internal routine which calculates the project Green's function subject to the Internal routine which calculates the project Green's function subject to the
@ -328,7 +328,7 @@ class SumkDFTTools(SumkDFT):
return G_proj return G_proj
def load_parproj(self,data_type=None): def load_parproj(self, data_type=None):
""" """
Internal routine which loads the n_parproj, proj_mat_all, rot_mat_all and 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. rot_mat_all_time_inv from parproj data from .h5 file.
@ -339,7 +339,7 @@ class SumkDFTTools(SumkDFT):
'band' - reads data converted by bands_convert() 'band' - reads data converted by bands_convert()
None - reads data converted by parproj_convert() None - reads data converted by parproj_convert()
""" """
#read in the projectors # read in the projectors
things_to_read = ['n_parproj', 'proj_mat_all'] things_to_read = ['n_parproj', 'proj_mat_all']
if data_type == 'band': if data_type == 'band':
subgroup_present, values_not_read = self.read_input_from_hdf( subgroup_present, values_not_read = self.read_input_from_hdf(
@ -352,7 +352,7 @@ class SumkDFTTools(SumkDFT):
if len(values_not_read) > 0 and mpi.is_master_node: if len(values_not_read) > 0 and mpi.is_master_node:
raise ValueError( raise ValueError(
'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read) 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
#read general data # read general data
things_to_read = ['rot_mat_all', 'rot_mat_all_time_inv'] things_to_read = ['rot_mat_all', 'rot_mat_all_time_inv']
subgroup_present, values_not_read = self.read_input_from_hdf( subgroup_present, values_not_read = self.read_input_from_hdf(
subgrp=self.parproj_data, things_to_read=things_to_read) subgrp=self.parproj_data, things_to_read=things_to_read)
@ -364,6 +364,7 @@ class SumkDFTTools(SumkDFT):
""" """
Calculates the band resolved density matrices (occupations) from the Matsubara Calculates the band resolved density matrices (occupations) from the Matsubara
frequency self-energy. frequency self-energy.
Parameters Parameters
---------- ----------
mu : double, optional mu : double, optional
@ -376,7 +377,8 @@ class SumkDFTTools(SumkDFT):
save_occ : boolean, optional save_occ : boolean, optional
If True, saves the band resolved density matrix in misc_data. If True, saves the band resolved density matrix in misc_data.
save_to_file : boolean, optional save_to_file : boolean, optional
If True, text files with the calculated data will be created. If True, text files with the calculated data will be created.\
Returns Returns
------- -------
occik : Dict of numpy arrays occik : Dict of numpy arrays
@ -387,7 +389,8 @@ class SumkDFTTools(SumkDFT):
mesh = self.Sigma_imp[0].mesh mesh = self.Sigma_imp[0].mesh
else: else:
mesh = self.mesh mesh = self.mesh
assert isinstance(mesh, MeshImFreq), "SumkDFT.mesh must be real if with_Sigma is True or mesh is not given" assert isinstance(
mesh, MeshImFreq), "SumkDFT.mesh must be real if with_Sigma is True or mesh is not given"
if mu is None: if mu is None:
mu = self.chemical_potential mu = self.chemical_potential
@ -395,20 +398,21 @@ class SumkDFTTools(SumkDFT):
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
occik = {} occik = {}
for sp in spn: for sp in spn:
#same format as gf.data ndarray # 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)] occik[sp] = [numpy.zeros([1, self.n_orbitals[ik, ntoi[sp]],
#calculate the occupations self.n_orbitals[ik, ntoi[sp]]], numpy.double) for ik in range(self.n_k)]
# calculate the occupations
ikarray = numpy.array(range(self.n_k)) ikarray = numpy.array(range(self.n_k))
for ik in range(self.n_k): for ik in range(self.n_k):
G_latt = self.lattice_gf( G_latt = self.lattice_gf(
ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc) ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc)
for bname, gf in G_latt: for bname, gf in G_latt:
occik[bname][ik][0,:,:] = gf.density().real occik[bname][ik][0, :, :] = gf.density().real
# Collect data from mpi: # Collect data from mpi:
for sp in spn: for sp in spn:
occik[sp] = mpi.all_reduce(occik[sp]) occik[sp] = mpi.all_reduce(occik[sp])
mpi.barrier() mpi.barrier()
#save to HDF5 file (if specified) # save to HDF5 file (if specified)
if save_occ and mpi.is_master_node(): if save_occ and mpi.is_master_node():
things_to_save_misc = ['occik'] things_to_save_misc = ['occik']
# Save it to the HDF: # Save it to the HDF:
@ -420,7 +424,6 @@ class SumkDFTTools(SumkDFT):
del ar del ar
return occik return occik
# Uses .data of only GfReFreq objects.
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): 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):
""" """
Calculates the correlated spectral function at the Fermi level (relating to the Fermi Calculates the correlated spectral function at the Fermi level (relating to the Fermi
@ -465,10 +468,10 @@ class SumkDFTTools(SumkDFT):
(Correlated) k-resolved spectral function projected to atoms and (Correlated) k-resolved spectral function projected to atoms and
resolved into orbital contributions. Empty if proj_type = None resolved into orbital contributions. Empty if proj_type = None
""" """
if(proj_type!=None): if (proj_type != None):
assert proj_type in ('wann'), "'proj_type' must be 'wann' if not None" assert proj_type in ('wann'), "'proj_type' must be 'wann' if not None"
#read in the energy contour energies and projectors # read in the energy contour energies and projectors
things_to_read = ['n_k','bmat','BZ_n_k','BZ_iknr','BZ_vkl', things_to_read = ['n_k', 'bmat', 'BZ_n_k', 'BZ_iknr', 'BZ_vkl',
'n_orbitals', 'proj_mat', 'hopping'] 'n_orbitals', 'proj_mat', 'hopping']
subgroup_present, values_not_read = self.read_input_from_hdf( subgroup_present, values_not_read = self.read_input_from_hdf(
subgrp=self.cont_data, things_to_read=things_to_read) subgrp=self.cont_data, things_to_read=things_to_read)
@ -479,51 +482,53 @@ class SumkDFTTools(SumkDFT):
if mu is None: if mu is None:
mu = self.chemical_potential mu = self.chemical_potential
if (with_Sigma): if (with_Sigma):
assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True" assert isinstance(
mesh=self.Sigma_imp[0].mesh self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
mesh = self.Sigma_imp[0].mesh
elif mesh is not None: elif mesh is not None:
assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq" assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
elif self.mesh is not None: elif self.mesh is not None:
assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq" assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
mesh=self.mesh mesh = self.mesh
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
else: else:
assert 0, "ReFreqMesh input required for calculations without real frequency self-energy" assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh)) mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
n_om = len(mesh) n_om = len(mesh)
om_minplot = mesh_val[0] - 0.001 om_minplot = mesh_val[0] - 0.001
om_maxplot = mesh_val[-1] + 0.001 om_maxplot = mesh_val[-1] + 0.001
#for Fermi Surface calculations # for Fermi Surface calculations
if FS: if FS:
dw = abs(mesh_val[1]-mesh_val[0]) dw = abs(mesh_val[1]-mesh_val[0])
#ensure that a few frequencies around the Fermi level are included # ensure that a few frequencies around the Fermi level are included
plot_range = [-2*dw, 2*dw] plot_range = [-2*dw, 2*dw]
mpi.report('Generated A(k,w) will be evaluted at closest frequency to 0.0 in given mesh ') mpi.report('Generated A(k,w) will be evaluted at closest frequency to 0.0 in given mesh ')
if plot_range is None: if plot_range is None:
n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]) 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)] mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
else: else:
om_minplot = plot_range[0] om_minplot = plot_range[0]
om_maxplot = plot_range[1] om_maxplot = plot_range[1]
n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]) 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)] mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
#\omega ~= 0.0 index for FS file # \omega ~= 0.0 index for FS file
abs_mesh_val = [abs(i) for i in mesh_val2] abs_mesh_val = [abs(i) for i in mesh_val2]
jw=[i for i in range(len(abs_mesh_val)) if abs_mesh_val[i] == numpy.min(abs_mesh_val[:])] 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 # calculate the spectral functions for the irreducible set of k-points
[Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh, \ [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh,
plot_shift=0.0, plot_range=plot_range, \ plot_shift=0.0, plot_range=plot_range,
shell_list=None, with_Sigma=with_Sigma, with_dc=with_dc, \ shell_list=None, with_Sigma=with_Sigma, with_dc=with_dc,
proj_type=proj_type) proj_type=proj_type)
if save_to_file and mpi.is_master_node(): if save_to_file and mpi.is_master_node():
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
vkc = numpy.zeros(3, float) vkc = numpy.zeros(3, float)
mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)] mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
if FS: if FS:
n_om = 1 n_om = 1
else: else:
@ -533,53 +538,54 @@ class SumkDFTTools(SumkDFT):
for iom in range(n_om): for iom in range(n_om):
if FS: if FS:
f = open('Akw_FS_' + sp + '.dat', 'w') f = open('Akw_FS_' + sp + '.dat', 'w')
jom=jw[0] jom = jw[0]
else: else:
f = open('Akw_omega_%s_%s.dat' % (iom, sp), 'w') f = open('Akw_omega_%s_%s.dat' % (iom, sp), 'w')
jom=iom jom = iom
f.write("#Spectral function evaluated at frequency = %s\n" %mesh_val2[jom]) f.write("#Spectral function evaluated at frequency = %s\n" % mesh_val2[jom])
for ik in range(self.BZ_n_k): for ik in range(self.BZ_n_k):
jk=self.BZ_iknr[ik] jk = self.BZ_iknr[ik]
vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:]) vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2], Akw[sp][jk, jom])) f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2], Akw[sp][jk, jom]))
f.close() f.close()
if (proj_type!=None): if (proj_type != None):
n_shells = len(pAkw[:]) n_shells = len(pAkw[:])
for iom in range(n_om): for iom in range(n_om):
for sp in spn: for sp in spn:
for ish in range(n_shells): for ish in range(n_shells):
if FS: if FS:
strng = 'Akw_FS' + '_' + proj_type + '_' + sp + '_proj' + str(ish) strng = 'Akw_FS' + '_' + proj_type + '_' + sp + '_proj' + str(ish)
jom=jw[0] jom = jw[0]
else: else:
strng = 'Akw_omega_' + str(iom) + '_' + proj_type + '_' + sp + '_proj' + str(ish) strng = 'Akw_omega_' + str(iom) + '_' + proj_type + \
jom=iom '_' + sp + '_proj' + str(ish)
jom = iom
f = open(strng + '.dat', 'w') f = open(strng + '.dat', 'w')
f.write("#Spectral function evaluated at frequency = %s\n" %mesh_val2[jom]) f.write("#Spectral function evaluated at frequency = %s\n" % mesh_val2[jom])
for ik in range(self.BZ_n_k): for ik in range(self.BZ_n_k):
jk=self.BZ_iknr[ik] jk = self.BZ_iknr[ik]
vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:]) vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2], f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2],
pAkw[ish][sp][jk, jom])) pAkw[ish][sp][jk, jom]))
f.close() f.close()
dim=len(pAkw_orb[ish][sp][0, 0, 0, :]) dim = len(pAkw_orb[ish][sp][0, 0, 0, :])
for i in range(dim): for i in range(dim):
for j in range(dim): for j in range(dim):
#For Elk with parproj - skip off-diagonal elements # For Elk with parproj - skip off-diagonal elements
if(proj_type=='elk') and (i!=j): continue if (proj_type == 'elk') and (i != j):
continue
strng2 = strng + '_' + str(i) + '_' + str(j) strng2 = strng + '_' + str(i) + '_' + str(j)
# Open file for storage: # Open file for storage:
f = open(strng2 + '.dat', 'w') f = open(strng2 + '.dat', 'w')
for ik in range(self.BZ_n_k): for ik in range(self.BZ_n_k):
jk=self.BZ_iknr[ik] jk = self.BZ_iknr[ik]
vkc[:] = numpy.matmul(self.bmat,self.BZ_vkl[ik,:]) vkc[:] = numpy.matmul(self.bmat, self.BZ_vkl[ik, :])
f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2], f.write("%s %s %s %s\n" % (vkc[0], vkc[1], vkc[2],
pAkw_orb[ish][sp][jk, jom, i, j])) pAkw_orb[ish][sp][jk, jom, i, j]))
f.close() f.close()
return Akw, pAkw, pAkw_orb return Akw, pAkw, pAkw_orb
# Uses .data of only GfReFreq objects.
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): 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):
""" """
Calculates the k-resolved spectral function A(k,w) (band structure) Calculates the k-resolved spectral function A(k,w) (band structure)
@ -627,37 +633,38 @@ class SumkDFTTools(SumkDFT):
resolved into orbital contributions. Empty if proj_type = None resolved into orbital contributions. Empty if proj_type = None
""" """
#initialisation # initialisation
if(proj_type!=None): if (proj_type != None):
assert proj_type in ('wann', 'wien2k'), "'proj_type' must be either 'wann', 'wien2k'" assert proj_type in ('wann', 'wien2k'), "'proj_type' must be either 'wann', 'wien2k'"
if(proj_type!='wann'): if (proj_type != 'wann'):
assert proj_type==self.dft_code, "proj_type must be from the corresponding dft inputs." assert proj_type == self.dft_code, "proj_type must be from the corresponding dft inputs."
things_to_read = ['n_k', 'n_orbitals', 'proj_mat', 'hopping'] things_to_read = ['n_k', 'n_orbitals', 'proj_mat', 'hopping']
subgroup_present, values_not_read = self.read_input_from_hdf( subgroup_present, values_not_read = self.read_input_from_hdf(
subgrp=self.bands_data, things_to_read=things_to_read) subgrp=self.bands_data, things_to_read=things_to_read)
if len(values_not_read) > 0 and mpi.is_master_node: if len(values_not_read) > 0 and mpi.is_master_node:
raise ValueError( raise ValueError(
'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read) 'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
if(proj_type=='wien2k'): if (proj_type == 'wien2k'):
self.load_parproj(data_type='band') self.load_parproj(data_type='band')
if mu is None: if mu is None:
mu = self.chemical_potential mu = self.chemical_potential
if (with_Sigma): if (with_Sigma):
assert isinstance(self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True" assert isinstance(
mesh=self.Sigma_imp[0].mesh self.Sigma_imp[0].mesh, MeshReFreq), "SumkDFT.mesh must be real if with_Sigma is True"
mesh = self.Sigma_imp[0].mesh
elif mesh is not None: elif mesh is not None:
assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq" assert isinstance(mesh, MeshReFreq), "mesh must be of form MeshReFreq"
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
elif self.mesh is not None: elif self.mesh is not None:
assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq" assert isinstance(self.mesh, MeshReFreq), "self.mesh must be of form MeshReFreq"
mesh=self.mesh mesh = self.mesh
if broadening is None: if broadening is None:
broadening=0.001 broadening = 0.001
else: else:
assert 0, "ReFreqMesh input required for calculations without real frequency self-energy" assert 0, "ReFreqMesh input required for calculations without real frequency self-energy"
mesh_val = numpy.linspace(mesh.w_min,mesh.w_max,len(mesh)) mesh_val = numpy.linspace(mesh.w_min, mesh.w_max, len(mesh))
n_om = len(mesh) n_om = len(mesh)
om_minplot = mesh_val[0] - 0.001 om_minplot = mesh_val[0] - 0.001
om_maxplot = mesh_val[-1] + 0.001 om_maxplot = mesh_val[-1] + 0.001
@ -667,31 +674,31 @@ class SumkDFTTools(SumkDFT):
else: else:
om_minplot = plot_range[0] om_minplot = plot_range[0]
om_maxplot = plot_range[1] om_maxplot = plot_range[1]
n_om = len(mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)]) n_om = len(mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)])
[Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh, \ [Akw, pAkw, pAkw_orb] = self.gen_Akw(mu=mu, broadening=broadening, mesh=mesh,
plot_shift=plot_shift, plot_range=plot_range, \ plot_shift=plot_shift, plot_range=plot_range,
shell_list=shell_list, with_Sigma=with_Sigma, with_dc=with_dc, \ shell_list=shell_list, with_Sigma=with_Sigma, with_dc=with_dc,
proj_type=proj_type) proj_type=proj_type)
if save_to_file and mpi.is_master_node(): if save_to_file and mpi.is_master_node():
mesh_val2 = mesh_val[(mesh_val > om_minplot)&(mesh_val < om_maxplot)] mesh_val2 = mesh_val[(mesh_val > om_minplot) & (mesh_val < om_maxplot)]
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
for sp in spn: for sp in spn:
# Open file for storage: # Open file for storage:
f = open('Akw_' + sp + '.dat', 'w') f = open('Akw_' + sp + '.dat', 'w')
for ik in range(self.n_k): for ik in range(self.n_k):
for iom in range(n_om): for iom in range(n_om):
f.write('%s %s %s\n' %(ik, mesh_val2[iom], Akw[sp][ik, iom])) f.write('%s %s %s\n' % (ik, mesh_val2[iom], Akw[sp][ik, iom]))
f.write('\n') f.write('\n')
f.close() f.close()
if (proj_type!=None): if (proj_type != None):
n_shells = len(pAkw[:]) n_shells = len(pAkw[:])
if shell_list==None: if shell_list == None:
shell_list=[ish for ish in range(n_shells)] shell_list = [ish for ish in range(n_shells)]
for sp in spn: for sp in spn:
for ish in range(n_shells): for ish in range(n_shells):
jsh=shell_list[ish] jsh = shell_list[ish]
f = open('Akw_' + proj_type + '_' + f = open('Akw_' + proj_type + '_' +
sp + '_proj' + str(jsh) + '.dat', 'w') sp + '_proj' + str(jsh) + '.dat', 'w')
for ik in range(self.n_k): for ik in range(self.n_k):
@ -700,12 +707,13 @@ class SumkDFTTools(SumkDFT):
ik, mesh_val2[iom], pAkw[ish][sp][ik, iom])) ik, mesh_val2[iom], pAkw[ish][sp][ik, iom]))
f.write('\n') f.write('\n')
f.close() f.close()
#get orbital dimension from the length of dimension of the array # get orbital dimension from the length of dimension of the array
dim=len(pAkw_orb[ish][sp][0, 0, 0, :]) dim = len(pAkw_orb[ish][sp][0, 0, 0, :])
for i in range(dim): for i in range(dim):
for j in range(dim): for j in range(dim):
#For Elk with parproj - skip off-diagonal elements # For Elk with parproj - skip off-diagonal elements
if(proj_type=='elk') and (i!=j): continue if(proj_type =='elk') and (i!=j):
continue
# Open file for storage: # Open file for storage:
f = open('Akw_' + proj_type + '_' + sp + '_proj' + str(jsh) f = open('Akw_' + proj_type + '_' + sp + '_proj' + str(jsh)
+ '_' + str(i) + '_' + str(j) + '.dat', 'w') + '_' + str(i) + '_' + str(j) + '.dat', 'w')