Add structure.as_ngkpt method

abipy/core/structure.py

@@ -172,6 +172,7 @@ def display_structure(obj, **kwargs):
 
 def get_structures_from_file(filepath: PathLike, index) -> list[Structure]:
     """
+    Read and return list of structures from filepath
     """
     #if index is None:
     #    index = -1
@@ -2409,16 +2410,31 @@ 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)
+    def as_ngkpt(self, ngkpt) -> np.ndarray:
+        """
+        Flexible API to compute the ABINIT variable ``ngkpt`` using different approaches.
+
+        The following cases are supported:
+            - If `ngkpt` is a 1D vector with 3 items, return `ngkpt` as-is.
+            - If `ngkpt` is a positive float, interpret it as `nksmall` (a scaling factor for the k-point grid).
+            - If `ngkpt` is a negative float, interpret it as the desired number of k-points per atom.
+
+        This method allows users to flexibly specify the k-point grid based on their preferred input format.
+        """
+        ngkpt = np.array(ngkpt)
+
+        if ngkpt.ndim == 1 and len(ngkpt) == 3:
+            return ngkpt
+
+        if (nksmall := float(ngkpt)) > 0:
+            return self.calc_ngkpt(nksmall)
+
+        if (kppa := -float(ngkpt)) > 0:
+            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"]
 
-    #    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"]
+        raise ValueError(f"Don't know how to convert {type(ngkpt)=}, {ngkpt=} to k-mesh!")
 
     def calc_shiftk(self, symprec=0.01, angle_tolerance=5) -> np.ndarray:
         """

abipy/core/tests/test_structure.py

@@ -159,6 +159,7 @@ def test_utils(self):
         assert len(si.abi_string)
         assert si.reciprocal_lattice == si.lattice.reciprocal_lattice
 
+        # Test k-sampling methods.
         kptbounds = si.calc_kptbounds()
         ksamp = si.calc_ksampling(nksmall=10)
 
@@ -168,6 +169,13 @@ def test_utils(self):
         self.assert_equal(ksamp.ngkpt, [10, 10, 10])
         self.assert_equal(ksamp.shiftk, shiftk)
 
+        # Test as_ngkpt API
+        self.assert_equal(si.as_ngkpt([1, 2, 3]), [1, 2, 3])
+        self.assert_equal(si.as_ngkpt(3), [3, 3, 3])
+        self.assert_equal(si.as_ngkpt(-1000), [8, 8, 8])
+        with self.assertRaises(ValueError):
+            si.as_ngkpt("foo")
+
         lif = Structure.from_abistring("""
 acell      7.7030079150    7.7030079150    7.7030079150 Angstrom
 rprim      0.0000000000    0.5000000000    0.5000000000

abipy/examples/flows/run_conducwork.py

@@ -15,8 +15,7 @@
 from abipy.abio.factories import conduc_kerange_from_inputs
 
 
-def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0),
-                   **variables):
+def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), **variables):
     """Build and return SCF input given the structure and pseudopotentials"""
 
     scf_inp = abilab.AbinitInput(structure, pseudos=pseudos)
@@ -31,8 +30,7 @@ def make_scf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0),
     return scf_inp
 
 
-def make_nscf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0),
-                    **variables):
+def make_nscf_input(structure, pseudos, ngkpt=(2,2,2), shiftk=(0,0,0), **variables):
     """Build and return NSCF input given the structure and pseudopotentials"""
     scf_inp = abilab.AbinitInput(structure, pseudos=pseudos)
 
@@ -76,14 +74,14 @@ def build_flow(options):
 
     # Kerange Variables
     nbr_proc = 4
-    ngqpt_fine = [16, 16, 16] # The sigma_ngkpt grid must be divisible by the qpt grid
+    ngqpt_fine = [16, 16, 16]         # The sigma_ngkpt grid must be divisible by the qpt grid
     sigma_ngkpt = [16, 16, 16]
-    einterp = [1, 5, 0, 0] # Star functions Interpolation
+    einterp = [1, 5, 0, 0]            # Star functions Interpolation
     sigma_erange = [-0.3, -0.3, "eV"] # Negative value for metals
 
     flow = flowtk.Flow(workdir=options.workdir)
 
-    # Create inputs Object
+    # Create inputs
     scf_input = make_scf_input(structure, pseudos,
                                tolvrs=1e-12,
                                ngkpt=ngkpt,
@@ -118,17 +116,18 @@ def build_flow(options):
                                        sigma_erange=sigma_erange,
                                        einterp=einterp,
                                        boxcutmin=boxcutmin, # 1.1 is the default value of the function
-                                       mixprec=mixprec # 1 is the default value of the function
+                                       mixprec=mixprec      # 1 is the default value of the function
                                        )
 
     # Here we can change multi to change the variable of a particular dataset
 
-    conduc_work = flowtk.ConducWork.from_phwork(phwork=ph_work, # Linking the DDB and DVDB via a |PhononWork|
-                                                multi=multi, # The multidataset object
+    conduc_work = flowtk.ConducWork.from_phwork(phwork=ph_work,    # Linking the DDB and DVDB via a PhononWork
+                                                multi=multi,       # The multidataset object
                                                 nbr_proc=nbr_proc, # Needed to parallelize the calculation
                                                 flow=flow,
                                                 with_kerange=True, # Using Kerange
-                                                omp_nbr_thread=1) # 1 is the default value of the function
+                                                omp_nbr_thread=1)  # 1 is the default value of the function
+
     # If you already have the DDB and DVDB, use from_filepath(ddb_path, dvdb_path, multi, ...) instead of from_phwork
 
     flow.register_work(conduc_work)
