Improve notebooks

abipy/dfpt/ddb.py

@@ -704,20 +704,41 @@ class DdbFile(TextFile, Has_Structure, NotebookWriter):
 
         nb.cells.extend([
             nbv.new_code_cell("ddb = abilab.abiopen('%s')" % self.filepath),
-            nbv.new_code_cell("print(ddb)"),
-            nbv.new_code_cell("display(ddb.header)"),
+            nbv.new_code_cell("units = 'eV'\nprint(ddb)"),
+            nbv.new_code_cell("# display(ddb.header)"),
+            nbv.new_markdown_cell("## Invoke `anaddb` to compute bands and dos"),
             nbv.new_code_cell("""\
-bstfile, phdosfile =  ddb.anaget_phbst_and_phdos_files(nqsmall=10, ndivsm=20, asr=2, chneut=1, dipdip=1, dos_method="tetra",
-                                   ngqpt=None, lo_to_splitting=False,
-                                   qptbounds=None, anaddb_kwargs=None)
-phbands, phdos = bstfile.phbands, phdosfile.phdos"""),
-            nbv.new_code_cell("fig = phbands.plot_with_phdos(phdos)"),
-            nbv.new_code_cell("fig = phbands.plot_fatbands(phdos_file=phdosfile)"),
-            nbv.new_code_cell("#fig = phbands.plot_colored_matched()"),
-            nbv.new_code_cell("#fig = phdosfile.plot_pjdos_type()"),
-            nbv.new_code_cell("#fig = phdosfile.plot_pjdos_redirs_type()"),
-            nbv.new_code_cell("#fig = phdosfile.plot_pjdos_redirs_site()"),
+bstfile, phdosfile =  ddb.anaget_phbst_and_phdos_files(nqsmall=10, ndivsm=20,
+    asr=2, chneut=1, dipdip=0, lo_to_splitting=False,
+    dos_method="tetra", ngqpt=None, qptbounds=None, verbose=0, anaddb_kwargs=None)
 
+phbands, phdos = bstfile.phbands, phdosfile.phdos"""),
+            nbv.new_markdown_cell("## q-point path"),
+            nbv.new_code_cell("fig = phbands.qpoints.plot()"),
+            nbv.new_markdown_cell("## Phonon bands with DOS"),
+            nbv.new_code_cell("fig = phbands.plot_with_phdos(phdos, units=units)"),
+            nbv.new_markdown_cell("## Phonon fatbands with DOS"),
+            nbv.new_code_cell("fig = phbands.plot_fatbands(phdos_file=phdosfile, units=units)"),
+            nbv.new_markdown_cell("## Distribution of phonon frequencies wrt mode index"),
+            nbv.new_code_cell("fig = phbands.boxplot(units=units)"),
+            nbv.new_markdown_cell("## Phonon band structure with different color for each line"),
+            nbv.new_code_cell("fig = phbands.plot_colored_matched(units=units)"),
+            nbv.new_markdown_cell("## Type-projected phonon DOS."),
+            nbv.new_code_cell("fig = phdosfile.plot_pjdos_type(units=units)"),
+            nbv.new_markdown_cell("## Type-projected phonon DOS decomposed along the three reduced directions"),
+            nbv.new_code_cell("fig = phdosfile.plot_pjdos_redirs_type(units=units)"),
+            nbv.new_code_cell("#fig = phdosfile.plot_pjdos_redirs_site(units=units)"),
+            nbv.new_markdown_cell("## Thermodinamic properties within the harmonic approximation"),
+            nbv.new_code_cell("fig = phdosfile.phdos.plot_harmonic_thermo(tstart=5, tstop=300)"),
+
+            nbv.new_markdown_cell("## Macroscopic dielectric tensor and Born effective charges"),
+            nbv.new_code_cell("""\
+if False:
+    emacro, becs = ddb.anaget_emacro_and_becs()
+    print(emacro)
+    print(becs)"""),
+
+            nbv.new_markdown_cell("## Call `anaddb` to compute phonon DOS with different BZ samplings"),
             nbv.new_code_cell("""\
 c = None
 if False:
@@ -729,15 +750,23 @@ if c is not None:
     c.plotter.ipw_select_plot()
     #for phdos in c.phdoses"""),
 
+            nbv.new_markdown_cell("## Analysis of the IFCs in real-space"),
             nbv.new_code_cell("""\
-if False:
-    emacro, becs = ddb.anaget_emacro_and_becs()
-    print(emacro)
-    print(becs)"""),
-
-            nbv.new_code_cell("""\
+ifc = None
 if False:
     ifc = ddb.anaget_ifc(ifcout=None, asr=2, chneut=1, dipdip=1, ngqpt=None, verbose=0, anaddb_kwargs=None)"""),
+
+            nbv.new_code_cell("""\
+if ifc is not None:
+    fig = ifc.plot_longitudinal_ifc(atom_indices=None, atom_element=None, neighbour_element=None, min_dist=None,
+                                    max_dist=None, ax=None)"""),
+            nbv.new_code_cell("""\
+if ifc is not None:
+    fig = ifc.plot_longitudinal_ifc_short_range(atom_indices=None, atom_element=None, neighbour_element=None)"""),
+            nbv.new_code_cell("""\
+if ifc is not None:
+    fig = ifc.plot_longitudinal_ifc_ewald(atom_indices=None, atom_element=None, neighbour_element=None,
+                                          min_dist=None, max_dist=None, ax=None)"""),
         ])
 
         return self._write_nb_nbpath(nb, nbpath)

