diff --git a/abipy/abilab.py b/abipy/abilab.py
index 97afc638d..754689327 100644
--- a/abipy/abilab.py
+++ b/abipy/abilab.py
@@ -423,7 +423,7 @@ def cmp_version(this: str, other: str, op: str = ">=") -> bool:
         Compare two version strings with the given operator ``op``
         >>> assert cmp_version("1.1.1", "1.1.0") and not cmp_version("1.1.1", "1.1.0", op="==")
         """
-        from pkg_resources import parse_version
+        from packaging.version import parse as parse_version
         from monty.operator import operator_from_str
         op = operator_from_str(op)
         return op(parse_version(this), parse_version(other))
diff --git a/abipy/abio/abivar_database/variables_abinit.py b/abipy/abio/abivar_database/variables_abinit.py
index b75d3c014..5e687ffbc 100644
--- a/abipy/abio/abivar_database/variables_abinit.py
+++ b/abipy/abio/abivar_database/variables_abinit.py
@@ -1603,6 +1603,39 @@
 """,
 ),
 
+Variable(
+    abivarname="chkparal",
+    varset="gstate",
+    vartype="integer",
+    topics=['parallelism_useful'],
+    dimensions="scalar",
+    defaultval=1,
+    mnemonics="CHecK whether the PARALelism is adequate",
+    added_in_version="before_v9",
+    text=r"""
+Not all parallelism types or level are allowed or simply relevant for the different [[optdriver]] values in ABINIT.
+It has been observed that some users do not understand well their relation. In particular, their expectation of the adequacy
+of some parallelism for some [[optdriver]] value was not correct, with a large loss of computing resources.
+Indeed, if the user does not sufficiently understand the parallelism in ABINIT, huge amount of ressources might be spend
+when they are booked for a run that cannot use these.
+Accordingly, the user might blame ABINIT for being slow while the user has simply not activated
+the relevant parallelism, or activated an irrelevant parallelism.
+
+However, if the user correctly understand the parallelism,
+it might be more convenient to leave in the input file irrelevant variables. This is especially the case for high-throughput calculations
+driven by workflows developed for earlier versions of ABINIT.
+
+The default value of [[chkparal]], will enforce some basic relevance of the input variables related to parallelism,
+thus hopefully preventing some users to loose computing power.
+
+The following relevances and adequacies are checked at present if [[chkparal]]=1 :
+the input variable [[autoparal]] is relevant only for [[optdriver]]=1 calculations (ground-state);
+the input variable [[paral_kgb]] is relevant only for [[optdriver]]=1 calculations (ground-state) or for [[optdriver]]=66 (Laczos-Sternheimer GW).
+
+The relevance of [[paral_atom]] or [[paral_rf]] or [[gwpara]] is not checked at present. The default values should not yield loss of computing power.
+""",
+),
+
 Variable(
     abivarname="chkprim",
     varset="gstate",
@@ -2475,7 +2508,7 @@
     characteristics=['[[ENERGY]]'],
     added_in_version="before_v9",
     text=r"""
-It is the value of the "scissors operator", the shift of conduction band
+It is the value of the "scissors operator" [[cite:Levine1989]], the shift of conduction band
 eigenvalues, used in response function calculations.
 Can be specified in Ha (the default), Ry, eV or Kelvin, since [[ecut]] has the
 [[ENERGY]] characteristics (1 Ha = 27.2113845 eV).
@@ -3346,6 +3379,9 @@
 Compute weight of configuration computed during CTQMC calculations.
 For example, for a calculation on $d$ orbitals, the calculations
 gives the weight of 0,1,2,3,4,5,6,7,8,9 and 10 electrons configurations.
+  * 0 --> Nothing done
+  * 1 --> Calculation of weight of configurations
+  * 2 --> Calculation of local magnetic susceptibilty
 """,
 ),
 
@@ -4118,7 +4154,7 @@
 zero-point renormalisation as well as temperature dependence.
 If set to 1 (default), use the non-adiabatic version of the Allen-Heine-Cardona equation to compute the
 zero-point renormalisation as well as temperature dependence.
-Note: The use of [[eph_ahc_type]]=0 is not recommanded in IR-active materials.
+Note: The use of [[eph_ahc_type]]=0 is not recommended in IR-active materials.
 """,
 ),
 
@@ -4353,6 +4389,7 @@
          Save results in GSTORE.nc file (requires netcdf library with MPI-IO support).
          The k-mesh must be equal to the one associated to the input WFK file, the q-mesh is specified
          by [[eph_ngqpt_fine]] (NB: the q-mesh must be a sub-mesh of the k-mesh or equal).
+* 14 --> Compute the molecular Berry curvature from GSTORE.nc. No support for metals or non-collinear magnetism yet. Reference: [[cite:Saparov2022]], [[cite:Coh2023]].
 * 15, -15 --> Write the average in r-space of the DFPT potentials to the V1QAVG.nc file.
               In the first case (+15) the q-points are specified via [[ph_nqpath]] and [[ph_qpath]]. The code assumes the
               input DVDB contains q-points in the IBZ and the potentials along the path
@@ -4469,7 +4506,8 @@
     mnemonics="EXPERTise of the USER",
     added_in_version="9.2.2",
     text=r"""
-If set to 0, the checking provided by ABINIT is maximum (default values of [[chkprim]], [[chksymbreak]], [[chksymtnons]], [[chkdilatmx]]).
+If set to 0, the checking provided by ABINIT is maximum (default values of [[chkprim]], [[chksymbreak]], [[chksymtnons]], [[chkdilatmx]],
+[[chkparal]]).
 If non-zero (up to three), the above-mentioned checking input variables are all disabled (set to zero)
 although it is still possible to activate particular tests by specifying input variables directly in the input file.
 In the future, the level three will always be the maximum allowed value, with all checks set to zero,
@@ -5961,7 +5999,7 @@
     produced in a previous ground state calculation
   * This variable or [[irdwfkfine]] is mandatory when [[bs_interp_mode]] == 1
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 
 **This variable is experimental. In development.**
 """,
@@ -6265,31 +6303,7 @@
   The [[KOKKOS]] GPU implementation can be used in conjuction with openMP threads
   on CPU (see [[gpu_kokkos_nthrd]]).
 
-For an expert use of ABINIT on [[GPU]], some additional keywords can be used. See
-[[gpu_use_nvtx]],[[gpu_nl_distrib]],[[gpu_nl_splitsize]].
-""",
-),
-
-Variable(
-    abivarname="gpu_use_nvtx",
-    varset="paral",
-    vartype="integer",
-    topics=['parallelism_expert'],
-    dimensions="scalar",
-    defaultval=0,
-    mnemonics="GPU: activate USE of NVTX tracing/profiling",
-    added_in_version="9.7.2",
-    text=r"""
-Only available if ABINIT executable has been linked with the [[NVTX]] (Nvidia) or
-[[ROCTX]] (AMD) library.
-When [[gpu_use_nvtx]] = 1 on Nvidia GPU(s),
-the use of [NVTX](https://github.com/NVIDIA/NVTX)
- tracing/profiling ((NVIDIA Tools eXtension Library) is activated.
- The trace can be used for instance with
- [NVidia NSight System](https://developer.nvidia.com/nsight-systems)
-When [[gpu_use_nvtx]] = 1 on AMD GPU(s),
-the use of [ROCTX](https://rocm.docs.amd.com/projects/roctracer/en/latest/roctracer_spec.html)
- tracing/profiling (ROCm Tools eXtension Library) is activated.
+For an expert use of ABINIT on [[GPU]], some additional keywords can be used. See [[gpu_nl_distrib]], [[gpu_nl_splitsize]].
 """,
 ),
 
@@ -7802,18 +7816,18 @@
 ),
 
 Variable(
-    abivarname="invol_blk_sliced",
+    abivarname="invovl_blksliced",
     varset="dev",
     vartype="integer",
     topics=['parallelism_expert'],
     dimensions="scalar",
-    mnemonics="INVerse OverLap: BLocK-diagonal matrix applied SLICED",
+    mnemonics="INVerse OVerLap: BLocK-diagonal matrix applied SLICED",
     added_in_version="9.7.2",
     defaultval=ValueWithConditions({'[[gpu_option]] > 0': '0', 'defaultval': 1}),
     text=r"""
 Only relevant if [[wfoptalg]] == 1 or 111 (WF optimization by Chebyshev filtering algorithm).
 In the Chebyshev-filtered subspace method (iterative diagonalization algorithm))
-one needs to apply the inverse of the overlap matrix. [[invol_blk_sliced]] allows one
+one needs to apply the inverse of the overlap matrix. [[invovl_blksliced]] allows one
 to choose between two variants, sliced (1) or non-sliced (0).
 Default value is different for an execution on GPU.
 """,
@@ -8222,7 +8236,7 @@
 When [[iscf]] < 0, the reading of a DEN file is always enforced.
 
 A non-zero value of **ird1den** is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8261,7 +8275,7 @@
   * if [[irdddk]] = 1: read first-order ddk wavefunctions from a disk file appended with _1WFx,
     produced in a previous response function calculation.
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8344,7 +8358,7 @@
 the charge neutrality sum rule is usually prohibitively large.
 
 A non-zero value of [[irdddb]] is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 
 The input variable [[getddb]] is an alternative to [[irdddb]], in the multidataset case.
 Note also that, starting Abinit v9, one can also use [[getddb_filepath]] to specify the path of the DDB file directly.
@@ -8372,7 +8386,7 @@
   * if [[irdddk]] = 1: read first-order ddk wavefunctions from a disk file appended with _1WFx,
     produced in a previous response function calculation
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8391,7 +8405,7 @@
 Alternative to [[getden_filepath]] and [[getden]].
 
 A non-zero value of [[irdden]] is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8412,7 +8426,7 @@
 Alternative to [[getkden]].
 
 A non-zero value of [[irdkden]] is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8492,7 +8506,7 @@
 [[irdsuscep]] or [[getsuscep]]) must be non-zero.
 
 A non-zero value of [[irdscr]] is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8513,7 +8527,7 @@
 (alternatively, [[irdscr]] or [[getscr]]) must be non-zero.
 
 A non-zero value of [[irdsuscep]] is treated in the same way as other "ird" variables.
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8577,7 +8591,7 @@
   * if [[irdddk]] = 1: read first-order ddk wavefunctions from a disk file appended with _1WFx,
     produced in a previous response function calculation
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -8614,7 +8628,7 @@
      produced in a previous ground state calculation
   * This variable or [[getwfkfine]] is mandatory when [[bs_interp_mode]] = 1
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 
 **This variable is experimental. In development.**
 """,
@@ -10358,7 +10372,7 @@
     requires="[[optdriver]] in [3,4,7,99]",
     added_in_version="before_v9",
     text=r"""
-The scissor operator energy added to the conductions states. In some cases,
+The scissor operator energy added to the conductions states [[cite:Levine1989]]. In some cases,
 it mimics a second iteration self-consistent GW calculation.
 """,
 ),
@@ -11015,6 +11029,7 @@
     will give a band structure with the first 8 bands converged within [[tolwfr]].
     For more complex systems and/or GS NSCF calculations with many empty states, one usually needs
     to increase [[nbdbuf]], let's say 10% of [[nband]].
+    This can be easily achieved by using a negative value e.g.: nbdbuf -10 means 10% of nband.
 
 
 This concept is useful in three situations: in non-self-consistent calculations, for the
@@ -12669,7 +12684,7 @@
 will be appended by _WFQ instead of _WFK
 Also, if 1 and a RF calculation is done, defines the wavevector of the perturbation.
 
-For further information about the *files file*, consult the [[help:abinit#files-file]].
+For further information about the naming of files in ABINIT, consult the [[help:abinit#files-file]].
 """,
 ),
 
@@ -13072,7 +13087,7 @@
     dimensions=[3, '[[natom]]'],
     defaultval=0.0,
     mnemonics="NUClear DIPole MOMents",
-    requires="[[usepaw]] = 1; [[pawcpxocc]] = 2; [[kptopt]] = 0 or 3",
+    requires="[[usepaw]] = 1; [[pawcpxocc]] = 2; [[optforces]] = 0; [[optstress]] = 0; [[kptopt]] = 0, 4, or 3",
     added_in_version="before_v9",
     text=r"""
 Places an array of nuclear magnetic dipole moments on the atomic
@@ -13095,15 +13110,15 @@
     dimensions="scalar",
     defaultval=0,
     mnemonics="NUClear site Electric Field Gradient",
-    requires="[[usepaw]] == 1, [[quadmom]]",
+    requires="[[usepaw]] == 1",
     added_in_version="before_v9",
     text=r"""
 If nonzero, calculate the electric field gradient at each atomic site in the unit cell.
