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

Merge pull request #185 from jkarp314/vectorize_omega

Browse Source 
Vectorize omega when calculating spectral function
This commit is contained in:
Alexander Hampel 2021-10-14 08:12:30 -04:00 committed by GitHub
commit 8fb173ef5a
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 52 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
python/triqs_dft_tools/sumk_dft_tools.py
View File

@ -153,9 +153,7 @@ class SumkDFTTools(SumkDFT):
# G_loc can now also be used to look at orbitally-resolved quantities
for ish in range(self.n_inequiv_shells):
for bname, gf in G_loc[self.inequiv_to_corr[ish]]: # loop over spins
for iom in range(n_om):
DOSproj[ish][bname][iom] -= gf.data[iom,
:, :].imag.trace() / numpy.pi
DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
DOSproj_orb[ish][bname][
:, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
@ -279,8 +277,7 @@ class SumkDFTTools(SumkDFT):
# G_loc can now also be used to look at orbitally-resolved quantities
for bname, gf in G_loc_all: # loop over spins
for iom in range(n_om):
DOSproj[bname][iom] -= gf.data[iom,:,:].imag.trace() / numpy.pi
DOSproj[bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
DOSproj_orb[bname][:,:,:] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
# Write to files
if save_to_file and mpi.is_master_node():
@ -392,10 +389,8 @@ class SumkDFTTools(SumkDFT):
G_latt_w *= self.bz_weights[ik]
# Non-projected DOS
for iom in range(n_om):
for bname, gf in G_latt_w:
DOS[bname][iom] -= gf.data[iom, :, :].imag.trace() / \
numpy.pi
for bname, gf in G_latt_w:
DOS[bname] -= gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
# Projected DOS:
for ish in range(self.n_shells):
@ -427,9 +422,7 @@ class SumkDFTTools(SumkDFT):
# G_loc can now also be used to look at orbitally-resolved quantities
for ish in range(self.n_shells):
for bname, gf in G_loc[ish]:
for iom in range(n_om):
DOSproj[ish][bname][iom] -= gf.data[iom,
:, :].imag.trace() / numpy.pi
DOSproj[ish][bname] = -gf.data.imag.trace(axis1=1, axis2=2) / numpy.pi
DOSproj_orb[ish][bname][
:, :, :] += (1.0j*(gf-gf.conjugate().transpose())/2.0/numpy.pi).data[:,:,:]
@ -857,15 +850,15 @@ class SumkDFTTools(SumkDFT):
if mu is None:
mu = self.chemical_potential
spn = self.spin_block_names[self.SO]
mesh = [x.real for x in self.Sigma_imp_w[0].mesh]
n_om = len(mesh)
mesh = numpy.array([x.real for x in self.Sigma_imp_w[0].mesh])
if plot_range is None:
om_minplot = mesh[0] - 0.001
om_maxplot = mesh[n_om - 1] + 0.001
om_maxplot = mesh[-1] + 0.001
else:
om_minplot = plot_range[0]
om_maxplot = plot_range[1]
n_om = len(mesh[(mesh > om_minplot)&(mesh < om_maxplot)])
if ishell is None:
Akw = {sp: numpy.zeros([self.n_k, n_om], numpy.float_)
@ -889,14 +882,11 @@ class SumkDFTTools(SumkDFT):
if ishell is None:
# Non-projected A(k,w)
for iom in range(n_om):
if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot):
for bname, gf in G_latt_w:
Akw[bname][ik, iom] += gf.data[iom, :,
:].imag.trace() / (-1.0 * numpy.pi)
# shift Akw for plotting stacked k-resolved eps(k)
# curves
Akw[bname][ik, iom] += ik * plot_shift
for bname, gf in G_latt_w:
Akw[bname][ik] = -gf.data[numpy.where((mesh > om_minplot)&(mesh < 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
else: # ishell not None
# Projected A(k,w):
@ -914,13 +904,9 @@ class SumkDFTTools(SumkDFT):
G_loc[bname] << self.rotloc(
ishell, gf, direction='toLocal', shells='all')
for iom in range(n_om):
if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot):
for ish in range(self.shells[ishell]['dim']):
for sp in spn:
Akw[sp][ish, ik, iom] = G_loc[sp].data[
iom, ish, ish].imag / (-1.0 * numpy.pi)
for ish in range(self.shells[ishell]['dim']):
for sp in spn:
Akw[sp][ish, ik] = -gf.data[numpy.where((mesh > om_minplot)&(mesh < om_maxplot))].imag.trace(axis1=1, axis2=2)/numpy.pi
# Collect data from mpi
for sp in spn:
Akw[sp] = mpi.all_reduce(mpi.world, Akw[sp], lambda x, y: x + y)

BIN
test/python/SrVO3_spectral.h5 Normal file
View File

Binary file not shown.

36
test/python/srvo3_spectral.py Normal file
View File

@ -0,0 +1,36 @@
from triqs_dft_tools.sumk_dft_tools import *
from triqs.utility.h5diff import h5diff
from h5 import HDFArchive
beta = 40
SK = SumkDFTTools(hdf_file='SrVO3_spectral.h5', use_dft_blocks=True)
if mpi.is_master_node():
with HDFArchive('SrVO3_Sigma.h5', 'a') as ar:
Sigma = ar['dmft_transp_input']['Sigma_w']
SK.chemical_potential = ar['dmft_transp_input']['chemical_potential']
SK.dc_imp = ar['dmft_transp_input']['dc_imp']
Sigma = mpi.bcast(Sigma)
SK.chemical_potential = mpi.bcast(SK.chemical_potential)
SK.dc_imp = mpi.bcast(SK.dc_imp)
SK.set_Sigma([Sigma])
dos_wannier = SK.dos_wannier_basis(broadening=0.01, with_Sigma=True, with_dc=True, save_to_file=False)
dos_parproj = SK.dos_parproj_basis(broadening=0.01, with_Sigma=True, with_dc=True, save_to_file=False)
spaghetti = SK.spaghettis(broadening=0.01, plot_shift=0.0, plot_range=(-1,1), ishell=None, save_to_file=False)
if mpi.is_master_node():
# with HDFArchive('srvo3_spectral.ref.h5', 'a') as ar:
# ar['dos_wannier'] = dos_wannier
# ar['dos_parproj'] = dos_parproj
# ar['spaghetti'] = spaghetti
with HDFArchive('srvo3_spectral.out.h5', 'a') as ar:
ar['dos_wannier'] = dos_wannier
ar['dos_parproj'] = dos_parproj
ar['spaghetti'] = spaghetti
h5diff('srvo3_spectral.out.h5', 'srvo3_spectral.ref.h5')

BIN
test/python/srvo3_spectral.ref.h5 Normal file
View File

Binary file not shown.