abipy/dfpt/gruneisen.py

@@ -205,13 +205,15 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
         return plotter
 
     @add_fig_kwargs
-    def plot_phbands_with_gruns(self, gamma_fact=1, alpha=0.6, with_doses="all", units="eV",
+    def plot_phbands_with_gruns(self, fill_with="gruns", gamma_fact=1, alpha=0.6, with_doses="all", units="eV",
                                 ylims=None, match_bands=False, **kwargs):
         """
         Plot the phonon bands corresponding to V0 (the central point) with markers
         showing the value and the sign of the Grunesein parameters.
 
         Args:
+            fill_with: Define the quantity used to plot stripes. "gruns" for Grunesein parameters,
+                "groupv" for phonon group velocities.
             gamma_fact: Scaling factor for Grunesein parameters.
                 Up triangle for positive values, down triangles for negative values.
             alpha: The alpha blending value for the markers between 0 (transparent) and 1 (opaque)
@@ -250,27 +252,35 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
         # Plot phonon bands.
         phbands.plot(ax=ax_bands, units=units, match_bands=match_bands, show=False)
 
-        # Plot gruneisen markers on top of band structure.
-        max_omega = np.abs(phbands.phfreqs).max()
+        if fill_with == "gruns":
             max_gamma = np.abs(phbands.grun_vals).max()
+        elif fill_with == "groupv":
+            # TODO: units?
+            dwdq_qpath = self.reader.read_value("gruns_dwdq_qpath")
+            groupv = np.linalg.norm(dwdq_qpath, axis=-1)
+            max_gamma = np.abs(groupv).max()
+        else:
+            raise ValueError("Unsupported fill_with: `%s`" % fill_with)
+
+        # Plot gruneisen markers on top of band structure.
         xvals = np.arange(len(phbands.phfreqs))
+        max_omega = np.abs(phbands.phfreqs).max()
+
         for nu in phbands.branches:
             omegas = phbands.phfreqs[:, nu].copy() * factor
-            sizes = phbands.grun_vals[:, nu].copy() * (gamma_fact * 0.02 * max_omega / max_gamma)
 
+            if fill_with == "gruns":
+                # Must handle positive-negative values
+                sizes = phbands.grun_vals[:, nu].copy() * (gamma_fact * 0.02 * max_omega / max_gamma)
                 yup = omegas + np.where(sizes >= 0, sizes, 0)
                 ydown = omegas + np.where(sizes < 0, sizes, 0)
-
                 ax_bands.fill_between(xvals, omegas, yup, alpha=alpha, facecolor="red")
                 ax_bands.fill_between(xvals, ydown, omegas, alpha=alpha, facecolor="blue")
 
-            # Use different symbols depending on the value of s. Cannot use negative s.
-            #xys = np.array([xos for xos in zip(xvals, omegas, sizes) if xos[2] >= 0]).T.copy()
-            #if xys.size:
-            #    ax_bands.scatter(xys[0], xys[1], s=xys[2], marker="^", label=" >0", color="blue", alpha=alpha)
-            #xys = np.array([xos for xos in zip(xvals, omegas, sizes) if xos[2] < 0]).T.copy()
-            #if xys.size:
-            #    ax_bands.scatter(xys[0], xys[1], s=np.abs(xys[2]), marker="v", label=" <0", color="blue", alpha=alpha)
+            elif fill_with == "groupv":
+                sizes = groupv[:, nu].copy() * (gamma_fact * 0.04 * max_omega / max_gamma)
+                ydown, yup = omegas - sizes / 2, omegas + sizes / 2
+                ax_bands.fill_between(xvals, ydown, yup, alpha=alpha, facecolor="red")
 
         set_axlims(ax_bands, ylims, "x")
 
@@ -304,18 +314,20 @@ class GrunsNcFile(AbinitNcFile, Has_Structure, NotebookWriter):
         nb.cells.extend([
             nbv.new_code_cell("ncfile = abilab.abiopen('%s')" % self.filepath),
             nbv.new_code_cell("print(ncfile)"),
+            nbv.new_code_cell("ncfile.structure"),
 
             nbv.new_code_cell("fig = ncfile.plot_doses()"),
             nbv.new_code_cell("fig = ncfile.plot_phbands_with_gruns()"),
 
+            #nbv.new_code_cell("fig = phbands_qpath_v0.plot_fatbands(phdos_file=phdosfile)"),
             nbv.new_code_cell("plotter = ncfile.get_plotter()\nprint(plotter)"),
             nbv.new_code_cell("df_phbands = plotter.get_phbands_frame()\ndisplay(df_phbands)"),
             nbv.new_code_cell("plotter.ipw_select_plot()"),
 
-            nbv.new_code_cell("gruns_data = ncfile.to_dataframe()"),
+            nbv.new_code_cell("gdata = ncfile.to_dataframe()\ngdata.describe()"),
             nbv.new_code_cell("""\
 #import df_widgets.seabornw as snsw
-#snsw.api_selector(gruns_data)"""),
+#snsw.api_selector(gdata)"""),
         ])
 
         return self._write_nb_nbpath(nb, nbpath)