-Using this option requires [[quadmom]] to be set as well.
 Values will be written to the main output file (search for Electric Field Gradient).
-If nucefg=1, only the quadrupole coupling in MHz and asymmetry are reported.
-If nucefg=2, the full electric field gradient tensors in atomic units are also given,
+If nucefg=1, the electric field gradient in both atomic units and SI units is given,
 showing separate contributions from the valence electrons, the ion cores, and the PAW reconstruction.
+If nucefg=2, the quadrupole couplings to the nuclear electric quadrupole moments are reported as well,
+and are based on the input values of [[quadmom]].
 If nucefg=3, then in addition to the nucefg=2 output, the EFGs are computed using an ionic point charge model.
 This is useful for comparing the accurate PAW-based results to those of simple ion-only models.
 Use of nucefg=3 requires that the variable [[ptcharge]] be set as well.
@@ -13799,6 +13814,7 @@
     characteristics=['[[DEVELOP]]'],
     requires="""[[usepaw]] == 1;
 [[usexcnhat]] == 0;
+[[pawxcdev]] == 0;
 [[paral_atom]] == 0;
 [[paral_kgb]] == 0;
 ([[kptopt]] == 3 or [[kptopt]] == 0) """,
@@ -13807,14 +13823,17 @@
 Compute quantities related to orbital magnetic moment. Typically used in the
 presence of a nonzero nuclear magnetic dipole moment, see [[nucdipmom]], to compute
 the nuclear magnetic shielding as measured in NMR. [[orbmag]]
-    is parallelized over k points only. The implementation follows the
-    theory outlined in [[cite:Zwanziger2023]].
-    The computed results are returned in the
-    standard output file, search for "Orbital magnetic moment". This calculation requires
-    both the ground state and DDK wavefunctions (see [[rfddk]]), and is triggered at the end of a
-    DDK calculation.
-
-* [[orbmag]] = 1: Compute orbital magnetization and Chern number
+is parallelized over k points only. The implementation follows the
+theory outlined in [[cite:Zwanziger2023]]. The computed results are returned in the
+standard output file, search for "Orbital magnetic moment". This calculation requires
+both the ground state and DDK wavefunctions (see [[rfddk]] or [[berryopt]]). The
+preferred way to use [[orbmag]] is at the end of a DFPT DDK calculation. Alternatively, it
+can be called in a ground state calculation if [[berryopt]] -2 has also been called,
+to generate discretized DDK wavefunctions. Note that convergence with kpt mesh is
+*much* faster using the DFPT approach, and the [[berryopt]] approach is not recommended
+unless a very specific ground state feature is also needed.
+
+* [[orbmag]] = 1: Compute orbital magnetization and Chern vector
 * [[orbmag]] = 2: Same as [[orbmag]] 1 but also print out values of each term making up total
 orbital magnetic moment and a band-by-band decomposition.
 """,
@@ -13883,8 +13902,11 @@
 Relevant only for PAW calculations.
 This keyword controls the parallel distribution of memory over atomic sites.
 Calculations are also distributed using the "kpt-band" communicator.
-Compatible with ground-state calculations and response function calculations
-""",
+Compatible with ground-state calculations and response function calculations.
+
+This parallelization concerns only a small part of the whole calculation in the sequential case.
+When using parallelism, it might be that this small part becomes predominant if [[paral_atom]] is not activated.
+"""
 ),
 
 Variable(
@@ -17255,15 +17277,14 @@
     dimensions=['[[ntypat]]'],
     defaultval=MultipleValue(number=None, value=0),
     mnemonics="QUADrupole MOMents",
-    requires="[[usepaw]] == 1 and [[nucefg]]>=1",
+    requires="[[usepaw]] == 1 and [[nucefg]]>1",
     added_in_version="before_v9",
     text=r"""
   * Array of quadrupole moments, in barns, of the nuclei. These values are used
   in conjunction with the electric field gradients computed with [[nucefg]] to
   calculate the quadrupole couplings in MHz, as well as the asymmetries. Note that
   the electric field gradient at a nuclear site is independent of the nuclear
-  quadrupole moment, thus the quadrupole moment of a nucleus can be input as 0,
-  and the option [[nucefg]] = 2 used to determine the electric field gradient at the site.
+  quadrupole moment, thus the quadrupole moment of a nucleus can be input as 0.
 """,
 ),
 
@@ -17856,7 +17877,7 @@
     electric field are different, one often does the calculation of derivatives in
     a separate dataset, followed by calculation of electric field response as well as phonon.
     The options 2 and 3 proves useful in that context; also, in case a scissor
-    shift is to be used, it is usually not applied for the $\,d / \,d k$ response).
+    shift is to be used, it is usually not applied for the $\,d / \,d k$ response [[cite:Levine1989]]).
 """,
 ),
 
@@ -18865,8 +18886,8 @@
     mnemonics="STRess FACTor",
     added_in_version="before_v9",
     text=r"""
-The stresses multiplied by [[strfact]] will be treated like forces in the
-process of optimization ([[ionmov]] = 2 or 22, non-zero [[optcell]]).
+The stresses (in atomic units) multiplied by [[strfact]] will be treated like forces (in atomic units) in the
+algorithms for optimization ([[ionmov]] = 2 or 22, non-zero [[optcell]]).
 For example, the stopping criterion defined by [[tolmxf]] relates to these
 scaled stresses.
 """,
@@ -20312,7 +20333,7 @@
     abivarname="usepotzero",
     varset="dev",
     vartype="integer",
-    topics=['Coulomb_useful'],
+    topics=['Coulomb_useful','Verification_useful'],
     dimensions="scalar",
     defaultval=0,
     mnemonics="USE POTential ZERO",
@@ -20324,7 +20345,7 @@
 
   * [[usepotzero]] = 1, the new convention: the all-electron physical potential is set to zero average value.
 
-  * [[usepotzero]] = 2, the PWscf convention: the potential of equivalent point charges is set to
+  * [[usepotzero]] = 2, the QE/PWscf convention: the potential of equivalent point charges is set to
   zero average value (convention also valid for NC pseudopotentials).
 """,
 ),
@@ -21252,18 +21273,18 @@
     text=r"""
 Selects a van-der-Waals density functional to apply the corresponding
 correction to the exchange-correlation energy. If set to zero, no correction
-will be applied.
+will be applied. For response functions, only [[vdw_xc]]=5, 6 and 7 have been implemented (obviously also 0)..
 Possible values are:
 
   * 0: no correction.
-  * 1: apply vdW-DF1 (DRSLL) from [[cite:Dion2004]].
-  * 2: apply vdw-DF2 (LMKLL) from [[cite:Lee2010]].
-  * 5: apply vdw-DFT-D2 as proposed by S. Grimme [[cite:Grimme2006]] (adding a semi-empirical dispersion potential). Available only for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence.
-  * 6: apply vdw-DFT-D3 as proposed by S. Grimme [[cite:Grimme2010]] (refined version of DFT-D2). Available only for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence and [[vdw_tol_3bt]] variable to include 3-body corrections.
-  * 7: apply vdw-DFT-D3(BJ) as proposed by Grimme (based on Becke-Jonhson method from [[cite:Becke2006]]). Available only for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence.
-  * 10: evaluate the vdW correlation energy from maximally localized Wannier functions, as proposed by P. L. Silvestrelli, also known as vdW-WF1 method [[cite:Silvestrelli2008]]. For details on this implementation please check [[cite:Espejo2012]]. The improvements introduced by Andrinopoulos _et al._ [[cite:Andrinopoulos2011]], namely the amalgamation procedure, splitting of p-like MLWFs into two s-like Wannier functions and fractional occupation of MLWFs are performed automatically.
-  * 11: evaluate the vdW correlation energy from maximally localized Wannier functions, as proposed by A. Ambrosetti and P. L. Silvestrelli, also known as vdW-WF2 method [[cite:Ambrosetti2012]].
-  * 14: apply DFT/vdW-QHO-WF method as proposed by Silvestrelli, which combines the quantum harmonic oscillator-model with localized Wannier functions [[cite:Silvestrelli2013]]. For periodic systems a supercell approach has to be used since **vdw_supercell** is not enabled in this case.
+  * 1: apply vdW-DF1 (DRSLL) from [[cite:Dion2004]]. Not available for response functions.
+  * 2: apply vdw-DF2 (LMKLL) from [[cite:Lee2010]]. Not available for response functions.
+  * 5: apply vdw-DFT-D2 as proposed by S. Grimme [[cite:Grimme2006]] (adding a semi-empirical dispersion potential). Available for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence.
+  * 6: apply vdw-DFT-D3 as proposed by S. Grimme [[cite:Grimme2010]] (refined version of DFT-D2). Available for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence and [[vdw_tol_3bt]] variable to include 3-body corrections.
+  * 7: apply vdw-DFT-D3(BJ) as proposed by Grimme (based on Becke-Jonhson method from [[cite:Becke2006]]). Available for ground-state calculations and response functions; see [[vdw_tol]] variable to control convergence.
+  * 10: evaluate the vdW correlation energy from maximally localized Wannier functions, as proposed by P. L. Silvestrelli, also known as vdW-WF1 method [[cite:Silvestrelli2008]]. For details on this implementation please check [[cite:Espejo2012]]. The improvements introduced by Andrinopoulos _et al._ [[cite:Andrinopoulos2011]], namely the amalgamation procedure, splitting of p-like MLWFs into two s-like Wannier functions and fractional occupation of MLWFs are performed automatically. Not available for response functions.
+  * 11: evaluate the vdW correlation energy from maximally localized Wannier functions, as proposed by A. Ambrosetti and P. L. Silvestrelli, also known as vdW-WF2 method [[cite:Ambrosetti2012]]. Not available for response functions.
+  * 14: apply DFT/vdW-QHO-WF method as proposed by Silvestrelli, which combines the quantum harmonic oscillator-model with localized Wannier functions [[cite:Silvestrelli2013]]. For periodic systems a supercell approach has to be used since **vdw_supercell** is not enabled in this case. Not available for response functions.
 
 For [[vdw_xc]] = 1 and [[vdw_xc]] = 2, the implementation follows the strategy
 devised in the article of Roman-Perez and Soler [[cite:Romanperez2009]].
@@ -22271,6 +22292,22 @@
      Find k-points inside (electron/hole) pockets according to the values specified by [[sigma_erange]].
      Write KERANGE.nc file with all the tables required by the code to automate NSCF band structure calculations
      inside the pocket(s) and electron lifetime computation in the EPH code when [[eph_task]] = -4.
+
+  * "wannier" --> Read WFK file and run Wannierization. It has the similar effect of
+      [[prtwant]] = 2, which uses the **ABINIT- Wannier90** interface. The difference is that with wfk_task "wannier",
+      the $\kk$-points in the full BZ is not necessary. Instead, the wavefunctions with the $\kk$-points not in
+      the IBZ will be reconstructed by symmetry. This functionality does not yet work with PAW when the wavefunction
+      is not already in full BZ.
+
+      ABINIT will produce the input files required by Wannier90 and it will run
+      Wannier90 to produce the Maximally-locallized Wannier functions (see [
+      http://www.wannier.org ](http://www.wannier.org) ).
+      !!! Notes
+
+          * The files that are created can also be used by Wannier90 in stand-alone mode.
+          * In order to use Wannier90 as a post-processing program for ABINIT you might have to
+            compile it with the appropriate flags (see ABINIT makefile). You might use ./configure --enable-wannier90
+          * There are some other variables related to the interface of Wannier90 and ABINIT. See [[varset:w90]].
 """,
 ),
 
@@ -23607,6 +23644,30 @@
   * **useextfpmd** = 4, the energy shift will be evaluated
 by making an integration of the trial potential over the real space and the
 contributions will be computed with integrals over the band number.
+
+  * **useextfpmd** = 5, the energy shift will be evaluated
+by making an integration of the trial potential over the real space and the
+contributions will be computed with discrete integer sums over the band indices.
+[[extfpmd_nband]] sets the last band index, which acts like [[nband]] for a
+conventional calculation.
+""",
+),
+
+Variable(
+    abivarname="extfpmd_nband",
+    varset="gstate",
+    vartype="integer",
+    topics=['ExtFPMD_basic'],
+    dimensions="scalar",
+    defaultval=0,
+    mnemonics="EXTended FPMD: Number of Bands",
+    added_in_version="10.1.0",
+    text=r"""
+Specifies the number of bands to use for extended FPMD contributions when using [[useextfpmd]] = 5.
+This acts like [[nband]] for a conventional calculation. **extfpmd_nband** must be sufficiently high
+so that its occupancy is close to zero. Extended FPMD contributions will be computed from [[nband]] to [[extfpmd_nband]].
+
+**extfpmd_nband** must be greater than [[nband]].
 """,
 ),
 
@@ -23817,6 +23878,26 @@
 """,
 ),
 
+Variable(
+    abivarname="gstore_gmode",
+    varset="eph",
+    vartype="string",
+    topics=['ElPhonInt_basic'],
+    dimensions="scalar",
+    defaultval="phonon",
+    mnemonics=r"GSTORE GMODE",
+    requires="[[optdriver]] == 7",
+    added_in_version="10.1.2",
+    text=r"""
+This input variable specifies the representation used to store the e-ph matrix elements in the GSTORE.nc file
+
+Possible values are:
+
+    "phonon" --> Store e-ph matrix elements in the phonon representation (collective displacement)
+    "atom" -->  Store e-ph matrix elements in the atom representation (displacement of a single atom along one of the reduced directions)
+""",
+),
+
 Variable(
     abivarname="gstore_brange",
     varset="eph",
@@ -23861,6 +23942,67 @@
 """,
 ),
 
