diff --git a/abipy/core/structure.py b/abipy/core/structure.py
index b72e2481c..63e449ea5 100644
--- a/abipy/core/structure.py
+++ b/abipy/core/structure.py
@@ -1328,7 +1328,6 @@ def check_image(structure, site):
 
         return True
 
-
     @lazy_property
     def hsym_kpath(self):
         """
@@ -2410,6 +2409,17 @@ def calc_ngkpt(self, nksmall) -> np.ndarray:
 
         return ngkpt
 
+    #def any_to_ngkpt(self, any_var) -> np.ndarray:
+    #    any_var = np.array(any_var)
+    #    if len(any_var) == 3:
+    #        return any_var
+    #    if (nksmall := float(any_var)) > 0:
+    #        return self.calc_ngkpt(nksmall)
+
+    #    import pymatgen.io.abinit.abiobjects as aobj
+    #    ksampling = aobj.KSampling.automatic_density(self, kppa, chksymbreak=0, shifts=(0,0,0))
+    #    return ksampling.to_abivars()["ngkpt"]
+
     def calc_shiftk(self, symprec=0.01, angle_tolerance=5) -> np.ndarray:
         """
         Find the values of ``shiftk`` and ``nshiftk`` appropriated for the sampling of the Brillouin zone.
diff --git a/abipy/dfpt/qha.py b/abipy/dfpt/qha.py
index 6e0474eb0..bb513c68b 100644
--- a/abipy/dfpt/qha.py
+++ b/abipy/dfpt/qha.py
@@ -100,7 +100,7 @@ def fit_energies(self, tstart=0, tstop=800, num=100):
         # list of minimum volumes and energies, one for each temperature
         min_volumes = np.array([fit.v0 for fit in fits])
         min_energies = np.array([fit.e0 for fit in fits])
-        F2D=  np.array([fit.b0 for fit in fits]) /min_volumes 
+        F2D=  np.array([fit.b0 for fit in fits]) /min_volumes
 
         return dict2namedtuple(tot_en=tot_en, fits=fits, min_en=min_energies, min_vol=min_volumes, temp=tmesh , F2D=F2D)
 