@@ -402,7 +414,7 @@ class GrunsReader(ETSF_Reader):
         grun_vals = self.read_value("gruns_gvals_qpath")
         freqs_vol = self.read_value("gruns_wvols_qpath") * abu.Ha_eV
         # TODO: Convert?
-        dwdw = self.read_value("gruns_dwdq_qpath")
+        dwdq_qpath = self.read_value("gruns_dwdq_qpath")
 
         amuz = self.read_amuz_dict()
         #print("amuz", amuz)

abipy/dfpt/phonons.py

@@ -1321,7 +1321,7 @@ class PhononBands(object):
     @add_fig_kwargs
     def boxplot(self, ax=None, units="eV", mode_range=None, swarm=False, **kwargs):
         """
-        Use seaborn to draw a box plot to show distributions of eigenvalues with respect to the band index.
+        Use seaborn to draw a box plot to show distributions of eigenvalues with respect to the mode index.
 
         Args:
             ax: matplotlib :class:`Axes` or None if a new figure should be created.
@@ -2651,7 +2651,8 @@ class PhononBandsPlotter(NotebookWriter):
         Return an ipython widget with controllers to select the plot.
         """
         def plot_callback(plot_type, units):
-            return getattr(self, plot_type)(units=units)
+            r = getattr(self, plot_type)(units=units, show=True)
+            if plot_type == "animate": return r
 
         import ipywidgets as ipw
         return ipw.interact_manual(
@@ -2858,7 +2859,7 @@ class PhononDosPlotter(NotebookWriter):
         Return an ipython widget with controllers to select the plot.
         """
         def plot_callback(plot_type, units):
-            return getattr(self, plot_type)(units=units)
+            getattr(self, plot_type)(units=units, show=True)
 
         import ipywidgets as ipw
         return ipw.interact_manual(
@@ -2873,8 +2874,8 @@ class PhononDosPlotter(NotebookWriter):
         from the phonon DOS within the harmonic approximation.
         """
         def plot_callback(tstart, tstop, num, units, formula_units):
-            return self.plot_harmonic_thermo(tstart=tstart, tstop=tstop, num=num,
-                                             units=units, formula_units=formula_units)
+            self.plot_harmonic_thermo(tstart=tstart, tstop=tstop, num=num,
+                                      units=units, formula_units=formula_units, show=True)
 
         import ipywidgets as ipw
         return ipw.interact_manual(
@@ -3188,6 +3189,7 @@ class InteratomicForceConstants(Has_Structure):
 
         ax.set_xlabel('Distance [Bohr]')
         ax.set_ylabel(r'IFC [Ha/Bohr$^2$]')
+        ax.grid(True)
 
         ax.plot(dist, filtered_ifc, **kwargs)
 
@@ -3198,6 +3200,7 @@ class InteratomicForceConstants(Has_Structure):
                               max_dist=None, ax=None, **kwargs):
         """
         Plots the total longitudinal ifcs in local coordinates, filtered according to the optional arguments.
+
         Args:
             atom_indices: a list of atom indices in the structure. Only neighbours of these atoms will be considered.
             atom_element: symbol of an element in the structure. Only neighbours of these atoms will be considered.
@@ -3206,6 +3209,7 @@ class InteratomicForceConstants(Has_Structure):
             max_dist: maximum distance between atoms and neighbours.
             ax: matplotlib :class:`Axes` or None if a new figure should be created.
             kwargs: kwargs passed to the matplotlib function 'plot'. Color defaults to blue, symbol to 'o' and lw to 0
+
         Returns:
             matplotlib figure
         """
@@ -3229,6 +3233,7 @@ class InteratomicForceConstants(Has_Structure):
             max_dist: maximum distance between atoms and neighbours.
             ax: matplotlib :class:`Axes` or None if a new figure should be created.
             kwargs: kwargs passed to the matplotlib function 'plot'. Color defaults to blue, symbol to 'o' and lw to 0
+
         Returns:
             matplotlib figure
         """
@@ -3247,6 +3252,7 @@ class InteratomicForceConstants(Has_Structure):
                                           min_dist=None, max_dist=None, ax=None, **kwargs):
         """
         Plots the Ewald part of the ifcs in local coordinates, filtered according to the optional arguments.
+
         Args:
             atom_indices: a list of atom indices in the structure. Only neighbours of these atoms will be considered.
             atom_element: symbol of an element in the structure. Only neighbours of these atoms will be considered.
@@ -3255,10 +3261,10 @@ class InteratomicForceConstants(Has_Structure):
             max_dist: maximum distance between atoms and neighbours.
             ax: matplotlib :class:`Axes` or None if a new figure should be created.
             kwargs: kwargs passed to the matplotlib function 'plot'. Color defaults to blue, symbol to 'o' and lw to 0
+
         Returns:
             matplotlib figure
         """
-
         if self.local_vectors is None:
             raise ValueError("Local coordinates are missing. Run anaddb with ifcana=1")
 

abipy/dfpt/tests/test_anaddbnc.py

@@ -16,8 +16,7 @@ class AnaddbNcFileTest(AbipyTest):
         anaddbnc_fname = abidata.ref_file("AlAs_nl_dte_anaddb.nc")
 
         with AnaddbNcFile(anaddbnc_fname) as anc:
-            repr(anc)
-            str(anc)
+            repr(anc); str(anc)
             assert anc.structure.formula == "Al1 As1"
             assert anc.becs is not None
             assert anc.emacro is not None

abipy/dfpt/tests/test_phonons.py

@@ -169,6 +169,26 @@ class PhbstFileTest(AbipyTest):
         if self.has_nbformat():
             ncfile.write_notebook(nbpath=self.get_tmpname(text=True))
 
+    def test_phbst_file_with_loto(self):
+        """Testing PHBST file with LOTO terms from anaddb.nc."""
+        with abilab.abiopen(abidata.ref_file("ZnSe_hex_886.out_PHBST.nc")) as ncfile:
+            phbands = ncfile.phbands
+
+        # Phonon frequencies with non analytical contributions, if calculated, are saved
+        # in the anaddb.nc file produced by anaddb. The results should be fetched from there
+        # and added to the phonon bands.
+        phbands.read_non_anal_from_file(abidata.ref_file("ZnSe_hex_886.anaddb.nc"))
+        nana = phbands.non_anal_ph
+        assert nana.structure == phbands.structure
+        print(nana.structure.reciprocal_lattice)
+        self.assert_almost_equal(nana.directions.flat,
+                [0.1234510847, -0.071274517, 0, 0.1646014463, 0, 0, 0, 0, 0.0751616546])
+        for direc in nana.directions:
+            assert nana.has_direction(direc, cartesian=True)
+
+        if self.has_matplotlib():
+            phbands.plot(title="ZnSe with LO-TO splitting", show=False)
+
 
 class PhononBandsPlotterTest(AbipyTest):
 
@@ -361,7 +381,7 @@ class NonAnalyticalPhTest(AbipyTest):
             phbands = ddb.anaget_phmodes_at_qpoint(qpoint=[0, 0, 0], lo_to_splitting=True)
 
             assert phbands.non_anal_ph is not None
-            str(phbands.non_anal_ph)
+            repr(phbands.non_anal_ph); str(phbands.non_anal_ph)
             assert phbands.structure == phbands.non_anal_ph.structure
             #assert phbands.non_anal_ph.has_direction(direction= cartesian=False)
 

abipy/electrons/ebands.py

@@ -1,5 +1,5 @@
 # coding: utf-8
-"""Classes for the analysis of electronic structures."""
+"""Classes to analyse electronic structures."""
 from __future__ import print_function, division, unicode_literals, absolute_import
 
 import sys
@@ -2517,12 +2517,13 @@ class ElectronBandsPlotter(NotebookWriter):
         Return an ipython widget with controllers to select the plot.
         """
         def plot_callback(plot_type, e0):
-            return getattr(self, plot_type)(e0=e0)
+            r = getattr(self, plot_type)(e0=e0, show=True)
+            if plot_type == "animate": return r
 
         import ipywidgets as ipw
         return ipw.interact_manual(
                 plot_callback,
-                plot_type=["combiplot", "gridplot", "boxplot", "combiboxplot"],
+                plot_type=["combiplot", "gridplot", "boxplot", "combiboxplot", "animate"],
                 e0=["fermie", "0.0"],
             )
 
@@ -3097,7 +3098,7 @@ class ElectronDosPlotter(NotebookWriter):
         Return an ipython widget with controllers to select the plot.
         """
         def plot_callback(plot_type, e0):
-            return getattr(self, plot_type)(e0=e0)
+            getattr(self, plot_type)(e0=e0, show=True)
 
         import ipywidgets as ipw
         return ipw.interact_manual(

abipy/electrons/fatbands.py

@@ -1279,6 +1279,11 @@ class FatBandsFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWrite
             # Compute PJDOS from self.
             pjdosfile = self
 
+        if not pjdosfile.ebands.kpoints.is_ibz:
+            cprint("DOS requires k-points in the IBZ but got pjdosfile: %s" % repr(pjdosfile), "yellow")
+            cprint("Returning None", "yellow")
+            return None
+
         if edos_kwargs is None: edos_kwargs = {}
 
         # Build plot grid.
@@ -1532,9 +1537,10 @@ class FatBandsFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWrite
         nb.cells.extend([
             nbv.new_code_cell("fbnc = abilab.abiopen('%s')\nprint(fbnc)" % self.filepath),
             nbv.new_code_cell("fbnc.structure"),
-            nbv.new_code_cell("fig = fbnc.ebands.kpoints.plot()"),
+            nbv.new_code_cell("# fig = fbnc.ebands.kpoints.plot()"),
             nbv.new_code_cell("xlims = (None, None)\nylims = (None, None)"),
-            #nbv.new_code_cell("fig = fbnc.ebands.plot(ylims=ylims)"),
+            nbv.new_code_cell("# fig = fbnc.ebands.plot(ylims=ylims)"),
+            nbv.new_code_cell("fig = fbnc.ebands.boxplot()"),
         ])
 
         if self.prtdos == 3: