Add gpath.py

abipy/core/release.py

@@ -75,6 +75,7 @@
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
     "Development Status :: 4 - Beta",
     "Intended Audience :: Science/Research",
     "License :: OSI Approved :: GNU General Public License v2 (GPLv2)",

abipy/dfpt/phtk.py

@@ -19,7 +19,7 @@
 def get_dyn_mat_eigenvec(phdispl, structure, amu=None, amu_symbol=None) -> np.ndarray:
     """
     Converts the phonon displacements to the orthonormal eigenvectors of the dynamical matrix.
-    Small discrepancies with the original values may be expected due to the different values 
+    Small discrepancies with the original values may be expected due to the different values
     of the atomic masses in abinit and pymatgen.
 
     .. note::
@@ -87,7 +87,7 @@ def match_eigenvectors(v1, v2) -> np.ndarray:
 class NonAnalyticalPh(Has_Structure):
     """
     Phonon data at gamma including non analytical contributions
-    Read from anaddb.nc
+    Usually read from anaddb.nc
     """
 
     def __init__(self, structure, directions, phfreqs, phdispl_cart, amu=None):

abipy/electrons/ebands.py

@@ -354,7 +354,7 @@ def from_file(cls, filepath: str) -> ElectronBands:
         else:
             raise NotImplementedError("ElectronBands can only be initialized from nc files")
 
-        assert new.__class__ == cls
+        assert new.__class__ is cls
         return new
 
     @classmethod
@@ -528,6 +528,7 @@ def __init__(self, structure, kpoints, eigens, fermie, occfacts, nelect, nspinor
         self._occfacts = np.atleast_3d(occfacts)
         assert self._eigens.shape == self._occfacts.shape
         self._linewidths = None
+
         if linewidths is not None:
             self._linewidths = np.reshape(linewidths, self._eigens.shape)
 
@@ -541,8 +542,11 @@ def __init__(self, structure, kpoints, eigens, fermie, occfacts, nelect, nspinor
             self.nband_sk.shape = (self.nsppol, self.nkpt)
 
         self.kpoints = kpoints
-        assert self.nkpt == len(self.kpoints)
-        assert isinstance(self.kpoints, KpointList)
+        if self.nkpt != len(self.kpoints):
+            raise ValueError(f"{self.nkpt=} != {len(self.kpoints)=}")
+
+        if not isinstance(self.kpoints, KpointList):
+            raise TypeError(f"{type(self.kpoints)=} is not an instance of KpointList")
 
         self.smearing = {} if smearing is None else smearing
         self.nelect = float(nelect)
@@ -560,7 +564,7 @@ def __init__(self, structure, kpoints, eigens, fermie, occfacts, nelect, nspinor
 
         if self.smearing and self.smearing.occopt == 1 and self.nsppol == 1 and self.nspden == 1:
             # This is the simplest case as the calculation has been done assuming a non-magnetic semiconductor
-            # so there must be a gap.
+            # so there must be a gap in single-particle theory.
             # On the other hand, in a NSCF run with a k-path it may happen that the HOMO energy
             # taken from the GS SCF run underestimates the real HOMO level.
             # For instance, the GS-SCF may have been done with a shifted k-mesh whereas
@@ -572,7 +576,7 @@ def __init__(self, structure, kpoints, eigens, fermie, occfacts, nelect, nspinor
 
         else:
             #print("This is a wannabe metal")
-            # This flag is true if the system must be metallic in a single particle theory
+            # This flag is true if the system must be metallic in single particle theory
             is_bloch_metal = self.nsppol == 1 and self.nspinor == 1 and self.nelect % 2 != 0
 
             self.is_metal = True
@@ -596,9 +600,11 @@ def structure(self) -> Structure:
 
     @lazy_property
     def _auto_klabels(self):
-        # Find the k-point names in the pymatgen database.
-        # We'll use _auto_klabels to label the point in the matplotlib plot
-        # if klabels are not specified by the user.
+        """
+        Find the k-point names in the pymatgen database.
+        We'll use _auto_klabels to label the point in the matplotlib plot
+        if klabels are not specified by the user.
+        """
 
         _auto_klabels = OrderedDict()
         # If the first or the last k-point are not recognized in findname_in_hsym_stars
@@ -856,20 +862,6 @@ def with_points_along_path(self, frac_bounds=None, knames=None, dist_tol=1e-12):
 
         return dict2namedtuple(ebands=new_ebands, ik_new2prev=ik_new2prev)
 
-    #def select_bands(self, bands, kinds=None):
-    #    """Build new ElectronBands object by selecting bands via band_slice (slice object)."""
-    #    bands = np.array(bands)
-    #    kinds = np.array(kinds) if kinds is not None else np.array(range(self.nkpt))
-    #    # This won't work because I need a KpointList object.
-    #    new_kpoints = self.kpoints[kinds]
-    #    new_eigens = self.eigens[:, kinds, bands].copy()
-    #    new_occfacts = self.occupation[:, kinds, bands].copy()
-    #    new_linewidths = None if not self.linewidths else self.linewidths[:, kinds, bands].copy()
-
-    #    return self.__class__(self.structure, new_kpoints, new_eigens, self.fermie, new_occfacts,
-    #                          self.nelect, self.nspinor, self.nspden,
-    #                          smearing=self.smearing, linewidths=new_linewidths)
-
     @classmethod
     def empty_with_ibz(cls, ngkpt, structure, fermie, nelect, nsppol, nspinor, nspden, mband,
                        shiftk=(0, 0, 0), kptopt=1,
@@ -1909,8 +1901,7 @@ def compare_gauss_edos(self, widths, step=0.1) -> ElectronDosPlotter:
         edos_plotter = ElectronDosPlotter()
         for width in widths:
             edos = self.get_edos(method="gaussian", step=0.1, width=width)
-            label = r"$\sigma = %s$ (eV)" % width
-            edos_plotter.add_edos(label, edos)
+            edos_plotter.add_edos(r"$\sigma = %s$ (eV)" % width, edos)
 
         return edos_plotter
 
@@ -2042,7 +2033,7 @@ def get_ejdos(self, spin, valence, conduction,
                         jdos += fact * gaussian(mesh, width, center=ec-ev)
 
         else:
-            raise NotImplementedError("Method %s is not supported" % str(method))
+            raise NotImplementedError(f"{method=} is not supported")
 
         return Function1D(mesh, jdos)
 
@@ -2130,6 +2121,7 @@ def apply_scissors(self, scissors):
         """
         if self.nsppol == 1 and not isinstance(scissors, Iterable):
             scissors = [scissors]