@@ -157,5 +156,3 @@ def main(options):
 
 if __name__ == "__main__":
     sys.exit(main())
-
-############################################

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.qha import VzsisaFlow
 
 
 def build_flow(options):
@@ -26,24 +26,25 @@ def build_flow(options):
     structure = abilab.Structure.from_file(abidata.cif_file("si.cif"))
     pseudos = abidata.pseudos("14si.pspnc")
 
-    # Build input for GS calculation.
+    # Select k-mesh for electrons and q-mesh for phonons.
     #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)
+    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=8, nband=4, tolvrs=1e-8, nstep=50) #, paral_kgb=0)
+    scf_input.set_vars(ecut=8, nband=4, tolvrs=1e-8, nstep=50)
     scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0])
 
-    return vZSISAFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt,
+    return VzsisaFlow.from_scf_input(options.workdir, scf_input, bo_scales, ph_scales, ngqpt,
                                      with_becs, with_quad, edos_ngkpt=None)
 
 

abipy/flowtk/qha.py

@@ -13,25 +13,25 @@
 from abipy.flowtk.flows import Flow
 
 
-class vZSISAFlow(Flow):
+class VzsisaFlow(Flow):
     """
-    Flow for QHA calculations with the vZSISA approach
+    Flow for QHA calculations with the VZSISA approach.
 
     .. rubric:: Inheritance Diagram
-    .. inheritance-diagram:: vZSISAFlow
+    .. 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:
+                       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:
+            bo_scales: List of scaling factors for the BO volumes.
+            ph_scales: List of scaling factors for the BO volumes.
             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`
@@ -45,8 +45,9 @@ def from_scf_input(cls, workdir, scf_input, bo_scales, ph_scales, ngqpt, with_be
         """
         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)
+        # 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)
 
         flow.register_work(work)
         return flow
@@ -61,47 +62,45 @@ def finalize(self):
 
         # 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]
+
+        entries, gsr_relax_volumes = [], []
+        for task in work.relax_tasks_vol:
+            with 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_relax_volumes.append(gsr.structure.volume)
+
+        data["gsr_entries"] = entries
+        data["gsr_relax_volumes"] = gsr_relax_volumes
+
         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
+        data["ddb_relax_volumes"] = [ph_work[0].input.structure.volume for ph_work in work.ph_works]
+
+        data["gsr_edos_path"] = [] if not work.edos_work else [task.gsr_path for task in work.edos_work]
 
         mjson_write(data, self.outdir.path_in("vzsisa.json"), indent=4)
 
         return super().finalize()
 
 
-class vZSISAWork(Work):
+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:: vZSISAWork
+    .. inheritance-diagram:: VzsisaWork
     """
 
     @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:
+                       optcell: int, ionmov: int, edos_ngkpt=None) -> VzsisaWork:
         """
         Build the work from an |AbinitInput| representing a GS-SCF calculation.
 
@@ -172,10 +171,16 @@ def on_all_ok(self):
             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, with_quad=self.with_quad,
                                                 ddk_tolerance=None)
 
+            # Reduce the number of files produced in the DFPT tasks to avoid possible disk quota issues.
+            prtvars = dict(prtwf=-1, prtden=0, prtpot=0)
+            for task in ph_work[1:]:
+                task.input.set_vars(**prtvars)
+
             self.flow.register_work(ph_work)
             self.ph_works.append(ph_work)
 

abipy/flowtk/tests/test_qha.py

@@ -1,19 +1,27 @@
 """Tests for qha module."""
 
 from abipy.core.testing import AbipyTest
-from abipy.flowtk.qha import QhaFlow
+from abipy.flowtk.qha import VzsisaFlow
 
 
-class TestQha(AbipyTest):
+class TestVzsisa(AbipyTest):
 
-    def test_qhaflow(self):
-        """Testing QhaFlow."""
+    def test_vzsisa_flow(self):
+        """Testing VzsisaFlow."""
         workdir = self.mkdtemp()
-        scf_input = self.get_gsinput_alas_ngkpt(ngkpt=[4, 4, 4])
-        #v0 = scf_input.structure.volume
-        #volumes = [0.98 * v0, v0, v0 * 1.02]
-        flow = QhaFlow.from_scf_input(workdir, scf_input, ngqpt=[2, 2, 2],
-                                      with_becs=False, edos_ngkpt=[4, 4, 4]) #, metadata={"mpi_id": 123})
+        ngkpt = [4, 4, 4]
+        ngqpt = [2, 2, 2]
+        scf_input = self.get_gsinput_alas_ngkpt(ngkpt=ngkpt)
+
+        bo_scales = [0.96, 0.98, 1, 1.02, 1.04]    # EinfVib4(S)
+        ph_scales = [1, 1.02, 1.04]                # EinfVib2(D)
+
+        with_becs = False
+        with_quad = False
+        #with_quad = not structure.has_zero_dynamical_quadrupoles
+
+        flow = VzsisaFlow.from_scf_input(workdir, scf_input, bo_scales, ph_scales, ngqpt,
+                                         with_becs, with_quad, edos_ngkpt=None)
 
         flow.allocate()
         flow.check_status()