Add test_pade

abipy/dfpt/vzsisa.py

@@ -27,32 +27,34 @@ def anaget_phdoses_with_gauss(nqsmall_or_qppa,
                               smearing_ev: float | None,
                               ddb_paths: list[PathLike],
                               anaget_kwargs: dict | None,
-                              verbose: int) -> tuple(list[str]):
+                              verbose: int) -> tuple(list[str]) -> tuple[list[str], list[str]]:
     """
+    Invoke anaddb to compute PHDOSes from a list of DDB filepaths with the gaussian method.
+
+    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.
+
     Args:
         nqsmall_or_qppa: Define the q-mesh for the computation of the PHDOS.
             if > 0, it is interpreted as nqsmall
             if < 0, it is interpreted as qppa.
-        smearing_ev: Gaussian smearing in eV
+        smearing_ev: Gaussian smearing in eV.
         ddb_paths: List of paths to the DDB files.
-        anaget_kwargs: Dict with extra arguments passed to anaget_phbst_and_phdos_files
+        anaget_kwargs: dictionary with extra arguments passed to anaget_phbst_and_phdos_files.
         verbose: Verbosity level.
 
-    Returns:
+    Returns: tuple with the list of filepaths for phdos and phbands.
     """
-    # 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, phbands_paths = [], []
 
     if smearing_ev is None:
         smearing_ev = 4.5e-6 * abu.Ha_eV
 
-    my_kwargs = dict(dos_method=f"gaussian:{smearing_ev} eV",
-                    return_input=False)
+    my_kwargs = dict(dos_method=f"gaussian:{smearing_ev} eV", return_input=False)
 
     if nqsmall_or_qppa > 0:
         my_kwargs["nqsmall"] = nqsmall_or_qppa