@@ -209,10 +209,10 @@ def get_thermal_expansion_coeff(self, tstart=0, tstop=800, num=100, tref=None, m
         eos_list = [ 'taylor', 'murnaghan', 'birch', 'birch_murnaghan','pourier_tarantola', 'vinet', 'antonschmidt']
 
         dt = f.temp[1] - f.temp[0]
-        if (method=="finite_difference"): 
+        if (method=="finite_difference"):
             alpha = (f.min_vol[2:] - f.min_vol[:-2]) / (2 * dt) / f.min_vol[1:-1]
         else :
-            if (self.eos_name in eos_list) : 
+            if (self.eos_name in eos_list) :
                 thermo = self.get_thermodynamic_properties(tstart=tstart, tstop=tstop, num=num)
                 entropy = thermo.entropy.T #* abu.e_Cb * abu.Avogadro
                 df_t = np.zeros((num,self.nvols))
@@ -227,10 +227,10 @@ def get_thermal_expansion_coeff(self, tstart=0, tstop=800, num=100, tref=None, m
                     param2[j] = np.array([3*param[j][0],2*param[j][1],param[j][2]])
 
                     p = np.poly1d(param2[j])
-                    d2f_t_v[j]= p(f.min_vol[j]) 
+                    d2f_t_v[j]= p(f.min_vol[j])
 
                 if tref is None:
-                    alpha= - 1/f.min_vol[1:-1] *d2f_t_v[1:-1] / f.F2D[1:-1] 
+                    alpha= - 1/f.min_vol[1:-1] *d2f_t_v[1:-1] / f.F2D[1:-1]
                 else :
                     alpha= - 1/f0.min_vol * d2f_t_v[1:-1] / f.F2D[1:-1]
             else :
@@ -789,7 +789,6 @@ def get_thermodynamic_properties(self, tstart=0, tstop=800, num=100):
         cv = np.zeros((self.nvols, num))
         free_energy = np.zeros((self.nvols, num))
         entropy = np.zeros((self.nvols, num))
-        internal_energy = np.zeros((self.nvols, num))
         zpe = np.zeros(self.nvols)
 
         for i, d in enumerate(self.doses):
diff --git a/abipy/dfpt/qha_aproximation.py b/abipy/dfpt/qha_aproximation.py
index 692670fc0..cf482f709 100644
--- a/abipy/dfpt/qha_aproximation.py
+++ b/abipy/dfpt/qha_aproximation.py
@@ -1,5 +1,8 @@
 # coding: utf-8
-"""Classes for V-ZSISA approximation."""
+"""
+Classes for V-ZSISA approximation.
+See PHYSICAL REVIEW B 110, 014103 (2024)
+"""
 from __future__ import annotations
 
 import numpy as np
@@ -10,10 +13,10 @@
 from pymatgen.analysis.eos import EOS
 from abipy.core.func1d import Function1D
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, set_grid_legend #, get_axarray_fig_plt
-from abipy.tools.typing import Figure
+from abipy.tools.typing import Figure, PathLike
 from abipy.electrons.gsr import GsrFile
 from abipy.dfpt.ddb import DdbFile
-from abipy.dfpt.phonons import PhdosFile
+from abipy.dfpt.phonons import PhdosFile, PhononDosPlotter
 
 
 class QHA_App:
@@ -24,8 +27,51 @@ class QHA_App:
     Does not include electronic entropic contributions for metals.
     """
 
+    def from_gsr_ddb_files(cls,
+                           gsr_paths: list[PathLike],
+                           ddb_paths: list[PathLike],
+                           verbose: int = 0,
+                           **extra_kwargs) -> QHA_App:
+        """
+        Creates an instance from a list of GSR files and a list of DDB files.
+        This is a simplified interface that computes the PHDOS.nc files automatically
+        from the DDB files by invoking anaddb
+
+        Args:
+            gsr_paths: list of paths to GSR files.
+            ddb_paths: list of paths to DDB files.
+        """
+        # From PHYSICAL REVIEW B110,014103(2024)
+        # The phonon density of states (PHDOS)was determined utilizing the Gaussian method,
+        # with a DOS smearing value set to 4.5×10−6 Hartree.
+        # Furthermore, a frequency grid step of 1.0×10−6 Hartree was employed for PHDOS calculations.
+        # These adjustments in numerical accuracy were imperative for versions of ABINIT preceding v9.10.
+        # Notably, in ABINIT versions v9.10 and after, these parameter values are preset as defaults for calculations.
+        phdos_paths = []
+        smear = "4.5 e−6 Hartree"
+        kwargs = dict(#nqsmall,
+                      qppa=100,
+                      dos_method=f"gaussian:{smear}",
+                      ndivsm=0,
+                      return_input=True,
+                     )
+
+        kwargs.update(extra_kwargs)
+        if verbose:
+            print(f"Calling anaget_phbst_and_phdos_files with {kwargs=}:")
+
+        for ddb_path in ddb_paths:
+            with DdbFile(ddb_path) as ddb:
+                with ddb.anaget_phbst_and_phdos_files(**kwargs) as g:
+                    if verbose:
+                        print(f"anaddb input file: {g.input=}")
+                    phbst_file, phdos_file = g[0], g[1]
+                    phdos_paths.append(phdos_file.filepath)
+
+        return cls.from_gsr_phdos_files(gsr_paths, phdos_paths)
+
     @classmethod
-    def from_files_app(cls, gsr_paths: list[str], phdos_paths: list[str]) -> QHA_App:
+    def from_gsr_phdos_files(cls, gsr_paths: list[str], phdos_paths: list[str]) -> QHA_App:
         """
         Creates an instance from a list of GSR files and a list of PHDOS.nc files.
 
@@ -66,8 +112,11 @@ def from_files_app(cls, gsr_paths: list[str], phdos_paths: list[str]) -> QHA_App
 
         return cls(structures, structures_from_phdos, index_list, phdoses, energies, pressures)
 
+    # FIXME: This to maintain previous interface. It should be removed
+    from_files_app = from_gsr_phdos_files
+
     @classmethod
-    def from_files_app_ddb(cls, ddb_paths: list[str], phdos_paths: list[str]) -> QHA_App:
+    def from_ddb_phdos_files(cls, ddb_paths: list[str], phdos_paths: list[str]) -> QHA_App:
         """
         Creates an instance from a list of DDB files and a list of PHDOS.nc files.
 
@@ -109,15 +158,18 @@ def from_files_app_ddb(cls, ddb_paths: list[str], phdos_paths: list[str]) -> QHA
 
         return cls(structures, structures_from_phdos, index_list, phdoses, energies, pressures)
 
+    # FIXME: This to maintain previous interface. It should be removed
+    from_files_app_ddb = from_ddb_phdos_files
+
     def __init__(self, structures, structures_from_phdos, index_list, phdoses, energies, pressures,
-                 eos_name: str='vinet', pressure: float=0.0):
+                 eos_name: str = 'vinet', pressure: float = 0.0):
         """
         Args:
-            structures: list of structures at different volumes.
-            structure_from_phdos:
-            index_list:
+            structures: list of structures at the different volumes for the BO energies.
+            structure_from_phdos: Structures used to compute phonon DOS
+            index_list: Index of each phdos in structures.
             phdoses: Phonon DOSes.
-            energies: list of SCF energies for the structures in eV.
+            energies: list of BO energies for the structures in eV.
             pressures: value of the pressure in GPa that will be considered in the p*V contribution to the energy.
             eos_name: string indicating the expression used to fit the energies. See pymatgen.analysis.eos.EOS.
             pressure:
@@ -137,7 +189,7 @@ def __init__(self, structures, structures_from_phdos, index_list, phdoses, energ
 
         self.angles_alpha = np.array([s.lattice.angles[0] for s in structures])
         self.angles_beta = np.array([s.lattice.angles[1] for s in structures])
-        self.angles_gama = np.array([s.lattice.angles[2] for s in structures])
+        self.angles_gamma = np.array([s.lattice.angles[2] for s in structures])
 
         self.phdoses = phdoses
         self.structures_from_phdos = np.array(structures_from_phdos)
@@ -168,6 +220,22 @@ def nvols(self) -> int:
         """Number of volumes"""
         return len(self.volumes_from_phdos)
 
+    def to_string(self, verbose: int = 0) -> str:
+        """String representation with verbosity level verbose:"""
+        lines = []
+        app = lines.append
+        #app(self.structure.to_string(verbose=verbose))
+        app(f"BO volumes: {self.volumes} Ang^3")
+        app(f"PHDOS volumes: {self.volumes_from_phdos} Ang^3")
+        app(f"eos_name: {self.eos_name}")
+        app(f"pressure: {self.pressure} GPa")
+        app(f"scale_points: {self.scale_points} GPa")
+
+        return "\n".join(lines)
+
+    def __str__(self) -> str:
+        return self.to_string()
+
     def set_eos(self, eos_name: str) -> None:
         """
         Set the EOS model used for the fit.
@@ -426,7 +494,7 @@ def vol_Einf_Vib4(self, tstart=0, tstop=1000, num=101) -> tuple:
         V0 = volumes0[iv0]
 
         tot_en = ( (volumes[np.newaxis, :].T - V0) * dfe_dV1 + 0.5 * ( volumes[np.newaxis, :].T - V0)**2*(dfe_dV2)
-                 + (volumes[np.newaxis, :].T - V0)**3 * dfe_dV3/6.0  + ( volumes[np.newaxis, :].T - V0)**4*(dfe_dV4/24.0)
+                 + (volumes[np.newaxis, :].T - V0)**3 * dfe_dV3/6.0 + ( volumes[np.newaxis, :].T - V0)**4*(dfe_dV4/24.0)
                  + fe_V0[:] + energies[np.newaxis, :].T)
 
         fits = [self.eos.fit(volumes, e) for e in tot_en.T]
@@ -463,6 +531,13 @@ def plot_energies(self, tstart=0, tstop=1000, num=1, ax=None, fontsize=10, **kwa
 
         return fig
 
+    def get_phdos_plotter(self) -> PhononDosPlotter:
+        """Build and return a PhononDosPlotter."""
+        plotter = PhononDosPlotter()
+        for volume, phdos in zip(self.volumes_from_phdos, self.phdoses):
+            plotter.add_phdos(f"V={volume:.2f}", phdos)
+        return plotter
+
     @add_fig_kwargs
     def plot_vol_vs_t(self, tstart=0, tstop=1000, num=101, fontsize=10, ax=None, **kwargs) -> Figure:
         """
@@ -613,28 +688,28 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
         dV = volumes[iv0+1]-volumes[iv0]
 
         if self.scale_points == "S":
-            vol,fits = self.vol_E2Vib1(num=num,tstop=tstop,tstart=tstart)
+            vol, fits = self.vol_E2Vib1(num=num, tstop=tstop, tstart=tstart)
             E2D = self.second_derivative_energy_v(self.volumes[self.iv0])
             #f = self.fit_tot_energies(0, 0, 1 ,self.energies[np.newaxis, :].T,self.volumes)
             #E2D = f.F2D
             if tref is None:
                 alpha_1 = - 1/vol[:] * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D
             else:
-                vol_ref,fits = self.vol_E2Vib1(num=1,tstop=tref,tstart=tref)
+                vol_ref, fits = self.vol_E2Vib1(num=1, tstop=tref, tstart=tref)
                 b0 = np.array([fit.b0 for fit in fits])
                 print("B (E2vib1)   @ ", tref, " K =", b0 * 160.21766208, "(GPa)" )
                 alpha_1 = - 1/vol_ref * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D
-            ax.plot(tmesh, alpha_1,color='b', lw=2, label="E2Vib1")
+            ax.plot(tmesh, alpha_1, color='b', lw=2, label="E2Vib1")
             data_to_save = np.column_stack((data_to_save,alpha_1))
             columns.append( 'E2vib1')
 
         if len(self.index_list) >= 2:
-            vol2, fits = self.vol_Einf_Vib1(num=num,tstop=tstop,tstart=tstart)
+            vol2, fits = self.vol_Einf_Vib1(num=num, tstop=tstop, tstart=tstart)
             E2D_V = self.second_derivative_energy_v(vol2)
             if tref is None:
                 alpha_2 = - 1/vol2[:] * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D_V[:]
             else:
-                vol2_ref, fits = self.vol_Einf_Vib1(num=1,tstop=tref,tstart=tref)
+                vol2_ref, fits = self.vol_Einf_Vib1(num=1, tstop=tref, tstart=tref)
                 b0 = np.array([fit.b0 for fit in fits])
                 print("B (Einfvib1) @ ", tref," K =", b0 * 160.21766208, "(GPa)" )
                 alpha_2 = - 1/vol2_ref * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D_V[:]
@@ -644,7 +719,7 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
             columns.append( 'Einfvib1')
 
         if len(self.index_list) >= 3:
-            vol3,fits = self.vol_Einf_Vib2(num=num,tstop=tstop,tstart=tstart)
+            vol3, fits = self.vol_Einf_Vib2(num=num, tstop=tstop, tstart=tstart)
             E2D_V = self.second_derivative_energy_v(vol3)
             dfe_dV2 = np.zeros( num)
             for i, e in enumerate(ph_energies.T):
@@ -655,23 +730,23 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
             if tref is None:
                 alpha_3 = - 1/vol3[:] * ds_dv / (E2D_V[:]+dfe_dV2[:])
             else:
-                vol3_ref,fits = self.vol_Einf_Vib2(num=1,tstop=tref,tstart=tref)
+                vol3_ref,fits = self.vol_Einf_Vib2(num=1, tstop=tref, tstart=tref)
                 b0 = np.array([fit.b0 for fit in fits])
                 #print("B (Einfvib2) @ ",tref," K =", b0*160.21766208, "(GPa)" )
                 alpha_3 = - 1/vol3_ref * ds_dv / (E2D_V[:]+dfe_dV2[:])
 
             ax.plot(tmesh, alpha_3,color='m', lw=2 ,  label=r"$E_\infty Vib2$")
-            data_to_save = np.column_stack((data_to_save,alpha_3))
+            data_to_save = np.column_stack((data_to_save, alpha_3))
             columns.append( 'Einfvib2')
 
         if len(self.index_list) == 5:
             alpha_qha = self.get_thermal_expansion_coeff(tstart, tstop, num, tref)
-            vol4,fits = self.vol_Einf_Vib4(num=num,tstop=tstop,tstart=tstart)
+            vol4,fits = self.vol_Einf_Vib4(num=num, tstop=tstop, tstart=tstart)
             E2D_V = self.second_derivative_energy_v(vol4)
 
-            d2fe_dV2 = np.zeros( num)
-            d3fe_dV3 = np.zeros( num)
-            d4fe_dV4 = np.zeros( num)
+            d2fe_dV2 = np.zeros(num)
+            d3fe_dV3 = np.zeros(num)
+            d4fe_dV4 = np.zeros(num)
             for i,e in enumerate(ph_energies.T):
                 d2fe_dV2[i] = (-e[4]+16*e[3]-30*e[2]+16*e[1]-e[0])/(12*dV**2)
                 d3fe_dV3[i] = (e[4]-2*e[3]+2*e[1]-e[0])/(2*dV**3)
@@ -685,7 +760,7 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
             if tref is None:
                 alpha_4 = - 1/vol4[:] * ds_dv / D2F
             else:
-                vol4_ref,fits = self.vol_Einf_Vib4(num=1,tstop=tref,tstart=tref)
+                vol4_ref,fits = self.vol_Einf_Vib4(num=1, tstop=tref, tstart=tref)
                 b0 = np.array([fit.b0 for fit in fits])
                 print("B (Einfvib4) @ ", tref, " K =", b0 * 160.21766208, "(GPa)" )
                 alpha_4 = - 1/vol4_ref * ds_dv / D2F
@@ -711,13 +786,13 @@ def get_abc(self, volumes, num=101) -> tuple:
             volumes:
         """
         param = np.zeros((num, 4))
-        param = np.polyfit(self.volumes, self.lattice_a , 3)
+        param = np.polyfit(self.volumes, self.lattice_a, 3)
         pa = np.poly1d(param)
         aa_qha = pa(volumes)
-        param = np.polyfit(self.volumes, self.lattice_b , 3)
+        param = np.polyfit(self.volumes, self.lattice_b, 3)
         pb = np.poly1d(param)
         bb_qha = pb(volumes)
-        param = np.polyfit(self.volumes, self.lattice_c , 3)
+        param = np.polyfit(self.volumes, self.lattice_c, 3)
         pc = np.poly1d(param)
         cc_qha = pc(volumes)
 
@@ -730,14 +805,14 @@ def get_angles(self, volumes, num=101) -> tuple:
         Args:
             num: int, optional Number of samples to generate.
         """
-        param = np.zeros((num,4))
+        param = np.zeros((num, 4))
         param = np.polyfit(self.volumes, self.angles_alpha, 3)
         pa = np.poly1d(param)
         gamma = pa(volumes)
-        param = np.polyfit(self.volumes, self.angles_beta , 3)
+        param = np.polyfit(self.volumes, self.angles_beta, 3)
         pb = np.poly1d(param)
         beta = pb(volumes)
-        param = np.polyfit(self.volumes, self.angles_gama , 3)
+        param = np.polyfit(self.volumes, self.angles_gamma, 3)
         pc = np.poly1d(param)
         alpha = pc(volumes)
 
@@ -770,29 +845,29 @@ def plot_thermal_expansion_coeff_abc(self, tstart=0, tstop=1000, num=101, tref=N
         columns = ['#Tmesh']
 
         if self.scale_points == "S":
-            vol2 ,fits = self.vol_E2Vib1(num=num, tstop=tstop, tstart=tstart)
+            vol2, fits = self.vol_E2Vib1(num=num, tstop=tstop, tstart=tstart)
             if tref is not None:
                 vol2_tref, fits = self.vol_E2Vib1(num=1, tstop=tref, tstart=tref)
 
         if len(self.index_list) == 2:
             vol, fits = self.vol_Einf_Vib1(num=num, tstop=tstop, tstart=tstart)
             if tref is not None:
-                vol_tref, fits  = self.vol_Einf_Vib1(num=1, tstop=tref, tstart=tref)
+                vol_tref, fits = self.vol_Einf_Vib1(num=1, tstop=tref, tstart=tref)
             method = r"$ (E_\infty Vib1)$"
 
         if len(self.index_list) == 3:
-            vol,fits = self.vol_Einf_Vib2(num=num, tstop=tstop, tstart=tstart)
+            vol, fits = self.vol_Einf_Vib2(num=num, tstop=tstop, tstart=tstart)
             if tref is not None:
                 vol_tref, fits = self.vol_Einf_Vib2(num=1, tstop=tref, tstart=tref)
             method = r"$ (E_\infty Vib2)$"
 
         if len(self.index_list) == 5:
             tot_en = self.energies_pdos[np.newaxis, :].T + ph_energies
-            f0 = self.fit_tot_energies(tstart, tstop, num,tot_en, self.volumes_from_phdos)
+            f0 = self.fit_tot_energies(tstart, tstop, num, tot_en, self.volumes_from_phdos)
             method = r"$ (E_\infty Vib4)$"
             vol, fits = self.vol_Einf_Vib4(num=num, tstop=tstop, tstart=tstart)
             if tref is not None:
-                vol_tref, fits = self.vol_Einf_Vib4(num=1,tstop=tref,tstart=tref)
+                vol_tref, fits = self.vol_Einf_Vib4(num=1, tstop=tref, tstart=tref)
 
         #alpha  = self.get_thermal_expansion_coeff(tstart, tstop, num, tref)
         tmesh = np.linspace(tstart, tstop, num)
@@ -839,9 +914,9 @@ def plot_thermal_expansion_coeff_abc(self, tstart=0, tstop=1000, num=101, tref=N
                 alpha2_b = (bb2[2:] - bb2[:-2]) / (2 * dt) / bb2_tref
                 alpha2_c = (cc2[2:] - cc2[:-2]) / (2 * dt) / cc2_tref
 
-            ax.plot(tmesh[1:-1] ,alpha2_a ,  linestyle='dashed' ,  color='r', lw=2 ,label = r"$\alpha_a$"" (E2vib1)")
-            ax.plot(tmesh[1:-1] ,alpha2_b ,  linestyle='dashed' ,  color='b', lw=2 ,label = r"$\alpha_b$"" (E2vib1)")
-            ax.plot(tmesh[1:-1] ,alpha2_c ,  linestyle='dashed' ,  color='m', lw=2 ,label = r"$\alpha_c$"" (E2vib1)")
+            ax.plot(tmesh[1:-1], alpha2_a, linestyle='dashed', color='r', lw=2, label=r"$\alpha_a$"" (E2vib1)")
+            ax.plot(tmesh[1:-1], alpha2_b, linestyle='dashed', color='b', lw=2, label=r"$\alpha_b$"" (E2vib1)")
+            ax.plot(tmesh[1:-1], alpha2_c, linestyle='dashed', color='m', lw=2, label=r"$\alpha_c$"" (E2vib1)")
             data_to_save = np.column_stack((data_to_save,alpha2_a,alpha2_b,alpha2_c))
             columns.append( 'E2vib1 (alpha_a,alpha_b,alpha_c)   ')
 
@@ -910,9 +985,9 @@ def plot_thermal_expansion_coeff_angles(self, tstart=0, tstop=1000, num=101, tre
         if tref is not None:
             alpha_tref, beta_tref, cc_tref = self.get_angles(vol_tref, num=1)
 
-        alpha_alpha = np.zeros( num-2)
-        alpha_beta = np.zeros( num-2)
-        alpha_gamma = np.zeros( num-2)
+        alpha_alpha = np.zeros(num-2)
+        alpha_beta = np.zeros(num-2)
+        alpha_gamma = np.zeros(num-2)
         if tref is None:
             alpha_alpha = (alpha[2:] - alpha[:-2]) / (2 * dt) / alpha[1:-1]
             alpha_beta = (beta[2:] - beta[:-2]) / (2 * dt) / beta[1:-1]
@@ -1008,7 +1083,7 @@ def plot_abc_vs_t(self, tstart=0, tstop=1000, num=101, lattice=None, tref=None,
 
         aa, bb, cc = self.get_abc(vol, num=num)
 
-        method_header = method+"  (a,b,c) |"
+        method_header = method + "  (a,b,c) |"
         data_to_save = np.column_stack((data_to_save, aa, bb, cc))
         columns.append(method_header)
 
@@ -1083,7 +1158,7 @@ def plot_angles_vs_t(self, tstart=0, tstop=1000, num=101, angle=None, tref=None,
 
         alpha, beta, gamma = self.get_angles(vol, num=num)
 
-        method_header = method+"  (alpha,beta,gamm) |"
+        method_header = method + "  (alpha,beta,gamm) |"
         data_to_save = np.column_stack((data_to_save, alpha, beta, gamma))
         columns.append(method_header)
 
@@ -1103,7 +1178,6 @@ def plot_angles_vs_t(self, tstart=0, tstop=1000, num=101, angle=None, tref=None,
                 ax.plot(tmesh, gamma2, linestyle='dashed', color='m', lw=2, label=r"$gamma(V(T))$""E2vib1")
 
         set_grid_legend(ax, fontsize, xlabel='T (K)', ylabel=None)
-
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
@@ -1152,11 +1226,16 @@ def fit_forth(self, tstart=0, tstop=1000, num=1, energy="energy", volumes="volum
             min_en[j] = p(min_vol[j])
             F2D_V[j] = p3(min_vol[j])
 
-        return dict2namedtuple(min_vol=min_vol, temp=tmesh , min_en=min_en , param=param , F2D_V=F2D_V) #, fits=fits)
+        return dict2namedtuple(min_vol=min_vol, temp=tmesh, min_en=min_en, param=param, F2D_V=F2D_V) #, fits=fits)
 
     def vol_E2Vib1_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
         """
-        Returns: Vol
+        Args:
+            tstart: The starting value (in Kelvin) of the temperature mesh.
+            tstop: The end value (in Kelvin) of the mesh.
+            num: int, optional. Number of samples to generate.
+
+        Returns: array with num volumes
         """
         volumes0 = self.volumes_from_phdos
         iv0 = self.iv0_vib
@@ -1184,7 +1263,13 @@ def vol_E2Vib1_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
 
     def vol_EinfVib1_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
         """
-        Returns: Vol
+
+        Args:
+            tstart: The starting value (in Kelvin) of the temperature mesh.
+            tstop: The end value (in Kelvin) of the mesh.
+            num: int, optional. Number of samples to generate.
+
+        Returns: array with num volumes
         """
         volumes0 = self.volumes_from_phdos
         iv0 = self.iv0_vib
@@ -1216,7 +1301,7 @@ def vol_Einf_Vib2_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
             tstop: The end value (in Kelvin) of the mesh.
             num: int, optional Number of samples to generate.
 
-        Returns: Vol
+        Returns: array with num volumes
         """
         tmesh = np.linspace(tstart, tstop, num)
 
@@ -1225,10 +1310,10 @@ def vol_Einf_Vib2_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
         energies = self.energies
         volumes0 = self.volumes_from_phdos
         volumes = self.volumes
-        vol = np.zeros( num)
-        dfe_dV = np.zeros( num)
-        d2fe_dV2 = np.zeros( num)
-        fe_V0 = np.zeros( num)
+        vol = np.zeros(num)
+        dfe_dV = np.zeros(num)
+        d2fe_dV2 = np.zeros(num)
+        fe_V0 = np.zeros(num)
         if len(self.index_list) == 5:
             iv0 = 2
         else:
@@ -1255,9 +1340,9 @@ def vol_Einf_Vib4_forth(self, tstart=0, tstop=1000, num=101) -> np.ndarray:
         Args:
             tstart: The starting value (in Kelvin) of the temperature mesh.
             tstop: The end value (in Kelvin) of the mesh.
-            num: int, optional Number of samples to generate.
+            num: int, optional. Number of samples to generate.
 
-        Returns: Vol
+        Returns: array with num volumes
         """
         tmesh = np.linspace(tstart, tstop, num)
 
@@ -1319,19 +1404,19 @@ def plot_vol_vs_t_4th(self, tstart=0, tstop=1000, num=101, ax=None, fontsize=10,
 
         if self.scale_points == "S":
             vol_4th = self.vol_E2Vib1_forth(num=num, tstop=tstop, tstart=tstart)
-            ax.plot(tmesh, vol_4th,color='b', lw=2, label="E2Vib1")
+            ax.plot(tmesh, vol_4th, color='b', lw=2, label="E2Vib1")
             data_to_save = np.column_stack((data_to_save, vol_4th))
             columns.append( 'E2vib1')
 
         if len(self.index_list) >= 2:
             vol2_4th = self.vol_EinfVib1_forth(num=num, tstop=tstop, tstart=tstart)
-            ax.plot(tmesh, vol2_4th,color='gold', lw=2 ,  label=r"$E_\infty Vib1$")
+            ax.plot(tmesh, vol2_4th, color='gold', lw=2, label=r"$E_\infty Vib1$")
             data_to_save = np.column_stack((data_to_save, vol2_4th))
             columns.append( 'Einfvib1')
 
         if len(self.index_list) >= 3:
             vol3_4th = self.vol_Einf_Vib2_forth(num=num,tstop=tstop,tstart=tstart)
-            ax.plot(tmesh, vol3_4th,color='m', lw=2 , label=r"$E_\infty Vib2$")
+            ax.plot(tmesh, vol3_4th, color='m', lw=2, label=r"$E_\infty Vib2$")
             data_to_save = np.column_stack((data_to_save, vol3_4th))
             columns.append( 'Einfvib2')
 
@@ -1340,12 +1425,12 @@ def plot_vol_vs_t_4th(self, tstart=0, tstop=1000, num=101, ax=None, fontsize=10,
             f0 = self.fit_forth( tstart, tstop, num ,tot_en,volumes)
             vol4_4th = self.vol_Einf_Vib4_forth(num=num, tstop=tstop, tstart=tstart)
             ax.plot(tmesh, vol4_4th,color='c', lw=2  ,label=r"$E_\infty Vib4$")
-            ax.plot(tmesh, f0.min_vol, color='k', linestyle='dashed', lw=1.5 ,label="QHA")
+            ax.plot(tmesh, f0.min_vol, color='k', linestyle='dashed', lw=1.5, label="QHA")
             data_to_save = np.column_stack((data_to_save, vol4_4th, f0.min_vol))
             columns.append( 'Einfvib4')
             columns.append( 'QHA')
 
-        ax.plot(0, self.volumes[self.iv0], color='g', lw=0, marker='o', ms=10,label="V0")
+        ax.plot(0, self.volumes[self.iv0], color='g', lw=0, marker='o', ms=10, label="V0")
         set_grid_legend(ax, fontsize, xlabel='T (K)', ylabel=r'V (${\AA}^3$)')
         ax.set_xlim(tstart, tstop)
 
@@ -1383,8 +1468,8 @@ def get_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=No
         gamma = np.zeros(num)
 
         for j in range(num):
-            param[j] = np.polyfit(self.volumes_from_phdos,df_t[j] , 3)
-            param2[j] = np.array([3*param[j][0],2*param[j][1],param[j][2]])
+            param[j] = np.polyfit(self.volumes_from_phdos, df_t[j], 3)
+            param2[j] = np.array([3*param[j][0], 2*param[j][1], param[j][2]])
             p = np.poly1d(param2[j])
             d2f_t_v[j] = p(f.min_vol[j])
 
@@ -1455,7 +1540,7 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
             if tref is None:
                 alpha_2 = - 1/vol2_4th[:] * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D_V[:]
             else:
-                vol2_4th_ref = self.vol_EinfVib1_forth(num=1,tstop=tref,tstart=tref)
+                vol2_4th_ref = self.vol_EinfVib1_forth(num=1, tstop=tref, tstart=tref)
                 alpha_2 = - 1/vol2_4th_ref * (df_t[:,iv1]-df_t[:,iv0])/(volumes[iv1]-volumes[iv0]) / E2D_V[:]
 
             ax.plot(tmesh, alpha_2,color='gold', lw=2 ,  label=r"$E_\infty Vib1$")
@@ -1465,7 +1550,7 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
         if len(self.index_list) >= 3:
             vol3_4th = self.vol_Einf_Vib2_forth(num=num, tstop=tstop, tstart=tstart)
             E2D_V = p3(vol3_4th)
-            dfe_dV2 = np.zeros( num)
+            dfe_dV2 = np.zeros(num)
             for i,e in enumerate(ph_energies.T):
                 dfe_dV2[i] = (e[iv0+2]-2.0*e[iv0+1]+e[iv0])/(dV)**2
 
@@ -1485,9 +1570,9 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
             vol4_4th = self.vol_Einf_Vib4_forth(num=num, tstop=tstop, tstart=tstart)
             E2D_V = p3(vol4_4th)
 
-            d2fe_dV2 = np.zeros( num)
-            d3fe_dV3 = np.zeros( num)
-            d4fe_dV4 = np.zeros( num)
+            d2fe_dV2 = np.zeros(num)
+            d3fe_dV3 = np.zeros(num)
+            d4fe_dV4 = np.zeros(num)
             for i,e in enumerate(ph_energies.T):
                 d2fe_dV2[i] = (-e[4]+16*e[3]-30*e[2]+16*e[1]-e[0])/(12*dV**2)
                 d3fe_dV3[i] = (e[4]-2*e[3]+2*e[1]-e[0])/(2*dV**3)
@@ -1508,10 +1593,10 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
             ax.plot(tmesh, alpha_4, color='c', linewidth=2, label=r"$E_\infty Vib4$")
 
             alpha_qha = self.get_thermal_expansion_coeff_4th(tstart, tstop, num, tref)
-            ax.plot(tmesh, alpha_qha, color='k',linestyle='dashed', lw=1.5 ,label="QHA")
+            ax.plot(tmesh, alpha_qha, color='k', linestyle='dashed', lw=1.5, label="QHA")
             data_to_save = np.column_stack((data_to_save, alpha_4, alpha_qha))
-            columns.append( 'Einfvib4')
-            columns.append( 'QHA')
+            columns.append('Einfvib4')
+            columns.append('QHA')
 
         set_grid_legend(ax, fontsize, xlabel='T (K)', ylabel=r'$\alpha$ (K$^{-1}$)')
         ax.set_xlim(tstart, tstop)
@@ -1562,7 +1647,7 @@ def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None, tref=No
 
         if len(self.index_list) == 5:
             tot_en = self.energies_pdos[np.newaxis, :].T + ph_energies
-            f0 = self.fit_forth( tstart, tstop, num ,tot_en,volumes)
+            f0 = self.fit_forth(tstart, tstop, num, tot_en, volumes)
             method = r"$ (E_\infty Vib4)$"
             vol = self.vol_Einf_Vib4_forth(num=num, tstop=tstop, tstart=tstart)
 
@@ -1573,11 +1658,11 @@ def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None, tref=No
         columns.append(method_header)
 
         if lattice is None or lattice == "a":
-            ax.plot(tmesh, aa , color='r', lw=2, label=r"$a(V(T))$" + method, **kwargs)
+            ax.plot(tmesh, aa, color='r', lw=2, label=r"$a(V(T))$" + method, **kwargs)
         if lattice is None or lattice == "b":
-            ax.plot(tmesh, bb , color='b', lw=2, label=r"$b(V(T))$" + method)
+            ax.plot(tmesh, bb, color='b', lw=2, label=r"$b(V(T))$" + method)
         if lattice is None or lattice == "c":
-            ax.plot(tmesh, cc , color='m', lw=2, label=r"$c(V(T))$" + method)
+            ax.plot(tmesh, cc, color='m', lw=2, label=r"$c(V(T))$" + method)
 
         if abs(abs(self.volumes[self.iv0]-volumes[iv0])-abs(volumes[iv1]-self.volumes[self.iv0])) < 1e-3 :
             if lattice is None or lattice == "a":
@@ -1636,13 +1721,13 @@ def plot_angles_vs_t_4th(self, tstart=0, tstop=1000, num=101, angle=None, tref=N
 
         if len(self.index_list) == 5:
             tot_en = self.energies_pdos[np.newaxis, :].T + ph_energies
-            f0 = self.fit_forth( tstart, tstop, num ,tot_en,volumes)
+            f0 = self.fit_forth(tstart, tstop, num, tot_en, volumes)
             method = r"$ (E_\infty Vib4)$"
             vol = self.vol_Einf_Vib4_forth(num=num, tstop=tstop, tstart=tstart)
 
         alpha, beta, gamma = self.get_angles(vol, num=num)
 
-        method_header = method+"  (alpha,beta,gamma) |"
+        method_header = method + "  (alpha,beta,gamma) |"
         data_to_save = np.column_stack((data_to_save, alpha, beta, gamma))
         columns.append(method_header)
 
@@ -1726,9 +1811,9 @@ def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tr
         if tref is not None:
             aa_tref, bb_tref, cc_tref = self.get_abc(vol_tref, num=1)
 
-        alpha_a = np.zeros( num-2)
-        alpha_b = np.zeros( num-2)
-        alpha_c = np.zeros( num-2)
+        alpha_a = np.zeros(num-2)
+        alpha_b = np.zeros(num-2)
+        alpha_c = np.zeros(num-2)
         if tref is None:
             alpha_a = (aa[2:] - aa[:-2]) / (2 * dt) / aa[1:-1]
             alpha_b = (bb[2:] - bb[:-2]) / (2 * dt) / bb[1:-1]
@@ -1742,7 +1827,7 @@ def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tr
         ax.plot(tmesh[1:-1], alpha_b, color='b', lw=2, label=r"$\alpha_b$" + method)
         ax.plot(tmesh[1:-1], alpha_c, color='m', lw=2, label=r"$\alpha_c$" + method)
 
-        method_header = method+"  (alpha_a,alpha_b,alpha_c) |"
+        method_header = method + "  (alpha_a,alpha_b,alpha_c) |"
         data_to_save = np.column_stack((data_to_save, alpha_a, alpha_b, alpha_c))
         columns.append(method_header)
 
@@ -1751,9 +1836,9 @@ def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tr
             if tref is not None:
                 aa2_tref, bb2_tref, cc2_tref = self.get_abc(vol2_tref, num=1)
 
-            alpha2_a = np.zeros( num-2)
-            alpha2_b = np.zeros( num-2)
-            alpha2_c = np.zeros( num-2)
+            alpha2_a = np.zeros(num-2)
+            alpha2_b = np.zeros(num-2)
+            alpha2_c = np.zeros(num-2)
             if tref is None:
                 alpha2_a = (aa2[2:] - aa2[:-2]) / (2 * dt) / aa2[1:-1]
                 alpha2_b = (bb2[2:] - bb2[:-2]) / (2 * dt) / bb2[1:-1]
diff --git a/abipy/dfpt/tests/test_qha_approximation.py b/abipy/dfpt/tests/test_qha_approximation.py
index 563007b93..0c9675f39 100644
--- a/abipy/dfpt/tests/test_qha_approximation.py
+++ b/abipy/dfpt/tests/test_qha_approximation.py
@@ -9,8 +9,8 @@
 
 class QhaTest(AbipyTest):
 
-    def test_v_ZSIZA(self):
-
+    def test_v_ZSISA(self):
+        """Testing v_SSISA postprocessing tools."""
         # Root points to the directory in the git submodule with the output results.
         root = os.path.join(abidata.dirpath, "data_v-ZSISA-QHA.git", "Si_v_ZSISA_approximation")
 
@@ -21,13 +21,17 @@ def test_v_ZSIZA(self):
 
         gsr_paths = [os.path.join(root, "scale_{:d}_GSR.nc".format(s)) for s in strains]
         ddb_paths = [os.path.join(root, "scale_{:d}_GSR_DDB".format(s)) for s in strains]
-        dos_paths = [os.path.join(root, "scale_{:d}_PHDOS.nc".format(s)) for s in strains2]
+        phdos_paths = [os.path.join(root, "scale_{:d}_PHDOS.nc".format(s)) for s in strains2]
 
-        qha = QHA_App.from_files_app_ddb(ddb_paths, dos_paths)
+        qha = QHA_App.from_ddb_phdos_files(ddb_paths, phdos_paths)
         tstart, tstop = 0, 800
 
         # Test basic properties and get methods of qha
         assert qha.nvols == 5
+        assert qha.eos_name == "vinet"
+        assert qha.scale_points == "D"
+        assert str(qha)
+        assert qha.to_string(verbose=1)
 
         data = qha.fit_tot_energies(tstart=0, tstop=0, num=1,
                                    tot_energies=qha.energies[np.newaxis, :].T, volumes=qha.volumes)
@@ -120,3 +124,6 @@ def test_v_ZSIZA(self):
             assert qha.plot_thermal_expansion_coeff_4th(tref=293, show=False)
             assert qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293, show=False)
             assert qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3, show=False)
+
+            plotter = qha.get_phdos_plotter()
+            plotter.combiplot(show=False)
diff --git a/abipy/electrons/gw.py b/abipy/electrons/gw.py
index 02128a0bb..c2f9a3a20 100644
--- a/abipy/electrons/gw.py
+++ b/abipy/electrons/gw.py
@@ -2464,12 +2464,21 @@ class GwRobotWithDisplacedAtom(SigresRobot):
     Specialized class to analyze GW or GWR calculations with displaced atom.
     """
     @classmethod
-    def from_displaced_atom(cls, site_index, reduced_dir, step_ang, gw_files):
+    def from_displaced_atom(cls, site_index, reduced_dir, step_ang, gw_files) -> GwRobotWithDisplacedAtom:
         """
-        Build an instance a list of SIG.nc or GWR.nc files.
+        Build an instance from a list of SIGRES.nc or GWR.nc files.
+
+        Args:
+            site_index: Index of the site that has been displaced.
+            reduced_dir: Reduced direction of the displacement.
+            step_ang: Step used to displace structures in Angstrom.
+            gw_files: List of paths to either SIGRES.nc or GWR.nc files.
+                Files are assumed to be ordered according to the displacement.
         """
         #print(f"{gw_files}")
         new = cls(*gw_files)
+        print("new = cls(*gw_files) done!")
+
         new.site_index = site_index
         new.reduced_dir = reduced_dir
         new.step_ang = step_ang
@@ -2478,11 +2487,13 @@ def from_displaced_atom(cls, site_index, reduced_dir, step_ang, gw_files):
         i0 = new.num_points // 2
         origin_structure = new[i0].structure
 
-        # Consistency check:
+        ####################
+        # Consistency check
+        ####################
+
         # 1) Make sure lattice parameters and all sites other than i0 are equal.
         fixed_site_indices = [i for i in range(len(origin_structure)) if i != site_index]
         if err_str := new.has_different_structures(site_indices=fixed_site_indices):
-            #diff_structures(structures, fmt="abivars", mode="table", headers=(), file=sys.stdout)
             raise ValueError(err_str)
 
         # TODO
@@ -2503,9 +2514,15 @@ def from_displaced_atom(cls, site_index, reduced_dir, step_ang, gw_files):
 
         return new
 
-    def get_dataframe_skb(self, spin, kpoint, band, with_params: bool=True) -> pd.DataFrame:
+    def get_dataframe_skb(self, spin, kpoint, band, with_params: bool = True) -> pd.DataFrame:
         """
-        Return pandas dataframe with the most important results.
+        Return a pandas dataframe with the most important results.
+
+        Args:
+            spin: Spin index.
+            kpoint: K-point in self-energy. Accepts |Kpoint|, vector or index.
+            band: band index.
+            with_params: True if metadata should be included.
         """
         # Create list of QPState for each file.
         qp_states = [ncfile.r.read_qplist_sk(spin, kpoint, band=band)[0] for ncfile in self.abifiles]
diff --git a/abipy/eph/varpeq.py b/abipy/eph/varpeq.py
index 372c43d81..7b3e6a2de 100644
--- a/abipy/eph/varpeq.py
+++ b/abipy/eph/varpeq.py
@@ -86,10 +86,10 @@
 
 @dataclasses.dataclass(kw_only=True)
 class Entry:
-    name: str
-    latex: str
-    info: str
-    utype: str
+    name: str   # Entry name
+    latex: str  # Latex label
+    info: str   # Description string
+    utype: str  # Unit type
     #color: str
 
 
@@ -771,7 +771,7 @@ def plot_bqnu_with_phbands(self, phbands_qpath,
 
         if not with_phdos:
             # Return immediately.
-            if with_title: 
+            if with_title:
                 fig.suptitle(self.get_title(with_gaps=True))
             return fig
 
@@ -1044,4 +1044,4 @@ def write_notebook(self, nbpath=None) -> str:
             #nbv.new_code_cell("ebands_plotter = robot.get_ebands_plotter()"),
         ])
 
-        return self._write_nb_nbpath(nb, nbpath)
\ No newline at end of file
+        return self._write_nb_nbpath(nb, nbpath)
diff --git a/abipy/examples/flows/_under_develop/run_qha.py b/abipy/examples/flows/run_qha_vzsisa.py
similarity index 64%
rename from abipy/examples/flows/_under_develop/run_qha.py
rename to abipy/examples/flows/run_qha_vzsisa.py
index 8c487ea30..bd9ff5a73 100755
--- a/abipy/examples/flows/_under_develop/run_qha.py
+++ b/abipy/examples/flows/run_qha_vzsisa.py
@@ -10,7 +10,7 @@
 import abipy.data as abidata
 
 from abipy import flowtk
-from abipy.flowtk.qha import QhaFlow
+from abipy.flowtk.qha import vZSISAFlow
 
 
 def build_flow(options):
@@ -19,7 +19,7 @@ def build_flow(options):
     """
     # Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_")
     if not options.workdir:
-        __file__ = os.path.join(os.getcwd(), "run_qha.py")
+        __file__ = os.path.join(os.getcwd(), "run_qha_vzsisa.py")
         options.workdir = os.path.basename(__file__).replace(".py", "").replace("run_", "flow_")
 
     # Initialize structure and pseudos
@@ -27,14 +27,24 @@ def build_flow(options):
     pseudos = abidata.pseudos("14si.pspnc")
 
     # Build input for GS calculation.
+    #ngkpt = [2, 2, 2]; ngqpt = [2, 2, 2]
+    #ngkpt = [4, 4, 4]; ngqpt = [4, 4, 4]
+    ngkpt = [2, 2, 2]; ngqpt = [1, 1, 1]
+    with_becs = False
+    with_quad = False
+    #with_quad = not structure.has_zero_dynamical_quadrupoles
+
+    #bo_scales = [0.96, 0.98, 1.0, 1.02, 1.04, 1.06]
+    #ph_scales = [0.98, 1.0, 1.02, 1.04, 1.06] # EinfVib4(D)
+    bo_scales = [0.96, 0.98, 1, 1.02, 1.04] # EinfVib4(S)
+    ph_scales = [1, 1.02, 1.04]         # EinfVib2(D)
+
     scf_input = abilab.AbinitInput(structure, pseudos)
-    scf_input.set_vars(ecut=12, nband=4, tolvrs=1e-8)
-    scf_input.set_kmesh(ngkpt=[2, 2, 2], shiftk=[0, 0, 0])
+    scf_input.set_vars(ecut=8, nband=4, tolvrs=1e-8, nstep=50) #, paral_kgb=0)
+    scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0])
 
-    with_becs = False
-    ngqpt = [2, 2, 2]
-    return QhaFlow.from_scf_input(options.workdir, scf_input, ngqpt, with_becs, eps=0.005,
-                                  ) # , edos_ngkpt=(4, 4, 4))
+    return vZSISAFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt,
+                                     with_becs, with_quad, edos_ngkpt=None)
 
 
 # This block generates the thumbnails in the Abipy gallery.
diff --git a/abipy/examples/plot/plot_qha_v-ZSISA.py b/abipy/examples/plot/plot_qha_v-ZSISA.py
index 9cae10c35..84a7306e7 100755
--- a/abipy/examples/plot/plot_qha_v-ZSISA.py
+++ b/abipy/examples/plot/plot_qha_v-ZSISA.py
@@ -79,3 +79,8 @@
 #%%
 qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3,
                          title="Angles as a function of T")
+
+
+#%%
+phdos_plotter = qha.get_phdos_plotter()
+phdos_plotter.combiplot()
diff --git a/abipy/flowtk/flows.py b/abipy/flowtk/flows.py
index f470612a5..156acf3c9 100644
--- a/abipy/flowtk/flows.py
+++ b/abipy/flowtk/flows.py
@@ -2105,7 +2105,7 @@ def register_work(self, work: Work, deps=None, manager=None, workdir=None) -> Wo
 
         return work
 
-    def register_work_from_cbk(self, cbk_name, cbk_data, deps, work_class, manager=None):
+    def register_work_from_cbk(self, cbk_name, cbk_data, deps, work_class, manager=None) -> Work:
         """
         Registers a callback function that will generate the :class:`Task` of the :class:`Work`.
 
diff --git a/abipy/flowtk/launcher.py b/abipy/flowtk/launcher.py
index f24f30add..c219f14a3 100644
--- a/abipy/flowtk/launcher.py
+++ b/abipy/flowtk/launcher.py
@@ -666,8 +666,7 @@ def _runem_all(self) -> None:
                   len(list(flow.iflat_tasks(status=flow.S_SUB))))
 
         if nqjobs >= self.max_njobs_inqueue:
-            print(f"Too many jobs in the queue: {nqjobs} >= {self.max_njobs_inqueue}.\n",
-                  "No job will be submitted.")
+            print(f"Too many jobs in the queue: {nqjobs} >= {self.max_njobs_inqueue}. No job will be submitted.")
             flow.check_status(show=False)
             return
 
diff --git a/abipy/flowtk/qha.py b/abipy/flowtk/qha.py
index aaebffaa5..2fc3dd0bd 100644
--- a/abipy/flowtk/qha.py
+++ b/abipy/flowtk/qha.py
@@ -6,31 +6,109 @@
 
 import numpy as np
 
-from abipy.tools.serialization import mjson_write # mjson_load, mjson_loads,
+from abipy.tools.serialization import mjson_write
 from abipy.dfpt.deformation_utils import generate_deformations
 from abipy.flowtk.works import Work, PhononWork
 from abipy.flowtk.tasks import RelaxTask
 from abipy.flowtk.flows import Flow
 
 
-class RelaxAndAddPhononWorks(Work):
+class vZSISAFlow(Flow):
+    """
+    Flow for QHA calculations with the vZSISA approach
+
+    .. rubric:: Inheritance Diagram
+    .. inheritance-diagram:: vZSISAFlow
+    """
+
+    @classmethod
+    def from_scf_input(cls, workdir, scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad,
+                       edos_ngkpt=None, manager=None) -> vZSISAFlow:
+        """
+        Build a flow for QHA calculations from an |AbinitInput| for GS-SCF calculation.
+
+        Args:
+            workdir: Working directory of the flow.
+            scf_input: |AbinitInput| for GS-SCF run used as template to generate the other inputs.
+            bo_scales:
+            ph_scales:
+            ngqpt: Three integers defining the q-mesh for phonon calculation.
+            with_becs: Activate calculation of Electric field and Born effective charges.
+            with_quad: Activate calculation of dynamical quadrupoles. Require `with_becs`
+                Note that only selected features are compatible with dynamical quadrupoles.
+                Please consult <https://docs.abinit.org/topics/longwave/>
+            edos_ngkpt: Three integers defining the the k-sampling for the computation of the
+                electron DOS with the relaxed structures. Useful for metals or small gap semiconductors
+                in which the electronic contribution should be included.
+                None disables the computation of the e-DOS.
+            manager: |TaskManager| instance. Use default if None.
+        """
+        flow = cls(workdir=workdir, manager=manager)
+
+        work = vZSISAWork.from_scf_input(scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad,
+                                         optcell=2, ionmov=2, edos_ngkpt=edos_ngkpt)
+
+        flow.register_work(work)
+        return flow
+
+    def finalize(self):
+        """
+        This method is called when the flow is completed.
+        It performs some basic post-processing of the results to facilitate further analysis.
+        """
+        work = self[0]
+        data = {}
+
+        # Build list of strings with path to the relevant output files ordered by V.
+        data["gsr_relax_paths"] = [task.gsr_path for task in work.relax_tasks_vol]
+        data["gsr_relax_volumes"] = [task.gsr_path for task in work.relax_tasks_vol]
+        data["ddb_paths"] = [ph_work.outdir.has_abiext("DDB") for ph_work in work.ph_works]
+        data["gsr_edos_path"] = []
+        if work.edos_work:
+            data["gsr_edos_paths"] = [task.gsr_path for task in work.edos_work]
+
+        #entries = []
+        #items = zip(self.relax_and_phonon_work.relax_tasks_vol, self.relax_and_phonon_work.ph_works_vol)
+        #for ivol, (relax_task, ph_work) in enumerate(items):
+        #    ddb_path = ph_work.outdir.has_abiext("DDB")
+        #    with relax_task.open_gsr() as gsr:
+        #        entries.append(dict(
+        #            volume=gsr.structure.volume,
+        #            energy_eV=float(gsr.energy),
+        #            pressure_GPa=float(gsr.pressure),
+        #            structure=gsr.structure,
+        #            gsr_edos_path=None,
+        #            ddb_path=ddb_path,
+        #        )
+
+        #data["entries_for_each_volume"] = entries
+
+        mjson_write(data, self.outdir.path_in("vzsisa.json"), indent=4)
+
+        return super().finalize()
+
+
+class vZSISAWork(Work):
     """
     This work performs a structural relaxation of the initial structure, then a set of distorted
     structures is genenerated and the relaxed structures are used
     to compute phonons, BECS and the dielectric tensor with DFPT.
 
     .. rubric:: Inheritance Diagram
-    .. inheritance-diagram:: RelaxAndAddPhononWorks
+    .. inheritance-diagram:: vZSISAWork
     """
 
     @classmethod
-    def from_scf_input(cls, scf_input, eps, ngqpt, with_becs: bool, optcell: int, ionmov: int,
-                       edos_ngkpt=None) -> RelaxAndAddPhononWorks:
+    def from_scf_input(cls, scf_input, bo_scales, ph_scales, ngqpt,
+                       with_becs: bool, with_quad: bool,
+                       optcell: int, ionmov: int, edos_ngkpt=None) -> vZSISAWork:
         """
         Build the work from an |AbinitInput| representing a GS-SCF calculation.
 
         Args:
             scf_input: |AbinitInput| for GS-SCF used as template to generate the other inputs.
+            bo_scales:
+            ph_scales:
             with_becs: Activate calculation of Electric field and Born effective charges.
             optcell, ionmov: Abinit input variables.
             edos_ngkpt: Three integers defining the the k-sampling for the computation of the
@@ -42,9 +120,11 @@ def from_scf_input(cls, scf_input, eps, ngqpt, with_becs: bool, optcell: int, io
 
         # Save attributes in work
         work.initial_scf_input = scf_input
-        work.eps = eps
+        work.bo_scales = np.array(bo_scales)
+        work.ph_scales = np.array(ph_scales)
         work.ngqpt = ngqpt
         work.with_becs = with_becs
+        work.with_quad = with_quad
         work.edos_ngkpt = edos_ngkpt if edos_ngkpt is None else np.reshape(edos_ngkpt, (3,))
 
         # Create input for relaxation and register the relaxation task.
@@ -55,6 +135,7 @@ def from_scf_input(cls, scf_input, eps, ngqpt, with_becs: bool, optcell: int, io
             relax_template.set_vars_ifnotin(ecutsm=0.5, dilatmx=1.05)
 
         work.initial_relax_task = work.register_relax_task(relax_template)
+
         return work
 
     def on_ok(self, sender):
@@ -63,14 +144,16 @@ def on_ok(self, sender):
         It executes on_all_ok when all tasks in self have reached `S_OK`.
         """
         if sender == self.initial_relax_task:
-            # Get relaxed structure and build deformed structures.
+            # Get relaxed structure and build new task for structural relaxation at fixed volume.
             relaxed_structure = sender.get_final_structure()
-            self.deformed_structures_dict = generate_deformations(relaxed_structure, eps=self.eps)
+            v0 = relaxed_structure.volume
 
-            # Add new tasks to relax the deformed structures at fixed cell.
-            for structure in self.deformed_structures_dict.values():
-                new_input = self.relax_template.new_with_structure(structure, optcell=0)
-                self.register_relax_task(new_input)
+            relax_template = self.relax_template
+            self.relax_tasks_vol = []
+            for bo_scale in self.bo_scales:
+                new_structure = relax_template.structure.scale_lattice(v0 * bo_scale)
+                task = self.register_relax_task(relax_template.new_with_structure(new_structure))
+                self.relax_tasks_vol.append(task)
 
             self.flow.allocate(build=True)
 
@@ -81,23 +164,23 @@ def on_all_ok(self):
         This callback is called when all tasks in the Work reach status `S_OK`.
         Here we add a new PhononWork for each volume using the relaxed structure.
         """
-        if self.edos_ngkpt is not None: self.edos_work = Work()
+        self.edos_work = None if self.edos_ngkpt is None else Work()
+        self.ph_works = []
 
         # Build phonon works for the different relaxed structures.
-        self.ph_works = []
-        start = 1
-        for task, def_name in zip(self[start:], self.deformed_structures_dict.keys(), strict=True):
+        for task, bo_scale in zip(self.relax_tasks_vol, self.bo_scales):
+            if all(abs(bo_scale - self.ph_scales)) > 1e-3: continue
             relaxed_structure = task.get_final_structure()
             scf_input = self.initial_scf_input.new_with_structure(relaxed_structure)
             ph_work = PhononWork.from_scf_input(scf_input, self.ngqpt, is_ngqpt=True, tolerance=None,
-                                                with_becs=self.with_becs, ddk_tolerance=None)
+                                                with_becs=self.with_becs, with_quad=self.with_quad,
+                                                ddk_tolerance=None)
 
-            ph_work.set_name(def_name)
-            self.ph_works.append(ph_work)
             self.flow.register_work(ph_work)
+            self.ph_works.append(ph_work)
 
             # Add electron DOS calculation.
-            if self.edos_ngkpt is not None:
+            if self.edos_work is not None:
                 edos_input = scf_input.make_edos_input(self.edos_ngkpt)
                 self.edos_work.register_nscf_task(edos_input, deps={ph_work[0]: "DEN"})
 
@@ -149,109 +232,3 @@ def on_all_ok(self):
         Implement the post-processing step at the end of the Work.
         """
         return super().on_all_ok()
-
-
-class QhaFlow(Flow):
-    """
-    Flow for QHA calculations.
-
-    .. rubric:: Inheritance Diagram
-    .. inheritance-diagram:: QhaFlow
-    """
-
-    @classmethod
-    def from_scf_input(cls, workdir, scf_input, ngqpt, with_becs, eps=0.005,
-                       edos_ngkpt=None, manager=None) -> QhaFlow:
-        """
-        Build a flow for QHA calculations from an |AbinitInput| for GS-SCF calculations.
-
-        Args:
-            workdir: Working directory of the flow.
-            scf_input: |AbinitInput| for GS-SCF run used as template to generate the other inputs.
-            ngqpt: Three integers defining the q-mesh for phonon calculation.
-            with_becs: Activate calculation of Electric field and Born effective charges.
-            eps:
-            edos_ngkpt: Three integers defining the the k-sampling for the computation of the
-                electron DOS with the relaxed structures. Useful for metals or small gap semiconductors
-                in which the electronic contribution should be included.
-                None disables the computation of the e-DOS.
-            manager: |TaskManager| instance. Use default if None.
-        """
-        flow = cls(workdir=workdir, manager=manager)
-
-        new_work = RelaxAndAddPhononWorks.from_scf_input(scf_input, eps, ngqpt,
-                                                         with_becs=with_becs, optcell=2, ionmov=2,
-                                                         edos_ngkpt=edos_ngkpt)
-        flow.register_work(new_work)
-        flow.relax_and_add_phonon_works = new_work
-
-        return flow
-
-    #def on_all_ok(self):
-    #    if self.on_all_ok_num_calls > 0: return True
-    #    self.on_all_ok_num_calls += 1
-
-    #    print(f"In on_all_ok with {self.on_all_ok_num_calls =}")
-
-    #    # implement_logic_to_create_new_work
-    #    #for work in self.relax_and_add_phonon_works.ph_works
-
-    #    """
-    #    self.register_work(work)
-    #    self.allocate()
-    #    self.build_and_pickle_dump()
-    #    """
-
-    #    #from abipy.dfpt.qha_approximation import QHA_App
-    #    #qha = QHA_App.from_files_app_ddb(cls, ddb_paths, phdos_paths)
-
-    #    # The scheduler will keep on running the flow.
-    #    return False
-
-    def finalize(self):
-        """
-        This method is called when the flow is completed.
-        It performs some basic post-processing of the results to facilitate further analysis.
-        """
-        data = {}
-        #work = self[0]
-        #data = dict(
-        #    initial_structure=work.initial_structure,
-        #    relaxed_structure=work.final_structure,
-        #    relaxed_structure_energy_eV=work.initial_structure,
-        #    initial_deformed_structures=work.initial_structure,
-        #    relaxed_deformed_structures_energy_ev=work.initial_structure,
-        #)
-
-        # QHA.from_files(cls, gsr_paths, phdos_paths)
-        # Build list of strings with path to the relevant output files ordered by V.
-        #data["gsr_relax_paths"] = [task.gsr_path for task in self.relax_and_phonon_work.relax_tasks_vol]
-        data["ddb_paths"] = [ph_work.outdir.has_abiext("DDB") for ph_work in self.relax_and_add_phonon_works.ph_works]
-        #data["gsr_edos_path"] = []
-        #if self.relax_and_phonon_work.edos_ngkpt is not None:
-        #    data["gsr_edos_paths"] = [task.gsr_path for task in self.relax_and_phonon_work.edos_work]
-
-        #entries = []
-        #items = zip(self.relax_and_phonon_work.relax_tasks_vol, self.relax_and_phonon_work.ph_works_vol)
-        #for ivol, (relax_task, ph_work) in enumerate(items):
-        #    ddb_path = ph_work.outdir.has_abiext("DDB")
-        #    with relax_task.open_gsr() as gsr:
-        #        entry = dict(
-        #            volume=gsr.structure.volume,
-        #            energy_eV=gsr.energy,
-        #            pressure_GPa=gsr.pressure,
-        #            structure=gsr.structure,
-        #            gsr_edos_path=None,
-        #            ddb_path=ddb_path,
-        #        )
-
-        #    if self.relax_and_phonon_work.edos_ngkpt is not None:
-        #        task = self.relax_and_phonon_work.edos_work[ivol]
-        #        entry["gsr_edos_path"] = task.gsr_path
-
-        #    entries.append(entry)
-        #data["entries_for_each_volume"] = entries
-
-        mjson_write(data, self.outdir.path_in("qha_data.json"), indent=4)
-
-        return super().finalize()
diff --git a/abipy/ml/extxyz_generator.py b/abipy/ml/extxyz_generator.py
index 650aa274d..c17c2741e 100644
--- a/abipy/ml/extxyz_generator.py
+++ b/abipy/ml/extxyz_generator.py
@@ -25,6 +25,9 @@
 
 class ExtxyzIOWriter:
     """
+    This object converts output files produced by ab-initio codes such as ABINIT or VASP
+    to extended XYZ files that can be used to train ML models.
+
     Example:
 
         # To find all the vasprun.xml files starting from a top-level directory, use:
@@ -113,6 +116,7 @@ def yield_atoms(self):
         Generate ASE atoms.
         """
         for filepath in self.filepaths:
+
             if self.ext == "vasprun.xml":
                 vasprun = Vasprun(filepath)
                 dirname = os.path.dirname(filepath)
diff --git a/abipy/wannier90/win.py b/abipy/wannier90/win.py
index 110f68d4e..f10969fe3 100644
--- a/abipy/wannier90/win.py
+++ b/abipy/wannier90/win.py
@@ -36,7 +36,7 @@ def structure2wannier90(structure, units="Bohr") -> str:
         fact = 1.0
         app("begin unit_cell_cart\nAng")
     else:
-        raise ValueError(f"Invalid {units =}")
+        raise ValueError(f"Invalid {units=}")
 
     for r in np.where(np.abs(structure.lattice.matrix) > 1e-8, structure.lattice.matrix, 0.0):
         r = r * fact
@@ -123,7 +123,7 @@ def __init__(self, structure, comment="", win_args=None, win_kwargs=None, spell_
         self._vars = OrderedDict(args)
 
     @property
-    def vars(self) ->dict:
+    def vars(self) -> dict:
         return self._vars
 
     # This stufff should be moved to the ABC