diff --git a/abipy/dfpt/deformation_utils.py b/abipy/dfpt/deformation_utils.py
index fb6a56a23..37e862e92 100644
--- a/abipy/dfpt/deformation_utils.py
+++ b/abipy/dfpt/deformation_utils.py
@@ -1,11 +1,14 @@
 # deformation_utils.py
+from __future__ import annotations
 
 import numpy as np
-from pymatgen.core import Structure, Lattice, Element
+import re
+
+from pymatgen.core import Lattice, Element
+#from abipy.core.structure import Structure
 from abipy.core.symmetries import AbinitSpaceGroup
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from abipy.abio.inputs import AbinitInput
-import re
 
 
 def generate_deformations_volumic(structure, eps_V=0.02, scales=None):
diff --git a/abipy/dfpt/tests/test_qha_approximation.py b/abipy/dfpt/tests/test_vzsisa.py
similarity index 100%
rename from abipy/dfpt/tests/test_qha_approximation.py
rename to abipy/dfpt/tests/test_vzsisa.py
diff --git a/abipy/dfpt/qha_aproximation.py b/abipy/dfpt/vzsisa.py
similarity index 98%
rename from abipy/dfpt/qha_aproximation.py
rename to abipy/dfpt/vzsisa.py
index cf482f709..3243cfb87 100644
--- a/abipy/dfpt/qha_aproximation.py
+++ b/abipy/dfpt/vzsisa.py
@@ -1,6 +1,6 @@
 # coding: utf-8
 """
-Classes for V-ZSISA approximation.
+Classes  for post-processing QHA results obtained with the V-ZSISA approximation.
 See PHYSICAL REVIEW B 110, 014103 (2024)
 """
 from __future__ import annotations
@@ -19,7 +19,7 @@
 from abipy.dfpt.phonons import PhdosFile, PhononDosPlotter
 
 
-class QHA_App:
+class Vzsisa:
     """
     Class for the approximations on QHA analysis.
     Provides some basic methods and plotting utils.
@@ -31,7 +31,7 @@ def from_gsr_ddb_files(cls,
                            gsr_paths: list[PathLike],
                            ddb_paths: list[PathLike],
                            verbose: int = 0,
-                           **extra_kwargs) -> QHA_App:
+                           **extra_kwargs) -> Vzsisa:
         """
         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
@@ -71,7 +71,7 @@ def from_gsr_ddb_files(cls,
         return cls.from_gsr_phdos_files(gsr_paths, phdos_paths)
 
     @classmethod
-    def from_gsr_phdos_files(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]) -> Vzsisa:
         """
         Creates an instance from a list of GSR files and a list of PHDOS.nc files.
 
@@ -112,11 +112,8 @@ def from_gsr_phdos_files(cls, gsr_paths: list[str], phdos_paths: list[str]) -> Q
 
         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_ddb_phdos_files(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]) -> Vzsisa:
         """
         Creates an instance from a list of DDB files and a list of PHDOS.nc files.
 
