Compute A2(E) and B2(E) using full BZ instead of IBZ

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):
  """

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.
 

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():

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):

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``).