+Variable(
+    abivarname="getabiwan_filepath",
+    varset="eph",
+    vartype="string",
+    topics=['ElPhonInt_basic'],
+    dimensions="scalar",
+    defaultval="None",
+    mnemonics="GET the ABIWAN.nc from FILEPATH",
+    added_in_version="10.1.1",
+    text=r"""
+This variable defines the path of the ABIWAN.nc file with the Wannier rotation matrices produced by ABINIT
+when wannier90 in invoked in library mode.
+See also [[getabiwan]].
+""",
+),
+
+Variable(
+    abivarname="getabiwan",
+    varset="eph",
+    vartype="int",
+    topics=['ElPhonInt_basic'],
+    dimensions="scalar",
+    defaultval="None",
+    mnemonics="GET the ABIWAN.nc from dataset",
+    added_in_version="10.1.1",
+    text=r"""
+This variable is similar in spirit to [[getabiwan_filepath]] but uses the dataset index
+instead of the filepath.
+""",
+),
+
+Variable(
+    abivarname="getgwan_filepath",
+    varset="eph",
+    vartype="string",
+    topics=['ElPhonInt_basic'],
+    dimensions="scalar",
+    defaultval="None",
+    mnemonics="GET the GWAN.nc from FILEPATH",
+    added_in_version="10.1.1",
+    text=r"""
+This variable defines the path of the GWAN.nc file with the e-ph matrix elements in the Wannier representation
+See also [[getgwan]].
+""",
+),
+
+Variable(
+    abivarname="getgwan",
+    varset="eph",
+    vartype="int",
+    topics=['ElPhonInt_basic'],
+    dimensions="scalar",
+    defaultval="None",
+    mnemonics="GET the GWAN.nc from dataset",
+    added_in_version="10.1.1",
+    text=r"""
+This variable is similar in spirit to [[getgwan_filepath]] but uses the dataset index
+instead of the filepath.
+""",
+),
+
 Variable(
     abivarname="gstore_erange",
     varset="eph",
@@ -23963,8 +24105,9 @@
 * EGW0 -->
 * G0EW -->
 * RPA_ENERGY --> Compute RPA correlation energy within the ACFDT framework.
-* CC4S -->
+* CC4S --> Read density from file, diagonalize the KS Hamiltonian and produce output files required by CC4S.
 * CC4S_FULL --> Same as CC4S but compute maximum number of eigenvectors according to [[ecut]]
+* CC4S_FROM_WFK --> Same as CC4S but read single particle orbitals from an external WFK file.
 
 !!! important
 
diff --git a/abipy/abio/abivar_database/variables_optic.py b/abipy/abio/abivar_database/variables_optic.py
index e2f8bfcbe..b71b91835 100644
--- a/abipy/abio/abivar_database/variables_optic.py
+++ b/abipy/abio/abivar_database/variables_optic.py
@@ -237,7 +237,7 @@
 This parameter provides a fixed shift to all the conduction bands. As
 LDA/GGA are known to underestimate the band-gap by a significant amount in
 some cases, in order to obtain a reasonable optical spectrum and make a realistic
-comparison with experiments one needs to correct for this.
+comparison with experiments one needs to correct for this [[cite:Levine1989]].
 The scissors shift is normally chosen to be the difference between the experimental and
 theoretical band-gap, and simply shifts the conduction bands. Alternatively, one may
 determine the self energy using the [[tutorial:gw1|GW approach]], in which case
diff --git a/abipy/abio/inputs.py b/abipy/abio/inputs.py
index f5a953559..a4ded6dc7 100644
--- a/abipy/abio/inputs.py
+++ b/abipy/abio/inputs.py
@@ -472,7 +472,7 @@ def enforce_znucl_and_typat(self, znucl, typat):
         if len(typat) != len(self.structure):
             raise ValueError("typat contains %d entries while it should be natom: %d" % (len(typat), len(self.structure)))
 
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
         if len(znucl) != ntypat:
             raise ValueError("znucl contains %d entries while it should be ntypat: %d" % (len(znucl), ntypat))
 
@@ -814,7 +814,7 @@ def escape(text):
                 for name in names:
                     value = vars[name]
                     if mnemonics and value is not None:
-                        print(f"{name=}")
+                        #print(f"{name=}")
                         app(escape("#### <" + var_database[name].mnemonics + ">"))
 
                     # Build variable, convert to string and append it
@@ -4053,6 +4053,7 @@ def to_string(self, sortmode=None, mode="text", verbose=0, files_file=False) ->
             #if mode == "html": vname = root + "#%s" % vname
             if mode == "html": vname = var_database[vname].html_link(label=vname)
             value = format_string_abivars(vname, value, "anaddb")
+            #print("vname:", vname, "value:", value)
             app(str(InputVariable(vname, value)))
 
         return "\n".join(lines) if mode == "text" else "\n".join(lines).replace("\n", "<br>")
diff --git a/abipy/core/kpoints.py b/abipy/core/kpoints.py
index 673c6ac44..de69b4878 100644
--- a/abipy/core/kpoints.py
+++ b/abipy/core/kpoints.py
@@ -101,7 +101,6 @@ def issamek(k1, k2, atol=None):
     """
     k1 = np.asarray(k1)
     k2 = np.asarray(k2)
-    #if k1.shape != k2.shape:
 
     return is_integer(k1 - k2, atol=atol)
 
@@ -409,32 +408,35 @@ def map_kpoints(other_kpoints, other_lattice, ref_lattice, ref_kpoints, ref_symr
 def kpoints_indices(frac_coords, ngkpt, check_mesh=0) -> np.ndarray:
     """
     This function is used when we need to insert k-dependent quantities in a (nx, ny, nz) array.
-    It compute the indices of the k-points assuming these points belong to
-    a mesh with ngkpt divisions.
+    It computes the indices of the k-points assuming these points belong to a mesh with ngkpt divisions.
 
     Args:
         frac_coords
         ngkpt:
-        check_mesh
+        check_mesh:
     """
+
     # Transforms kpt in its corresponding reduced number in the interval [0,1[
     k_indices = [np.round((kpt % 1) * ngkpt) for kpt in frac_coords]
-
     k_indices = np.array(k_indices, dtype=int)
-    #k_indices = np.array(list(map(tuple, k_indices)))
-    #k_indices = list(map(tuple, k_indices))
 
+    # Debug secction.
     if check_mesh:
         print(f"kpoints_indices: Testing whether k-points belong to the {ngkpt =} mesh")
         ierr = 0
         for kpt, inds in zip(frac_coords, k_indices):
-            #print("kpt:", kpt, "inds:", inds)
+            if check_mesh > 1: print("kpt:", kpt, "inds:", inds)
             same_k = np.array((inds[0]/ngkpt[0], inds[1]/ngkpt[1], inds[2]/ngkpt[2]))
             if not issamek(kpt, same_k):
                 ierr += 1; print(kpt, "-->", same_k)
         if ierr:
             raise ValueError("Wrong mapping")
-        print("Check succesful!")
+
+        #for kpt, inds in zip(frac_coords, k_indices):
+        #    if np.any(inds >= ngkpt):
+        #        raise ValueError(f"inds >= nkgpt for {kpt=}, {np.round(kpt % 1)=} {inds=})")
+
+        print("Check succesfull!")
 
     return k_indices
 
@@ -526,6 +528,7 @@ def kpath_from_bounds_and_ndivsm(bounds, ndivsm, structure):
         for j in range(ndivs[i]):
             p = bounds[i] + j * (bounds[i + 1] - bounds[i]) / ndivs[i]
             path.append(p)
+
     path.append(bounds[-1])
 
     return np.array(path)
@@ -1400,6 +1403,8 @@ def lines(self):
             for line in self.lines:
                 vals_on_line = eigens[spin, line, band]
         """
+        if len(self) < 2:
+            return tuple()
         prev = self.versors[0]
         lines = [[0]]
 
diff --git a/abipy/core/skw.py b/abipy/core/skw.py
index 39f33c26b..2a2a3dc19 100644
--- a/abipy/core/skw.py
+++ b/abipy/core/skw.py
@@ -890,8 +890,8 @@ def __init__(self, lpratio, kpts, eigens, fermie, nelect, cell, symrel, has_timr
         cprint("FIT vs input data: Mean Absolute Error= %.3e (meV)" % mae, color="red" if warn else "green")
         if warn:
             # Issue warning if error too large.
-            cprint("Large error in SKW interpolation!", "red")
-            cprint("MAE:", mae, "[meV]", "red")
+            cprint("Large error in SKW interpolation!", color="red")
+            cprint(f"MAE: {mae} [meV]", color="red")
 
         self.mae = mae
 
diff --git a/abipy/core/structure.py b/abipy/core/structure.py
index d65982af6..328280f5e 100644
--- a/abipy/core/structure.py
+++ b/abipy/core/structure.py
@@ -1272,6 +1272,50 @@ def print_neighbors(self, radius=2.0) -> None:
                 print("\t", repr(s), " at distance", dist)
             print("")
 
+    @lazy_property
+    def has_zero_dynamical_quadrupoles(self):
+        """
+        Dynamical quadrupoles are nonzero in all noncentrosymmetric crystals,
+        but also in centrosymmetric ones if one or more atoms are placed at noncentrosymmetric sites.
+        """
+        def create_image(s1, s2):
+            """
+            Creates the image of s2 through s1
+            This image is a fictitious atom in the structure
+            """
+            image = PeriodicSite.from_dict(s2.as_dict())
+            image.coords = s1.coords - (s2.coords - s1.coords)
+            image.frac_coords = s1.frac_coords - (s2.frac_coords - s1.frac_coords)
+
+            return image
+
+        def check_image(structure, site):
+            """
+            Checks if a fictitious site is an image of a site of the structure
+            """
+            for site in structure.sites:
+                if site.is_periodic_image(site):
+                    return True
+
+            return False
+
+        # If the centrosymmetry is broken at a given atomic site of the given structure,
+        # returns False. Else, return True
+        sites = self.sites
+
+        for s1 in sites:
+            for s2 in sites:
+                # Do not take s1 = s2 into account
+                if s2 != s1:
+                    # Create the image of s2 through s1 (fictitious atom)
+                    image = create_image(s1, s2)
+                    is_image = check_image(self, image)
+                    if not is_image:
+                        return False
+
+        return True
+
+
     @lazy_property
     def hsym_kpath(self):
         """
@@ -1892,8 +1936,7 @@ def get_smallest_supercell(self, qpoint, max_supercell):
         Returns: the scaling matrix of the supercell
         """
         if np.allclose(qpoint, 0.0):
-            scale_matrix = np.eye(3, 3)
-            return scale_matrix
+            return np.eye(3, 3)
 
         l = max_supercell
 
@@ -1966,11 +2009,12 @@ def get_smallest_supercell(self, qpoint, max_supercell):
             raise ValueError('max_supercell is not large enough for this q-point')
 
         # Fortran 2 python!!!
-        return scale_matrix.T
+        return scale_matrix.T.copy()
 
-    def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
+    def make_doped_supercells(self, scaling_matrix, replaced_atom, dopant_atom):
         """
         Returns a list doped supercell structures, one for each non-equivalent site of the replaced atom.
+
         Args:
             scaling_matrix: A scaling matrix for transforming the lattice vectors.
                 Has to be all integers. Several options are possible:
@@ -1985,18 +2029,18 @@ def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
             dopant_atom : Symbol of the dopant_atom (ex: 'Eu')
         """
         ### list of positions of non-equivalent sites for the replaced atom. ###
-        irred=self.spget_equivalent_atoms().eqmap # mapping from inequivalent sites to atoms sites
-        positions=self.get_symbol2indices()[replaced_atom] # get indices of the replaced atom
+        irred = self.spget_equivalent_atoms().eqmap # mapping from inequivalent sites to atoms sites
+        positions = self.get_symbol2indices()[replaced_atom] # get indices of the replaced atom
 
-        index_non_eq_sites=[]
+        index_non_eq_sites = []
         for pos in positions:
             if len(irred[pos]) != 0:
                  index_non_eq_sites.append(irred[pos][0])
 
-        doped_supercell=self.copy()
+        doped_supercell = self.copy()
         doped_supercell.make_supercell(scaling_matrix)
 
-        doped_structure_list=[]
+        doped_structure_list = []
 
         for index in index_non_eq_sites:
             final_structure=doped_supercell.copy()
@@ -2005,8 +2049,6 @@ def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
 
         return doped_structure_list
 
-
-
     def get_trans_vect(self, scale_matrix):
         """
         Returns the translation vectors for a given scale matrix.
@@ -2038,8 +2080,7 @@ def range_vec(i):
 
         # find the translation vectors (in terms of the initial lattice vectors)
         # that are inside the unit cell defined by the scale matrix
-        # we're using a slightly offset interval from 0 to 1 to avoid numerical
-        # precision issues
+        # we're using a slightly offset interval from 0 to 1 to avoid numerical precision issues
         #print(scale_matrix)
         inv_matrix = np.linalg.inv(scale_matrix)
 
@@ -2106,8 +2147,8 @@ def write_vib_file(self, xyz_file, qpoint, displ, do_real=True, frac_coords=True
     def frozen_2phonon(self, qpoint, displ1, displ2, eta=1, frac_coords=False, scale_matrix=None, max_supercell=None):
         """
         Creates the supercell needed for a given qpoint and adds the displacements.
-        The displacements are normalized so that the largest atomic displacement will correspond to the
-        value of eta in Angstrom.
+        The displacements are normalized so that the largest atomic displacement will correspond
+        to the value of eta in Angstrom.
 
         Args:
             qpoint: q vector in reduced coordinate in reciprocal space.
@@ -2182,15 +2223,14 @@ def frozen_phonon(self, qpoint, displ, eta=1, frac_coords=False, scale_matrix=No
         value of eta in Angstrom.
 
         Args:
-            qpoint: q vector in reduced coordinate in reciprocal space.
+            qpoint: q-vector in reduced coordinate in reciprocal space.
             displ: displacement in real space of the atoms.
-            eta: pre-factor multiplying the displacement. Gives the value in Angstrom of the
-                largest displacement.
+            eta: pre-factor multiplying the displacement. Gives the value in Angstrom of the largest displacement.
             frac_coords: whether the displacements are given in fractional or cartesian coordinates
             scale_matrix: the scaling matrix of the supercell. If None a scaling matrix suitable for
                 the qpoint will be determined.
             max_supercell: mandatory if scale_matrix is None, ignored otherwise. Defines the largest
-                supercell in the search for a scaling matrix suitable for the q point.
+                supercell in the search for a scaling matrix suitable for the input q-point.
 
         Returns:
             A namedtuple with a Structure with the displaced atoms, a numpy array containing the
@@ -2199,35 +2239,34 @@ def frozen_phonon(self, qpoint, displ, eta=1, frac_coords=False, scale_matrix=No
 
         if scale_matrix is None:
             if max_supercell is None:
-                raise ValueError("If scale_matrix is not provided, please provide max_supercell!")
-
+                raise ValueError("If scale_matrix is not provided in input, please provide max_supercell!")
             scale_matrix = self.get_smallest_supercell(qpoint, max_supercell=max_supercell)
 
         scale_matrix = np.array(scale_matrix, np.int16)
         if scale_matrix.shape != (3, 3):
             scale_matrix = np.array(scale_matrix * np.eye(3), np.int16)
-        print("scale_matrix:", scale_matrix)
+        #print("scale_matrix:\n", scale_matrix)
 
         old_lattice = self._lattice
         new_lattice = Lattice(np.dot(scale_matrix, old_lattice.matrix))
 
         tvects = self.get_trans_vect(scale_matrix)
-        print("tvects", tvects)
+        #print("tvects\n", tvects)
 
         if frac_coords:
-            displ = np.array((old_lattice.get_cartesian_coords(d) for d in displ))
+            displ = np.array([old_lattice.get_cartesian_coords(d) for d in displ])
         else:
             displ = np.array(displ)
-        # from here displ are in cartesian coordinates
 
+        # from here on, displ is in cartesian coordinates.
         displ = eta * displ / np.linalg.norm(displ, axis=1).max()
 
         new_displ = np.zeros(3, dtype=float)
         new_sites = []
         displ_list = []
-        for at, site in enumerate(self):
+        for iat, site in enumerate(self):
             for t in tvects:
-                new_displ[:] = np.real(np.exp(2*1j*np.pi*(np.dot(qpoint,t)))*displ[at,:])
+                new_displ[:] = np.real(np.exp(2*1j*np.pi * (np.dot(qpoint, t))) * displ[iat,:])
                 displ_list.append(list(new_displ))
 
                 coords = site.coords + old_lattice.get_cartesian_coords(t) + new_displ
@@ -2699,7 +2738,7 @@ def structure2siesta(structure: Structure, verbose=0) -> str:
     lines = []
     app = lines.append
     app("NumberOfAtoms %d" % len(structure))
-    app("NumberOfSpecies %d" % structure.ntypesp)
+    app("NumberOfSpecies %d" % structure.n_elems)
 
     if verbose:
         app("# The species number followed by the atomic number, and then by the desired label")
diff --git a/abipy/core/testing.py b/abipy/core/testing.py
index d44f423c4..c863a4471 100644
--- a/abipy/core/testing.py
+++ b/abipy/core/testing.py
@@ -39,7 +39,7 @@ def cmp_version(this: str, other: str, op: str = ">=") -> bool:
     Compare two version strings with the given operator ``op``
     >>> assert cmp_version("1.1.1", "1.1.0") and not cmp_version("1.1.1", "1.1.0", op="==")
     """
-    from pkg_resources import parse_version
+    from packaging.version import parse as parse_version
     from monty.operator import operator_from_str
     op = operator_from_str(op)
     return op(parse_version(this), parse_version(other))
diff --git a/abipy/core/tests/test_structure.py b/abipy/core/tests/test_structure.py
index 4f5473ea3..e908e87bf 100644
--- a/abipy/core/tests/test_structure.py
+++ b/abipy/core/tests/test_structure.py
@@ -397,7 +397,7 @@ def test_frozen_phonon_methods(self):
         f2p_data = structure.frozen_2phonon(qpoint, 0.05 * displ, 0.02*displ2, eta=0.5, frac_coords=False,
                                             max_supercell=mx_sc, scale_matrix=scale_matrix)
 
-        d_tot = 0.05*displ+0.02*displ2
+        d_tot = 0.05 *displ + 0.02 * displ2
         max_displ = np.linalg.norm(d_tot, axis=1).max()
         self.assert_almost_equal(f2p_data.structure[0].coords,
                                     structure[0].coords + 0.5*d_tot[0]/max_displ)
diff --git a/abipy/dfpt/ddb.py b/abipy/dfpt/ddb.py
index 3dd80cecc..74fe7943d 100644
--- a/abipy/dfpt/ddb.py
+++ b/abipy/dfpt/ddb.py
@@ -1966,9 +1966,9 @@ def _run_anaddb_task(self, anaddb_input, mpi_procs, workdir, manager, verbose):
 
         return task
 
-    def write(self, filepath, filter_blocks=None) -> None:
+    def write(self, filepath: str, filter_blocks=None) -> None:
         """
-        Writes the DDB file in filepath. Requires the blocks data.
+        Writes the DDB file to filepath. Requires the blocks data.
         Only the information stored in self.header.lines and in self.blocks
         are be used to produce the file
         """
@@ -2087,8 +2087,8 @@ def remove_block(self, dord, qpt=None, qpt3=None):
 
     def get_2nd_ord_dict(self) -> dict:
         """
-        Generates an ordered dictionary with the second order derivative of the form
-        {qpt: {(idir1, ipert1, idir2, ipert2): complex value}}.
+        Generates an ordered dictionary with the second order derivatives in the form
+            {qpt: {(idir1, ipert1, idir2, ipert2): complex value}}.
 
         Returns:
             OrderedDict: a dictionary with all the elements of a dynamical matrix
@@ -2106,8 +2106,7 @@ def get_2nd_ord_dict(self) -> dict:
 
     def set_2nd_ord_data(self, data, replace=True) -> None:
         """
-        Insert the blocks corresponding to the data provided for the second
-        order perturbations.
+        Insert the blocks corresponding to the data provided for the second order perturbations.
 
         Args:
             data: a dict of the form {qpt: {(idir1, ipert1, idir2, ipert2): complex value}}.
@@ -2226,7 +2225,7 @@ def write_notebook(self, nbpath=None) -> str:
 
 class Becs(Has_Structure, MSONable):
     """
-    This object stores the Born effective charges and provides simple tools for data analysis.
+    This object stores the Born effective charges and provides tools for data analysis.
     """
 
     @pmg_serialize
@@ -2741,8 +2740,6 @@ def plot_e0w_qdirs(self, qdirs=None, w_min=0, w_max=None, gamma_ev=1e-4, num=500
             with_phfreqs: True to show phonon frequencies with dots.
             ax: |matplotlib-Axes| or None if a new figure should be created.
             fontsize: Legend and label fontsize.
-
-        Return: |matplotlib-Figure|
         """
         wmesh = self._get_wmesh(gamma_ev, num, units, w_min, w_max)
         t = np.zeros((num, 3, 3), dtype=complex)
@@ -2830,7 +2827,7 @@ def _add_phfreqs_plotly(self, fig, rcd, units, with_phfreqs) -> None:
     def reflectivity(self, qdir, w, gamma_ev=1e-4, units='eV'):
         """
         Calculates the reflectivity from the dielectric tensor along the specified direction
-        according to eq. (58) in :cite:`Gonze1997` PRB55, 10355 (1997).
+        according to Eq. (58) in :cite:`Gonze1997` PRB55, 10355 (1997).
 
         Args:
             qdir: a list with three components defining the direction.
@@ -2841,8 +2838,7 @@ def reflectivity(self, qdir, w, gamma_ev=1e-4, units='eV'):
             units: string specifying the units used for phonon frequencies. Possible values in
                 ("eV", "meV", "Ha", "cm-1", "Thz"). Case-insensitive.
 
-        Returns:
-            the value of the reflectivity.
+        Returns: the value of the reflectivity.
         """
 
         qdir = np.array(qdir) / np.linalg.norm(qdir)
@@ -2874,8 +2870,6 @@ def plot_reflectivity(self, qdirs=None, w_min=0, w_max=None, gamma_ev=1e-4, num=
             with_phfreqs: True to show phonon frequencies with dots.
             ax: |matplotlib-Axes| or None if a new figure should be created.
             fontsize: Legend and label fontsize.
-
-        Return: |matplotlib-Figure|
         """
         wmesh = self._get_wmesh(gamma_ev, num, units, w_min, w_max)
         t = np.zeros((num, 3, 3), dtype=complex)
diff --git a/abipy/dfpt/phonons.py b/abipy/dfpt/phonons.py
index 6fe50ed0d..64205b8f0 100644
--- a/abipy/dfpt/phonons.py
+++ b/abipy/dfpt/phonons.py
@@ -1123,7 +1123,8 @@ def plot(self, ax=None, units="eV", qlabels=None, branch_range=None, match_bands
         self.plot_ax(ax, branch_range, units=units, match_bands=match_bands, **kwargs)
 
         if points is not None:
-            ax.scatter(points.x, np.array(points.y), s=np.abs(points.s), marker="o", c="b")
+            ax.scatter(points.x, np.array(points.y), s=np.abs(points.s),
+                        marker=points.marker, c=points.color, alpha=points.alpha)
 
         if temp is not None:
             # Scatter plot with Bose-Einstein occupation factors for T = temp
@@ -1626,7 +1627,7 @@ def plot_fatbands(self, use_eigvec=True, units="eV", colormap="jet", phdos_file=
         """
         lw = kwargs.pop("lw", 2)
         factor = abu.phfactor_ev2units(units)
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
 
         # Prepare PJDOS.
         close_phdos_file = False
@@ -1767,7 +1768,7 @@ def plotly_fatbands(self, use_eigvec=True, units="eV", colormap="G10", phdos_fil
         """
         lw = kwargs.pop("lw", 2)
         factor = abu.phfactor_ev2units(units)
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
 
         # Prepare PJDOS.
         close_phdos_file = False
@@ -2058,7 +2059,7 @@ def plot_phdispl(self, qpoint, cart_dir=None, use_reduced_coords=False, ax=None,
 
         ax, fig, plt = get_ax_fig_plt(ax=ax)
         cmap = plt.get_cmap(colormap)
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
 
         if is_non_analytical_direction:
             ax.set_title("q-direction = %s" % repr(qpoint), fontsize=fontsize)
@@ -2940,7 +2941,7 @@ def __init__(self, filepath: str):
             path: path to the file
         """
         super().__init__(filepath)
-        self.reader = PHBST_Reader(filepath)
+        self.reader = self.r = PHBST_Reader(filepath)
 
         # Initialize Phonon bands and add metadata from ncfile
         self._phbands = PhononBands.from_file(filepath)
@@ -2982,7 +2983,7 @@ def phbands(self) -> PhononBands:
 
     def close(self) -> None:
         """Close the file."""
-        self.reader.close()
+        self.r.close()
 
     @lazy_property
     def params(self) -> dict:
@@ -3623,12 +3624,16 @@ def __init__(self, filepath: str):
         # Open the file, read data and create objects.
         super().__init__(filepath)
 
-        self.reader = r = PhdosReader(filepath)
-        self.wmesh = r.wmesh
+        self.reader = self.r = PhdosReader(filepath)
+        self.wmesh = self.r.wmesh
 
     def close(self) -> None:
         """Close the file."""
-        self.reader.close()
+        self.r.close()
+
+    @lazy_property
+    def qptrlatt(self):
+        return self.r.read_value("qptrlatt")
 
     @lazy_property
     def params(self) -> dict:
@@ -3666,12 +3671,12 @@ def to_string(self, verbose: int = 0) -> str:
     @lazy_property
     def structure(self) -> Structure:
         """|Structure| object."""
-        return self.reader.structure
+        return self.r.structure
 
     @lazy_property
     def phdos(self) -> PhononDos:
         """|PhononDos| object."""
-        return self.reader.read_phdos()
+        return self.r.read_phdos()
 
     @lazy_property
     def pjdos_symbol(self):
@@ -3680,7 +3685,7 @@ def pjdos_symbol(self):
         where PhononDos is the contribution to the total DOS summed over atoms
         with chemical symbol `symbol`.
         """
-        return self.reader.read_pjdos_symbol_dict()
+        return self.r.read_pjdos_symbol_dict()
 
     @lazy_property
     def msqd_dos(self):
@@ -3688,7 +3693,7 @@ def msqd_dos(self):
         |MsqDos| object with Mean square displacement tensor in cartesian coords.
         Allows one to calculate Debye Waller factors by integration with 1/omega and the Bose-Einstein factor.
         """
-        return self.reader.read_msq_dos()
+        return self.r.read_msq_dos()
 
     @add_fig_kwargs
     def plot_pjdos_type(self, units="eV", stacked=True, colormap="jet", alpha=0.7, exchange_xy=False,
@@ -3839,7 +3844,7 @@ def plot_pjdos_cartdirs_type(self, units="eV", stacked=True, colormap="jet", alp
         Returns: |matplotlib-Figure|.
         """
         lw = kwargs.pop("lw", 2)
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
         factor = abu.phfactor_ev2units(units)
 
         # Three rows for each direction.
@@ -3851,7 +3856,7 @@ def plot_pjdos_cartdirs_type(self, units="eV", stacked=True, colormap="jet", alp
         cmap = plt.get_cmap(colormap)
 
         # symbol --> [three, number_of_frequencies] in cart dirs
-        pjdos_symbol_rc = self.reader.read_pjdos_symbol_xyz_dict()
+        pjdos_symbol_rc = self.r.read_pjdos_symbol_xyz_dict()
 
         xx = self.phdos.mesh * factor
         for idir, ax in enumerate(ax_list):
@@ -3865,7 +3870,7 @@ def plot_pjdos_cartdirs_type(self, units="eV", stacked=True, colormap="jet", alp
 
             # Plot Type projected DOSes along cartesian direction idir
             cumulative = np.zeros(len(self.wmesh))
-            for itype, symbol in enumerate(self.reader.chemical_symbols):
+            for itype, symbol in enumerate(self.r.chemical_symbols):
                 color = cmap(float(itype) / max(1, ntypat - 1))
                 yy = pjdos_symbol_rc[symbol][idir] / factor
 
@@ -3908,7 +3913,7 @@ def plot_pjdos_cartdirs_site(self, view="inequivalent", units="eV", stacked=True
         """
         # Define num_plots and ax2atom depending on view.
         factor = abu.phfactor_ev2units(units)
-        #natom, ntypat = len(self.structure), self.structure.ntypesp
+        #natom, ntypat = len(self.structure), self.structure.n_elems
         lw = kwargs.pop("lw", 2)
 
         # Select atoms.
@@ -3923,7 +3928,7 @@ def plot_pjdos_cartdirs_site(self, view="inequivalent", units="eV", stacked=True
         cmap = plt.get_cmap(colormap)
 
         # [natom, three, nomega] array with PH-DOS projected over atoms and cartesian directions
-        pjdos_atdir = self.reader.read_pjdos_atdir()
+        pjdos_atdir = self.r.read_pjdos_atdir()
 
         xx = self.phdos.mesh * factor
         for idir, ax in enumerate(ax_list):
@@ -4019,7 +4024,7 @@ def to_pymatgen(self) -> PmgCompletePhononDos:
         total_dos = self.phdos.to_pymatgen()
 
         # [natom, three, nomega] array with PH-DOS projected over atoms and cartesian directions"""
-        pjdos_atdir = self.reader.read_pjdos_atdir()
+        pjdos_atdir = self.r.read_pjdos_atdir()
 
         factor = abu.phfactor_ev2units("thz")
         summed_pjdos = np.sum(pjdos_atdir, axis=1) / factor
diff --git a/abipy/electrons/ebands.py b/abipy/electrons/ebands.py
index 48877ab78..9d0fb1e5e 100644
--- a/abipy/electrons/ebands.py
+++ b/abipy/electrons/ebands.py
@@ -185,7 +185,7 @@ def num_leadblanks(string):
             return _TIPS
 
 
-class ElectronTransition(object):
+class ElectronTransition:
     """
     This object describes an electronic transition between two single-particle states.
     """
@@ -516,8 +516,7 @@ def __init__(self, structure, kpoints, eigens, fermie, occfacts, nelect, nspinor
         Args:
             structure: |Structure| object.
             kpoints: |KpointList| instance.
-            eigens: Array-like object with the eigenvalues (eV) stored as [s, k, b]
-                where s: spin , k: kpoint, b: band index
+            eigens: Array-like object with the eigenvalues (eV) stored as [s, k, b] where s: spin , k: kpoint, b: band index
             fermie: Fermi level in eV.
             occfacts: Occupation factors (same shape as eigens)
             nelect: Number of valence electrons in the unit cell.
@@ -798,8 +797,7 @@ def with_points_along_path(self, frac_bounds=None, knames=None, dist_tol=1e-12):
             frac_bounds: [M, 3] array  with the vertexes of the k-path in reduced coordinates.
                 If None, the k-path is automatically selected from the structure.
             knames: List of strings with the k-point labels defining the k-path. It has precedence over frac_bounds.
-            dist_tol: A point is considered to be on the path if its distance from the line
-                is less than dist_tol.
+            dist_tol: A point is considered to be on the path if its distance from the line is less than dist_tol.
 
         Return:
             namedtuple with the following attributes::
@@ -2216,7 +2214,8 @@ def plot(self, spin=None, band_range=None, klabels=None, e0="fermie", ax=None, y
                 self.plot_ax(ax, e0, spin=spin, band=band, **opts)
 
         if points is not None:
-            ax.scatter(points.x, np.array(points.y) - e0, s=np.abs(points.s), marker="o", c="b")
+            ax.scatter(points.x, np.array(points.y) - e0, s=np.abs(points.s),
+                        marker=points.marker, c=points.color, alpha=points.alpha)
 
         if with_gaps and (self.mband > self.nspinor * self.nelect // 2):
             # Show fundamental and direct gaps for each spin.
@@ -2435,9 +2434,63 @@ def plotly(self, spin=None, band_range=None, klabels=None, e0="fermie", fig=None
 
         return fig
 
+    @add_fig_kwargs
+    def plot_split(self, ylims_list: list,
+                   hspace: float = 0.1, ratio: float = 0.5, **kwargs) -> Figure:
+        """
+        Plot electronic bands with broken y-axis.
+        Based on: https://matplotlib.org/stable/gallery/subplots_axes_and_figures/broken_axis.html
 
-    # TODO: Is this really useful?
+        Args:
+            ylims_list: List of tuples. Each tuple defines the y-limits for the subplots.
+            hspace: space between Axes.
+            ratio: proportion of vertical to horizontal extent of the slanted line.
+            kwargs: keyword arguments passed to self.plot method.
+        """
+        import matplotlib.pyplot as plt
+        num_axis = len(ylims_list)
+
+        fig, ax_list = plt.subplots(num_axis, 1, sharex=True)
+        ax_list = np.flip(ax_list)
+        # Adjust space between Axes
+        fig.subplots_adjust(hspace=hspace)
+
+        for ix, (ax, ylims) in enumerate(zip(ax_list, ylims_list)):
+            # Plot the same data on all Axes
+            self.plot(ax=ax, show=False, **kwargs)
+            # Zoom-in / limit the view to different portions of the data.
+            ax.set_ylim(ylims[0], ylims[1])
+
+            # Hide the spines between axis.
+            if ix == 0 or (num_axis > 2 and ix != num_axis - 1):
+                ax.spines.top.set_visible(False)
+                ax.xaxis.tick_bottom()
+                set_visible(ax, False, *["title"])
+
+            if ix == num_axis - 1 or (num_axis > 2 and ix != num_axis - 1):
+                ax.spines.bottom.set_visible(False)
+                ax.xaxis.tick_top()
+                ax.tick_params(labeltop=False)  # don't put tick labels at the top
+                set_visible(ax, False, *["legend", "xlabel", "yabel"])
+
+        # Now, let's turn towards the cut-out slanted lines.
+        # We create line objects in axes coordinates, in which (0,0), (0,1),
+        # (1,0), and (1,1) are the four corners of the Axes.
+        # The slanted lines themselves are markers at those locations, such that the
+        # lines keep their angle and position, independent of the Axes size or scale
+        # Finally, we need to disable clipping.
+        plt_kws = dict(marker=[(-1, -ratio), (1, ratio)], markersize=12,
+                       linestyle="none", color='k', mec='k', mew=1, clip_on=False)
+
+        for ix, ax in enumerate(ax_list):
+            if ix == 0:
+                ax.plot([0, 1], [1, 1], transform=ax.transAxes, **plt_kws)
+            else:
+                ax.plot([0, 1], [0, 0], transform=ax.transAxes, **plt_kws)
+
+        return fig
 
+    # TODO: Is this really useful?
     @add_fig_kwargs
     def plot_scatter3d(self, band, spin=0, e0="fermie", colormap="jet", ax=None, **kwargs) -> Figure:
         r"""
@@ -2512,6 +2565,7 @@ def decorate_plotly(self, fig, **kwargs) -> None:
         """
         Add q-labels and unit name to figure ``fig``.
         Use units="" to add k-labels without unit name.
+
         Args:
             klabels:
             klabel_size:
@@ -2539,7 +2593,7 @@ def add_fundgap_span(self, ax_or_axlist, spin, span_dir="v", fontsize=8, **kwarg
 
         Args:
             ax_or_axlist:
-            spin
+            spin:
             spand_dir:
         """
         ks_lumo = self.lumos[spin]
@@ -2588,7 +2642,7 @@ def get_e0(self, e0):
 
     def plot_ax(self, ax, e0, spin=None, band=None, **kwargs) -> list:
         """
-        Helper function to plot the energies for (spin, band) on the axis ax with matplotlib..
+        Helper function to plot the energies for (spin, band) on the axis ax with matplotlib.
 
         Args:
             ax: |matplotlib-Axes|.
@@ -2720,7 +2774,6 @@ def plot_with_edos(self, edos, klabels=None, ax_list=None, e0="fermie", points=N
                 * ``edos_fermie``: Use the Fermi energy computed from the DOS to define the zero of energy in both subplots.
                 *  Number e.g ``e0 = 0.5``: shift all eigenvalues to have zero energy at 0.5 eV
                 *  None: Don't shift energies, equivalent to ``e0 = 0``
-
             points: Marker object with the position and the size of the marker.
                 Used for plotting purpose e.g. QP energies, energy derivatives...
             with_gaps: True to add markers and arrows showing the fundamental and the direct gap.
@@ -2799,7 +2852,6 @@ def plotly_with_edos(self, edos, klabels=None, fig=None, band_rcd=None, dos_rcd=
                 * ``edos_fermie``: Use the Fermi energy computed from the DOS to define the zero of energy in both subplots.
                 *  Number e.g ``e0 = 0.5``: shift all eigenvalues to have zero energy at 0.5 eV
                 *  None: Don't shift energies, equivalent to ``e0 = 0``
-
             points: Marker object with the position and the size of the marker.
                 Used for plotting purpose e.g. QP energies, energy derivatives...
             with_gaps: True to add markers and arrows showing the fundamental and the direct gap.
diff --git a/abipy/eph/gkq.py b/abipy/eph/gkq.py
index 3f3de7c65..423cc7874 100644
--- a/abipy/eph/gkq.py
+++ b/abipy/eph/gkq.py
@@ -198,7 +198,7 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
 
         # (nsppol, nkpt, 3*natom, mband, mband) real array.
         absg = np.abs(self.read_all_gkq(mode="phonon")) * abu.Ha_meV
-        absg = absg[:,ik]
+        absgk = absg[:,ik].copy()
         absg_unsym = absg[:,ik].copy()
 
         # Average over phonons.
@@ -208,11 +208,10 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
         natom3 = len(self.structure) * 3
         phfreqs_ha = self.phfreqs_ha
         nband = self.ebands.nband
-
-        absg_sym = np.empty_like(absg)
-        g2_mn = np.empty((nband,nband), dtype=float)
+        absg_sym = np.zeros_like(absgk)
 
         for spin in range(self.ebands.nsppol):
+            g2_mn = np.zeros((nband,nband), dtype=float)
             for nu in range(natom3):
                 w_1 = phfreqs_ha[nu]
                 g2_mn[:], nn = 0.0, 0
@@ -220,12 +219,12 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
                     w_2 = phfreqs_ha[mu]
                     if abs(w_1 - w_2) >= eps_ha: continue
                     nn += 1
-                    g2_mn += absg[spin,mu,:,:] ** 2
+                    g2_mn += absgk[spin,mu,:,:] ** 2
                 absg_sym[spin,nu,:,:] = np.sqrt(g2_mn / nn)
 
         # Average over k electrons.
-        absg = absg_sym
-        g2_nu = np.empty((natom3), dtype=float)
+        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):
@@ -235,23 +234,22 @@ def get_absg_kpoint(self, kpoint) -> tuple[np.ndarray, np.ndarray, int, Kpoint]:
                         w_2 = self.ebands.eigens[spin, ik, pbnd]
                         if abs(w_2 - w_1) >= eps_ev: continue
                         nn += 1
-                        g2_nu += absg[spin,:,jbnd,pbnd] ** 2
-                    absg_sym[spin,:,ibnd,jbnd] = np.sqrt(g2_nu / nn)
-
+                        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
-        absg = absg_sym
+        absgk = absg_sym.copy()
         for spin in range(self.ebands.nsppol):
-            for jbnd in range(nband):
-                for ibnd in range(nband):
-                    w_1 = eigens_kq[spin, ik, ibnd]
+            for ibnd in range(nband):
+                for jbnd in range(nband):
+                    w_1 = eigens_kq[spin, ik, jbnd]
                     g2_nu[:], nn = 0.0, 0
                     for pbnd in range(nband):
                         w_2 = eigens_kq[spin, ik, pbnd]
                         if abs(w_2 - w_1) >= eps_ev: continue
                         nn += 1
-                        g2_nu += absg[spin,:,jbnd,pbnd] ** 2
-                    absg_sym[spin,:,ibnd,jbnd] = np.sqrt(g2_nu / nn)
+                        g2_nu += absgk[spin,:,pbnd,ibnd] ** 2
+                    absg_sym[spin,:,jbnd,ibnd] = np.sqrt(g2_nu / nn)
 
         return absg_sym, absg_unsym, ik, kpoint
 
diff --git a/abipy/eph/varpeq.py b/abipy/eph/varpeq.py
index 8f1002522..80247b8f7 100644
--- a/abipy/eph/varpeq.py
+++ b/abipy/eph/varpeq.py
@@ -8,7 +8,6 @@
 
 import dataclasses
 import numpy as np
-#import pandas as pd
 import abipy.core.abinit_units as abu
 
 from collections import defaultdict
@@ -30,6 +29,7 @@
 from abipy.abio.robots import Robot
 from abipy.eph.common import BaseEphReader
 
+
 ITER_LABELS = [
     r'$E_{pol}$',
     r'$E_{el}$',
@@ -48,6 +48,7 @@ class VarpeqFile(AbinitNcFile, Has_Header, Has_Structure, Has_ElectronBands, Not
 
     .. code-block:: python
 
+        from abipy.eph.varpeq import VarpeqFile
         with VarpeqFile("out_VARPEQ.nc") as varpeq:
             print(varpeq)
             varpeq.plot_scf_cycle()
@@ -109,7 +110,6 @@ def to_string(self, verbose=0) -> str:
         app(self.filestat(as_string=True))
         app("")
         app(self.structure.to_string(verbose=verbose, title="Structure"))
-
         app("")
         app(self.ebands.to_string(with_structure=False, verbose=verbose, title="Electronic Bands"))
         #if verbose > 1:
@@ -211,8 +211,10 @@ def write_notebook(self, nbpath=None) -> str:
 class Polaron:
     """
     This object stores the polaron coefficients A_kn, B_qnu for a given spin.
-    Provides methods to plot |A_nk|^2 or |B_qnu|^2 together with band structures.
+    Provides methods to plot |A_nk|^2 or |B_qnu|^2 together with band structures
+    (fatbands-like plots)
     """
+
     spin: int          # Spin index.
     nb: int            # Number of bands.
     nk: int            # Number of k-points (including filtering if any)
@@ -254,7 +256,9 @@ def __str__(self) -> str:
         return self.to_string()
 
     def to_string(self, verbose=0) -> str:
-        """String representation with verbosiy level ``verbose``."""
+        """
+        String representation with verbosiy level verbose.
+        """
         lines = []; app = lines.append
 
         app(marquee(f"Ank for spin: {self.spin}", mark="="))
@@ -277,7 +281,7 @@ def to_string(self, verbose=0) -> str:
 
         return "\n".join(lines)
 
-    def get_a2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGridInterpolator:
+    def get_a2_interpolator(self, method: str, check_mesh: int=0) -> BzRegularGridInterpolator:
         """
         Build and return an interpolator for |A_nk|^2
 
@@ -285,7 +289,7 @@ def get_a2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGrid
             method: String defining the interpolation method.
             check_mesh: Check whether k-points belong to the mesh ff !=0.
         """
-        # Neeed to know the size of the k-mesh.
+        # Need to know the size of the k-mesh.
         ksampling = self.ebands.kpoints.ksampling
         ngkpt, shifts = ksampling.mpdivs, ksampling.shifts
 
@@ -303,10 +307,9 @@ def get_a2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGrid
                 ix, iy, iz = k_inds
                 a_data[ib, ix, iy, iz] = a_cplx
 
-        return BzRegularGridInterpolator(self.structure, shifts, np.abs(a_data) ** 2,
-                                         method=method)
+        return BzRegularGridInterpolator(self.structure, shifts, np.abs(a_data) ** 2, method=method)
 
-    def get_b2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGridInterpolator:
+    def get_b2_interpolator(self, method: str, check_mesh: int=0) -> BzRegularGridInterpolator:
         """
         Build and return an interpolator for |B_qnu|^2.
 
@@ -314,7 +317,7 @@ def get_b2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGrid
             method: String defining the interpolation method.
             check_mesh: Check whether k-points belong to the mesh ff !=0.
         """
-        # Neeed to know the size of the q-mesh (always Gamma-centered)
+        # Need to know the shape of the q-mesh (always Gamma-centered)
         ngqpt, shifts = self.varpeq.r.ngqpt, [0, 0, 0]
         q_indices = kpoints_indices(self.qpoints, ngqpt, check_mesh=check_mesh)
 
@@ -326,8 +329,7 @@ def get_b2_interpolator(self, method: str, check_mesh: int = 0) -> BzRegularGrid
                 ix, iy, iz = q_inds
                 b_data[nu, ix, iy, iz] = b_cplx
 
-        return BzRegularGridInterpolator(self.structure, shifts, np.abs(b_data) ** 2,
-                                         method=method)
+        return BzRegularGridInterpolator(self.structure, shifts, np.abs(b_data) ** 2, method=method)
 
     @add_fig_kwargs
     def plot_bz_sampling(self, what="kpoints", fold=False,
@@ -341,26 +343,31 @@ def plot_bz_sampling(self, what="kpoints", fold=False,
                 Defaults to False.
         """
         bz_points = dict(kpoints=self.kpoints, qpoints=self.qpoints)[what]
-        kws = dict(ax=ax, pmg_path=pmg_path, with_labels=with_labels, fold=fold,
-                   kpoints=bz_points)
+        kws = dict(ax=ax, pmg_path=pmg_path, with_labels=with_labels, fold=fold, kpoints=bz_points)
 
         return self.structure.plot_bz(show=False, **kws)
 
     @add_fig_kwargs
-    def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
-                             ax=None, scale=50, fontsize=12, **kwargs) -> Figure:
+    def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None, nksmall: int = 20,
+                             lpratio: int = 5, step: float = 0.1, width: float = 0.2, method: str = "linear",
+                             ax_list=None, ylims=None, scale=10, fontsize=8, **kwargs) -> Figure:
         """
         Plot electronic energies with markers whose size is proportional to |A_nk|^2.
 
         Args:
-            ebands: ElectronBands or Abipy file providing an electronic band structure along a path.
-            ax: |matplotlib-Axes| or None if a new figure should be created.
+            ebands_kpath: ElectronBands or Abipy file providing an electronic band structure along a path.
+            ebands_kmesh: ElectronBands or Abipy file providing an electronic band structure in the IBZ.
+            method=Interpolation method.
+            ax_list: List of |matplotlib-Axes| or None if a new figure should be created.
             scale: Scaling factor for |A_nk|^2.
+            ylims: Set the data limits for the y-axis. Accept tuple e.g. ``(left, right)``
+                   or scalar e.g. ``left``. If left (right) is None, default values are used.
+            fontsize: fontsize for legends and titles
         """
         ebands_kpath = ElectronBands.as_ebands(ebands_kpath)
 
         # Interpolate A_nk
-        a2_interp = self.get_a2_interpolator("linear")
+        a2_interp = self.get_a2_interpolator(method)
 
         # DEBUG SECTION
         #ref_kn = np.abs(self.a_kn) ** 2
@@ -382,7 +389,7 @@ def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
         points = Marker(x, y, s)
 
         nrows, ncols = 1, 2
-        ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
+        ax_list, fig, plt = get_axarray_fig_plt(ax_list, nrows=nrows, ncols=ncols,
                                                 sharex=False, sharey=True, squeeze=False)
         ax_list = ax_list.ravel()
 
@@ -390,34 +397,24 @@ def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
         ebands_kpath.plot(ax=ax, points=points, show=False)
 
         vertices_names = [(k.frac_coords, k.name) for k in ebands_kpath.kpoints]
-        lpratio = 5
-        kmesh = np.array([12, 12, 12]) * 3
-        step = 0.1
-        width = 0.2
 
         if ebands_kmesh is None:
-            # Compute ebands_kmesh with Star-function interpolation.
+            # Compute ebands_kmesh with star-function interpolation.
+            kmesh = self.structure.calc_ngkpt(nksmall)
             r = self.ebands.interpolate(lpratio=lpratio, vertices_names=vertices_names, kmesh=kmesh)
             ebands_kmesh = r.ebands_kmesh
 
         # Get electronic DOS.
         edos = ebands_kmesh.get_edos(step=step, width=width)
-
         mesh = edos.spin_dos[self.spin].mesh
         ank_dos = np.zeros(len(mesh))
         e0 = self.ebands.fermie
 
-        ymin -= 0.5 * abs(ymin)
-        ymin -= e0
-        ymax += 0.5 * abs(ymax)
-        ymax -= e0
-        #ymax = None
-
         #################
         # Compute Ank DOS
         #################
         # NB: This is just to sketch the ideas. I don't think the present version
-        # is correct as only k --> -k symmetry can be used.
+        # is correct as only the k --> -k symmetry can be used.
 
         for ik, kpoint in enumerate(ebands_kmesh.kpoints):
             weight = kpoint.weight
@@ -428,30 +425,64 @@ def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
 
         ax = ax_list[1]
         edos.plot_ax(ax, e0, spin=self.spin, exchange_xy=True, label="eDOS(E)")
+        ax.set_xlabel("arbitrary units", fontsize=fontsize)
 
         ank_dos = Function1D(mesh, ank_dos)
+        print("A2(E) integrates to:", ank_dos.integral_value, " Ideally, it should be 1.")
         ank_dos.plot_ax(ax, exchange_xy=True, label=r"$A^2$(E)")
         ax.grid(True)
         ax.legend(loc="best", shadow=True, fontsize=fontsize)
-        print("A2(E) integrates to", ank_dos.integral_value)
+
+        if ylims is None:
+            ymin -= 0.5 * abs(ymin)
+            ymin -= e0
+            ymax += 0.5 * abs(ymax)
+            ymax -= e0
+            ylims = [ymin, ymax]
 
         for ax in ax_list:
-            ax.set_ylim(ymin, ymax)
+            set_axlims(ax, ylims, "y")
 
         return fig
 
     @add_fig_kwargs
-    def plot_bqnu_with_phbands(self, phbands_qpath, ax=None, scale=10, **kwargs) -> Figure:
+    def plot_bqnu_with_ddb(self, ddb, with_phdos=True, anaget_kwargs=None, **kwargs) -> Figure:
+        """
+        High-level interface to plot phonon energies with markers whose size is proportional to |B_qnu|^2.
+        Similar to plot_bqnu_with_phbands but this function receives in input a DdbFile or a
+        path to a ddbfile and automates the computation of the phonon bands by invoking anaddb.
+
+        Args:
+            ddb:
+            anaget_kwargs:
+        """
+        from abipy.dfpt.ddb import DdbFile
+        ddb = DdbFile.as_ddb(ddb)
+        anaget_kwargs = {} if anaget_kwargs is None else anaget_kwargs
+
+        with ddb.anaget_phbst_and_phdos_files(**anaget_kwargs) as g:
+            phbst_file, phdos_file = g[0], g[1]
+            phbands_qpath = phbst_file.phbands
+            return self.plot_bqnu_with_phbands(phbands_qpath,
+                                               phdos_file=phdos_file if with_phdos else None,
+                                               **kwargs)
+
+    @add_fig_kwargs
+    def plot_bqnu_with_phbands(self, phbands_qpath, phdos_file=None,
+                               method="linear",
+                               ax=None, scale=10, **kwargs) -> Figure:
         """
         Plot phonon energies with markers whose size is proportional to |B_qnu|^2.
 
         Args:
             phbands_qpath: PhononBands or Abipy file providing a phonon band structure.
+            phdos_file:
+            method=Interpolation method.
             ax: |matplotlib-Axes| or None if a new figure should be created.
             scale: Scaling factor for |B_qnu|^2.
         """
         phbands_qpath = PhononBands.as_phbands(phbands_qpath)
-        b2_interp = self.get_b2_interpolator("linear")
+        b2_interp = self.get_b2_interpolator(method)
 
         # DEBUG SECTION
         #ref_qnu = np.abs(self.b_qnu) ** 2
@@ -464,13 +495,42 @@ def plot_bqnu_with_phbands(self, phbands_qpath, ax=None, scale=10, **kwargs) ->
         for iq, qpoint in enumerate(phbands_qpath.qpoints):
             omegas_nu = phbands_qpath.phfreqs[iq,:]
             b2_nu = b2_interp.eval_kpoint(qpoint.frac_coords)
-            assert len(omegas_nu) == len(b2_nu)
             for w, b2 in zip(omegas_nu, b2_nu):
                 x.append(iq); y.append(w); s.append(scale * b2)
+
         points = Marker(x, y, s)
 
-        ax, fig, plt = get_ax_fig_plt(ax=ax)
-        phbands_qpath.plot(ax=ax, points=points, show=False)
+        nrows, ncols = 1, 1 if phdos_file is None else 2
+        ax_list, fig, plt = get_axarray_fig_plt(None, nrows=nrows, ncols=ncols,
+                                                sharex=False, sharey=True, squeeze=False)
+        ax_list = ax_list.ravel()
+
+        phbands_qpath.plot(ax=ax_list[0], points=points, show=False)
+
+        if phdos_file is None:
+            return fig
+
+        # Add phdos and |B_qn| dos
+        # mesh is given in eV, values are in states/eV.
+        #phbands_qmesh = phdos_file.phbands
+        phdos = phdos_file.phdos
+        ngqpt = np.diagonal(phdos_file.qptrlatt)
+        mesh = phdos.mesh
+        bqnu_dos = np.zeros(len(mesh))
+
+        from abipy.core.kpoints import Ktables
+        qtabs = Ktables(self.structure, ngqpt, is_shift=False, has_timrev=True)
+        print(qtabs)
+
+        width = 0.2
+        for iq, qpoint in enumerate(phbands_qmesh.qpoints):
+            weight = qpoint.weight
+            #enes_n = phbands_qmesh.eigens[self.spin, iq, self.bstart:self.bstop]
+            b2_nu = b2_interp.eval_kpoint(qpoint.frac_coords)
+            for w, b2 in zip(enes_n, b2_nu):
+                bqnu_dos += weight * b2 * gaussian(mesh, width, center=w)
+
+        #ank_dos = Function1D(mesh, ank_dos)
 
         return fig
 
diff --git a/abipy/flowtk/psrepos.py b/abipy/flowtk/psrepos.py
index a0839ff80..268a6d614 100644
--- a/abipy/flowtk/psrepos.py
+++ b/abipy/flowtk/psrepos.py
@@ -429,10 +429,10 @@ class JthRepo(PseudosRepo):
 
     @classmethod
     def from_abinit_website(cls, xc_name: str, relativity_type: str, version: str) -> JthRepo:
-
         # https://www.abinit.org/ATOMICDATA/JTH-LDA-atomicdata.tar.gz
         # ATOMPAW-LDA-JTHv0.4
-        url = f"https://www.abinit.org/ATOMICDATA/JTH-{xc_name}-atomicdata.tar.gz"
+        #url = f"https://www.abinit.org/ATOMICDATA/JTH-{xc_name}-atomicdata.tar.gz"
+        url = f"https://abinit.github.io/abinit_web/ATOMICDATA/JTH-{xc_name}-atomicdata.tar.gz"  # TODO New url
         return cls(cls.ps_generator, xc_name, relativity_type, cls.project_name, version, url)
 
     @property
diff --git a/abipy/flowtk/tasks.py b/abipy/flowtk/tasks.py
index d1d5f16d3..025e09da0 100644
--- a/abipy/flowtk/tasks.py
+++ b/abipy/flowtk/tasks.py
@@ -631,8 +631,17 @@ def from_user_config(cls) -> TaskManager:
         #    return _USER_CONFIG_TASKMANAGER
 
         # Try in the current directory then in user configuration directory.
-        path = os.path.join(os.getcwd(), cls.YAML_FILE)
-        if not os.path.exists(path):
+        try:
+            cwd = os.getcwd()
+        except FileNotFoundError:
+            # This error is triggered when running pytes with workers.
+            cwd = None
+
+        if cwd is not None:
+            path = os.path.join(cwd, cls.YAML_FILE)
+            if not os.path.exists(path):
+                path = os.path.join(cls.USER_CONFIG_DIR, cls.YAML_FILE)
+        else:
             path = os.path.join(cls.USER_CONFIG_DIR, cls.YAML_FILE)
 
         if not os.path.exists(path):
@@ -1263,7 +1272,7 @@ def version_ge(self, version_string):
 
     def compare_version(self, version_string, op):
         """Compare Abinit version to `version_string` with operator `op`"""
-        from pkg_resources import parse_version
+        from packaging.version import parse as parse_version
         from monty.operator import operator_from_str
         op = operator_from_str(op)
         return op(parse_version(self.version), parse_version(version_string))
diff --git a/abipy/ml/aseml.py b/abipy/ml/aseml.py
index 823bdc154..d5d57e570 100644
--- a/abipy/ml/aseml.py
+++ b/abipy/ml/aseml.py
@@ -46,7 +46,7 @@
 from abipy.tools.iotools import workdir_with_prefix, PythonScript, yaml_safe_load_path
 from abipy.tools.typing import Figure, PathLike
 from abipy.tools.printing import print_dataframe
-from abipy.tools.serialization import HasPickleIO
+from abipy.tools.serialization import HasPickleIO, mjson_write
 from abipy.tools.context_managers import Timer
 from abipy.tools.parallel import get_max_nprocs # , pool_nprocs_pmode
 from abipy.abio.enums import StrEnum, EnumMixin
@@ -3195,13 +3195,16 @@ def __init__(self, atoms, nn_name, verbose, workdir, prefix=None):
         self.verbose = verbose
 
     def run(self):
+        """Run the calculation."""
         calc = CalcBuilder(self.nn_name).get_calculator()
         self.atoms.calc = calc
         res = AseResults.from_atoms(self.atoms)
         print(res.to_string(verbose=self.verbose))
 
-        # Write json file with GS results.
-        from abipy.tools.serialization import mjson_write
+        # Write json file GS results.
+        # To read the dictionary from json use:
+        #   from abipy.tools.serialization import mjson_load
+        #   data = mjson_load(self.workdir / "gs.json")
         data = dict(
             structure=Structure.as_structure(self.atoms),
             ene=res.ene,
@@ -3210,10 +3213,6 @@ def run(self):
         )
         mjson_write(data, self.workdir / "gs.json", indent=4)
 
-        # To read the dictionary from json use:
-        #from abipy.tools.serialization import mjson_load
-        #data = mjson_load(self.workdir / "gs.json")
-
         # Write ASE trajectory file with results.
         with open(self.workdir / "gs.traj", "wb") as fd:
             write_traj(fd, [self.atoms])
@@ -3221,6 +3220,82 @@ def run(self):
         return 0
 
 
+class FrozenPhononMl(MlBase):
+    """
+    Frozen-phonon calculations with ML potential.
+    """
+
+    @classmethod
+    def from_ddb_file(cls, ddb_filepath, qpoint, eta_list, nn_name, verbose, workdir, prefix=None, **anaddb_kwargs):
+        """
+        """
+        from abipy.dfpt.ddb import DdbFile
+        with DdbFile(ddb_filepath) as ddb:
+            # Call anaddb to get all phonon modes for this q-point.
+            phbands = ddb.anaget_phmodes_at_qpoint(qpoint=qpoint, verbose=verbose, **anaddb_kwargs)
+
+        return cls(ddb.structure, qpoint, phbands.phdispl_cart, eta_list, nn_name, verbose, workdir, prefix=prefix)
+
+    def __init__(self, structure, qpoint, phdispl_cart, eta_list, nn_name, verbose, workdir, prefix=None):
+        """
+        Args:
+            qpoint: q-vector in reduced coordinate in reciprocal space.
+            displ_cart: displacement of the atoms in real space .
+            eta: pre-factor multiplying the displacement. Gives the value in Angstrom of the largest displacement.
+        """
+        super().__init__(workdir, prefix)
+        self.initial_structure = structure
+        natom = len(structure)
+        self.nn_name = nn_name
+        self.verbose = verbose
+        # TODO: Should check that qpoint is [1/Nx, 1/Ny, 1/Nz]
+        self.qpoint = np.array(qpoint)
+        self.phdispl_cart = phdispl_cart
+        #self.phdispl_cart = np.reshape(phdispl_cart, (-1, 3*natom, 3*natom))
+        self.eta_list = np.array(eta_list)
+
+    def run(self):
+        """Run the calculation."""
+        calc = CalcBuilder(self.nn_name).get_calculator()
+
+        max_sc = np.ones(3, dtype=int)
+        for i, qf in enumerate(self.qpoint):
+            if qf != 0: max_sc[i] = np.round(1 / qf)
+        print(f"{max_sc =}")
+
+        for displ_cart in self.phdispl_cart:
+            for eta in self.eta_list:
+                print(f"{eta=}")
+                print(f"{displ_cart.shape=}")
+                print(f"{displ_cart=}")
+                scell = self.initial_structure.frozen_phonon(self.qpoint, displ_cart, eta=eta, frac_coords=True, max_supercell=max_sc)
+                print(scell.scale_matrix)
+                print(scell.structure)
+
+                scell.structure.to_ase_atoms()
+                atoms.calc = calc
+                res = AseResults.from_atoms(atoms)
+                #print(res.to_string(verbose=self.verbose))
+
+        # Write json file GS results.
+        # To read the dictionary from json use:
+        #   from abipy.tools.serialization import mjson_load
+        #   data = mjson_load(self.workdir / "gs.json")
+        #data = dict(
+        #    structure=Structure.as_structure(self.atoms),
+        #    ene=res.ene,
+        #    stress=res.stress,
+        #    forces=res.forces,
+        #)
+        #mjson_write(data, self.workdir / "gs.json", indent=4)
+
+        ## Write ASE trajectory file with results.
+        #with open(self.workdir / "gs.traj", "wb") as fd:
+        #    write_traj(fd, [self.atoms])
+
+        return 0
+
+
 class MlCompareNNs(MlBase):
     """
     Compare energies, forces and stresses obtained with different ML potentials.
diff --git a/abipy/scripts/abicomp.py b/abipy/scripts/abicomp.py
index 3e2b9010e..8025c4156 100755
--- a/abipy/scripts/abicomp.py
+++ b/abipy/scripts/abicomp.py
@@ -16,7 +16,7 @@
 from monty.functools import prof_main
 from monty.termcolor import cprint
 from abipy import abilab
-from abipy.tools.plotting import get_ax_fig_plt, GenericDataFilesPlotter, FilesPlotter
+from abipy.tools.plotting import get_ax_fig_plt, GenericDataFilesPlotter, FilesPlotter, Exposer
 
 
 def remove_disordered(structures, paths):
@@ -594,10 +594,10 @@ def abicomp_abiwan(options):
 
 def abicomp_abiwan_ebands(options):
     """
-    Compare Wannier-interpolated band structure with ab-initio results.
+    Compare Wannier-interpolated band structure with ab-initio data.
     """
     if len(options.paths) != 2:
-        raise ValueError("Two arguments with ABIWAN.nc and nc file with ElectronBands are required.")
+        raise ValueError("Two arguments with ABIWAN.nc and netcdf file with ElectronBands are required!")
     from abipy.wannier90 import AbiwanFile
     abiwan_path, ebands_path = options.paths[0], options.paths[1]
     if not abiwan_path.endswith("ABIWAN.nc"):
@@ -605,8 +605,11 @@ def abicomp_abiwan_ebands(options):
 
     abiwan = AbiwanFile(abiwan_path)
     print(abiwan)
-    abiwan.hwan.plot()
-    abiwan.plot_with_ebands(ebands_path)
+
+    with Exposer.as_exposer("mpl") as e:
+        plot_kws = dict(show=False)
+        e(abiwan.hwan.plot(**plot_kws))
+        e(abiwan.plot_with_ebands(ebands_path, **plot_kws))
 
     return 0
 
diff --git a/abipy/scripts/abiml.py b/abipy/scripts/abiml.py
index 90e1cfb5f..2b656bbb5 100755
--- a/abipy/scripts/abiml.py
+++ b/abipy/scripts/abiml.py
@@ -17,6 +17,7 @@
 from abipy.core.structure import Structure
 from abipy.tools.printing import print_dataframe
 
+
 ASE_OPTIMIZERS = aseml.ase_optimizer_cls("__all__")
 
 DEFAULT_NN = "chgnet"
@@ -721,6 +722,28 @@ def gs(ctx, filepath, nn_name,
     return 0
 
 
+@main.command()
+@herald
+@click.pass_context
+@click.argument("filepath", type=str)
+@click.option("--qpoint", "-q", nargs=3, type=float, help="q-point in reduced coordinates.")
+@add_nn_name_opt
+@add_workdir_verbose_opts
+@click.option('--config', default='abiml_phfrozen.yml', type=click.Path(), callback=set_default, is_eager=True, expose_value=False)
+def phddb_frozen(ctx, filepath, qpoint, nn_name,
+       workdir, verbose,
+       ):
+    """
+    Frozen-phonon calculation with ML potential.
+    """
+    qpoint = [0, 0, 0]
+    eta_list = [1, 2]
+    frozen_ph = aseml.FrozenPhononMl.from_ddb_file(filepath, qpoint, eta_list, nn_name, verbose, workdir, prefix="_abiml_phfrozen")
+    frozen_ph.run()
+
+    return 0
+
+
 @main.command()
 @herald
 @click.pass_context
diff --git a/abipy/scripts/abistruct.py b/abipy/scripts/abistruct.py
index 23b095405..ca42cb154 100755
--- a/abipy/scripts/abistruct.py
+++ b/abipy/scripts/abistruct.py
@@ -575,6 +575,10 @@ def show_examples_and_exit(err_msg=None, error_code=1):
     if options.verbose > 2:
         print(options)
 
+    #structure = abilab.Structure.from_file(options.filepath)
+    #print("structure.has_zero_dynamical_quadrupoles:", structure.has_zero_dynamical_quadrupoles)
+    #sys.exit(0)
+
     if options.command == "spglib":
         structure = abilab.Structure.from_file(options.filepath)
         print(structure.spget_summary(symprec=options.symprec, angle_tolerance=options.angle_tolerance,
diff --git a/abipy/tools/plotting.py b/abipy/tools/plotting.py
index 69af0fe9d..d259b6b13 100644
--- a/abipy/tools/plotting.py
+++ b/abipy/tools/plotting.py
@@ -14,6 +14,7 @@
 import functools
 import numpy as np
 import pandas as pd
+import matplotlib.collections as mcoll
 
 from collections import namedtuple, OrderedDict
 from typing import Any, Callable, Iterator
@@ -23,8 +24,6 @@
 from abipy.tools.iotools import dataframe_from_filepath
 from abipy.tools.typing import Figure, Axes, VectorLike
 from abipy.tools.numtools import data_from_cplx_mode
-
-import matplotlib.collections as mcoll
 from plotly.tools import mpl_to_plotly
 
 __all__ = [
@@ -80,7 +79,9 @@ def __init__(self, filepaths: list[str]):
 
     @add_fig_kwargs
     def plot(self, **kwargs) -> Figure:
-        """Loop through the PNG files and display them in subplots."""
+        """
+        Loop through the PNG files and display them in subplots.
+        """
         # Build grid of plots.
         num_plots, ncols, nrows = len(self.filepaths), 1, 1
         if num_plots > 1:
@@ -942,9 +943,8 @@ def plot(self, cplx_mode="abs", colormap="jet", fontsize=8, **kwargs) -> Figure:
         return fig
 
 
-#TODO use object and introduce c for color, client code should be able to customize it.
-# Rename it to ScatterData
-class Marker(namedtuple("Marker", "x y s")):
+#TODO Rename it to ScatterData?
+class Marker:
     """
     Stores the position and the size of the marker.
     A marker is a list of tuple(x, y, s) where x, and y are the position
@@ -955,49 +955,25 @@ class Marker(namedtuple("Marker", "x y s")):
 
         x, y, s = [1, 2, 3], [4, 5, 6], [0.1, 0.2, -0.3]
         marker = Marker(x, y, s)
-        marker.extend((x, y, s))
-
     """
-    def __new__(cls, *xys):
-        """Extends the base class adding consistency check."""
-        if not xys:
-            xys = ([], [], [])
-            return super().__new__(cls, *xys)
-
-        if len(xys) != 3:
-            raise TypeError("Expecting 3 entries in xys got %d" % len(xys))
 
-        x = np.asarray(xys[0])
-        y = np.asarray(xys[1])
-        s = np.asarray(xys[2])
-        xys = (x, y, s)
+    def __init__(self, x, y, s, marker: str = "o", color: str = "y", alpha: float = 1.0):
+        self.x, self.y, self.s = np.array(x), np.array(y), np.array(s)
 
-        for s in xys[-1]:
-            if np.iscomplex(s):
-                raise ValueError("Found ambiguous complex entry %s" % str(s))
+        if len(self.x) != len(self.y):
+            raise ValueError("len(self.x) != len(self.y)")
+        if len(self.y) != len(self.s):
+            raise ValueError("len(self.y) != len(self.s)")
 
-        return super().__new__(cls, *xys)
+        self.marker = marker
+        self.color = color
+        self.alpha = alpha
 
     def __bool__(self):
         return bool(len(self.s))
 
     __nonzero__ = __bool__
 
-    def extend(self, xys):
-        """
-        Extend the marker values.
-        """
-        if len(xys) != 3:
-            raise TypeError("Expecting 3 entries in xys got %d" % len(xys))
-
-        self.x.extend(xys[0])
-        self.y.extend(xys[1])
-        self.s.extend(xys[2])
-
-        lens = np.array((len(self.x), len(self.y), len(self.s)))
-        if np.any(lens != lens[0]):
-            raise TypeError("x, y, s vectors should have same lengths but got %s" % str(lens))
-
     def posneg_marker(self) -> tuple[Marker, Marker]:
         """
         Split data into two sets: the first one contains all the points with positive size.
@@ -2492,10 +2468,8 @@ def plotly_points(points, lattice=None, coords_are_cartesian=False, fold=False,
     from pymatgen.electronic_structure.plotter import fold_point
     vecs = []
     for p in points:
-
         if fold:
             p = fold_point(p, lattice, coords_are_cartesian=coords_are_cartesian)
-
         elif not coords_are_cartesian:
             p = lattice.get_cartesian_coords(p)
 
@@ -2666,8 +2640,11 @@ def add_colorscale_dropwdowns(fig):
 
     return fig
 
-def mpl_to_ply(fig, latex=False):
-    """Nasty workaround for plotly latex rendering in legend/breaking exception"""
+
+def mpl_to_ply(fig: Figure, latex: bool= False):
+    """
+    Nasty workaround for plotly latex rendering in legend/breaking exception
+    """
     if is_plotly_figure(fig):
         return fig
 
@@ -2748,3 +2725,46 @@ def parse_latex(label):
         trace.name = new_label
 
     return plotly_fig
+
+
+class PolyfitPlotter:
+    """
+    Fit data with polynomals of different degrees and visualize the results.
+    """
+
+    def __init__(self, xs, ys):
+        self.xs, self.ys = xs, ys
+
+    @add_fig_kwargs
+    def plot(self, deg_list: list[int], num=100, ax=None, xlabel=None, ylabel=None,
+             fontsize=8, **kwargs) -> Figure:
+        """
+        Args:
+            deg_list: List with degrees of the fitting polynomial.
+            num: Number of samples to generate. Default is 100. Must be non-negative.
+            ax: |matplotlib-Axes| or None if a new figure should be created.
+            fontsize: Legend fontsize.
+        """
+        xs, ys = self.xs, self.ys
+        ax, fig, plt = get_ax_fig_plt(ax=ax)
+
+        for i, deg in enumerate(deg_list):
+            # Fit a ndeg polynomial to the data points and get the polynomial function.
+            coefficients = np.polyfit(xs, ys, deg)
+            polynomial = np.poly1d(coefficients)
+            #print("Coefficients:", coefficients); print("Polynomial:", polynomial)
+
+            if i == 0:
+                # Plot the original data points
+                ax.scatter(xs, ys, color='red', marker="o", label='Data Points')
+
+            # Generate (x, y) values for plotting the fit
+            x_fit = np.linspace(min(xs), max(xs), num)
+            y_fit = polynomial(x_fit)
+            ax.plot(x_fit, y_fit, label=f"{deg}-order fit")
+
+        if xlabel is not None: ax.set_xlabel(xlabel)
+        if ylabel is not None: ax.set_ylabel(ylabel)
+        ax.legend(loc="best", fontsize=fontsize, shadow=True)
+
+        return fig
diff --git a/abipy/tools/tests/test_plotting.py b/abipy/tools/tests/test_plotting.py
index 5582f8de6..236cdd037 100644
--- a/abipy/tools/tests/test_plotting.py
+++ b/abipy/tools/tests/test_plotting.py
@@ -102,11 +102,10 @@ def test_array_plotter(self):
 
     def test_marker(self):
         """Testing Marker."""
-        marker = Marker()
-        assert not marker
-        assert not marker.x
+
         x, y, s = [1, 2, 3], [4, 5, 6], [0.1, 0.2, -0.3]
-        marker.extend((x, y, s))
+        marker = Marker(x, y, s)
+
         assert marker
         self.assert_equal(marker.x, x)
         self.assert_equal(marker.y, y)
@@ -117,13 +116,6 @@ def test_marker(self):
         assert all(s < 0 for s in neg_mark.s)
         assert len(neg_mark.s) == 1 and neg_mark.s[0] == -0.3
 
-        with self.assertRaises(TypeError):
-            marker.extend((x, y))
-
-        with self.assertRaises(TypeError):
-            x.append(-1)
-            marker.extend((x, y, s))
-
     def test_plot_cell_tools(self):
         """Testing plot_unit_cell."""
         lattice = abilab.Lattice.hexagonal(a=2, c=4)
diff --git a/requirements.txt b/requirements.txt
index 9e7217df1..0a8e05f69 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
+packaging
 monty
 tabulate
 apscheduler<=3.10.4
diff --git a/setup.py b/setup.py
index 2a2899edb..5a41561d4 100755
--- a/setup.py
+++ b/setup.py
@@ -145,6 +145,7 @@ def cleanup():
 
 install_requires = [
     "monty",
+    "packaging",
     "tabulate",
     "apscheduler<=3.10.4",
     "pydispatcher>=2.0.5",