diff --git a/abipy/core/kpoints.py b/abipy/core/kpoints.py
index 7c7675ab6..20c802947 100644
--- a/abipy/core/kpoints.py
+++ b/abipy/core/kpoints.py
@@ -176,7 +176,7 @@ def kmesh_from_mpdivs(mpdivs, shifts, pbc=False, order="bz"):
         mpdivs: The three MP divisions.
         shifts: Array-like object with the MP shift.
         pbc: If True, periodic images of the k-points will be included i.e. closed mesh.
-        order: "unit_cell" if the kpoint coordinates must be in [0,1)
+        order: "unit_cell" if the kpoint coordinates must be in [0, 1)
                "bz" if the kpoint coordinates must be in [-1/2, +1/2)
     """
     shifts = np.reshape(shifts, (-1, 3))
@@ -220,6 +220,7 @@ def map_grid2ibz(structure, ibz, ngkpt, has_timrev, pbc=False):
     if abispg is None:
         raise ValueError("Structure does not contain Abinit spacegroup info!")
 
+
     # Extract rotations in reciprocal space (FM part).
     symrec_fm = [o.rot_g for o in abispg.fm_symmops]
 
@@ -242,27 +243,15 @@ def map_grid2ibz(structure, ibz, ngkpt, has_timrev, pbc=False):
 
     if np.any(bzgrid2ibz == -1):
         #for ik_bz, ik_ibz in enumerate(self.bzgrid2ibz): print(ik_bz, ">>>", ik_ibz)
-        msg = "Found %s/%s invalid entries in bzgrid2ibz array" % ((bzgrid2ibz == -1).sum(), bzgrid2ibz.size)
-        msg += "This can happen if there an inconsistency between the input IBZ and ngkpt"
-        msg += "ngkpt: %s, has_timrev: %s" % (str(ngkpt), has_timrev)
+        msg =  " Found %s/%s invalid entries in bzgrid2ibz array\n" % ((bzgrid2ibz == -1).sum(), bzgrid2ibz.size)
+        msg += " This can happen if there is an inconsistency between the input IBZ and ngkpt\n"
+        msg += " ngkpt: %s, has_timrev: %s\n" % (str(ngkpt), has_timrev)
+        msg += f" {abispg=}\n"
         raise ValueError(msg)
 
     bz2ibz = bzgrid2ibz.flatten()
     return bz2ibz
 
-    """
-    for ik_bz, kbz in enumerate(bz):
-        found = False
-        for ik_ibz, kibz in enumerate(ibz):
-            if found: break
-            for symmop in structure.spacegroup:
-                krot = symmop.rotate_k(kibz)
-                if issamek(krot, kbz):
-                    bz2ibz[ik_bz] = ik_ibz
-                    found = True
-                    break
-    """
-
 
 def has_timrev_from_kptopt(kptopt):
     """
diff --git a/abipy/core/structure.py b/abipy/core/structure.py
index 42921bfa1..7233114e2 100644
--- a/abipy/core/structure.py
+++ b/abipy/core/structure.py
@@ -20,7 +20,7 @@
 from monty.functools import lazy_property
 from monty.string import is_string, marquee, list_strings
 from monty.termcolor import cprint
-from monty.dev import deprecated
+#from monty.dev import deprecated
 from pymatgen.core.structure import Structure as pmg_Structure
 from pymatgen.core.sites import PeriodicSite
 from pymatgen.core.lattice import Lattice
@@ -29,7 +29,7 @@
 from abipy.core.symmetries import AbinitSpaceGroup
 from abipy.tools.plotting import add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt, add_plotly_fig_kwargs
 from abipy.iotools import as_etsfreader, Visualizer