@@ -64,9 +66,6 @@ def anaget_phdoses_with_gauss(nqsmall_or_qppa,
     if anaget_kwargs:
         my_kwargs.update(anaget_kwargs)
 
-    #if verbose:
-    #    print(f"Calling anaget_phbst_and_phdos_files with {my_kwargs=}:")
-
     from abipy.dfpt.ddb import DdbRobot
     with DdbRobot.from_files(ddb_paths) as robot:
         r = robot.anaget_phonon_plotters(**my_kwargs)
@@ -94,7 +93,9 @@ class Vzsisa(HasPickleIO):
 
     @classmethod
     def from_phonopy_files(cls, filepaths: list, bo_energies, phdos_kwargs=None):
-
+        """
+        Build an instance from a phonopy calculation.
+        """
         import phonopy
         bo_structures, ph_structures, phdoses = [], [], []
 
@@ -139,7 +140,7 @@ def from_json_file(cls,
                        smearing_ev: float | None = None,
                        verbose: int = 0) -> Vzsisa:
         """
-        Build an instance from a json file `filepath` typically produced by an Abipy flow.
+        Build an instance from a json file `filepath` typically produced by an AbiPy flow.
         For the meaning of the other arguments see from_gsr_ddb_paths.
         """
         data = mjson_load(filepath)
@@ -155,9 +156,9 @@ def from_gsr_ddb_paths(cls,
                            smearing_ev: float | None = None,
                            verbose: int = 0) -> Vzsisa:
         """
-        Creates an instance from a list of GSR files and a list of DDB files.
+        Create 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
+        from the DDB files by invoking anaddb.
 
         Args:
             nqsmall_or_qppa: Define the q-mesh for the computation of the PHDOS.
@@ -239,7 +240,7 @@ def __init__(self,
         """
         Args:
             bo_structures: list of structures at the different volumes for the BO bo_energies.
-            structure_from_phdos: Structures used to compute phonon DOS
+            ph_structures: Structures used to compute phonon DOSes.
             bo_energies: list of BO energies for the structures in eV.
             phdoses: Phonon DOSes.
             edoses: Electron DOSes. None if electronic contribution should not be considered.
@@ -364,7 +365,7 @@ def set_eos(self, eos_name: str) -> None:
     def get_gibbs_phvol(self, tstart, tstop, num) -> np.array:
         """
         Compute the Gibbs free energy in eV for all the ph_volumes and the list of temperatures
-        specified in input. Returnd array of shape [ph_nvol, num]
+        specified in input. Return array of shape [ph_nvol, num]
 
         Args:
             tstart: The starting value (in Kelvin) of the temperature mesh.
@@ -657,7 +658,7 @@ def get_phbands_plotter(self) -> PhononBandsPlotter:
     def get_edos_plotter(self) -> ElectronDosPlotter:
         """Build and return a ElectronDosPlotter with electron doses indexed by volume"""
         if self.edoses[0] is None:
-            raise ValueError("edoes_list is not available")
+            raise ValueError("edoses_list is not available")
         plotter = ElectronDosPlotter()
         for volume, edos in zip(self.ph_volumes, self.edoses, strict=True):
             plotter.add_edos(f"V={volume:.2f} ų", edos)
@@ -856,7 +857,7 @@ def plot_thermal_expansion_coeff(self, tstart=0, tstop=1000, num=101, tref=None,
         data = {"tmesh": tmesh}
         iv0 = self.iv0_vib
         iv1 = self.iv1_vib
-        dV = ph_volumes[iv0+1]-ph_volumes[iv0]
+        dV = ph_volumes[iv0+1] - ph_volumes[iv0]
 
         if self.scale_points == "S":
             vols, fits = self.vol_E2Vib1(num=num, tstop=tstop, tstart=tstart)
@@ -1799,7 +1800,7 @@ def plot_abc_vs_t_4th(self, tstart=0, tstop=1000, num=101, lattice=None, tref=No
         if lattice is None or lattice == "c":
             ax.plot(tmesh, cc, color='m', lw=2, label=r"$c(V(T))$" + method)
 
-        if abs(abs(self.bo_volumes[self.iv0] - ph_volumes[iv0])-abs(ph_volumes[iv1]-self.bo_volumes[self.iv0])) < 1e-3 :
+        if abs(abs(self.bo_volumes[self.iv0] - ph_volumes[iv0])-abs(ph_volumes[iv1]-self.bo_volumes[self.iv0])) < 1e-3:
             if lattice is None or lattice == "a":
                 ax.plot(tmesh, aa2, linestyle='dashed', color='r', lw=2, label=r"$a(V(T))$""E2vib1")
             if lattice is None or lattice == "b":

abipy/electrons/gwr.py

@@ -175,6 +175,7 @@ class PadeData:
 
 @dataclasses.dataclass(kw_only=True)
 class SigmaTauFit:
+    """Stores the fit for Sigma(i tau)"""
 
     tau_mesh: np.ndarray
     values: np.ndarray
@@ -220,6 +221,7 @@ def __init__(self, *args, **kwargs):
 
     def tau_fit(self, first, last, xs) -> tuple[np.ndarray, complex, float]:
         """
+        Performs the exponential fit in imaginary time.
         """
         mp_taus = self.c_tau.mesh
         vals_mptaus = self.c_tau.values

abipy/examples/flows/run_qha_vzsisa.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 r"""
-Flow for QHA calculations v-ZSISA-QHA
-=====================================
+Flow for v-ZSISA-QHA calculations
+================================
 
 See [Phys. Rev. B 110, 014103](https://doi.org/10.1103/PhysRevB.110.014103)
 """
@@ -30,22 +30,25 @@ def build_flow(options):
     from abipy.flowtk.psrepos import get_oncvpsp_pseudos
     pseudos = get_oncvpsp_pseudos(xc_name="PBEsol", version="0.4")
 
-    # Select k-mesh for electrons and q-mesh for phonons.
-    ngkpt = [2, 2, 2]; ngqpt = [1, 1, 1]
-    #ngkpt = [4, 4, 4]; ngqpt = [4, 4, 4]
+    # Select k-mesh for electrons and q-mesh for phonons. NB: ngqpt should be a divisor of ngkpt.
+    # These values are clearly UNDERCONVERGED, just for testing.
+    ngkpt = [2, 2, 2]
+    ngqpt = [1, 1, 1]
 
+    # BECS are not needed for Si. quadrupoles require pseudos without non-linear core correction.
     with_becs = False
     with_quad = False
     #with_quad = not structure.has_zero_dynamical_quadrupoles
 
+    # List of volumetric scaling factors for the BO energies and the phonon part.
     #bo_vol_scales = [0.96, 0.98, 1.0, 1.02, 1.04, 1.06]
     #ph_vol_scales = [0.98, 1.0, 1.02, 1.04, 1.06] # EinfVib4(D)
     bo_vol_scales = [0.96, 0.98, 1, 1.02, 1.04]    # EinfVib4(S)
     ph_vol_scales = [1, 1.02, 1.04]                # EinfVib2(D)
 
     scf_input = abilab.AbinitInput(structure, pseudos)
 
-    # Set other important variables
+    # Set other important variables assuming a spin unpolarized semiconductor (nsppol 1, nspinor 1).
     # All the other DFPT runs will inherit these parameters.
     scf_input.set_vars(
         nband=scf_input.num_valence_electrons // 2,
@@ -54,7 +57,7 @@ def build_flow(options):
         nstep=100,
         ecutsm=1.0,
         tolvrs=1.0e-8,   # SCF stopping criterion (modify default)
-        #tolvrs=1.0e-18,   # SCF stopping criterion (modify default)
+        #tolvrs=1.0e-18, # SCF stopping criterion (modify default)
     )
 
     scf_input.set_kmesh(ngkpt=ngkpt, shiftk=[0, 0, 0])

abipy/flowtk/vzsisa.py

@@ -83,7 +83,7 @@ def finalize(self):
         data = {"bo_vol_scales": work.bo_vol_scales, "ph_vol_scales": work.ph_vol_scales}
         data["initial_structure"] = work.initial_scf_input.structure
 
-        # Build list of strings with path to the relevant output files ordered by V.
+        # Build list of strings with path to the relevant output files ordered by volume.
         data["gsr_relax_paths"] = [task.gsr_path for task in work.relax_tasks_vol]
 
         gsr_relax_entries, gsr_relax_volumes = [], []
@@ -208,7 +208,7 @@ def on_all_ok(self):
             self.flow.register_work(ph_work)
             self.ph_works.append(ph_work)
 
-            # Add task for electron DOS calculation to edos_work
+            # Add task for electron DOS calculation to edos_work.
             if self.edos_work is not None:
                 edos_input = scf_input.make_edos_input(self.edos_ngkpt, prtwf=-1)
                 self.edos_work.register_nscf_task(edos_input, deps={ph_work[0]: "DEN"})

abipy/tools/pade.py

@@ -1,13 +1,14 @@
 # coding: utf-8
 """
 Functions to perform analytic continuation with Pade'
-Some of these routines have been directly translated from the Fortran versions
+Some of these routines have been directly translated from the Fortran version
 implemented in ABINIT.
 """
 from __future__ import annotations
 
 import numpy as np
 
+
 class SigmaPade:
     """
     High-level interface to perform the analytic continuation of the self-energy with the Pade' method.
@@ -24,7 +25,8 @@ def __init__(self, zs, f_zs):
 
     def eval(self, z_evals: np.ndarray) -> tuple[np.ndarray, np.np.ndarray]:
         """
-        Compute the Pade fit for a list of zz points
+        Compute the Pade fit for a list of zz points.
+        Return f(z_evals) and f'(z_evals)
         """
         nn = len(z_evals)
         sws = np.zeros(nn, dtype=complex)
@@ -37,8 +39,9 @@ def eval(self, z_evals: np.ndarray) -> tuple[np.ndarray, np.np.ndarray]:
 
         return sws, dsdws
 
-    def _eval_one(self, z_eval) -> tuple[np.ndarray, np.ndarray]:
+    def _eval_one(self, z_eval) -> tuple:
         """
+        Pade for a single point z_eval.
         """
         # if z_eval is in 2 or 3 quadrant, avoid the branch cut in the complex plane using Sigma(-iw) = Sigma(iw)*.
         # See also sigma_pade_eval in m_dyson_solver.F90
@@ -146,6 +149,7 @@ def dpade(zs: np.ndarray, f_zs: np.ndarray, z_eval: complex) -> complex:
     # print("Bz(n):", Bz[n])
     # if np.isclose(Bz[n].real, 0.0) and np.isclose(Bz[n].imag, 0.0):
     #     print("Warning: Bz(n) is close to zero:",
+    return dpade_value
 
 
 def calculate_pade_a(zs: np.ndarray, f_zs: np.ndarray) -> np.ndarray:

abipy/tools/parallel.py

@@ -1,5 +1,5 @@
 """
-Tools used to parallelize sections of python code.
+Tools used to parallelize sections of python code with multiprocessing or threads.
 """
 from __future__ import annotations
 
@@ -45,7 +45,6 @@ def pool_nprocs_pmode(nprocs: int | None, pmode: str):
     """
     from multiprocessing.pool import ThreadPool
     from multiprocessing import Pool
-
     max_nprocs = get_max_nprocs()
 
     if pmode == "seq":

abipy/tools/tests/test_pade.py

@@ -0,0 +1,34 @@
+"""Tests for frozen_phonons"""
+import os
+import numpy as np
+#import abipy.data as abidata
+
+from abipy.core.testing import AbipyTest
+from abipy.tools.pade import SigmaPade
+
+
+class PadeTest(AbipyTest):
+
+    def test_sigma_pade(self):
+        """Testing SigmaPade interface."""
+
+        # List of points and values for the Pade'
+        zs = np.array([0.3204639033968435j, 1.159257347625358j, 2.7851769489568112j,
+                       6.6074592214004655j, 16.793387237006716j, 50.739714260170196j],
+                       dtype=complex)
+
+        f_zs = np.array([(-3.4591316381250774-0.09890936491537991j), (-3.4538645067311458-0.24793469903457657j),
+                        (-3.4173978321863663-0.6621397224074206j), (-3.257949585597619-1.324699039193591j),
+                        (-2.523177170083569-2.4583423787866163j), (-0.8792031880464662-2.2566436943461583j)],
+                        dtype=complex)
+
+        pade = SigmaPade(zs, f_zs)
+
+        # Evaluate Pade' as w_vals.
+        w_vals = np.array([-0.49705069294512505, -0.397050692945125, -0.297050692945125, -0.197050692945125])
+        this_f, this_fp = pade.eval(w_vals)
+
+        ref_f = [-3.36245542+0.0371258j, -3.38196076+0.03677179j, -3.40149899+0.036502j, -3.42108152+0.03632633j]
+        self.assert_almost_equal(this_f, ref_f, decimal=7)
+        ref_fp = [-0.19492091-0.00392504j, -0.19520083-0.00313668j, -0.19558294-0.00224237j, -0.1960902-0.00125779j]
+        self.assert_almost_equal(this_fp, ref_fp, decimal=7)