@@ -158,9 +155,6 @@ def from_ddb_phdos_files(cls, ddb_paths: list[str], phdos_paths: list[str]) -> Q
 
         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):
         """
@@ -529,6 +523,8 @@ def plot_energies(self, tstart=0, tstop=1000, num=1, ax=None, fontsize=10, **kwa
         ax.plot(f.min_vol, f.min_en - self.energies[self.iv0], color='r', linestyle='dashed' , lw=1, marker='o', ms=5)
         set_grid_legend(ax, fontsize, xlabel=r'V (${\AA}^3$)', ylabel='E (eV)', legend=False)
 
+        fig.suptitle("Energies as a function of volume for different T")
+
         return fig
 
     def get_phdos_plotter(self) -> PhononDosPlotter:
@@ -605,6 +601,8 @@ def plot_vol_vs_t(self, tstart=0, tstop=1000, num=101, fontsize=10, ax=None, **k
         set_grid_legend(ax, fontsize, xlabel='T (K)', ylabel=r'V (${\AA}^3$)')
         ax.set_xlim(tstart, tstop)
 
+        fig.suptitle("Volume as a function of T")
+
         return fig
 
     def get_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None) -> Function1D:
@@ -775,6 +773,8 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Volumetric thermal expansion coefficient as a function of T")
+
         return fig
 
     def get_abc(self, volumes, num=101) -> tuple:
@@ -924,6 +924,8 @@ def plot_thermal_expansion_coeff_abc(self, tstart=0, tstop=1000, num=101, tref=N
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Thermal expansion coefficient as a function of T")
+
         return fig
 
     @add_fig_kwargs
@@ -1107,6 +1109,8 @@ def plot_abc_vs_t(self, tstart=0, tstop=1000, num=101, lattice=None, tref=None,
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Lattice as a function of T")
+
         return fig
 
     @add_fig_kwargs
@@ -1181,6 +1185,8 @@ def plot_angles_vs_t(self, tstart=0, tstop=1000, num=101, angle=None, tref=None,
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Angles as a function of T")
+
         return fig
 
     def fit_forth(self, tstart=0, tstop=1000, num=1, energy="energy", volumes="volumes"):
@@ -1434,6 +1440,8 @@ def plot_vol_vs_t_4th(self, tstart=0, tstop=1000, num=101, ax=None, fontsize=10,
         set_grid_legend(ax, fontsize, xlabel='T (K)', ylabel=r'V (${\AA}^3$)')
         ax.set_xlim(tstart, tstop)
 
+        fig.suptitle("Volume as a function of T")
+
         return fig
 
     def get_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=None) -> np.ndarray:
@@ -1558,7 +1566,7 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
             ds_dv = ds_dv + (df_t[:,iv0+2]-2*df_t[:,iv0+1]+df_t[:,iv0])/dV**2 * (vol3_4th[:]-volumes[iv0+1])
             if tref is None:
                 alpha_3 = - 1/vol3_4th[:] * ds_dv / (E2D_V[:]+dfe_dV2[:])
-            else :
+            else:
                 vol3_4th_ref = self.vol_Einf_Vib2_forth(num=1, tstop=tref, tstart=tref)
                 alpha_3 = - 1/vol3_4th_ref * ds_dv / (E2D_V[:]+dfe_dV2[:])
 
@@ -1602,6 +1610,8 @@ def plot_thermal_expansion_coeff_4th(self, tstart=0, tstop=1000, num=101, tref=N
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Volumetric thermal expansion coefficient as a function of T")
+
         return fig
 
     @add_fig_kwargs
@@ -1676,6 +1686,8 @@ def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None, tref=No
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Lattice as a function of T")
+
         return fig
 
     @add_fig_kwargs
@@ -1750,6 +1762,8 @@ def plot_angles_vs_t_4th(self, tstart=0, tstop=1000, num=101, angle=None, tref=N
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Angles as a function of T")
+
         return fig
 
     @add_fig_kwargs
@@ -1858,6 +1872,8 @@ def plot_thermal_expansion_coeff_abc_4th(self, tstart=0, tstop=1000, num=101, tr
         ax.set_xlim(tstart, tstop)
         ax.get_yaxis().get_major_formatter().set_powerlimits((0, 0))
 
+        fig.suptitle("Thermal expansion coefficient as a function of T")
+
         return fig
 
     @add_fig_kwargs
diff --git a/abipy/examples/flows/run_qha_vzsisa.py b/abipy/examples/flows/run_qha_vzsisa.py
index 0dc6b98b0..b0442e209 100755
--- a/abipy/examples/flows/run_qha_vzsisa.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 VzsisaFlow
+from abipy.flowtk.vzsisa import VzsisaFlow
 
 
 def build_flow(options):
diff --git a/abipy/examples/plot/plot_qha_v-ZSISA.py b/abipy/examples/plot/plot_qha_vzsisa.py
similarity index 51%
rename from abipy/examples/plot/plot_qha_v-ZSISA.py
rename to abipy/examples/plot/plot_qha_vzsisa.py
index 84a7306e7..bc27d85fa 100755
--- a/abipy/examples/plot/plot_qha_v-ZSISA.py
+++ b/abipy/examples/plot/plot_qha_vzsisa.py
@@ -9,7 +9,7 @@
 import os
 import abipy.data as abidata
 
-from abipy.dfpt.qha_aproximation import QHA_App
+from abipy.dfpt.vzsisa import Vzsisa
 
 # 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")
@@ -23,64 +23,64 @@
 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]
 
-qha = QHA_App.from_files_app_ddb(ddb_paths, dos_paths)
+qha = Vzsisa.from_files_app_ddb(ddb_paths, dos_paths)
 tstart, tstop = 0, 800
 
 #%%
-qha.plot_energies(tstop=tstop, tstart=tstart, num=11,
-                  title="Energies as a function of volume for different T")
+# Plot Energies as a function of volume for different T
+qha.plot_energies(tstop=tstop, tstart=tstart, num=11)
 
 #%%
-qha.plot_vol_vs_t(tstop=tstop, tstart=tstart, num=101,
-                  title="Volume as a function of T")
+# Plot Volume as a function of T
+qha.plot_vol_vs_t(tstop=tstop, tstart=tstart, num=101)
 
 #%%
-qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101,
-                  title="Lattice as a function of T")
-
+# Plot Lattice as a function of T
+qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101)
 
 #%%
-qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, lattice="b",
-                  title="Lattice as a function of T")
+# Plot Lattice as a function of T")
+qha.plot_abc_vs_t(tstop=tstop, tstart=tstart, num=101, lattice="b")
 
 #%%
-qha.plot_thermal_expansion_coeff(tstop=tstop, tstart=tstart ,num=101,
-                                 title="Volumetric thermal expansion coefficient as a function of T")
+# Plot Volumetric thermal expansion coefficient as a function of T
+qha.plot_thermal_expansion_coeff(tstop=tstop, tstart=tstart ,num=101)
 
 #%%
-qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart ,num=101,
-                                     title="Thermal expansion coefficient as a function of T")
+# Plot Thermal expansion coefficient as a function of T
+qha.plot_thermal_expansion_coeff_abc(tstop=tstop, tstart=tstart ,num=101)
 
 #%%
-qha.plot_angles_vs_t(tstop=tstop, tstart=tstart, num=101,
-                     title="Angles as a function of T")
+# Plot Angles as a function of T
+qha.plot_angles_vs_t(tstop=tstop, tstart=tstart, num=101)
 
 #%%
-# 4th order polinomial
-qha.plot_vol_vs_t_4th(tstop=tstop, tstart=tstart, num=101,
-                      title="Volume as a function of T")
+#
+# Plot Volume as a function of T. 4th order polinomial
+qha.plot_vol_vs_t_4th(tstop=tstop, tstart=tstart, num=101)
 
 #%%
-qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, num=101, lattice="a",
-                      title="Lattice as a function of T")
+# Plot Lattice as a function of T. 4th order polinomial
+qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart, num=101, lattice="a")
 
 #%%
-qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart,
-                      title="Lattice as a function of T")
+# Plot Lattice as a function of T. 4th order polinomial
+qha.plot_abc_vs_t_4th(tstop=tstop, tstart=tstart)
 
 #%%
-qha.plot_thermal_expansion_coeff_4th(tref=293,
-                                     title="Volumetric thermal expansion coefficient as a function of T")
+# Plot Volumetric thermal expansion coefficient as a function of T
+qha.plot_thermal_expansion_coeff_4th(tref=293)
 
 #%%
-qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293,
-                                         title="Thermal expansion coefficient as a function of T")
+# Plot Thermal expansion coefficient as a function of T
+qha.plot_thermal_expansion_coeff_abc_4th(tstop=tstop, tstart=tstart ,num=101, tref=293)
 
 #%%
-qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3,
-                         title="Angles as a function of T")
+# Plot Angles as a function of T.
+qha.plot_angles_vs_t_4th(tstop=tstop, tstart=tstart, num=101, angle=3)
 
 
 #%%
+# Create plotter to plot all the phonon DOS.
 phdos_plotter = qha.get_phdos_plotter()
 phdos_plotter.combiplot()
diff --git a/abipy/flowtk/psrepos.py b/abipy/flowtk/psrepos.py
index 268a6d614..5953929f0 100644
--- a/abipy/flowtk/psrepos.py
+++ b/abipy/flowtk/psrepos.py
@@ -56,7 +56,13 @@ def get_oncvpsp_pseudos(xc_name: str, version: str, relativity_type: str = "SR",
     ps_generator, project_name = OncvpspRepo.ps_generator, OncvpspRepo.project_name
     repo_name = f"{ps_generator}-{xc_name}-{relativity_type}-{project_name}v{version}"
 
-    return get_repo_from_name(repo_name).get_pseudos(accuracy)
+    repo = get_repo_from_name(repo_name)
+
+    # Istall it at runtime if needed.
+    if not repo.is_installed():
+        repo.install()
+
+    return repo.get_pseudos(accuracy)
 
 
 # Installation directory.
diff --git a/abipy/flowtk/tests/test_qha.py b/abipy/flowtk/tests/test_vzsisa.py
similarity index 100%
rename from abipy/flowtk/tests/test_qha.py
rename to abipy/flowtk/tests/test_vzsisa.py
diff --git a/abipy/flowtk/qha.py b/abipy/flowtk/vzsisa.py
similarity index 83%
rename from abipy/flowtk/qha.py
rename to abipy/flowtk/vzsisa.py
index 0b773e90c..28c28ee4e 100644
--- a/abipy/flowtk/qha.py
+++ b/abipy/flowtk/vzsisa.py
@@ -7,7 +7,6 @@
 import numpy as np
 
 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
@@ -45,9 +44,8 @@ def from_scf_input(cls, workdir, scf_input, bo_scales, ph_scales, ngqpt, with_be
         """
         flow = cls(workdir=workdir, manager=manager)
 
