""" For epw related utilities """
import numpy as np
import os
from EPWpy.utilities.printing import *
from EPWpy.plotting import plot_polaron_matplotlib as plot_pl
from EPWpy.plotting import plot_wannier_matplotlib as plot_wann
from EPWpy.utilities import kramers_kronig as kk
import re
[docs]
class EPWProperties:
"""
Class representing computed properties from EPW (Electron-Phonon Wannier) calculations.
This class stores and organizes EPW-calculated quantities such as
electron-phonon coupling matrices and polaronic displacements.
It is typically used to read, process, and analyze EPW output data.
Attributes
----------
file : str
Name of the EPW output file containing computed properties.
system : str
Name or prefix of the system being studied.
prefix : str
Alias for the system, used to construct file paths.
energy_range : float, optional
Energy range for properties extraction or analysis (e.g., around Fermi level).
epw_params : dict, optional
Dictionary of EPW parameters used for calculations (can include grid size,
smearing, interpolation settings, etc.).
epw_fold : str
Name of the EPW output folder. Default is 'epw'.
epw_file : str
Full path to the EPW output file, constructed from system, folder, and filename.
Notes
-----
- This class is intended for post-processing EPW output.
- The primary properties handled include:
- `dtau` : polaronic displacements.
- `gkk` : electron-phonon matrix elements.
Examples
--------
>>> props = EPWProperties(file='epw2.out', system='graphene', energy_range=0.5)
>>> print(props.epw_file)
'graphene/epw/epw1.out'
"""
def __init__(self, file: str,
system: str,
energy_range: float = None,
epw_fold: str = 'epw',
epw_params: dict = None, version: float = 5.9):
"""
Initialize the EPWProperties object.
Parameters
----------
file : str
Name of the EPW output file.
system : str
Prefix or name of the system being studied.
energy_range : float, optional
Energy range to consider for the properties.
epw_fold : str, optional
Name of the EPW output folder. Default is 'epw'.
epw_params : dict, optional
Dictionary of EPW calculation parameters.
"""
self.file = file
self.system = system
self.prefix = self.system
self.energy_range = energy_range
self.epw_params = epw_params
self.epw_fold = epw_fold
self.epw_file = f'{self.system}/{self.epw_fold}/{self.file}'
self.epw_ver = version
@property
def dtau(self):
"""
Polaron displacement from EPW.
Returns
-------
array-like or None
Returns the polaron displacement if the EPW file exists;
otherwise, returns None.
"""
try:
return read_dtau(self.epw_file)
except FileNotFoundError:
return None
@property
def psir_plrn(self):
"""
Polaron wavefunction from EPW.
Returns
-------
array-like or None
Returns the polaron wavefunction if the EPW file exists;
otherwise, returns None.
"""
try:
return read_psir_plrn(self.epw_file)
except FileNotFoundError:
return None
@property
def gkk(self):
"""
Electron-phonon matrix elements.
The returned array has dimensions g(ibnd, jbnd, imode, ik, iq).
Returns
-------
array-like or None
Electron-phonon matrix if the EPW file exists;
otherwise, prints a message and returns None.
"""
try:
return get_gkk(self.epw_file, self.epw_ver)
except FileNotFoundError:
print('EPW file not found')
return None
@property
def ibte_mobility_e(self):
_,ibte_mob = read_mobility(f'{self.epw_file}',float(self.epw_params['temps']))
return(ibte_mob)
@property
def serta_mobility_e(self):
serta_mob,_ = read_mobility(f'{self.epw_file}',float(self.epw_params['temps']))
return(serta_mob)
@property
def ibte_mobility_h(self):
_,ibte_mob = read_mobility(f'{self.epw_file}',float(self.epw_params['temps']),typ='Hole')
return(ibte_mob)
@property
def serta_mobility_h(self):
serta_mob,_ = read_mobility(f'{self.epw_file}',float(self.epw_params['temps']),typ = 'Hole')
return(serta_mob)
@property
def Eform(self):
Ef = read_plrn(f'{self.epw_file}')
return(Ef)
@property
def relaxation_time_e(self):
"""
get electron relaxation time
"""
if (self.epw_ver >= 6.0):
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_taucb_0.fmt'))
else:
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_taucb.fmt'))
@property
def relaxation_time_h(self):
"""
get hole relaxation time
"""
if (self.epw_ver >= 6.0):
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_tau_0.fmt'))
else:
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_tau.fmt'))
@property
def relaxation_time_imp(self):
"""
get ionized impurity relaxation time
"""
if (self.epw_ver >= 6.0):
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_tau_0.fmt'))
else:
return(read_inv_taucb(f'{self.system}/{self.epw_fold}/inv_tau.fmt'))
[docs]
def eps1(self, eps0=1.0, leng=100):
T=self.temp
file = f'{self.system}/{self.epw_fold}'
if('lindabs' in self.epw_params.keys()):
_,eps_real,_,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K.dat',eps0,leng)
if('loptabs' in self.epw_params.keys()):
_,eps_real,_,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K_'+self.epw_params['meshnum']+'.dat',eps0, leng)
return(eps_real)
[docs]
def nr(self, eps0=1, leng=100):
file = f'{self.system}/{self.epw_fold}'
T=self.temp#default_epw_input['temps']
if('lindabs' in self.epw_params.keys()):
nr,eps_real,_,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K.dat',eps0,leng)
if('loptabs' in self.epw_params.keys()):
nr,eps_real,_,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K_'+self.default_epw_input['meshnum']+'.dat',eps0,leng)
return(nr)
[docs]
def eps2(self, eps0=1, leng=100):
file = f'{self.system}/{self.epw_fold}'
T=self.temp#default_epw_input['temps']
if('lindabs' in self.epw_params.keys()):
_,_,eps2,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K.dat', eps0,leng)
if('loptabs' in self.epw_params.keys()):
_,_,eps2,_ = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K_'+self.default_epw_input['meshnum']+'.dat',eps,leng)
return(eps2)
[docs]
def omega(self, eps0=1, leng=100):
file = f'{self.system}/{self.epw_fold}'
T=self.temp#default_epw_input['temps']
if('lindabs' in self.epw_params.keys()):
_,_,_,omega = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K.dat',eps0,leng)
if('loptabs' in self.epw_params.keys()):
_,_,_,omega = kk.calc_epr_real(f'{file}/epsilon2_indabs_'+str(T)+'K_'+self.default_epw_input['meshnum']+'.dat',eps0,leng)
return(omega)
[docs]
def optical_properties(self, eps: int = 1, leng: int = None, Temps= None):
"""
Extract and return optical properties from EPW output.
This method parses the dielectric response (ε₂) data obtained
from EPW calculations. It reads the direct, indirect, and
quasi-degenerate perturbation theory (QDPT) contributions to
the imaginary part of the dielectric function from the
`epsilon2` output files.
Parameters
----------
eps : int, optional
Index specifying which dielectric tensor component to extract
(default: ``1``).
leng : int, optional
Number of frequency points to read. If ``None``, uses the
full data length (default: ``None``).
Temps : list or None, optional
List of temperatures at which the optical properties are computed.
If not provided, attempts to read from ``self.epw_params['temps']``.
Returns
-------
dict
Dictionary containing dielectric data:
- ``'direct'`` : dict
Direct electronic transitions contribution.
- ``'indirect'`` : dict
Phonon-assisted (indirect) transitions contribution.
- ``'qdpt'`` : dict
QDPT-corrected optical response.
Notes
-----
- The method automatically detects available temperatures from
EPW input parameters.
- Currently, only ε₂ (imaginary part) is parsed; ε₁ and derived
quantities (e.g., refractive index) can be added later.
Example
-------
>>> epw = EPWProperties(file='epw2.out', system='Si', epw_fold='epw')
>>> data = epw.optical_properties()
>>> direct_eps2 = data['direct']
"""
file = f'{self.system}/{self.epw_fold}'
if Temps is None:
T = self.epw_params['temps']
print('Available temps:', T)
try:
arr = [float(x) for x in T.split()]
except AttributeError:
arr = np.array([T,0])
TT = arr[0]
direct, indirect, qdpt = read_epsilon2_files(f'{file}')
for key in direct.keys():
raw_eps2 = direct[key]
if leng is None:
leng = len(raw_eps2[:,0])
#omega = self.omega(leng = leng)
#eps1 = self.eps1(leng = leng)
#eps2 = self.eps2(leng = leng)
#nr = self.nr(leng = leng)
return {'direct':direct,
'indirect':indirect,
'qdpt':qdpt}
[docs]
def wannier_files(self):
"""
List all Wannier cube files in the EPW output folder.
This method scans the EPW output folder for files with the `.cube`
extension, which are typically Wannier-related outputs for visualization.
Returns
-------
list of str
List of cube filenames present in the folder.
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> cube_files = props.wannier_files()
Present Wannier files are:
graphene_epw1.cube
graphene_epw2.cube
"""
folder = f'{self.system}/{self.epw_fold}'
cube_files = [f for f in os.listdir(folder) if f.endswith(".cube")]
print('Present Wannier files are:')
for files in cube_files:
print(files)
return cube_files
[docs]
def polaron(self, psir_file='psir_plrn.xsf', dtau_file='dtau_disp.plrn'):
"""
Retrieve polaron properties from EPW output files.
Reads the polaron wavefunction and displacement from EPW-generated
files and returns them in a dictionary along with formation energy.
Parameters
----------
psir_file : str, optional
Filename for the polaron wavefunction. Default is 'psir_plrn.xsf'.
dtau_file : str, optional
Filename for the polaron displacement. Default is 'dtau_disp.plrn'.
Returns
-------
dict
Dictionary containing:
- 'dtau': polaronic displacement
- 'psir_plrn': polaron wavefunction
- 'Eform': polaron formation energy
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> polaron_data = props.polaron()
>>> print(polaron_data['dtau'])
"""
self.epw_file = f'{self.system}/{self.epw_fold}/{psir_file}'
psir_plrn = self.psir_plrn
self.epw_file = f'{self.system}/{self.epw_fold}/{dtau_file}'
dtau = self.dtau
self.epw_file = f'{self.system}/{self.epw_fold}/{self.file}'
return {'dtau': dtau,
'psir_plrn': psir_plrn,
'Eform': self.Eform}
[docs]
def transport(self, tau=False):
"""
Retrieve transport properties from EPW calculations.
Returns a dictionary of electron and hole mobilities, as well as
relaxation times for both carrier types. Values are available for
both semi-classical (SERTA) and fully iterated (IBTE) Boltzmann transport
approaches.
Parameters
----------
tau : bool, optional
If True, may be used in future to scale properties by relaxation time
(currently not implemented). Default is False.
Returns
-------
dict
Dictionary containing transport properties:
- 'serta_mobility_elec': electron mobility (SERTA)
- 'serta_mobility_hole': hole mobility (SERTA)
- 'ibte_mobility_elec': electron mobility (IBTE)
- 'ibte_mobility_hole': hole mobility (IBTE)
- 'relaxation_time_elec': electron relaxation time
- 'relaxation_time_hole': hole relaxation time
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> transport_data = props.transport()
>>> print(transport_data['serta_mobility_elec'])
"""
return {'serta_mobility_elec': self.serta_mobility_e,
'serta_mobility_hole': self.serta_mobility_h,
'ibte_mobility_elec': self.ibte_mobility_e,
'ibte_mobility_hole': self.ibte_mobility_h,
'relaxation_time_elec': self.relaxation_time_e,
'relaxation_time_hole': self.relaxation_time_h}
[docs]
def elph(self):
"""
Retrieve electron-phonon (el-ph) data from EPW output.
Returns a dictionary containing the electron-phonon matrix and
derived properties such as the number of k-points, q-points, modes, and bands.
Returns
-------
dict
Dictionary with keys:
- 'gkk' : electron-phonon matrix elements g(ibnd, jbnd, imode, ik, iq)
- 'nk' : number of k-points
- 'nq' : number of q-points
- 'nmode' : number of phonon modes
- 'nband' : number of electronic bands
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> elph_data = props.elph()
>>> print(elph_data['nk'])
"""
return {'gkk': self.gkk,
'nk': len(self.gkk[0, 0, :, 0, 0]),
'nq': len(self.gkk[0, 0, 0, :, 0]),
'nmode': len(self.gkk[0, 0, 0, 0, :]),
'nband': len(self.gkk[:, 0, 0, 0, 0])
}
[docs]
def plot_polaron(self, psir_file='psir_plrn.xsf', view=None):
"""
Plot the polaron wavefunction using Matplotlib.
Parameters
----------
psir_file : str, optional
Filename of the polaron wavefunction. Default is 'psir_plrn.xsf'.
view : dict, optional
Dictionary of plotting parameters. Can include 'iso_level' to set
the isosurface threshold. Default is {}.
Returns
-------
matplotlib.figure.Figure
Matplotlib figure object showing the polaron isosurface.
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> fig = props.plot_polaron(iso_level=0.01)
"""
if view is None:
view = {}
Data = read_psir_plrn(f'{self.system}/{self.epw_fold}/{psir_file}')
iso_level = 0.005
if 'iso_level' in view.keys():
iso_level = view['iso_level']
return plot_pl.plot_psir_plrn_matplotlib(
Data, view=view, iso_level=iso_level*np.nanmax(Data['Density'])
)
[docs]
def plot_wannier(self, wannier_file=None, view=None, bond_length=5.0):
"""
Plot Wannier functions (cube files) using Matplotlib.
Parameters
----------
wannier_file : str, optional
Filename of the cube file to plot. If None, the first available
Wannier cube file in the folder is used.
view : dict, optional
Dictionary of plotting parameters (e.g., camera angles). Default is {}.
bond_length : float, optional
Bond length cutoff used to read the cube file. Default is 5.0 Å.
Returns
-------
matplotlib.figure.Figure
Matplotlib figure object showing the Wannier isosurfaces.
Example
-------
>>> props = EPWProperties(file='epw1.out', system='graphene')
>>> fig = props.plot_wannier(bond_length=4.0)
"""
if view is None:
view = {}
if wannier_file is None:
wannier_file = self.wannier_files()[0]
print(f'Plotting {wannier_file} since no file was provided')
Data = read_cube_file(f'{self.system}/{self.epw_fold}/{wannier_file}', bond_length)
return plot_wann.plot_isosurfaces_matplotlib(Data, view=view)
[docs]
def read_epsilon2_files(folder_path):
"""
Reads all epsilon2 files in the given folder and categorizes them into:
- direct: 'dirabs'
- indirect: 'indabs' without number suffix
- qdpt: 'indabs' with number suffix (e.g., _12, _5, etc.)
Returns:
direct: dict of {filename: data}
indirect: dict of {filename: data}
qdpt: dict of {filename: data}
"""
direct = {}
indirect = {}
qdpt = {}
dirnum = 0
indnum = 0
qdptnum = 0
for fname in os.listdir(folder_path):
if fname.startswith("epsilon2") and fname.endswith(".dat"):
path = os.path.join(folder_path, fname)
try:
data = np.loadtxt(path)
except Exception as e:
print(f"Error reading {fname}: {e}")
continue
if "dirabs" in fname:
direct[f'{dirnum}'] = data
dirnum +=1
elif "indabs" in fname:
# QDPT files have a number suffix before '.dat'
if re.search(r'_?\d+\.dat$', fname):
qdpt[f'{qdptnum}'] = data
qdptnum +=1
else:
indirect[f'{indnum}'] = data
indnum +=1
return direct, indirect, qdpt
[docs]
def read_inv_taucb(filein):
"""
reads inverse taucb
"""
inv_tau=np.loadtxt(f'{filein}')
inv_tau[:,3]=inv_tau[:,3]*13.6
inv_tau[:,4]=inv_tau[:,4]*20670.6944033
return(inv_tau)
[docs]
def read_plrn(file1):
with open(file1, 'r') as f:
for line in f:
if( len(line.split())>0 ):
if( line.split()[0] == 'Formation'):
Eform=float(line.split()[3])
return(Eform)
[docs]
def sig2eig(filename, outfile, eigs=11):
Data = []
eig = []
with open (f'{filename}','r') as f:
for line in f:
Data.append(line.split()[-1])
eig.append(int(line.split()[0]))
t = 0
eigs = max(eig)
with open(f'{outfile}','w') as f:
f.write(' \n')
f.write(' \n')
for data in Data:
if (t == eigs):
f.write(' \n')
t = 0
f.write(f'{data}\n')
t += 1
return(Data)
[docs]
def get_connections(x, y, z, bond_leng=3.5, min_leng=0.4):
connections = list()
distance = []
for i in range(len(x)):
for j in range(len(x)):
dist = np.sqrt((x[i]-x[j])**2+(y[i]-y[j])**2+(z[i]-z[j])**2)
if ((dist < bond_leng) & (dist > min_leng)):
connections.append((i, j))
distance.append(dist)
return(connections, distance)
[docs]
def read_dtau(filein):
"""
reads the dtau
"""
t = 0
data=[]
x=[]
y=[]
z=[]
u=[]
v=[]
w=[]
ilength = np.zeros((3,1),dtype=int)
nline = 1000
mat = []
with open(filein, 'r') as f:
for line in f:
if (t == 10):
nline = int(line.split()[0])
print(nline)
if ((t>10) & (t<=10+nline)):
mat.append(int(line.split()[0]))
x.append(float(line.split()[1]))
y.append(float(line.split()[2]))
z.append(float(line.split()[3]))
u.append(float(line.split()[4]))
v.append(float(line.split()[5]))
w.append(float(line.split()[6]))
t +=1
x = np.array(x)
y = np.array(y)
z = np.array(z)
u = np.array(u)
v = np.array(v)
w = np.array(w)
mat = np.array(mat)
connections,_ = get_connections(x,y,z)
connections = np.array(connections)
return(x,y,z,u,v,w,connections,mat)
[docs]
def read_psir_plrn(filein, bond_leng=2.5):
"""
reads the psir_plrn
"""
t = 0
data=[]
x=[]
y=[]
z=[]
u=[]
v=[]
w=[]
ilength = np.zeros((3,1),dtype=int)
lattice_vec = np.zeros((3,3),dtype = float)
nline = 1000
mat = []
grid = list()
Density=[]
base =1e12
base2=2e12
mat = []
with open(filein, 'r') as f:
for line in f:
if ((t > 5) & (t<9)):
print(line.split())
lattice_vec[t-6,:] = np.array(line.split()[:]).astype(float)
if (t == 10):
nline = int(line.split()[0])
print(nline)
if ((t>10) & (t<=10+nline)):
# print(line.split()[0])
mat.append(int(line.split()[0]))
x.append(float(line.split()[1]))
y.append(float(line.split()[2]))
z.append(float(line.split()[3]))
u.append(float(line.split()[4]))
v.append(float(line.split()[5]))
w.append(float(line.split()[6]))
if (t == 10+nline+6):
base = 10+nline+10
for d in line.split():
grid.append(int(d))
base2 = grid[0]*grid[1]*grid[2]
if ((t > base) & (t < base+base2-1)):
try:
for data in line.split():
Density.append(float(data))
except ValueError:
break
#dd +=1
t +=1
Dense = np.zeros((grid[0],grid[1],grid[2]),dtype = float)
x = np.array(x)
y = np.array(y)
z = np.array(z)
u = np.array(u)
v = np.array(v)
w = np.array(w)
mat = np.array(mat)
Grid = np.zeros((grid[0],grid[1],grid[2],3),dtype = float)
xx,yy,zz = np.mgrid[0:1:int(grid[0])*1j,0:1:grid[1]*1j,0:1:grid[2]*1j]
t = 0
pts = []
#print(np.shape(Dense))
#print(len(Density))
for i in range(grid[0]):
for j in range(grid[1]):
for k in range(grid[2]):
Grid[i,j,k,:] = lattice_vec[0,:]*xx[i,j,k]+lattice_vec[1,:]*yy[i,j,k]+lattice_vec[2,:]*zz[i,j,k]
Dense[i,j,k] = Density[t]
pts.append(Grid[i,j,k,:])
t +=1
connections,_ = get_connections(x,y,z, bond_leng = bond_leng)
Data = {'mat':mat,
'x':x,
'y':y,
'z':z,
'u':u,
'v':v,
'w':w,
'Dense':Dense,
'Grid':Grid,
'pts':pts,
'Density':Density,
'connections':np.array(connections)}
return(Data)
[docs]
def read_wfc(filein):
"""
reads the psir_plrn
"""
t = 0
data=[]
x=[]
y=[]
z=[]
u=[]
v=[]
w=[]
ilength = np.zeros((3,1),dtype=int)
nline = 1000
mat = []
with open(filein, 'r') as f:
for line in f:
if (t == 10):
nline = int(line.split()[0])
print(nline)
if((t>10) & (t<=10+nline)):
mat.append(int(line.split()[0]))
x.append(float(line.split()[1]))
y.append(float(line.split()[2]))
z.append(float(line.split()[3]))
u.append(float(line.split()[4]))
v.append(float(line.split()[5]))
w.append(float(line.split()[6]))
t +=1
x = np.array(x)
y = np.array(y)
z = np.array(z)
u = np.array(u)
v = np.array(v)
w = np.array(w)
mat = np.array(mat)
connections,_ = get_connections(x,y,z)
return(x,y,z,u,v,w,np.array(connections),mat)
[docs]
def read_cube_file(filein, bond_leng=3.5):
"""
Read a .cube file and return the scalar field data, grid dimensions,
atomic positions, origin, and axis vectors.
"""
Data = {}
with open(filein, 'r') as f:
lines = f.readlines()
# Skip the first two comment lines
i = 2
num_atoms, origin_x, origin_y, origin_z = map(float, lines[i].split()) # Number of atoms and origin
num_atoms = int(num_atoms)
i += 1
# Read the grid dimensions and voxel information
grid_info = []
for _ in range(3):
grid_info.append(list(map(float, lines[i].split())))
i += 1
# Extract grid dimensions and axis vectors
n_points_x, n_points_y, n_points_z = abs(int(grid_info[0][0])), abs(int(grid_info[1][0])), abs(int(grid_info[2][0]))
spacing_x, spacing_y, spacing_z = grid_info[0][1], grid_info[1][1], grid_info[2][1]
axis_x = np.array(grid_info[0][1:4]) # Axis vector for X
axis_y = np.array(grid_info[1][1:4]) # Axis vector for Y
axis_z = np.array(grid_info[2][1:4]) # Axis vector for Z
# Atoms data
atomic_positions = []
for _ in range(num_atoms):
parts = lines[i].split()
atom_number = int(parts[0]) # Atom number
charge = float(parts[1]) # Atom charge
x, y, z = map(float, parts[2:]) # Atomic position
atomic_positions.append((atom_number, charge, x, y, z))
i += 1
# Bonds
connections,_ = get_connections(np.array(atomic_positions)[:,2],np.array(atomic_positions)[:,3],np.array(atomic_positions)[:,4], bond_leng)
connections = np.array(connections)
# Read the scalar field data
scalar_data = []
expected_num_values = n_points_x * n_points_y * n_points_z
# Collect scalar data in a list of lists
while i < len(lines):
line = lines[i].split()
if len(line) > 0:
scalar_data.append(line)
i += 1
# Flatten the scalar data and ensure correct number of values
scalar_data = [float(value) for sublist in scalar_data for value in sublist]
# Check if the scalar data contains the expected number of values
scalar_data_size = len(scalar_data)
if scalar_data_size != expected_num_values:
print(f"Warning: Expected {expected_num_values} scalar data points, but found {scalar_data_size}.")
# Pad the data with NaN if it is smaller than expected
if scalar_data_size < expected_num_values:
padding = [np.nan] * (expected_num_values - scalar_data_size)
scalar_data.extend(padding)
# Truncate if there are more values than expected
elif scalar_data_size > expected_num_values:
scalar_data = scalar_data[:expected_num_values]
# Reshape the scalar data into a 3D grid
scalar_data = np.array(scalar_data).reshape((n_points_z, n_points_y, n_points_x))
# Return all extracted data
Data = {'scalar_data': scalar_data,
'grid_shape': (n_points_x, n_points_y, n_points_z),
'spacing': (spacing_x, spacing_y, spacing_z),
'origin': (origin_x, origin_y, origin_z),
'atomic_positions': atomic_positions,
'connections': connections,
'axis_x': axis_x, 'axis_y': axis_y, 'axis_z': axis_z}
return Data
[docs]
def read_gkk(filn, nbnd, nq, nk, nmode, version=5.9):
""" Read the g matrix """
G=np.zeros((nbnd,nbnd,nmode,nk,nq),dtype=float)
t=0
iq=0
with open(filn,'r') as f:
k=0
q=0
flag = 0
read=0
for line in f:
if ('ik ' in line):
k +=1
if (k == nk+1):
k = 1
if ('iq ' in line):
q +=1
if (q == nq+1):
q = 1
if (('-----' in line) & (k > 0)):
flag = 1
if (('-----' in line) & (flag >= 2)):
read = 0
flag = 0
if ((flag >= 2) & (len(line.split())>0)):
#print(line.split())
try:
ibnd=int(line.split()[0])-1
jbnd=int(line.split()[1])-1
imode=int(line.split()[2])-1
if (version >= 6.0):
G[ibnd,jbnd,imode,k-1,q-1] = float(line.split()[-3])
else:
G[ibnd,jbnd,imode,k-1,q-1] = float(line.split()[-1])
flag +=1
except ValueError:
flag = 0
if (flag >= 1):
flag +=1
f.close()
return(G)
[docs]
def get_gkk(epw_file, version=5.9):
nq = []
nk = []
q_arr = []
k_arr = []
g_mat = []
band1=[]
band2=[]
mode = []
kk = []
qq = []
t = 0
iq = 0
ik = 0
with open(epw_file, 'r') as f:
for line in f:
if ('iq' in line):
Arr = line.split()[2:]
t = 1
if (Arr in q_arr):
pass
else:
q_arr.append(Arr)
iq +=1
if ('ik' in line):
Arr = line.split()[2:]
if (Arr in k_arr):
pass
else:
k_arr.append(Arr)
ik +=1
t = t*2
if ((t > 16) & ('-----------' in line)):
t = 0
if (t > 16):
kk.append(ik)
qq.append(iq)
band1.append(int(line.split()[0]))
band2.append(int(line.split()[1]))
mode.append(int(line.split()[2]))
g_mat.append(float(line.split()[-1]))
f.close()
g = read_gkk(epw_file, max(band1), len(q_arr), len(k_arr), max(mode), version=version)
return(g)
[docs]
def read_ibndmin(file):
""" Reads the minimum bands used"""
with open(file, 'r') as f:
for line in f:
if('ibndmin' in line):
ibndmin = int(line.split()[2])
f.close()
return(ibndmin)
[docs]
def read_ibndmax(file):
""" Reads the maximum number of bands """
with open(file, 'r') as f:
for line in f:
if('ibndmax' in line):
ibndmax = int(line.split()[2])
f.close()
return(ibndmax)
[docs]
def read_mobility(file, T, typ='Elec'):
read_mob_serta = False
read_mob_ibte = False
read_serta = False
read_ibte = False
serta_mob=np.zeros((3,3),dtype=float)
ibte_mob=np.zeros((3,3),dtype=float)
temp_found = False
t = 0
with open(f'{file}', 'r') as f:
for line in f:
if ((len(line.split()) != 0) & (read_mob_serta == True)):
try:
if((float(line.split()[0]) == T) & (t == 0)):
temp_found = True
if(temp_found == True):
serta_mob[t,:] = np.array(line.split()).astype(float)[-3:]
t +=1
if(t == 3):
read_mob_serta = False
read_serta = False
t = 0
temp_found = False
except ValueError:
pass
if ((len(line.split()) != 0) & (read_mob_ibte == True)):
try:
if((float(line.split()[0]) == T) & (t == 0)):
temp_found = True
if(temp_found == True):
ibte_mob[t,:] = np.array(line.split()).astype(float)[-3:]
t +=1
if(t == 3):
read_mob_ibte = False
t = 0
temp_found = False
except ValueError:
pass
if ('SERTA' in line):
read_serta = True
if ('Iteration number' in line):
read_ibte = True
if (f'Drift {typ} mobility' in line):
if(read_serta == True) :
read_mob_serta = True
elif(read_ibte == True):
read_mob_ibte = True
f.close()
return(serta_mob, ibte_mob)
[docs]
def read_grr(filn, nbnd, nq, nmode):
pass
[docs]
def read_dtau_plrn(filename):
"""
For reading dtau.plrn for further postprocessing
"""
file = open(filename, 'r');
control = ( file.readline() ).split();
control = np.array(control, dtype=int);
#
dim = [control[0], control[1], control[2]]
dim.append(control[-1] // 3)
dim.append(3)
#
disp = np.zeros((dim),dtype=complex);
#
for ix in range(control[0]):
for iy in range(control[1]):
for iz in range(control[2]):
for iat in range(dim[-2]):
for icart in range(3):
temp = ( file.readline() ).split()
disp[ix,iy,iz,iat,icart] = float(temp[-2]) + 1j * float(temp[-1])
#
file.close();
return disp
[docs]
def generate_dtau_disp(dtau, direction='x'):
"""
Shift the polaron along the [100] direction
"""
shift_x = np.zeros((np.shape(dtau)),dtype=complex);
#
params = np.shape(dtau)
#
idx = np.zeros((params[0]), dtype=int)
idy = np.zeros((params[0]), dtype=int)
idz = np.zeros((params[0]), dtype=int)
idx[0:params[0] - 1] = np.arange(1, params[0], dtype=int); idx[-1] = 0
idy[0:params[1] - 1] = np.arange(1, params[1], dtype=int); idy[-1] = 0
idz[0:params[2] - 1] = np.arange(1, params[2], dtype=int); idz[-1] = 0
#
for ix in range(params[0]):
for iy in range(params[1]):
for iz in range(params[2]):
for iat in range(params[3]):
for icart in range(3):
if direction == 'x':
shift_x[idx[ix], iy, iz, iat, icart] = dtau[ix, iy, iz, iat, icart]
elif direction == 'y':
shift_x[ix, idy[iy], iz, iat, icart] = dtau[ix, iy, iz, iat, icart]
elif direction == 'z':
shift_x[ix, iy, idz[iz], iat, icart] = dtau[ix, iy, iz, iat, icart]
else:
raise ValueError("Only shift directions x, y, z are implemented.")
#
return shift_x
[docs]
def write_dtau_disp(initial, final, ntau, fileout):
"""
Generate the dtau_disp.plrn_* files
"""
params = np.shape(final)
midpoints = ntau;
scaling = np.linspace(0,1,ntau,endpoint=True)
shift = final - initial
for i in range(midpoints):
file = open(fileout + str(i+1), 'w');
file.write( '{:10d}'.format(params[0]) + '{:10d}'.format(params[1]) + '{:10d}'.format(params[2]) );
file.write( '{:10d}'.format(params[0] * params[1] * params[2]) );
file.write( '{:10d}'.format(params[3] * 3) + '\n');
for ix in range(params[0]):
for iy in range(params[1]):
for iz in range(params[2]):
iunit = ix * params[1] * params[2] + iy * params[2] + iz + 1
for iat in range(params[3]):
for icart in range(3):
ifree = iat * 3 + icart + 1
disp_value = shift[ix, iy, iz, iat, icart] * scaling[i]+\
initial[ix, iy, iz, iat, icart];
file.write( '{:18.10E}'.format(np.real(disp_value)) + '{:18.10E}'.format(np.imag(disp_value)) +'\n' )
file.close();
if __name__=="__main__":
#from plot_g import *
import matplotlib.pyplot as plt
G=get_gkk('notebooks_basic/si/epw/epw1.out')
print(np.shape(G))
#plot_gkk(4,4,G[:,:,:,0,:])
#plt.show()
# print(np.shape(G))
# print(G)