-
+from abipy.tools.typing import Figure
 
 
 __all__ = [
@@ -170,14 +170,48 @@ def display_structure(obj, **kwargs):
     return nbjsmol_display(structure.to(fmt="cif"), ext=".cif", **kwargs)
 
 
-class Structure(pmg_Structure, NotebookWriter):
+def get_structures_from_file(filepath: PathLike, index) -> list[Structure]:
     """
-    Extends :class:`pymatgen.core.structure.Structure` with Abinit-specific methods.
+    """
+    #if index is None:
+    #    index = -1
+
+    if filepath.endswith(".traj"):
+        # ASE trajectory file.
+        from ase.atoms import Atoms
+        from ase.io import read
+        atoms_list = read(filepath, index=index)
+        if not isinstance(atoms_list, list):
+            assert isinstance(atoms_list, Atoms)
+            atoms_list = [atoms_list]
+
+    # TODO: HIST.nc and XDATCAR
+    #elif filepath.endswith("_HIST.nc"):
+    #    from abipy.dynamics.hist import HistFile
+    #    with HistFile(filepath) as hist:
+    #        return hist.read_structures(index)
+
+    #elif filepath.endswith("XDATCAR"):
+
+    else:
+        raise NotImplementedError(f"Don't know how to extract structures with index from {filepath=}")
+
+    return [Structure.as_structure(atoms) for atoms in atoms_list]
 
-    A jupyter_ notebook documenting the usage of this object is available at
-    <https://nbviewer.jupyter.org/github/abinit/abitutorials/blob/master/abitutorials/structure.ipynb>
 
-    For the pymatgen project see :cite:`Ong2013`.
+def get_first_and_last_structure_from_file(filepath: PathLike) -> tuple[Structure]:
+    """
+    Helper function to read the first and the last structure from a trajectory file.
+    Simplified wrapper around get_structures_from_file.
+    """
+    first_structure = get_structures_from_file(filepath, index=0)[0]
+    last_structure = get_structures_from_file(filepath, index=-1)[0]
+    return first_structure, last_structure
+
+
+class Structure(pmg_Structure, NotebookWriter):
+    """
+    Extends :class:`pymatgen.core.structure.Structure` with Abinit-specific methods.
 
     .. rubric:: Inheritance Diagram
     .. inheritance-diagram:: Structure
@@ -277,8 +311,7 @@ def from_file(cls, filepath: str, primitive: bool = False, sort: bool = False) -
             if closeit: ncfile.close()
 
         elif filepath.endswith(".abi") or filepath.endswith(".in"):
-            # Abinit input file.
-            # Here I assume that the input file contains a single structure.
+            # Abinit input file. Here I assume that the input file contains a single structure.
             from abipy.abio.abivars import AbinitInputFile
             return AbinitInputFile.from_file(filepath).structure
 
@@ -291,6 +324,7 @@ def from_file(cls, filepath: str, primitive: bool = False, sort: bool = False) -
                 #print("initial_structures:\n", out.initial_structures, "\nfinal_structures:\n", out.final_structures)
                 if out.final_structures: return out.final_structure
                 if out.initial_structures: return out.initial_structure
+
             raise ValueError("Cannot find structure in Abinit output file `%s`" % filepath)
 
         elif filepath.endswith(".abivars") or filepath.endswith(".ucell"):
@@ -323,8 +357,7 @@ def from_file(cls, filepath: str, primitive: bool = False, sort: bool = False) -
                 from ase.io import read
             except ImportError as exc:
                 raise RuntimeError("ase is required to read xyz files. Use `pip install ase`") from exc
-            atoms = read(filepath)
-            return cls.as_structure(atoms)
+            return cls.as_structure(read(filepath))
 
         else:
             # Invoke pymatgen and change class. Note that AbinitSpacegroup is missing here.
@@ -389,7 +422,6 @@ def from_ase_atoms(cls, atoms) -> Structure:
     @property
     def n_elems(self) -> int:
         """Number of types of atoms."""
-
         return len(self.types_of_species)
 
     def to_ase_atoms(self, calc=None):
@@ -508,7 +540,6 @@ def zincblende(cls, a, species, units="ang", **kwargs) -> Structure:
             kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
 
         Example::
-
             Structure.zincblende(a, ["Zn", "S"])
         """
         a = pmg_units.Length(a, units).to("ang")
@@ -532,7 +563,6 @@ def rocksalt(cls, a, species, units="ang", **kwargs) -> Structure:
             kwargs: All keyword arguments accepted by :class:`pymatgen.Structure`
 
         Example::
-
             Structure.rocksalt(a, ["Na", "Cl"])
         """
         a = pmg_units.Length(a, units).to("ang")
@@ -570,10 +600,9 @@ def ABO3(cls, a, species, units="ang", **kwargs) -> Structure:
     @classmethod
     def from_abistring(cls, string: str) -> Structure:
         """
-        Initialize Structure from string with Abinit input variables.
+        Initialize Structure from a string with Abinit input variables.
         """
         from abipy.abio.abivars import AbinitInputFile, structure_from_abistruct_fmt
-
         if "xred_symbols" not in string:
             # Standard (verbose) input file with znucl, typat etc.
             return AbinitInputFile.from_string(string).structure
@@ -679,12 +708,10 @@ def write_cif_with_spglib_symms(self, filename, symprec=1e-3, angle_tolerance=5.
             symprec (float): If not none, finds the symmetry of the structure
                 and writes the cif with symmetry information. Passes symprec
                 to the SpacegroupAnalyzer.
-            significant_figures (int): Specifies precision for formatting of floats.
-                Defaults to 8.
+            significant_figures (int): Specifies precision for formatting of floats. Defaults to 8.
             angle_tolerance (float): Angle tolerance for symmetry finding. Passes
-                angle_tolerance to the SpacegroupAnalyzer. Used only if symprec
-                is not None.
-            ret_string: True to return string
+                angle_tolerance to the SpacegroupAnalyzer. Used only if symprec is not None.
+            ret_string: True to return string.
         """
         from pymatgen.io.cif import CifWriter
         cif_str = str(CifWriter(self,
@@ -722,7 +749,7 @@ def to_abivars(self, enforce_znucl=None, enforce_typat=None, **kwargs) -> dict:
 
     @property
     def latex_formula(self) -> str:
-        """LaTeX formatted formula. E.g., Fe2O3 is transformed to Fe$_{2}$O$_{3}$."""
+        """LaTeX formatted formula: e.g., Fe2O3 is transformed to Fe$_{2}$O$_{3}$."""
         from pymatgen.util.string import latexify
         return latexify(self.formula)
 
@@ -743,7 +770,6 @@ def get_abi_string(self, fmt: str = "abinit_input") -> str:
         fmt="abicell" is the lightweight format that uses `xred_symbols`
         This format can be used to include the structure in the input via the structure "abivars:FILE" syntax.
         """
-
         lines = []
         app = lines.append
 
@@ -811,8 +837,8 @@ def get_panel(self, with_inputs=True, **kwargs):
     def get_conventional_standard_structure(self, international_monoclinic=True,
                                            symprec=1e-3, angle_tolerance=5) -> Structure:
         """
-        Gives a structure with a conventional cell according to certain
-        standards. The standards are defined in :cite:`Setyawan2010`
+        Gives a structure with a conventional cell according to certain standards.
+        The standards are defined in :cite:`Setyawan2010`
         They basically enforce as much as possible norm(a1) < norm(a2) < norm(a3)
 
         Returns: The structure in a conventional standardized cell
@@ -822,8 +848,7 @@ def get_conventional_standard_structure(self, international_monoclinic=True,
         return self.__class__.as_structure(new)
 
     def abi_primitive(self, symprec=1e-3, angle_tolerance=5, no_idealize=0) -> Structure:
-        #TODO: this should be moved to pymatgen in the get_refined_structure or so ...
-        # to be considered in February 2016
+        #TODO: this should be moved to pymatgen in the get_refined_structure or so ... to be considered in February 2016
         import spglib
         from pymatgen.io.ase import AseAtomsAdaptor
         try:
@@ -943,11 +968,6 @@ def abi_sanitize(self, symprec=1e-3, angle_tolerance=5,
 
         return self.__class__.as_structure(structure)
 
-    #def relax(self, **kwargs) -> Structure:
-    #    new = super().relax(**kwargs)
-    #    new.__class__ = self.__class__
-    #    return new
-
     def get_oxi_state_decorated(self, **kwargs) -> Structure:
         """
         Use :class:`pymatgen.analysis.bond_valence.BVAnalyzer` to estimate oxidation states
@@ -976,24 +996,10 @@ def reciprocal_lattice(self):
         """
         return self._lattice.reciprocal_lattice
 
-    @deprecated(message="lattice_vectors is deprecated and will be removed in abipy 1.0")
-    def lattice_vectors(self, space="r"):
-        """
-        Returns the vectors of the unit cell in Angstrom.
-
-        Args:
-            space: "r" for real space vectors, "g" for reciprocal space basis vectors.
-        """
-        if space.lower() == "r":
-            return self.lattice.matrix
-        if space.lower() == "g":
-            return self.lattice.reciprocal_lattice.matrix
-        raise ValueError("Wrong value for space: %s " % str(space))
-
     def spget_lattice_type(self, symprec=1e-3, angle_tolerance=5) -> str:
         """
         Call spglib to get the lattice for the structure, e.g., (triclinic,
-        orthorhombic, cubic, etc.).This is the same than the
+        orthorhombic, cubic, etc.). This is the same than the
         crystal system with the exception of the hexagonal/rhombohedral lattice
 
         Args:
@@ -1627,7 +1633,7 @@ def plotly_bz(self, fig=None, pmg_path=True, with_labels=True, **kwargs):
 
     @add_fig_kwargs
     def plot_xrd(self, wavelength="CuKa", symprec=0, debye_waller_factors=None,
-                 two_theta_range=(0, 90), annotate_peaks=True, ax=None, **kwargs):
+                 two_theta_range=(0, 90), annotate_peaks=True, ax=None, **kwargs) -> Figure:
         """
         Use pymatgen :class:`XRDCalculator` to show the XRD plot.
 
diff --git a/abipy/core/tests/test_structure.py b/abipy/core/tests/test_structure.py
index 4a03934ef..35a204d29 100644
--- a/abipy/core/tests/test_structure.py
+++ b/abipy/core/tests/test_structure.py
@@ -224,9 +224,9 @@ def test_utils(self):
 
         # Temporarily disables ad webserver is down.
         #if self.is_url_reachable("www.crystallography.net"):
-        #    mgb2_cod = Structure.from_cod_id(1526507, primitive=True)
-        #    assert mgb2_cod.formula == "Mg1 B2"
-        #    assert mgb2_cod.spget_lattice_type() == "hexagonal"
+        mgb2_cod = Structure.from_cod_id(1526507, primitive=True)
+        assert mgb2_cod.formula == "Mg1 B2"
+        assert mgb2_cod.spget_lattice_type() == "hexagonal"
 
         mgb2 = abidata.structure_from_ucell("MgB2")
         if self.has_ase():
diff --git a/abipy/dynamics/hist.py b/abipy/dynamics/hist.py
index d2204d12a..7b27dd1c8 100644
--- a/abipy/dynamics/hist.py
+++ b/abipy/dynamics/hist.py
@@ -59,6 +59,8 @@ def params(self) -> dict:
     def __str__(self) -> str:
         return self.to_string()
 
+    #def read_structures(self, index: str):
+
     # TODO: Add more metadata.
     #@lazy_property
     #def nsppol(self):
diff --git a/abipy/electrons/ebands.py b/abipy/electrons/ebands.py
index ab7aacd08..59412e7d5 100644
--- a/abipy/electrons/ebands.py
+++ b/abipy/electrons/ebands.py
@@ -28,7 +28,7 @@
 from abipy.core.func1d import Function1D
 from abipy.core.mixins import Has_Structure, NotebookWriter
 from abipy.core.kpoints import (Kpoint, KpointList, Kpath, IrredZone, KSamplingInfo, KpointsReaderMixin,
-    Ktables, has_timrev_from_kptopt, map_grid2ibz) #, kmesh_from_mpdivs)
+    Ktables, has_timrev_from_kptopt, map_grid2ibz, kmesh_from_mpdivs)
 from abipy.core.structure import Structure
 from abipy.iotools import ETSF_Reader
 from abipy.tools import duck
@@ -971,6 +971,29 @@ def has_timrev(self) -> bool:
         """True if time-reversal symmetry is used in the BZ sampling."""
         return has_timrev_from_kptopt(self.kptopt)
 
+    def get_bz2ibz_bz_points(self, require_gamma_centered=False):
+        """
+        Return mapping bz2ibz and the list of k-points in the BZ.
+
+        Args:
+            require_gamma_centered: True if the k-mesh should be Gamma-centered
+        """
+        err_msg = self.isnot_ibz_sampling(require_gamma_centered=require_gamma_centered)
+        if err_msg:
+            raise TypeError(err_msg)
+
+        if not self.kpoints.is_mpmesh:
+            raise ValueError("Only Monkhorst-Pack meshes are supported")
+
+        ngkpt, shifts = self.kpoints.mpdivs_shifts
+        print(f"{ngkpt = }, {shifts =}")
+        # TODO: Handle shifts
+
+        bz2ibz = map_grid2ibz(self.structure, self.kpoints.frac_coords, ngkpt, self.has_timrev)
+        bz_kpoints = kmesh_from_mpdivs(ngkpt, shifts)
+
+        return bz2ibz, bz_kpoints
+
     def isnot_ibz_sampling(self, require_gamma_centered=False) -> str:
         """
         Test whether the k-points in the band structure represent an IBZ with an associated k-mesh
@@ -5414,7 +5437,7 @@ def xcrysden_view(self):  # pragma: no cover
         from abipy.iotools.visualizer import Xcrysden
         return Xcrysden(tmp_filepath)()
 
-    def to_bxsf(self, filepath, unit="eV") -> str:
+    def to_bxsf(self, filepath: str, unit: str="eV") -> str:
         """
         Export the full band structure to ``filepath`` in BXSF format
         suitable for the visualization of the Fermi surface with xcrysden_ (use ``xcrysden --bxsf FILE``).
diff --git a/abipy/eph/varpeq.py b/abipy/eph/varpeq.py
index d91ffb486..4e9125e73 100644
--- a/abipy/eph/varpeq.py
+++ b/abipy/eph/varpeq.py
@@ -17,7 +17,7 @@
 from monty.termcolor import cprint
 from abipy.core.func1d import Function1D
 from abipy.core.structure import Structure
-from abipy.core.kpoints import kpoints_indices
+from abipy.core.kpoints import kpoints_indices, map_grid2ibz, kmesh_from_mpdivs
 from abipy.core.mixins import AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter
 from abipy.tools.typing import PathLike
 from abipy.tools.plotting import (add_fig_kwargs, get_ax_fig_plt, get_axarray_fig_plt, set_axlims, set_visible,
@@ -28,6 +28,7 @@
 from abipy.dfpt.ddb import DdbFile
 from abipy.tools.typing import Figure
 from abipy.tools.numtools import BzRegularGridInterpolator, gaussian
+from abipy.iotools import bxsf_write
 from abipy.abio.robots import Robot
 from abipy.eph.common import BaseEphReader
 
@@ -103,6 +104,8 @@ class Entry:
 # Convert to dictionary: name --> Entry
 _ALL_ENTRIES = {e.name: e for e in _ALL_ENTRIES}
 
+# TODO:
+# Handle multiple states.
 
 class VarpeqFile(AbinitNcFile, Has_Structure, Has_ElectronBands, NotebookWriter):
     """
@@ -196,7 +199,7 @@ def to_string(self, verbose=0) -> str:
     def get_last_iteration_dict_ev(self, spin: int) -> dict:
         """
         Return dictionary mapping the latex label to the value of the last iteration
-        for the given spin index. All energies are in eV.
+        for the given spin. All energies are in eV.
         """
         nstep2cv_spin = self.r.read_value('nstep2cv')
         iter_rec_spin = self.r.read_value('iter_rec')
@@ -252,8 +255,7 @@ def plot_scf_cycle(self, ax_mat=None, fontsize=8, **kwargs) -> Figure:
                         ax.set_yscale("log")
 
                 ax.set_xlim(1, nstep2cv)
-                set_grid_legend(ax, fontsize, xlabel="Iteration",
-                                title=self.get_title_spin(spin),
+                set_grid_legend(ax, fontsize, xlabel="Iteration", title=self.get_title_spin(spin),
                                 ylabel="Energy (eV)" if iax == 0 else r"$|\Delta|$ Energy (eV)",
                 )
 
@@ -284,19 +286,18 @@ def write_notebook(self, nbpath=None) -> str:
 class Polaron:
     """
     This object stores the polaron coefficients A_kn, B_qnu for a given spin.
-    Provides methods to plot |A_nk|^2 or |B_qnu|^2 together with the band structures
-    (fatbands-like plots).
+    Provides methods to plot |A_nk|^2 or |B_qnu|^2 together with the band structures (fatbands-like plots).
     """
 
     spin: int          # Spin index.
-    nb: int            # Number of bands.
-    nk: int            # Number of k-points (including filtering if any)
-    nq: int            # Number of q-points (including filtering if any)
+    nb: int            # Number of bands in A_kn,
+    nk: int            # Number of k-points in A_kn, (including filtering if any)
+    nq: int            # Number of q-points in B_qnu (including filtering if any)
     bstart: int        # First band starts at bstart
-    bstop: int
+    bstop: int         # Last band (python convention)
 
-    kpoints: np.ndarray
-    qpoints: np.ndarray
+    kpoints: np.ndarray   # Reduced coordinates of the k-points
+    qpoints: np.ndarray   # Reduced coordinates of the q-points
     a_kn: np.ndarray
     b_qnu: np.ndarray
     varpeq: VarpeqFile
@@ -313,8 +314,8 @@ def from_varpeq(cls, varpeq: VarpeqFile, spin: int) -> Polaron:
         qpoints = r.read_value("qpts_spin")[spin, :nq]
         a_kn = r.read_value("a_spin", cmode="c")[spin, :nk, :nb]
         b_qnu = r.read_value("b_spin", cmode="c")[spin, :nq]
-
         data = locals()
+
         return cls(**{k: data[k] for k in [field.name for field in dataclasses.fields(Polaron)]})
 
     @property
@@ -365,69 +366,130 @@ def get_title(self, with_gaps: bool=True) -> str:
             gaps_string = self.varpeq.ebands.get_gaps_string()
             return f"{pre}{varpeq.r.varpeq_pkind}, {gaps_string}"
 
-    def get_a2_interpolator(self, check_mesh: int=0) -> BzRegularGridInterpolator:
+    @lazy_property
+    def ngkpt_and_shifts(self):
         """
-        Build and return an interpolator for |A_nk|^2
-
-        Args:
-            check_mesh: Check whether k-points belong to the mesh if != 0.
+        Return k-mesh divisions and shifts.
         """
-        # Need to know the size of the k-mesh.
         ksampling = self.ebands.kpoints.ksampling
         ngkpt, shifts = ksampling.mpdivs, ksampling.shifts
 
         if ngkpt is None:
             raise ValueError("Non diagonal k-meshes are not supported!")
+
         if len(shifts) > 1:
             raise ValueError("Multiple k-shifts are not supported!")
 
-        k_indices = kpoints_indices(self.kpoints, ngkpt, check_mesh=check_mesh)
+        return ngkpt, shifts
 
+    def insert_a_inbox(self, fill_value=None):
+        """
+        """
+        # Need to know the size of the k-mesh.
+        ngkpt, shifts = self.ngkpt_and_shifts
+        k_indices = kpoints_indices(self.kpoints, ngkpt)
         nx, ny, nz = ngkpt
-        a_data = np.empty((self.nb, nx, ny, nz), dtype=complex)
+
+        shape = (self.nb, nx, ny, nz)
+        if fill_value is None:
+            a_data = np.empty(shape, dtype=complex)
+        else:
+            a_data = np.full(shape, fill_value, dtype=complex)
+
         for ib in range(self.nb):
             for a_cplx, k_inds in zip(self.a_kn[:,ib], k_indices):
                 ix, iy, iz = k_inds
                 a_data[ib, ix, iy, iz] = a_cplx
 
-        return BzRegularGridInterpolator(self.structure, shifts, np.abs(a_data) ** 2, method="linear")
+        return a_data, ngkpt, shifts
 
-    def get_b2_interpolator(self, check_mesh: int=0) -> BzRegularGridInterpolator:
+    def insert_b_inbox(self, fill_value=None):
         """
-        Build and return an interpolator for |B_qnu|^2.
-
-        Args
-            check_mesh: Check whether k-points belong to the mesh if != 0.
         """
         # Need to know the shape of the q-mesh (always Gamma-centered)
         ngqpt, shifts = self.varpeq.r.ngqpt, [0, 0, 0]
-        q_indices = kpoints_indices(self.qpoints, ngqpt, check_mesh=check_mesh)
+        q_indices = kpoints_indices(self.qpoints, ngqpt)
 
         natom3 = 3 * len(self.structure)
         nx, ny, nz = ngqpt
-        b_data = np.empty((natom3, nx, ny, nz), dtype=complex)
+
+        shape = (natom3, nx, ny, nz)
+        if fill_value is None:
+            b_data = np.empty(shape, dtype=complex)
+        else:
+            b_data = np.full(shape, fill_value, dtype=complex)
+
         for nu in range(natom3):
             for b_cplx, q_inds in zip(self.b_qnu[:,nu], q_indices):
                 ix, iy, iz = q_inds
                 b_data[nu, ix, iy, iz] = b_cplx
 
+        return b_data, ngqpt, shifts
+
+    def get_a2_interpolator(self) -> BzRegularGridInterpolator:
+        """
+        Build and return an interpolator for |A_nk|^2
+        """
+        a_data, ngkpt, shifts = self.insert_a_inbox()
+
+        return BzRegularGridInterpolator(self.structure, shifts, np.abs(a_data) ** 2, method="linear")
+
+    def get_b2_interpolator(self) -> BzRegularGridInterpolator:
+        """
+        Build and return an interpolator for |B_qnu|^2.
+        """
+        b_data, ngqpt, shifts = self.insert_b_inbox()
+
         return BzRegularGridInterpolator(self.structure, shifts, np.abs(b_data) ** 2, method="linear")
 
-    @add_fig_kwargs
-    def plot_bz_sampling(self, what="kpoints", fold=False,
-                        ax=None, pmg_path=True, with_labels=True, **kwargs) -> Figure:
+    def write_a2_bxsf(self, filepath: PathLike) -> None:
         """
-        Plots a 3D representation of the Brillouin zone with the sampling.
+        Export |A_nk|^2 in BXSF format suitable for visualization with xcrysden (use ``xcrysden --bxsf FILE``).
+        Require gamma-centered k-mesh.
 
         Args:
-            what: "kpoints" or "qpoints"
-            fold: whether the points should be folded inside the first Brillouin Zone.
-                Defaults to False.
+            filepath: BXSF filename.
+        """
+        # NB: kmesh must be gamma-centered, multiple shifts are not supported.
+        # Init values with 0. This is relevant only if kfiltering is being used.
+        a_data, ngkpt, shifts = self.insert_a_inbox(fill_value=0.0)
+        nkbz = np.product(ngkpt)
+        a2_data = np.abs(a_data) ** 2
+        fermie = a2_data.mean()
+
+        bxsf_write(filepath, self.structure, 1, num_states, ngkpt, a2_data, fermie, unit="Ha")
+
+    def write_b2_bxsf(self, filepath: PathLike) -> None:
         """
-        bz_points = dict(kpoints=self.kpoints, qpoints=self.qpoints)[what]
-        kws = dict(ax=ax, pmg_path=pmg_path, with_labels=with_labels, fold=fold, kpoints=bz_points)
+        Export |B_qnu|^2 in BXSF format suitable for visualization with xcrysden (use ``xcrysden --bxsf FILE``).
 
-        return self.structure.plot_bz(show=False, **kws)
+        Args:
+            filepath: BXSF filename.
+        """
+        # NB: qmesh must be gamma-centered, multiple shifts are not supported.
+        # Init values with 0. This is relevant only if kfiltering is being used.
+        b_data, ngqpt, shifts = self.insert_b_inbox(fill_value=0.0)
+        nqbz = np.product(nqkpt)
+        b2_data = np.abs(b_data) ** 2
+        fermie = b2_data.mean()
+
+        bxsf_write(filepath, self.structure, 1, num_states, ngqpt, b2_data, fermie, unit="Ha")
+
+    #@add_fig_kwargs
+    #def plot_bz_sampling(self, what="kpoints", fold=False,
+    #                     ax=None, pmg_path=True, with_labels=True, **kwargs) -> Figure:
+    #    """
+    #    Plots a 3D representation of the Brillouin zone with the sampling.
+
+    #    Args:
+    #        what: "kpoints" or "qpoints"
+    #        fold: whether the points should be folded inside the first Brillouin Zone.
+    #            Defaults to False.
+    #    """
+    #    bz_points = dict(kpoints=self.kpoints, qpoints=self.qpoints)[what]
+    #    kws = dict(ax=ax, pmg_path=pmg_path, with_labels=with_labels, fold=fold, kpoints=bz_points)
+
+    #    return self.structure.plot_bz(show=False, **kws)
 
     @add_fig_kwargs
     def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
@@ -493,39 +555,53 @@ def plot_ank_with_ebands(self, ebands_kpath, ebands_kmesh=None,
 
         if ebands_kmesh is None:
             print(f"Computing ebands_kmesh with star-function interpolation and {nksmall=}")
-            kmesh = self.structure.calc_ngkpt(nksmall)
-            r = self.ebands.interpolate(lpratio=lpratio, vertices_names=vertices_names, kmesh=kmesh)
+            this_ngkpt = self.structure.calc_ngkpt(nksmall)
+            r = self.ebands.interpolate(lpratio=lpratio, vertices_names=vertices_names, kmesh=this_ngkpt)
             ebands_kmesh = r.ebands_kmesh
 
-        # Get electronic DOS.
+        # Get electronic DOS from ebands_kmesh.
         edos_kws = dict(method=method, step=step, width=width)
         edos = ebands_kmesh.get_edos(**edos_kws)
         mesh = edos.spin_dos[self.spin].mesh
         ank_dos = np.zeros(len(mesh))
         e0 = self.ebands.fermie
 
-        #################
-        # Compute Ank DOS
-        #################
+        ##################
+        # Compute A_nk DOS
+        ##################
+        # FIXME
         # NB: This is just to sketch the ideas. I don't think the present version
         # is correct as only the k --> -k symmetry can be used.
 
-        for ik, kpoint in enumerate(ebands_kmesh.kpoints):
+        for ik_ibz, kpoint in enumerate(ebands_kmesh.kpoints):
             weight = kpoint.weight
-            enes_n = ebands_kmesh.eigens[self.spin, ik, self.bstart:self.bstop]
+            enes_n = ebands_kmesh.eigens[self.spin, ik_ibz, self.bstart:self.bstop]
             a2_n = a2_interp.eval_kpoint(kpoint)
             for e, a2 in zip(enes_n, a2_n):
                 ank_dos += weight * a2 * gaussian(mesh, width, center=e-e0)
 
+        # TODO New version using the BZ. Requires new VARPEQ.nc file with all symmetries
+        # The A_nk do not necessarily have the symmetry of the lattice so we have to loop over the full BZ.
+        # Get mapping BZ --> IBZ needed to obtain the KS eigenvalues e_nk from the IBZ for the DOS
+        """
+        bz2ibz, bz_kpoints = ebands_kmesh.get_bz2ibz_bz_points()
+
+        for ik_ibz, kpoint in zip(bz2ibz, bz_kpoints):
+            enes_n = ebands_kmesh.eigens[self.spin, ik_ibz, self.bstart:self.bstop]
+            a2_n = a2_interp.eval_kpoint(kpoint)
+            for e, a2 in zip(enes_n, a2_n):
+                ank_dos += a2 * gaussian(mesh, width, center=e-e0)
+        ank_dos /= np.product(ngkpt)
+        """
+
         ank_dos = Function1D(mesh, ank_dos)
         print("A2(E) integrates to:", ank_dos.integral_value, " Ideally, it should be 1.")
 
-        edos_opts = {"color": "black",} if self.spin == 0 else {"color": "red",  }
+        edos_opts = {"color": "black",} if self.spin == 0 else {"color": "red"}
 
         ax = ax_list[1]
         edos.plot_ax(ax, e0, spin=self.spin, normalize=normalize, exchange_xy=True, label="eDOS(E)", **edos_opts)
         ank_dos.plot_ax(ax, exchange_xy=True, normalize=normalize, label=r"$A^2$(E)", color=color)
-
         set_grid_legend(ax, fontsize, xlabel="arb. unit")
 
         if ylims is None:
@@ -586,7 +662,7 @@ def plot_bqnu_with_phbands(self, phbands_qpath, phdos_file=None,
             with_title: True to add title with chemical formula and gaps.
             ax_list: List of |matplotlib-Axes| or None if a new figure should be created.
             scale: Scaling factor for |B_qnu|^2.
-            fontsize: fontsize for legends and titles
+            fontsize: fontsize for legends and titles.
         """
         with_phdos = phdos_file is not None  and ddb is not None
         nrows, ncols, gridspec_kw = 1, 1, None
@@ -619,6 +695,10 @@ def plot_bqnu_with_phbands(self, phbands_qpath, phdos_file=None,
                 fig.suptitle(self.get_title(with_gaps=True))
             return fig
 
+        ##################
+        # Compute B_qnu DOS
+        ##################
+
         # Add phdos and |B_qn| dos. Mesh is given in eV, values are in states/eV.
         phdos = phdos_file.phdos
         ngqpt = np.diagonal(phdos_file.qptrlatt)
@@ -629,6 +709,9 @@ def plot_bqnu_with_phbands(self, phbands_qpath, phdos_file=None,
         anaddb_kwargs = {} if anaddb_kwargs is None else anaddb_kwargs
         phbands_qmesh = ddb.anaget_phmodes_at_qpoints(ngqpt=ngqpt, verbose=verbose, **anaddb_kwargs)
 
+        # FIXME
+        # NB: This is just to sketch the ideas. I don't think the present version
+        # is correct as only the k --> -k symmetry can be used.
         for iq, qpoint in enumerate(phbands_qmesh.qpoints):
             weight = qpoint.weight
             freqs_nu = phbands_qmesh.phfreqs[iq]
@@ -636,6 +719,21 @@ def plot_bqnu_with_phbands(self, phbands_qpath, phdos_file=None,
             for w, b2 in zip(freqs_nu, b2_nu):
                 bqnu_dos += weight * b2 * gaussian(mesh, width, center=w)
 
+        # TODO New version using the BZ. Requires new VARPEQ.nc file with all symmetries
+        # The B_qnu do not necessarily have the symmetry of the lattice so we have to loop over the full BZ.
+        # Get mapping BZ --> IBZ needed to obtain the KS eigenvalues e_nk from the IBZ for the DOS
+        """
+        shifts = [0, 0, 0]
+        bz_qpoints = kmesh_from_mpdivs(ngqpt, shifts)
+        bz2ibz = map_grid2ibz(self.structure, phbands_qmesh.qpoints.frac_coords, ngqpt, has_timrev)
+        for iq_ibz, qpoint in zip(bz2ibz, bz_qpoints):
+            freqs_nu = phbands_qmesh.phfreqs[iq_ibz]
+            b2_nu = b2_interp.eval_kpoint(qpoint)
+            for w, b2 in zip(freqs_nu, b2_nu):
+                bqnu_dos += weight * b2 * gaussian(mesh, width, center=w)
+        bqnu_dos /= np.product(ngqpt)
+        """
+
         bqnu_dos = Function1D(mesh, bqnu_dos)
 
         ax = ax_list[1]
diff --git a/abipy/iotools/xsf.py b/abipy/iotools/xsf.py
index 845bf5fbc..c90579157 100644
--- a/abipy/iotools/xsf.py
+++ b/abipy/iotools/xsf.py
@@ -170,7 +170,7 @@ def xsf_write_data(file, structure, data, add_replicas=True, cplx_mode=None,
     fwrite('END_BLOCK_DATAGRID_3D\n')
 
 
-def bxsf_write(file, structure, nsppol, nband, ndivs, ucdata_sbk, fermie, unit="eV") -> None:
+def bxsf_write(file, structure, nsppol, nband, ngkpt, ucdata_sbk, fermie, unit="eV") -> None:
     """
     Write band structure data in the Xcrysden format (XSF)
 
@@ -179,9 +179,9 @@ def bxsf_write(file, structure, nsppol, nband, ndivs, ucdata_sbk, fermie, unit="
         structure: :class:`Structure` object.
         nsppol: Number of spins.
         nband: Number of bands.
-        ndivs: Number of divisions of the full k-mesh.
-        ucdata_sbk: Array [nsppol, nband, ndivs[0], ndivs[1], mpdvis[2]] with energies
-            in the unic cell mesh in unit `unit`.
+        ngkpt: Number of divisions of the full k-mesh.
+        ucdata_sbk: Array [nsppol, nband, ngkpt[0], ngkpt[1], ngkpt[2]] with energies
+            in the unit cell mesh in unit `unit`.
         fermie: Fermi energy.
         unit=Unit of input `ucdata_sbk` and `fermie`. Energies will be converted to Hartree before writing.
 
@@ -201,11 +201,11 @@ def bxsf_write(file, structure, nsppol, nband, ndivs, ucdata_sbk, fermie, unit="
     ucdata_sbk = EnergyArray(ucdata_sbk, unit).to("Ha")
     fermie = Energy(fermie, unit).to("Ha")
 
-    ucdata_sbk = np.reshape(ucdata_sbk, (nsppol, nband, np.prod(ndivs)))
+    ucdata_sbk = np.reshape(ucdata_sbk, (nsppol, nband, np.prod(ngkpt)))
 
     close_it = False
     if not hasattr(file, "write"):
-        file = open(file, mode="w")
+        file = open(str(file), mode="wt")
         close_it = True
 
     fw = file.write
@@ -223,8 +223,8 @@ def bxsf_write(file, structure, nsppol, nband, ndivs, ucdata_sbk, fermie, unit="
     fw(' BEGIN_BANDGRID_3D\n')
 
     fw(str(nsppol * nband) + "\n")  # Number of bands written.
-    fw("%d %d %d\n" % tuple(ndivs)) # Number of division in the full BZ mesh.
-    fw("0 0 0\n")                   # Unshifted meshes are not supported.
+    fw("%d %d %d\n" % tuple(ngkpt)) # Number of division in the full BZ mesh.
+    fw("0 0 0\n")                   # NB: Unshifted meshes are not supported.
 
     # Reciprocal lattice vectors in Ang^{-1}
     gcell = structure.lattice.reciprocal_lattice.matrix
diff --git a/abipy/lumi/tests/test_lineshape.py b/abipy/lumi/tests/test_lineshape.py
index 71a362c92..2d4578915 100644
--- a/abipy/lumi/tests/test_lineshape.py
+++ b/abipy/lumi/tests/test_lineshape.py
@@ -29,7 +29,7 @@ def test_deltaSCF(self):
                 [3,3,3], # test same size dSCF scell than phonon scell
                 [4,4,4], # test same size dSCF scell smaller phonon scell
             ]
-        qpts_s=[kmesh_from_mpdivs(mpdivs=qgrid,shifts=[0,0,0],order="unit_cell") for qgrid in qgrids]
+        qpts_s=[kmesh_from_mpdivs(mpdivs=qgrid, shifts=[0,0,0], order="unit_cell") for qgrid in qgrids]
 
         ddb_pristine_s=[ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts) for qpts in qpts_s]
 
diff --git a/abipy/ml/aseml.py b/abipy/ml/aseml.py
index f0ffaee34..f50f1dd2f 100644
--- a/abipy/ml/aseml.py
+++ b/abipy/ml/aseml.py
@@ -1535,7 +1535,15 @@ def get_calculator(self, with_delta: bool = True, reset: bool = False) -> Calcul
         #if reset: self.reset()
         self.reset()
 
-        if self.nn_type == "m3gnet":
+        if self.nn_type == "emt":
+            # Set the calculator (EMT for testing purposes)
+            from ase.calculators.emt import EMT
+            class MyEMTCalculator(_MyCalculator, EMTCalculator):
+                """Add abi_forces and abi_stress"""
+
+            return MyEMTCalculator(**self.calc_kwargs)
+
+        elif self.nn_type == "m3gnet":
             # m3gnet legacy version.
             if self._model is None:
                 silence_tensorflow()
diff --git a/abipy/scripts/abiview.py b/abipy/scripts/abiview.py
index a56df820e..7f1a7694d 100755
--- a/abipy/scripts/abiview.py
+++ b/abipy/scripts/abiview.py
@@ -528,6 +528,26 @@ def abiview_lobster(options) -> int:
     return 0
 
 
+def abiview_xrd_traj(options) -> int:
+    """
+    Compare XRD spectra using the first and the last structure read from a trajectory file.
+    """
+    from abipy.core.structure import get_first_and_last_structure_from_file
+    structures = get_first_and_last_structure_from_file(options.filepath)
+
+    dfs = abilab.dataframes_from_structures(structures, index=["first", "last"])
+    abilab.print_dataframe(dfs.lattice, title="Lattice parameters:")
+    if options.verbose:
+        abilab.print_dataframe(dfs.coords, title="Atomic positions (columns give the site index):")
+
+    from pymatgen.analysis.diffraction.xrd import XRDCalculator
+    xrd = XRDCalculator(wavelength=options.wavelength, symprec=options.symprec)
+    two_theta_range = tuple(float(t) for t in options.two_theta_range)
+    xrd.plot_structures(structures, two_theta_range=two_theta_range, fontsize=6,
+                        annotate_peaks=not options.no_annotate_peaks, tight_layout=True)
+    return 0
+
+
 def get_epilog() -> str:
     return """\
 Usage example:
@@ -797,6 +817,21 @@ def add_args(p, *args):
         help=abiview_lobster.__doc__)
     p_lobster.add_argument("--prefix", type=str, default="", help="Prefix for lobster output files. Default: ''")
 
+    # Subparser for xrd.
+    p_xrd = subparsers.add_parser('xrd_traj', parents=[copts_parser],
+        help="Compare X-ray diffraction for the first and the last structure in a trajectory file.")
+    p_xrd.add_argument("-w", "--wavelength", default="CuKa", type=str, help=(
+        "The wavelength can be specified as a string. It must be one of the "
+        "supported definitions in the WAVELENGTHS dict declared in pymatgen/analysis/diffraction/xrd.py."
+        "Defaults to 'CuKa', i.e, Cu K_alpha radiation."))
+    p_xrd.add_argument("-s", "--symprec", default=0, type=float, help=(
+        "Symmetry precision for structure refinement. "
+        "If set to 0, no refinement is done. Otherwise, refinement is performed using spglib with provided precision."))
+    p_xrd.add_argument("-t", "--two-theta-range", default=(0, 90), nargs=2, help=(
+        "Tuple for range of two_thetas to calculate in degrees. Defaults to (0, 90)."))
+    p_xrd.add_argument("-nap", "--no-annotate-peaks", default=False, action="store_true",
+        help="Whether to annotate the peaks with plane information.")
+
     # Subparser for denpot command.
     p_denpot = subparsers.add_parser('denpot', parents=[copts_parser], help=abiview_denpot.__doc__)
     p_denpot.add_argument("-a", "--appname", type=str, default="vesta",