Add plot examples for gpath

abipy/eph/gpath.py

@@ -114,9 +114,6 @@ def to_string(self, verbose: int = 0) -> str:
         if self.r.eph_fix_korq == "q":
             app(self.ebands_k.to_string(with_structure=False, verbose=verbose, title="Electronic Bands (k)"))
 
-        #app(f"gstore_cplex: {self.r.cplex}")
-        #app(f"gstore_qptopt: {self.r.qptopt}")
-
         return "\n".join(lines)
 
     @staticmethod
@@ -134,16 +131,26 @@ def _get_band_range(self, band_range):
         return (self.r.bstart, self.r.bstop) if band_range is None else band_range
 
     @add_fig_kwargs
-    def plot_g_qpath(self, band_range=None, which_g="avg", with_qexp: int = 0, scale=1, gmax_mev=250,
-                     ph_modes=None, with_phbands=True, with_ebands=False,
-                     ax_mat=None, fontsize=8, **kwargs) -> Figure:
+    def plot_g_qpath(self,
+                     band_range=None,
+                     which_g="avg",
+                     with_qexp: int = 0,
+                     scale=1,
+                     gmax_mev=250,
+                     ph_modes=None,
+                     with_phbands=True,
+                     with_ebands=False,
+                     ax_mat=None,
+                     fontsize=8,
+                     **kwargs) -> Figure:
         """
-        Plot the averaged |g(k,q)| in meV units along the q-path
+        Plot the averaged |g(k,q)| in meV units along the q-path.
 
         Args:
             band_range: Band range that will be averaged over (python convention).
-            which_g: "avg" to plot the symmetrized |g|, "raw" for unsymmetrized |g|."all" for both.
-            with_qexp: Multiply |g(q)| by |q|^{with_qexp}.
+                If None all bands are considered.
+            which_g: "avg" to plot the symmetrized |g|, "raw" for unsymmetrized |g|. "all" for both.
+            with_qexp: Multiply |g(k, q)| by |q|^{with_qexp}.
             scale: Scaling factor for the marker size used when with_phbands is True.
             gmax_mev: Show results up to gmax in meV.
             ph_modes: List of ph branch indices to show (start from 0). If None all modes are shown.
@@ -221,8 +228,16 @@ def plot_g_qpath(self, band_range=None, which_g="avg", with_qexp: int = 0, scale
         return fig
 
     @add_fig_kwargs
-    def plot_g_kpath(self, band_range=None, which_g="avg", scale=1, gmax_mev=250, ph_modes=None,
-                    with_ebands=True, ax_mat=None, fontsize=8, **kwargs) -> Figure:
+    def plot_g_kpath(self,
+                     band_range=None,
+                     which_g="avg",
+                     scale=1,
+                     gmax_mev=250,
+                     ph_modes=None,
+                     with_ebands=True,
+                     ax_mat=None,
+                     fontsize=8,
+                     **kwargs) -> Figure:
         """
         Plot the averaged |g(k,q)| in meV units along the k-path
 
@@ -339,8 +354,9 @@ def __init__(self, filepath: PathLike):
         self.eph_fix_korq = self.read_string("eph_fix_korq")
         if self.eph_fix_korq not in {"k", "q"}:
             raise ValueError(f"Invalid value for {self.eph_fix_korq=}")
+
         self.eph_fix_wavec = self.read_value("eph_fix_wavevec")
-        self.dbdb_add_lr = self.read_value("dvdb_add_lr")
+        self.dvdb_add_lr = self.read_value("dvdb_add_lr")
         #self.used_ftinterp = self.read_value("used_ftinterp")
         #self.completed = self.read_value("gstore_completed")
 
@@ -621,10 +637,14 @@ class GpathRobot(Robot, RobotWithEbands):
     EXT = "GPATH"
 
     @add_fig_kwargs
-    def plot_g_qpath(self, which_g="avg", gmax_mev=250, ph_modes=None,
-                    colormap="jet", **kwargs) -> Figure:
+    def plot_g_qpath(self,
+                     which_g="avg",
+                     gmax_mev=250,
+                     ph_modes=None,
+                     colormap="jet",
+                     **kwargs) -> Figure:
         """
-        Compare the g-matrix along a q-path.
+        Compare the g-matrix stored in the Robot along a q-path.
 
         Args
             which_g: "avg" to plot the symmetrized |g|, "raw" for unsymmetrized |g|."all" for both.
@@ -677,6 +697,7 @@ def plot_g_qpath(self, which_g="avg", gmax_mev=250, ph_modes=None,
 
     #@add_fig_kwargs
     #def plot_g_kpath(self, **kwargs) --> Figure
+    #   """Compare the g-matrix stored in the Robot along a q-path."""
 
     def yield_figs(self, **kwargs):  # pragma: no cover
         """