+
         if self.nsppol == 2 and len(scissors) != 2:
             raise ValueError("Expecting two scissors operators for spin up and down")
 
@@ -2745,7 +2737,6 @@ def plotly_traces(self, fig, e0, rcd=None, spin=None, band=None, showlegend=Fals
                     fig.add_scatter(x=xx, y=yy + w, mode='lines', line=lw_opts, name='',
                                     showlegend=False, fill='tonexty', row=ply_row, col=ply_col)
 
-
     def _make_ticks_and_labels(self, klabels):
         """Return ticks and labels from the mapping qlabels."""
         if klabels is not None:
@@ -3189,16 +3180,14 @@ def to_bxsf(self, filepath):
         suitable for the visualization of isosurfaces with xcrysden_ (xcrysden --bxsf FILE).
         Require k-points in IBZ and gamma-centered k-mesh.
         """
-        err_msg = self.isnot_ibz_sampling(require_gamma_centered=True)
-        if err_msg:
+        if err_msg := self.isnot_ibz_sampling(require_gamma_centered=True):
             raise ValueError(err_msg)
 
         self.get_ebands3d().to_bxsf(filepath)
 
     #@memoized_method(maxsize=5, typed=False)
     def get_ebands3d(self):
-        err_msg = self.isnot_ibz_sampling()
-        if err_msg:
+        if err_msg := self.isnot_ibz_sampling():
             raise ValueError(err_msg)
 
         return ElectronBands3D(self.structure, self.kpoints, self.has_timrev, self.eigens, self.fermie)
@@ -3222,8 +3211,7 @@ def derivatives(self, spin, band, order=1, acc=4):
             #ebranch = 0.5 * units.Ha_to_eV * np.array([(k.norm * units.bohr_to_ang)**2 for k in self.kpoints])
 
             # Compute derivatives by finite differences.
-            ders_onlines = self.kpoints.finite_diff(ebranch, order=order, acc=acc)
-            return ders_onlines
+            return self.kpoints.finite_diff(ebranch, order=order, acc=acc)
 
         else:
             raise NotImplementedError("Derivatives on homogeneous k-meshes are not supported yet")
@@ -3348,8 +3336,7 @@ def interpolate(self, lpratio=5, knames=None, vertices_names=None, line_density=
                     interpolator: |SkwInterpolator| object.
         """
         # Get symmetries from abinit spacegroup (read from file).
-        abispg = self.structure.abi_spacegroup
-        if abispg is None:
+        if (abispg := self.structure.abi_spacegroup) is None:
             abispg = self.structure.spgset_abi_spacegroup(has_timerev=self.has_timrev)
 
         fm_symrel = [s for (s, afm) in zip(abispg.symrel, abispg.symafm) if afm == 1]
@@ -3439,8 +3426,7 @@ def dataframe_from_ebands(ebands_objects, index=None, with_spglib=True) -> pd.Da
     # Use OrderedDict to have columns ordered nicely.
     odict_list = [(ebands.get_dict4pandas(with_spglib=with_spglib)) for ebands in ebands_list]
 
-    return pd.DataFrame(odict_list, index=index,
-                        columns=list(odict_list[0].keys()) if odict_list else None)
+    return pd.DataFrame(odict_list, index=index, columns=list(odict_list[0].keys()) if odict_list else None)
 
 
 class ElectronBandsPlotter(NotebookWriter):
@@ -3528,8 +3514,7 @@ def get_ebands_frame(self, with_spglib=True) -> pd.DataFrame:
         Build a |pandas-DataFrame| with the most important results available in the band structures.
         Useful to analyze band-gaps.
         """
-        return dataframe_from_ebands(list(self.ebands_dict.values()),
-                                     index=list(self.ebands_dict.keys()), with_spglib=with_spglib)
+        return dataframe_from_ebands(list(self.ebands_dict.values()), index=list(self.ebands_dict.keys()), with_spglib=with_spglib)
 
     @property
     def ebands_list(self) -> List[ElectronBands]:
@@ -5950,7 +5935,6 @@ def gridplot_with_hue(self, hue, ylims=None, fontsize=8,
 
 
 
-
 def find_yaml_section_in_lines(lines, tag):
 
     magic = f"--- !{tag}"

abipy/eph/gkq.py

@@ -186,8 +186,11 @@ def read_all_gkq(self, mode: str = "phonon") -> np.ndarray:
 
         return gkq_nu
 
-    def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
+    def get_absg_kpoint(self, kpoint, eps_mev: float=0.01) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
         """
+        Args:
+            kpoint: |Kpoint| object or list/tuple with reduced coordinates or integer with the index
+            eps_mev: Tolerance in mev used to detect degeneracies
         """
         if duck.is_intlike(kpoint):
             ik = kpoint
@@ -196,22 +199,26 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
             kpoint = Kpoint.as_kpoint(kpoint, self.structure.reciprocal_lattice)
             ik = self.kpoints.index(kpoint)
 
-        # (nsppol, nkpt, 3*natom, mband, mband) real array.
-        absg = np.abs(self.read_all_gkq(mode="phonon")) * abu.Ha_meV
-        absgk = absg[:,ik].copy()
-        absg_unsym = absg[:,ik].copy()
-
-        # Average over phonons.
-        eps_ha = 0.01 / abu.Ha_meV
+        eps_ha = eps_mev / abu.Ha_meV
         eps_ev = eps_ha * abu.Ha_eV
 
+        nsppol = self.ebands.nsppol
         natom3 = len(self.structure) * 3
+        nb = self.ebands.nband
+
         phfreqs_ha = self.phfreqs_ha
-        nband = self.ebands.nband
+        eigens_k = self.ebands.eigens
+        eigens_kq = self.eigens_kq
+
+        # (nsppol, nkpt, 3*natom, mband, mband) real array.
+        absg = np.abs(self.read_all_gkq(mode="phonon")) * abu.Ha_meV
+        absgk = absg[:,ik].copy()
+        absg_unsym = absg[:,ik].copy()
         absg_sym = np.zeros_like(absgk)
 
-        for spin in range(self.ebands.nsppol):
-            g2_mn = np.zeros((nband,nband), dtype=float)
+        # Average over phonons.
+        for spin in range(nsppol):
+            g2_mn = np.zeros((nb, nb), dtype=float)
             for nu in range(natom3):
                 w_1 = phfreqs_ha[nu]
                 g2_mn[:], nn = 0.0, 0
@@ -225,26 +232,27 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
         # Average over k electrons.
         absg = absg_sym.copy()
         g2_nu = np.zeros((natom3), dtype=float)
-        for spin in range(self.ebands.nsppol):
-            for jbnd in range(nband):
-                for ibnd in range(nband):
-                    w_1 = self.ebands.eigens[spin, ik, ibnd]
+        for spin in range(nsppol):
+            for jbnd in range(nb):
+                for ibnd in range(nb):
+                    w_1 = eigens_k[spin, ik, ibnd]
                     g2_nu[:], nn = 0.0, 0
-                    for pbnd in range(nband):
-                        w_2 = self.ebands.eigens[spin, ik, pbnd]
+                    for pbnd in range(nb):
+                        w_2 = eigens_k[spin, ik, pbnd]
                         if abs(w_2 - w_1) >= eps_ev: continue
                         nn += 1
+                        # MG FIXME: Why absgk and not absg here as done below for k+q?
                         g2_nu += absgk[spin,:,jbnd,pbnd] ** 2
                     absg_sym[spin,:,jbnd,ibnd] = np.sqrt(g2_nu / nn)
+
         # Average over k+q electrons.
-        eigens_kq = self.eigens_kq
         absgk = absg_sym.copy()
-        for spin in range(self.ebands.nsppol):
-            for ibnd in range(nband):
-                for jbnd in range(nband):
+        for spin in range(nsppol):
+            for ibnd in range(nb):
+                for jbnd in range(nb):
                     w_1 = eigens_kq[spin, ik, jbnd]
                     g2_nu[:], nn = 0.0, 0
-                    for pbnd in range(nband):
+                    for pbnd in range(nb):
                         w_2 = eigens_kq[spin, ik, pbnd]
                         if abs(w_2 - w_1) >= eps_ev: continue
                         nn += 1