From f4d0a33b21c624e08398bd403206b20a6665d620 Mon Sep 17 00:00:00 2001
From: Matteo Giantomassi <gmatteo@gmail.com>
Date: Sat, 1 Feb 2025 11:46:38 +0100
Subject: [PATCH] More f-strings

---
 abipy/abio/inputs.py                   | 60 ++++++++++++--------------
 abipy/abio/robots.py                   | 12 +++---
 abipy/dfpt/gruneisen.py                | 22 +++++-----
 abipy/dfpt/phtk.py                     |  6 ++-
 abipy/dfpt/raman.py                    | 11 ++---
 abipy/dfpt/vzsisa.py                   | 17 ++------
 abipy/dynamics/analyzer.py             | 20 ++++-----
 abipy/dynamics/cpx.py                  |  3 +-
 abipy/dynamics/hist.py                 | 14 +++---
 abipy/electrons/arpes.py               | 23 ++++------
 abipy/electrons/bse.py                 | 56 +++++++++---------------
 abipy/electrons/fatbands.py            | 52 +++++++++++-----------
 abipy/electrons/fold2bloch.py          | 33 +++++++-------
 abipy/electrons/gwr.py                 |  8 ++--
 abipy/electrons/lruj.py                |  6 ---
 abipy/electrons/psps.py                | 31 ++++++-------
 abipy/electrons/scissors.py            | 21 +++++----
 abipy/electrons/tests/test_psps.py     | 12 ++++--
 abipy/eph/tests/test_varpeq.py         | 59 +++++++++++++++----------
 abipy/eph/varpeq.py                    |  5 ++-
 abipy/examples/flows/run_conducwork.py |  2 +-
 abipy/examples/flows/run_qha_2d.py     |  8 ++--
 abipy/examples/plot/plot_qha_2d.py     |  8 ++--
 abipy/examples/plot/plot_qha_vzsisa.py | 29 ++++++-------
 abipy/flowtk/vzsisa.py                 |  2 +-
 abipy/iotools/cube.py                  |  2 +-
 abipy/iotools/visualizer.py            | 14 +++---
 abipy/tools/plotting.py                | 25 +++--------
 abipy/wannier90/abiwan.py              | 10 ++---
 abipy/wannier90/wout.py                |  3 +-
 docs/api/abipy.tools.rst               |  8 ++++
 31 files changed, 279 insertions(+), 303 deletions(-)

diff --git a/abipy/abio/inputs.py b/abipy/abio/inputs.py
index 7924e1934..a059c0afc 100644
--- a/abipy/abio/inputs.py
+++ b/abipy/abio/inputs.py
@@ -119,8 +119,6 @@
 #        raise ValueError("Don't know how to reallocate variable %s" % str(name))
 
 
-
-
 class AbstractInput(MutableMapping, metaclass=abc.ABCMeta):
     """
     Abstract class defining the methods that must be implemented by Input classes.
@@ -421,7 +419,8 @@ def __init__(self, structure, pseudos,
 
         if pseudo_dir:
             pseudo_dir = os.path.abspath(pseudo_dir)
-            if not os.path.exists(pseudo_dir): raise self.Error("Directory `%s` does not exist" % pseudo_dir)
+            if not os.path.exists(pseudo_dir):
+                raise self.Error(f"Directory `{pseudo_dir}` does not exist")
             pseudos = [os.path.join(pseudo_dir, p) for p in list_strings(pseudos)]
 
         try:
@@ -610,7 +609,7 @@ def _check_varname(self, key: str) -> None:
     input.set_spell_check(False)
 
 to disable spell checking. Perhaps the internal database is not in synch
-with the Abinit version you are using. Please contact the AbiPy developers.""" % key)
+with the Abinit version you are using? Please contact the AbiPy developers.""" % key)
 
     @property
     def runlevel(self):
@@ -850,7 +849,7 @@ def escape(text):
                     app(str(InputVariable(vname, value)))
 
         else:
-            raise ValueError("Unsupported value for sortmode %s" % str(sortmode))
+            raise ValueError(f"Unsupported value for {sortmode=}")
 
         def to_html(string):
             string = string.replace("\n", "<br>")
@@ -1804,8 +1803,8 @@ def make_edos_input(self, ngkpt, shiftk=(0, 0, 0), tolwfr=1e-20, nscf_nband=None
                             prtwf=-1,
                             comment="Input file for electron DOS calculation from a GS SCF input (NSCF on kmesh)",
                             )
-        edos_input.set_kmesh(ngkpt, shiftk)
 
+        edos_input.set_kmesh(ngkpt, shiftk)
         edos_input.set_vars(**extra_abivars)
 
         return edos_input
@@ -1952,7 +1951,8 @@ def make_ph_inputs_qpoint(self, qpt, tolerance=None,
 
         return ph_inputs
 
-    def make_ddkpert_input(self, perturbation, kptopt=2, only_vk=False, use_symmetries=False, tolerance=None, manager=None) -> AbinitInput:
+    def make_ddkpert_input(self, perturbation, kptopt=2, only_vk=False,
+                           use_symmetries=False, tolerance=None, manager=None) -> AbinitInput:
         """
         Returns |AbinitInput| for the calculation of an electric field perturbation.
         This function should be called with an input that represents a GS run and
