2021-08-04 18:56:09 +02:00
##########################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
2023-03-15 20:33:49 +01:00
# Copyright (c) 2022-2023 Simons Foundation
2021-08-04 18:56:09 +02:00
#
# 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/>.
#
2023-03-15 20:33:49 +01:00
# Authors: M. Aichhorn, S. Beck, A. Hampel, L. Pourovskii, V. Vildosola
2021-08-04 18:56:09 +02:00
##########################################################################
import sys
import numpy
2023-01-18 20:27:19 +01:00
from warnings import warn
2021-08-04 18:56:09 +02:00
from triqs . gf import *
import triqs . utility . mpi as mpi
from . symmetry import *
import scipy . constants as cst
import os . path
__all__ = [ ' transport_distribution ' , ' conductivity_and_seebeck ' , ' write_output_to_hdf ' ,
2023-10-04 14:12:51 +02:00
' init_spectroscopy ' , ' transport_function ' , ' raman_vertex ' ]
2021-08-04 18:56:09 +02:00
# ----------------- helper functions -----------------------
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def read_transport_input_from_hdf ( sum_k ) :
r """
Reads the data for transport calculations from the hdf5 archive .
Parameters
- - - - - - - - - -
sum_k : sum_k object
triqs SumkDFT object
Returns
- - - - - - -
sum_k : sum_k object
triqs SumkDFT object
"""
2023-03-15 20:35:39 +01:00
assert sum_k . dft_code in (
' wien2k ' , ' elk ' , ' w90 ' ) , " read_transport_input_from_hdf() is only implemented for wien2k and elk inputs "
2021-08-04 18:56:09 +02:00
2023-03-15 20:35:39 +01:00
if sum_k . dft_code in ( ' wien2k ' , ' elk ' ) :
thingstoread = [ ' band_window_optics ' , ' velocities_k ' ]
else :
thingstoread = [ ' band_window_optics ' ]
2021-08-04 18:56:09 +02:00
2023-03-15 20:35:39 +01:00
sum_k . read_input_from_hdf ( subgrp = sum_k . transp_data , things_to_read = thingstoread )
if ( sum_k . dft_code == " wien2k " ) :
thingstoread = [ ' band_window ' , ' lattice_angles ' , ' lattice_constants ' ,
' lattice_type ' , ' n_symmetries ' , ' rot_symmetries ' ]
elif ( sum_k . dft_code == " elk " ) :
thingstoread = [ ' band_window ' , ' n_symmetries ' , ' rot_symmetries ' ,
' cell_vol ' ]
elif ( sum_k . dft_code == ' w90 ' ) :
thingstoread = [ ' band_window ' , ' n_symmetries ' , ' rot_symmetries ' ]
2021-08-04 18:56:09 +02:00
sum_k . read_input_from_hdf ( subgrp = sum_k . misc_data , things_to_read = thingstoread )
2023-03-15 20:35:39 +01:00
if ( sum_k . dft_code == " wien2k " ) :
sum_k . cell_vol = cellvolume (
sum_k . lattice_type , sum_k . lattice_constants , sum_k . lattice_angles ) [ 1 ]
2021-08-04 18:56:09 +02:00
return sum_k
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def write_output_to_hdf ( sum_k , things_to_save , subgrp = ' user_data ' ) :
r """
Saves data from a list into the HDF file . Prints a warning if a requested data is not found in SumkDFT object .
Parameters
- - - - - - - - - -
hdf_file : hdf5 archive
hd5 file
things_to_save : list of strings
List of datasets to be saved into the hdf5 file .
subgrp : string , optional
Name of hdf5 file subgroup in which the data are to be stored .
"""
if not ( mpi . is_master_node ( ) ) :
return # do nothing on nodes
with HDFArchive ( sum_k . hdf_file , ' a ' ) as ar :
2023-03-15 20:35:39 +01:00
if not subgrp in ar :
ar . create_group ( subgrp )
2021-08-04 18:56:09 +02:00
for it , val in things_to_save . items ( ) :
2023-03-15 20:35:39 +01:00
if it in [ " gf_struct_sumk " , " gf_struct_solver " ,
" solver_to_sumk " , " sumk_to_solver " , " solver_to_sumk_block " ] :
2021-08-04 18:56:09 +02:00
warn ( " It is not recommended to save ' {} ' individually. Save ' block_structure ' instead. " . format ( it ) )
ar [ subgrp ] [ it ] = val
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def cellvolume ( lattice_type , lattice_constants , latticeangle ) :
r """
Determines the conventional und primitive unit cell volumes .
Parameters
- - - - - - - - - -
lattice_type : string
Lattice type according to the Wien2k convention ( P , F , B , R , H , CXY , CYZ , CXZ ) .
lattice_constants : list of double
Lattice constants ( a , b , c ) .
lattice angles : list of double
Lattice angles ( : math : ` \alpha , \beta , \gamma ` ) .
Returns
- - - - - - -
vol_c : double
Conventional unit cell volume .
vol_p : double
Primitive unit cell volume .
"""
a = lattice_constants [ 0 ]
b = lattice_constants [ 1 ]
c = lattice_constants [ 2 ]
c_al = numpy . cos ( latticeangle [ 0 ] )
c_be = numpy . cos ( latticeangle [ 1 ] )
c_ga = numpy . cos ( latticeangle [ 2 ] )
vol_c = a * b * c * \
numpy . sqrt ( 1 + 2 * c_al * c_be * c_ga -
c_al * * 2 - c_be * * 2 - c_ga * * 2 )
det = { " P " : 1 , " F " : 4 , " B " : 2 , " R " : 3 ,
" H " : 1 , " CXY " : 2 , " CYZ " : 2 , " CXZ " : 2 }
vol_p = vol_c / det [ lattice_type ]
return vol_c , vol_p
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def fermi_dis ( w , beta , der = 0 ) :
r """
Fermi distribution .
. . math : :
f ( x ) = 1 / ( e ^ x + 1 ) .
Parameters
- - - - - - - - - -
w : double
frequency
beta : double
inverse temperature
der : integer
order of derivative
Returns
- - - - - - -
f : double
"""
2023-10-24 22:14:20 +02:00
exponent = numpy . longdouble ( w * beta )
2021-08-04 18:56:09 +02:00
fermi = 1.0 / ( numpy . exp ( exponent ) + 1 )
if der == 0 :
return fermi
elif der == 1 :
return - beta * fermi * * 2 * numpy . exp ( exponent )
else :
2023-10-24 22:14:20 +02:00
raise ValueError ( ' higher order of derivative than 1 not implemented ' )
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
2023-01-18 20:27:19 +01:00
def recompute_w90_input_on_different_mesh ( sum_k , seedname , nk_optics , pathname = ' ./ ' , calc_velocity = False , calc_inverse_mass = False , oc_select = ' both ' , oc_basis = ' h ' ) :
2021-08-04 18:56:09 +02:00
r """
Recomputes dft_input objects on a finer mesh using WannierBerri and Wannier90 input .
Parameters
- - - - - - - - - -
sum_k : sum_k object
triqs SumkDFT object
seedname : string
Wannier90 seedname
nk_optics : single integer / float or three integers
if single integer given , mesh is [ nk_optics , nk_optics , nk_optics ]
elif single float given , mesh is ceiling of * sum_k . kpts * nk_optics
elif three integers given , mesh is nk_optics
pathname : string , optional , default = ' ./ '
location of Wannier90 data
calc_velocity : boolean , optional , default = False
whether the velocity ( first derivative of H ( k ) ) is computed
calc_inverse_mass : boolean , optional , default = False
whether the inverse effective mass ( second derivative of H ( k ) ) is computed
2023-01-18 20:27:19 +01:00
oc_select : string , optional , default = ' both '
select contributions for optical conductivity from [ ' intra ' , ' inter ' , ' both ' ]
oc_basis : string , optional , default = ' h '
gauge choice options ' h ' for Hamiltonian / band and ' w ' for Wannier basis
2021-08-04 18:56:09 +02:00
Returns
- - - - - - -
sum_k : sum_k object
triqs SumkDFT object
things_to_store : dictionary
dictionary of datasets to be temporarily overwritten
"""
2023-01-18 20:27:19 +01:00
mpi . report ( ' Starting Wannier interpolation... ' )
2021-08-04 18:56:09 +02:00
BOHRTOANG = cst . physical_constants [ ' Bohr radius ' ] [ 0 ] / cst . angstrom
HARTREETOEV = cst . physical_constants [ ' Hartree energy ' ] [ 0 ] / cst . eV
n_inequiv_spin_blocks = sum_k . SP + 1 - sum_k . SO
# set-up k mesh depending on input shape
# read in transport input and some checks
2023-03-15 20:33:49 +01:00
read_transport_input_from_hdf ( sum_k )
2021-08-04 18:56:09 +02:00
# first check for right formatting of sum_k.nk_optics
2023-03-15 20:35:39 +01:00
assert len ( nk_optics ) in [ 1 , 3 ] , ' " nk_optics " must be given as three integers or one float '
if len ( nk_optics ) == 1 :
assert numpy . array ( list ( nk_optics ) ) . dtype in (
int , float ) , ' " nk_optics " single value must be float or integer '
if len ( nk_optics ) == 3 :
assert numpy . array ( list ( nk_optics ) ) . dtype == int , ' " nk_optics " mesh must be integers '
2021-08-04 18:56:09 +02:00
if len ( nk_optics ) == 1 :
interpolate_factor = nk_optics [ 0 ]
2023-03-15 20:35:39 +01:00
nk_x , nk_y , nk_z = list (
map ( lambda i : int ( numpy . ceil ( interpolate_factor * len ( set ( sum_k . kpts [ : , i ] ) ) ) ) , range ( 3 ) ) )
2021-08-04 18:56:09 +02:00
else :
nk_x , nk_y , nk_z = nk_optics
2023-03-15 20:33:49 +01:00
# check for spin calculation (not supported)
assert sum_k . SP == 0 , ' spin dependent transport calculations are not supported. '
2021-08-04 18:56:09 +02:00
n_orb = numpy . max ( [ sum_k . n_orbitals [ ik ] [ 0 ] for ik in range ( sum_k . n_k ) ] )
# temporarily recompute the following quantities on a different mesh
things_to_modify = { ' bz_weights ' : None , ' hopping ' : None , ' kpt_weights ' : None , ' kpts ' : None ,
' n_k ' : None , ' n_orbitals ' : None , ' proj_mat ' : None , ' band_window ' : None , ' band_window_optics ' : None }
things_to_store = dict . fromkeys ( things_to_modify , None )
# initialize variables
n_kpts = nk_x * nk_y * nk_z
kpts = numpy . zeros ( ( n_kpts , 3 ) )
hopping = numpy . zeros ( ( n_kpts , 1 , n_orb , n_orb ) , dtype = complex )
2023-03-15 20:35:39 +01:00
proj_mat = numpy . zeros ( numpy . shape (
hopping [ : , 0 , 0 , 0 ] ) + numpy . shape ( sum_k . proj_mat [ 0 , : ] ) , dtype = complex )
2021-08-04 18:56:09 +02:00
cell_volume = kpts = None
2023-03-15 20:35:39 +01:00
if calc_velocity :
velocities_k = None
if calc_inverse_mass :
inverse_mass = None
2021-08-04 18:56:09 +02:00
if mpi . is_master_node ( ) :
# try wannierberri import
try :
import wannierberri as wb
except ImportError :
2024-07-17 16:30:09 +02:00
print ( ' ImportError: WannierBerri needs to be installed to run optics calculations with Wannier90 ' )
2021-08-04 18:56:09 +02:00
try :
mpi . MPI . COMM_WORLD . Abort ( 1 )
except :
sys . exit ( )
# initialize WannierBerri system
2023-03-15 20:35:39 +01:00
shift_gamma = numpy . array ( [ 0.0 , 0.0 , 0.0 ] )
# wberri = wb.System_w90(pathname + seedname, berry=True, fft='numpy')
2021-08-04 18:56:09 +02:00
# WannierBerri uses python multiprocessing which might conflict with mpi.
# if there's a segfault, uncomment the following line
2024-08-15 20:12:14 +02:00
wberri = wb . System_w90 ( pathname + seedname , berry = True , fftlib = ' numpy ' , npar = 16 )
2021-08-04 18:56:09 +02:00
grid = wb . Grid ( wberri , NKdiv = 1 , NKFFT = [ nk_x , nk_y , nk_z ] )
2024-07-17 16:30:09 +02:00
dataK = wb . data_K . Data_K_R ( wberri , dK = shift_gamma , grid = grid , fftlib = ' numpy ' )
2021-08-04 18:56:09 +02:00
2023-03-15 20:33:49 +01:00
assert dataK . HH_K . shape == hopping [ : , 0 , : , : ] . shape , ' wberri / wannier Hamiltonian has different number of orbitals than SumK object. Disentanglement is not supported as of now. '
2021-08-04 18:56:09 +02:00
# read in hoppings and proj_mat
2023-01-18 20:27:19 +01:00
if oc_basis == ' h ' :
2023-03-15 20:35:39 +01:00
hopping [ : , 0 , range ( hopping . shape [ 2 ] ) , range ( hopping . shape [ 3 ] ) ] = dataK . E_K
2023-01-18 20:27:19 +01:00
elif oc_basis == ' w ' :
2023-03-15 20:35:39 +01:00
hopping [ : , 0 , : , : ] = dataK . HH_K
2023-01-18 20:27:19 +01:00
fake_proj_mat = numpy . zeros ( numpy . shape ( dataK . UU_K ) , dtype = complex )
2023-03-15 20:35:39 +01:00
fake_proj_mat [ : , range ( numpy . shape ( fake_proj_mat ) [ 1 ] ) , range (
numpy . shape ( fake_proj_mat ) [ 2 ] ) ] = 1. + 1 j * 0.0
2023-01-18 20:27:19 +01:00
2021-08-04 18:56:09 +02:00
for isp in range ( n_inequiv_spin_blocks ) :
iorb = 0
for icrsh in range ( sum_k . n_corr_shells ) :
dim = sum_k . corr_shells [ icrsh ] [ ' dim ' ]
2023-01-18 20:27:19 +01:00
if oc_basis == ' h ' :
2023-03-15 20:35:39 +01:00
proj_mat [ : , isp , icrsh , 0 : dim , : ] = dataK . UU_K [ : , iorb : iorb + dim , : ]
2023-01-18 20:27:19 +01:00
elif oc_basis == ' w ' :
2023-03-15 20:35:39 +01:00
proj_mat [ : , isp , icrsh , 0 : dim , : ] = fake_proj_mat [ : , iorb : iorb + dim , : ]
2021-08-04 18:56:09 +02:00
iorb + = dim
if calc_velocity :
2023-01-18 20:27:19 +01:00
# velocity: [k x n_orb x n_orb x R]
def _commutator ( A , B ) :
term1 = numpy . einsum ( ' kmo, kona -> kmna ' , A , B )
term2 = numpy . einsum ( ' kmoa, kon -> kmna ' , B , A )
return term1 - term2
# in the band basis
# vh_alpha = Hhbar_alpha + i [Hh, Ahbar_alpha]
if oc_basis == ' h ' :
# first term
Hhbar_alpha = dataK . Xbar ( ' Ham ' , 1 )
# second term
2023-03-15 20:35:39 +01:00
c_Hh_Ahbar_alpha = _commutator ( hopping [ : , 0 , : , : ] , dataK . Xbar ( ' AA ' ) )
2023-01-18 20:27:19 +01:00
velocities_k = ( Hhbar_alpha + 1 j * c_Hh_Ahbar_alpha ) / HARTREETOEV / BOHRTOANG
# split into diag and offdiag elements, corresponding to intra- and interband contributions
v_diag = numpy . zeros ( numpy . shape ( velocities_k ) , dtype = complex )
2023-03-15 20:35:39 +01:00
v_diag [ : , range ( numpy . shape ( velocities_k ) [ 1 ] ) ,
range ( numpy . shape ( velocities_k ) [ 2 ] ) , : ] = velocities_k [ : , range ( numpy . shape ( velocities_k ) [ 1 ] ) ,
range ( numpy . shape ( velocities_k ) [ 2 ] ) , : ]
2023-01-18 20:27:19 +01:00
v_offdiag = velocities_k . copy ( )
2023-03-15 20:35:39 +01:00
v_offdiag [ : , range ( numpy . shape ( velocities_k ) [ 1 ] ) , range (
numpy . shape ( velocities_k ) [ 2 ] ) , : ] = 0. + 1 j * 0.0
2023-01-18 20:27:19 +01:00
if oc_select == ' intra ' :
velocities_k = v_diag
elif oc_select == ' inter ' :
velocities_k = v_offdiag
elif oc_select == ' both ' :
velocities_k = v_diag + v_offdiag
# in the orbital basis
# vw_alpha = Hw_alpha + i [Hw, Aw_alpha]
elif oc_basis == ' w ' :
2024-08-15 20:12:14 +02:00
# first term, taken from
# github.com/wannier-berri/wannier-berri/blob/2d3982331c02775f5ee033c664849d5f2d41d0c1/wannierberri/data_K.py#L687
Hw_alpha = wb . data_K . Data_K_R . _R_to_k_H ( dataK , dataK . Ham_R , der = 1 , hermitian = True )
2023-01-18 20:27:19 +01:00
# second term
2024-08-15 20:12:14 +02:00
Aw_alpha = dataK . fft_R_to_k ( dataK . get_R_mat ( ' AA ' ) , hermitian = True )
2023-03-15 20:35:39 +01:00
c_Hw_Aw_alpha = _commutator ( hopping [ : , 0 , : , : ] , Aw_alpha )
2023-01-18 20:27:19 +01:00
velocities_k = ( Hw_alpha + 1 j * c_Hw_Aw_alpha ) / HARTREETOEV / BOHRTOANG
2021-08-04 18:56:09 +02:00
if calc_inverse_mass :
2023-10-04 13:55:56 +02:00
# in the band basis
# ToDo: change units of inverse_mass consistently
if oc_basis == ' h ' :
inverse_mass = dataK . Xbar ( ' Ham ' , 2 )
# in the orbital basis
elif oc_basis == ' w ' :
Hw_alphabeta_R = dataK . Ham_R . copy ( )
for i in range ( 2 ) :
shape_cR = numpy . shape ( dataK . cRvec_wcc )
Hw_alphabeta_R = 1 j * Hw_alphabeta_R . reshape ( ( Hw_alphabeta_R . shape ) + ( 1 , ) ) * dataK . cRvec_wcc . reshape (
( shape_cR [ 0 ] , shape_cR [ 1 ] , dataK . system . nRvec ) + ( 1 , ) * len ( Hw_alphabeta_R . shape [ 3 : ] ) + ( 3 , ) )
inverse_mass = dataK . fft_R_to_k ( Hw_alphabeta_R , hermitean = False ) [ dataK . select_K ]
2021-08-04 18:56:09 +02:00
# read in rest from dataK
cell_volume = dataK . cell_volume / BOHRTOANG * * 3
kpts = dataK . kpoints_all
# broadcast everything
2022-11-02 15:48:44 +01:00
sum_k . cell_vol = mpi . bcast ( cell_volume )
2021-08-04 18:56:09 +02:00
kpts = mpi . bcast ( kpts )
hopping = mpi . bcast ( hopping )
proj_mat = mpi . bcast ( proj_mat )
2023-03-15 20:35:39 +01:00
if calc_velocity :
velocities_k = mpi . bcast ( velocities_k )
if calc_inverse_mass :
inverse_mass = mpi . bcast ( inverse_mass )
2021-08-04 18:56:09 +02:00
# write interpolated sumk quantities into "things_to_modify"
things_to_modify [ ' n_k ' ] = n_kpts
things_to_modify [ ' n_orbitals ' ] = numpy . full ( ( n_kpts , 1 ) , n_orb )
for key in [ ' bz_weights ' , ' kpt_weights ' ] :
things_to_modify [ key ] = numpy . full ( n_kpts , 1 / n_kpts )
n_inequiv_spin_blocks = sum_k . SP + 1 - sum_k . SO
for key in [ ' band_window ' , ' band_window_optics ' ] :
2023-03-15 20:35:39 +01:00
things_to_modify [ key ] = [ numpy . full ( ( n_kpts , 2 ) , sum_k . band_window [ isp ] [ 0 ] )
for isp in range ( n_inequiv_spin_blocks ) ]
2021-08-04 18:56:09 +02:00
things_to_modify [ ' kpts ' ] = kpts
things_to_modify [ ' hopping ' ] = hopping
things_to_modify [ ' proj_mat ' ] = proj_mat
if calc_velocity :
sum_k . velocities_k = None
things_to_modify [ ' velocities_k ' ] = velocities_k
if calc_inverse_mass :
sum_k . inverse_mass = None
things_to_modify [ ' inverse_mass ' ] = inverse_mass
# now save previous sum_k instances into "things_to_store" and overwrite
# TODO: decide whether this should be undone after the run
for key in things_to_modify :
things_to_store [ key ] = getattr ( sum_k , key )
setattr ( sum_k , key , things_to_modify [ key ] )
# write velocities to file
if calc_velocity :
if mpi . is_master_node ( ) :
ar = HDFArchive ( sum_k . hdf_file , ' a ' )
ar [ ' dft_transp_input ' ] [ ' velocities_k ' ] = velocities_k
if calc_inverse_mass :
if mpi . is_master_node ( ) :
ar = HDFArchive ( sum_k . hdf_file , ' a ' )
ar [ ' dft_transp_input ' ] [ ' inverse_mass ' ] = inverse_mass
2022-11-02 15:48:44 +01:00
return sum_k , things_to_store
2021-08-04 18:56:09 +02:00
# ----------------- transport -----------------------
2023-10-04 14:12:51 +02:00
def raman_vertex ( sumk , ik , direction , code , isp , options = None ) :
r """
Reads all necessary quantities for transport calculations from transport subgroup of the hdf5 archive .
Performs checks on input . Uses interpolation if code = wannier90 .
Parameters
- - - - - - - - - -
sum_k : sum_k object
triqs SumkDFT object
ik : integer
the ik index in which will be calculate the raman vertex
direction : string
direction to calculate the Raman vertex .
code : string
DFT code from which velocities are being read . Options : ' wien2k ' , ' wannier90 '
isp : integer
spin index
options : dictionary , optional
additional keywords necessary in case a new direction wants to be implemented
Returns
- - - - - - -
raman_vertex : 2 D array
numpy array size [ n_orb , n_orb ] of complex type
"""
if code in ( ' wien2k ' ) :
assert 0 , ' Raman for wien2k not yet implemented ' #ToDo
dir_names = [ ' xx ' , ' yy ' , ' zz ' , ' B2g ' , ' B1g ' , ' A1g ' , ' Eg ' ]
dir_array = [ [ [ 1. , 0. , 0. ] , [ 0. , 0. , 0. ] , [ 0. , 0. , 0. ] ] ,
[ [ 0. , 0. , 0. ] , [ 0. , 1. , 0. ] , [ 0. , 0. , 0. ] ] ,
[ [ 0. , 0. , 0. ] , [ 0. , 0. , 0. ] , [ 0. , 0. , 1. ] ] ,
[ [ 0. , 1. , 0. ] , [ 0. , 0. , 0. ] , [ 0. , 0. , 0. ] ] ,
[ [ 0.5 , 0. , 0. ] , [ 0. , - 0.5 , 0. ] , [ 0. , 0. , 0. ] ] ,
[ [ 0. , 0. , 0. ] , [ 0. , 0. , 0. ] , [ 0. , 0. , 1. ] ] ,
[ [ 0. , 0. , 1. ] , [ 0. , 0. , 0. ] , [ 0. , 0. , 0. ] ] ]
# Load custom directions
if " custom_dir " in options :
assert isinstance ( options [ " custom_dir " ] , dict ) , " raman_vertex: in options, custom_dir must be a dictionary "
for dire in options [ " custom_dir " ] :
assert numpy . shape ( numpy . array ( options [ " custom_dir " ] [ dire ] ) ) == ( 3 , 3 ) , " raman_vertex: custom_dir must have shape 3x3 (a numpy array or a nested list) "
if dire in dir_names :
if direction == dire and ik == 0 : #reports only when ik==0
mpi . report ( " Warning: the direction %s was already loaded and will be replace with the one provided in custom_dir " % dire )
idir = dir_names . index ( dire )
dir_array [ idir ] = options [ " custom_dir " ] [ dire ]
else :
dir_names . append ( dire )
dir_array . append ( options [ " custom_dir " ] [ dire ] )
dir_array = [ numpy . array ( el , dtype = numpy . float_ ) for el in dir_array ] # convert list to numpy array
if direction in dir_names :
idir = dir_names . index ( direction )
else :
assert 0 , ' raman_vertex: direction %s is not supported by default. Try to add it by raman_vertex options custom_dir. ' % direction
n_bands = sumk . n_orbitals [ ik ] [ isp ]
ram_vert = numpy . zeros ( ( n_bands , n_bands ) , dtype = numpy . complex_ )
for i in range ( n_bands ) :
for j in range ( n_bands ) :
ram_vert [ i , j ] = numpy . dot ( dir_array [ idir ] , sumk . inverse_mass [ ik , i , j , : , : ] ) . trace ( )
return ram_vert
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def init_spectroscopy ( sum_k , code = ' wien2k ' , w90_params = { } ) :
r """
Reads all necessary quantities for transport calculations from transport subgroup of the hdf5 archive .
Performs checks on input . Uses interpolation if code = wannier90 .
Parameters
- - - - - - - - - -
sum_k : sum_k object
triqs SumkDFT object
code : string
DFT code from which velocities are being read . Options : ' wien2k ' , ' wannier90 '
w90_params : dictionary , optional
additional keywords necessary in case code == ' wannier90 '
Returns
- - - - - - -
sum_k : sum_k object
triqs SumkDFT object , interpolated
"""
2023-01-18 20:27:19 +01:00
mpi . report ( ' Initializing optical conductivity... ' )
2021-08-04 18:56:09 +02:00
# up and down are equivalent if SP = 0
2023-01-18 20:27:19 +01:00
n_inequiv_spin_blocks = sum_k . SP + 1 - sum_k . SO
2021-08-04 18:56:09 +02:00
# ----------------- set-up input from DFT -----------------------
2022-11-02 15:48:44 +01:00
if code in ( ' wien2k ' , ' elk ' ) :
2021-08-04 18:56:09 +02:00
# Check if wien converter was called and read transport subgroup form
# hdf file
if mpi . is_master_node ( ) :
ar = HDFArchive ( sum_k . hdf_file , ' r ' )
if not ( sum_k . transp_data in ar ) :
2023-03-15 20:35:39 +01:00
raise IOError (
" transport_distribution: No %s subgroup in hdf file found! Call convert_transp_input first. " % sum_k . transp_data )
2021-08-04 18:56:09 +02:00
# check if outputs file was converted
if not ( ' n_symmetries ' in ar [ ' dft_misc_input ' ] ) :
2023-03-15 20:35:39 +01:00
raise IOError (
" transport_distribution: n_symmetries missing. Check if case.outputs file is present and call convert_misc_input() or convert_dft_input(). " )
2021-08-04 18:56:09 +02:00
sum_k = read_transport_input_from_hdf ( sum_k )
elif code in ( ' wannier90 ' ) :
required_entries = [ ' seedname ' , ' nk_optics ' ]
assert all ( entry in w90_params for entry in required_entries ) , ' Please provide additional keywords " seedname " and " nk_optics " for " code = " wannier90 " " '
# check if spin-unpolarized
assert n_inequiv_spin_blocks == 1 , " Spin-polarized optical conductivity calculations not implemented with Wannier90 "
# check some of the input
pathname = w90_params [ ' pathname ' ] if ' pathname ' in w90_params else ' ./ '
2023-03-15 20:35:39 +01:00
assert all ( isinstance ( name , str ) for name in [
' seedname ' , ' pathname ' ] ) , f ' Check pathname { w90_params [ " pathname " ] } and seedname { w90_params [ " seedname " ] } '
2021-08-04 18:56:09 +02:00
for file_ending in [ ' .wout ' , ' _hr.dat ' , ' .chk ' , ' .mmn ' , ' .eig ' ] :
filename = [ pathname , w90_params [ ' seedname ' ] , file_ending ]
2023-03-15 20:35:39 +01:00
assert os . path . isfile (
' ' . join ( filename ) ) , f ' Filename { " " . join ( filename ) } does not exist! '
2021-08-04 18:56:09 +02:00
calc_velocity = w90_params [ ' calc_velocity ' ] if ' calc_velocity ' in w90_params else True
calc_inverse_mass = w90_params [ ' calc_inverse_mass ' ] if ' calc_inverse_mass ' in w90_params else False
2023-03-15 20:35:39 +01:00
assert all ( isinstance ( name , bool ) for name in [
calc_velocity , calc_inverse_mass ] ) , f ' Parameter { calc_velocity } or { calc_inverse_mass } not bool! '
2021-08-04 18:56:09 +02:00
2023-01-18 20:27:19 +01:00
# select contributions to be used
oc_select = w90_params [ ' oc_select ' ] if ' oc_select ' in w90_params else ' both '
2023-03-15 20:35:39 +01:00
assert oc_select in [ ' intra ' , ' inter ' ,
' both ' ] , ' " oc_select " needs to be either [ " intra " , " inter " , " both " ] '
2023-01-18 20:27:19 +01:00
# gauge choice options 'h' for Hamiltonian and 'w' for Wannier
oc_basis = w90_params [ ' oc_basis ' ] if ' oc_basis ' in w90_params else ' h '
assert oc_basis in [ ' h ' , ' w ' ] , ' " oc_basis " needs to be either [ " h " , " w " ] '
# finally, make sure oc_select is 'both' for oc_basis = 'w'
if oc_basis == ' w ' and oc_select != ' both ' :
warn ( f ' " oc_select " must be " both " for " oc_basis " = " w " ! ' )
oc_select = ' both '
# further checks for calc_inverse_mass
if calc_inverse_mass :
assert oc_basis == ' h ' , ' " calc_inverse_mass " only implemented for " oc_basis " == " h " '
assert oc_select == ' both ' , ' " oc_select " not implemented for " calc_inverse_mass " '
# print some information
mpi . report ( f ' { " Basis choice [h (Hamiltonian), w (Wannier)]: " : <60s } { oc_basis } ' )
mpi . report ( f ' { " Contributions from [intra(-band), inter(-band), both]: " : <60s } { oc_select } ' )
2021-08-04 18:56:09 +02:00
# recompute sum_k instances on denser grid
2022-11-02 15:48:44 +01:00
sum_k , _ = recompute_w90_input_on_different_mesh ( sum_k , w90_params [ ' seedname ' ] , nk_optics = w90_params [ ' nk_optics ' ] , pathname = pathname ,
2023-03-15 20:35:39 +01:00
calc_velocity = calc_velocity , calc_inverse_mass = calc_inverse_mass , oc_select = oc_select , oc_basis = oc_basis )
2021-08-04 18:56:09 +02:00
# k-dependent-projections.
# to be checked. But this should be obsolete atm, works for both cases
# k_dep_projection is nowhere used
# assert sum_k.k_dep_projection == 0, "transport_distribution: k dependent projection is not implemented!"
2022-11-02 15:48:44 +01:00
return sum_k
2021-08-04 18:56:09 +02:00
# Uses .data of only GfReFreq objects.
2023-03-15 20:35:39 +01:00
2022-11-02 15:48:44 +01:00
def transport_distribution ( sum_k , beta , directions = [ ' xx ' ] , energy_window = None , Om_mesh = [ 0.0 ] , with_Sigma = False , n_om = None , broadening = 0.0 , code = ' wien2k ' ) :
2021-08-04 18:56:09 +02:00
r """
Calculates the transport distribution
. . math : :
\Gamma_ { \alpha \beta } \left ( \omega + \Omega / 2 , \omega - \Omega / 2 \right ) = \frac { 1 } { V } \sum_k Tr \left ( v_ { k , \alpha } A_ { k } ( \omega + \Omega / 2 ) v_ { k , \beta } A_ { k } \left ( \omega - \Omega / 2 \right ) \right )
in the direction : math : ` \alpha \beta ` . The velocities : math : ` v_ { k } ` are read from the transport subgroup of the hdf5 archive .
Parameters
- - - - - - - - - -
sum_k : sum_k object
triqs SumkDFT object
beta : double
Inverse temperature : math : ` \beta ` .
directions : list of string , optional
: math : ` \alpha \beta ` e . g . : [ ' xx ' , ' yy ' , ' zz ' , ' xy ' , ' xz ' , ' yz ' ] .
energy_window : list of double , optional
Specifies the upper and lower limit of the frequency integration for : math : ` \Omega = 0.0 ` . The window is automatically enlarged by the largest : math : ` \Omega ` value ,
hence the integration is performed in the interval [ energy_window [ 0 ] - max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) ] .
Om_mesh : list of double , optional
: math : ` \Omega ` frequency mesh of the optical conductivity . For the conductivity and the Seebeck coefficient : math : ` \Omega = 0.0 ` has to be
part of the mesh . In the current version Om_mesh is repined to the mesh provided by the self - energy ! The actual mesh is printed on the screen and given as output .
with_Sigma : boolean , optional
Determines whether the calculation is performed with or without self energy . If this parameter is set to False the self energy is set to zero ( i . e . the DFT band
structure : math : ` A ( k , \omega ) ` is used ) . Note : For with_Sigma = False it is necessary to specify the parameters energy_window , n_om and broadening .
n_om : integer , optional
Number of equidistant frequency points in the interval [ energy_window [ 0 ] - max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) ] . This parameters is only used if
with_Sigma = False .
broadening : double , optional
Lorentzian broadening . It is necessary to specify the boradening if with_Sigma = False , otherwise this parameter can be set to 0.0 .
code : string
DFT code from which velocities are being read . Options : ' wien2k ' , ' wannier90 '
Returns
- - - - - - -
Gamma_w : dictionary of double matrices
transport distribution function in each direction , frequency given by Om_mesh_out and omega
omega : list of double
omega vector
Om_mesh_out : list of double
frequency mesh of the optical conductivity recomputed on the mesh provided by the self energy
"""
2023-01-18 20:27:19 +01:00
mpi . report ( ' Computing transport distribution... ' )
2021-08-04 18:56:09 +02:00
n_inequiv_spin_blocks = sum_k . SP + 1 - sum_k . SO
# up and down are equivalent if SP = 0
# positive om_mesh
2023-03-15 20:35:39 +01:00
assert all (
Om > = 0.0 for Om in Om_mesh ) , " transport_distribution: Om_mesh should not contain negative values! "
2021-08-04 18:56:09 +02:00
# Check if energy_window is sufficiently large and correct
if ( energy_window [ 0 ] > = energy_window [ 1 ] or energy_window [ 0 ] > = 0 or energy_window [ 1 ] < = 0 ) :
assert 0 , " transport_distribution: energy_window wrong! "
if ( abs ( fermi_dis ( energy_window [ 0 ] , beta ) * fermi_dis ( - energy_window [ 0 ] , beta ) ) > 1e-5
or abs ( fermi_dis ( energy_window [ 1 ] , beta ) * fermi_dis ( - energy_window [ 1 ] , beta ) ) > 1e-5 ) :
mpi . report (
" \n #################################################################### " )
mpi . report (
" transport_distribution: WARNING - energy window might be too narrow! " )
mpi . report (
" #################################################################### \n " )
# ----------------- calculate A(k,w) -----------------------
# Define mesh for Green's function and in the specified energy window
if ( with_Sigma == True ) :
omega = numpy . array ( [ round ( x . real , 12 )
2023-03-15 20:35:39 +01:00
for x in sum_k . Sigma_imp [ 0 ] . mesh ] )
2021-08-04 18:56:09 +02:00
mesh = None
mu = sum_k . chemical_potential
n_om = len ( omega )
mpi . report ( " Using omega mesh provided by Sigma! " )
if energy_window :
# Find according window in Sigma mesh
ioffset = numpy . sum (
omega < energy_window [ 0 ] - max ( Om_mesh ) )
2023-03-15 20:35:39 +01:00
omega = omega [ numpy . logical_and (
omega > = energy_window [ 0 ] - max ( Om_mesh ) , omega < = energy_window [ 1 ] + max ( Om_mesh ) ) ]
2021-08-04 18:56:09 +02:00
n_om = len ( omega )
# Truncate Sigma to given omega window
# In the future there should be an option in gf to manipulate the mesh (e.g. truncate) directly.
# For now we stick with this:
for icrsh in range ( sum_k . n_corr_shells ) :
Sigma_save = sum_k . Sigma_imp [ icrsh ] . copy ( )
spn = sum_k . spin_block_names [ sum_k . corr_shells [ icrsh ] [ ' SO ' ] ]
2023-03-15 20:35:39 +01:00
def glist ( ) : return [ GfReFreq ( target_shape = ( block_dim , block_dim ) , window = ( omega [
0 ] , omega [ - 1 ] ) , n_points = n_om ) for block , block_dim in sum_k . gf_struct_sumk [ icrsh ] ]
2021-08-04 18:56:09 +02:00
sum_k . Sigma_imp [ icrsh ] = BlockGf (
name_list = spn , block_list = glist ( ) , make_copies = False )
for i , g in sum_k . Sigma_imp [ icrsh ] :
for iL in g . indices [ 0 ] :
for iR in g . indices [ 0 ] :
for iom in range ( n_om ) :
g . data [ iom , int ( iL ) , int ( iR ) ] = Sigma_save [
i ] . data [ ioffset + iom , int ( iL ) , int ( iR ) ]
else :
assert n_om is not None , " transport_distribution: Number of omega points (n_om) needed to calculate transport distribution! "
assert energy_window is not None , " transport_distribution: Energy window needed to calculate transport distribution! "
assert broadening != 0.0 and broadening is not None , " transport_distribution: Broadening necessary to calculate transport distribution! "
omega = numpy . linspace (
energy_window [ 0 ] - max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) , n_om )
2022-11-02 15:48:44 +01:00
mesh = MeshReFreq ( energy_window [ 0 ] -
2023-03-15 20:35:39 +01:00
max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) , n_om )
2021-08-04 18:56:09 +02:00
mu = 0.0
dir_to_int = { ' x ' : 0 , ' y ' : 1 , ' z ' : 2 }
# Define mesh for optic conductivity
d_omega = round ( numpy . abs ( omega [ 0 ] - omega [ 1 ] ) , 12 )
iOm_mesh = numpy . array ( [ round ( ( Om / d_omega ) , 0 ) for Om in Om_mesh ] )
temp_Om_mesh = iOm_mesh * d_omega
if mpi . is_master_node ( ) :
print ( " Chemical potential: " , mu )
2023-03-15 20:35:39 +01:00
print ( " Using n_om = %s points in the energy_window [ %s , %s ] " % (
n_om , omega [ 0 ] , omega [ - 1 ] ) , end = ' ' )
2021-08-04 18:56:09 +02:00
print ( " where the omega vector is: " )
print ( omega )
print ( " Calculation requested for Omega mesh: " , numpy . array ( Om_mesh ) )
print ( " Omega mesh automatically repined to: " , temp_Om_mesh )
2023-03-15 20:35:39 +01:00
Gamma_w = { direction : numpy . zeros ( ( len ( temp_Om_mesh ) , n_om ) ,
2024-06-21 17:47:39 +02:00
dtype = numpy . float64 ) for direction in directions }
2021-08-04 18:56:09 +02:00
# Sum over all k-points
ikarray = numpy . array ( list ( range ( sum_k . n_k ) ) )
for ik in mpi . slice_array ( ikarray ) :
# Calculate G_w for ik and initialize A_kw
G_w = sum_k . lattice_gf ( ik , mu , broadening = broadening , mesh = mesh , with_Sigma = with_Sigma )
2024-06-21 17:47:39 +02:00
A_kw = [ numpy . zeros ( ( sum_k . n_orbitals [ ik ] [ isp ] , sum_k . n_orbitals [ ik ] [ isp ] , n_om ) , dtype = numpy . complex128 )
2021-08-04 18:56:09 +02:00
for isp in range ( n_inequiv_spin_blocks ) ]
for isp in range ( n_inequiv_spin_blocks ) :
# copy data from G_w (swapaxes is used to have omega in the 3rd
# dimension)
A_kw [ isp ] = copy . deepcopy ( G_w [ sum_k . spin_block_names [ sum_k . SO ] [
isp ] ] . data . swapaxes ( 0 , 1 ) . swapaxes ( 1 , 2 ) )
# calculate A(k,w) for each frequency
for iw in range ( n_om ) :
A_kw [ isp ] [ : , : , iw ] = - 1.0 / ( 2.0 * numpy . pi * 1 j ) * (
A_kw [ isp ] [ : , : , iw ] - numpy . conjugate ( numpy . transpose ( A_kw [ isp ] [ : , : , iw ] ) ) )
2023-03-15 20:35:39 +01:00
# Akw_write[ik] = A_kw[isp].copy() * sum_k.bz_weights[ik]
2021-08-04 18:56:09 +02:00
b_min = max ( sum_k . band_window [ isp ] [
ik , 0 ] , sum_k . band_window_optics [ isp ] [ ik , 0 ] )
b_max = min ( sum_k . band_window [ isp ] [
ik , 1 ] , sum_k . band_window_optics [ isp ] [ ik , 1 ] )
A_i = slice (
b_min - sum_k . band_window [ isp ] [ ik , 0 ] , b_max - sum_k . band_window [ isp ] [ ik , 0 ] + 1 )
v_i = slice ( b_min - sum_k . band_window_optics [ isp ] [
ik , 0 ] , b_max - sum_k . band_window_optics [ isp ] [ ik , 0 ] + 1 )
# loop over all symmetries
for R in sum_k . rot_symmetries :
# get transformed velocity under symmetry R
if code in ( ' wien2k ' ) :
vel_R = copy . deepcopy ( sum_k . velocities_k [ isp ] [ ik ] )
elif code in ( ' wannier90 ' ) :
vel_R = copy . deepcopy ( sum_k . velocities_k [ ik ] )
for nu1 in range ( sum_k . band_window_optics [ isp ] [ ik , 1 ] - sum_k . band_window_optics [ isp ] [ ik , 0 ] + 1 ) :
for nu2 in range ( sum_k . band_window_optics [ isp ] [ ik , 1 ] - sum_k . band_window_optics [ isp ] [ ik , 0 ] + 1 ) :
vel_R [ nu1 ] [ nu2 ] [ : ] = numpy . dot (
R , vel_R [ nu1 ] [ nu2 ] [ : ] )
# calculate Gamma_w for each direction from the velocities
# vel_R and the spectral function A_kw
for direction in directions :
for iw in range ( n_om ) :
for iq in range ( len ( temp_Om_mesh ) ) :
2023-03-15 20:35:39 +01:00
if ( iw + iOm_mesh [ iq ] > = n_om or omega [ iw ] < - temp_Om_mesh [ iq ] + energy_window [ 0 ] or omega [ iw ] > temp_Om_mesh [ iq ] + energy_window [ 1 ] ) :
2021-08-04 18:56:09 +02:00
continue
Gamma_w [ direction ] [ iq , iw ] + = ( numpy . dot ( numpy . dot ( numpy . dot ( vel_R [ v_i , v_i , dir_to_int [ direction [ 0 ] ] ] , A_kw [ isp ] [ A_i , A_i , int ( iw + iOm_mesh [ iq ] ) ] ) ,
2023-03-15 20:35:39 +01:00
vel_R [ v_i , v_i , dir_to_int [ direction [ 1 ] ] ] ) , A_kw [ isp ] [ A_i , A_i , iw ] ) . trace ( ) . real * sum_k . bz_weights [ ik ] )
2021-08-04 18:56:09 +02:00
for direction in directions :
2023-06-02 16:51:48 +02:00
Gamma_w [ direction ] = ( mpi . all_reduce ( Gamma_w [ direction ] ) / sum_k . cell_vol / sum_k . n_symmetries )
2021-08-04 18:56:09 +02:00
return Gamma_w , omega , temp_Om_mesh
2023-03-15 20:35:39 +01:00
2023-01-18 20:27:19 +01:00
def transport_function ( beta , directions , hopping , velocities , energy_window , n_om , rot_symmetries ) :
2021-08-04 18:56:09 +02:00
r """
Calculates the transport function
. . math : :
\Phi_ \alpha \beta ( \omega ) = \sum_k v_ { k , \alpha } v_ { k , \beta } \delta ( \omega - \varepsilon )
in the direction : math : ` \alpha \beta ` .
Parameters
- - - - - - - - - -
beta : double
Inverse temperature : math : ` \beta ` .
directions : list of string , optional
: math : ` \alpha \beta ` e . g . : [ ' xx ' , ' yy ' , ' zz ' , ' xy ' , ' xz ' , ' yz ' ] .
hopping : double array
Hamiltonian in band basis : math : ` \epsilon ( k ) `
veolcities : complex array
matrix elements derivative of Hamiltonian : math : ` \frac { d \epsilon ( k ) } { dk } `
energy_window : list of double
Specifies the upper and lower limit of the frequency integration for : math : ` \Omega = 0.0 ` . The window is automatically enlarged by the largest : math : ` \Omega ` value ,
hence the integration is performed in the interval [ energy_window [ 0 ] - max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) ] .
n_om : integer
Number of equidistant frequency points in the interval [ energy_window [ 0 ] - max ( Om_mesh ) , energy_window [ 1 ] + max ( Om_mesh ) ] . This parameters is only used if
with_Sigma = False .
2023-01-18 20:27:19 +01:00
rot_symmetries : list of 3 x 3 matrices
rotational symmetries to restore the full FBZ
2021-08-04 18:56:09 +02:00
Returns
- - - - - - -
transp_func : dictionary of double array
transport function in each direction , frequencies given by energy_window
"""
2023-01-18 20:27:19 +01:00
mpi . report ( ' Computing transport function... ' )
# check that velocities are computed on the FBZ
2023-03-15 20:35:39 +01:00
assert numpy . shape ( rot_symmetries ) [
0 ] == 1 , ' Using symmetries currently not implemented for transport function. '
2023-01-18 20:27:19 +01:00
2021-08-04 18:56:09 +02:00
dir_to_int = { ' x ' : 0 , ' y ' : 1 , ' z ' : 2 }
tol = 1 / beta
orb_1 , orb_2 = velocities . shape [ 1 : 3 ]
ws = numpy . linspace ( energy_window [ 0 ] , energy_window [ 1 ] , n_om )
transp_func = { direction : numpy . zeros ( shape = ( ws . shape [ 0 ] ) ) for direction in directions }
for ct , w in enumerate ( ws ) :
2023-03-15 20:35:39 +01:00
idx = numpy . where ( numpy . abs ( hopping [ : , 0 , range ( orb_1 ) , range ( orb_2 ) ] . real - w ) < = tol )
fermi_wg = fermi_dis ( hopping [ : , 0 , range ( orb_1 ) , range ( orb_2 ) ]
[ idx ] . real - w , beta , 1 ) / fermi_dis ( 0. , beta , 1 )
2021-08-04 18:56:09 +02:00
for direction in directions :
dir_a , dir_b = [ dir_to_int [ x ] for x in direction ]
2023-03-15 20:35:39 +01:00
matrix_product = numpy . einsum (
' kmn, kno -> kmo ' , velocities [ : , : , : , dir_a ] , velocities [ : , : , : , dir_b ] )
transp_func [ direction ] [ ct ] = numpy . sum (
fermi_wg * matrix_product [ : , range ( orb_1 ) , range ( orb_2 ) ] [ idx ] ) . real
2021-08-04 18:56:09 +02:00
return transp_func
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def transport_coefficient ( Gamma_w , omega , Om_mesh , spin_polarization , direction , iq , n , beta , method = None ) :
r """
Calculates the transport coefficient A_n in a given direction for a given : math : ` \Omega ` . The required members ( Gamma_w , directions , Om_mesh ) have to be obtained first
by calling the function : meth : ` transport_distribution < dft . sumk_dft_tools . SumkDFTTools . transport_distribution > ` . For n > 0 A is set to NaN if : math : ` \Omega ` is not 0.0 .
Parameters
- - - - - - - - - -
Gamma_w : dictionary of double matrices
transport distribution function in each direction , frequency given by Om_mesh_out and omega
omega : list of double
omega vector
Om_mesh : list of double
frequency mesh of the optical conductivity recomputed on the mesh provided by the self energy
spin_polarization : integer
Boolean - type integer whether system is spin - polarized or not
direction : string
: math : ` \alpha \beta ` e . g . : ' xx ' , ' yy ' , ' zz ' , ' xy ' , ' xz ' , ' yz ' .
iq : integer
Index of : math : ` \Omega ` point in the member Om_mesh .
n : integer
Number of the desired moment of the transport distribution .
beta : double
Inverse temperature : math : ` \beta ` .
method : string
2024-07-01 12:50:33 +02:00
Integration method : cubic spline and scipy . integrate . quad ( ' quad ' ) , simpson rule ( ' simpson ' ) , trapezoidal rule ( ' trapz ' ) , rectangular integration ( otherwise )
2021-08-04 18:56:09 +02:00
Note that the sampling points of the the self - energy are used !
Returns
- - - - - - -
A : double
Transport coefficient .
"""
2023-03-15 20:33:49 +01:00
from scipy . interpolate import interp1d
2024-07-01 12:50:33 +02:00
from scipy . integrate import simpson , quad
2023-03-15 20:33:49 +01:00
2021-08-04 18:56:09 +02:00
if not ( mpi . is_master_node ( ) ) :
2023-06-02 16:51:48 +02:00
return None
2021-08-04 18:56:09 +02:00
if ( Om_mesh [ iq ] == 0.0 or n == 0.0 ) :
A = 0.0
# setup the integrand
if ( Om_mesh [ iq ] == 0.0 ) :
2023-03-15 20:35:39 +01:00
A_int = Gamma_w [ direction ] [ iq ] * \
( fermi_dis ( omega , beta ) * fermi_dis ( - omega , beta ) ) * ( omega * beta ) * * n
2021-08-04 18:56:09 +02:00
elif ( n == 0.0 ) :
2023-03-15 20:35:39 +01:00
A_int = Gamma_w [ direction ] [ iq ] * ( fermi_dis ( omega , beta ) -
fermi_dis ( omega + Om_mesh [ iq ] , beta ) ) / ( Om_mesh [ iq ] * beta )
2021-08-04 18:56:09 +02:00
# w-integration
if method == ' quad ' :
# quad on interpolated w-points with cubic spline
A_int_interp = interp1d ( omega , A_int , kind = ' cubic ' )
A = quad ( A_int_interp , min ( omega ) , max ( omega ) ,
epsabs = 1.0e-12 , epsrel = 1.0e-12 , limit = 500 )
A = A [ 0 ]
2024-07-01 12:50:33 +02:00
elif method == ' simpson ' :
2021-08-04 18:56:09 +02:00
# simpson rule for w-grid
2024-07-01 12:50:33 +02:00
A = simpson ( A_int , omega )
2021-08-04 18:56:09 +02:00
elif method == ' trapz ' :
# trapezoidal rule for w-grid
A = numpy . trapz ( A_int , omega )
else :
# rectangular integration for w-grid (orignal implementation)
d_w = omega [ 1 ] - omega [ 0 ]
for iw in range ( Gamma_w [ direction ] . shape [ 1 ] ) :
A + = A_int [ iw ] * d_w
A = A * numpy . pi * ( 2.0 - spin_polarization )
else :
A = numpy . nan
return A
2023-03-15 20:35:39 +01:00
2021-08-04 18:56:09 +02:00
def conductivity_and_seebeck ( Gamma_w , omega , Om_mesh , SP , directions , beta , method = None ) :
r """
Calculates the Seebeck coefficient and the optical conductivity by calling
: meth : ` transport_coefficient < dft . sumk_dft_tools . SumkDFTTools . transport_coefficient > ` .
The required members ( Gamma_w , directions , Om_mesh ) have to be obtained first by calling the function
: meth : ` transport_distribution < dft . sumk_dft_tools . SumkDFTTools . transport_distribution > ` .
Parameters
- - - - - - - - - -
Gamma_w : dictionary of double matrices
transport distribution function in each direction , frequency given by Om_mesh_out and omega
omega : list of double
omega vector
Om_mesh : list of double
frequency mesh of the optical conductivity recomputed on the mesh provided by the self energy
spin_polarization : integer
Boolean - type integer whether system is spin - polarized or not
directions : list of string , optional
: math : ` \alpha \beta ` e . g . : [ ' xx ' , ' yy ' , ' zz ' , ' xy ' , ' xz ' , ' yz ' ] .
beta : double
Inverse temperature : math : ` \beta ` .
method : string
2024-07-01 12:50:33 +02:00
Integration method : cubic spline and scipy . integrate . quad ( ' quad ' ) , simpson rule ( ' simpson ' ) , trapezoidal rule ( ' trapz ' ) , rectangular integration ( otherwise )
2021-08-04 18:56:09 +02:00
Note that the sampling points of the the self - energy are used !
Returns
- - - - - - -
optic_cond : dictionary of double vectors
Optical conductivity in each direction and frequency given by Om_mesh .
seebeck : dictionary of double
Seebeck coefficient in each direction . If zero is not present in Om_mesh the Seebeck coefficient is set to NaN .
kappa : dictionary of double .
thermal conductivity in each direction . If zero is not present in Om_mesh the thermal conductivity is set to NaN
"""
2023-01-18 20:27:19 +01:00
mpi . report ( ' Computing optical conductivity and kinetic coefficients... ' )
2021-08-04 18:56:09 +02:00
if not ( mpi . is_master_node ( ) ) :
2023-06-02 16:51:48 +02:00
return None , None , None
2021-08-04 18:56:09 +02:00
n_q = Gamma_w [ directions [ 0 ] ] . shape [ 0 ]
# initialization
A0 = { direction : numpy . full ( ( n_q , ) , numpy . nan ) for direction in directions }
A1 = { direction : numpy . full ( ( n_q , ) , numpy . nan ) for direction in directions }
A2 = { direction : numpy . full ( ( n_q , ) , numpy . nan ) for direction in directions }
optic_cond = { direction : numpy . full ( ( n_q , ) , numpy . nan ) for direction in directions }
seebeck = { direction : numpy . nan for direction in directions }
kappa = { direction : numpy . nan for direction in directions }
for direction in directions :
for iq in range ( n_q ) :
2023-03-15 20:35:39 +01:00
A0 [ direction ] [ iq ] = transport_coefficient (
Gamma_w , omega , Om_mesh , SP , direction , iq = iq , n = 0 , beta = beta , method = method )
A1 [ direction ] [ iq ] = transport_coefficient (
Gamma_w , omega , Om_mesh , SP , direction , iq = iq , n = 1 , beta = beta , method = method )
A2 [ direction ] [ iq ] = transport_coefficient (
Gamma_w , omega , Om_mesh , SP , direction , iq = iq , n = 2 , beta = beta , method = method )
print ( " A_0 in direction %s for Omega = %.2f %e a.u. " %
( direction , Om_mesh [ iq ] , A0 [ direction ] [ iq ] ) )
print ( " A_1 in direction %s for Omega = %.2f %e a.u. " %
( direction , Om_mesh [ iq ] , A1 [ direction ] [ iq ] ) )
print ( " A_2 in direction %s for Omega = %.2f %e a.u. " %
( direction , Om_mesh [ iq ] , A2 [ direction ] [ iq ] ) )
2021-08-04 18:56:09 +02:00
if ~ numpy . isnan ( A1 [ direction ] [ iq ] ) :
# Seebeck and kappa are overwritten if there is more than one Omega =
# 0 in Om_mesh
seebeck [ direction ] = - A1 [ direction ] [ iq ] / A0 [ direction ] [ iq ] * 86.17
2023-03-15 20:35:39 +01:00
kappa [ direction ] = A2 [ direction ] [ iq ] - \
A1 [ direction ] [ iq ] * A1 [ direction ] [ iq ] / A0 [ direction ] [ iq ]
2021-08-04 18:56:09 +02:00
kappa [ direction ] * = 293178.0
# factor for optical conductivity: hbar * velocity_Hartree_to_SI * volume_Hartree_to_SI * m_to_cm * 10^-4 final unit
2023-03-15 20:35:39 +01:00
convert_to_SI = cst . hbar * ( cst . c * cst . fine_structure ) * * 2 * \
( 1 / cst . physical_constants [ ' Bohr radius ' ] [ 0 ] ) * * 3 * 1e-6
2021-08-04 18:56:09 +02:00
optic_cond [ direction ] = beta * convert_to_SI * A0 [ direction ]
for iq in range ( n_q ) :
2023-03-15 20:35:39 +01:00
print ( " Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1 " %
( direction , Om_mesh [ iq ] , optic_cond [ direction ] [ iq ] ) )
2021-08-04 18:56:09 +02:00
if not ( numpy . isnan ( A1 [ direction ] [ iq ] ) ) :
2023-03-15 20:35:39 +01:00
print ( " Seebeck in direction %s for Omega = 0.00 %f x 10^(-6) V/K " %
( direction , seebeck [ direction ] ) )
print ( " kappa in direction %s for Omega = 0.00 %f W/(m * K) " %
( direction , kappa [ direction ] ) )
2021-08-04 18:56:09 +02:00
return optic_cond , seebeck , kappa