@@ -703,4 +724,4 @@ def write_notebook(self, nbpath=None) -> str:
         #nb.cells.extend(self.get_baserobot_code_cells())
         #nb.cells.extend(self.get_ebands_code_cells())
 
-        return self._write_nb_nbpath(nb, nbpath)
+        return self._write_nb_nbpath(nb, nbpath)

abipy/eph/tests/test_gpath.py

@@ -0,0 +1,84 @@
+"""Tests for varpeq module."""
+import abipy.data as abidata
+
+from abipy.core.testing import AbipyTest
+from abipy.eph.gpath import GpathFile, GpathRobot
+
+
+class GpathTest(AbipyTest):
+
+    def test_gpath_file_fixed_k(self):
+        """Testing GPATH.nc file with fixed k."""
+        with GpathFile(abidata.ref_file("teph4zpr_9o_DS1_GPATH.nc")) as gpath:
+            # g(k,q) at fixed k along q-path.
+            repr(gpath)
+            str(gpath)
+            assert gpath.to_string(verbose=2)
+            assert not gpath.params
+            assert gpath.structure.formula == "Mg1 O1" and len(gpath.structure) == 2
+            assert gpath.r.eph_fix_korq == "k"
+            assert gpath.r.bstart == 5 and gpath.r.bstop == 8
+            assert gpath.r.nk_path == 1
+            assert gpath.r.nq_path == 42
+            assert len(gpath.ebands_k.kpoints) == 1
+            assert len(gpath.ebands_kq.kpoints) == 42
+            assert len(gpath.phbands.qpoints) == 42
+            self.assert_equal(gpath.r.eph_fix_wavec, (0, 0, 0))
+
+            if self.has_matplotlib():
+                assert gpath.plot_g_qpath(show=False)
+                with self.assertRaises(ValueError):
+                    gpath.plot_g_kpath(show=False)
+
+            # Test jupyter notebook creation
+            if self.has_nbformat():
+                gpath.write_notebook(nbpath=self.get_tmpname(text=True))
+
+    def test_gpath_file_fixed_q(self):
+        """Testing GPATH.nc file with fixed q."""
+        with GpathFile(abidata.ref_file("teph4zpr_9o_DS2_GPATH.nc")) as gpath:
+            # g(k,q) at fixed q along k-path.
+            repr(gpath)
+            str(gpath)
+            assert gpath.to_string(verbose=2)
+            assert not gpath.params
+            assert gpath.structure.formula == "Mg1 O1" and len(gpath.structure) == 2
+            assert gpath.r.eph_fix_korq == "q"
+            assert gpath.r.bstart == 5 and gpath.r.bstop == 8
+            assert gpath.r.nk_path == 42
+            assert gpath.r.nq_path == 1
+            assert len(gpath.ebands_k.kpoints) == 42
+            #assert len(gpath.ebands_kq.kpoints) == 1
+            assert len(gpath.phbands.qpoints) == 1
+            self.assert_equal(gpath.r.eph_fix_wavec, (0.11, 0, 0))
+
+            if self.has_matplotlib():
+                assert gpath.plot_g_kpath(show=False)
+                with self.assertRaises(ValueError):
+                    gpath.plot_g_qpath(show=False)
+
+            # Test jupyter notebook creation
+            if self.has_nbformat():
+                gpath.write_notebook(nbpath=self.get_tmpname(text=True))
+
+    def test_gpath_robot(self):
+        """Testing GpathRobot."""
+        files = abidata.ref_files(
+            "abinitio_qpath_V1QAVG.nc",
+            "interpolated_qpath_V1QAVG.nc",
+        )
+
+        with GpathRobot() as robot:
+            robot.add_file("one", abidata.ref_file("teph4zpr_9o_DS1_GPATH.nc"))
+            robot.add_file("two", abidata.ref_file("teph4zpr_9o_DS1_GPATH.nc"))
+            assert len(robot) == 2
+            repr(robot); str(robot)
+            robot.to_string(verbose=2)
+
+            # Test matplotlib methods
+            if self.has_matplotlib():
+                assert robot.plot_g_qpath(show=False)
+
+            # Test jupyter notebook creation
+            if self.has_nbformat():
+                robot.write_notebook(nbpath=self.get_tmpname(text=True))