@@ -2657,16 +2657,19 @@ def make_eph_transport_input(self, ddb_ngqpt, sigma_erange, tmesh, eph_ngqpt_fin
 
     #    return new
 
-    def make_gwr_qprange_input(self, gwr_ntau, nband, ecuteps, gw_qprange=0, gwr_task=GWR_TASK.G0W0,
+    def make_gwr_qprange_input(self, gwr_ntau: int, nband: int, ecuteps: float,
+                               gw_qprange: int = 0,
+                               gwr_task=GWR_TASK.G0W0,
                                **kwargs) -> AbinitInput:
         """
         Build and return an input file to compute QP corrections with the GWR code.
 
         Args:
             gwr_ntau: Number of minimax points.
-            nband: Number of bands in Green's function
-            ecuteps: Cutoff energy for chi0
-            gw_qprange = 0 to compute the QP corrections only for the fundamental and the direct gap.
+            nband: Number of bands in Green's function.
+            ecuteps: Cutoff energy for chi0 in Ha.
+            gw_qprange: 0 to compute the QP corrections only for the fundamental and the direct gap.
+                For other values see Abinit docs.
             gwr_task: String defining the GWR task
         """
         new = self.new_with_vars(
@@ -2725,7 +2728,7 @@ def abiget_dryrun_task(self, workdir, manager, **extra_vars):
 
     def abiget_spacegroup(self, tolsym=None, retdict=False, workdir=None, manager=None, verbose=0):
         """
-        This function invokes Abinit to get the space group (as detected by Abinit, not by spglib)
+        This function invokes Abinit to get the space group as detected by Abinit, and not by spglib.
         It should be called with an input file that contains all the mandatory variables required by ABINIT.
 
         Args:
@@ -2765,18 +2768,6 @@ def abiget_dims_spginfo(self, workdir=None, manager=None, verbose=0) -> tuple[di
         except Exception as exc:
             self._handle_task_exception(task, exc)
 
-    #def get_qptopt(self) -> int:
-    #    """Helper function that returns"""
-    #    if "qptopt" in self:
-    #        return self["qptopt"]
-
-    #    kptopt = self.get("kptopt", 1)
-    #    nspinor = self.get("nspinor", 1)
-    #    nsppol = self.get("nsppol", 1)
-    #    nspden = self.get("nspden", 1)
-    #    qptopt = 1
-    #    return qptopt
-
     def abiget_ibz(self, ngkpt=None, shiftk=None, kptopt=None, workdir=None, manager=None, verbose=0):
         """
         This function computes the list of points in the IBZ with the corresponding weights.
@@ -2860,7 +2851,7 @@ def abiget_scr_ibz(self, ngkpt=None, shiftk=None, kptopt=None, workdir=None, man
         if verbose:
             print("Computing qptdms with input:\n", str(inp))
 
-        # Build a Task to run Abinit in a shell subprocess
+        # Build a Task to run Abinit in a shell subprocess.
         task = AbinitTask.temp_shell_task(inp, workdir=workdir, manager=manager)
         task.start_and_wait(autoparal=False)
 
@@ -3068,7 +3059,6 @@ def abiget_irred_dteperts(self, ngkpt=None, shiftk=None, kptopt=None, ixc=None,
 
             [{'idir': 1, 'ipert': 4, 'qpt': [0.0, 0.0, 0.0]},
              {'idir': 2, 'ipert': 4, 'qpt': [0.0, 0.0, 0.0]}]
-
         """
         dteperts_vars = dict(d3e_pert1_phon=1 if phonon_pert else 0,  # phonon-type perturbation
                              d3e_pert2_phon=0,
@@ -3142,7 +3132,7 @@ def pop_par_vars(self, all=False) -> dict:
         """
         Remove all the variables associated to parallelism from the input file.
         Useful in case of a restart when we need to remove the parallel variables
-        before rerunning autoparal
+        before rerunning with autoparal.
         """
         parvars = ['npkpt', 'npfft', 'npband', 'npspinor', 'npimage']
         if all:
@@ -3465,7 +3455,11 @@ def has_same_structures(self) -> bool:
     def __str__(self) -> str:
         return self.to_string()
 
-    def to_string(self, mode="text", verbose=0, with_pseudos=True, files_file=False) -> str:
+    def to_string(self,
+                  mode: str = "text",
+                  verbose: int =0,
+                  with_pseudos: bool = True,
+                  files_file: bool = False) -> str:
         """
         String representation i.e. the ABINIT input file.
 
@@ -3803,7 +3797,7 @@ def modes_at_qpoint(cls, structure, qpoint, asr=2, chneut=1, dipdip=1, dipquad=1
         qpoint = qpoint.frac_coords if hasattr(qpoint, "frac_coords") else np.array(qpoint)
 
         if len(qpoint) != 3:
-            raise ValueError("Wrong q-point %s" % qpoint)
+            raise ValueError(f"Wrong {qpoint=}")
 
         return cls.modes_at_qpoints(structure=structure, qpoints=[qpoint], asr=asr, chneut=chneut, dipdip=dipdip,
                                     dipquad=dipquad, quadquad=quadquad,
@@ -3942,7 +3936,7 @@ def phbands_and_dos(cls, structure, ngqpt, nqsmall, qppa=None, ndivsm=20, line_d
 
                 dossmear = Energy(float(value), eunit).to("Ha")
         else:
-            raise NotImplementedError("Wrong value for dos_method: %s" % str(dos_method))
+            raise NotImplementedError(f"Wrong value for {dos_method=}")
 
         new = cls(structure, comment="ANADDB input for phonon bands and DOS generated by AbiPy" if not comment else comment,
                   anaddb_args=anaddb_args, anaddb_kwargs=anaddb_kwargs, spell_check=spell_check)
@@ -4240,7 +4234,7 @@ def to_string(self, sortmode=None, mode="text", verbose=0, files_file=False) ->
             # alphabetical order.
             keys = sorted(self.keys())
         else:
-            raise ValueError("Unsupported value for sortmode %s" % str(sortmode))
+            raise ValueError(f"Unsupported value for {sortmode=}")
 
         if mode == "html":
             var_database = get_anaddb_variables()
@@ -4448,7 +4442,7 @@ def as_dict(self) -> dict:
             value = self.vars.get(name)
             if value is None: value = self.get_default(name)
             if value is None:
-                raise self.Error("Variable %s is missing" % name)
+                raise self.Error(f"Variable {name=} is missing")
 
             var = self._NAME2VAR[name]
             grp = var.group
@@ -4465,7 +4459,7 @@ def to_string(self, verbose: int = 0, files_file: bool = False) -> str:
             value = self.vars.get(name)
             if value is None: value = self.get_default(name)
             if value is None:
-                raise self.Error("Variable %s is missing" % name)
+                raise self.Error(f"Variable {name=} is missing")
 
             # One line per variable --> valperline set to None
             variable = InputVariable("", value, valperline=None)
diff --git a/abipy/abio/robots.py b/abipy/abio/robots.py
index e8cdee48e..3ff12aec5 100644
--- a/abipy/abio/robots.py
+++ b/abipy/abio/robots.py
@@ -1,7 +1,9 @@
 # coding: utf-8
 """
-This module defines the Robot BaseClass. Robots operates on multiple files and provide helper
-functions to plot the data e.g. convergence studies and to build pandas dataframes from the output files.
+This module defines the Robot BaseClass.
+Robots operates on multiple files and provide helper
+functions to plot the data e.g. convergence studies
+and to build pandas dataframes from the output files.
 """
 from __future__ import annotations
 
@@ -150,7 +152,7 @@ def _open_files_in_dir(cls, top: str, walk: bool) -> list:
         Open files in directory tree starting from `top`. Return list of Abinit files.
         """
         if not os.path.isdir(top):
-            raise ValueError("%s: no such directory" % str(top))
+            raise ValueError(f"{top=}: no such directory")
         from abipy.abilab import abiopen
         items = []
         if walk:
@@ -431,7 +433,7 @@ def add_file(self, label, abifile, filter_abifile=None) -> None:
             self._do_close[abifile.filepath] = True
 
         if label in self._abifiles:
-            raise ValueError("label %s is already present!" % label)
+            raise ValueError(f"{label=} is already present!")
 
         self._abifiles[label] = abifile
 
@@ -773,7 +775,7 @@ def is_sortable(self, aname: str, raise_exc: bool = False) -> bool:
     %s
 
 Note that this list is automatically generated.
-Not all entries are sortable (Please select number-like quantities)""" % (self.__class__.__name__, aname, str(attrs)))
+Not all entries are sortable (please select number-like quantities)""" % (self.__class__.__name__, aname, str(attrs)))
 
     def _sortby_labelfile_list(self, labelfile_list, func_or_string, reverse=False, unpack=False):
         """
diff --git a/abipy/dfpt/gruneisen.py b/abipy/dfpt/gruneisen.py
index 9db80be8d..bd21ceec8 100644
--- a/abipy/dfpt/gruneisen.py
+++ b/abipy/dfpt/gruneisen.py
@@ -10,7 +10,6 @@
 
 from functools import lru_cache
 from collections import OrderedDict
-#from typing import List
 from monty.string import marquee, list_strings
 from monty.termcolor import cprint
 from monty.collections import AttrDict
@@ -77,7 +76,7 @@ def __str__(self):
         return self.to_string()
 
     def to_string(self, verbose: int = 0) -> str:
-        """String representation."""
+        """String representation with verbosite level `verbose`."""
         lines = []; app = lines.append
 
         app(marquee("File Info", mark="="))
@@ -185,7 +184,6 @@ def to_dataframe(self) -> pd.DataFrame:
         freq            Phonon frequency in eV.
         ==============  ==========================
         """
-
         grun_vals = self.gvals_qibz
         nqibz, natom3 = grun_vals.shape
         phfreqs = self.reader.rootgrp.variables["gruns_wvols_qibz"][:, self.iv0, :] * abu.Ha_eV
@@ -323,7 +321,7 @@ def plot_phbands_with_gruns(self, fill_with="gruns", gamma_fact=1, alpha=0.6, wi
             max_gamma = np.abs(self.grun_vals_finite_differences(match_eigv=True)).max()
             values = self.grun_vals_finite_differences(match_eigv=True)
         else:
-            raise ValueError("Unsupported fill_with: `%s`" % fill_with)
+            raise ValueError(f"Unsupported {fill_with=}")
 
         # Plot gruneisen markers on top of band structure.
         xvals = np.arange(len(phbands.phfreqs))
@@ -399,7 +397,7 @@ def plot_gruns_scatter(self, values="gruns", ax=None, units="eV", cmap="rainbow"
         elif values == "gruns_fd":
             y = self.gvals_qibz_finite_differences(match_eigv=True)
         else:
-            raise ValueError("Unsupported values: `%s`" % values)
+            raise ValueError(f"Unsupported {values=}")
 
         w = self.wvols_qibz[:, self.iv0, :] * abu.phfactor_ev2units(units)
 
@@ -481,7 +479,7 @@ def plot_gruns_bs(self, values="gruns", ax=None, branch_range=None, qlabels=None
         elif values == "gruns_fd":
             y = self.split_gruns_finite_differences(match_eigv=True)
         else:
-            raise ValueError("Unsupported values: `%s`" % values)
+            raise ValueError(f"Unsupported {values=}")
 
         phbands = self.phbands_qpath_vol[self.iv0]
 
@@ -673,9 +671,11 @@ def acoustic_debye_temp(self):
         return self.phdos.get_acoustic_debye_temp(len(self.structure))
 
     @classmethod
-    def from_ddb_list(cls, ddb_list, nqsmall=10, qppa=None, ndivsm=20, line_density=None, asr=2, chneut=1, dipdip=1,
-                     dos_method="tetra", lo_to_splitting="automatic", ngqpt=None, qptbounds=None, anaddb_kwargs=None,
-                     verbose=0, mpi_procs=1, workdir=None, manager=None):
+    def from_ddb_list(cls, ddb_list, nqsmall=10, qppa=None, ndivsm=20, line_density=None,
+                      asr=2, chneut=1, dipdip=1,
+                      dos_method="tetra", lo_to_splitting="automatic",
+                      ngqpt=None, qptbounds=None, anaddb_kwargs=None,
+                      verbose=0, mpi_procs=1, workdir=None, manager=None):
         """
         Execute anaddb to compute generate the object from a list of ddbs.
 
@@ -748,9 +748,7 @@ def from_ddb_list(cls, ddb_list, nqsmall=10, qppa=None, ndivsm=20, line_density=
         if not report.run_completed:
             raise ddb0.AnaddbError(task=task, report=report)
 
-        gruns = cls.from_file(os.path.join(task.workdir, "run.abo_GRUNS.nc"))
-
-        return gruns
+        return cls.from_file(os.path.join(task.workdir, "run.abo_GRUNS.nc"))
 
     @lru_cache()
     def grun_vals_finite_differences(self, match_eigv=True):
diff --git a/abipy/dfpt/phtk.py b/abipy/dfpt/phtk.py
index 98df80412..7c1efc151 100644
--- a/abipy/dfpt/phtk.py
+++ b/abipy/dfpt/phtk.py
@@ -199,9 +199,11 @@ def has_direction(self, direction, cartesian=False) -> bool:
             return False
 
 
-def open_file_phononwebsite(filename, port=8000,
+def open_file_phononwebsite(filename,
+                            port=8000,
                             website="http://henriquemiranda.github.io/phononwebsite",
-                            host="localhost", browser=None): # pragma: no cover
+                            host="localhost",
+                            browser=None): # pragma: no cover
     """
     Take a file, detect the type and open it on the phonon website
     Based on a similar function in <https://github.com/henriquemiranda/phononwebsite/phononweb.py>
diff --git a/abipy/dfpt/raman.py b/abipy/dfpt/raman.py
index 2e86f2c0c..0cff9f5dc 100644
--- a/abipy/dfpt/raman.py
+++ b/abipy/dfpt/raman.py
@@ -155,7 +155,6 @@ def _get_prefactor(w, temp, laser_freq) -> np.array:
         Returns:
             An array with shape (n modes) with the coefficient for the Raman intensities.
         """
-
         c = np.zeros_like(w)
         ind = np.where(w > 1e-5)
 
@@ -165,7 +164,7 @@ def _get_prefactor(w, temp, laser_freq) -> np.array:
 
         return c
 
-    def _get_lorentz_freqs_and_factor(self, intensity, non_anal_dir, min_freq, max_freq, 
+    def _get_lorentz_freqs_and_factor(self, intensity, non_anal_dir, min_freq, max_freq,
                                       num, width, units) -> tuple:
         """
         Helper method to get the list of frequencies and the main spread factors to
@@ -187,7 +186,6 @@ def _get_lorentz_freqs_and_factor(self, intensity, non_anal_dir, min_freq, max_f
             Tuple with list of "num" frequencies in eV and factors for the Lorentz broadening
             with shape (n modes, num).
         """
-
         if non_anal_dir is None:
             w = self.phfreqs
         else:
@@ -287,7 +285,7 @@ def get_lorentz_intensity(self, temp, laser_freq, width, non_anal_dir=None, min_
 
             return li_func
 
-    def get_powder_intensity(self, temp, laser_freq, non_anal_dir=None, 
+    def get_powder_intensity(self, temp, laser_freq, non_anal_dir=None,
                              relative=False, units="eV") -> PowderIntensity:
         """
         Calculates the Raman intensities in arbitrary units for each mode integrated over all possible
@@ -309,7 +307,6 @@ def get_powder_intensity(self, temp, laser_freq, non_anal_dir=None,
             A PowderIntensity with the parallel, perpendicular and total components of the powder
             intensities. Each one is an array with length n modes.
         """
-
         if non_anal_dir is None:
             w = self.phfreqs
             sus = self.susceptibility
@@ -369,7 +366,6 @@ def get_powder_lorentz_intensity(self, temp, laser_freq, width, non_anal_dir=Non
             A PowderIntensity with the parallel, perpendicular and total components of the powder
             intensities. Each one is a Function1D with "num" points.
         """
-
         pi = self.get_powder_intensity(temp=temp, laser_freq=laser_freq, non_anal_dir=non_anal_dir, units=units)
 
         freqs, lorentz = self._get_lorentz_freqs_and_factor(intensity=pi.tot, non_anal_dir=non_anal_dir, min_freq=min_freq,
@@ -386,7 +382,7 @@ def get_powder_lorentz_intensity(self, temp, laser_freq, width, non_anal_dir=Non
 
     @add_fig_kwargs
     def plot_intensity(self, temp, laser_freq, width, value, non_anal_dir=None, min_freq=None, max_freq=None,
-                       num=1000, relative=False, units="eV", ax=None, 
+                       num=1000, relative=False, units="eV", ax=None,
                        plot_phfreqs=False, **kwargs) -> Figure:
         """
         Plot one representation of the broadened Raman intensities.
@@ -420,7 +416,6 @@ def plot_intensity(self, temp, laser_freq, width, value, non_anal_dir=None, min_
 
         Returns: |matplotlib-Figure|
         """
-
         ax, fig, plt = get_ax_fig_plt(ax=ax)
 
         if width:
diff --git a/abipy/dfpt/vzsisa.py b/abipy/dfpt/vzsisa.py
index b7dd7c3f5..ce6251cea 100644
--- a/abipy/dfpt/vzsisa.py
+++ b/abipy/dfpt/vzsisa.py
@@ -27,7 +27,7 @@ def anaget_phdoses_with_gauss(nqsmall_or_qppa,
                               smearing_ev: float | None,
                               ddb_paths: list[PathLike],
                               anaget_kwargs: dict | None,
-                              verbose: int) -> tuple(list[str]) -> tuple[list[str], list[str]]:
+                              verbose: int) -> tuple[list[str], list[str]]:
     """
     Invoke anaddb to compute PHDOSes from a list of DDB filepaths with the gaussian method.
 
@@ -71,17 +71,6 @@ def anaget_phdoses_with_gauss(nqsmall_or_qppa,
         r = robot.anaget_phonon_plotters(**my_kwargs)
         return r.phdos_paths, r.phbands_paths
 
-    #for ddb_path in ddb_paths:
-    #    with DdbFile(ddb_path) as ddb:
-    #        with ddb.anaget_phbst_and_phdos_files(**my_kwargs) as g:
-    #            if verbose:
-    #                print(f"anaddb input file: {str(g.input)=}")
-    #            phbst_file, phdos_file = g[0], g[1]
-    #            phdos_paths.append(phdos_file.filepath)
-    #            phbands_paths.append(phbst_file.filepath)
-
-    #return phdos_paths, phbands_paths
-
 
 class Vzsisa(HasPickleIO):
     """
@@ -656,7 +645,7 @@ def get_phbands_plotter(self) -> PhononBandsPlotter:
         return plotter
 
     def get_edos_plotter(self) -> ElectronDosPlotter:
-        """Build and return a ElectronDosPlotter with electron doses indexed by volume"""
+        """Build and return a ElectronDosPlotter with electron DOSEs indexed by volume"""
         if self.edoses[0] is None:
             raise ValueError("edoses_list is not available")
         plotter = ElectronDosPlotter()
@@ -1876,7 +1865,7 @@ def plot_angles_vs_t_4th(self, tstart=0, tstop=1000, num=101, angle=None, tref=N
         if angle is None or angle == 3:
             ax.plot(tmesh, gamma, color='m', lw=2, label=r"$gamma(V(T))$" + method)
 
-        if abs(abs(self.bo_volumes[self.iv0]- ph_volumes[iv0])-abs(ph_volumes[iv1]-self.bo_volumes[self.iv0])) < 1e-3:
+        if abs(abs(self.bo_volumes[self.iv0]- ph_volumes[iv0]) - abs(ph_volumes[iv1]-self.bo_volumes[self.iv0])) < 1e-3:
             if angle is None or angle == 1:
                 ax.plot(tmesh, alpha2, linestyle='dashed', color='r', lw=2, label=r"$alpha(V(T))$""E2vib1")
             if angle is None or angle == 2:
diff --git a/abipy/dynamics/analyzer.py b/abipy/dynamics/analyzer.py
index 4324c1479..4da596da4 100644
--- a/abipy/dynamics/analyzer.py
+++ b/abipy/dynamics/analyzer.py
@@ -621,7 +621,7 @@ def latex_formula_n_temp(self) -> str:
 
     @property
     def latex_avg_volume(self) -> str:
-        return "V$_{\mathrm{ave}}$ = " + f"{self.avg_volume:.2f}" + '$\mathrm{{\AA}^3}$'
+        return r"V$_{\mathrm{ave}}$ = " + f"{self.avg_volume:.2f}" + r'$\mathrm{{\AA}^3}$'
 
     @property
     def avg_volume(self) -> float:
@@ -819,7 +819,7 @@ def plot_sqdt_symbols(self, symbols, t0: float = 0.0, atom_inds=None, with_dw=0,
 
         ax.legend(fontsize=fontsize, loc="upper left")
         ax.set_xlabel('t (ps)', fontsize=fontsize)
-        ax.set_ylabel('mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
+        ax.set_ylabel(r'mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
         set_axlims(ax, xlims, "x")
         #set_ticks_fontsize(ax, fontsize)
         set_logscale(ax, xy_log)
@@ -857,14 +857,14 @@ def plot_sqdt_symbols_tmax(self, symbols, tmax: float, atom_inds=None, nprocs=No
             ts = self.times[t_start:] - self.times[t_start]
 
             ax.plot(ts, msd_t,
-                    label=symbol + " <msd($t$, $t_0$)>$\{$t_0$\}$, $t$ = [0, " + str(int(self.times[index_tmax])) + " ps]",
+                    label=symbol + r" <msd($t$, $t_0$)>$\{$t_0$\}$, $t$ = [0, " + str(int(self.times[index_tmax])) + " ps]",
                     color=self.color_symbol[symbol],
                     )
 
         set_axlims(ax, xlims, "x")
         ax.legend(fontsize=fontsize, loc="upper left")
         ax.set_xlabel('t (ps)', fontsize=fontsize)
-        ax.set_ylabel('average mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
+        ax.set_ylabel(r'average mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
         #set_ticks_fontsize(ax, fontsize)
         set_logscale(ax, xy_log)
         ax.add_artist(AnchoredText(f"{self.latex_formula_n_temp}\n{self.latex_avg_volume}",
@@ -1001,7 +1001,7 @@ def plot(self, ax=None, xy_log=None, fontsize=8, xlims=None, **kwargs) -> Figure
 
         ax, fig, plt = get_ax_fig_plt(ax=ax)
         ax.plot(ts, self.msd_t,
-                label=self.symbol + " <msd($t, t_0$)>$\{$t_0$\}$, t = [0, " + str(int(self.times[index_tmax])) + " ps]",
+                label=self.symbol + r" <msd($t, t_0$)>$\{$t_0$\}$, t = [0, " + str(int(self.times[index_tmax])) + " ps]",
                 color=self.mda.color_symbol[self.symbol],
                 )
 
@@ -1191,7 +1191,6 @@ def plot(self, sharex=True, sharey=True, fontsize=8, **kwargs) -> Figure:
         return fig
 
 
-
 @dataclasses.dataclass(kw_only=True)
 class SigmaBerend:
     """
@@ -1231,7 +1230,7 @@ def plot(self, fontsize=8, ax_list=None, **kwargs) -> Figure:
                               self.latex_formula + ', '+ 'T = %4.0f' % self.temperature + 'K')
             ax.legend(fontsize=fontsize, loc="lower right")
             ax.set_xlabel('N. of data in block', fontsize=fontsize)
-            ax.set_ylabel('$\sigma$ ($\AA^2$)', fontsize=fontsize)
+            ax.set_ylabel(r'$\sigma$ ($\AA^2$)', fontsize=fontsize)
             ax.grid(True)
 
         fig.suptitle("Variance of correlated data as function of block number")
@@ -1444,7 +1443,7 @@ def expose(self, exposer="mpl", **kwargs):
 
 class MultiMdAnalyzer(HasPickleIO):
     """
-    High-level interface to analyze multiple MD trajectories
+    High-level interface to analyze multiple MD trajectories.
     """
 
     @classmethod
@@ -1618,7 +1617,7 @@ def plot_sqdt_symbols(self, symbols, t0: float = 0.0,
             set_axlims(ax, xlims, "x")
             ax.legend(fontsize=fontsize, loc="upper left")
             ax.set_xlabel('t (ps)', fontsize=fontsize)
-            ax.set_ylabel('average mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
+            ax.set_ylabel(r'average mean square displacement ($\mathrm{{\AA}^2}$)', fontsize=fontsize)
             #set_ticks_fontsize(ax, fontsize)
             set_logscale(ax, xy_log)
 
@@ -1982,7 +1981,8 @@ def add_entry_from_file(self, filepath: PathLike, key: str, mpl_style=None) -> N
         self.append(ArrheniusEntry.from_file(filepath, key, mpl_style))
 
     @add_fig_kwargs
-    def plot(self, thinvt_arange=None, what="diffusion", ncar=None, colormap="jet", with_t=True, text=None,
+    def plot(self, thinvt_arange=None, what="diffusion", ncar=None,
+             colormap="jet", with_t=True, text=None,
              ax=None, fontsize=8, xlims=None, ylims=None, **kwargs) -> Figure:
         """
         Arrhenius plot.
diff --git a/abipy/dynamics/cpx.py b/abipy/dynamics/cpx.py
index 76e729a4c..f1625ff41 100644
--- a/abipy/dynamics/cpx.py
+++ b/abipy/dynamics/cpx.py
@@ -158,7 +158,7 @@ def times(self) -> np.ndarray:
         """Array with times in ps units."""
         return np.array(self.df["tps(ps)"].values, dtype=float)
 
-    def to_string(self, verbose=0) -> str:
+    def to_string(self, verbose: int = 0) -> str:
         """String representation with verbosity level verbose."""
         lines = []; app = lines.append
         app(self.df.describe(percentiles=None).to_string())
@@ -419,4 +419,3 @@ def downsample_xyz(input_xyz, take_every, output_xyz, skip_head=None, verbose=1)
 
     if verbose: print(f"Wrote {count=} configurations to {output_xyz=} with {take_every=} and {skip_head=}")
     return count
-
diff --git a/abipy/dynamics/hist.py b/abipy/dynamics/hist.py
index a565592d5..d90f00cab 100644
--- a/abipy/dynamics/hist.py
+++ b/abipy/dynamics/hist.py
@@ -1,5 +1,7 @@
 # coding: utf-8
-"""History file with structural relaxation results."""
+"""
+Interface with the HIST.nc file containing ABINIT structural relaxation results or MD
+"""
 from __future__ import annotations
 
 import os
@@ -206,7 +208,7 @@ def write_xdatcar(self, filepath="XDATCAR", groupby_type=True, overwrite=False,
         from pymatgen.io.vasp.outputs import Xdatcar
 
         if filepath is not None and os.path.exists(filepath) and not overwrite:
-            raise RuntimeError("Cannot overwrite pre-existing file `%s`" % filepath)
+            raise RuntimeError(f"Cannot overwrite pre-existing file: {filepath}")
 
         if filepath is None:
             import tempfile
@@ -219,7 +221,7 @@ def write_xdatcar(self, filepath="XDATCAR", groupby_type=True, overwrite=False,
         ntypat = self.r.read_dimvalue("ntypat")
         num_pseudos = self.r.read_dimvalue("npsp")
         if num_pseudos != ntypat:
-            raise NotImplementedError("Alchemical mixing is not supported, num_pseudos != ntypat")
+            raise NotImplementedError("Alchemical mixing is not supported, {num_pseudos=} != {ntypat=}")
         #print("znucl:", znucl, "\ntypat:", typat)
 
         symb2pos = OrderedDict()
@@ -280,7 +282,7 @@ def visualize(self, appname="ovito", to_unit_cell=False):  # pragma: no cover
         from abipy.iotools import Visualizer
         visu = Visualizer.from_name(appname)
         if visu.name != "ovito":
-            raise NotImplementedError("visualizer: %s" % visu.name)
+            raise NotImplementedError(f"{visu.name=} is not supported")
 
         filepath = self.write_xdatcar(filepath=None, groupby_type=True, to_unit_cell=to_unit_cell)
 
@@ -486,7 +488,7 @@ def plotly_traces(self, fig, what, rcd=None, fontsize=8, showlegend=False, **kwa
             fig.layout['yaxis%u' % rcd.iax].title.text = 'F stats (eV/A)'
 
         else:
-            raise ValueError("Invalid value for what: `%s`" % str(what))
+            raise ValueError(f"Invalid value for {what=}")
 
         fig.layout.legend.font.size = fontsize
 
@@ -928,7 +930,7 @@ def read_all_structures(self) -> list[Structure]:
         num_pseudos = self.read_dimvalue("npsp")
         ntypat = self.read_dimvalue("ntypat")
         if num_pseudos != ntypat:
-            raise NotImplementedError("Alchemical mixing is not supported, num_pseudos != ntypat")
+            raise NotImplementedError("Alchemical mixing is not supported, {num_pseudos=} != {ntypat=}")
 
         znucl, typat = self.read_value("znucl"), self.read_value("typat").astype(int)
         #print(znucl.dtype, typat)
diff --git a/abipy/electrons/arpes.py b/abipy/electrons/arpes.py
index 5e6a1abd7..e0f8fe102 100644
--- a/abipy/electrons/arpes.py
+++ b/abipy/electrons/arpes.py
@@ -7,22 +7,17 @@
 from scipy.interpolate import UnivariateSpline
 from monty.collections import dict2namedtuple
 from abipy.core.mixins import Has_Structure, Has_ElectronBands, NotebookWriter
+from abipy.tools.typing import Figure
 from abipy.electrons import ElectronBands
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_ax3d_fig_plt, get_axarray_fig_plt #set_axlims,
 
 
 class ArpesPlotter(Has_Structure, Has_ElectronBands, NotebookWriter):
     """
-    Usage example:
-
-    .. code-block:: python
-
-        with abilab.abiopen("foo_ABIWAN.nc") as abiwan:
-            print(abiwan)
-
     .. rubric:: Inheritance Diagram
     .. inheritance-diagram:: ArpesPlotter
     """
+
     @classmethod
     def model_from_ebands(cls, ebands, tmesh=(0, 300, 600), poorman_polaron=False):
         ebands = ElectronBands.as_ebands(ebands)
@@ -101,10 +96,10 @@ def ebands(self):
         """|ElectronBands| object."""
         return self._ebands
 
-    def __str__(self):
+    def __str__(self) -> str:
         return self.to_string()
 
-    def to_string(self, verbose=0):
+    def to_string(self, verbose: int = 0) -> str:
         """String representation with verbosity level `verbose`."""
         lines = []; app = lines.append
         app(self.structure.to_string(verbose=verbose, title="Structure"))
@@ -148,7 +143,7 @@ def with_points_along_path(self, frac_bounds=None, knames=None, dist_tol=1e-12):
 
     #    return self.__class__(new_ebands, new_aw, aw_meshes, self.tmesh)
 
-    def get_emesh_eminmax(self, estep):
+    def get_emesh_eminmax(self, estep: float) -> np.ndarray:
         """Compute linear mesh covering entire energy range."""
         emin = self.ebands.enemin()
         emin -= 0.1 * abs(emin)
@@ -189,7 +184,7 @@ def get_atw(self, wmesh, spin, ikpt, band_inds, temp_inds):
 
     @add_fig_kwargs
     def plot_ekmap_temps(self, temp_inds=None, spins=None, estep=0.02, with_colorbar=True,
-                        ylims=None, fontsize=8, **kwargs):
+                        ylims=None, fontsize=8, **kwargs) -> Figure:
         """
         Plot (k, e) color maps for different temperatures.
 
@@ -220,7 +215,7 @@ def plot_ekmap_temps(self, temp_inds=None, spins=None, estep=0.02, with_colorbar
         return fig
 
     @add_fig_kwargs
-    def plot_ekmap_itemp(self, itemp=0, spins=None, estep=0.02, ax=None, ylims=None, with_colorbar=True, **kwargs):
+    def plot_ekmap_itemp(self, itemp=0, spins=None, estep=0.02, ax=None, ylims=None, with_colorbar=True, **kwargs) -> Figure:
         """
         Plot (k, e) color map for given temperature.
 
@@ -279,7 +274,7 @@ def plot_ekmap_itemp(self, itemp=0, spins=None, estep=0.02, ax=None, ylims=None,
 
     @add_fig_kwargs
     def plot_ak_vs_temp(self, temp_inds=None, spins=None, band_inds=None, kpt_inds=None,
-                        apad=1.0, estep=0.02, colormap="jet", fontsize=8, **kwargs):
+                        apad=1.0, estep=0.02, colormap="jet", fontsize=8, **kwargs) -> Figure:
         """
 
         Args:
@@ -388,7 +383,7 @@ def plot_ak_vs_temp(self, temp_inds=None, spins=None, band_inds=None, kpt_inds=N
     #    return fig
 
     @add_fig_kwargs
-    def plot_3dlines(self, itemp=0, estep=0.02, spins=None, band_inds=None, ax=None, **kwargs):
+    def plot_3dlines(self, itemp=0, estep=0.02, spins=None, band_inds=None, ax=None, **kwargs) -> Figure:
         ax, fig, plt = get_ax3d_fig_plt(ax=ax)
 
         xs, emin, emax = self.get_emesh_eminmax(estep)
diff --git a/abipy/electrons/bse.py b/abipy/electrons/bse.py
index 38c191e85..df6e63cb3 100644
--- a/abipy/electrons/bse.py
+++ b/abipy/electrons/bse.py
@@ -33,7 +33,7 @@
 
 
 # Deprecated: should be rewritten from scratch.
-class _DielectricTensor(object):
+class _DielectricTensor:
     """
     This object stores the frequency-dependent macroscopic dielectric tensor
     obtained from the dielectric functions for different q-directions.
@@ -152,7 +152,7 @@ def plot_ax(self, ax, what, red_coords, *args, **kwargs) -> list:
         if duck.is_intlike(what):
             f = self.to_func1d(red_coords)[int(what)]
         else:
-            raise ValueError("Don't know how to handle %s" % str(what))
+            raise ValueError(f"Don't know how to handle {what=}")
 
         return f.plot_ax(ax, *args, **kwargs)
 
@@ -301,7 +301,7 @@ def plot_ax(self, ax, qpoint=None, **kwargs) -> list:
             f = self.emacro_avg
 
         else:
-            raise ValueError("Don't know how to handle %s" % str(qpoint))
+            raise ValueError(f"Don't know how to handle {qpoint=}")
 
         return f.plot_ax(ax, **kwargs)
 
@@ -326,7 +326,7 @@ def from_file(cls, filepath: str) -> MdfFile:
 
     def __init__(self, filepath: str):
         super().__init__(filepath)
-        self.reader = MdfReader(filepath)
+        self.r = MdfReader(filepath)
 
         # TODO Add electron Bands.
         #self._ebands = r.read_ebands()
@@ -351,32 +351,32 @@ def to_string(self, verbose=0) -> str:
 
     def close(self) -> None:
         """Close the file."""
-        self.reader.close()
+        self.r.close()
 
     @lazy_property
     def structure(self) -> Structure:
         """|Structure| object."""
-        return self.reader.read_structure()
+        return self.r.read_structure()
 
     @lazy_property
     def exc_mdf(self):
         "Excitonic macroscopic dieletric function."""
-        return self.reader.read_exc_mdf()
+        return self.r.read_exc_mdf()
 
     @lazy_property
     def rpanlf_mdf(self):
         """RPA dielectric function without local-field effects."""
-        return self.reader.read_rpanlf_mdf()
+        return self.r.read_rpanlf_mdf()
 
     @lazy_property
     def gwnlf_mdf(self):
         """RPA-GW dielectric function without local-field effects."""
-        return self.reader.read_gwnlf_mdf()
+        return self.r.read_gwnlf_mdf()
 
     @property
     def qpoints(self):
         """List of q-points."""
-        return self.reader.qpoints
+        return self.r.qpoints
 
     @property
     def qfrac_coords(self) -> np.ndarray:
@@ -389,7 +389,7 @@ def params(self) -> dict:
         Dictionary with the parameters that are usually tested for convergence.
         Used to build Pandas dataframes in Robots.
         """
-        return self.reader.read_params()
+        return self.r.read_params()
 
     def get_mdf(self, mdf_type="exc"):
         """"
@@ -552,7 +552,7 @@ def read_gwnlf_mdf(self):
 
 class MdfPlotter:
     """
-    Class for plotting Macroscopic dielectric functions.
+    Class for plotting macroscopic dielectric functions.
 
     Usage example:
 
@@ -580,7 +580,7 @@ def add_mdf(self, label, mdf):
         self._mdfs[label] = mdf
 
     @add_fig_kwargs
-    def plot(self, ax=None, cplx_mode="Im", qpoint=None, xlims=None, ylims=None, 
+    def plot(self, ax=None, cplx_mode="Im", qpoint=None, xlims=None, ylims=None,
              fontsize=8, **kwargs) -> Figure:
         """
         Get a matplotlib plot showing the MDFs.
@@ -621,21 +621,6 @@ def plot(self, ax=None, cplx_mode="Im", qpoint=None, xlims=None, ylims=None,
 
         return fig
 
-    #def ipw_plot(self) # pragma: no cover
-    #    """
-    #    Return an ipython widget with controllers to select the plot.
-    #    """
-    #    def plot_callback(plot_type, qpoint):
-    #        if qpoint == "None": qpoint = None
-    #        return self.plot(cplx_type=cplx_type, qpoint=qpoint)
-
-    #    import ipywidgets as ipw
-    #    return ipw.interact_manual(
-    #            plot_callback,
-    #            cplx_type=["re", "im", "abs"],
-    #            qpoint=["None"] + list(range(self.,
-    #        )
-
 
 class MultipleMdfPlotter:
     """
@@ -711,7 +696,7 @@ def add_mdf_file(self, label, obj):
                 self._mdfs[label][mdf_type] = obj.get_mdf(mdf_type=mdf_type)
 
     @add_fig_kwargs
-    def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None, 
+    def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None,
              fontsize=8, **kwargs) -> Figure:
         """
         Plot all macroscopic dielectric functions (MDF) stored in the plotter
@@ -736,7 +721,7 @@ def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None,
             qpoints = self._get_qpoints()
             ncols, nrows = 2, len(qpoints)
         else:
-            raise ValueError("Invalid value of qview: %s" % str(qview))
+            raise ValueError(f"Invalid value of {qview=}")
 
         ax_mat, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
                                                sharex=True, sharey=True, squeeze=False)
@@ -747,6 +732,7 @@ def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None,
                                    fontsize=fontsize, with_legend=True, show=False)
             self.plot_mdftype_cplx(mdf_type, "Im", ax=ax_mat[0, 1], xlims=xlims, ylims=ylims,
                                    fontsize=fontsize, with_legend=False, show=False)
+
         elif qview == "all":
             # Plot MDF(q)
             nqpt = len(qpoints)
@@ -758,14 +744,14 @@ def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None,
                     fontsize=fontsize, with_legend=False, with_xlabel=islast, with_ylabel=islast, show=False)
 
         else:
-            raise ValueError("Invalid value of qview: `%s`" % str(qview))
+            raise ValueError(f"Invalid value of {qview=}")
 
         #ax_mat[0, 0].legend(loc="best", fontsize=fontsize, shadow=True)
 
         return fig
 
     #@add_fig_kwargs
-    #def plot_mdftypes(self, qview="avg", xlims=None, ylims=None, **kwargs):
+    #def plot_mdf_types(self, qview="avg", xlims=None, ylims=None, **kwargs):
     #    """
     #    Args:
     #        qview:
@@ -808,7 +794,7 @@ def plot(self, mdf_type="exc", qview="avg", xlims=None, ylims=None,
 
     @add_fig_kwargs
     def plot_mdftype_cplx(self, mdf_type, cplx_mode, qpoint=None, ax=None, xlims=None, ylims=None,
-                          with_legend=True, with_xlabel=True, with_ylabel=True, 
+                          with_legend=True, with_xlabel=True, with_ylabel=True,
                           fontsize=8, **kwargs) -> Figure:
         """
         Helper function to plot data corresponds to ``mdf_type``, ``cplx_mode``, ``qpoint``.
@@ -1028,7 +1014,7 @@ def _from_cart_to_red(cartesian_tensor,lattice):
 
 
 # TODO Remove
-class _Tensor(object):
+class _Tensor:
     """Representation of a 3x3 tensor"""
 
     def __init__(self, red_tensor, lattice, space="r"):
@@ -1049,7 +1035,7 @@ def __init__(self, red_tensor, lattice, space="r"):
         elif space == "r":
             self._is_real_space = True
         else:
-            raise ValueError("space should be either 'g' or 'r'")
+            raise ValueError(f"space should be either 'g' or 'r' but got {space=}")
 
     def __eq__(self, other):
         if other is None: return False
diff --git a/abipy/electrons/fatbands.py b/abipy/electrons/fatbands.py
index d21e99203..fd6e61436 100644
--- a/abipy/electrons/fatbands.py
+++ b/abipy/electrons/fatbands.py
@@ -21,6 +21,9 @@
 
 
 def gaussians_dos(dos, mesh, width, values, energies, weights):
+    """
+    Accumulate dos with the Gaussian method.
+    """
     assert len(dos) == len(mesh) and len(values) == len(energies) == len(weights)
     for vw, e, w in zip(values * weights, energies, weights):
         dos += vw * gaussian(mesh, width, center=e)
@@ -97,7 +100,6 @@ def __init__(self, filepath: str):
         self.r = self.reader = r = ElectronsReader(filepath)
 
         # Initialize the electron bands from file
-        self._ebands = r.read_ebands()
         self.natom = len(self.structure)
 
         # Read metadata so that we know how to handle the content of the file.
@@ -195,7 +197,7 @@ def _read_wal_sbk(self, key="dos_fractions"):
         # To keep it simple, the code always operate on an array dimensioned with the total number of atoms
         # Entries that are not computed are set to zero and a warning is issued.
         if self.prtdos != 3:
-            raise RuntimeError("The file does not contain L-DOS since prtdos=%i" % self.prtdos)
+            raise RuntimeError(f"The file does not contain L-DOS since {self.prtdos=}")
 
         wshape = (self.natom, self.mbesslang, self.nsppol, self.mband, self.nkpt)
 
@@ -240,9 +242,9 @@ def _read_walm_sbk(self, key="dos_fraction_m"):
         # To keep it simple, the code always operate on an array dimensioned with the total number of atoms
         # Entries that are not computed are set to zero and a warning is issued.
         if self.prtdos != 3:
-            raise RuntimeError("The file does not contain L-DOS since prtdos=%i" % self.prtdos)
+            raise RuntimeError(f"The file does not contain L-DOS since {self.prtdos=}")
         if self.prtdosm == 0:
-            raise RuntimeError("The file does not contain LM-DOS since prtdosm=%i" % self.prtdosm)
+            raise RuntimeError("The file does not contain LM-DOS since {self.prtdosm=}")
 
         wshape = (self.natom, self.mbesslang**2, self.nsppol, self.mband, self.nkpt)
 
@@ -275,10 +277,10 @@ def _read_walm_sbk(self, key="dos_fraction_m"):
 
         return walm_sbk
 
-    @property
+    @lazy_property
     def ebands(self) -> ElectronBands:
         """|ElectronBands| object."""
-        return self._ebands
+        return self.r.read_ebands()
 
     @property
     def structure(self) -> Structure:
@@ -295,7 +297,7 @@ def close(self) -> None:
         """Called at the end of the ``with`` context manager."""
         return self.r.close()
 
-    def __str__(self):
+    def __str__(self) -> str:
         """String representation"""
         return self.to_string()
 
@@ -412,7 +414,7 @@ def get_spilling(self, spin=None, band=None):
 
         return 1.0 - sp
 
-    def eb_plotax_kwargs(self, spin):
+    def eb_plotax_kwargs(self, spin: int) -> dict:
         """
         Dictionary with the options passed to ``ebands.plot_ax``
         when plotting a band line with spin index ``spin``.
@@ -426,7 +428,7 @@ def eb_plotax_kwargs(self, spin):
             #klabel_size=self.klabel_size,
         )
 
-    def eb_plotly_kwargs(self, spin):
+    def eb_plotly_kwargs(self, spin: int) -> tuple[dict]:
         """
         Dictionary with the options passed to ``ebands.plot_ax``
         when plotting a band line with spin index ``spin``.
@@ -476,7 +478,7 @@ def plot_fatbands_siteview(self, e0="fermie", view="inequivalent", fact=1.0, fon
             num_plots = len(ea.irred_pos)
             ax2iatom = ea.irred_pos
         else:
-            raise ValueError("Wrong value for view: %s" % str(view))
+            raise ValueError(f"Wrong value for {view=}")
 
         # Build plot grid.
         ncols, nrows = 1, 1
@@ -1063,8 +1065,8 @@ def plot_spilling(self, e0="fermie", fact=1.0, ax_list=None, ylims=None, blist=N
     #    return fig
 
     # TODO: THIS CODE IS STILL UNDER DEVELOPMENT
-    #def nelect_in_spheres(self, start_energy=None, stop_energy=None,
-    #                     method="gaussian", step=0.1, width=0.2):
+    #def get_nelect_in_spheres(self, start_energy=None, stop_energy=None,
+    #                          method="gaussian", step=0.1, width=0.2):
     #    """
     #    Print the number of electrons inside each atom-centered sphere.
     #    Note that this is a very crude estimate of the charge density distribution.
@@ -1105,7 +1107,7 @@ def plot_spilling(self, e0="fermie", fact=1.0, ax_list=None, ylims=None, blist=N
     #            nel_spin[spin] = edos[spin].integral(start=start_spin[spin], stop=stop_spin[spin])
     #        print("iatom", iatm, "site", site, nel_spin)
 
-    def get_dos_integrator(self, method, step, width):
+    def get_dos_integrator(self, method: str, step: float, width: float):
         """
         FatBandsFile can use differerent integrators that are cached in self._cached_dos_integrators
         """
@@ -1164,7 +1166,7 @@ def plot_pjdos_lview(self, e0="fermie", lmax=None, method="gaussian", step=0.1,
             msg += ("Error while trying to compute the DOS.\n"
                     "Verify that the k-points form a homogenous sampling of the BZ.\n"
                     "Returning None\n")
-            cprint(msg, "red")
+            cprint(msg, color="red")
             return None
 
         # Get energy mesh from total DOS and define the zero of energy
@@ -1322,7 +1324,7 @@ def plotly_pjdos_lview(self, e0="fermie", lmax=None, method="gaussian", step=0.1
             msg += ("Error while trying to compute the DOS.\n"
                     "Verify that the k-points form a homogenous sampling of the BZ.\n"
                     "Returning None\n")
-            cprint(msg, "red")
+            cprint(msg, color="red")
             return None
 
         # Get energy mesh from total DOS and define the zero of energy
@@ -1493,7 +1495,7 @@ def plot_pjdos_typeview(self, e0="fermie", lmax=None, l_list=None, method="gauss
             msg += ("Error while trying to compute the DOS.\n"
                     "Verify that the k-points form a homogenous sampling of the BZ.\n"
                     "Returning None\n")
-            cprint(msg, "red")
+            cprint(msg, color="red")
             return None
 
         # Get energy mesh from total DOS and define the zero of energy
@@ -1652,7 +1654,7 @@ def plotly_pjdos_typeview(self, e0="fermie", lmax=None, method="gaussian", step=
             msg += ("Error while trying to compute the DOS.\n"
                     "Verify that the k-points form a homogenous sampling of the BZ.\n"
                     "Returning None\n")
-            cprint(msg, "red")
+            cprint(msg, color="red")
             return None
 
         # Get energy mesh from total DOS and define the zero of energy
@@ -1815,8 +1817,8 @@ def plot_fatbands_with_pjdos(self, e0="fermie", fact=1.0, lmax=None, blist=None,
             pjdosfile = self
 
         if not pjdosfile.ebands.kpoints.is_ibz:
-            cprint("DOS requires k-points in the IBZ but got pjdosfile: %s" % repr(pjdosfile), "yellow")
-            cprint("Returning None", "yellow")
+            cprint("DOS requires k-points in the IBZ but got pjdosfile: %s" % repr(pjdosfile), color="yellow")
+            cprint("Returning None", color="yellow")
             return None
 
         if edos_kwargs is None: edos_kwargs = {}
@@ -1862,7 +1864,7 @@ def plot_fatbands_with_pjdos(self, e0="fermie", fact=1.0, lmax=None, blist=None,
                                           with_info=False, with_spin_sign=False, show=False, ylims=ylims,
                                           **edos_kwargs)
         else:
-            raise ValueError("Don't know how to handle view=%s" % str(view))
+            raise ValueError(f"Don't know how to handle {view=}")
 
         # Remove labels from DOS plots.
         for ax in pjdos_axmat.ravel():
@@ -1912,8 +1914,8 @@ def plotly_fatbands_with_pjdos(self, e0="fermie", fact=1.0, lmax=None, blist=Non
             pjdosfile = self
 
         if not pjdosfile.ebands.kpoints.is_ibz:
-            cprint("DOS requires k-points in the IBZ but got pjdosfile: %s" % repr(pjdosfile), "yellow")
-            cprint("Returning None", "yellow")
+            cprint("DOS requires k-points in the IBZ but got pjdosfile: %s" % repr(pjdosfile), color="yellow")
+            cprint("Returning None", color="yellow")
             return None
 
         if edos_kwargs is None: edos_kwargs = {}
@@ -1941,7 +1943,7 @@ def plotly_fatbands_with_pjdos(self, e0="fermie", fact=1.0, lmax=None, blist=Non
                                             combined_spins=False, fontsize=fontsize, with_info=False,
                                             with_spin_sign=False, ylims=ylims, band_and_dos=1, show=False, **edos_kwargs)
         else:
-            raise ValueError("Don't know how to handle view=%s" % str(view))
+            raise ValueError(f"Don't know how to handle {view=}")
 
         if closeit: pjdosfile.close()
         return fig
@@ -2009,7 +2011,7 @@ def plot_pawdos_terms(self, lmax=None, method="gaussian", step=0.1, width=0.2, x
                                 pawt1dos_al[iatom, l, spin] += weight * gs * pawt1_wal_sbk[iatom, l, spin, band, k]
 
         else:
-            raise ValueError("Method %s is not supported" % method)
+            raise ValueError(f"{method=} is not supported")
 
         # TOT = PW + AE - PS
         pwdos_al = totdos_al - paw1dos_al + pawt1dos_al
@@ -2275,7 +2277,7 @@ def symbols_lso(self):
                 symbols_lso[symbol] = lso
 
         else:
-            raise ValueError("Method %s is not supported" % self.method)
+            raise ValueError(f"{self.method=} is not supported")
 
         return symbols_lso
 
diff --git a/abipy/electrons/fold2bloch.py b/abipy/electrons/fold2bloch.py
index d981e886f..8998c677d 100644
--- a/abipy/electrons/fold2bloch.py
+++ b/abipy/electrons/fold2bloch.py
@@ -35,13 +35,14 @@ class Fold2BlochNcfile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBand
     .. rubric:: Inheritance Diagram
     .. inheritance-diagram:: Fold2BlochNcfile
     """
+
     @classmethod
-    def from_wfkpath(cls, wfkpath, folds, workdir=None, manager=None, mpi_procs=1, verbose=0) -> Fold2BlochNcfile:
+    def from_wfkpath(cls, wfk_path, folds, workdir=None, manager=None, mpi_procs=1, verbose=0) -> Fold2BlochNcfile:
         """
         Run fold2bloch in workdir.
 
         Args:
-            wfkpath:
+            wfk_path:
             folds:
             workdir: Working directory of the fake task used to compute the ibz. Use None for temporary dir.
             manager: :class:`TaskManager` of the task. If None, the manager is initialized from the config file.
@@ -55,9 +56,9 @@ def from_wfkpath(cls, wfkpath, folds, workdir=None, manager=None, mpi_procs=1, v
 
         # Create temporary directory and link to the WFK file
         workdir = get_workdir(workdir)
-        wfkpath = os.path.abspath(wfkpath)
-        link = os.path.join(workdir, os.path.basename(wfkpath))
-        os.symlink(wfkpath, link)
+        wfk_path = os.path.abspath(wfk_path)
+        link = os.path.join(workdir, os.path.basename(wfk_path))
+        os.symlink(wfk_path, link)
 
         # Run fold2bloch
         ncpath = fold2bloch.unfold(link, folds, workdir=workdir)
@@ -66,7 +67,7 @@ def from_wfkpath(cls, wfkpath, folds, workdir=None, manager=None, mpi_procs=1, v
 
     def __init__(self, filepath: str):
         super().__init__(filepath)
-        self.reader = ElectronsReader(filepath)
+        self.r = ElectronsReader(filepath)
 
         # Initialize the electron bands from file.
         # Spectral weights are dimensioned with `nss`
@@ -74,9 +75,8 @@ def __init__(self, filepath: str):
         # nctkarr_t("reduced_coordinates_of_unfolded_kpoints", "dp", "number_of_reduced_dimensions, nk_unfolded")
         # nctkarr_t("unfolded_eigenvalues", "dp", "max_number_of_states, nk_unfolded, number_of_spins")
         # nctkarr_t("spectral_weights", "dp", "max_number_of_states, nk_unfolded, nsppol_times_nspinor")
-        self._ebands = self.reader.read_ebands()
         self.nss = max(self.nsppol, self.nspinor)
-        self.fold_matrix = self.reader.read_value("fold_matrix")
+        self.fold_matrix = self.r.read_value("fold_matrix")
 
         # Compute direct lattice of the primitive cell from fold_matrix.
         if is_diagonal(self.fold_matrix):
@@ -86,15 +86,15 @@ def __init__(self, filepath: str):
             raise NotImplementedError("non diagonal fold_matrix: %s" % str(self.fold_matrix))
 
         # Read fold2bloch output data.
-        self.uf_kfrac_coords = self.reader.read_value("reduced_coordinates_of_unfolded_kpoints")
+        self.uf_kfrac_coords = self.r.read_value("reduced_coordinates_of_unfolded_kpoints")
         self.uf_kpoints = KpointList(self.pc_lattice.reciprocal_lattice, self.uf_kfrac_coords)
         self.uf_nkpt = len(self.uf_kpoints)
 
     def __str__(self) -> str:
         return self.to_string()
 
-    def to_string(self, verbose=0) -> str:
-        """String representation."""
+    def to_string(self, verbose: int = 0) -> str:
+        """String representation with verbosity level `verbose`."""
         lines = []; app = lines.append
 
         app(marquee("File Info", mark="="))
@@ -121,10 +121,10 @@ def to_string(self, verbose=0) -> str:
 
         return "\n".join(lines)
 
-    @property
+    @lazy_property
     def ebands(self) -> ElectronBands:
         """|ElectronBands| object with folded band energies."""
-        return self._ebands
+        return self.r.read_ebands()
 
     @property
     def structure(self) -> Structure:
@@ -133,7 +133,7 @@ def structure(self) -> Structure:
 
     def close(self) -> None:
         """Close the file."""
-        self.reader.close()
+        self.r.close()
 
     @lazy_property
     def params(self) -> dict:
@@ -145,13 +145,13 @@ def params(self) -> dict:
     def uf_eigens(self) -> np.ndarray:
         """[nsppol, nk_unfolded, nband] |numpy-array| with unfolded eigenvalues in eV."""
         # nctkarr_t("unfolded_eigenvalues", "dp", "max_number_of_states, nk_unfolded, number_of_spins")
-        return self.reader.read_value("unfolded_eigenvalues") * units.Ha_to_eV
+        return self.r.read_value("unfolded_eigenvalues") * units.Ha_to_eV
 
     @lazy_property
     def uf_weights(self) -> np.ndarray:
         """[nss, nk_unfolded, nband] array with spectral weights. nss = max(nspinor, nsppol)."""
         # nctkarr_t("spectral_weights", "dp", "max_number_of_states, nk_unfolded, nsppol_times_nspinor")
-        return self.reader.read_value("spectral_weights")
+        return self.r.read_value("spectral_weights")
 
     def get_spectral_functions(self, step=0.01, width=0.02):
         """
@@ -215,6 +215,7 @@ def plot_unfolded(self, kbounds, klabels, ylims=None, dist_tol=1e-12, verbose=0,
         if len(p.ikfound) == 0:
             cprint("Warning: find_points_along_path returned zero points along the path. Try to increase dist_tol.", "yellow")
             return None
+
         if verbose:
             uf_frac_coords = np.reshape([k.frac_coords for k in self.uf_kpoints], (-1, 3))
             fcoords = uf_frac_coords[p.ikfound]
diff --git a/abipy/electrons/gwr.py b/abipy/electrons/gwr.py
index 2e54c5dca..13b4cceab 100644
--- a/abipy/electrons/gwr.py
+++ b/abipy/electrons/gwr.py
@@ -476,7 +476,7 @@ def ks_dirgaps(self) -> np.ndarray:
     @lazy_property
     def qpz0_dirgaps(self) -> np.ndarray:
        """
-       QP direct gaps in eV computed with the Z factor at the KS energy
+       QP direct gaps in eV computed with the renormalization Z factor at the KS energy
        Shape: [nsppol, nkcalc]
        """
        return self.r.read_value("qpz_gaps") * abu.Ha_eV
@@ -523,14 +523,14 @@ def find_qpkinds(self, qp_kpoints) -> _MyQpkindsList:
                 qp_kpoints = [self.sigma_kpoints[ikc] for ikc in ik_list]
                 items = qp_kpoints, ik_list
         else:
-            raise TypeError("Don't know how to interpret `%s`" % (type(qp_kpoints)))
+            raise TypeError(f"Don't know how to interpret {type(qp_kpoints)}")
 
         # Check indices
         errors = []
         eapp = errors.append
         for ikc in items[1]:
             if ikc >= self.nkcalc:
-                eapp("K-point index %d >= nkcalc %d, check input qp_kpoints" % (ikc, self.nkcalc))
+                eapp(f"K-point index {ikc} >= {self.nkcalc=}. Please check input qp_kpoints")
         if errors:
             raise ValueError("\n".join(errors))
 
@@ -577,7 +577,7 @@ def to_string(self, verbose: int = 0) -> str:
 
         # GWR section.
         app(marquee("GWR parameters", mark="="))
-        app("gwr_task: %s" % self.r.gwr_task)
+        app(f"gwr_task: {self.r.gwr_task}")
 
         if self.r.gwr_task == GWR_TASK.RPA_ENERGY:
             pass
diff --git a/abipy/electrons/lruj.py b/abipy/electrons/lruj.py
index 4a0b73b38..42bb53de6 100644
--- a/abipy/electrons/lruj.py
+++ b/abipy/electrons/lruj.py
@@ -44,8 +44,6 @@
 #        self.r = r = EtsfReader(filepath)
 
 
-#===============================================================================================================
-#===============================================================================================================
 @dataclasses.dataclass(kw_only=True)
 class LrujResults:
     """
@@ -186,8 +184,6 @@ def find(header, dtype=None):
         _data = locals()
         return cls(**{k: _data[k] for k in [field.name for field in dataclasses.fields(cls)]})
 
-#===============================================================================================================
-#===============================================================================================================
     @add_fig_kwargs
     def plot(self, ax=None, degrees="all", inset=True, insetdegree=1, insetlocale="lower left",
              ptcolor0='k', ptcolor='k', gradcolor1='#3575D5',gradcolor2='#FDAE7B',
@@ -295,8 +291,6 @@ def dfvalue(keywo):
 
 
 
-#===============================================================================================================
-#===============================================================================================================
 class LrujAnalyzer:
     """
     Analyzes multiple sets of LRUJ files.
diff --git a/abipy/electrons/psps.py b/abipy/electrons/psps.py
index e12a1eb41..e7d5f1384 100644
--- a/abipy/electrons/psps.py
+++ b/abipy/electrons/psps.py
@@ -8,7 +8,7 @@
 import numpy as np
 import pandas as pd
 
-#from typing import List, Any
+#from typing import Any
 from collections import OrderedDict
 from monty.bisect import find_gt
 from monty.string import marquee # list_strings,
@@ -82,7 +82,7 @@ class PspsFile(AbinitNcFile, NotebookWriter):
 
     .. code-block:: python
 
-        with PspsFile("foo_PSPS.nc") as psps:
+        with PspsFile("out_PSPS.nc") as psps:
             psps.plot_tcore_rspace()
     """
     linestyles_der = ["-", "--", '-.', ':', ":", ":"]
@@ -101,7 +101,7 @@ def from_abinit_run(cls, pseudo, ecut, vloc_rcut=None, workdir=None) -> PspsFile
         abinit with cutoff energy `ecut`.
 
         Args:
-            vloc_rcut: Radial cutoff in Bohr (Abinit input variable)
+            vloc_rcut: Radial cutoff in Bohr (Abinit input variable).
         """
         from abipy.flowtk.pseudos import Pseudo
         pseudo = Pseudo.as_pseudo(pseudo)
@@ -126,14 +126,13 @@ def from_abinit_run(cls, pseudo, ecut, vloc_rcut=None, workdir=None) -> PspsFile
 
     def __init__(self, filepath: str):
         super().__init__(filepath)
-        self.reader = PspsReader(filepath)
-
+        self.r = PspsReader(filepath)
         # TODO
-        #self.ecut = self.reader.read_value("ecut")
+        #self.ecut = self.r.read_value("ecut")
 
     def close(self) -> None:
         """Close the file."""
-        self.reader.close()
+        self.r.close()
 
     @lazy_property
     def params(self) -> dict:
@@ -198,7 +197,7 @@ def plot_tcore_rspace(self, ax=None, ders=(0, 1, 2, 3), scale=1.0, rmax=3.0, **k
         """
         if not isinstance(ders, (list, tuple)): ders = [ders]
 
-        rmeshes, coresd = self.reader.read_coresd(rmax=rmax)
+        rmeshes, coresd = self.r.read_coresd(rmax=rmax)
         ax, fig, plt = get_ax_fig_plt(ax=ax)
 
         for rmesh, mcores in zip(rmeshes, coresd):
@@ -221,8 +220,7 @@ def plot_tcore_rspace(self, ax=None, ders=(0, 1, 2, 3), scale=1.0, rmax=3.0, **k
         return fig
 
     @add_fig_kwargs
-    def plot_tcore_qspace(self, ax=None, ders=(0,), with_fact=True,
-                          with_qn=0, scale=1.0, **kwargs) -> Figure:
+    def plot_tcore_qspace(self, ax=None, ders=(0,), with_fact=True, with_qn=0, scale=1.0, **kwargs) -> Figure:
         """
         Plot the model core charge in q space.
 
@@ -241,7 +239,7 @@ def plot_tcore_qspace(self, ax=None, ders=(0,), with_fact=True,
         color = kwargs.pop("color", "black")
         linewidth = kwargs.pop("linewidth", _LW)
 
-        qmesh, tcore_spl = self.reader.read_tcorespl()
+        qmesh, tcore_spl = self.r.read_tcorespl()
         ecuts = 2 * (np.pi * qmesh)**2
 
         lines = []
@@ -273,8 +271,7 @@ def plot_tcore_qspace(self, ax=None, ders=(0,), with_fact=True,
         return fig
 
     @add_fig_kwargs
-    def plot_q2vq(self, ax=None, ders=(0,), with_qn=0, with_fact=True,
-                  scale=None, **kwargs) -> Figure:
+    def plot_q2vq(self, ax=None, ders=(0,), with_qn=0, with_fact=True, scale=None, **kwargs) -> Figure:
         r"""
         Plot the local part of the pseudopotential in q space.
 
@@ -297,7 +294,7 @@ def plot_q2vq(self, ax=None, ders=(0,), with_qn=0, with_fact=True,
         color = kwargs.pop("color", "black")
         linewidth = kwargs.pop("linewidth", _LW)
 
-        qmesh, vlspl = self.reader.read_vlspl()
+        qmesh, vlspl = self.r.read_vlspl()
         ecuts = 2 * (np.pi * qmesh)**2
         for atype, vl_atype in enumerate(vlspl):
             for der, values in enumerate(vl_atype):
@@ -315,7 +312,7 @@ def plot_q2vq(self, ax=None, ders=(0,), with_qn=0, with_fact=True,
                 #    ax.plot(ecuts, yvals, color=color, linewidth=linewidth, label="q*f(q) x %2.f" % fact)
 
                 #if der == 0:
-                #   z_val = self.reader.zion_typat[atype]
+                #   z_val = self.r.zion_typat[atype]
                 #   ax.axhline(y=-z_val / np.pi, linewidth=1, color='k', linestyle="dashed")
 
         ax.grid(True)
@@ -355,8 +352,8 @@ def plot_ffspl(self, ax=None, ecut_ffnl=None, ders=(0,), l_select=None,
         l_seen = set()
 
         # vlspl has shape [ntypat, 2, mqgrid_vl]
-        qmesh, vlspl = self.reader.read_vlspl()
-        all_projs = self.reader.read_projectors()
+        qmesh, vlspl = self.r.read_vlspl()
+        all_projs = self.r.read_projectors()
 
         for itypat, projs_type in enumerate(all_projs):
             # Loop over the projectors for this atom type.
diff --git a/abipy/electrons/scissors.py b/abipy/electrons/scissors.py
index 2189497cb..eb6b9149f 100644
--- a/abipy/electrons/scissors.py
+++ b/abipy/electrons/scissors.py
@@ -1,5 +1,7 @@
 # coding: utf-8
 """Scissors operator."""
+from __future__ import annotations
+
 import os
 import numpy as np
 import pickle
@@ -7,6 +9,7 @@
 from collections import OrderedDict
 from monty.collections import AttrDict
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt
+from abipy.tools.typing import Figure
 
 
 __all__ = [
@@ -19,7 +22,7 @@ class ScissorsError(Exception):
     """Base class for the exceptions raised by :class:`Scissors`"""
 
 
-class Scissors(object):
+class Scissors:
     """
     This object represents an energy-dependent scissors operator.
     The operator is defined by a list of domains (energy intervals)
@@ -119,7 +122,7 @@ def apply(self, eig):
         raise self.Error("Cannot find location of eigenvalue %s in domains:\n%s" % (eig, domains))
 
 
-class ScissorsBuilder(object):
+class ScissorsBuilder:
     """
     This object facilitates the creation of :class:`Scissors` instances.
 
@@ -141,7 +144,7 @@ class ScissorsBuilder(object):
     """
 
     @classmethod
-    def from_file(cls, filepath):
+    def from_file(cls, filepath: str):
         """
         Generate object from (SIGRES.nc) file. Main entry point for client code.
         """
@@ -150,7 +153,7 @@ def from_file(cls, filepath):
             return cls(qps_spin=ncfile.qplist_spin, sigres_ebands=ncfile.ebands)
 
     @classmethod
-    def pickle_load(cls, filepath):
+    def pickle_load(cls, filepath: str):
         """Load the object from a pickle file."""
         with open(filepath, "rb") as fh:
             d = AttrDict(pickle.load(fh))
@@ -202,17 +205,17 @@ def __init__(self, qps_spin, sigres_ebands):
         self.build()
 
     @property
-    def nsppol(self):
+    def nsppol(self) -> int:
         """Number of spins."""
         return len(self._qps_spin)
 
     @property
-    def e0min(self):
+    def e0min(self) -> float:
         """Minimum KS energy in eV (takes into account spin)"""
         return self._e0min
 
     @property
-    def e0max(self):
+    def e0max(self) -> float:
         """Maximum KS energy in eV (takes into account spin)"""
         return self._e0max
 
@@ -275,7 +278,7 @@ def build(self, domains_spin=None, bounds_spin=None, k=3):
         return domains_spin
 
     @add_fig_kwargs
-    def plot_qpe_vs_e0(self, with_fields="all", **kwargs):
+    def plot_qpe_vs_e0(self, with_fields="all", **kwargs) -> Figure:
         """Plot the quasiparticle corrections as function of the KS energy."""
         ax_list = None
         for spin, qps in enumerate(self._qps_spin):
@@ -286,7 +289,7 @@ def plot_qpe_vs_e0(self, with_fields="all", **kwargs):
         return fig
 
     @add_fig_kwargs
-    def plot_fit(self, ax=None, fontsize=8, **kwargs):
+    def plot_fit(self, ax=None, fontsize=8, **kwargs) -> Figure:
         """
         Compare fit functions with input quasi-particle corrections.
 
diff --git a/abipy/electrons/tests/test_psps.py b/abipy/electrons/tests/test_psps.py
index d6df52ef5..427dad69e 100644
--- a/abipy/electrons/tests/test_psps.py
+++ b/abipy/electrons/tests/test_psps.py
@@ -10,19 +10,19 @@
 class PspsFileTestCase(AbipyTest):
 
     def test_psps_nc_silicon(self):
-        """Test PSPS.nc file with Ga.oncvpsp"""
+        """Testing PSPS.nc file with Ga.oncvpsp"""
         pseudo = abidata.pseudo("Ga.oncvpsp")
 
         with pseudo.open_pspsfile(ecut=10) as psps:
             repr(psps); print(psps)
-            r = psps.reader
+            r = psps.r
             assert r.usepaw == 0 and r.ntypat == 1
             assert not psps.params
 
             robot = PspsRobot.from_files([psps.filepath])
             repr(psps); print(psps)
 
-            all_projs = psps.reader.read_projectors()
+            all_projs = psps.r.read_projectors()
             for itypat, p_list in enumerate(all_projs):
                 for p in p_list:
                     assert p.to_string(verbose=1)
@@ -30,7 +30,11 @@ def test_psps_nc_silicon(self):
             if self.has_matplotlib():
                 # psps plots.
                 assert psps.plot(what="all", with_qn=True, show=False)
-
+                assert psps.plot_tcore_rspace(ax=None, ders=(0, 1, 2, 3), scale=1.0, rmax=3.0, show=False)
+                assert psps.plot_tcore_qspace(ax=None, ders=(0,), with_fact=True, with_qn=0, scale=1.0, show=False)
+                assert psps.plot_q2vq(ax=None, ders=(0,), with_qn=0, with_fact=True, scale=None, show=False)
+                assert psps.plot_ffspl(ax=None, ecut_ffnl=None, ders=(0,), l_select=None,
+                                       with_qn=0, with_fact=False, scale=None, show=False)
                 # robot plots.
                 assert robot.plot_tcore_rspace(ders=(0, 1, 2, 3), with_qn=0, scale=None, fontsize=8, show=False)
                 assert robot.plot_tcore_qspace(ders=(0, 1), with_qn=0, scale=None, fontsize=8, show=False)
diff --git a/abipy/eph/tests/test_varpeq.py b/abipy/eph/tests/test_varpeq.py
index a7b6f62ea..cff9679df 100644
--- a/abipy/eph/tests/test_varpeq.py
+++ b/abipy/eph/tests/test_varpeq.py
@@ -5,29 +5,42 @@
 from abipy.eph.varpeq import VarpeqFile
 
 
-#class VarpeqTest(AbipyTest):
-#
-#    def test_varpeq_file(self):
-#        """Testing VarpeqFile."""
-#        with VarpeqFile(abidata.ref_file("abinitio_qpath_V1QAVG.nc")) as varpeq:
-#            repr(varpeq)
-#            str(varpeq)
-#            assert varpeq.to_string(verbose=2)
-#            assert not varpeq.params
-#            assert varpeq.structure.formula == "Ga1 P1" and len(varpeq.structure) == 2
-#
-#            #print(varpeq.ebands.kpoints.ksampling)
-#            for polaron in varpeq.polaron_spin:
-#                #print(polaron)
-#                #polaron.plot_bz_sampling()
-#                #polaron.plot_bqnu_with_phbands(phbands)
-#                polaron.plot_ank_with_ebands(ebands)
-#
-#                if self.has_matplotlib():
-#
-#            # Test jupyter notebook creation
-#            if self.has_nbformat():
-#                varpeq.write_notebook(nbpath=self.get_tmpname(text=True))
+class VarpeqTest(AbipyTest):
+
+    def test_varpeq_file(self):
+        """Testing VarpeqFile."""
+        varpeq_filepath = "/Users/giantomassi/git_repos/abinit/_build/tests/POLARON/varpeq6/out_VARPEQ.nc"
+        with VarpeqFile(varpeq_filepath) as varpeq:
+            repr(varpeq)
+            str(varpeq)
+            assert varpeq.to_string(verbose=2)
+            assert varpeq.structure.formula == "Li1 F1" and len(varpeq.structure) == 2
+            params = varpeq.params
+            #assert params["nkbz"] ==
+            #assert params["ngkpt"] ==
+
+            #print(varpeq.ebands.kpoints.ksampling)
+            for polaron in varpeq.polaron_spin:
+                print(polaron)
+                #assert polaron.spin == 0
+                #assert polaron.nstates == 0
+                #assert polaron.nb == 0
+                #assert polaron.nk == 0
+                #assert polaron.nq == 0
+                #assert polaron.bstart == 0
+                #assert polaron.bstop == 0
+                df = polaron.get_final_results_df(with_params=True)
+                #polaron.write_a2_bxsf(self, filepath: PathLike, fill_value: float = 0.0) -> None:
+                #polaron.write_b2_bxsf(self, filepath: PathLike, fill_value: float = 0.0) -> None:
+                if self.has_matplotlib():
+                    polaron.plot_scf_cycle(show=False)
+                    #polaron.plot_ank_with_ebands(ebands_kpath, ebands_kmesh=None)
+                    #polaron.plot_bqnu_with_ddb("in_DDB", with_phdos=True)
+                    #polaron.plot_bqnu_with_phbands(phbands_qpath)
+
+            # Test jupyter notebook creation
+            #if self.has_nbformat():
+            #    varpeq.write_notebook(nbpath=self.get_tmpname(text=True))
 
 
 #class VarpeqRobotTest(AbipyTest):
diff --git a/abipy/eph/varpeq.py b/abipy/eph/varpeq.py
index 1049ebf41..4b5c1f433 100644
--- a/abipy/eph/varpeq.py
+++ b/abipy/eph/varpeq.py
@@ -106,6 +106,7 @@ class Entry:
 # Convert to dictionary: name --> Entry
 _ALL_ENTRIES = {e.name: e for e in _ALL_ENTRIES}
 
+
 class VarpeqFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter):
     """
     This file stores the results of a VARPEQ calculations: SCF cycle, A_nk, B_qnu coefficients
@@ -118,7 +119,9 @@ class VarpeqFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter)
         from abipy.eph.varpeq import VarpeqFile
         with VarpeqFile("out_VARPEQ.nc") as varpeq:
             print(varpeq)
-            varpeq.plot_scf_cycle()
+            for polaron in varpeq.polaron_spin:
+                print(polaron)
+                polaron.plot_scf_cycle()
 
     .. rubric:: Inheritance Diagram
     .. inheritance-diagram:: VarpeqFile
diff --git a/abipy/examples/flows/run_conducwork.py b/abipy/examples/flows/run_conducwork.py
index 604e5aa26..ff31a5022 100755
--- a/abipy/examples/flows/run_conducwork.py
+++ b/abipy/examples/flows/run_conducwork.py
@@ -15,7 +15,7 @@
 from abipy.abio.factories import conduc_kerange_from_inputs
 
 
-def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), **variables):
+def make_scf_input(structure, pseudos, ngkpt=(2, 2, 2), shiftk=(0, 0, 0), **variables):
     """Build and return SCF input given the structure and pseudopotentials"""
 
     scf_inp = abilab.AbinitInput(structure, pseudos=pseudos)
diff --git a/abipy/examples/flows/run_qha_2d.py b/abipy/examples/flows/run_qha_2d.py
index 59fca9ad0..6f0f4e0e4 100755
--- a/abipy/examples/flows/run_qha_2d.py
+++ b/abipy/examples/flows/run_qha_2d.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 r"""
-Flow for quasi-harmonic calculations under development
-======================================================
+Flow for QHA calculations with 2 DOFs
+=====================================
 Warning: This code is still under development.
 """
 import sys
@@ -16,7 +16,7 @@
 
 def build_flow(options):
     """
-    Create a `QhaFlow` for quasi-harmonic calculations.
+    Create a `Qha2dFlow` for quasi-harmonic calculations with 2 DOFs
     """
     # Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
     if not options.workdir:
@@ -43,7 +43,7 @@ def build_flow(options):
    0.0000000000    0.0000000000    9.7234377918
 """)
 
-    # Use NC PBE pseudos from pseudodojo v0.4
+    # Use NC PBEsol pseudos from pseudodojo v0.4
     from abipy.flowtk.psrepos import get_oncvpsp_pseudos
     pseudos = get_oncvpsp_pseudos(xc_name="PBEsol", version="0.4")
 
diff --git a/abipy/examples/plot/plot_qha_2d.py b/abipy/examples/plot/plot_qha_2d.py
index 9793274bf..ec9dffcbd 100755
--- a/abipy/examples/plot/plot_qha_2d.py
+++ b/abipy/examples/plot/plot_qha_2d.py
@@ -34,14 +34,16 @@
 
 qha = QHA_2D.from_files(gsr_paths, phdos_paths, bo_strains_ac, phdos_strains_ac, gsr_file="DDB")
 
+tstart, tstop, num = 0, 1000, 101
+
 #%%
 qha.plot_energies()
 
 #%%
-qha.plot_free_energies(tstop=500, tstart=0, num=6)
+qha.plot_free_energies(tstart=tstart, tstop=500, num=6)
 
 #%%
-qha.plot_thermal_expansion(tstop=1000, tstart=0, num=101)
+qha.plot_thermal_expansion(tstart=tstart, tstop=tstop, num=num)
 
 #%%
-qha.plot_lattice(tstop=1000, tstart=0, num=101)
+qha.plot_lattice(tstart=tstart, tstop=tstop, num=num)
diff --git a/abipy/examples/plot/plot_qha_vzsisa.py b/abipy/examples/plot/plot_qha_vzsisa.py
index ad4bfb3bd..e13351c88 100755
--- a/abipy/examples/plot/plot_qha_vzsisa.py
+++ b/abipy/examples/plot/plot_qha_vzsisa.py
@@ -24,48 +24,48 @@
 phdos_paths = [os.path.join(root, "scale_{:d}_PHDOS.nc".format(s)) for s in strains2]
 
 qha = Vzsisa.from_ddb_phdos_files(ddb_paths, phdos_paths)
-tstart, tstop = 0, 800
+tstart, tstop, num = 0, 800, 101
 
 #%%
 # Plot BO Energies as a function of volume for different T
-qha.plot_bo_energies(tstop=tstop, tstart=tstart, num=11)
+qha.plot_bo_energies(tstart=tstart, tstop=tstop, num=11)
 
 #%%
 # Plot Volume as a function of T
-qha.plot_vol_vs_t(tstop=tstop, tstart=tstart, num=101)
+qha.plot_vol_vs_t(tstart=tstart, tstop=tstop, num=num)
 
 #%%
 # Plot Lattice as a function of T
-qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101)
+qha.plot_abc_vs_t(tstart=tstart, tstop=tstop, num=num)
 
 #%%
-# Plot Lattice as a function of T")
-qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, lattice="b")
+# Plot Lattice as a function of T
+qha.plot_abc_vs_t(tstart=tstart, tstop=tstop, num=num, lattice="b")
 
 #%%
 # Plot Volumetric thermal expansion coefficient as a function of T
-qha.plot_thermal_expansion_coeff(tstop=tstop, tstart=tstart ,num=101)
+qha.plot_thermal_expansion_coeff(tstart=tstart, tstop=tstop, num=num)
 
 #%%
 # Plot Thermal expansion coefficient as a function of T
-qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart ,num=101)
+qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart, num=num)
 
 #%%
 # Plot Angles as a function of T
-qha.plot_angles_vs_t(tstop=tstop, tstart=tstart, num=101)
+qha.plot_angles_vs_t(tstart=tstart, tstop=tstop, num=num)
 
 #%%
 #
 # Plot Volume as a function of T. 4th order polinomial
-qha.plot_vol_vs_t_4th(tstop=tstop, tstart=tstart, num=101)
+qha.plot_vol_vs_t_4th(tstart=tstart, tstop=tstop, num=num)
 
 #%%
 # Plot Lattice as a function of T. 4th order polinomial
-qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, num=101, lattice="a")
+qha.plot_abc_vs_t_4th(tstart=tstart, tstop=tstop, num=num, lattice="a")
 
 #%%
 # Plot Lattice as a function of T. 4th order polinomial
-qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart)
+qha.plot_abc_vs_t_4th(tstart=tstart, tstop=tstop)
 
 #%%
 # Plot Volumetric thermal expansion coefficient as a function of T
@@ -73,12 +73,11 @@
 
 #%%
 # Plot Thermal expansion coefficient as a function of T
-qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293)
+qha.plot_thermal_expansion_coeff_abc_4th(tstart=tstart, tstop=tstop, num=num, tref=293)
 
 #%%
 # Plot Angles as a function of T.
-qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3)
-
+qha.plot_angles_vs_t_4th(tstart=tstart, tstop=tstop, num=num, angle=3)
 
 #%%
 # Create plotter to plot all the phonon DOS.
diff --git a/abipy/flowtk/vzsisa.py b/abipy/flowtk/vzsisa.py
index 30f3fc792..52b3f918f 100644
--- a/abipy/flowtk/vzsisa.py
+++ b/abipy/flowtk/vzsisa.py
@@ -54,7 +54,7 @@ def from_scf_input(cls,
                 This option is the recommended one if the k-path contains two consecutive high symmetry k-points
                 that are very close as ndivsm > 0 may produce a very large number of wavevectors.
                 if 0, deactivate band structure calculation for electrons.
-            edos_ngkpt: Three integers defining the the k-sampling for the computation of the
+            edos_ngkpt: Three integers defining the k-sampling for the computation of the
                 electron DOS with the relaxed structures. Useful for metals or small gap semiconductors
                 in which the electronic contribution should be included. None disables the computation of the e-DOS.
             manager: |TaskManager| instance. Use default if None.
diff --git a/abipy/iotools/cube.py b/abipy/iotools/cube.py
index c67c94d7d..a0410be58 100644
--- a/abipy/iotools/cube.py
+++ b/abipy/iotools/cube.py
@@ -75,7 +75,7 @@ def cube_read_structure_mesh_data(filepath: str) -> tuple:
                 data[ii // (ny * nz), (ii // nz) % ny, ii % nz] = float(val)
                 ii += 1
         data = data / (bohr_to_angstrom ** 3)
-        if ii != nx*ny*nz:
+        if ii != nx * ny * nz:
             raise ValueError('Wrong number of data points ...')
         from abipy.core.structure import Structure
         structure = Structure.from_sites(sites=sites)
diff --git a/abipy/iotools/visualizer.py b/abipy/iotools/visualizer.py
index 27f329bf3..0930931b0 100644
--- a/abipy/iotools/visualizer.py
+++ b/abipy/iotools/visualizer.py
@@ -21,7 +21,7 @@ def is_macosx() -> bool:
     return "darwin" in sys.platform
 
 
-def find_loc(app_name) -> str:
+def find_loc(app_name: str) -> str:
     """
     Returns the location of the application from its name. None if not found.
     """
@@ -31,7 +31,7 @@ def find_loc(app_name) -> str:
     return path
 
 
-def _find_loc(app_name) -> str:  # pragma: no cover
+def _find_loc(app_name: str) -> str:  # pragma: no cover
     # Try command line version
     path = which(app_name)
     if path is not None: return path
@@ -92,11 +92,11 @@ def __str__(self):
         return "%s: %s, is_macosx_app %s, filepath: %s" % (
             self.__class__.__name__, self.binpath, self.is_macosx_app, self.filepath)
 
-    def __call__(self):  # pragma: no cover
+    def __call__(self) -> int:  # pragma: no cover
         """
         Call the visualizer in a subprocess.
 
-        Returns: exit status of the subprocess.
+        Return: exit status of the subprocess.
         """
         from subprocess import call
         if not self.is_macosx_app:
@@ -160,12 +160,12 @@ def from_file(cls, filepath: str):
         # Get the file extension.
         root, ext = os.path.splitext(filepath)
         if not ext:
-            raise ValueError("Cannot detect file extension in %s " % filepath)
+            raise ValueError(f"Cannot detect file extension in {filepath}")
 
         avail_visus = cls.get_available(ext=ext)
 
         if not avail_visus:
-            raise cls.Error("Cannot find available visualizer for extension %s " % ext)
+            raise cls.Error(f"Cannot find available visualizer for extension {ext}")
 
         return avail_visus[0](filepath)
 
@@ -180,7 +180,7 @@ def from_name(cls, appname):
         for visu in cls.__subclasses__():
             if visu.name == appname: return visu
 
-        raise cls.Error("`%s` is not among the list of supported visualizers" % appname)
+        raise cls.Error(f"{appname=} is not among the list of supported visualizers")
 
     @classmethod
     def all_visunames(cls):
diff --git a/abipy/tools/plotting.py b/abipy/tools/plotting.py
index 32f8844e0..205887df9 100644
--- a/abipy/tools/plotting.py
+++ b/abipy/tools/plotting.py
@@ -751,7 +751,7 @@ def __init__(self, xkey: str, xs: VectorLike, yvals_dict: dict[str, VectorLike],
             if not duck.is_listlike(ytols): ytols = [ytols]
 
             if any(yt <= 0 for yt in ytols):
-                raise ValueError(f"tolerances cannot be negative: {ytols}")
+                raise ValueError(f"tolerances cannot be negative: {ytols=}")
 
             # Sort input tolerances just to be on the safe side.
             self.ytols_dict[ykey] = np.sort(np.array(ytols))[::-1]
@@ -841,7 +841,7 @@ def get_dataframe_ykey(self, ykey: str) -> pd.DataFrame:
             #rows.append(dict(ytol=ytol, ix=ix, xx=xx), xx_best=xx_best, ykey=ykey)
         return pd.DataFrame(rows)
 
-    def to_string(self, verbose=0) -> str:
+    def to_string(self, verbose: int = 0) -> str:
         """
         String representation with verbosity level `verbose`.
         """
@@ -1131,9 +1131,9 @@ class Exposer:
     """
 
     @classmethod
-    def as_exposer(cls, exposer, **kwargs) -> Exposer:
+    def as_exposer(cls, exposer: str | Exposer, **kwargs) -> Exposer:
         """
-        Return an instance of Exposer, usually from a string with then name.
+        Return an instance of Exposer, usually from a string with the name.
 
         Args:
             exposer: "mpl" for MplExposer, "panel" for PanelExposer.
@@ -1741,9 +1741,9 @@ def __init__(self, py_row: int, py_col: int, nrows: int, ncols: int):
         self.iax = 1 + self.py_col + self.py_row * self.ncols
         # Note that plotly col and row start from 1.
         if nrows == 1 and ncols == 1:
-            self.ply_row, self.ply_col = (None, None)
+            self.ply_row, self.ply_col = None, None
         else:
-            self.ply_row, self.ply_col = (self.py_row + 1, self.py_col + 1)
+            self.ply_row, self.ply_col = self.py_row + 1, self.py_col + 1
 
     def __str__(self) -> str:
         lines = []
@@ -2888,16 +2888,3 @@ def plot(self, deg_list: list[int],
         ax.legend(loc="best", fontsize=fontsize, shadow=True)
 
         return fig
-
-
-#class PolyExtrapolator:
-#    """
-#    Fit data with polynomals, extrapolate to zero and visualize the results.
-#    """
-#    def __init__(self, xs, ys):
-#        self.xs, self.ys = np.array(xs), np.array(ys)
-#
-#    def extrapolate_to_zero(self, deg: int):
-
-
-
diff --git a/abipy/wannier90/abiwan.py b/abipy/wannier90/abiwan.py
index 27e697f0f..1a4aedd4d 100644
--- a/abipy/wannier90/abiwan.py
+++ b/abipy/wannier90/abiwan.py
@@ -45,7 +45,7 @@ def from_file(cls, filepath: str) -> AbiwanFile:
 
     def __init__(self, filepath: str):
         super().__init__(filepath)
-        self.reader = self.r = AbiwanReader(filepath)
+        self.r = AbiwanReader(filepath)
 
         # Number of bands actually used to construct the Wannier functions
         self.num_bands_spin = self.r.read_value("num_bands")
@@ -80,7 +80,7 @@ def bands_in(self) -> np.ndarray:
     def lwindow(self) -> np.ndarray:
         """
         [nsppol, nkpt, max_num_bands] array. Only if disentanglement.
-        True if this band at this k-point lies within the outer window
+        True if this band at this k-point lies within the outer window.
         """
         return self.r.read_value("lwindow_int").astype(bool)
 
@@ -450,7 +450,7 @@ def __init__(self, structure, nwan_spin, spin_vmatrix, spin_rmn, irvec, ndegen):
         self.nband = nwan_spin[0]
         #self.nelect
 
-    def eval_sk(self, spin, kpt, der1=None, der2=None) -> np.ndarray:
+    def eval_sk(self, spin: int, kpt, der1=None, der2=None) -> np.ndarray:
         """
         Interpolate eigenvalues for all bands at a given (spin, k-point).
         Optionally compute gradients and Hessian matrices.
@@ -465,7 +465,7 @@ def eval_sk(self, spin, kpt, der1=None, der2=None) -> np.ndarray:
             oeigs[nband]
         """
         if der1 is not None or der2 is not None:
-            raise NotImplementedError("Derivatives")
+            raise NotImplementedError("Derivatives are not coded")
 
         # O_ij(k) = sum_R e^{+ik.R}*O_ij(R)
         j2pi = 2.0j * np.pi
@@ -538,7 +538,7 @@ class AbiwanRobot(Robot, RobotWithEbands):
     """
     EXT = "ABIWAN"
 
-    def get_dataframe(self, with_geo=True, abspath=False, funcs=None, **kwargs) -> pd.DataFrame:
+    def get_dataframe(self, with_geo: bool = True, abspath: bool = False, funcs=None, **kwargs) -> pd.DataFrame:
         """
         Return a |pandas-DataFrame| with the most important Wannier90 results and the filenames as index.
 
diff --git a/abipy/wannier90/wout.py b/abipy/wannier90/wout.py
index 7d4b8d723..6ca701690 100644
--- a/abipy/wannier90/wout.py
+++ b/abipy/wannier90/wout.py
@@ -7,6 +7,7 @@
 
 from collections import OrderedDict
 from monty.string import marquee
+from monty.termcolor import cprint
 from abipy.core.mixins import BaseFile, Has_Structure, NotebookWriter
 from abipy.core.structure import Structure
 from abipy.tools.plotting import add_fig_kwargs, get_axarray_fig_plt
@@ -345,7 +346,7 @@ def plot_centers_spread(self, fontsize=8, **kwargs) -> Figure:
         Returns: |matplotlib-Figure|
         """
         if self._parse_iterations() != 0:
-            print("Wout files does not contain Wannierization cycles. Returning None")
+            cprint("Wout files does not contain Wannierization cycles. Returning None", color="red")
             return None
 
         # Build grid of plots.
diff --git a/docs/api/abipy.tools.rst b/docs/api/abipy.tools.rst
index 4224736f8..383c549e3 100644
--- a/docs/api/abipy.tools.rst
+++ b/docs/api/abipy.tools.rst
@@ -105,6 +105,14 @@ abipy.tools.numtools module
    :undoc-members:
    :show-inheritance:
 
+abipy.tools.pade module
+-----------------------
+
+.. automodule:: abipy.tools.pade
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 abipy.tools.parallel module
 ---------------------------