diff --git a/CHANGELOG.md b/CHANGELOG.md index f0de9d7c2..c765bf1fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ we hit release version 1.0.0. ## [0.15.0] - YYYY-MM-DD ### Added +- added `vectorsSileSiesta` to read vibra eigenmode output - added `dihedral` to `Geometry`, #773 - ability to retain sub-classes through `.new` calls - added `Listify` to ensure arguments behaves as *iterables* diff --git a/src/sisl/io/siesta/binaries.py b/src/sisl/io/siesta/binaries.py index 2f1c1a57b..4a63b412c 100644 --- a/src/sisl/io/siesta/binaries.py +++ b/src/sisl/io/siesta/binaries.py @@ -1919,7 +1919,16 @@ def yield_eigenstate(self): # The loop in which the generator was used has been broken. self._close_wfsx() - def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron: + @deprecate_argument( + "ktol", + "atol", + "use atol instead of ktol", + "0.15", + "0.16", + ) + def read_eigenstate( + self, k=(0, 0, 0), spin: int = 0, atol: float = 1e-4 + ) -> EigenstateElectron: """Reads a specific eigenstate from the file. This method iterates over the states until it finds a match. Do not call @@ -1929,10 +1938,10 @@ def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron: ---------- k: array-like of shape (3,), optional The k point of the state you want to find. - spin: integer, optional + spin: The spin index of the state you want to find. Only meaningful for polarized calculations. - ktol: float, optional + atol: The threshold value for considering two k-points the same (i.e. to match the query k point with the states k point). @@ -1944,16 +1953,21 @@ def read_eigenstate(self, k=(0, 0, 0), spin=0, ktol=1e-4) -> EigenstateElectron: ------- EigenstateElectron or None: If found, the state that was queried. - If not found, returns `None`. NOTE this may change to an exception in the future + + Raises + ------ + LookupError : + in case the requested k-point can not be found in the file. """ # Iterate over all eigenstates in the file for state in self.yield_eigenstate(): - if state.info.get("spin", 0) == spin and np.allclose( - state.info["k"], k, atol=ktol - ): + info = state.info + if info.get("spin", 0) == spin and np.allclose(info["k"], k, atol=atol): # This is the state that the user requested return state - return None + raise LookupError( + f"{self.__class__.__name__}.read_eigenstate could not find k-point: {k!s} eigenstate" + ) def read_info(self): """Reads the information for all the k points contained in the file diff --git a/src/sisl/io/siesta/tests/test_vibra.py b/src/sisl/io/siesta/tests/test_vibra.py index fba6da026..5b372a9bd 100644 --- a/src/sisl/io/siesta/tests/test_vibra.py +++ b/src/sisl/io/siesta/tests/test_vibra.py @@ -5,6 +5,7 @@ import os.path as osp +import numpy as np import pytest import sisl @@ -26,15 +27,31 @@ def test_eigenmode_read(sisl_files): nlines = len(blines) iklast = [sum(int(l.split()[0]) for l in blines[: i + 1]) for i in range(nlines)] klast = [list(map(float, l.split()[1:4])) for l in blines] - print(iklast, klast) # yield modes nmodes = 0 nmodes_total = 0 + ks = [] for state in vectors.yield_eigenmode(): nmodes += 1 nmodes_total += len(state) + ks.append(state.info["k"]) - print(nmodes, nmodes_total) assert nmodes == nk assert nmodes_total == geometry.na * 3 * nk + + bz = vectors.read_brillouinzone() + assert len(bz) == nk + assert np.allclose(bz.k, ks) + + +def test_eigenmode_values(sisl_files): + fdf = sisl.get_sile(sisl_files(_dir, "h_chain_vibra.fdf")) + vectors = sisl.io.siesta.vectorsSileSiesta( + sisl_files(_dir, "h_chain_vibra.vectors"), + geometry=fdf.read_geometry(), + ) + + mode = vectors.read_eigenmode() + assert np.allclose(mode.state[0, 0:3], [-0.7071, -0.2400e-4, -0.2400e-4]) + assert np.allclose(mode.state[0, 3:7], [0.7071, 0.2400e-4, 0.2400e-4]) diff --git a/src/sisl/io/siesta/vibra.py b/src/sisl/io/siesta/vibra.py index 56fc7fd29..e782e9a04 100644 --- a/src/sisl/io/siesta/vibra.py +++ b/src/sisl/io/siesta/vibra.py @@ -4,6 +4,7 @@ from __future__ import annotations import re +from typing import Iterator import numpy as np @@ -11,6 +12,7 @@ from sisl import Geometry, Lattice from sisl._internal import set_module from sisl.messages import warn +from sisl.physics.brillouinzone import BrillouinZone from sisl.physics.phonon import EigenmodePhonon from sisl.unit.siesta import unit_convert @@ -30,11 +32,10 @@ class vectorsSileSiesta(SileSiesta): Parameters ---------- parent : obj, optional - a parent may contain a DynamicalMatrix, or Geometry + a parent may contain a DynamicalMatrix, Geometry or Lattice geometry : Geometry, optional a geometry contains a cell with corresponding lattice vectors used to convert k [1/Ang] -> [b], and the atomic masses - """ def _setup(self, *args, **kwargs): @@ -62,7 +63,7 @@ def _setup(self, *args, **kwargs): if geometry is not None and lattice is None: lattice = geometry.lattice - lattice = kwargs.get("lattice", kwargs.get("sc", lattice)) + lattice = kwargs.get("lattice", lattice) if lattice is None and geometry is not None: raise ValueError( f"{self.__class__.__name__}(geometry=Geometry, lattice=None) is not an allowed argument combination." @@ -82,7 +83,7 @@ def conv(k): if not np.allclose(k, 0.0): warn( f"{self.__class__.__name__} cannot convert stored k-points from 1/Ang to reduced coordinates. " - "Please ensure parent=DynamicalMatrix, geometry=Geometry, or lattice=Lattice to ensure reduced k." + "Please ensure parent=DynamicalMatrix, geometry=Geometry, or lattice=Lattice to calculate reduced k." ) return k / _Bohr2Ang @@ -93,7 +94,7 @@ def conv(k): self._convert_k = conv - def _open(self, rewind=False): + def _open(self, rewind: bool = False): """Open the file Here we also initialize some variables to keep track of the state of the read. @@ -101,6 +102,7 @@ def _open(self, rewind=False): super()._open() if rewind: self.seek(0) + # Here we initialize the variables that will keep track of the state of the read. # The process for identification is done on this basis: # _ik is the current (Python) index for the k-point to be read @@ -110,27 +112,14 @@ def _open(self, rewind=False): self._state = -1 self._ik = 0 - def _setup_parsing(self, close=True): - """Gets all the things needed to parse the wfsx file. - - Parameters - ----------- - close: bool, optional - Whether the file unit should be closed afterwards. - """ - self._open() - # Read the sizes relevant to the file. - # We also read whether there's only gamma point information or there are multiple points - if self._state == -1: - self._sizes = self._r_next_sizes() - self._state = 0 - - if close: - self._close() + # read in sizes + self._sizes = self._r_next_sizes() + self._state = 0 - def _r_next_sizes(self): + def _r_next_sizes(self) -> None: """Determine the dimension of the data stored in the vectors file""" pos = self.fh.tell() + # Skip past header self.readline() # Empty line self.readline() # K point @@ -161,7 +150,7 @@ def _r_next_sizes(self): # Rewind self.seek(pos) - def _r_next_eigenmode(self): + def _r_next_eigenmode(self) -> EigenmodePhonon: """Reads the next phonon eigenmode in the vectors file. Returns @@ -172,14 +161,14 @@ def _r_next_eigenmode(self): # Skip empty line at the head of the file self.readline() - k = _a.asarrayd(list(map(float, self.readline().split()[2:]))) + k = _a.fromiterd(map(float, self.readline().split()[2:])) if len(k) == 0: - raise GeneratorExit + return None # Read first eigenvector index # Determine number of atoms (i.e. rows per mode) - mode = np.empty((self._nmodes, 3, self._natoms), dtype=np.complex128) - c = np.empty(self._nmodes, dtype=np.float64) + mode = _a.emptyz([self._nmodes, self._natoms, 3]) + c = _a.emptyd(self._nmodes) for imode in range(self._nmodes): self.readline() @@ -191,19 +180,19 @@ def _r_next_eigenmode(self): self.readline() for iatom in range(self._natoms): line = self.readline() - mode[imode, :, iatom].real = list(map(float, line.split())) + mode[imode, iatom].real = list(map(float, line.split())) # Read imaginary part of eigenmode # Skip eigenmode header self.readline() for iatom in range(self._natoms): line = self.readline() - mode[imode, :, iatom].imag = list(map(float, line.split())) + mode[imode, iatom].imag = list(map(float, line.split())) info = dict(k=self._convert_k(k), parent=self._parent, gauge="orbital") return EigenmodePhonon(mode.reshape(self._nmodes, -1), c * _cm1_eV, **info) - def yield_eigenmode(self): + def yield_eigenmode(self) -> Iterator[EigenmodePhonon]: """Iterates over the modes in the vectors file Yields @@ -211,19 +200,16 @@ def yield_eigenmode(self): EigenmodePhonon """ # Open file and get parsing information - self._setup_parsing(close=False) - - try: - # Iterate over all eigenmodes in the WFSX file, yielding control to the caller at - # each iteration. - while True: - yield self._r_next_eigenmode() - # We ran out of eigenmodes - except GeneratorExit: - # The loop in which the generator was used has been broken. - pass - - def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon: + self._open() + + # Iterate over all eigenmodes in the vectors file, yielding control to the caller at + # each iteration. + mode = self._r_next_eigenmode() + while mode is not None: + yield mode + mode = self._r_next_eigenmode() + + def read_eigenmode(self, k=(0, 0, 0), atol: float = 1e-4) -> EigenmodePhonon: """Reads a specific eigenmode from the file. This method iterates over the modes until it finds a match. Do not call @@ -233,7 +219,7 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon: ---------- k: array-like of shape (3,), optional The k point of the mode you want to find. - ktol: + atol: The threshold value for considering two k-points the same (i.e. to match the query k point with the modes k point). @@ -245,14 +231,25 @@ def read_eigenmode(self, k=(0, 0, 0), ktol: float = 1e-4) -> EigenmodePhonon: ------- EigenmodePhonon or None: If found, the mode that was queried. - If not found, returns `None`. NOTE this may change to an exception in the future + + Raises + ------ + LookupError : + in case the requested k-point can not be found in the file. """ # Iterate over all eigenmodes in the file for mode in self.yield_eigenmode(): - if np.allclose(mode.info["k"], k, atol=ktol): + if np.allclose(mode.info["k"], k, atol=atol): # This is the mode that the user requested return mode - return None + raise LookupError( + f"{self.__class__.__name__}.read_eigenmode could not find k-point: {k!s} eigenmode" + ) + + def read_brillouinzone(self) -> BrillouinZone: + """Read the brillouin zone object (only containing the k-points)""" + k = [mode.info["k"] for mode in self.yield_eigenmode()] + return BrillouinZone(self._parent, k=k) add_sile("vectors", vectorsSileSiesta, gzip=True)