abipy/examples/plot/plot_gpath.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+r"""
+e-ph matrix along a path in the BZ
+==================================
+
+This example shows how to use the GPATH.nc file produced by ABINIT
+to plot e-ph matrix along a path in the BZ.
+"""
+import os
+import abipy.data as abidata
+
+from abipy.eph.gpath import GpathFile
+
+#%%
+# Here we use one of the GPATH files shipped with abipy
+# with the e-ph matrix g(k,q) with q along a path and k fixed.
+# Replace the call to abidata.ref_file("") with the path to your GPATH.nc
+
+gpath_q = GpathFile(abidata.ref_file("teph4zpr_9o_DS1_GPATH.nc"))
+print(gpath_q)
+
+#%%
+# To plot the e-ph matrix elements as a function of q.
+gpath_q.plot_g_qpath(band_range=None,
+                     which_g="avg",
+                     with_qexp=0,
+                     scale=1,
+                     gmax_mev=250,
+                     ph_modes=None,
+                     with_phbands=True,
+                     with_ebands=False,
+                     )
+
+#%%
+# For the meaning of the different arguments, see the docstring
+print(gpath_q.plot_g_qpath.__doc__)
+
+#%%
+# Here we read another file
+# with the e-ph matrix g(k,q) with k along a path and q fixed.
+
+gpath_k = GpathFile(abidata.ref_file("teph4zpr_9o_DS2_GPATH.nc"))
+print(gpath_k)
+
+#%%
+# To plot the e-ph matrix elements as a function of q.
+gpath_k.plot_g_kpath(band_range=None,
+                     which_g="avg",
+                     scale=1,
+                     gmax_mev=250,
+                     ph_modes=None,
+                     with_ebands=True,
+                     )
+
+#%%
+# For the meaning of the different arguments, see the docstring
+print(gpath_k.plot_g_kpath.__doc__)

abipy/flowtk/abiphonopy.py

@@ -9,6 +9,7 @@
 from phonopy.interface.vasp import read_vasp_from_strings
 from phonopy.interface.abinit import parse_set_of_forces
 
+from abipy.abio.inputs import AbinitInput
 from abipy.core.structure import Structure
 from abipy.flowtk.works import Work
 
@@ -65,17 +66,19 @@ class PhonopyWork(Work):
     """
 
     @classmethod
-    def from_gs_input(cls, gs_inp, scdims, phonopy_kwargs=None, displ_kwargs=None) -> PhonopyWork:
+    def from_gs_input(cls,
+                      gs_inp: AbinitInput,
+                      scdims,
+                      phonopy_kwargs=None,
+                      displ_kwargs=None) -> PhonopyWork:
         """
-        Build the work from an :class:`AbinitInput` object representing a GS calculations.
+        Build the work from an AbinitInput object representing a GS calculations.
 
         Args:
-            gs_inp::class:`AbinitInput` object representing a GS calculation in the initial unit cell.
+            gs_inp: AbinitInput object representing a GS calculation in the initial unit cell.
             scdims: Number of unit cell replicas along the three reduced directions.
             phonopy_kwargs: (Optional) dictionary with arguments passed to Phonopy constructor.
             displ_kwargs: (Optional) dictionary with arguments passed to generate_displacements.
-
-        Return: PhonopyWork instance.
         """
         new = cls()
         new.phonopy_tasks, new.bec_tasks = [], []
@@ -117,13 +120,11 @@ def on_all_ok(self):
         phonon = self.phonon
 
         # Write POSCAR with initial unit cell.
-        structure = structure_from_atoms(phonon.get_primitive())
+        structure = structure_from_atoms(phonon.primitive)
         structure.to(filename=self.outdir.path_in("POSCAR"))
 
         # Write yaml file with displacements.
-        #supercell = phonon.get_supercell()
         supercell = phonon.supercell
-        #displacements = phonon.get_displacements()
         displacements = phonon.displacements
         file_IO.write_disp_yaml(displacements, supercell, # directions=directions,
                                 filename=self.outdir.path_in('disp.yaml'))
@@ -200,9 +201,14 @@ class PhonopyGruneisenWork(Work):
         numpy arrays with the number of cells in the supercell along the three reduced directions.
     """
     @classmethod
