Add run_qha_vzsisa.py flow script

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.

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):

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)

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]

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)
+        return self._write_nb_nbpath(nb, nbpath)

abipy/examples/flows/_under_develop/run_qha.py → 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,22 +19,32 @@ 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
     structure = abilab.Structure.from_file(abidata.cif_file("si.cif"))
     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.

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()

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`.
 

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