-        # optcell = 3: constant-volume optimization of cell geometry
         work = VzsisaWork.from_scf_input(scf_input, bo_scales, ph_scales, ngqpt, with_becs, with_quad,
-                                         optcell=3, ionmov=2, edos_ngkpt=edos_ngkpt)
+                                         ionmov=2, edos_ngkpt=edos_ngkpt)
 
         flow.register_work(work)
         return flow
@@ -100,7 +98,7 @@ class VzsisaWork(Work):
     @classmethod
     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:
+                       ionmov: int, edos_ngkpt=None) -> VzsisaWork:
         """
         Build the work from an |AbinitInput| representing a GS-SCF calculation.
 
@@ -109,7 +107,6 @@ def from_scf_input(cls, scf_input, bo_scales, ph_scales, ngqpt,
             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
                 electron DOS with the relaxed structures. Useful for metals
                 in which the electronic contribution should be included.
@@ -129,9 +126,9 @@ def from_scf_input(cls, scf_input, bo_scales, ph_scales, ngqpt,
         # Create input for relaxation and register the relaxation task.
         work.relax_template = relax_template = scf_input.deepcopy()
         relax_template.pop_tolerances()
-        relax_template.set_vars(tolvrs=1e-8, toldff=1.e-6, optcell=optcell, ionmov=ionmov)
-        if optcell is not None and optcell != 0:
-            relax_template.set_vars_ifnotin(ecutsm=0.5, dilatmx=1.05)
+        # optcell = 3: constant-volume optimization of cell geometry
+        relax_template.set_vars(optcell=3, ionmov=ionmov, tolvrs=1e-8, toldff=1.e-6)
+        relax_template.set_vars_ifnotin(ecutsm=0.5, dilatmx=1.05)
 
         work.initial_relax_task = work.register_relax_task(relax_template)
 
@@ -193,47 +190,3 @@ def on_all_ok(self):
         self.flow.allocate(build=True)
 
         return super().on_all_ok()
-
-
-class ThermalRelaxTask(RelaxTask):
-
-    def _on_ok(self):
-        results = super()._on_ok()
-        # Check for convergence.
-        #if not self.collinear_done:
-        #    self.input.set_vars(strtarget=strtarget)
-        #    self.finalized = False
-        #    self.restart()
-
-        return results
-
-
-class ThermalRelaxWork(Work):
-    """
-    .. rubric:: Inheritance Diagram
-    .. inheritance-diagram:: ThermalRelaxWork
-    """
-
-    @classmethod
-    def from_relax_input(cls, relax_input, qha, temperatures, pressures):
-        """
-        """
-        work = cls()
-
-        work.temperatures = temperatures
-        work.pressures = pressures
-        #work.qha = qha
-
-        for pressure in pressures:
-            for temperature in temperatures:
-                strtarget = qha.get_strtarget(temperature, pressure)
-                new_input = relax_input.new_with_vars(strtarget=strtarget)
-                work.register_relax_task(new_input)
-
-        return work
-
-    def on_all_ok(self):
-        """
-        Implement the post-processing step at the end of the Work.
-        """
-        return super().on_all_ok()
diff --git a/abipy/tools/plotting.py b/abipy/tools/plotting.py
index ae564a6d4..32f8844e0 100644
--- a/abipy/tools/plotting.py
+++ b/abipy/tools/plotting.py
@@ -67,8 +67,7 @@
 
 
 def add_fig_kwargs(func):
-    """Decorator that adds keyword arguments for functions returning matplotlib
-    figures.
+    """Decorator that adds keyword arguments for functions returning matplotlib figures.
 
     The function should return either a matplotlib figure or None to signal
     some sort of error/unexpected event.
diff --git a/docs/api/dfpt_api.rst b/docs/api/dfpt_api.rst
index db2916d8d..8d06b7ab7 100644
--- a/docs/api/dfpt_api.rst
+++ b/docs/api/dfpt_api.rst
@@ -101,6 +101,15 @@ dfpt Package
    :undoc-members:
    :show-inheritance:
 
+:mod:`vzsisa` Module
+--------------------
+
+.. automodule:: abipy.dfpt.vzsisa
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 :mod:`qha_2D` Module
 --------------------
 
diff --git a/docs/api/flowtk_api.rst b/docs/api/flowtk_api.rst
index dc95a2642..ee5d4585c 100644
--- a/docs/api/flowtk_api.rst
+++ b/docs/api/flowtk_api.rst
@@ -173,10 +173,18 @@ flowtk Package
    :undoc-members:
    :show-inheritance:
 
-:mod:`qha` Module
------------------
+:mod:`vzsisa` Module
+--------------------
+
+.. automodule:: abipy.flowtk.vzsisa
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+:mod:`zsisa` Module
+-------------------
 
-.. automodule:: abipy.flowtk.qha
+.. automodule:: abipy.flowtk.zsisa
    :members:
    :undoc-members:
    :show-inheritance: