import numpy as np
[docs]
class PHProperties:
"""
A class for parsing and analyzing phonon (PH) calculation properties.
This class provides access to various physical quantities derived from
Quantum ESPRESSO (QE) `ph.x` output files, such as dielectric constants,
Fermi level, band structures, and lattice vectors. It serves as a post-processing
utility to extract, store, and manipulate phonon-related data.
Parameters
----------
file : str
Path to the phonon output file generated by QE (`ph.out`).
energy_range : float, optional
Energy window (in eV) around the Fermi level for evaluating band-related
quantities. Defaults to ``None``.
Attributes
----------
eps : float
Static dielectric constant of the system.
nbnd_energy : int
Number of bands within the specified energy range.
lattice_vec : ndarray
Lattice vectors extracted from the associated PW calculation.
fermi_level : float
Fermi energy or the highest occupied band energy (in eV).
atomic_positions : ndarray
Atomic positions in Cartesian or fractional coordinates.
band_scf : ndarray
Self-consistent field (SCF) band structure data array.
Examples
--------
>>> ph = PHProperties(file='Si.ph.out', energy_range=5.0)
>>> print(ph.file)
'Si.ph.out'
>>> print(ph.energy_range)
5.0
"""
def __init__(self, file: str,
energy_range: float = None,
version: float = 7.4
):
self.file = file
self.energy_range = energy_range
self.version = version
[docs]
def ph_properties(self):
"""
Returns various phonon related outputs
"""
return({'q_points':self.qpoints,
'output_file':self.file,
'dielectric':self.eps_ph})
[docs]
def plot_pl_vector(self,mode):
"""
performs a quiver plot of pl vectors
"""
pass
@property
def eps_ph(self):
"""
Returns high-frequency dielectric matrix from DFPT
"""
dielectric=get_dielectric_matrix(self.file)
return(dielectric)
@property
def qpoints(self):
"""
Returns q-points for Ph calculation
"""
return(get_qpoints(self.file))
@property
def gkk(self):
"""
Returns the electron–phonon coupling matrix g(ibnd, jbnd, imode, ik, iq)
extracted from a Quantum ESPRESSO `ph.x` output file.
The matrix has dimensions:
(nbnd_i, nbnd_j, nmodes, nkpoints, nqpoints)
where:
- ibnd, jbnd: band indices
- imode: phonon mode index
- ik, iq: k-point and q-point indices
Returns
-------
np.ndarray
5D array of electron–phonon coupling strengths in units consistent with the QE output.
"""
return(get_gkk(self.file))
@property
def dynmat(self):
"""
Returns the dynamical matrix
"""
pass
@property
def polarization_vector(self):
"""
Compute and return the polarization vectors from the dynamical matrix file.
This property reads the `dynmat.axsf` file generated by `dynmat.x`
and extracts the atomic displacement (polarization) vectors corresponding
to phonon modes at the Γ-point or other q-points, depending on the
underlying calculation.
Returns
-------
numpy.ndarray or None
A 3D NumPy array of shape ``(nmode, natoms, 6)`` containing the
polarization vectors and atomic coordinates for each phonon mode.
Returns ``None`` if no polarization data are found.
Notes
-----
Each entry in the last dimension (6 columns) corresponds typically to:
`[x, y, z, dx, dy, dz]`, where `(x, y, z)` are atomic positions and
`(dx, dy, dz)` are displacement components for a given mode.
"""
return(get_polarization_vec(self.dynmat_file, self.atoms, self.natoms))
[docs]
def get_polarization_vec(dynmat_file, atoms, natoms):
"""
Reads dynmat.axsf file to return polarization
vectors
"""
matrix = []
t = 0
with open(dynmat_file,'r') as f:
for line in f:
#print(line.split())
if (line.split()[0] == 'ANIMSTEPS'):
nmode = int(line.split()[1])
if (line.split()[0] == 'PRIMCOORD'):
check1 = True
t = 0
if (line.split()[0] in atoms):
matrix.append(line.split()[1:])
if (len(matrix) > 0):
matrix1 = np.array(matrix).astype(float).reshape(nmode, natoms, 6)
return(matrix1)
else:
return(None)
[docs]
def get_dielectric_matrix(ph_file):
T = 0
i = 0
dielectric = np.zeros((3,3), dtype = float)
with open(ph_file, 'r') as f:
for line in f:
if ('Dielectric constant in cartesian axis' in line):
T = 1
if ((T > 1) & (i < 3) & (len(line.split()) > 0)):
dielectric[i,:] = np.array(line.split()[1:4]).astype(float)
i +=1
T = 2*T
f.close()
return(dielectric)
[docs]
def get_qpoints(ph_file):
T = 0
i = 0
nqs = 0
with open(ph_file, 'r') as f:
for line in f:
if ('Dynamical matrices for' in line):
T = 1
if (T == 2):
nqs = int(line.split()[1])
ph_qpoints = np.zeros((nqs,3),dtype=float)
if ((T > 4) & (i < nqs) & (len(line.split()) > 0)):
ph_qpoints[i,:3] = np.array(line.split()[1:4]).astype(float)
i +=1
T = 2*T
f.close()
return(ph_qpoints)
[docs]
def get_gkk(ph_file):
nq = []
nk = []
q_arr = []
k_arr = []
g_mat = []
band1=[]
band2=[]
mode = []
kk = []
qq = []
t = 0
iq = 0
ik = 0
with open(ph_file, 'r') as f:
for line in f:
if ('q coord' in line):
Arr = line.split()[2:]
t = 1
if (Arr in q_arr):
pass
else:
q_arr.append(Arr)
iq +=1
if ('k coord' 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]))
if (self.version < 7.5):
g_mat.append(float(line.split()[-3]))
else:
g_mat.append(float(line.split()[-3]))
f.close()
g = np.zeros((max(band1),max(band2),max(mode),len(k_arr),len(q_arr)),dtype = float)
t = 0
for i in range(len(g_mat)):
g[band1[i]-1,band2[i]-1,mode[i]-1,kk[i]-1,qq[i]-1] = g_mat[i]
return(g)
if __name__=="__main__":
GKK = get_gkk('notebooks_basic/si/ph/ph_g.out')
for i in range(4):
for j in range(4):
for k in range(6):
print(GKK[i,j,k,0,0])