-    def from_gs_input(cls, gs_inp, voldelta, scdims, phonopy_kwargs=None, displ_kwargs=None) -> PhonopyGruneisenWork:
+    def from_gs_input(cls,
+                      gs_inp: AbinitInput,
+                      voldelta,
+                      scdims,
+                      phonopy_kwargs=None,
+                      displ_kwargs=None) -> PhonopyGruneisenWork:
         """
-        Build the work from an :class:`AbinitInput` object representing a GS calculations.
+        Build the work from an AbinitInput object representing a GS calculations.
 
         Args:
             gs_inp: :class:`AbinitInput` object representing a GS calculation in the initial unit cell.
@@ -211,8 +217,6 @@ def from_gs_input(cls, gs_inp, voldelta, scdims, phonopy_kwargs=None, displ_kwar
             scdims: Number of unit cell replicas along the three reduced directions.
             phonopy_kwargs: (Optional) dictionary with arguments passed to Phonopy constructor.
             displ_kwargs: (Optional) dictionary with arguments passed to generate_displacements.
-
-        Return: `PhonopyGruneisenWork` instance.
         """
         new = cls()
 
@@ -253,7 +257,7 @@ def on_all_ok(self):
         self.add_phonopy_works_and_build()
         return super().on_all_ok()
 
-    def add_phonopy_works_and_build(self):
+    def add_phonopy_works_and_build(self) -> None:
         """
         Get relaxed structures from the tasks, build Phonopy works with supercells
         constructed from the new structures, add them to the flow and build new directories.

abipy/ml/aseml.py

@@ -74,21 +74,12 @@ class RX_MODE(EnumMixin, StrEnum):  # StrEnum added in 3.11
     cell = "cell"
 
 
-def to_ase_atoms(structure: PmgStructure, calc=None) -> Atoms:
-    """Convert pymatgen structure to ASE atoms. Optionally, attach a calculator."""
-    structure = Structure.as_structure(structure)
-    atoms = AseAtomsAdaptor.get_atoms(structure)
-    if calc:
-        atoms.calc = calc
-    return atoms
-
-
 def get_atoms(obj: Any) -> Atoms:
     """Return ASE Atoms from object."""
     if isinstance(obj, str):
-        return to_ase_atoms(Structure.from_file(obj))
+        return Structure.from_file(obj).to_ase_atoms()
     if isinstance(obj, PmgStructure):
-        return to_ase_atoms(obj)
+        return Structure.as_structure(obj).to_ase_atoms()
     if isinstance(obj, Atoms):
         return obj
     raise TypeError(f"Don't know how to construct Atoms object from {type(obj)}")
@@ -100,7 +91,7 @@ def abisanitize_atoms(atoms: Atoms, **kwargs) -> Atoms:
     """
     structure = Structure.as_structure(atoms)
     new_structure = structure.abi_sanitize(**kwargs)
-    return to_ase_atoms(get_atoms(new_structure), calc=atoms.calc)
+    return new_structure.to_ase_atoms(calc=atoms.calc)
 
 
 def fix_atoms(atoms: Atoms,
@@ -3867,7 +3858,7 @@ def run_nn_name_set_name(self, nn_name: str, set_name: str) -> dict:
                     for volume in volumes:
                         #ase = structure.get_ase().copy()
                         #ase.set_cell(ase.get_cell() * float(scale_factor)**(1 / 3), scale_atoms=True)
-                        atoms = to_ase_atoms(Structure.as_structure(v0_atoms).scale_lattice(volume))
+                        atoms = Structure.as_structure(v0_atoms).scale_lattice(volume).to_ase_atoms()
                         r = AseResults.from_atoms(atoms, calc=calc)
                         energies.append(r.ene)
                         stresses.append(r.stress.tolist())

abipy/ml/ml_phonopy.py

@@ -13,6 +13,7 @@
 from monty.termcolor import cprint
 from ase.calculators.calculator import Calculator
 from ase.atoms import Atoms
+from ase.filters import FrechetCellFilter # ExpCellFilter,
 from pymatgen.io.phonopy import get_phonopy_structure
 from abipy.core.structure import Structure
 from abipy.dfpt.ddb import DdbFile
@@ -572,7 +573,6 @@ def run(self) -> None:
         atoms = self.initial_atoms.copy()
         atoms.calc = calculator
 
-        from ase.filters import ExpCellFilter
         relax_mode = RX_MODE.cell
         print(f"Relaxing initial atoms with {relax_mode=}.")
         relax_kws = dict(optimizer=self.optimizer,
@@ -598,7 +598,7 @@ def run(self) -> None:
             new_vol = v0 * bo_vol_scale
             new_atoms = new_structure.scale_lattice(new_vol).to_ase_atoms()
             new_atoms.calc = calculator
-            new_ucf = ExpCellFilter(new_atoms, constant_volume=True)
+            new_ucf = FrechetCellFilter(new_atoms, constant_volume=True)
 
             vol_workdir = self.get_path(f"VOL_{ivol}", info=f"Directory with phonopy results for {ivol=}, {new_vol=}")
             os.mkdir(vol_workdir)