From 65d537933526aa5f41bf68df4c66bad928df272b Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Thu, 6 Jun 2024 19:16:24 +0800 Subject: [PATCH] Add types annotations for `core.interface` (#3822) * tweak type and docstring move dunder methods to the top add more types and tweaks relocate more dunder methods to top more types and format tweaks fix type error add types for composition help fix https://github.com/materialsproject/pymatgen/pull/3792#discussion_r1584003750 reverse compare order for readability Revert "reverse compare order for readability" This reverts commit 05ea23a35ffaf03dc70e9dcbdbf302ac69765603. Revert "help fix https://github.com/materialsproject/pymatgen/pull/3792#discussion_r1584003750" This reverts commit cae7aed4464f4f54c1b2f76d4f12c7a07a49d04e. add types for `core.bonds` finish `core.ion` add some types revert non-interface changes * use math.gcd over gcd * use more specific types for ClassVar * more: use more specific types for ClassVar * fix some type errors and comment tweaks * fix mypy errors * enable types: more type errors to fix * fix type errors * fix type errors * fix mypy errors doesn't show locally * revert change in test and convert rotation_axis and plane to tuple * cast plane type to tuple * remove `del` of var name * add and update new type `Tuple3Ints = tuple[int, int, int]` * relocate `Tuple4Ints` to `core.interface` * relocate `Tuple4Ints` * use `Tuple3Floats` * revert usage of assert_allclose * use more meaningful types * fix replacement * Revert "fix replacement" This reverts commit 6b9589ad681c1e53ab62cf5763cb090f94133541. * revert type aliases in docstring --- pymatgen/analysis/diffraction/tem.py | 57 +- pymatgen/analysis/graphs.py | 9 +- .../interfaces/coherent_interfaces.py | 5 +- .../analysis/interfaces/substrate_analyzer.py | 5 +- pymatgen/analysis/local_env.py | 3 +- pymatgen/analysis/surface_analysis.py | 4 +- pymatgen/core/composition.py | 2 +- pymatgen/core/interface.py | 1205 +++++++++-------- pymatgen/core/surface.py | 13 +- pymatgen/core/xcfunc.py | 10 +- pymatgen/electronic_structure/dos.py | 6 +- pymatgen/ext/optimade.py | 2 +- pymatgen/io/abinit/abiobjects.py | 32 +- pymatgen/io/abinit/pseudos.py | 8 +- pymatgen/io/adf.py | 2 +- pymatgen/io/aims/inputs.py | 6 +- pymatgen/io/aims/parsers.py | 5 +- pymatgen/io/cp2k/inputs.py | 4 +- pymatgen/io/fiesta.py | 6 +- pymatgen/io/lobster/inputs.py | 3 +- pymatgen/io/nwchem.py | 4 +- pymatgen/io/pwmat/outputs.py | 4 +- pymatgen/io/pwscf.py | 2 +- pymatgen/io/res.py | 4 +- pymatgen/io/vasp/inputs.py | 18 +- pymatgen/io/vasp/outputs.py | 9 +- pymatgen/io/vasp/sets.py | 6 +- pymatgen/phonon/bandstructure.py | 4 +- pymatgen/symmetry/groups.py | 8 +- pymatgen/util/typing.py | 9 +- pymatgen/vis/structure_vtk.py | 2 +- tests/core/test_interface.py | 6 +- 32 files changed, 796 insertions(+), 667 deletions(-) diff --git a/pymatgen/analysis/diffraction/tem.py b/pymatgen/analysis/diffraction/tem.py index 3cf05361a95..7c7a7584102 100644 --- a/pymatgen/analysis/diffraction/tem.py +++ b/pymatgen/analysis/diffraction/tem.py @@ -15,6 +15,7 @@ from pymatgen.analysis.diffraction.core import AbstractDiffractionPatternCalculator from pymatgen.symmetry.analyzer import SpacegroupAnalyzer from pymatgen.util.string import latexify_spacegroup, unicodeify_spacegroup +from pymatgen.util.typing import Tuple3Ints if TYPE_CHECKING: from numpy.typing import NDArray @@ -47,7 +48,7 @@ def __init__( self, symprec: float | None = None, voltage: float = 200, - beam_direction: tuple[int, int, int] = (0, 0, 1), + beam_direction: Tuple3Ints = (0, 0, 1), camera_length: int = 160, debye_waller_factors: dict[str, float] | None = None, cs: float = 1, @@ -104,9 +105,7 @@ def generate_points(coord_left: int = -10, coord_right: int = 10) -> np.ndarray: points_matrix = (np.ravel(points[i]) for i in range(3)) return np.vstack(list(points_matrix)).transpose() - def zone_axis_filter( - self, points: list[tuple[int, int, int]] | np.ndarray, laue_zone: int = 0 - ) -> list[tuple[int, int, int]]: + def zone_axis_filter(self, points: list[Tuple3Ints] | np.ndarray, laue_zone: int = 0) -> list[Tuple3Ints]: """Filter out all points that exist within the specified Laue zone according to the zone axis rule. Args: @@ -122,11 +121,11 @@ def zone_axis_filter( return [] filtered = np.where(np.dot(np.array(self.beam_direction), np.transpose(points)) == laue_zone) result = points[filtered] # type: ignore - return cast(list[tuple[int, int, int]], [tuple(x) for x in result.tolist()]) + return cast(list[Tuple3Ints], [tuple(x) for x in result.tolist()]) def get_interplanar_spacings( - self, structure: Structure, points: list[tuple[int, int, int]] | np.ndarray - ) -> dict[tuple[int, int, int], float]: + self, structure: Structure, points: list[Tuple3Ints] | np.ndarray + ) -> dict[Tuple3Ints, float]: """ Args: structure (Structure): the input structure. @@ -142,9 +141,7 @@ def get_interplanar_spacings( interplanar_spacings_val = np.array([structure.lattice.d_hkl(x) for x in points_filtered]) return dict(zip(points_filtered, interplanar_spacings_val)) - def bragg_angles( - self, interplanar_spacings: dict[tuple[int, int, int], float] - ) -> dict[tuple[int, int, int], float]: + def bragg_angles(self, interplanar_spacings: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, float]: """Get the Bragg angles for every hkl point passed in (where n = 1). Args: @@ -158,7 +155,7 @@ def bragg_angles( bragg_angles_val = np.arcsin(self.wavelength_rel() / (2 * interplanar_spacings_val)) return dict(zip(plane, bragg_angles_val)) - def get_s2(self, bragg_angles: dict[tuple[int, int, int], float]) -> dict[tuple[int, int, int], float]: + def get_s2(self, bragg_angles: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, float]: """ Calculates the s squared parameter (= square of sin theta over lambda) for each hkl plane. @@ -175,8 +172,8 @@ def get_s2(self, bragg_angles: dict[tuple[int, int, int], float]) -> dict[tuple[ return dict(zip(plane, s2_val)) def x_ray_factors( - self, structure: Structure, bragg_angles: dict[tuple[int, int, int], float] - ) -> dict[str, dict[tuple[int, int, int], float]]: + self, structure: Structure, bragg_angles: dict[Tuple3Ints, float] + ) -> dict[str, dict[Tuple3Ints, float]]: """ Calculates x-ray factors, which are required to calculate atomic scattering factors. Method partially inspired by the equivalent process in the xrd module. @@ -205,8 +202,8 @@ def x_ray_factors( return x_ray_factors def electron_scattering_factors( - self, structure: Structure, bragg_angles: dict[tuple[int, int, int], float] - ) -> dict[str, dict[tuple[int, int, int], float]]: + self, structure: Structure, bragg_angles: dict[Tuple3Ints, float] + ) -> dict[str, dict[Tuple3Ints, float]]: """ Calculates atomic scattering factors for electrons using the Mott-Bethe formula (1st order Born approximation). @@ -232,8 +229,8 @@ def electron_scattering_factors( return electron_scattering_factors def cell_scattering_factors( - self, structure: Structure, bragg_angles: dict[tuple[int, int, int], float] - ) -> dict[tuple[int, int, int], int]: + self, structure: Structure, bragg_angles: dict[Tuple3Ints, float] + ) -> dict[Tuple3Ints, int]: """ Calculates the scattering factor for the whole cell. @@ -258,9 +255,7 @@ def cell_scattering_factors( scattering_factor_curr = 0 return cell_scattering_factors - def cell_intensity( - self, structure: Structure, bragg_angles: dict[tuple[int, int, int], float] - ) -> dict[tuple[int, int, int], float]: + def cell_intensity(self, structure: Structure, bragg_angles: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, float]: """ Calculates cell intensity for each hkl plane. For simplicity's sake, take I = |F|**2. @@ -317,8 +312,8 @@ def get_pattern( return pd.DataFrame(rows, columns=field_names) def normalized_cell_intensity( - self, structure: Structure, bragg_angles: dict[tuple[int, int, int], float] - ) -> dict[tuple[int, int, int], float]: + self, structure: Structure, bragg_angles: dict[Tuple3Ints, float] + ) -> dict[Tuple3Ints, float]: """ Normalizes the cell_intensity dict to 1, for use in plotting. @@ -340,8 +335,8 @@ def normalized_cell_intensity( def is_parallel( self, structure: Structure, - plane: tuple[int, int, int], - other_plane: tuple[int, int, int], + plane: Tuple3Ints, + other_plane: Tuple3Ints, ) -> bool: """ Checks if two hkl planes are parallel in reciprocal space. @@ -357,7 +352,7 @@ def is_parallel( phi = self.get_interplanar_angle(structure, plane, other_plane) return phi in (180, 0) or np.isnan(phi) - def get_first_point(self, structure: Structure, points: list) -> dict[tuple[int, int, int], float]: + def get_first_point(self, structure: Structure, points: list) -> dict[Tuple3Ints, float]: """Get the first point to be plotted in the 2D DP, corresponding to maximum d/minimum R. Args: @@ -378,7 +373,7 @@ def get_first_point(self, structure: Structure, points: list) -> dict[tuple[int, return {max_d_plane: max_d} @staticmethod - def get_interplanar_angle(structure: Structure, p1: tuple[int, int, int], p2: tuple[int, int, int]) -> float: + def get_interplanar_angle(structure: Structure, p1: Tuple3Ints, p2: Tuple3Ints) -> float: """Get the interplanar angle (in degrees) between the normal of two crystal planes. Formulas from International Tables for Crystallography Volume C pp. 2-9. @@ -432,9 +427,9 @@ def get_interplanar_angle(structure: Structure, p1: tuple[int, int, int], p2: tu @staticmethod def get_plot_coeffs( - p1: tuple[int, int, int], - p2: tuple[int, int, int], - p3: tuple[int, int, int], + p1: Tuple3Ints, + p2: Tuple3Ints, + p3: Tuple3Ints, ) -> np.ndarray: """ Calculates coefficients of the vector addition required to generate positions for each DP point @@ -454,7 +449,7 @@ def get_plot_coeffs( x = np.dot(a_pinv, b) return np.ravel(x) - def get_positions(self, structure: Structure, points: list) -> dict[tuple[int, int, int], np.ndarray]: + def get_positions(self, structure: Structure, points: list) -> dict[Tuple3Ints, np.ndarray]: """ Calculates all the positions of each hkl point in the 2D diffraction pattern by vector addition. Distance in centimeters. @@ -524,7 +519,7 @@ def tem_dots(self, structure: Structure, points) -> list: class dot(NamedTuple): position: NDArray - hkl: tuple[int, int, int] + hkl: Tuple3Ints intensity: float film_radius: float d_spacing: float diff --git a/pymatgen/analysis/graphs.py b/pymatgen/analysis/graphs.py index f760f4239b6..7f403a7b5c4 100644 --- a/pymatgen/analysis/graphs.py +++ b/pymatgen/analysis/graphs.py @@ -43,6 +43,7 @@ from pymatgen.analysis.local_env import NearNeighbors from pymatgen.core import Species + from pymatgen.util.typing import Tuple3Ints logger = logging.getLogger(__name__) @@ -58,7 +59,7 @@ class ConnectedSite(NamedTuple): site: PeriodicSite - jimage: tuple[int, int, int] + jimage: Tuple3Ints index: Any # TODO: use more specific type weight: float dist: float @@ -338,8 +339,8 @@ def add_edge( self, from_index: int, to_index: int, - from_jimage: tuple[int, int, int] = (0, 0, 0), - to_jimage: tuple[int, int, int] | None = None, + from_jimage: Tuple3Ints = (0, 0, 0), + to_jimage: Tuple3Ints | None = None, weight: float | None = None, warn_duplicates: bool = True, edge_properties: dict | None = None, @@ -756,7 +757,7 @@ def map_indices(grp: Molecule) -> dict[int, int]: warn_duplicates=False, ) - def get_connected_sites(self, n: int, jimage: tuple[int, int, int] = (0, 0, 0)) -> list[ConnectedSite]: + def get_connected_sites(self, n: int, jimage: Tuple3Ints = (0, 0, 0)) -> list[ConnectedSite]: """Get a named tuple of neighbors of site n: periodic_site, jimage, index, weight. Index is the index of the corresponding site diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index c7fb9b08736..3e145fc5618 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -18,6 +18,7 @@ from collections.abc import Iterator, Sequence from pymatgen.core import Structure + from pymatgen.util.typing import Tuple3Ints class CoherentInterfaceBuilder: @@ -30,8 +31,8 @@ def __init__( self, substrate_structure: Structure, film_structure: Structure, - film_miller: tuple[int, int, int], - substrate_miller: tuple[int, int, int], + film_miller: Tuple3Ints, + substrate_miller: Tuple3Ints, zslgen: ZSLGenerator | None = None, ): """ diff --git a/pymatgen/analysis/interfaces/substrate_analyzer.py b/pymatgen/analysis/interfaces/substrate_analyzer.py index 49d816861e9..4d6f3ecf431 100644 --- a/pymatgen/analysis/interfaces/substrate_analyzer.py +++ b/pymatgen/analysis/interfaces/substrate_analyzer.py @@ -14,6 +14,7 @@ from typing_extensions import Self from pymatgen.core import Structure + from pymatgen.util.typing import Tuple3Ints @dataclass @@ -24,8 +25,8 @@ class SubstrateMatch(ZSLMatch): energy if provided, and the elastic energy. """ - film_miller: tuple[int, int, int] - substrate_miller: tuple[int, int, int] + film_miller: Tuple3Ints + substrate_miller: Tuple3Ints strain: Strain von_mises_strain: float ground_state_energy: float diff --git a/pymatgen/analysis/local_env.py b/pymatgen/analysis/local_env.py index 8c571a6ae61..26fbc93755c 100644 --- a/pymatgen/analysis/local_env.py +++ b/pymatgen/analysis/local_env.py @@ -38,6 +38,7 @@ from typing_extensions import Self from pymatgen.core.composition import SpeciesLike + from pymatgen.util.typing import Tuple3Ints __author__ = "Shyue Ping Ong, Geoffroy Hautier, Sai Jayaraman, " @@ -540,7 +541,7 @@ def _get_nn_shell_info( return list(all_sites.values()) @staticmethod - def _get_image(structure: Structure, site: Site) -> tuple[int, int, int]: + def _get_image(structure: Structure, site: Site) -> Tuple3Ints: """Private convenience method for get_nn_info, gives lattice image from provided PeriodicSite and Structure. diff --git a/pymatgen/analysis/surface_analysis.py b/pymatgen/analysis/surface_analysis.py index d114cb93bd3..c251e44c479 100644 --- a/pymatgen/analysis/surface_analysis.py +++ b/pymatgen/analysis/surface_analysis.py @@ -58,6 +58,8 @@ if TYPE_CHECKING: from typing_extensions import Self + from pymatgen.util.typing import Tuple3Ints + EV_PER_ANG2_TO_JOULES_PER_M2 = 16.0217656 __author__ = "Richard Tran" @@ -566,7 +568,7 @@ def area_frac_vs_chempot_plot( all_chempots = np.linspace(min(chempot_range), max(chempot_range), increments) # initialize a dictionary of lists of fractional areas for each hkl - hkl_area_dict: dict[tuple[int, int, int], list[float]] = {} + hkl_area_dict: dict[Tuple3Ints, list[float]] = {} for hkl in self.all_slab_entries: hkl_area_dict[hkl] = [] diff --git a/pymatgen/core/composition.py b/pymatgen/core/composition.py index c8446dd71fc..b15775eeb30 100644 --- a/pymatgen/core/composition.py +++ b/pymatgen/core/composition.py @@ -81,7 +81,7 @@ class Composition(collections.abc.Hashable, collections.abc.Mapping, MSONable, S # Special formula handling for peroxides and certain elements. This is so # that formula output does not write LiO instead of Li2O2 for example. - special_formulas: ClassVar = dict( + special_formulas: ClassVar[dict[str, str]] = dict( LiO="Li2O2", NaO="Na2O2", KO="K2O2", diff --git a/pymatgen/core/interface.py b/pymatgen/core/interface.py index 35e3cfa94be..d7d5c8eccf3 100644 --- a/pymatgen/core/interface.py +++ b/pymatgen/core/interface.py @@ -1,14 +1,16 @@ -"""This module provides classes to store, generate, and manipulate material interfaces, including grain boundaries.""" +"""This module provides classes to store, generate, +and manipulate interfaces, including grain boundaries. +""" from __future__ import annotations import logging +import math import warnings from fractions import Fraction from functools import reduce from itertools import chain, combinations, product -from math import cos, floor, gcd -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal, Union, cast import numpy as np from monty.fractions import lcm @@ -22,18 +24,21 @@ from pymatgen.core.structure import Structure from pymatgen.core.surface import Slab from pymatgen.symmetry.analyzer import SpacegroupAnalyzer +from pymatgen.util.typing import Tuple3Ints if TYPE_CHECKING: from collections.abc import Sequence - from typing import Any + from typing import Any, Callable - from numpy.typing import ArrayLike + from numpy.typing import ArrayLike, NDArray from typing_extensions import Self - from pymatgen.util.typing import CompositionLike, Vector3D + from pymatgen.core import Element + from pymatgen.util.typing import CompositionLike, Matrix3D, MillerIndex, Tuple3Floats, Vector3D + +Tuple4Ints = tuple[int, int, int, int] +logger = logging.getLogger(__name__) -# This module implements representations of grain boundaries, as well as -# algorithms for generating them. __author__ = "Xiang-Guo Li" __copyright__ = "Copyright 2018, The Materials Virtual Lab" @@ -42,18 +47,16 @@ __email__ = "xil110@ucsd.edu" __date__ = "7/30/18" -logger = logging.getLogger(__name__) - class GrainBoundary(Structure): """ - Subclass of Structure representing a GrainBoundary (GB) object. Implements additional - attributes pertaining to gbs, but the init method does not actually implement any + Representation of grain boundary (GB). Implements additional + attributes pertaining to GBs, but the init method does not actually implement any algorithm that creates a GB. This is a DUMMY class who's init method only holds information about the GB. Also has additional methods that returns other information about a GB such as sigma value. - Note that all gbs have the GB surface normal oriented in the c-direction. This means + Note that all GBs have their surface normal oriented in the c-direction. This means the lattice vectors a and b are in the GB surface plane (at least for one grain) and the c vector is out of the surface plane (though not necessarily perpendicular to the surface). @@ -64,10 +67,10 @@ def __init__( lattice: np.ndarray | Lattice, species: Sequence[CompositionLike], coords: Sequence[ArrayLike], - rotation_axis: Vector3D, + rotation_axis: Tuple3Ints | Tuple4Ints, rotation_angle: float, - gb_plane: Vector3D, - join_plane: Vector3D, + gb_plane: Tuple3Ints, + join_plane: Tuple3Ints, init_cell: Structure, vacuum_thickness: float, ab_shift: tuple[float, float], @@ -77,9 +80,7 @@ def __init__( coords_are_cartesian: bool = False, properties: dict | None = None, ) -> None: - """ - Makes a GB structure, a structure object with additional information - and methods pertaining to gbs. + """A Structure with additional information and methods pertaining to GBs. Args: lattice (Lattice/3x3 array): The lattice, either as an instance or @@ -139,9 +140,35 @@ def __init__( properties=properties, ) + def __str__(self) -> str: + comp = self.composition + outs = [ + f"Gb Summary ({comp.formula})", + f"Reduced Formula: {comp.reduced_formula}", + f"Rotation axis: {self.rotation_axis}", + f"Rotation angle: {self.rotation_angle}", + f"GB plane: {self.gb_plane}", + f"Join plane: {self.join_plane}", + f"vacuum thickness: {self.vacuum_thickness}", + f"ab_shift: {self.ab_shift}", + ] + + def to_str(number: float, rjust: int = 10) -> str: + """Convert a float to string and right justify.""" + return (f"{number:0.6f}").rjust(rjust) + + outs += ( + f"abc : {' '.join(to_str(i) for i in self.lattice.abc)}", + f"angles: {' '.join(to_str(i) for i in self.lattice.angles)}", + f"Sites ({len(self)})", + ) + for idx, site in enumerate(self, start=1): + outs.append(f"{idx} {site.species_string} {' '.join(to_str(coord, 12) for coord in site.frac_coords)}") + return "\n".join(outs) + def copy(self) -> Self: # type: ignore[override] - """Make a copy of the GrainBoundary object.""" - return GrainBoundary( + """Make a copy of the GrainBoundary.""" + return type(self)( self.lattice, self.species_and_occu, self.frac_coords, @@ -156,8 +183,12 @@ def copy(self) -> Self: # type: ignore[override] self.oriented_unit_cell, ) - def get_sorted_structure(self, key=None, reverse=False): - """Get a sorted copy of the structure. The parameters have the same + def get_sorted_structure( + self, + key: Callable | None = None, + reverse: bool = False, + ) -> Self: + """Get a sorted copy of the Structure. The parameters have the same meaning as in list.sort. By default, sites are sorted by the electronegativity of the species. Note that Slab has to override this because of the different __init__ args. @@ -171,7 +202,7 @@ def get_sorted_structure(self, key=None, reverse=False): """ sites = sorted(self, key=key, reverse=reverse) struct = Structure.from_sites(sites) - return GrainBoundary( + return type(self)( struct.lattice, struct.species_and_occu, struct.frac_coords, @@ -194,16 +225,14 @@ def sigma(self) -> int: @property def sigma_from_site_prop(self) -> int: """The sigma value of the GB from site properties. - If the GB structure merge some atoms due to the atoms too closer with + If the GB structure merge some atoms due to the atoms too close with each other, this property will not work. """ - n_coi = 0 if None in self.site_properties["grain_label"]: - raise RuntimeError("Site were merged, this property do not work") - for tag in self.site_properties["grain_label"]: - if "incident" in tag: - n_coi += 1 - return int(round(len(self) / n_coi)) + raise ValueError("Sites were merged, this property does not work") + + n_coi = sum("incident" in tag for tag in self.site_properties["grain_label"]) + return round(len(self) / n_coi) @property def top_grain(self) -> Structure: @@ -225,59 +254,35 @@ def bottom_grain(self) -> Structure: @property def coincidents(self) -> list[Site]: - """The a list of coincident sites.""" + """A list of coincident sites.""" coincident_sites = [] for idx, tag in enumerate(self.site_properties["grain_label"]): if "incident" in tag: coincident_sites.append(self.sites[idx]) return coincident_sites - def __str__(self) -> str: - comp = self.composition - outs = [ - f"Gb Summary ({comp.formula})", - f"Reduced Formula: {comp.reduced_formula}", - f"Rotation axis: {self.rotation_axis}", - f"Rotation angle: {self.rotation_angle}", - f"GB plane: {self.gb_plane}", - f"Join plane: {self.join_plane}", - f"vacuum thickness: {self.vacuum_thickness}", - f"ab_shift: {self.ab_shift}", - ] - - def to_str(x, rjust=10): - return (f"{x:0.6f}").rjust(rjust) - - outs += ( - f"abc : {' '.join(to_str(i) for i in self.lattice.abc)}", - f"angles: {' '.join(to_str(i) for i in self.lattice.angles)}", - f"Sites ({len(self)})", - ) - for idx, site in enumerate(self, start=1): - outs.append(f"{idx} {site.species_string} {' '.join(to_str(coord, 12) for coord in site.frac_coords)}") - return "\n".join(outs) - - def as_dict(self): + def as_dict(self) -> dict: # type: ignore[override] """ Returns: Dictionary representation of GrainBoundary object. """ - dct = super().as_dict() - dct["@module"] = type(self).__module__ - dct["@class"] = type(self).__name__ - dct["init_cell"] = self.init_cell.as_dict() - dct["rotation_axis"] = self.rotation_axis - dct["rotation_angle"] = self.rotation_angle - dct["gb_plane"] = self.gb_plane - dct["join_plane"] = self.join_plane - dct["vacuum_thickness"] = self.vacuum_thickness - dct["ab_shift"] = self.ab_shift - dct["oriented_unit_cell"] = self.oriented_unit_cell.as_dict() - return dct + return { + **super().as_dict(), + "@module": type(self).__module__, + "@class": type(self).__name__, + "init_cell": self.init_cell.as_dict(), + "rotation_axis": self.rotation_axis, + "rotation_angle": self.rotation_angle, + "gb_plane": self.gb_plane, + "join_plane": self.join_plane, + "vacuum_thickness": self.vacuum_thickness, + "ab_shift": self.ab_shift, + "oriented_unit_cell": self.oriented_unit_cell.as_dict(), + } @classmethod - def from_dict(cls, dct: dict) -> GrainBoundary: # type: ignore[override] - """Generate a GrainBoundary object from a dictionary created by as_dict(). + def from_dict(cls, dct: dict) -> Self: # type: ignore[override] + """Generate GrainBoundary from a dict created by as_dict(). Args: dct: dict @@ -289,7 +294,7 @@ def from_dict(cls, dct: dict) -> GrainBoundary: # type: ignore[override] sites = [PeriodicSite.from_dict(site_dict, lattice) for site_dict in dct["sites"]] struct = Structure.from_sites(sites) - return GrainBoundary( + return cls( lattice=lattice, species=struct.species_and_occu, coords=struct.frac_coords, @@ -307,11 +312,10 @@ def from_dict(cls, dct: dict) -> GrainBoundary: # type: ignore[override] class GrainBoundaryGenerator: """ - This class is to generate grain boundaries (GBs) from bulk - conventional cell (fcc, bcc can from the primitive cell), and works for Cubic, - Tetragonal, Orthorhombic, Rhombohedral, and Hexagonal systems. - It generate GBs from given parameters, which includes - GB plane, rotation axis, rotation angle. + Generate grain boundaries (GBs) from bulk conventional cell (FCC, BCC can + from the primitive cell), and works for Cubic, Tetragonal, Orthorhombic, + Rhombohedral, and Hexagonal systems. It generate GBs from given parameters, + which includes GB plane, rotation axis, rotation angle. This class works for any general GB, including twist, tilt and mixed GBs. The three parameters, rotation axis, GB plane and rotation angle, are @@ -321,10 +325,16 @@ class GrainBoundaryGenerator: rotation angles for a specific sigma value. The same sigma value (with rotation axis fixed) can correspond to multiple rotation angles. + Users can use structure matcher in pymatgen to get rid of the redundant structures. """ - def __init__(self, initial_structure: Structure, symprec: float = 0.1, angle_tolerance: float = 1) -> None: + def __init__( + self, + initial_structure: Structure, + symprec: float = 0.1, + angle_tolerance: float = 1.0, + ) -> None: """ Args: initial_structure (Structure): Initial input structure. It can @@ -344,8 +354,9 @@ def __init__(self, initial_structure: Structure, symprec: float = 0.1, angle_tol """ analyzer = SpacegroupAnalyzer(initial_structure, symprec, angle_tolerance) self.lat_type = analyzer.get_lattice_type()[0] + + # Use the conventional cell for tetragonal if self.lat_type == "t": - # need to use the conventional cell for tetragonal initial_structure = analyzer.get_conventional_standard_structure() a, b, c = initial_structure.lattice.abc # c axis of tetragonal structure not in the third direction @@ -356,6 +367,7 @@ def __init__(self, initial_structure: Structure, symprec: float = 0.1, angle_tol # b == c, rotate a to the third direction else: initial_structure.make_supercell([[0, 1, 0], [0, 0, 1], [1, 0, 0]]) + elif self.lat_type == "h": alpha, beta, gamma = initial_structure.lattice.angles # c axis is not in the third direction @@ -366,38 +378,40 @@ def __init__(self, initial_structure: Structure, symprec: float = 0.1, angle_tol # beta = 120 or 60, rotate c, a to a, b vectors elif abs(beta - 90) > angle_tolerance: initial_structure.make_supercell([[0, 0, 1], [1, 0, 0], [0, 1, 0]]) + + # Use primitive cell for rhombohedra elif self.lat_type == "r": - # need to use primitive cell for rhombohedra initial_structure = analyzer.get_primitive_standard_structure() + + # Use the conventional cell for orthorhombic elif self.lat_type == "o": - # need to use the conventional cell for orthorhombic initial_structure = analyzer.get_conventional_standard_structure() + self.initial_structure = initial_structure def gb_from_parameters( self, - rotation_axis, - rotation_angle, - expand_times=4, - vacuum_thickness=0.0, + rotation_axis: Tuple3Ints, + rotation_angle: float, + expand_times: int = 4, + vacuum_thickness: float = 0.0, ab_shift: tuple[float, float] = (0, 0), - normal=False, - ratio=None, - plane=None, - max_search=20, - tol_coi=1.0e-8, - rm_ratio=0.7, - quick_gen=False, + normal: bool = False, + ratio: list[int] | None = None, + plane: Tuple3Ints | None = None, + max_search: int = 20, + tol_coi: float = 1.0e-8, + rm_ratio: float = 0.7, + quick_gen: bool = False, ) -> GrainBoundary: """ Args: - rotation_axis (list): Rotation axis of GB in the form of a list of integer - e.g.: [1, 1, 0] + rotation_axis (tuple of 3): Rotation axis of GB e.g.: (1, 1, 0). rotation_angle (float, in unit of degree): rotation angle used to generate GB. Make sure the angle is accurate enough. You can use the enum* functions in this class to extract the accurate angle. e.g.: The rotation angle of sigma 3 twist GB with the rotation axis - [1, 1, 1] and GB plane (1, 1, 1) can be 60 degree. + (1, 1, 1) and GB plane (1, 1, 1) can be 60 degree. If you do not know the rotation angle, but know the sigma value, we have provide the function get_rotation_angle_from_sigma which is able to return all the rotation angles of sigma value you provided. @@ -407,12 +421,11 @@ def gb_from_parameters( vacuum_thickness (float, in angstrom): The thickness of vacuum that you want to insert between two grains of the GB. Default to 0. ab_shift (list of float, in unit of a, b vectors of Gb): in plane shift of two grains - normal (logic): + normal (bool): determine if need to require the c axis of top grain (first transformation matrix) perpendicular to the surface or not. default to false. - ratio (list of integers): - lattice axial ratio. + ratio (list[int]): lattice axial ratio. For cubic system, ratio is not needed. For tetragonal system, ratio = [mu, mv], list of two integers, that is, mu/mv = c2/a2. If it is irrational, set it to none. @@ -427,9 +440,8 @@ def gb_from_parameters( This code also supplies a class method to generate the ratio from the structure (get_ratio). User can also make their own approximation and input the ratio directly. - plane (list): Grain boundary plane in the form of a list of integers - e.g.: [1, 2, 3]. If none, we set it as twist GB. The plane will be perpendicular - to the rotation axis. + plane (tuple of 3): Grain boundary plane. If none, we set it as twist GB. + The plane will be perpendicular to the rotation axis. max_search (int): max search for the GB lattice vectors that give the smallest GB lattice. If normal is true, also max search the GB c vector that perpendicular to the plane. For complex GB, if you want to speed up, you can reduce this value. @@ -446,89 +458,107 @@ def gb_from_parameters( find the smallest cell. Returns: - Grain boundary structure (GB object). + GrainBoundary object """ - lat_type = self.lat_type - # if the initial structure is primitive cell in cubic system, + lat_type = self.lat_type.lower() + # If the initial structure is primitive cell in cubic system, # calculate the transformation matrix from its conventional cell - # to primitive cell, basically for bcc and fcc systems. + # to primitive cell, basically for BCC and FCC systems. trans_cry = np.eye(3) if lat_type == "c": analyzer = SpacegroupAnalyzer(self.initial_structure) convention_cell = analyzer.get_conventional_standard_structure() vol_ratio = self.initial_structure.volume / convention_cell.volume - # bcc primitive cell, belong to cubic system + # BCC primitive cell, belong to cubic system if abs(vol_ratio - 0.5) < 1.0e-3: trans_cry = np.array([[0.5, 0.5, -0.5], [-0.5, 0.5, 0.5], [0.5, -0.5, 0.5]]) logger.info("Make sure this is for cubic with bcc primitive cell") - # fcc primitive cell, belong to cubic system + # FCC primitive cell, belong to cubic system elif abs(vol_ratio - 0.25) < 1.0e-3: trans_cry = np.array([[0.5, 0.5, 0], [0, 0.5, 0.5], [0.5, 0, 0.5]]) logger.info("Make sure this is for cubic with fcc primitive cell") else: logger.info("Make sure this is for cubic with conventional cell") + elif lat_type == "t": logger.info("Make sure this is for tetragonal system") if ratio is None: logger.info("Make sure this is for irrational c2/a2") elif len(ratio) != 2: raise RuntimeError("Tetragonal system needs correct c2/a2 ratio") + elif lat_type == "o": logger.info("Make sure this is for orthorhombic system") if ratio is None: raise RuntimeError("CSL does not exist if all axial ratios are irrational for an orthorhombic system") if len(ratio) != 3: raise RuntimeError("Orthorhombic system needs correct c2:b2:a2 ratio") + elif lat_type == "h": logger.info("Make sure this is for hexagonal system") if ratio is None: logger.info("Make sure this is for irrational c2/a2") elif len(ratio) != 2: raise RuntimeError("Hexagonal system needs correct c2/a2 ratio") + elif lat_type == "r": logger.info("Make sure this is for rhombohedral system") if ratio is None: logger.info("Make sure this is for irrational (1+2*cos(alpha)/cos(alpha) ratio") elif len(ratio) != 2: raise RuntimeError("Rhombohedral system needs correct (1+2*cos(alpha)/cos(alpha) ratio") + else: raise RuntimeError( "Lattice type not implemented. This code works for cubic, " "tetragonal, orthorhombic, rhombohedral, hexagonal systems" ) - # transform four index notation to three index notation for hexagonal and rhombohedral + # Transform four index notation to three index notation for hexagonal and rhombohedral if len(rotation_axis) == 4: u1 = rotation_axis[0] v1 = rotation_axis[1] w1 = rotation_axis[3] - if lat_type.lower() == "h": + if lat_type == "h": u = 2 * u1 + v1 v = 2 * v1 + u1 w = w1 - rotation_axis = [u, v, w] - elif lat_type.lower() == "r": + _rotation_axis: Tuple3Ints | Tuple4Ints = (u, v, w) + elif lat_type == "r": u = 2 * u1 + v1 + w1 v = v1 + w1 - u1 w = w1 - 2 * v1 - u1 - rotation_axis = [u, v, w] - - # make sure gcd(rotation_axis)==1 - if reduce(gcd, rotation_axis) != 1: - rotation_axis = [int(round(x / reduce(gcd, rotation_axis))) for x in rotation_axis] - # transform four index notation to three index notation for plane - if plane is not None and len(plane) == 4: - u1, v1, w1 = plane[0], plane[1], plane[3] - plane = [u1, v1, w1] - # set the plane for grain boundary when plane is None. + _rotation_axis = (u, v, w) + else: + _rotation_axis = cast(Tuple4Ints, tuple(rotation_axis)) + + elif len(rotation_axis) == 3: + _rotation_axis = cast(Tuple3Ints, tuple(rotation_axis)) + + else: + raise ValueError("Invalid length of rotation axis.") + + # Make sure math.gcd(rotation_axis) == 1 + if reduce(math.gcd, _rotation_axis) != 1: + _rotation_axis = tuple(round(x / reduce(math.gcd, _rotation_axis)) for x in _rotation_axis) # type: ignore[assignment] + + # Transform four index notation to three index notation for plane + if plane is not None: + if len(plane) == 4: + u1, v1, w1 = plane[0], plane[1], plane[3] + plane = (u1, v1, w1) + elif len(plane) == 3: + plane = cast(Tuple3Ints, tuple(plane)) + + # Set the plane for grain boundary when plane is None if plane is None: - if lat_type.lower() == "c": - plane = rotation_axis + if lat_type == "c" and len(_rotation_axis) == 3: + _plane = _rotation_axis else: - if lat_type.lower() == "h": + if lat_type == "h": c2_a2_ratio = 1.0 if ratio is None else ratio[0] / ratio[1] metric = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, c2_a2_ratio]]) - elif lat_type.lower() == "r": + elif lat_type == "r": cos_alpha = 0.5 if ratio is None else 1.0 / (ratio[0] / ratio[1] - 2) metric = np.array( [ @@ -537,48 +567,51 @@ def gb_from_parameters( [cos_alpha, cos_alpha, 1], ] ) - elif lat_type.lower() == "t": + elif lat_type == "t": c2_a2_ratio = 1.0 if ratio is None else ratio[0] / ratio[1] metric = np.array([[1, 0, 0], [0, 1, 0], [0, 0, c2_a2_ratio]]) - elif lat_type.lower() == "o": + elif lat_type == "o" and ratio is not None: for idx in range(3): if ratio[idx] is None: ratio[idx] = 1 metric = np.array([[1, 0, 0], [0, ratio[1] / ratio[2], 0], [0, 0, ratio[0] / ratio[2]]]) else: - raise RuntimeError("Lattice type has not implemented.") + raise RuntimeError("Lattice type is not implemented.") - plane = np.matmul(rotation_axis, metric) - fractions = [Fraction(x).limit_denominator() for x in plane] + _plane = np.matmul(_rotation_axis, metric) + fractions = [Fraction(x).limit_denominator() for x in _plane] least_mul = reduce(lcm, [fraction.denominator for fraction in fractions]) - plane = [int(round(x * least_mul)) for x in plane] + _plane = cast(Tuple3Ints, tuple(round(x * least_mul) for x in _plane)) + + else: + _plane = plane - if reduce(gcd, plane) != 1: - index = reduce(gcd, plane) - plane = [int(round(x / index)) for x in plane] + if reduce(math.gcd, _plane) != 1: + index = reduce(math.gcd, _plane) + _plane = cast(Tuple3Ints, tuple(round(x / index) for x in _plane)) t1, t2 = self.get_trans_mat( - r_axis=rotation_axis, + r_axis=_rotation_axis, angle=rotation_angle, normal=normal, trans_cry=trans_cry, lat_type=lat_type, ratio=ratio, - surface=plane, + surface=_plane, max_search=max_search, quick_gen=quick_gen, ) - # find the join_plane - if lat_type.lower() != "c": - if lat_type.lower() == "h": + # Find the join_plane + if lat_type != "c": + if lat_type == "h": if ratio is None: mu, mv = [1, 1] else: mu, mv = ratio trans_cry1 = np.array([[1, 0, 0], [-0.5, np.sqrt(3.0) / 2.0, 0], [0, 0, np.sqrt(mu / mv)]]) - elif lat_type.lower() == "r": + elif lat_type == "r": if ratio is None: c2_a2_ratio = 1.0 else: @@ -593,14 +626,14 @@ def gb_from_parameters( ) else: - if lat_type.lower() == "t": + if lat_type == "t": if ratio is None: mu, mv = [1, 1] else: mu, mv = ratio lam = mv - elif lat_type.lower() == "o": + elif lat_type == "o" and ratio is not None: new_ratio = [1 if v is None else v for v in ratio] mu, lam, mv = new_ratio @@ -614,12 +647,12 @@ def gb_from_parameters( grain_matrix = np.dot(t2, trans_cry1) plane_init = np.cross(grain_matrix[0], grain_matrix[1]) - if lat_type.lower() != "c": + if lat_type != "c": plane_init = np.dot(plane_init, trans_cry1.T) join_plane = self.vec_to_surface(plane_init) parent_structure = self.initial_structure.copy() - # calculate the bond_length in bulk system. + # Calculate the bond_length in bulk system if len(parent_structure) == 1: temp_str = parent_structure.copy() temp_str.make_supercell([1, 1, 2]) @@ -628,19 +661,19 @@ def gb_from_parameters( distance = parent_structure.distance_matrix bond_length = np.min(distance[np.nonzero(distance)]) - # top grain + # Top grain top_grain = fix_pbc(parent_structure * t1) - # obtain the smallest oriented cell + # Obtain the smallest oriented cell if normal and not quick_gen: t_temp = self.get_trans_mat( - r_axis=rotation_axis, + r_axis=_rotation_axis, angle=rotation_angle, normal=False, trans_cry=trans_cry, lat_type=lat_type, ratio=ratio, - surface=plane, + surface=_plane, max_search=max_search, ) oriented_unit_cell = fix_pbc(parent_structure * t_temp[0]) @@ -654,10 +687,10 @@ def gb_from_parameters( oriented_unit_cell = top_grain.copy() unit_ab_adjust = 0.0 - # bottom grain, using top grain's lattice matrix + # Bottom grain, using top grain's lattice matrix bottom_grain = fix_pbc(parent_structure * t2, top_grain.lattice.matrix) - # label both grains with 'top','bottom','top_incident','bottom_incident' + # Label both grains with 'top', 'bottom', 'top_incident', 'bottom_incident' n_sites = len(top_grain) t_and_b = Structure( top_grain.lattice, @@ -694,38 +727,38 @@ def gb_from_parameters( site_properties={"grain_label": bottom_labels}, ) - # expand both grains + # Expand both grains top_grain.make_supercell([1, 1, expand_times]) bottom_grain.make_supercell([1, 1, expand_times]) top_grain = fix_pbc(top_grain) bottom_grain = fix_pbc(bottom_grain) - # determine the top-grain location. + # Determine the top-grain location edge_b = 1.0 - max(bottom_grain.frac_coords[:, 2]) edge_t = 1.0 - max(top_grain.frac_coords[:, 2]) c_adjust = (edge_t - edge_b) / 2.0 - # construct all species + # Construct all species all_species = [] all_species.extend([site.specie for site in bottom_grain]) all_species.extend([site.specie for site in top_grain]) half_lattice = top_grain.lattice - # calculate translation vector, perpendicular to the plane + # Calculate translation vector, perpendicular to the plane normal_v_plane = np.cross(half_lattice.matrix[0], half_lattice.matrix[1]) unit_normal_v = normal_v_plane / np.linalg.norm(normal_v_plane) translation_v = unit_normal_v * vacuum_thickness - # construct the final lattice + # Construct the final lattice whole_matrix_no_vac = np.array(half_lattice.matrix) whole_matrix_no_vac[2] = half_lattice.matrix[2] * 2 whole_matrix_with_vac = whole_matrix_no_vac.copy() whole_matrix_with_vac[2] = whole_matrix_no_vac[2] + translation_v * 2 whole_lat = Lattice(whole_matrix_with_vac) - # construct the coords, move top grain with translation_v + # Construct the coords, move top grain with translation_v all_coords = [] - grain_labels = bottom_grain.site_properties["grain_label"] + top_grain.site_properties["grain_label"] + grain_labels = bottom_grain.site_properties["grain_label"] + top_grain.site_properties["grain_label"] # type: ignore[operator] for site in bottom_grain: all_coords.append(site.coords) for site in top_grain: @@ -745,7 +778,7 @@ def gb_from_parameters( coords_are_cartesian=True, site_properties={"grain_label": grain_labels}, ) - # merge closer atoms. extract near GB atoms. + # Merge closer atoms. extract near GB atoms. cos_c_norm_plane = np.dot(unit_normal_v, whole_matrix_with_vac[2]) / whole_lat.c range_c_len = abs(bond_length / cos_c_norm_plane / whole_lat.c) sites_near_gb = [] @@ -762,18 +795,18 @@ def gb_from_parameters( if len(sites_near_gb) >= 1: s_near_gb = Structure.from_sites(sites_near_gb) s_near_gb.merge_sites(tol=bond_length * rm_ratio, mode="delete") - all_sites = sites_away_gb + s_near_gb.sites # type: ignore + all_sites = sites_away_gb + s_near_gb.sites # type: ignore[operator] gb_with_vac = Structure.from_sites(all_sites) - # move coordinates into the periodic cell. + # Move coordinates into the periodic cell gb_with_vac = fix_pbc(gb_with_vac, whole_lat.matrix) return GrainBoundary( whole_lat, gb_with_vac.species, gb_with_vac.cart_coords, # type: ignore[arg-type] - rotation_axis, + _rotation_axis, rotation_angle, - plane, + _plane, join_plane, self.initial_structure, vacuum_thickness, @@ -783,9 +816,12 @@ def gb_from_parameters( coords_are_cartesian=True, ) - def get_ratio(self, max_denominator=5, index_none=None): - """ - find the axial ratio needed for GB generator input. + def get_ratio( + self, + max_denominator: int = 5, + index_none: int | None = None, + ) -> list[int] | None: + """Find the axial ratio needed for GB generator input. Args: max_denominator (int): the maximum denominator for @@ -798,22 +834,25 @@ def get_ratio(self, max_denominator=5, index_none=None): """ structure = self.initial_structure lat_type = self.lat_type - if lat_type in ("t", "h"): - # For tetragonal and hexagonal system, ratio = c2 / a2. + + # For tetragonal and hexagonal systems, ratio = c2 / a2 + if lat_type in {"t", "h"}: a, _, c = structure.lattice.lengths if c > a: frac = Fraction(c**2 / a**2).limit_denominator(max_denominator) - ratio = [frac.numerator, frac.denominator] + ratio: list[int | None] = [frac.numerator, frac.denominator] else: frac = Fraction(a**2 / c**2).limit_denominator(max_denominator) ratio = [frac.denominator, frac.numerator] + + # For rhombohedral system, ratio = (1 + 2 * cos(alpha)) / cos(alpha) elif lat_type == "r": - # For rhombohedral system, ratio = (1 + 2 * cos(alpha)) / cos(alpha). - cos_alpha = cos(structure.lattice.alpha / 180 * np.pi) + cos_alpha = math.cos(structure.lattice.alpha / 180 * np.pi) frac = Fraction((1 + 2 * cos_alpha) / cos_alpha).limit_denominator(max_denominator) ratio = [frac.numerator, frac.denominator] + + # For orthorhombic system, ratio = c2:b2:a2. If irrational for one axis, set it to None elif lat_type == "o": - # For orthorhombic system, ratio = c2:b2:a2.If irrational for one axis, set it to None. ratio = [None] * 3 lat = (structure.lattice.c, structure.lattice.b, structure.lattice.a) index = [0, 1, 2] @@ -824,8 +863,9 @@ def get_ratio(self, max_denominator=5, index_none=None): frac2 = Fraction(lat[index[1]] ** 2 / lat[min_index] ** 2).limit_denominator(max_denominator) com_lcm = lcm(frac1.denominator, frac2.denominator) ratio[min_index] = com_lcm - ratio[index[0]] = frac1.numerator * int(round(com_lcm / frac1.denominator)) - ratio[index[1]] = frac2.numerator * int(round(com_lcm / frac2.denominator)) + ratio[index[0]] = frac1.numerator * round(com_lcm / frac1.denominator) + ratio[index[1]] = frac2.numerator * round(com_lcm / frac2.denominator) + else: index.pop(index_none) if lat[index[0]] > lat[index[1]]: @@ -836,24 +876,26 @@ def get_ratio(self, max_denominator=5, index_none=None): frac = Fraction(lat[index[1]] ** 2 / lat[index[0]] ** 2).limit_denominator(max_denominator) ratio[index[1]] = frac.numerator ratio[index[0]] = frac.denominator + + # Cubic system does not need axial ratio elif lat_type == "c": - # Cubic system does not need axial ratio. return None + else: raise RuntimeError("Lattice type not implemented.") - return ratio + return cast(list[int], ratio) @staticmethod def get_trans_mat( - r_axis, - angle, - normal=False, - trans_cry=None, - lat_type="c", - ratio=None, - surface=None, - max_search=20, - quick_gen=False, + r_axis: Tuple3Ints | Tuple4Ints, + angle: float, + normal: bool = False, + trans_cry: NDArray | None = None, + lat_type: str = "c", + ratio: list[int] | None = None, + surface: Tuple3Ints | Tuple4Ints | None = None, + max_search: int = 20, + quick_gen: bool = False, ): """Find the two transformation matrix for each grain from given rotation axis, GB plane, rotation angle and corresponding ratio (see explanation for ratio @@ -863,10 +905,10 @@ def get_trans_mat( The algorithm for this code is from reference, Acta Cryst, A32,783(1976). Args: - r_axis (list of 3 integers, e.g. u, v, w or 4 integers, e.g. u, v, t, w for hex/rho system only): the - rotation axis of the grain boundary. + r_axis ((u, v, w) or (u, v, t, w) for hex/rho systems): + the rotation axis of the grain boundary. angle (float, in unit of degree): the rotation angle of the grain boundary - normal (logic): determine if need to require the c axis of one grain associated with + normal (bool): determine if need to require the c axis of one grain associated with the first transformation matrix perpendicular to the surface or not. default to false. trans_cry (np.array): shape 3x3. If the structure given are primitive cell in cubic system, e.g. @@ -878,7 +920,7 @@ def get_trans_mat( 'o' or 'O': orthorhombic system 'h' or 'H': hexagonal system 'r' or 'R': rhombohedral system - ratio (list of integers): lattice axial ratio. + ratio (list[int]): lattice axial ratio. For cubic system, ratio is not needed. For tetragonal system, ratio = [mu, mv], list of two integers, that is, mu/mv = c2/a2. If it is irrational, set it to none. @@ -889,9 +931,9 @@ def get_trans_mat( If irrational, set it to None. For hexagonal system, ratio = [mu, mv], list of two integers, that is, mu/mv = c2/a2. If it is irrational, set it to none. - surface (list of 3 integers, e.g. h, k, l or 4 integers, e.g. h, k, i, l for hex/rho system only): The - miller index of grain boundary plane, with the format of [h,k,l] if surface is not given, the default - is perpendicular to r_axis, which is a twist grain boundary. + surface ((h, k, l) or (h, k, i, l) for hex/rho systems): The + miller index of grain boundary plane, with the format of (h, k, l) if surface + is not given, the default is perpendicular to r_axis, which is a twist grain boundary. max_search (int): max search for the GB lattice vectors that give the smallest GB lattice. If normal is true, also max search the GB c vector that perpendicular to the plane. @@ -902,42 +944,47 @@ def get_trans_mat( t1 (3 by 3 integer array): The transformation array for one grain. t2 (3 by 3 integer array): The transformation array for the other grain """ - if trans_cry is None: - trans_cry = np.eye(3) - # transform four index notation to three index notation + trans_cry = np.eye(3) if trans_cry is None else trans_cry + lat_type = lat_type.lower() + + # Transform four index notation to three index notation if len(r_axis) == 4: u1 = r_axis[0] v1 = r_axis[1] w1 = r_axis[3] - if lat_type.lower() == "h": + if lat_type == "h": u = 2 * u1 + v1 v = 2 * v1 + u1 w = w1 - r_axis = [u, v, w] - elif lat_type.lower() == "r": + r_axis = (u, v, w) + elif lat_type == "r": u = 2 * u1 + v1 + w1 v = v1 + w1 - u1 w = w1 - 2 * v1 - u1 - r_axis = [u, v, w] + r_axis = (u, v, w) - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + # Make sure gcd(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast( + Union[Tuple3Ints, Tuple4Ints], + tuple(round(x / reduce(math.gcd, r_axis)) for x in r_axis), + ) if surface is not None and len(surface) == 4: u1 = surface[0] v1 = surface[1] w1 = surface[3] - surface = [u1, v1, w1] - # set the surface for grain boundary. + surface = (u1, v1, w1) + + # Set the surface for grain boundary. if surface is None: - if lat_type.lower() == "c": + if lat_type == "c": surface = r_axis else: - if lat_type.lower() == "h": + if lat_type == "h": c2_a2_ratio = 1.0 if ratio is None else ratio[0] / ratio[1] metric = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, c2_a2_ratio]]) - elif lat_type.lower() == "r": + elif lat_type == "r": cos_alpha = 0.5 if ratio is None else 1.0 / (ratio[0] / ratio[1] - 2) metric = np.array( [ @@ -946,12 +993,13 @@ def get_trans_mat( [cos_alpha, cos_alpha, 1], ] ) - elif lat_type.lower() == "t": + elif lat_type == "t": c2_a2_ratio = 1.0 if ratio is None else ratio[0] / ratio[1] metric = np.array([[1, 0, 0], [0, 1, 0], [0, 0, c2_a2_ratio]]) - elif lat_type.lower() == "o": + elif lat_type == "o": + assert ratio is not None, "Invalid ratio for orthorhombic system" for idx in range(3): - if ratio[idx] is None: + if ratio is not None and ratio[idx] is None: ratio[idx] = 1 metric = np.array( [ @@ -966,28 +1014,28 @@ def get_trans_mat( surface = np.matmul(r_axis, metric) fractions = [Fraction(x).limit_denominator() for x in surface] least_mul = reduce(lcm, [fraction.denominator for fraction in fractions]) - surface = [int(round(x * least_mul)) for x in surface] + surface = cast(Union[Tuple3Ints, Tuple4Ints], tuple(round(x * least_mul) for x in surface)) - if reduce(gcd, surface) != 1: - index = reduce(gcd, surface) - surface = [int(round(x / index)) for x in surface] + if reduce(math.gcd, surface) != 1: + index = reduce(math.gcd, surface) + surface = cast(Union[Tuple3Ints, Tuple4Ints], tuple(round(x / index) for x in surface)) lam = None - if lat_type.lower() == "h": - # set the value for u,v,w,mu,mv,m,n,d,x - # check the reference for the meaning of these parameters - u, v, w = r_axis - # make sure mu, mv are coprime integers. + if lat_type == "h": + # Set the values for u, v, w, mu, mv, m, n, d, x + # Check the reference for the meaning of these parameters + u, v, w = cast(Tuple3Ints, r_axis) + # Make sure mu, mv are coprime integers if ratio is None: - mu, mv = [1, 1] + mu, mv = (1, 1) if w != 0 and (u != 0 or (v != 0)): raise RuntimeError("For irrational c2/a2, CSL only exist for [0,0,1] or [u,v,0] and m = 0") else: mu, mv = ratio - if gcd(mu, mv) != 1: - temp = gcd(mu, mv) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) + if math.gcd(mu, mv) != 1: + temp = math.gcd(mu, mv) + mu = round(mu / temp) + mv = round(mv / temp) d = (u**2 + v**2 - u * v) * mv + w**2 * mu if abs(angle - 180.0) < 1.0e0: m = 0 @@ -999,7 +1047,7 @@ def get_trans_mat( m = fraction.denominator n = fraction.numerator - # construct the rotation matrix, check reference for details + # Construct the rotation matrix, check reference for details r_list = [ (u**2 * mv - v**2 * mv - w**2 * mu) * n**2 + 2 * w * mu * m * n + 3 * mu * m**2, (2 * v - u) * u * mv * n**2 - 4 * w * mu * m * n, @@ -1026,27 +1074,27 @@ def get_trans_mat( m = -1 * m F = 3 * mu * m**2 + d * n**2 all_list = r_list + r_list_inv + [F] - com_fac = reduce(gcd, all_list) + com_fac = reduce(math.gcd, all_list) sigma = F / com_fac r_matrix = (np.array(r_list) / com_fac / sigma).reshape(3, 3) - elif lat_type.lower() == "r": - # set the value for u,v,w,mu,mv,m,n,d - # check the reference for the meaning of these parameters - u, v, w = r_axis - # make sure mu, mv are coprime integers. + elif lat_type == "r": + # Set the values for u, v, w, mu, mv ,m, n, d + # Check the reference for the meaning of these parameters + u, v, w = cast(Tuple3Ints, r_axis) + # make sure mu, mv are coprime integers if ratio is None: - mu, mv = [1, 1] + mu, mv = (1, 1) if u + v + w != 0 and (u != v or u != w): raise RuntimeError( "For irrational ratio_alpha, CSL only exist for [1,1,1] or [u, v, -(u+v)] and m =0" ) else: mu, mv = ratio - if gcd(mu, mv) != 1: - temp = gcd(mu, mv) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) + if math.gcd(mu, mv) != 1: + temp = math.gcd(mu, mv) + mu = round(mu / temp) + mv = round(mv / temp) d = (u**2 + v**2 + w**2) * (mu - 2 * mv) + 2 * mv * (v * w + w * u + u * v) if abs(angle - 180.0) < 1.0e0: m = 0 @@ -1056,7 +1104,7 @@ def get_trans_mat( m = fraction.denominator n = fraction.numerator - # construct the rotation matrix, check reference for details + # Construct the rotation matrix, check reference for details r_list = [ (mu - 2 * mv) * (u**2 - v**2 - w**2) * n**2 + 2 * mv * (v - w) * m * n @@ -1101,27 +1149,27 @@ def get_trans_mat( m = -1 * m F = mu * m**2 + d * n**2 all_list = r_list_inv + r_list + [F] - com_fac = reduce(gcd, all_list) + com_fac = reduce(math.gcd, all_list) sigma = F / com_fac r_matrix = (np.array(r_list) / com_fac / sigma).reshape(3, 3) else: - u, v, w = r_axis - mu = mv = None - if lat_type.lower() == "c": + u, v, w = cast(Tuple3Ints, r_axis) + mu = mv = None # type: ignore[assignment] + if lat_type == "c": mu = 1 lam = 1 mv = 1 - elif lat_type.lower() == "t": + elif lat_type == "t": if ratio is None: - mu, mv = [1, 1] + mu, mv = (1, 1) if w != 0 and (u != 0 or (v != 0)): raise RuntimeError("For irrational c2/a2, CSL only exist for [0,0,1] or [u,v,0] and m = 0") else: mu, mv = ratio lam = mv - elif lat_type.lower() == "o": - if None in ratio: + elif lat_type == "o": + if ratio is not None and None in ratio: mu, lam, mv = ratio non_none = [i for i in ratio if i is not None] if len(non_none) < 2: @@ -1146,7 +1194,7 @@ def get_trans_mat( if u != 0 and (w != 0 or (v != 0)): raise RuntimeError("For irrational a2, CSL only exist for [1,0,0] or [0,v,w] and m = 0") else: - mu, lam, mv = ratio + mu, lam, mv = cast(list[int], ratio) if u == 0 and v == 0: mu = 1 if u == 0 and w == 0: @@ -1154,12 +1202,15 @@ def get_trans_mat( if v == 0 and w == 0: mv = 1 - # make sure mu, lambda, mv are coprime integers. - if reduce(gcd, [mu, lam, mv]) != 1: - temp = reduce(gcd, [mu, lam, mv]) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) - lam = int(round(lam / temp)) + # Make sure mu, lambda, mv are coprime integers + assert mu is not None + assert lam is not None + assert mv is not None + if reduce(math.gcd, [mu, lam, mv]) != 1: + temp = cast(int, reduce(math.gcd, [mu, lam, mv])) + mu = round(mu / temp) + mv = round(mv / temp) + lam = round(lam / temp) d = (mv * u**2 + lam * v**2) * mv + w**2 * mu * mv if abs(angle - 180.0) < 1.0e0: m = 0 @@ -1194,25 +1245,29 @@ def get_trans_mat( m = -1 * m F = mu * lam * m**2 + d * n**2 all_list = r_list + r_list_inv + [F] - com_fac = reduce(gcd, all_list) + com_fac = reduce(math.gcd, all_list) sigma = F / com_fac r_matrix = (np.array(r_list) / com_fac / sigma).reshape(3, 3) if sigma > 1000: raise RuntimeError("Sigma >1000 too large. Are you sure what you are doing, Please check the GB if exist") - # transform surface, r_axis, r_matrix in terms of primitive lattice + # Transform surface, r_axis, r_matrix in terms of primitive lattice surface = np.matmul(surface, np.transpose(trans_cry)) + assert surface is not None fractions = [Fraction(x).limit_denominator() for x in surface] least_mul = reduce(lcm, [fraction.denominator for fraction in fractions]) - surface = [int(round(x * least_mul)) for x in surface] - if reduce(gcd, surface) != 1: - index = reduce(gcd, surface) - surface = [int(round(x / index)) for x in surface] + surface = cast(Tuple3Ints, tuple(round(x * least_mul) for x in surface)) + if reduce(math.gcd, surface) != 1: + index = reduce(math.gcd, surface) + surface = cast(Tuple3Ints, tuple(round(x / index) for x in surface)) r_axis = np.rint(np.matmul(r_axis, np.linalg.inv(trans_cry))).astype(int) - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + if reduce(math.gcd, r_axis) != 1: + r_axis = cast( + Tuple3Ints, + tuple(round(x / reduce(math.gcd, r_axis)) for x in r_axis), + ) r_matrix = np.dot(np.dot(np.linalg.inv(trans_cry.T), r_matrix), trans_cry.T) - # set one vector of the basis to the rotation axis direction, and + # Set one vector of the basis to the rotation axis direction, and # obtain the corresponding transform matrix eye = np.eye(3, dtype=int) hh = kk = ll = None @@ -1225,7 +1280,7 @@ def get_trans_mat( trans = eye.T new_rot = np.array(r_matrix) - # with the rotation matrix to construct the CSL lattice, check reference for details + # With the rotation matrix to construct the CSL lattice, check reference for details fractions = [Fraction(x).limit_denominator() for x in new_rot[:, kk]] least_mul = reduce(lcm, [fraction.denominator for fraction in fractions]) scale = np.zeros((3, 3)) @@ -1242,22 +1297,22 @@ def get_trans_mat( if n_final is None: raise RuntimeError("Something is wrong. Check if this GB exists or not") scale[kk, ll] = n_final - # each row of mat_csl is the CSL lattice vector + # Each row of mat_csl is the CSL lattice vector csl_init = np.rint(np.dot(np.dot(r_matrix, trans), scale)).astype(int).T if abs(r_axis[hh]) > 1: csl_init = GrainBoundaryGenerator.reduce_mat(np.array(csl_init), r_axis[hh], r_matrix) csl = np.rint(Lattice(csl_init).get_niggli_reduced_lattice().matrix).astype(int) - # find the best slab supercell in terms of the conventional cell from the csl lattice, + # Find the best slab supercell in terms of the conventional cell from the csl lattice, # which is the transformation matrix - # now trans_cry is the transformation matrix from crystal to Cartesian coordinates. + # Now trans_cry is the transformation matrix from crystal to Cartesian coordinates. # for cubic, do not need to change. - if lat_type.lower() != "c": - if lat_type.lower() == "h": - trans_cry = np.array([[1, 0, 0], [-0.5, np.sqrt(3.0) / 2.0, 0], [0, 0, np.sqrt(mu / mv)]]) - elif lat_type.lower() == "r": - c2_a2_ratio = 1.0 if ratio is None else 3.0 / (2 - 6 * mv / mu) + if lat_type != "c": + if lat_type == "h": + trans_cry = np.array([[1, 0, 0], [-0.5, np.sqrt(3.0) / 2.0, 0], [0, 0, np.sqrt(mu / mv)]]) # type: ignore[operator] + elif lat_type == "r": + c2_a2_ratio = 1.0 if ratio is None else 3.0 / (2 - 6 * mv / mu) # type: ignore[operator] trans_cry = np.array( [ [0.5, np.sqrt(3.0) / 6.0, 1.0 / 3 * np.sqrt(c2_a2_ratio)], @@ -1266,7 +1321,7 @@ def get_trans_mat( ] ) else: - trans_cry = np.array([[1, 0, 0], [0, np.sqrt(lam / mv), 0], [0, 0, np.sqrt(mu / mv)]]) + trans_cry = np.array([[1, 0, 0], [0, np.sqrt(lam / mv), 0], [0, 0, np.sqrt(mu / mv)]]) # type: ignore[operator] t1_final = GrainBoundaryGenerator.slab_from_csl( csl, surface, normal, trans_cry, max_search=max_search, quick_gen=quick_gen ) @@ -1274,15 +1329,17 @@ def get_trans_mat( return t1_final, t2_final @staticmethod - def enum_sigma_cubic(cutoff, r_axis): + def enum_sigma_cubic( + cutoff: int, + r_axis: Tuple3Ints, + ) -> dict[int, list[float]]: """Find all possible sigma values and corresponding rotation angles within a sigma value cutoff with known rotation axis in cubic system. The algorithm for this code is from reference, Acta Cryst, A40,108(1984). Args: cutoff (int): the cutoff of sigma values. - r_axis (list of 3 integers, e.g. u, v, w): - the rotation axis of the grain boundary, with the format of [u,v,w]. + r_axis ((u, v, w)): the rotation axis of the grain boundary. Returns: dict: sigmas dictionary with keys as the possible integer sigma values @@ -1295,14 +1352,13 @@ def enum_sigma_cubic(cutoff, r_axis): you need to analyze the symmetry of the structure. Different angles may result in equivalent microstructures. """ - sigmas = {} - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + # Make sure math.gcd(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast(Tuple3Ints, tuple(round(x / reduce(math.gcd, r_axis)) for x in r_axis)) - # count the number of odds in r_axis + # Count the number of odds in r_axis odd_r = len(list(filter(lambda x: x % 2 == 1, r_axis))) - # Compute the max n we need to enumerate. + # Compute the max n we need to enumerate if odd_r == 3: a_max = 4 elif odd_r == 0: @@ -1310,14 +1366,15 @@ def enum_sigma_cubic(cutoff, r_axis): else: a_max = 2 n_max = int(np.sqrt(cutoff * a_max / sum(np.array(r_axis) ** 2))) - # enumerate all possible n, m to give possible sigmas within the cutoff. + # Enumerate all possible n, m to give possible sigmas within the cutoff + sigmas: dict[int, list[float]] = {} for n_loop in range(1, n_max + 1): n = n_loop m_max = int(np.sqrt(cutoff * a_max - n**2 * sum(np.array(r_axis) ** 2))) for m in range(m_max + 1): - if gcd(m, n) == 1 or m == 0: + if math.gcd(m, n) == 1 or m == 0: n = 1 if m == 0 else n_loop - # construct the quadruple [m, U,V,W], count the number of odds in + # Construct the quadruple [m, U,V,W], count the number of odds in # quadruple to determine the parameter a, refer to the reference quadruple = [m] + [x * n for x in r_axis] odd_qua = len(list(filter(lambda x: x % 2 == 1, quadruple))) @@ -1327,7 +1384,7 @@ def enum_sigma_cubic(cutoff, r_axis): a = 2 else: a = 1 - sigma = int(round((m**2 + n**2 * sum(np.array(r_axis) ** 2)) / a)) + sigma = round((m**2 + n**2 * sum(np.array(r_axis) ** 2)) / a) if 1 < sigma <= cutoff: if sigma not in list(sigmas): if m == 0: @@ -1345,34 +1402,40 @@ def enum_sigma_cubic(cutoff, r_axis): return sigmas @staticmethod - def enum_sigma_hex(cutoff, r_axis, c2_a2_ratio): + def enum_sigma_hex( + cutoff: int, + r_axis: Tuple3Ints | Tuple4Ints, + c2_a2_ratio: tuple[int, int], + ) -> dict[int, list[float]]: """Find all possible sigma values and corresponding rotation angles within a sigma value cutoff with known rotation axis in hexagonal system. The algorithm for this code is from reference, Acta Cryst, A38,550(1982). Args: cutoff (int): the cutoff of sigma values. - r_axis (list of 3 integers, e.g. u, v, w or 4 integers, e.g. u, v, t, w): the rotation axis of the grain - boundary. - c2_a2_ratio (list of two integers, e.g. mu, mv): mu/mv is the square of the hexagonal axial ratio, + r_axis ((u, v, w) or (u, v, t, w)): the rotation axis of the grain boundary. + c2_a2_ratio ((mu, mv)): mu/mv is the square of the hexagonal axial ratio, which is rational number. If irrational, set c2_a2_ratio = None Returns: dict: sigmas dictionary with keys as the possible integer sigma values and values as list of the possible rotation angles to the corresponding sigma values. e.g. the format as - {sigma1: [angle11,angle12,...], sigma2: [angle21, angle22,...],...} + {sigma1: [angle11, angle12, ...], sigma2: [angle21, angle22, ...], ...} Note: the angles are the rotation angles of one grain respect to the other grain. When generate the microstructures of the grain boundary using these angles, you need to analyze the symmetry of the structure. Different angles may result in equivalent microstructures. """ - sigmas = {} - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] - # transform four index notation to three index notation + # Make sure math.gcd(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast( + Union[Tuple3Ints, Tuple4Ints], + tuple([round(x / reduce(math.gcd, r_axis)) for x in r_axis]), + ) + + # Transform four index notation to three index notation if len(r_axis) == 4: u1 = r_axis[0] v1 = r_axis[1] @@ -1383,33 +1446,34 @@ def enum_sigma_hex(cutoff, r_axis, c2_a2_ratio): else: u, v, w = r_axis - # make sure mu, mv are coprime integers. + # Make sure mu, mv are coprime integers if c2_a2_ratio is None: mu, mv = [1, 1] if w != 0 and (u != 0 or (v != 0)): raise RuntimeError("For irrational c2/a2, CSL only exist for [0,0,1] or [u,v,0] and m = 0") else: mu, mv = c2_a2_ratio - if gcd(mu, mv) != 1: - temp = gcd(mu, mv) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) + if math.gcd(mu, mv) != 1: + temp = math.gcd(mu, mv) + mu = round(mu / temp) + mv = round(mv / temp) - # refer to the meaning of d in reference + # Refer to the meaning of d in reference d = (u**2 + v**2 - u * v) * mv + w**2 * mu - # Compute the max n we need to enumerate. + # Compute the max n we need to enumerate n_max = int(np.sqrt((cutoff * 12 * mu * mv) / abs(d))) - # Enumerate all possible n, m to give possible sigmas within the cutoff. + # Enumerate all possible n, m to give possible sigmas within the cutoff + sigmas: dict[int, list[float]] = {} for n in range(1, n_max + 1): if (c2_a2_ratio is None) and w == 0: m_max = 0 else: m_max = int(np.sqrt((cutoff * 12 * mu * mv - n**2 * d) / (3 * mu))) for m in range(m_max + 1): - if gcd(m, n) == 1 or m == 0: - # construct the rotation matrix, refer to the reference + if math.gcd(m, n) == 1 or m == 0: + # Construct the rotation matrix, refer to the reference R_list = [ (u**2 * mv - v**2 * mv - w**2 * mu) * n**2 + 2 * w * mu * m * n + 3 * mu * m**2, (2 * v - u) * u * mv * n**2 - 4 * w * mu * m * n, @@ -1422,7 +1486,7 @@ def enum_sigma_hex(cutoff, r_axis, c2_a2_ratio): (w**2 * mu - u**2 * mv - v**2 * mv + u * v * mv) * n**2 + 3 * mu * m**2, ] m = -1 * m - # inverse of the rotation matrix + # Inverse of the rotation matrix R_list_inv = [ (u**2 * mv - v**2 * mv - w**2 * mu) * n**2 + 2 * w * mu * m * n + 3 * mu * m**2, (2 * v - u) * u * mv * n**2 - 4 * w * mu * m * n, @@ -1439,8 +1503,8 @@ def enum_sigma_hex(cutoff, r_axis, c2_a2_ratio): all_list = R_list_inv + R_list + [F] # Compute the max common factors for the elements of the rotation matrix # and its inverse. - com_fac = reduce(gcd, all_list) - sigma = int(round((3 * mu * m**2 + d * n**2) / com_fac)) + com_fac = reduce(math.gcd, all_list) + sigma = round((3 * mu * m**2 + d * n**2) / com_fac) if 1 < sigma <= cutoff: if sigma not in list(sigmas): angle = 180.0 if m == 0 else 2 * np.arctan(n / m * np.sqrt(d / 3.0 / mu)) / np.pi * 180 @@ -1454,18 +1518,19 @@ def enum_sigma_hex(cutoff, r_axis, c2_a2_ratio): return sigmas @staticmethod - def enum_sigma_rho(cutoff, r_axis, ratio_alpha): + def enum_sigma_rho( + cutoff: int, + r_axis: Tuple3Ints | Tuple4Ints, + ratio_alpha: tuple[int, int], + ) -> dict[int, list[float]]: """Find all possible sigma values and corresponding rotation angles within a sigma value cutoff with known rotation axis in rhombohedral system. The algorithm for this code is from reference, Acta Cryst, A45,505(1989). Args: cutoff (int): the cutoff of sigma values. - r_axis (list[int]): of 3 integers, e.g. u, v, w - or 4 integers, e.g. u, v, t, w): - the rotation axis of the grain boundary, with the format of [u,v,w] - or Weber indices [u, v, t, w]. - ratio_alpha (list of two integers, e.g. mu, mv): + r_axis ((u, v, w) or (u, v, t, w)): the rotation axis of the grain boundary. + ratio_alpha (tuple of two integers, e.g. mu, mv): mu/mv is the ratio of (1+2*cos(alpha))/cos(alpha) with rational number. If irrational, set ratio_alpha = None. @@ -1479,8 +1544,7 @@ def enum_sigma_rho(cutoff, r_axis, ratio_alpha): angles, you need to analyze the symmetry of the structure. Different angles may result in equivalent microstructures. """ - sigmas = {} - # transform four index notation to three index notation + # Transform four index notation to three index notation if len(r_axis) == 4: u1 = r_axis[0] v1 = r_axis[1] @@ -1488,37 +1552,40 @@ def enum_sigma_rho(cutoff, r_axis, ratio_alpha): u = 2 * u1 + v1 + w1 v = v1 + w1 - u1 w = w1 - 2 * v1 - u1 - r_axis = [u, v, w] - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + r_axis = (u, v, w) + + # Make sure math.(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast(Tuple3Ints, tuple([round(x / reduce(math.gcd, r_axis)) for x in r_axis])) u, v, w = r_axis - # make sure mu, mv are coprime integers. + + # Make sure mu, mv are coprime integers if ratio_alpha is None: mu, mv = [1, 1] if u + v + w != 0 and (u != v or u != w): raise RuntimeError("For irrational ratio_alpha, CSL only exist for [1,1,1] or [u, v, -(u+v)] and m =0") else: mu, mv = ratio_alpha - if gcd(mu, mv) != 1: - temp = gcd(mu, mv) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) + if math.gcd(mu, mv) != 1: + temp = math.gcd(mu, mv) + mu = round(mu / temp) + mv = round(mv / temp) - # refer to the meaning of d in reference + # Refer to the meaning of d in reference d = (u**2 + v**2 + w**2) * (mu - 2 * mv) + 2 * mv * (v * w + w * u + u * v) - # Compute the max n we need to enumerate. + # Compute the max n we need to enumerate n_max = int(np.sqrt((cutoff * abs(4 * mu * (mu - 3 * mv))) / abs(d))) - # Enumerate all possible n, m to give possible sigmas within the cutoff. + # Enumerate all possible n, m to give possible sigmas within the cutoff + sigmas: dict[int, list[float]] = {} for n in range(1, n_max + 1): if ratio_alpha is None and u + v + w == 0: m_max = 0 else: m_max = int(np.sqrt((cutoff * abs(4 * mu * (mu - 3 * mv)) - n**2 * d) / (mu))) for m in range(m_max + 1): - if gcd(m, n) == 1 or m == 0: - # construct the rotation matrix, refer to the reference + if math.gcd(m, n) == 1 or m == 0: + # Construct the rotation matrix, refer to the reference R_list = [ (mu - 2 * mv) * (u**2 - v**2 - w**2) * n**2 + 2 * mv * (v - w) * m * n @@ -1540,7 +1607,7 @@ def enum_sigma_rho(cutoff, r_axis, ratio_alpha): + mu * m**2, ] m = -1 * m - # inverse of the rotation matrix + # Inverse of the rotation matrix R_list_inv = [ (mu - 2 * mv) * (u**2 - v**2 - w**2) * n**2 + 2 * mv * (v - w) * m * n @@ -1565,8 +1632,8 @@ def enum_sigma_rho(cutoff, r_axis, ratio_alpha): F = mu * m**2 + d * n**2 all_list = R_list_inv + R_list + [F] # Compute the max common factors for the elements of the rotation matrix and its inverse. - com_fac = reduce(gcd, all_list) - sigma = int(round(abs(F / com_fac))) + com_fac = reduce(math.gcd, all_list) + sigma = round(abs(F / com_fac)) if 1 < sigma <= cutoff: if sigma not in list(sigmas): angle = 180.0 if m == 0 else 2 * np.arctan(n / m * np.sqrt(d / mu)) / np.pi * 180 @@ -1580,61 +1647,63 @@ def enum_sigma_rho(cutoff, r_axis, ratio_alpha): return sigmas @staticmethod - def enum_sigma_tet(cutoff, r_axis, c2_a2_ratio): + def enum_sigma_tet( + cutoff: int, + r_axis: Tuple3Ints, + c2_a2_ratio: tuple[int, int], + ) -> dict[int, list[float]]: """Find all possible sigma values and corresponding rotation angles within a sigma value cutoff with known rotation axis in tetragonal system. The algorithm for this code is from reference, Acta Cryst, B46,117(1990). Args: cutoff (int): the cutoff of sigma values. - r_axis (list of 3 integers, e.g. u, v, w): - the rotation axis of the grain boundary, with the format of [u,v,w]. - c2_a2_ratio (list of two integers, e.g. mu, mv): - mu/mv is the square of the tetragonal axial ratio with rational number. - if irrational, set c2_a2_ratio = None + r_axis ((u, v, w)): the rotation axis of the grain boundary. + c2_a2_ratio ((mu, mv)): mu/mv is the square of the tetragonal axial + ratio with rational number. If irrational, set c2_a2_ratio = None. Returns: dict: sigmas dictionary with keys as the possible integer sigma values and values as list of the possible rotation angles to the corresponding sigma values. e.g. the format as - {sigma1: [angle11,angle12,...], sigma2: [angle21, angle22,...],...} + {sigma1: [angle11, angle12, ...], sigma2: [angle21, angle22, ...], ...} Note: the angles are the rotation angle of one grain respect to the other grain. When generate the microstructure of the grain boundary using these angles, you need to analyze the symmetry of the structure. Different angles may result in equivalent microstructures. """ - sigmas = {} - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + # Make sure math.gcd(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast(Tuple3Ints, tuple(round(x / reduce(math.gcd, r_axis)) for x in r_axis)) u, v, w = r_axis - # make sure mu, mv are coprime integers. + # Make sure mu, mv are coprime integers if c2_a2_ratio is None: mu, mv = [1, 1] if w != 0 and (u != 0 or (v != 0)): raise RuntimeError("For irrational c2/a2, CSL only exist for [0,0,1] or [u,v,0] and m = 0") else: mu, mv = c2_a2_ratio - if gcd(mu, mv) != 1: - temp = gcd(mu, mv) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) + if math.gcd(mu, mv) != 1: + temp = math.gcd(mu, mv) + mu = round(mu / temp) + mv = round(mv / temp) - # refer to the meaning of d in reference + # Refer to the meaning of d in reference d = (u**2 + v**2) * mv + w**2 * mu - # Compute the max n we need to enumerate. + # Compute the max n we need to enumerate n_max = int(np.sqrt((cutoff * 4 * mu * mv) / d)) - # Enumerate all possible n, m to give possible sigmas within the cutoff. + # Enumerate all possible n, m to give possible sigmas within the cutoff + sigmas: dict[int, list[float]] = {} for n in range(1, n_max + 1): m_max = 0 if c2_a2_ratio is None and w == 0 else int(np.sqrt((cutoff * 4 * mu * mv - n**2 * d) / mu)) for m in range(m_max + 1): - if gcd(m, n) == 1 or m == 0: - # construct the rotation matrix, refer to the reference + if math.gcd(m, n) == 1 or m == 0: + # Construct the rotation matrix, refer to the reference R_list = [ (u**2 * mv - v**2 * mv - w**2 * mu) * n**2 + mu * m**2, 2 * v * u * mv * n**2 - 2 * w * mu * m * n, @@ -1647,7 +1716,7 @@ def enum_sigma_tet(cutoff, r_axis, c2_a2_ratio): (w**2 * mu - u**2 * mv - v**2 * mv) * n**2 + mu * m**2, ] m = -1 * m - # inverse of rotation matrix + # Inverse of rotation matrix R_list_inv = [ (u**2 * mv - v**2 * mv - w**2 * mu) * n**2 + mu * m**2, 2 * v * u * mv * n**2 - 2 * w * mu * m * n, @@ -1663,9 +1732,9 @@ def enum_sigma_tet(cutoff, r_axis, c2_a2_ratio): F = mu * m**2 + d * n**2 all_list = R_list + R_list_inv + [F] # Compute the max common factors for the elements of the rotation matrix - # and its inverse. - com_fac = reduce(gcd, all_list) - sigma = int(round((mu * m**2 + d * n**2) / com_fac)) + # and its inverse + com_fac = reduce(math.gcd, all_list) + sigma = round((mu * m**2 + d * n**2) / com_fac) if 1 < sigma <= cutoff: if sigma not in list(sigmas): angle = 180.0 if m == 0 else 2 * np.arctan(n / m * np.sqrt(d / mu)) / np.pi * 180 @@ -1680,16 +1749,20 @@ def enum_sigma_tet(cutoff, r_axis, c2_a2_ratio): return sigmas @staticmethod - def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): + def enum_sigma_ort( + cutoff: int, + r_axis: Tuple3Ints, + c2_b2_a2_ratio: Tuple3Floats, + ) -> dict[int, list[float]]: """Find all possible sigma values and corresponding rotation angles within a sigma value cutoff with known rotation axis in orthorhombic system. - The algorithm for this code is from reference, Scipta Metallurgica 27, 291(1992). + + Reference: Scipta Metallurgica 27, 291(1992). Args: cutoff (int): the cutoff of sigma values. - r_axis (list of 3 integers, e.g. u, v, w): - the rotation axis of the grain boundary, with the format of [u,v,w]. - c2_b2_a2_ratio (list of 3 integers, e.g. mu,lambda, mv): + r_axis ((u, v, w)): the rotation axis of the grain boundary. + c2_b2_a2_ratio ((mu, lambda, mv)): mu:lam:mv is the square of the orthorhombic axial ratio with rational numbers. If irrational for one axis, set it to None. e.g. mu:lam:mv = c2,None,a2, means b2 is irrational. @@ -1705,23 +1778,23 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): angles, you need to analyze the symmetry of the structure. Different angles may result in equivalent microstructures. """ - sigmas = {} - # make sure gcd(r_axis)==1 - if reduce(gcd, r_axis) != 1: - r_axis = [int(round(x / reduce(gcd, r_axis))) for x in r_axis] + # Make sure math.gcd(r_axis) == 1 + if reduce(math.gcd, r_axis) != 1: + r_axis = cast(Tuple3Ints, tuple(round(x / reduce(math.gcd, r_axis)) for x in r_axis)) u, v, w = r_axis - # make sure mu, lambda, mv are coprime integers. + + # Make sure mu, lambda, mv are coprime integers if None in c2_b2_a2_ratio: mu, lam, mv = c2_b2_a2_ratio non_none = [i for i in c2_b2_a2_ratio if i is not None] if len(non_none) < 2: raise RuntimeError("No CSL exist for two irrational numbers") non1, non2 = non_none - if reduce(gcd, non_none) != 1: - temp = reduce(gcd, non_none) - non1 = int(round(non1 / temp)) - non2 = int(round(non2 / temp)) + if reduce(math.gcd, non_none) != 1: # type: ignore[arg-type] + temp = reduce(math.gcd, non_none) # type: ignore[arg-type] + non1 = round(non1 / temp) + non2 = round(non2 / temp) if mu is None: lam = non1 mv = non2 @@ -1742,23 +1815,24 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): raise RuntimeError("For irrational a2, CSL only exist for [1,0,0] or [0,v,w] and m = 0") else: mu, lam, mv = c2_b2_a2_ratio - if reduce(gcd, c2_b2_a2_ratio) != 1: - temp = reduce(gcd, c2_b2_a2_ratio) - mu = int(round(mu / temp)) - mv = int(round(mv / temp)) - lam = int(round(lam / temp)) + if reduce(math.gcd, c2_b2_a2_ratio) != 1: # type: ignore[arg-type] + temp = reduce(math.gcd, c2_b2_a2_ratio) # type: ignore[arg-type] + mu = round(mu / temp) + mv = round(mv / temp) + lam = round(lam / temp) if u == 0 and v == 0: mu = 1 if u == 0 and w == 0: lam = 1 if v == 0 and w == 0: mv = 1 - # refer to the meaning of d in reference + # Refer to the reference for the meaning of d d = (mv * u**2 + lam * v**2) * mv + w**2 * mu * mv - # Compute the max n we need to enumerate. + # Compute the max n we need to enumerate n_max = int(np.sqrt((cutoff * 4 * mu * mv * mv * lam) / d)) - # Enumerate all possible n, m to give possible sigmas within the cutoff. + # Enumerate all possible n, m to give possible sigmas within the cutoff + sigmas: dict[int, list[float]] = {} for n in range(1, n_max + 1): mu_temp, lam_temp, mv_temp = c2_b2_a2_ratio if (mu_temp is None and w == 0) or (lam_temp is None and v == 0) or (mv_temp is None and u == 0): @@ -1766,8 +1840,8 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): else: m_max = int(np.sqrt((cutoff * 4 * mu * mv * lam * mv - n**2 * d) / mu / lam)) for m in range(m_max + 1): - if gcd(m, n) == 1 or m == 0: - # construct the rotation matrix, refer to the reference + if math.gcd(m, n) == 1 or m == 0: + # Construct the rotation matrix, refer to the reference R_list = [ (u**2 * mv * mv - lam * v**2 * mv - w**2 * mu * mv) * n**2 + lam * mu * m**2, 2 * lam * (v * u * mv * n**2 - w * mu * m * n), @@ -1780,7 +1854,7 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): (w**2 * mu * mv - u**2 * mv * mv - v**2 * mv * lam) * n**2 + lam * mu * m**2, ] m = -1 * m - # inverse of rotation matrix + # Inverse of rotation matrix R_list_inv = [ (u**2 * mv * mv - lam * v**2 * mv - w**2 * mu * mv) * n**2 + lam * mu * m**2, 2 * lam * (v * u * mv * n**2 - w * mu * m * n), @@ -1796,9 +1870,9 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): F = mu * lam * m**2 + d * n**2 all_list = R_list + R_list_inv + [F] # Compute the max common factors for the elements of the rotation matrix - # and its inverse. - com_fac = reduce(gcd, all_list) - sigma = int(round((mu * lam * m**2 + d * n**2) / com_fac)) + # and its inverse. + com_fac = reduce(math.gcd, all_list) # type: ignore[arg-type] + sigma = round((mu * lam * m**2 + d * n**2) / com_fac) if 1 < sigma <= cutoff: if sigma not in list(sigmas): angle = 180.0 if m == 0 else 2 * np.arctan(n / m * np.sqrt(d / mu / lam)) / np.pi * 180 @@ -1813,26 +1887,30 @@ def enum_sigma_ort(cutoff, r_axis, c2_b2_a2_ratio): return sigmas @staticmethod - def enum_possible_plane_cubic(plane_cutoff, r_axis, r_angle): + def enum_possible_plane_cubic( + plane_cutoff: int, + r_axis: Tuple3Ints, + r_angle: float, + ) -> dict[Literal["Twist", "Symmetric tilt", "Normal tilt", "Mixed"], list[list]]: """Find all possible plane combinations for GBs given a rotation axis and angle for - cubic system, and classify them to different categories, including 'Twist', - 'Symmetric tilt', 'Normal tilt', 'Mixed' GBs. + cubic system, and classify them to different categories, including "Twist", + "Symmetric tilt", "Normal tilt", "Mixed" GBs. Args: plane_cutoff (int): the cutoff of plane miller index. - r_axis (list of 3 integers, e.g. u, v, w): - the rotation axis of the grain boundary, with the format of [u,v,w]. + r_axis ((u, v, w])): the rotation axis of the grain boundary. r_angle (float): rotation angle of the GBs. Returns: - dict: all combinations with keys as GB type, e.g. 'Twist','Symmetric tilt',etc. + dict: all combinations with keys as GB type, e.g. "Twist", "Symmetric tilt", etc. and values as the combination of the two plane miller index (GB plane and joining plane). """ - all_combinations = {} - all_combinations["Symmetric tilt"] = [] - all_combinations["Twist"] = [] - all_combinations["Normal tilt"] = [] - all_combinations["Mixed"] = [] + all_combinations: dict[Literal["Twist", "Symmetric tilt", "Normal tilt", "Mixed"], list[list]] = { + "Symmetric tilt": [], + "Twist": [], + "Normal tilt": [], + "Mixed": [], + } sym_plane = symm_group_cubic([[1, 0, 0], [1, 1, 0]]) j = np.arange(0, plane_cutoff + 1) combination = [] @@ -1851,7 +1929,7 @@ def enum_possible_plane_cubic(plane_cutoff, r_axis, r_angle): miller = np.array(combination) miller = miller[np.argsort(np.linalg.norm(miller, axis=1))] for val in miller: - if reduce(gcd, val) == 1: + if reduce(math.gcd, val) == 1: matrix = GrainBoundaryGenerator.get_trans_mat(r_axis, r_angle, surface=val, quick_gen=True) vec = np.cross(matrix[1][0], matrix[1][1]) miller2 = GrainBoundaryGenerator.vec_to_surface(vec) @@ -1878,12 +1956,17 @@ def enum_possible_plane_cubic(plane_cutoff, r_axis, r_angle): return all_combinations @staticmethod - def get_rotation_angle_from_sigma(sigma, r_axis, lat_type="C", ratio=None): - """Find all possible rotation angle for the given sigma value. + def get_rotation_angle_from_sigma( + sigma: int, + r_axis: Tuple3Ints | Tuple4Ints, + lat_type: str = "c", + ratio: tuple[int, int] | Tuple3Ints | None = None, + ) -> list[float]: + """Find all possible rotation angles for the given sigma value. Args: sigma (int): sigma value provided - r_axis (list of 3 integers, e.g. u, v, w or 4 integers, e.g. u, v, t, w for hex/rho system only): the + r_axis ((u, v, w) or (u, v, t, w) for hex/rho systems): the rotation axis of the grain boundary. lat_type (str): one character to specify the lattice type. Defaults to 'c' for cubic. 'c' or 'C': cubic system @@ -1891,17 +1974,17 @@ def get_rotation_angle_from_sigma(sigma, r_axis, lat_type="C", ratio=None): 'o' or 'O': orthorhombic system 'h' or 'H': hexagonal system 'r' or 'R': rhombohedral system - ratio (list of integers): lattice axial ratio. + ratio (tuple[int, ...]): lattice axial ratio. For cubic system, ratio is not needed. - For tetragonal system, ratio = [mu, mv], list of two integers, + For tetragonal system, ratio = (mu, mv), tuple of two integers, that is, mu/mv = c2/a2. If it is irrational, set it to none. - For orthorhombic system, ratio = [mu, lam, mv], list of 3 integers, + For orthorhombic system, ratio = (mu, lam, mv), tuple of 3 integers, that is, mu:lam:mv = c2:b2:a2. If irrational for one axis, set it to None. e.g. mu:lam:mv = c2,None,a2, means b2 is irrational. - For rhombohedral system, ratio = [mu, mv], list of two integers, + For rhombohedral system, ratio = (mu, mv), tuple of two integers, that is, mu/mv is the ratio of (1+2*cos(alpha)/cos(alpha). If irrational, set it to None. - For hexagonal system, ratio = [mu, mv], list of two integers, + For hexagonal system, ratio = (mu, mv), tuple of two integers, that is, mu/mv = c2/a2. If it is irrational, set it to none. Returns: @@ -1909,35 +1992,58 @@ def get_rotation_angle_from_sigma(sigma, r_axis, lat_type="C", ratio=None): If the sigma value is not correct, return the rotation angle corresponding to the correct possible sigma value right smaller than the wrong sigma value provided. """ - if lat_type.lower() == "c": + lat_type = lat_type.lower() + + # Check r_axis length + if lat_type in {"c", "t"}: + assert len(r_axis) == 3, "r_axis length incompatible with selected lattice system" + + # Check lattice axial ratio length + if lat_type == "o" and (ratio is None or len(ratio) != 3): + raise RuntimeError("Orthorhombic system needs correct c2:b2:a2 ratio") + + if lat_type in {"t", "h"} and (ratio is None or len(ratio) != 2): + raise RuntimeError("System needs correct c2/a2 ratio.") + + if lat_type == "r" and (ratio is None or len(ratio) != 2): + raise RuntimeError("Rhombohedral system needs correct (1+2*cos(alpha)/cos(alpha) ratio") + + if lat_type == "c": logger.info("Make sure this is for cubic system") - sigma_dict = GrainBoundaryGenerator.enum_sigma_cubic(cutoff=sigma, r_axis=r_axis) - elif lat_type.lower() == "t": + sigma_dict = GrainBoundaryGenerator.enum_sigma_cubic(cutoff=sigma, r_axis=cast(Tuple3Ints, r_axis)) + + elif lat_type == "t": logger.info("Make sure this is for tetragonal system") if ratio is None: logger.info("Make sure this is for irrational c2/a2 ratio") - elif len(ratio) != 2: - raise RuntimeError("Tetragonal system needs correct c2/a2 ratio") - sigma_dict = GrainBoundaryGenerator.enum_sigma_tet(cutoff=sigma, r_axis=r_axis, c2_a2_ratio=ratio) - elif lat_type.lower() == "o": + sigma_dict = GrainBoundaryGenerator.enum_sigma_tet( + cutoff=sigma, r_axis=cast(Tuple3Ints, r_axis), c2_a2_ratio=cast(tuple[int, int], ratio) + ) + + elif lat_type == "o": logger.info("Make sure this is for orthorhombic system") - if len(ratio) != 3: - raise RuntimeError("Orthorhombic system needs correct c2:b2:a2 ratio") - sigma_dict = GrainBoundaryGenerator.enum_sigma_ort(cutoff=sigma, r_axis=r_axis, c2_b2_a2_ratio=ratio) - elif lat_type.lower() == "h": + sigma_dict = GrainBoundaryGenerator.enum_sigma_ort( + cutoff=sigma, + r_axis=cast(Tuple3Ints, r_axis), + c2_b2_a2_ratio=cast(Tuple3Ints, ratio), + ) + + elif lat_type == "h": logger.info("Make sure this is for hexagonal system") if ratio is None: logger.info("Make sure this is for irrational c2/a2 ratio") - elif len(ratio) != 2: - raise RuntimeError("Hexagonal system needs correct c2/a2 ratio") - sigma_dict = GrainBoundaryGenerator.enum_sigma_hex(cutoff=sigma, r_axis=r_axis, c2_a2_ratio=ratio) - elif lat_type.lower() == "r": + sigma_dict = GrainBoundaryGenerator.enum_sigma_hex( + cutoff=sigma, r_axis=r_axis, c2_a2_ratio=cast(tuple[int, int], ratio) + ) + + elif lat_type == "r": logger.info("Make sure this is for rhombohedral system") if ratio is None: logger.info("Make sure this is for irrational (1+2*cos(alpha)/cos(alpha) ratio") - elif len(ratio) != 2: - raise RuntimeError("Rhombohedral system needs correct (1+2*cos(alpha)/cos(alpha) ratio") - sigma_dict = GrainBoundaryGenerator.enum_sigma_rho(cutoff=sigma, r_axis=r_axis, ratio_alpha=ratio) + sigma_dict = GrainBoundaryGenerator.enum_sigma_rho( + cutoff=sigma, r_axis=cast(Tuple3Ints, r_axis), ratio_alpha=cast(tuple[int, int], ratio) + ) + else: raise RuntimeError("Lattice type not implemented") @@ -1957,61 +2063,64 @@ def get_rotation_angle_from_sigma(sigma, r_axis, lat_type="C", ratio=None): return rotation_angles @staticmethod - def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=False): - """ - By linear operation of csl lattice vectors to get the best corresponding + def slab_from_csl( + csl: NDArray, + surface: Tuple3Ints, + normal: bool, + trans_cry: NDArray, + max_search: int = 20, + quick_gen: bool = False, + ) -> Matrix3D: + """By linear operation of csl lattice vectors to get the best corresponding slab lattice. That is the area of a,b vectors (within the surface plane) is the smallest, the c vector first, has shortest length perpendicular to surface [h,k,l], second, has shortest length itself. Args: - csl (3 by 3 integer array): - input csl lattice. - surface (list of 3 integers, e.g. h, k, l): - the miller index of the surface, with the format of [h,k,l] - normal (logic): - determine if the c vector needs to perpendicular to surface - trans_cry (3 by 3 array): - transform matrix from crystal system to orthogonal system + csl (3 by 3 integer array): input csl lattice. + surface ((h, k, l)): the miller index of the surface. + determine if the c vector needs to perpendicular to surface. + normal (bool): search the GB c vector that is perpendicular to the plane. + trans_cry (3 by 3 array): transform matrix from crystal system to orthogonal system. max_search (int): max search for the GB lattice vectors that give the smallest GB - lattice. If normal is true, also max search the GB c vector that perpendicular + lattice. If normal is True, also max search the GB c vector that is perpendicular to the plane. quick_gen (bool): whether to quickly generate a supercell, no need to find the smallest cell if set to true. Returns: - t_matrix: a slab lattice ( 3 by 3 integer array): + t_matrix: a slab lattice (3 by 3 integer array) """ - # set the transform matrix in real space + # Set the transform matrix in real space trans = trans_cry - # transform matrix in reciprocal space + # Transform matrix in reciprocal space ctrans = np.linalg.inv(trans.T) t_matrix = csl.copy() - # vectors constructed from csl that perpendicular to surface + # Vectors constructed from csl that perpendicular to surface ab_vector = [] - # obtain the miller index of surface in terms of csl. + # Obtain the miller index of surface in terms of csl miller = np.matmul(surface, csl.T) - if reduce(gcd, miller) != 1: - miller = [int(round(x / reduce(gcd, miller))) for x in miller] + if reduce(math.gcd, miller) != 1: + miller = [round(x / reduce(math.gcd, miller)) for x in miller] miller_nonzero = [] - # quickly generate a supercell, normal is not work in this way + # Quickly generate a supercell, normal is not working in this way if quick_gen: scale_factor = [] eye = np.eye(3, dtype=int) - for ii, jj in enumerate(miller): - if jj == 0: - scale_factor.append(eye[ii]) + for i, j in enumerate(miller): + if j == 0: + scale_factor.append(eye[i]) else: - miller_nonzero.append(ii) + miller_nonzero.append(i) if len(scale_factor) < 2: index_len = len(miller_nonzero) - for ii in range(index_len): - for jj in range(ii + 1, index_len): - lcm_miller = lcm(miller[miller_nonzero[ii]], miller[miller_nonzero[jj]]) + for i in range(index_len): + for j in range(i + 1, index_len): + lcm_miller = lcm(miller[miller_nonzero[i]], miller[miller_nonzero[j]]) scl_factor = [0, 0, 0] - scl_factor[miller_nonzero[ii]] = -int(round(lcm_miller / miller[miller_nonzero[ii]])) - scl_factor[miller_nonzero[jj]] = int(round(lcm_miller / miller[miller_nonzero[jj]])) + scl_factor[miller_nonzero[i]] = -round(lcm_miller / miller[miller_nonzero[i]]) + scl_factor[miller_nonzero[j]] = round(lcm_miller / miller[miller_nonzero[j]]) scale_factor.append(scl_factor) if len(scale_factor) == 2: break @@ -2023,22 +2132,22 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals return t_matrix c_index = 0 - for ii, jj in enumerate(miller): - if jj == 0: - ab_vector.append(csl[ii]) + for i, j in enumerate(miller): + if j == 0: + ab_vector.append(csl[i]) else: - c_index = ii - miller_nonzero.append(jj) + c_index = i + miller_nonzero.append(j) if len(miller_nonzero) > 1: t_matrix[2] = csl[c_index] index_len = len(miller_nonzero) lcm_miller = [] - for ii in range(index_len): - for jj in range(ii + 1, index_len): - com_gcd = gcd(miller_nonzero[ii], miller_nonzero[jj]) - mil1 = int(round(miller_nonzero[ii] / com_gcd)) - mil2 = int(round(miller_nonzero[jj] / com_gcd)) + for i in range(index_len): + for j in range(i + 1, index_len): + com_gcd = math.gcd(miller_nonzero[i], miller_nonzero[j]) + mil1 = round(miller_nonzero[i] / com_gcd) + mil2 = round(miller_nonzero[j] / com_gcd) lcm_miller.append(max(abs(mil1), abs(mil2))) lcm_sorted = sorted(lcm_miller) max_j = lcm_sorted[0] if index_len == 2 else lcm_sorted[1] @@ -2050,13 +2159,13 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals return t_matrix max_j = abs(miller_nonzero[0]) max_j = min(max_j, max_search) - # area of a, b vectors + # Area of a, b vectors area = None - # length of c vector + # Length of c vector c_norm = np.linalg.norm(np.matmul(t_matrix[2], trans)) # c vector length along the direction perpendicular to surface c_length = np.abs(np.dot(t_matrix[2], surface)) - # check if the init c vector perpendicular to the surface + # Check if the init c vector perpendicular to the surface if normal: c_cross = np.cross(np.matmul(t_matrix[2], trans), np.matmul(surface, ctrans)) normal_init = np.linalg.norm(c_cross) < 1e-8 @@ -2077,8 +2186,9 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals new_i = list(ii).copy() new_i[np.nonzero(ii)[0][0]] = -1 * new_i[np.nonzero(ii)[0][0]] combination.append(new_i) - for ii in combination: - if reduce(gcd, ii) == 1: + + for ii in combination: # type: ignore[assignment] + if reduce(math.gcd, ii) == 1: temp = np.dot(np.array(ii), csl) if abs(np.dot(temp, surface) - 0) < 1.0e-8: ab_vector.append(temp) @@ -2125,8 +2235,9 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals new_i = list(ii).copy() new_i[np.nonzero(ii)[0][0]] = -1 * new_i[np.nonzero(ii)[0][0]] combination.append(new_i) - for ii in combination: - if reduce(gcd, ii) == 1: + + for ii in combination: # type: ignore[assignment] + if reduce(math.gcd, ii) == 1: temp = np.dot(np.array(ii), csl) if abs(np.dot(temp, surface) - 0) > 1.0e-8: c_cross = np.cross(np.matmul(temp, trans), np.matmul(surface, ctrans)) @@ -2144,7 +2255,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals if normal_init: logger.info("Found perpendicular c vector") - # find the best a, b vectors with their formed area smallest and average norm of a,b smallest. + # Find the best a, b vectors with their formed area smallest and average norm of a,b smallest ab_norm = None for ii in combinations(ab_vector, 2): area_temp = np.linalg.norm(np.cross(np.matmul(ii[0], trans), np.matmul(ii[1], trans))) @@ -2162,7 +2273,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals area = area_temp ab_norm = ab_norm_temp - # make sure we have a left-handed crystallographic system + # Make sure we have a left-handed crystallographic system if np.linalg.det(np.matmul(t_matrix, trans)) < 0: t_matrix *= -1 @@ -2171,9 +2282,8 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals return t_matrix @staticmethod - def reduce_mat(mat, mag, r_matrix): - """ - Reduce integer array mat's determinant mag times by linear combination + def reduce_mat(mat: NDArray, mag: int, r_matrix: NDArray) -> NDArray: + """Reduce integer array mat's determinant mag times by linear combination of its row vectors, so that the new array after rotation (r_matrix) is still an integer array. @@ -2185,21 +2295,21 @@ def reduce_mat(mat, mag, r_matrix): Returns: the reduced integer array """ - max_j = abs(int(round(np.linalg.det(mat) / mag))) + max_j = abs(round(np.linalg.det(mat) / mag)) reduced = False for h in range(3): - kk = h + 1 if h + 1 < 3 else abs(2 - h) - ll = h + 2 if h + 2 < 3 else abs(1 - h) + kk = h + 1 if h < 2 else abs(2 - h) + ll = h + 2 if h < 1 else abs(1 - h) jj = np.arange(-max_j, max_j + 1) for j1, j2 in product(jj, repeat=2): temp = mat[h] + j1 * mat[kk] + j2 * mat[ll] if all(np.round(x, 5).is_integer() for x in list(temp / mag)): mat_copy = mat.copy() - mat_copy[h] = np.array([int(round(ele / mag)) for ele in temp]) + mat_copy[h] = np.array([round(ele / mag) for ele in temp]) new_mat = np.dot(mat_copy, np.linalg.inv(r_matrix.T)) if all(np.round(x, 5).is_integer() for x in list(np.ravel(new_mat))): reduced = True - mat[h] = np.array([int(round(ele / mag)) for ele in temp]) + mat[h] = np.array([round(ele / mag) for ele in temp]) break if reduced: break @@ -2209,23 +2319,23 @@ def reduce_mat(mat, mag, r_matrix): return mat @staticmethod - def vec_to_surface(vec): - """ - Transform a float vector to a surface miller index with integers. + def vec_to_surface(vec: Vector3D) -> MillerIndex: + """Transform a float vector to a surface miller index with integers. Args: - vec (1 by 3 array float vector): input float vector + vec (Vector3D): input float vector Returns: the surface miller index of the input vector. """ - miller = [None] * 3 + miller: list[None | int] = [None] * 3 index = [] for idx, value in enumerate(vec): if abs(value) < 1.0e-8: miller[idx] = 0 else: index.append(idx) + if len(index) == 1: miller[index[0]] = 1 else: @@ -2241,14 +2351,13 @@ def vec_to_surface(vec): else: com_lcm = lcm(frac[0].denominator, frac[1].denominator) miller[true_index] = com_lcm - miller[index[0]] = frac[0].numerator * int(round(com_lcm / frac[0].denominator)) - miller[index[1]] = frac[1].numerator * int(round(com_lcm / frac[1].denominator)) - return miller + miller[index[0]] = frac[0].numerator * round(com_lcm / frac[0].denominator) + miller[index[1]] = frac[1].numerator * round(com_lcm / frac[1].denominator) + return cast(Tuple3Ints, miller) -def fix_pbc(structure, matrix=None): - """ - Set all frac_coords of the input structure within [0,1]. +def fix_pbc(structure: Structure, matrix: NDArray = None) -> Structure: + """Wrap all frac_coords of the input structure within [0, 1]. Args: structure (pymatgen structure object): input structure @@ -2267,7 +2376,7 @@ def fix_pbc(structure, matrix=None): spec.append(site.specie) coord = np.array(site.frac_coords) for i in range(3): - coord[i] -= floor(coord[i]) + coord[i] -= math.floor(coord[i]) if np.allclose(coord[i], 1) or np.allclose(coord[i], 0): coord[i] = 0 else: @@ -2277,9 +2386,8 @@ def fix_pbc(structure, matrix=None): return Structure(latte, spec, coords, site_properties=structure.site_properties) -def symm_group_cubic(mat): - """ - Obtain cubic symmetric equivalents of the list of vectors. +def symm_group_cubic(mat: NDArray) -> list: + """Obtain cubic symmetric equivalents of the list of vectors. Args: mat (n by 3 array/matrix): lattice matrix @@ -2323,25 +2431,23 @@ def symm_group_cubic(mat): class Interface(Structure): - """This class stores data for defining an interface between two structures. - It is a subclass of pymatgen Structure. - """ + """Store data for defining an interface between two Structures.""" def __init__( self, - lattice, - species, - coords, - site_properties, - validate_proximity=False, - to_unit_cell=False, - coords_are_cartesian=False, + lattice: Lattice | NDArray, + species: list[Any], + coords: NDArray, + site_properties: dict[str, Any], + validate_proximity: bool = False, + to_unit_cell: bool = False, + coords_are_cartesian: bool = False, in_plane_offset: tuple[float, float] = (0, 0), gap: float = 0, vacuum_over_film: float = 0, interface_properties: dict | None = None, ) -> None: - """Make an interface structure, a structure object with additional information + """Make an Interface, a Structure with additional information and methods pertaining to interfaces. Args: @@ -2364,7 +2470,8 @@ def __init__( each species. validate_proximity (bool): Whether to check if there are sites that are less than 0.01 Ang apart. Defaults to False. - to_unit_cell (bool): Whether to translate sites into the unit cell. Defaults to False. + to_unit_cell (bool): Whether to translate sites into the unit cell. + Defaults to False. coords_are_cartesian (bool): Set to True if you are providing coordinates in Cartesian coordinates. Defaults to False. site_properties (dict): Properties associated with the sites as a @@ -2376,7 +2483,7 @@ def __init__( gap: gap between substrate and film in Angstroms; zero corresponds to the original distance between substrate and film sites vacuum_over_film: vacuum space above the film in Angstroms. Defaults to 0. - interface_properties: properties associated with the interface. Defaults to None. + interface_properties: properties associated with the Interface. Defaults to None. """ assert ( "interface_label" in site_properties @@ -2401,9 +2508,7 @@ def __init__( @property def in_plane_offset(self) -> np.ndarray: - """The shift between the film and substrate in fractional - coordinates. - """ + """The shift between the film and substrate in fractional coordinates.""" return self._in_plane_offset @in_plane_offset.setter @@ -2458,7 +2563,7 @@ def substrate_sites(self) -> list[Site]: @property def substrate(self) -> Structure: - """A pymatgen Structure for just the substrate.""" + """A Structure for just the substrate.""" return Structure.from_sites(self.substrate_sites) @property @@ -2473,15 +2578,19 @@ def film_sites(self) -> list[Site]: @property def film(self) -> Structure: - """A pymatgen Structure for just the film.""" + """A Structure for just the film.""" return Structure.from_sites(self.film_sites) def copy(self) -> Self: # type: ignore[override] """Make a copy of the Interface.""" - return Interface.from_dict(self.as_dict()) + return type(self).from_dict(self.as_dict()) - def get_sorted_structure(self, key=None, reverse=False) -> Structure: - """Get a sorted structure for the interface. The parameters have the same + def get_sorted_structure( + self, + key: Callable | None = None, + reverse: bool = False, + ) -> Structure: + """Get a sorted structure for the Interface. The parameters have the same meaning as in list.sort. By default, sites are sorted by the electronegativity of the species. @@ -2496,7 +2605,10 @@ def get_sorted_structure(self, key=None, reverse=False) -> Structure: struct_copy.sort(key=key, reverse=reverse) return struct_copy - def get_shifts_based_on_adsorbate_sites(self, tolerance: float = 0.1) -> list[tuple[float, float]]: + def get_shifts_based_on_adsorbate_sites( + self, + tolerance: float = 0.1, + ) -> list[tuple[float, float]]: """Compute possible in-plane shifts based on an adsorbate site algorithm. Args: @@ -2510,7 +2622,8 @@ def get_shifts_based_on_adsorbate_sites(self, tolerance: float = 0.1) -> list[tu substrate.lattice.inv_matrix, ) - # Film gets forced into substrate lattice anyways, so shifts can be computed in fractional coords + # Film gets forced into substrate lattice anyways, so shifts can be + # computed in fractional coords film_surface_sites = np.dot( list(chain.from_iterable(AdsorbateSiteFinder(film).find_adsorption_sites().values())), film.lattice.inv_matrix, @@ -2522,7 +2635,7 @@ def get_shifts_based_on_adsorbate_sites(self, tolerance: float = 0.1) -> list[tu ] ) - def _base_round(x, base=0.05): + def _base_round(x, base: float = 0.05): return base * (np.array(x) / base).round() # Round shifts to tolerance @@ -2560,8 +2673,11 @@ def substrate_layers(self) -> int: return count_layers(self.substrate, sorted_element_list[0][0]) def _update_c(self, new_c: float) -> None: - """Modifies the c-direction of the lattice without changing the site Cartesian coordinates - Be careful you can mess up the interface by setting a c-length that can't accommodate all the sites. + """Modify the c-direction of the lattice without changing the site + Cartesian coordinates. + + WARNING: you can mess up the Interface by setting a c-length that + can't accommodate all the sites. """ if new_c <= 0: raise ValueError("New c-length must be greater than 0") @@ -2574,14 +2690,15 @@ def _update_c(self, new_c: float) -> None: site._lattice = new_lattice # Update the lattice site.coords = c_coords # Put back into original Cartesian space - def as_dict(self): + def as_dict(self) -> dict: # type: ignore[override] """MSONable dict.""" - dct = super().as_dict() - dct["in_plane_offset"] = self.in_plane_offset.tolist() - dct["gap"] = self.gap - dct["vacuum_over_film"] = self.vacuum_over_film - dct["interface_properties"] = self.interface_properties - return dct + return { + **super().as_dict(), + "in_plane_offset": self.in_plane_offset.tolist(), + "gap": self.gap, + "vacuum_over_film": self.vacuum_over_film, + "interface_properties": self.interface_properties, + } @classmethod def from_dict(cls, dct: dict) -> Self: # type: ignore[override] @@ -2602,7 +2719,7 @@ def from_dict(cls, dct: dict) -> Self: # type: ignore[override] "vacuum_over_film": dct.get("vacuum_over_film"), "interface_properties": dct.get("interface_properties"), } - return Interface( + return cls( lattice=lattice, species=struct.species_and_occu, coords=struct.frac_coords, @@ -2621,26 +2738,26 @@ def from_slabs( interface_properties: dict | None = None, center_slab: bool = True, ) -> Self: - """Make an interface structure by merging a substrate and film slabs + """Make an Interface by merging a substrate and film slabs The film a- and b-vectors will be forced to be the substrate slab's a- and b-vectors. For now, it's suggested to use a factory method that will ensure the - appropriate interface structure is already met. + appropriate Interface is already met. Args: - substrate_slab (Slab): slab for the substrate - film_slab (Slab): slab for the film - in_plane_offset (tuple): fractional shift in plane for the film with respect to the substrate. - For example, (0.5, 0.5) will shift the film by half the substrate's a- and b-vectors. - Defaults to (0, 0). + substrate_slab (Slab): slab for the substrate. + film_slab (Slab): slab for the film. + in_plane_offset (tuple): fractional shift in plane for the film with + respect to the substrate. For example, (0.5, 0.5) will shift the + film by half the substrate's a- and b-vectors. Defaults to (0, 0). gap (float): gap between substrate and film in Angstroms. Defaults to 1.6. - vacuum_over_film (float): vacuum space above the film in Angstroms. Defaults to 0. - interface_properties (dict): misc properties to assign to the interface. Defaults to None. + vacuum_over_film (float): vacuum space above the film in Angstroms. + Defaults to 0. + interface_properties (dict): misc properties to assign to the Interface. + Defaults to None. center_slab (bool): center the slab. Defaults to True. """ - interface_properties = interface_properties or {} - # Ensure c-axis is orthogonal to a/b plane if isinstance(substrate_slab, Slab): substrate_slab = substrate_slab.get_orthogonal_c_slab() @@ -2667,7 +2784,7 @@ def from_slabs( film_max_c = np.max(film_coords[:, 2]) * film_slab.lattice.c min_height = np.abs(film_max_c - film_min_c) + np.abs(sub_max_c - sub_min_c) - # construct new lattice + # Construct new lattice abc = substrate_slab.lattice.abc[:2] + (min_height + gap + vacuum_over_film,) angles = substrate_slab.lattice.angles lattice = Lattice.from_parameters(*abc, *angles) @@ -2715,7 +2832,7 @@ def from_slabs( in_plane_offset=in_plane_offset, gap=gap, vacuum_over_film=vacuum_over_film, - interface_properties=interface_properties, + interface_properties=interface_properties or {}, ) iface.sort() @@ -2723,7 +2840,7 @@ def from_slabs( def label_termination(slab: Structure) -> str: - """Labels the slab surface termination.""" + """Label the slab surface termination.""" frac_coords = slab.frac_coords n = len(frac_coords) @@ -2762,8 +2879,8 @@ def label_termination(slab: Structure) -> str: return f"{form}_{sp_symbol}_{len(top_plane)}" -def count_layers(struct: Structure, el=None) -> int: - """Counts the number of 'layers' along the c-axis.""" +def count_layers(struct: Structure, el: Element | None = None) -> int: + """Count the number of layers along the c-axis.""" el = el or struct.elements[0] frac_coords = [site.frac_coords for site in struct if site.species_string == str(el)] n_el_sites = len(frac_coords) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 61bf6dcd196..7aedc4e9014 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -35,6 +35,7 @@ from pymatgen.symmetry.analyzer import SpacegroupAnalyzer from pymatgen.util.coord import in_coord_list from pymatgen.util.due import Doi, due +from pymatgen.util.typing import Tuple3Ints if TYPE_CHECKING: from collections.abc import Sequence @@ -939,7 +940,7 @@ def __init__( def reduce_vector(vector: MillerIndex) -> MillerIndex: """Helper function to reduce vectors.""" divisor = abs(reduce(gcd, vector)) # type: ignore[arg-type] - return cast(tuple[int, int, int], tuple(int(idx / divisor) for idx in vector)) + return cast(Tuple3Ints, tuple(int(idx / divisor) for idx in vector)) def add_site_types() -> None: """Add Wyckoff symbols and equivalent sites to the initial structure.""" @@ -1985,7 +1986,7 @@ def get_symmetrically_equivalent_miller_indices( else: symm_ops = structure.lattice.get_recp_symmetry_operation() - equivalent_millers: list[tuple[int, int, int]] = [_miller_index] + equivalent_millers: list[Tuple3Ints] = [_miller_index] for miller in itertools.product(idx_range, idx_range, idx_range): if miller == _miller_index: continue @@ -2038,7 +2039,7 @@ def get_symmetrically_distinct_miller_indices( spg_analyzer = SpacegroupAnalyzer(structure) if spg_analyzer.get_crystal_system() == "trigonal": transf = spg_analyzer.get_conventional_to_primitive_transformation_matrix() - miller_list: list[tuple[int, int, int]] = [hkl_transformation(transf, hkl) for hkl in conv_hkl_list] + miller_list: list[Tuple3Ints] = [hkl_transformation(transf, hkl) for hkl in conv_hkl_list] prim_structure = SpacegroupAnalyzer(structure).get_primitive_standard_structure() symm_ops = prim_structure.lattice.get_recp_symmetry_operation() @@ -2051,7 +2052,7 @@ def get_symmetrically_distinct_miller_indices( for idx, miller in enumerate(miller_list): denom = abs(reduce(gcd, miller)) # type: ignore[arg-type] - miller = cast(tuple[int, int, int], tuple(int(idx / denom) for idx in miller)) + miller = cast(Tuple3Ints, tuple(int(idx / denom) for idx in miller)) if not _is_in_miller_family(miller, unique_millers, symm_ops): if spg_analyzer.get_crystal_system() == "trigonal": # Now we find the distinct primitive hkls using @@ -2091,7 +2092,7 @@ def _is_in_miller_family( def hkl_transformation( transf: np.ndarray, miller_index: MillerIndex, -) -> tuple[int, int, int]: +) -> Tuple3Ints: """Transform the Miller index from setting A to B with a transformation matrix. Args: @@ -2125,7 +2126,7 @@ def miller_index_from_sites( coords_are_cartesian: bool = True, round_dp: int = 4, verbose: bool = True, -) -> tuple[int, int, int]: +) -> Tuple3Ints: """Get the Miller index of a plane, determined by a given set of coordinates. A minimum of 3 sets of coordinates are required. If more than 3 diff --git a/pymatgen/core/xcfunc.py b/pymatgen/core/xcfunc.py index 4264175cddf..a17d9c939fb 100644 --- a/pymatgen/core/xcfunc.py +++ b/pymatgen/core/xcfunc.py @@ -81,8 +81,8 @@ class type_name(NamedTuple): name: str xcf = LibxcFunc - defined_aliases: ClassVar = { - # (x, c) --> type_name + defined_aliases: ClassVar[dict[tuple[LibxcFunc, LibxcFunc], type_name]] = { + # (X, C) --> type_name # LDAs (xcf.LDA_X, xcf.LDA_C_PW): type_name("LDA", "PW"), # ixc 7 (xcf.LDA_X, xcf.LDA_C_PW_MOD): type_name("LDA", "PW_MOD"), @@ -101,13 +101,11 @@ class type_name(NamedTuple): (xcf.GGA_X_B88, xcf.GGA_C_LYP): type_name("GGA", "BLYP"), } - del type_name - # Correspondence between Abinit ixc notation and libxc notation. # see: http://www.abinit.org/doc/helpfiles/for-v7.8/input_variables/varbas.html#ixc # and 42_libpaw/m_pawpsp.F90 for the implementation. # Fortunately, all the other cases are handled with libxc. - abinitixc_to_libxc: ClassVar = { + abinitixc_to_libxc: ClassVar[dict[int, dict[str, LibxcFunc]]] = { 1: {"xc": xcf.LDA_XC_TETER93}, 2: {"x": xcf.LDA_X, "c": xcf.LDA_C_PZ}, # PZ 001009 4: {"x": xcf.LDA_X, "c": xcf.LDA_C_WIGNER}, # W @@ -118,8 +116,6 @@ class type_name(NamedTuple): 15: {"x": xcf.GGA_X_RPBE, "c": xcf.GGA_C_PBE}, # RPBE } - del xcf - def __init__( self, xc: LibxcFunc | None = None, diff --git a/pymatgen/electronic_structure/dos.py b/pymatgen/electronic_structure/dos.py index 77e7cf874a7..be074213ca6 100644 --- a/pymatgen/electronic_structure/dos.py +++ b/pymatgen/electronic_structure/dos.py @@ -24,7 +24,7 @@ from typing_extensions import Self from pymatgen.core.sites import PeriodicSite - from pymatgen.util.typing import SpeciesLike + from pymatgen.util.typing import SpeciesLike, Tuple3Floats class DOS(Spectrum): @@ -256,9 +256,7 @@ def get_interpolated_value(self, energy: float) -> dict[Spin, float]: energies[spin] = get_linear_interpolated_value(self.energies, self.densities[spin], energy) return energies - def get_interpolated_gap( - self, tol: float = 0.001, abs_tol: bool = False, spin: Spin | None = None - ) -> tuple[float, float, float]: + def get_interpolated_gap(self, tol: float = 0.001, abs_tol: bool = False, spin: Spin | None = None) -> Tuple3Floats: """Expects a DOS object and finds the gap. Args: diff --git a/pymatgen/ext/optimade.py b/pymatgen/ext/optimade.py index 0144d0253b8..6fbd1d18a05 100644 --- a/pymatgen/ext/optimade.py +++ b/pymatgen/ext/optimade.py @@ -58,7 +58,7 @@ class OptimadeRester: # regenerate on-demand from official providers.json using OptimadeRester.refresh_aliases() # these aliases are provided as a convenient shortcut for users of the OptimadeRester class - aliases: ClassVar = { + aliases: ClassVar[dict[str, str]] = { "aflow": "http://aflow.org/API/optimade/", "alexandria": "https://alexandria.icams.rub.de/pbe", "alexandria.pbe": "https://alexandria.icams.rub.de/pbe", diff --git a/pymatgen/io/abinit/abiobjects.py b/pymatgen/io/abinit/abiobjects.py index 47af77e2651..8be6d95af2d 100644 --- a/pymatgen/io/abinit/abiobjects.py +++ b/pymatgen/io/abinit/abiobjects.py @@ -17,7 +17,7 @@ from pymatgen.core import ArrayWithUnit, Lattice, Species, Structure, units if TYPE_CHECKING: - from typing import ClassVar + from typing import Any, ClassVar from typing_extensions import Self @@ -404,7 +404,7 @@ class Smearing(AbivarAble, MSONable): """ # Map string_mode to occopt - _mode2occopt: ClassVar = dict( + _mode2occopt: ClassVar[dict[str, int]] = dict( nosmearing=1, fermi_dirac=3, marzari4=4, @@ -509,8 +509,8 @@ def from_dict(cls, dct: dict) -> Self: class ElectronsAlgorithm(dict, AbivarAble, MSONable): """Variables controlling the SCF/NSCF algorithm.""" - # None indicates that we use abinit defaults. - _DEFAULT: ClassVar = dict( + # None indicates that we use abinit defaults + _DEFAULT: ClassVar[dict[str, int | None]] = dict( iprcell=None, iscf=None, diemac=None, @@ -1035,7 +1035,7 @@ class RelaxationMethod(AbivarAble, MSONable): The set of variables are constructed in to_abivars depending on ionmov and optcell. """ - _default_vars: ClassVar = dict( + _default_vars: ClassVar[dict[str, Any]] = dict( ionmov=MANDATORY, optcell=MANDATORY, ntime=80, @@ -1337,10 +1337,14 @@ class Screening(AbivarAble): """ # Approximations used for W - _WTYPES: ClassVar = dict(RPA=0) + _WTYPES: ClassVar[dict[str, int]] = dict(RPA=0) - # Self-consistecy modes - _SC_MODES: ClassVar = dict(one_shot=0, energy_only=1, wavefunctions=2) + # Self-consistency modes + _SC_MODES: ClassVar[dict[str, int]] = dict( + one_shot=0, + energy_only=1, + wavefunctions=2, + ) def __init__( self, @@ -1424,7 +1428,7 @@ def to_abivars(self): class SelfEnergy(AbivarAble): """Define the parameters used for the computation of the self-energy.""" - _SIGMA_TYPES: ClassVar = dict( + _SIGMA_TYPES: ClassVar[dict[str, int]] = dict( gw=0, hartree_fock=5, sex=6, @@ -1433,7 +1437,7 @@ class SelfEnergy(AbivarAble): model_gw_cd=9, ) - _SC_MODES: ClassVar = dict( + _SC_MODES: ClassVar[dict[str, int]] = dict( one_shot=0, energy_only=1, wavefunctions=2, @@ -1558,14 +1562,18 @@ class ExcHamiltonian(AbivarAble): """Contain parameters for the solution of the Bethe-Salpeter equation.""" # Types of excitonic Hamiltonian. - _EXC_TYPES: ClassVar = dict( + _EXC_TYPES: ClassVar[dict[str, int]] = dict( TDA=0, # Tamm-Dancoff approximation. coupling=1, # Calculation with coupling. ) # Algorithms used to compute the macroscopic dielectric function # and/or the exciton wavefunctions. - _ALGO2VAR: ClassVar = dict(direct_diago=1, haydock=2, cg=3) + _ALGO2VAR: ClassVar[dict[str, int]] = dict( + direct_diago=1, + haydock=2, + cg=3, + ) # Options specifying the treatment of the Coulomb term. _COULOMB_MODES = ("diago", "full", "model_df") diff --git a/pymatgen/io/abinit/pseudos.py b/pymatgen/io/abinit/pseudos.py index 40e49fca580..a3f29c60940 100644 --- a/pymatgen/io/abinit/pseudos.py +++ b/pymatgen/io/abinit/pseudos.py @@ -675,7 +675,7 @@ def _int_from_str(string): class NcAbinitHeader(AbinitHeader): """The abinit header found in the NC pseudopotential files.""" - _VARS: ClassVar = dict( + _VARS: ClassVar[dict[str, tuple]] = dict( zatom=(None, _int_from_str), zion=(None, float), pspdat=(None, float), @@ -866,7 +866,7 @@ def tm_header(filename, ppdesc): class PawAbinitHeader(AbinitHeader): """The abinit header found in the PAW pseudopotential files.""" - _VARS: ClassVar = dict( + _VARS: ClassVar[dict[str, tuple]] = dict( zatom=(None, _int_from_str), zion=(None, float), pspdat=(None, float), @@ -1012,7 +1012,7 @@ class ppdesc(NamedTuple): format: None # TODO Recheck - _PSPCODES: ClassVar = { + _PSPCODES: ClassVar[dict[int, ppdesc]] = { 1: ppdesc(1, "TM", "NC", None), 2: ppdesc(2, "GTH", "NC", None), 3: ppdesc(3, "HGH", "NC", None), @@ -1024,8 +1024,6 @@ class ppdesc(NamedTuple): 10: ppdesc(10, "HGHK", "NC", None), } - del ppdesc - # renumber functionals from oncvpsp todo confirm that 3 is 2 # _FUNCTIONALS = {1: {'n': 4, 'name': 'Wigner'}, # 2: {'n': 5, 'name': 'HL'}, diff --git a/pymatgen/io/adf.py b/pymatgen/io/adf.py index e61d7122a8a..e302530505c 100644 --- a/pymatgen/io/adf.py +++ b/pymatgen/io/adf.py @@ -370,7 +370,7 @@ class AdfTask(MSONable): ADF does not support calculating force/gradient. """ - operations: ClassVar = dict( + operations: ClassVar[dict[str, str]] = dict( energy="Evaluate the single point energy.", optimize="Minimize the energy by varying the molecular structure.", frequencies="Compute second derivatives and print out an analysis of molecular vibrations.", diff --git a/pymatgen/io/aims/inputs.py b/pymatgen/io/aims/inputs.py index 64877178e73..fda6aba8a75 100644 --- a/pymatgen/io/aims/inputs.py +++ b/pymatgen/io/aims/inputs.py @@ -25,6 +25,8 @@ from typing_extensions import Self + from pymatgen.util.typing import Tuple3Floats, Tuple3Ints + __author__ = "Thomas A. R. Purcell" __version__ = "1.0" __email__ = "purcellt@arizona.edu" @@ -252,9 +254,9 @@ class AimsCube(MSONable): """ type: str = field(default_factory=str) - origin: Sequence[float] | tuple[float, float, float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + origin: Sequence[float] | Tuple3Floats = field(default_factory=lambda: [0.0, 0.0, 0.0]) edges: Sequence[Sequence[float]] = field(default_factory=lambda: 0.1 * np.eye(3)) - points: Sequence[int] | tuple[int, int, int] = field(default_factory=lambda: [0, 0, 0]) + points: Sequence[int] | Tuple3Ints = field(default_factory=lambda: [0, 0, 0]) format: str = "cube" spin_state: int | None = None kpoint: int | None = None diff --git a/pymatgen/io/aims/parsers.py b/pymatgen/io/aims/parsers.py index 44ff922a026..f06f1307893 100644 --- a/pymatgen/io/aims/parsers.py +++ b/pymatgen/io/aims/parsers.py @@ -11,6 +11,7 @@ from pymatgen.core import Lattice, Molecule, Structure from pymatgen.core.tensors import Tensor +from pymatgen.util.typing import Tuple3Floats if TYPE_CHECKING: from collections.abc import Generator, Sequence @@ -567,9 +568,9 @@ def _parse_lattice_atom_pos( elif "atom " in line: line_split = line.split() species.append(line_split[4]) - coords.append(cast(tuple[float, float, float], tuple(float(inp) for inp in line_split[1:4]))) + coords.append(cast(Tuple3Floats, tuple(float(inp) for inp in line_split[1:4]))) elif "velocity " in line: - velocities.append(cast(tuple[float, float, float], tuple(float(inp) for inp in line.split()[1:4]))) + velocities.append(cast(Tuple3Floats, tuple(float(inp) for inp in line.split()[1:4]))) lattice = Lattice(lattice_vectors) if len(lattice_vectors) == 3 else None return species, coords, velocities, lattice diff --git a/pymatgen/io/cp2k/inputs.py b/pymatgen/io/cp2k/inputs.py index 081f16525b8..f6c47644dc2 100644 --- a/pymatgen/io/cp2k/inputs.py +++ b/pymatgen/io/cp2k/inputs.py @@ -52,7 +52,7 @@ from pymatgen.core.lattice import Lattice from pymatgen.core.structure import Molecule, Structure - from pymatgen.util.typing import Kpoint + from pymatgen.util.typing import Kpoint, Tuple3Ints __author__ = "Nicholas Winner" __version__ = "2.0" @@ -1786,7 +1786,7 @@ def f3(x): nel_beta = [] n_alpha = [] n_beta = [] - unpaired_orbital: tuple[int, int, int] = (0, 0, 0) + unpaired_orbital: Tuple3Ints = (0, 0, 0) while tmp: tmp2 = -min((esv[0][2], tmp)) if tmp > 0 else min((f2(esv[0][1]) - esv[0][2], -tmp)) l_alpha.append(esv[0][1]) diff --git a/pymatgen/io/fiesta.py b/pymatgen/io/fiesta.py index f4f98bc8af7..c98b3590f34 100644 --- a/pymatgen/io/fiesta.py +++ b/pymatgen/io/fiesta.py @@ -25,6 +25,8 @@ from typing_extensions import Self + from pymatgen.util.typing import Tuple3Ints + __author__ = "ndardenne" __copyright__ = "Copyright 2012, The Materials Project" __version__ = "0.1" @@ -103,9 +105,7 @@ class FiestaRun(MSONable): otherwise it breaks. """ - def __init__( - self, folder: str | None = None, grid: tuple[int, int, int] = (2, 2, 2), log_file: str = "log" - ) -> None: + def __init__(self, folder: str | None = None, grid: Tuple3Ints = (2, 2, 2), log_file: str = "log") -> None: """ Args: folder: Folder to look for runs. diff --git a/pymatgen/io/lobster/inputs.py b/pymatgen/io/lobster/inputs.py index 3c89e7aaab1..654dfd287c6 100644 --- a/pymatgen/io/lobster/inputs.py +++ b/pymatgen/io/lobster/inputs.py @@ -35,6 +35,7 @@ from typing_extensions import Self from pymatgen.core.composition import Composition + from pymatgen.util.typing import Tuple3Ints __author__ = "Janine George, Marco Esters" @@ -437,7 +438,7 @@ def write_KPOINTS( reciprocal_density: int = 100, isym: int = -1, from_grid: bool = False, - input_grid: tuple[int, int, int] = (5, 5, 5), + input_grid: Tuple3Ints = (5, 5, 5), line_mode: bool = True, kpoints_line_density: int = 20, symprec: float = 0.01, diff --git a/pymatgen/io/nwchem.py b/pymatgen/io/nwchem.py index 5063405f0e2..e070934bcc5 100644 --- a/pymatgen/io/nwchem.py +++ b/pymatgen/io/nwchem.py @@ -48,7 +48,7 @@ class NwTask(MSONable): """Base task for Nwchem.""" - theories: ClassVar = { + theories: ClassVar[dict[str, str]] = { "g3gn": "some description", "scf": "Hartree-Fock", "dft": "DFT", @@ -69,7 +69,7 @@ class NwTask(MSONable): "tddft": "Time Dependent DFT", } - operations: ClassVar = { + operations: ClassVar[dict[str, str]] = { "energy": "Evaluate the single point energy.", "gradient": "Evaluate the derivative of the energy with respect to nuclear coordinates.", "optimize": "Minimize the energy by varying the molecular structure.", diff --git a/pymatgen/io/pwmat/outputs.py b/pymatgen/io/pwmat/outputs.py index 519132139c6..8013be25dc5 100644 --- a/pymatgen/io/pwmat/outputs.py +++ b/pymatgen/io/pwmat/outputs.py @@ -13,7 +13,7 @@ if TYPE_CHECKING: from pymatgen.core import Structure - from pymatgen.util.typing import PathLike + from pymatgen.util.typing import PathLike, Tuple3Ints __author__ = "Hanyu Liu" @@ -190,7 +190,7 @@ def __init__(self, filename: PathLike): self._eigenvalues = self._parse_eigen() self._kpts, self._kpts_weight, self._hsps = self._parse_kpt() - def _parse_band(self) -> tuple[int, int, int]: + def _parse_band(self) -> Tuple3Ints: """Parse REPORT file to obtain spin switches, the number of kpoints and the number of bands. diff --git a/pymatgen/io/pwscf.py b/pymatgen/io/pwscf.py index 721f5331e9f..49f281569d2 100644 --- a/pymatgen/io/pwscf.py +++ b/pymatgen/io/pwscf.py @@ -501,7 +501,7 @@ class PWInputError(BaseException): class PWOutput: """Parser for PWSCF output file.""" - patterns: ClassVar = dict( + patterns: ClassVar[dict[str, str]] = dict( energies=r"total energy\s+=\s+([\d\.\-]+)\sRy", ecut=r"kinetic\-energy cutoff\s+=\s+([\d\.\-]+)\s+Ry", lattice_type=r"bravais\-lattice index\s+=\s+(\d+)", diff --git a/pymatgen/io/res.py b/pymatgen/io/res.py index d2859d24a1d..0c0c2e38842 100644 --- a/pymatgen/io/res.py +++ b/pymatgen/io/res.py @@ -30,7 +30,7 @@ from typing_extensions import Self - from pymatgen.util.typing import Vector3D + from pymatgen.util.typing import Tuple3Ints, Vector3D @dataclass(frozen=True) @@ -488,7 +488,7 @@ def get_cut_grid_gmax_fsbc(self) -> tuple[float, float, float, str] | None: self._raise_or_none(ResParseError("Could not find line with cut-off energy.")) return None - def get_mpgrid_offset_nkpts_spacing(self) -> tuple[tuple[int, int, int], Vector3D, int, float] | None: + def get_mpgrid_offset_nkpts_spacing(self) -> tuple[Tuple3Ints, Vector3D, int, float] | None: """ Retrieves the MP grid, the grid offsets, number of kpoints, and maximum kpoint spacing. diff --git a/pymatgen/io/vasp/inputs.py b/pymatgen/io/vasp/inputs.py index ab9bf4f4abb..f570e28b2d5 100644 --- a/pymatgen/io/vasp/inputs.py +++ b/pymatgen/io/vasp/inputs.py @@ -37,7 +37,7 @@ from pymatgen.electronic_structure.core import Magmom from pymatgen.util.io_utils import clean_lines from pymatgen.util.string import str_delimited -from pymatgen.util.typing import Kpoint, Vector3D +from pymatgen.util.typing import Kpoint, Tuple3Floats, Tuple3Ints, Vector3D if TYPE_CHECKING: from collections.abc import Iterator @@ -1194,7 +1194,7 @@ def kpts(self) -> Sequence[Kpoint]: return cast(Sequence[Kpoint], list(map(tuple, self._kpts))) # type: ignore[arg-type] if all(isinstance(point, (int, float)) for point in self._kpts) and len(self._kpts) == 3: - return [cast(tuple[float, float, float], tuple(self._kpts))] + return [cast(Kpoint, tuple(self._kpts))] raise ValueError(f"Invalid Kpoint {self._kpts}.") @@ -1316,9 +1316,7 @@ def automatic_density(cls, structure: Structure, kppa: float, force_gamma: bool ngrid = kppa / len(structure) mult: float = (ngrid * lengths[0] * lengths[1] * lengths[2]) ** (1 / 3) - num_div: tuple[int, int, int] = cast( - tuple[int, int, int], [math.floor(max(mult / length, 1)) for length in lengths] - ) + num_div: Tuple3Ints = cast(Tuple3Ints, [math.floor(max(mult / length, 1)) for length in lengths]) is_hexagonal: bool = lattice.is_hexagonal() is_face_centered: bool = structure.get_space_group_info()[0][0] == "F" @@ -1371,7 +1369,7 @@ def automatic_gamma_density(cls, structure: Structure, kppa: float) -> Self: comment, n_kpts, style, - [cast(tuple[int, int, int], tuple(n_div))], + [cast(Tuple3Ints, tuple(n_div))], (0, 0, 0), ) @@ -1424,7 +1422,7 @@ def automatic_density_by_lengths( lattice = structure.lattice abc = lattice.abc - num_div: tuple[int, int, int] = tuple(np.ceil(ld / abc[idx]) for idx, ld in enumerate(length_densities)) + num_div: Tuple3Ints = tuple(np.ceil(ld / abc[idx]) for idx, ld in enumerate(length_densities)) is_hexagonal: bool = lattice.is_hexagonal() is_face_centered: bool = structure.get_space_group_info()[0][0] == "F" @@ -1528,7 +1526,7 @@ def from_str(cls, string: str) -> Self: _kpt: list[float] = [float(i) for i in lines[3].split()] if len(_kpt) != 3: raise ValueError("Invalid Kpoint length.") - kpt: tuple[float, float, float] = cast(tuple[float, float, float], tuple(_kpt)) + kpt: Tuple3Floats = cast(Tuple3Floats, tuple(_kpt)) kpts_shift: Vector3D = (0, 0, 0) if len(lines) > 4 and coord_pattern.match(lines[4]): @@ -1568,7 +1566,7 @@ def from_str(cls, string: str) -> Self: "Cartesian" if lines[3].lower()[0] in "ck" else "Reciprocal" ) _style = cls.supported_modes.Line_mode - _kpts: list[tuple[float, float, float]] = [] + _kpts: list[Kpoint] = [] labels = [] patt = re.compile(r"([e0-9.\-]+)\s+([e0-9.\-]+)\s+([e0-9.\-]+)\s*!*\s*(.*)") for idx in range(4, len(lines)): @@ -1597,7 +1595,7 @@ def from_str(cls, string: str) -> Self: for idx in range(3, 3 + num_kpts): tokens = lines[idx].split() - kpts.append(cast(tuple[float, float, float], tuple(float(i) for i in tokens[:3]))) + kpts.append(cast(Kpoint, tuple(float(i) for i in tokens[:3]))) kpts_weights.append(float(tokens[3])) if len(tokens) > 4: labels.append(tokens[4]) diff --git a/pymatgen/io/vasp/outputs.py b/pymatgen/io/vasp/outputs.py index 6d1e6a824ec..da8ed8efec1 100644 --- a/pymatgen/io/vasp/outputs.py +++ b/pymatgen/io/vasp/outputs.py @@ -42,6 +42,7 @@ from pymatgen.io.wannier90 import Unk from pymatgen.util.io_utils import clean_lines, micro_pyawk from pymatgen.util.num import make_symmetric_matrix_from_upper_tri +from pymatgen.util.typing import Kpoint, Tuple3Floats, Vector3D if TYPE_CHECKING: from typing import Any, Callable, Literal @@ -1411,7 +1412,7 @@ def parse_atomic_symbol(symbol: str) -> str: return [parse_atomic_symbol(sym) for sym in atomic_symbols], potcar_symbols @staticmethod - def _parse_kpoints(elem: XML_Element) -> tuple[Kpoints, list[tuple[float, float, float]], list[float]]: + def _parse_kpoints(elem: XML_Element) -> tuple[Kpoints, list[Tuple3Floats], list[float]]: """Parse Kpoints.""" e = elem if elem.find("generation") is None else elem.find("generation") kpoint = Kpoints("Kpoints from vasprun.xml") @@ -1423,10 +1424,10 @@ def _parse_kpoints(elem: XML_Element) -> tuple[Kpoints, list[tuple[float, float, if name == "divisions": kpoint.kpts = [ - cast(tuple[int, int, int], tuple(int(i) for i in tokens)), + cast(Kpoint, tuple(int(i) for i in tokens)), ] elif name == "usershift": - kpoint.kpts_shift = cast(tuple[float, float, float], tuple(float(i) for i in tokens)) + kpoint.kpts_shift = cast(Vector3D, tuple(float(i) for i in tokens)) elif name in {"genvec1", "genvec2", "genvec3", "shift"}: setattr(kpoint, name, [float(i) for i in tokens]) @@ -1435,7 +1436,7 @@ def _parse_kpoints(elem: XML_Element) -> tuple[Kpoints, list[tuple[float, float, for va in elem.findall("varray"): name = va.attrib["name"] if name == "kpointlist": - actual_kpoints = cast(list[tuple[float, float, float]], list(map(tuple, _parse_vasp_array(va)))) + actual_kpoints = cast(list[Tuple3Floats], list(map(tuple, _parse_vasp_array(va)))) elif name == "weights": weights = [i[0] for i in _parse_vasp_array(va)] elem.clear() diff --git a/pymatgen/io/vasp/sets.py b/pymatgen/io/vasp/sets.py index b2e383401fe..be252b4a5d2 100644 --- a/pymatgen/io/vasp/sets.py +++ b/pymatgen/io/vasp/sets.py @@ -60,7 +60,7 @@ from typing_extensions import Self - from pymatgen.util.typing import Vector3D + from pymatgen.util.typing import Tuple3Ints, Vector3D UserPotcarFunctional = Union[ Literal["PBE", "PBE_52", "PBE_54", "LDA", "LDA_52", "LDA_54", "PW91", "LDA_US", "PW91_US"], None @@ -1791,7 +1791,7 @@ class MPSOCSet(VaspInputSet): **kwargs: Keywords supported by VaspInputSet. """ - saxis: tuple[int, int, int] = (0, 0, 1) + saxis: Tuple3Ints = (0, 0, 1) nbands_factor: float = 1.2 lepsilon: bool = False lcalcpol: bool = False @@ -2978,7 +2978,7 @@ class MPAbsorptionSet(VaspInputSet): copy_wavecar: bool = True nbands_factor: float = 2 reciprocal_density: float = 400 - nkred: tuple[int, int, int] | None = None + nkred: Tuple3Ints | None = None nedos: int = 2001 inherit_incar: bool = True force_gamma: bool = True diff --git a/pymatgen/phonon/bandstructure.py b/pymatgen/phonon/bandstructure.py index f6bf5155321..d22c535acba 100644 --- a/pymatgen/phonon/bandstructure.py +++ b/pymatgen/phonon/bandstructure.py @@ -20,8 +20,10 @@ from numpy.typing import ArrayLike from typing_extensions import Self + from pymatgen.util.typing import Tuple3Ints -def get_reasonable_repetitions(n_atoms: int) -> tuple[int, int, int]: + +def get_reasonable_repetitions(n_atoms: int) -> Tuple3Ints: """Choose the number of repetitions in a supercell according to the number of atoms in the system. """ diff --git a/pymatgen/symmetry/groups.py b/pymatgen/symmetry/groups.py index e68fe464ce1..d9fbc66cfdb 100644 --- a/pymatgen/symmetry/groups.py +++ b/pymatgen/symmetry/groups.py @@ -180,7 +180,7 @@ class SpaceGroup(SymmetryGroup): """ SYMM_OPS = loadfn(os.path.join(os.path.dirname(__file__), "symm_ops.json")) - SG_SYMBOLS: ClassVar[set] = set(SYMM_DATA["space_group_encoding"]) + SG_SYMBOLS: ClassVar[set[str]] = set(SYMM_DATA["space_group_encoding"]) for op in SYMM_OPS: op["hermann_mauguin"] = re.sub(r" ", "", op["hermann_mauguin"]) op["universal_h_m"] = re.sub(r" ", "", op["universal_h_m"]) @@ -191,8 +191,10 @@ class SpaceGroup(SymmetryGroup): # POINT_GROUP_ENC = SYMM_DATA["point_group_encoding"] sg_encoding = SYMM_DATA["space_group_encoding"] abbrev_sg_mapping = SYMM_DATA["abbreviated_spacegroup_symbols"] - translations: ClassVar = {k: Fraction(v) for k, v in SYMM_DATA["translations"].items()} - full_sg_mapping: ClassVar = {v["full_symbol"]: k for k, v in SYMM_DATA["space_group_encoding"].items()} + translations: ClassVar[dict[str, Fraction]] = {k: Fraction(v) for k, v in SYMM_DATA["translations"].items()} + full_sg_mapping: ClassVar[dict[str, str]] = { + v["full_symbol"]: k for k, v in SYMM_DATA["space_group_encoding"].items() + } def __init__(self, int_symbol: str) -> None: """Initialize a Space Group from its full or abbreviated international diff --git a/pymatgen/util/typing.py b/pymatgen/util/typing.py index 52b293d0ff5..7f498c5e93f 100644 --- a/pymatgen/util/typing.py +++ b/pymatgen/util/typing.py @@ -18,6 +18,9 @@ from pymatgen.entries.computed_entries import ComputedEntry, ComputedStructureEntry, GibbsComputedStructureEntry from pymatgen.entries.exp_entries import ExpEntry +# Commonly used composite types +Tuple3Ints = tuple[int, int, int] +Tuple3Floats = tuple[float, float, float] PathLike = Union[str, OsPathLike] PbcLike = tuple[bool, bool, bool] @@ -42,13 +45,13 @@ "GibbsComputedStructureEntry", ] -Vector3D = tuple[float, float, float] +Vector3D = Tuple3Floats Matrix3D = tuple[Vector3D, Vector3D, Vector3D] SitePropsType = Union[list[dict[Any, Sequence[Any]]], dict[Any, Sequence[Any]]] # Types specific to io.vasp -Kpoint = Union[tuple[float, float, float], tuple[int,]] +Kpoint = Union[Tuple3Floats, tuple[int,]] # Miller index -MillerIndex = tuple[int, int, int] +MillerIndex = Tuple3Ints diff --git a/pymatgen/vis/structure_vtk.py b/pymatgen/vis/structure_vtk.py index cddd1b98460..b4e75e40f40 100644 --- a/pymatgen/vis/structure_vtk.py +++ b/pymatgen/vis/structure_vtk.py @@ -881,7 +881,7 @@ def make_movie(structures, output_filename="movie.mp4", zoom=1.0, fps=20, bitrat class MultiStructuresVis(StructureVis): """Visualization for multiple structures.""" - DEFAULT_ANIMATED_MOVIE_OPTIONS: ClassVar = dict( + DEFAULT_ANIMATED_MOVIE_OPTIONS: ClassVar[dict[str, str | float]] = dict( time_between_frames=0.1, looping_type="restart", number_of_loops=1, diff --git a/tests/core/test_interface.py b/tests/core/test_interface.py index 3d24180cd1a..22cd178dd5b 100644 --- a/tests/core/test_interface.py +++ b/tests/core/test_interface.py @@ -38,11 +38,11 @@ def setUp(self): def test_init(self): assert self.Cu_GB1.rotation_angle == approx(123.74898859588858) assert self.Cu_GB1.vacuum_thickness == approx(1.5) - assert self.Cu_GB2.rotation_axis == [1, 2, 3] + assert self.Cu_GB2.rotation_axis == (1, 2, 3) assert_allclose(self.Cu_GB1.ab_shift, [0.0, 0.0]) assert_allclose(self.Cu_GB2.ab_shift, [0.2, 0.2]) - assert self.Cu_GB1.gb_plane == [1, 3, 1] - assert self.Cu_GB2.gb_plane == [1, 2, 3] + assert self.Cu_GB1.gb_plane == (1, 3, 1) + assert self.Cu_GB2.gb_plane == (1, 2, 3) assert_allclose(self.Cu_GB1.init_cell.lattice.matrix, self.Cu_conv.lattice.matrix) def test_copy(self):