From 65e3ec86a19903272e24b7036924d70403bcd5ed Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 12 Apr 2024 10:33:15 +0800 Subject: [PATCH 01/38] put `is_symmetric/polar` next to `property` --- pymatgen/core/surface.py | 86 ++++++++++++++++++++-------------------- 1 file changed, 43 insertions(+), 43 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 03af6ad463d..61719169ed0 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -62,12 +62,12 @@ class Slab(Structure): """Class to hold information for a Slab, with additional - attributes pertaining to slabs, but the init method does not - actually create a slab. Also has additional methods that returns other information - about a Slab such as the surface area, normal, and atom adsorption. + attributes pertaining to slabs, but does not actually create a slab. + Also has additional methods for a Slab such as the surface area, + normal, and adsorbate atoms. - Note that all Slabs have the surface normal oriented perpendicular to the a - and b lattice vectors. This means the lattice vectors a and b are in the + Note that all Slabs have the surface normal oriented perpendicular to the + a and b lattice vectors. This means the lattice vectors a and b are in the surface plane and the c vector is out of the surface plane (though not necessarily perpendicular to the surface). """ @@ -214,6 +214,44 @@ def surface_area(self) -> float: matrix = self.lattice.matrix return np.linalg.norm(np.cross(matrix[0], matrix[1])) + def is_symmetric(self, symprec: float = 0.1) -> bool: + """Check if Slab is symmetric, i.e., contains inversion, mirror on (hkl) plane, + or screw axis (rotation and translation) about [hkl]. + + Args: + symprec (float): Symmetry precision used for SpaceGroup analyzer. + + Returns: + bool: Whether surfaces are symmetric. + """ + spg_analyzer = SpacegroupAnalyzer(self, symprec=symprec) + symm_ops = spg_analyzer.get_point_group_operations() + + # Check for inversion symmetry. Or if sites from surface (a) can be translated + # to surface (b) along the [hkl]-axis, surfaces are symmetric. Or because the + # two surfaces of our slabs are always parallel to the (hkl) plane, + # any operation where there's an (hkl) mirror plane has surface symmetry + return ( + spg_analyzer.is_laue() + or any(op.translation_vector[2] != 0 for op in symm_ops) + or any(np.all(op.rotation_matrix[2] == np.array([0, 0, -1])) for op in symm_ops) + ) + + def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: + """Check if the Slab is polar by computing the normalized dipole per unit area. + Normalized dipole per unit area is used as it is more reliable than + using the absolute value, which varies with surface area. + + Note that the Slab must be oxidation state decorated for this to work properly. + Otherwise, the Slab will always have a dipole moment of 0. + + Args: + tol_dipole_per_unit_area (float): A tolerance above which the Slab is + considered polar. + """ + dip_per_unit_area = self.dipole / self.surface_area + return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area + @classmethod def from_dict(cls, dct: dict[str, Any]) -> Self: # type: ignore[override] """ @@ -279,44 +317,6 @@ def copy(self, site_properties: dict[str, Any] | None = None) -> Slab: # type: reorient_lattice=self.reorient_lattice, ) - def is_symmetric(self, symprec: float = 0.1) -> bool: - """Check if Slab is symmetric, i.e., contains inversion, mirror on (hkl) plane, - or screw axis (rotation and translation) about [hkl]. - - Args: - symprec (float): Symmetry precision used for SpaceGroup analyzer. - - Returns: - bool: Whether surfaces are symmetric. - """ - spg_analyzer = SpacegroupAnalyzer(self, symprec=symprec) - symm_ops = spg_analyzer.get_point_group_operations() - - # Check for inversion symmetry. Or if sites from surface (a) can be translated - # to surface (b) along the [hkl]-axis, surfaces are symmetric. Or because the - # two surfaces of our slabs are always parallel to the (hkl) plane, - # any operation where there's an (hkl) mirror plane has surface symmetry - return ( - spg_analyzer.is_laue() - or any(op.translation_vector[2] != 0 for op in symm_ops) - or any(np.all(op.rotation_matrix[2] == np.array([0, 0, -1])) for op in symm_ops) - ) - - def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: - """Check if the Slab is polar by computing the normalized dipole per unit area. - Normalized dipole per unit area is used as it is more reliable than - using the absolute value, which varies with surface area. - - Note that the Slab must be oxidation state decorated for this to work properly. - Otherwise, the Slab will always have a dipole moment of 0. - - Args: - tol_dipole_per_unit_area (float): A tolerance above which the Slab is - considered polar. - """ - dip_per_unit_area = self.dipole / self.surface_area - return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area - def get_surface_sites(self, tag: bool = False) -> dict[str, list]: """Returns the surface sites and their indices in a dictionary. Useful for analysis involving broken bonds and for finding adsorption sites. From 90ef691ced44528893bac96d065fd9ffd84a6fe5 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 12 Apr 2024 10:48:32 +0800 Subject: [PATCH 02/38] BREAKING: reverse `get_slab` shift direction --- pymatgen/core/surface.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 61719169ed0..da3890c3a8d 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -259,7 +259,7 @@ def from_dict(cls, dct: dict[str, Any]) -> Self: # type: ignore[override] dct: dict. Returns: - Creates slab from dict. + Slab: Created from dict. """ lattice = Lattice.from_dict(dct["lattice"]) sites = [PeriodicSite.from_dict(sd, lattice) for sd in dct["sites"]] @@ -329,8 +329,8 @@ def get_surface_sites(self, tag: bool = False) -> dict[str, list]: site as well. This will only work for single-element systems for now. Args: - tag (bool): Option to adds site attribute "is_surfsite" (bool) - to all sites of slab. Defaults to False + tag (bool): Add attribute "is_surf_site" (bool) + to all sites of the Slab. Defaults to False. Returns: A dictionary grouping sites on top and bottom of the slab together. @@ -638,8 +638,8 @@ def symmetrically_add_atom( specie: str | Element | Species | None = None, coords_are_cartesian: bool = False, ) -> None: - """Add a species at a specified site in a slab. Will also add an - equivalent site on the other side of the slab to maintain symmetry. + """Add a species at a selected site in a Slab. Will also add an + equivalent site on the other side to maintain symmetry. TODO (@DanielYang59): use "site" over "point" as arg name for consistency @@ -649,8 +649,6 @@ def symmetrically_add_atom( specie: Deprecated argument name in #3691. Use 'species' instead. coords_are_cartesian (bool): If the site is in Cartesian coordinates. """ - # For now just use the species of the surface atom as the element to add - # Check if deprecated argument is used if specie is not None: warnings.warn("The argument 'specie' is deprecated. Use 'species' instead.", DeprecationWarning) @@ -1075,7 +1073,8 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No intended for other generation methods. Args: - shift (float): The shift value along the lattice c direction in Angstrom. + shift (float): The shift value along the lattice c direction + in fractional coordinates. tol (float): Tolerance to determine primitive cell. energy (float): The energy to assign to the slab. @@ -1104,10 +1103,8 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No species = self.oriented_unit_cell.species_and_occu # Shift all atoms - # DEBUG(@DanielYang59): shift value in Angstrom inconsistent with frac_coordis frac_coords = self.oriented_unit_cell.frac_coords - # DEBUG(@DanielYang59): suspicious shift direction (positive for downwards shift) - frac_coords = np.array(frac_coords) + np.array([0, 0, -shift])[None, :] + frac_coords = np.array(frac_coords) + np.array([0, 0, shift])[None, :] frac_coords -= np.floor(frac_coords) # wrap frac_coords to the [0, 1) range # Scale down z-coordinate by the number of layers From 185a79c475f06ba1275eeb2af442fec3ebcc8a36 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 12 Apr 2024 10:56:04 +0800 Subject: [PATCH 03/38] clarify `shift` for single atom --- pymatgen/core/surface.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index da3890c3a8d..4c72ca1ec58 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1228,8 +1228,7 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Clustering does not work when there is only one atom if n_atoms == 1: - # TODO (@DanielYang59): why this magic number 0.5? - shift = frac_coords[0][2] + 0.5 + shift = -frac_coords[0][2] + 0.5 # shift to center return [shift - math.floor(shift)] # Compute a Cartesian z-coordinate distance matrix From f1d54332a3f1137d94ed677686b3cc20a553391d Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 12 Apr 2024 11:05:41 +0800 Subject: [PATCH 04/38] use `isclose` to determine overlap atoms --- pymatgen/core/surface.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 4c72ca1ec58..89013b13b26 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1269,7 +1269,9 @@ def gen_possible_shifts(ftol: float) -> list[float]: return sorted(shifts) - def get_z_ranges(bonds: dict[tuple[Species | Element, Species | Element], float]) -> list[tuple[float, float]]: + def get_z_ranges( + bonds: dict[tuple[Species | Element, Species | Element], float], tol: float + ) -> list[tuple[float, float]]: """Collect occupied z ranges where each z_range is a (lower_z, upper_z) tuple. This method examines all sites in the oriented unit cell (OUC) @@ -1279,7 +1281,7 @@ def get_z_ranges(bonds: dict[tuple[Species | Element, Species | Element], float] Args: bonds (dict): A {(species1, species2): max_bond_dist} dict. - tol (float): Fractional tolerance for determine overlapping positions. + tol (float): Tolerance for determine overlapping positions in Angstrom. """ # Sanitize species in dict keys bonds = {(get_el_sp(s1), get_el_sp(s2)): dist for (s1, s2), dist in bonds.items()} @@ -1302,15 +1304,13 @@ def get_z_ranges(bonds: dict[tuple[Species | Element, Species | Element], float] z_ranges.extend([(0, z_range[1]), (z_range[0] + 1, 1)]) # Neglect overlapping positions - elif z_range[0] != z_range[1]: - # TODO (@DanielYang59): use the following for equality check - # elif not isclose(z_range[0], z_range[1], abs_tol=tol): + elif not isclose(z_range[0], z_range[1], abs_tol=tol): z_ranges.append(z_range) return z_ranges # Get occupied z_ranges - z_ranges = [] if bonds is None else get_z_ranges(bonds) + z_ranges = [] if bonds is None else get_z_ranges(bonds, ftol) slabs = [] for shift in gen_possible_shifts(ftol=ftol): From 86bbc342cbdf66d8991f2cf71067ce945c21a6c8 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 12 Apr 2024 11:29:12 +0800 Subject: [PATCH 05/38] correct tolerance unit --- pymatgen/core/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 89013b13b26..50b8d126a44 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1281,7 +1281,7 @@ def get_z_ranges( Args: bonds (dict): A {(species1, species2): max_bond_dist} dict. - tol (float): Tolerance for determine overlapping positions in Angstrom. + tol (float): Fractional tolerance for determine overlapping positions. """ # Sanitize species in dict keys bonds = {(get_el_sp(s1), get_el_sp(s2)): dist for (s1, s2), dist in bonds.items()} @@ -1310,7 +1310,7 @@ def get_z_ranges( return z_ranges # Get occupied z_ranges - z_ranges = [] if bonds is None else get_z_ranges(bonds, ftol) + z_ranges = [] if bonds is None else get_z_ranges(bonds, tol=tol) slabs = [] for shift in gen_possible_shifts(ftol=ftol): From a7f7dc0e7f8b49f670904aa2541d3923925327ab Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sat, 13 Apr 2024 10:39:04 +0800 Subject: [PATCH 06/38] Revert "put `is_symmetric/polar` next to `property`" This reverts commit 65e3ec86a19903272e24b7036924d70403bcd5ed. --- pymatgen/core/surface.py | 86 ++++++++++++++++++++-------------------- 1 file changed, 43 insertions(+), 43 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 50b8d126a44..bb93efa05fb 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -62,12 +62,12 @@ class Slab(Structure): """Class to hold information for a Slab, with additional - attributes pertaining to slabs, but does not actually create a slab. - Also has additional methods for a Slab such as the surface area, - normal, and adsorbate atoms. + attributes pertaining to slabs, but the init method does not + actually create a slab. Also has additional methods that returns other information + about a Slab such as the surface area, normal, and atom adsorption. - Note that all Slabs have the surface normal oriented perpendicular to the - a and b lattice vectors. This means the lattice vectors a and b are in the + Note that all Slabs have the surface normal oriented perpendicular to the a + and b lattice vectors. This means the lattice vectors a and b are in the surface plane and the c vector is out of the surface plane (though not necessarily perpendicular to the surface). """ @@ -214,44 +214,6 @@ def surface_area(self) -> float: matrix = self.lattice.matrix return np.linalg.norm(np.cross(matrix[0], matrix[1])) - def is_symmetric(self, symprec: float = 0.1) -> bool: - """Check if Slab is symmetric, i.e., contains inversion, mirror on (hkl) plane, - or screw axis (rotation and translation) about [hkl]. - - Args: - symprec (float): Symmetry precision used for SpaceGroup analyzer. - - Returns: - bool: Whether surfaces are symmetric. - """ - spg_analyzer = SpacegroupAnalyzer(self, symprec=symprec) - symm_ops = spg_analyzer.get_point_group_operations() - - # Check for inversion symmetry. Or if sites from surface (a) can be translated - # to surface (b) along the [hkl]-axis, surfaces are symmetric. Or because the - # two surfaces of our slabs are always parallel to the (hkl) plane, - # any operation where there's an (hkl) mirror plane has surface symmetry - return ( - spg_analyzer.is_laue() - or any(op.translation_vector[2] != 0 for op in symm_ops) - or any(np.all(op.rotation_matrix[2] == np.array([0, 0, -1])) for op in symm_ops) - ) - - def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: - """Check if the Slab is polar by computing the normalized dipole per unit area. - Normalized dipole per unit area is used as it is more reliable than - using the absolute value, which varies with surface area. - - Note that the Slab must be oxidation state decorated for this to work properly. - Otherwise, the Slab will always have a dipole moment of 0. - - Args: - tol_dipole_per_unit_area (float): A tolerance above which the Slab is - considered polar. - """ - dip_per_unit_area = self.dipole / self.surface_area - return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area - @classmethod def from_dict(cls, dct: dict[str, Any]) -> Self: # type: ignore[override] """ @@ -317,6 +279,44 @@ def copy(self, site_properties: dict[str, Any] | None = None) -> Slab: # type: reorient_lattice=self.reorient_lattice, ) + def is_symmetric(self, symprec: float = 0.1) -> bool: + """Check if Slab is symmetric, i.e., contains inversion, mirror on (hkl) plane, + or screw axis (rotation and translation) about [hkl]. + + Args: + symprec (float): Symmetry precision used for SpaceGroup analyzer. + + Returns: + bool: Whether surfaces are symmetric. + """ + spg_analyzer = SpacegroupAnalyzer(self, symprec=symprec) + symm_ops = spg_analyzer.get_point_group_operations() + + # Check for inversion symmetry. Or if sites from surface (a) can be translated + # to surface (b) along the [hkl]-axis, surfaces are symmetric. Or because the + # two surfaces of our slabs are always parallel to the (hkl) plane, + # any operation where there's an (hkl) mirror plane has surface symmetry + return ( + spg_analyzer.is_laue() + or any(op.translation_vector[2] != 0 for op in symm_ops) + or any(np.all(op.rotation_matrix[2] == np.array([0, 0, -1])) for op in symm_ops) + ) + + def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: + """Check if the Slab is polar by computing the normalized dipole per unit area. + Normalized dipole per unit area is used as it is more reliable than + using the absolute value, which varies with surface area. + + Note that the Slab must be oxidation state decorated for this to work properly. + Otherwise, the Slab will always have a dipole moment of 0. + + Args: + tol_dipole_per_unit_area (float): A tolerance above which the Slab is + considered polar. + """ + dip_per_unit_area = self.dipole / self.surface_area + return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area + def get_surface_sites(self, tag: bool = False) -> dict[str, list]: """Returns the surface sites and their indices in a dictionary. Useful for analysis involving broken bonds and for finding adsorption sites. From 8c78e433fd8f646531c6e9c8c7cbc1e87b863027 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sat, 13 Apr 2024 10:39:42 +0800 Subject: [PATCH 07/38] reapply docstring changes --- pymatgen/core/surface.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index bb93efa05fb..733a859e820 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -62,12 +62,12 @@ class Slab(Structure): """Class to hold information for a Slab, with additional - attributes pertaining to slabs, but the init method does not - actually create a slab. Also has additional methods that returns other information - about a Slab such as the surface area, normal, and atom adsorption. + attributes pertaining to slabs, but does not actually create a slab. + Also has additional methods for a Slab such as the surface area, + normal, and adsorbate atoms. - Note that all Slabs have the surface normal oriented perpendicular to the a - and b lattice vectors. This means the lattice vectors a and b are in the + Note that all Slabs have the surface normal oriented perpendicular to the + a and b lattice vectors. This means the lattice vectors a and b are in the surface plane and the c vector is out of the surface plane (though not necessarily perpendicular to the surface). """ From c7420cf0288f16218a53745f6e8a4bd46bcd8818 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 22 Apr 2024 19:00:45 +0800 Subject: [PATCH 08/38] use `isclose` to check close --- pymatgen/core/surface.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 733a859e820..1418aceadb7 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -405,8 +405,6 @@ def get_symmetric_site( symmetric properties of a slab when creating adsorbed structures or symmetric reconstructions. - TODO (@DanielYang59): use "site" over "point" as arg name for consistency - Args: point (ArrayLike): Fractional coordinate of the original site. cartesian (bool): Use Cartesian coordinates. @@ -423,7 +421,7 @@ def get_symmetric_site( for op in ops: slab = self.copy() site_other = op.operate(point) - if f"{site_other[2]:.6f}" == f"{point[2]:.6f}": + if isclose(site_other[2], point[2], abs_tol=1e-6): continue # Add dummy sites to check if the overall structure is symmetric @@ -641,8 +639,6 @@ def symmetrically_add_atom( """Add a species at a selected site in a Slab. Will also add an equivalent site on the other side to maintain symmetry. - TODO (@DanielYang59): use "site" over "point" as arg name for consistency - Args: species (str | Element | Species): The species to add. point (ArrayLike): The coordinate of the target site. From b4627e9b8d407d68f6151718e0d53b3c3d02f247 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 22 Apr 2024 19:23:59 +0800 Subject: [PATCH 09/38] fix typo in tag --- pymatgen/core/surface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 1418aceadb7..e373975d447 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1456,7 +1456,7 @@ def move_to_other_side( # Calculate the moving distance as the fractional height # of the Slab inside the cell - # DEBUG(@DanielYang59): the use actually sizes for slab/vac + # DEBUG(@DanielYang59): use actual sizes for slab/vac # instead of the input arg (min_slab/vac_size) n_layers_slab: int = math.ceil(self.min_slab_size / height) n_layers_vac: int = math.ceil(self.min_vac_size / height) From 4be466cc2ed36baa41b47c473dddbefb716e037f Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Wed, 1 May 2024 09:53:29 +0800 Subject: [PATCH 10/38] remove debug tag --- pymatgen/core/surface.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 47a816d47d2..ca8fe64b1a5 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1684,8 +1684,6 @@ def generate_all_slabs( def get_d(slab: Slab) -> float: """Determine the z-spacing between the bottom two layers for a Slab. - - TODO (@DanielYang59): this should be private/internal to ReconstructionGenerator? """ # Sort all sites by z-coordinates sorted_sites = sorted(slab, key=lambda site: site.frac_coords[2]) @@ -1693,7 +1691,6 @@ def get_d(slab: Slab) -> float: distance = None for site, next_site in zip(sorted_sites, sorted_sites[1:]): if not isclose(site.frac_coords[2], next_site.frac_coords[2], abs_tol=1e-6): - # DEBUG (@DanielYang59): code will break if no distinguishable layers found distance = next_site.frac_coords[2] - site.frac_coords[2] break From 16b1999b2efeaf98216ddbf4427d60d12f846769 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 1 May 2024 01:53:58 +0000 Subject: [PATCH 11/38] pre-commit auto-fixes --- pymatgen/core/surface.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index ca8fe64b1a5..f7c7e55d589 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1683,8 +1683,7 @@ def generate_all_slabs( def get_d(slab: Slab) -> float: - """Determine the z-spacing between the bottom two layers for a Slab. - """ + """Determine the z-spacing between the bottom two layers for a Slab.""" # Sort all sites by z-coordinates sorted_sites = sorted(slab, key=lambda site: site.frac_coords[2]) From c7f0d322032a47fa4a526490d1e30730449bf7f6 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 3 May 2024 11:15:06 +0800 Subject: [PATCH 12/38] make `center_slab` very simple --- pymatgen/core/surface.py | 62 ++++++++++++---------------------------- 1 file changed, 18 insertions(+), 44 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index f7c7e55d589..4a1af38aab9 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -730,23 +730,12 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]: def center_slab(slab: Slab) -> Slab: - """Relocate the Slab to the center such that its center - (the slab region) is close to z=0.5. + """Relocate the Slab region to the center. This makes it easier to find surface sites and apply operations like doping. - There are two possible cases: - - 1. When the slab region is completely positioned between - two vacuum layers in the cell but is not centered, we simply - shift the Slab to the center along z-axis. - - 2. If the Slab completely resides outside the cell either - from the bottom or the top, we iterate through all sites that - spill over and shift all sites such that it is now - on the other side. An edge case being, either the top - of the Slab is at z = 0 or the bottom is at z = 1. + TODO (@DanielYang59): need to check if the Slab is continuous TODO (@DanielYang59): this should be a method for `Slab`? @@ -756,34 +745,19 @@ def center_slab(slab: Slab) -> Slab: Returns: Slab: The centered Slab. """ - # Get all site indices - all_indices = list(range(len(slab))) - - # Get a reasonable cutoff radius to sample neighbors - bond_dists = sorted(nn[1] for nn in slab.get_neighbors(slab[0], 10) if nn[1] > 0) - # TODO (@DanielYang59): magic number for cutoff radius (would 3 be too large?) - cutoff_radius = bond_dists[0] * 3 - - # TODO (@DanielYang59): do we need the following complex method? - # Why don't we just calculate the center of the Slab and move it to z=0.5? - # Before moving we need to ensure there is only one Slab layer though - - # If structure is case 2, shift all the sites - # to the other side until it is case 1 - for site in slab: # DEBUG (@DanielYang59): Slab position changes during loop? - # DEBUG (@DanielYang59): sites below z=0 is not considered (only check coord > c) - if any(nn[1] >= slab.lattice.c for nn in slab.get_neighbors(site, cutoff_radius)): - # TODO (@DanielYang59): the magic offset "0.05" seems unnecessary, - # as the Slab would be centered later anyway - shift = 1 - site.frac_coords[2] + 0.05 - slab.translate_sites(all_indices, [0, 0, shift]) - - # Now the slab is case 1, move it to the center - weights = [site.species.weight for site in slab] - center_of_mass = np.average(slab.frac_coords, weights=weights, axis=0) + # Calculate shift vector + center_of_mass = np.average( + slab.frac_coords, + weights=[site.species.weight for site in slab], + axis=0, + ) shift = 0.5 - center_of_mass[2] - slab.translate_sites(all_indices, [0, 0, shift]) + # Move the slab to the center + slab.translate_sites( + indices=list(range(len(slab))), # all sites + vector=[0, 0, shift], + ) return slab @@ -799,7 +773,7 @@ def get_slab_regions( TODO (@DanielYang59): this should be a method for `Slab`? - TODO (@DanielYang59): maybe project all z coordinates to 1D? + TODO (@DanielYang59): project all z coordinates to 1D? Args: slab (Slab): The Slab to analyse. @@ -816,10 +790,10 @@ def get_slab_regions( neighbors = slab.get_neighbors(site, blength) for nn in neighbors: # TODO (@DanielYang59): use z coordinate (z<0) to check - # if a Slab is contiguous is suspicious (Slab could locate + # whether a Slab being continuous is suspicious (Slab could locate # entirely below z=0) - # Find sites with z < 0 (sites noncontiguous within cell) + # Find sites with z < 0 (sites non-continuous within cell) if nn[0].frac_coords[2] < 0: frac_coords.append(nn[0].frac_coords[2]) indices.append(nn[-2]) @@ -827,7 +801,7 @@ def get_slab_regions( if nn[-2] not in all_indices: all_indices.append(nn[-2]) - # If slab is noncontiguous + # If slab is non-continuous if frac_coords: # Locate the lowest site within the upper Slab last_fcoords = [] @@ -841,7 +815,7 @@ def get_slab_regions( frac_coords, indices = [], [] for nn in neighbors: if 1 > nn[0].frac_coords[2] > 0 and nn[0].frac_coords[2] < site.frac_coords[2]: - # Sites are noncontiguous within cell + # Sites are non-continuous within cell frac_coords.append(nn[0].frac_coords[2]) indices.append(nn[-2]) if nn[-2] not in all_indices: From f45c3861895c255e61fde27f3a3bd05f506d9686 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Fri, 3 May 2024 12:05:00 +0800 Subject: [PATCH 13/38] Failed: make center_slab and get_slab_regions methods for `Slab` --- pymatgen/core/surface.py | 191 +++++++++++++++++++++---------------- tests/core/test_surface.py | 3 + 2 files changed, 111 insertions(+), 83 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 4a1af38aab9..0770f489878 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -317,6 +317,103 @@ def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: dip_per_unit_area = self.dipole / self.surface_area return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area + def center_slab(self) -> Slab: + """Relocate the Slab region to the center. + + This makes it easier to find surface sites and apply + operations like doping. + + TODO (@DanielYang59): need to check if the Slab is continuous + + Returns: + Slab: The centered Slab. + """ + # Calculate shift vector + center_of_mass = np.average( + self.frac_coords, + weights=[site.species.weight for site in self], + axis=0, + ) + shift = 0.5 - center_of_mass[2] + + # Move the slab to the center + self.translate_sites( + indices=list(range(len(self))), # all sites + vector=[0, 0, shift], + ) + + return self + + def get_slab_regions( + self, + blength: float = 3.5, + ) -> list[tuple[float, float]]: + """Find the z-ranges for the slab region. + + Useful for discerning where the slab ends and vacuum begins + if the slab is not fully within the cell. + + TODO (@DanielYang59): project all z coordinates to 1D? + + Args: + blength (float): The bond length between atoms in Angstrom. + You generally want this value to be larger than the actual + bond length in order to find atoms that are part of the slab. + """ + frac_coords: list = [] # TODO (@DanielYang59): zip site and coords? + indices: list = [] + + all_indices: list = [] + + for site in self: + neighbors = self.get_neighbors(site, blength) + for nn in neighbors: + # TODO (@DanielYang59): use z coordinate (z<0) to check + # whether a Slab being continuous is suspicious (Slab could locate + # entirely below z=0) + + # Find sites with z < 0 (sites non-continuous within cell) + if nn[0].frac_coords[2] < 0: + frac_coords.append(nn[0].frac_coords[2]) + indices.append(nn[-2]) + + if nn[-2] not in all_indices: + all_indices.append(nn[-2]) + + # If slab is non-continuous + if frac_coords: + # Locate the lowest site within the upper Slab + last_fcoords = [] + last_indices = [] + while frac_coords: + last_fcoords = copy.copy(frac_coords) + last_indices = copy.copy(indices) + + site = self[indices[frac_coords.index(min(frac_coords))]] + neighbors = self.get_neighbors(site, blength, include_index=True, include_image=True) + frac_coords, indices = [], [] + for nn in neighbors: + if 1 > nn[0].frac_coords[2] > 0 and nn[0].frac_coords[2] < site.frac_coords[2]: + # Sites are non-continuous within cell + frac_coords.append(nn[0].frac_coords[2]) + indices.append(nn[-2]) + if nn[-2] not in all_indices: + all_indices.append(nn[-2]) + + # Locate the highest site within the lower Slab + upper_fcoords: list = [] + for site in self: + if all(nn.index not in all_indices for nn in self.get_neighbors(site, blength)): + upper_fcoords.append(site.frac_coords[2]) + coords: list = copy.copy(frac_coords) if frac_coords else copy.copy(last_fcoords) + min_top = self[last_indices[coords.index(min(coords))]].frac_coords[2] + return [(0, max(upper_fcoords)), (min_top, 1)] + + # If the entire slab region is within the cell, just + # set the range as the highest and lowest site in the Slab + sorted_sites = sorted(self, key=lambda site: site.frac_coords[2]) + return [(sorted_sites[0].frac_coords[2], sorted_sites[-1].frac_coords[2])] + def get_surface_sites(self, tag: bool = False) -> dict[str, list]: """Returns the surface sites and their indices in a dictionary. Useful for analysis involving broken bonds and for finding adsorption sites. @@ -732,12 +829,7 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]: def center_slab(slab: Slab) -> Slab: """Relocate the Slab region to the center. - This makes it easier to find surface sites and apply - operations like doping. - - TODO (@DanielYang59): need to check if the Slab is continuous - - TODO (@DanielYang59): this should be a method for `Slab`? + Deprecated, please use the center_slab method of Slab directly. Args: slab (Slab): The Slab to center. @@ -745,95 +837,28 @@ def center_slab(slab: Slab) -> Slab: Returns: Slab: The centered Slab. """ - # Calculate shift vector - center_of_mass = np.average( - slab.frac_coords, - weights=[site.species.weight for site in slab], - axis=0, - ) - shift = 0.5 - center_of_mass[2] + warnings.warn("Please use the center_slab method of Slab directly.", DeprecationWarning) - # Move the slab to the center - slab.translate_sites( - indices=list(range(len(slab))), # all sites - vector=[0, 0, shift], - ) - - return slab + return slab.center_slab() -def get_slab_regions( - slab: Slab, - blength: float = 3.5, -) -> list[tuple[float, float]]: +def get_slab_regions(slab: Slab, blength: float = 3.5) -> list[tuple[float, float]]: """Find the z-ranges for the slab region. + Deprecated, please use the get_slab_regions method of Slab directly. + Useful for discerning where the slab ends and vacuum begins if the slab is not fully within the cell. - TODO (@DanielYang59): this should be a method for `Slab`? - - TODO (@DanielYang59): project all z coordinates to 1D? - Args: slab (Slab): The Slab to analyse. blength (float): The bond length between atoms in Angstrom. You generally want this value to be larger than the actual bond length in order to find atoms that are part of the slab. """ - frac_coords: list = [] # TODO (@DanielYang59): zip site and coords? - indices: list = [] - - all_indices: list = [] - - for site in slab: - neighbors = slab.get_neighbors(site, blength) - for nn in neighbors: - # TODO (@DanielYang59): use z coordinate (z<0) to check - # whether a Slab being continuous is suspicious (Slab could locate - # entirely below z=0) - - # Find sites with z < 0 (sites non-continuous within cell) - if nn[0].frac_coords[2] < 0: - frac_coords.append(nn[0].frac_coords[2]) - indices.append(nn[-2]) - - if nn[-2] not in all_indices: - all_indices.append(nn[-2]) - - # If slab is non-continuous - if frac_coords: - # Locate the lowest site within the upper Slab - last_fcoords = [] - last_indices = [] - while frac_coords: - last_fcoords = copy.copy(frac_coords) - last_indices = copy.copy(indices) - - site = slab[indices[frac_coords.index(min(frac_coords))]] - neighbors = slab.get_neighbors(site, blength, include_index=True, include_image=True) - frac_coords, indices = [], [] - for nn in neighbors: - if 1 > nn[0].frac_coords[2] > 0 and nn[0].frac_coords[2] < site.frac_coords[2]: - # Sites are non-continuous within cell - frac_coords.append(nn[0].frac_coords[2]) - indices.append(nn[-2]) - if nn[-2] not in all_indices: - all_indices.append(nn[-2]) + warnings.warn("Please use the get_slab_regions method of Slab directly.", DeprecationWarning) - # Locate the highest site within the lower Slab - upper_fcoords: list = [] - for site in slab: - if all(nn.index not in all_indices for nn in slab.get_neighbors(site, blength)): - upper_fcoords.append(site.frac_coords[2]) - coords: list = copy.copy(frac_coords) if frac_coords else copy.copy(last_fcoords) - min_top = slab[last_indices[coords.index(min(coords))]].frac_coords[2] - return [(0, max(upper_fcoords)), (min_top, 1)] - - # If the entire slab region is within the cell, just - # set the range as the highest and lowest site in the Slab - sorted_sites = sorted(slab, key=lambda site: site.frac_coords[2]) - return [(sorted_sites[0].frac_coords[2], sorted_sites[-1].frac_coords[2])] + return slab.get_slab_regions(blength=blength) class SlabGenerator: @@ -1100,6 +1125,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No props = {k: v * n_layers_slab for k, v in props.items()} # Generate Slab + # DEBUG (@DanielYang): Structure created instead of a Slab slab = Structure(new_lattice, species * n_layers_slab, all_coords, site_properties=props) # (Optionally) Post-process the Slab @@ -1117,8 +1143,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # Center the slab layer around the vacuum if self.center_slab: - c_center = np.mean([coord[2] for coord in slab.frac_coords]) - slab.translate_sites(list(range(len(slab))), [0, 0, 0.5 - c_center]) + slab = center_slab(slab) # Reduce to primitive cell if self.primitive: @@ -1604,7 +1629,7 @@ def generate_all_slabs( will be in min_slab_size/math.ceil(self._proj_height/dhkl) multiples of oriented unit cells. """ - all_slabs = [] + all_slabs: list[Slab] = [] for miller in get_symmetrically_distinct_miller_indices(structure, max_index): gen = SlabGenerator( diff --git a/tests/core/test_surface.py b/tests/core/test_surface.py index 0fc58e56284..ae7345c5349 100644 --- a/tests/core/test_surface.py +++ b/tests/core/test_surface.py @@ -580,15 +580,18 @@ def test_move_to_other_side(self): def test_bonds_broken(self): # Querying the Materials Project database for Si struct = self.get_structure("Si") + # Conventional unit cell is supplied to ensure miller indices # correspond to usual crystallographic definitions conv_bulk = SpacegroupAnalyzer(struct).get_conventional_standard_structure() slab_gen = SlabGenerator(conv_bulk, [1, 1, 1], 10, 10, center_slab=True) + # Setting a generous estimate for max_broken_bonds # so that all terminations are generated. These slabs # are ordered by ascending number of bonds broken # which is assigned to Slab.energy slabs = slab_gen.get_slabs(bonds={("Si", "Si"): 2.40}, max_broken_bonds=30) + # Looking at the two slabs generated in VESTA, we # expect 2 and 6 bonds broken so we check for this. # Number of broken bonds are floats due to primitive From 7d0871ca0e58da0d319adddd1071e6a858512d92 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sun, 12 May 2024 18:36:07 +0800 Subject: [PATCH 14/38] revert deprecation of `center_slab` function --- pymatgen/core/surface.py | 54 ++++++++++++++++------------------------ 1 file changed, 21 insertions(+), 33 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 0f7cd60520a..93f6d139245 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -317,33 +317,6 @@ def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: dip_per_unit_area = self.dipole / self.surface_area return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area - def center_slab(self) -> Slab: - """Relocate the Slab region to the center. - - This makes it easier to find surface sites and apply - operations like doping. - - TODO (@DanielYang59): need to check if the Slab is continuous - - Returns: - Slab: The centered Slab. - """ - # Calculate shift vector - center_of_mass = np.average( - self.frac_coords, - weights=[site.species.weight for site in self], - axis=0, - ) - shift = 0.5 - center_of_mass[2] - - # Move the slab to the center - self.translate_sites( - indices=list(range(len(self))), # all sites - vector=[0, 0, shift], - ) - - return self - def get_slab_regions( self, blength: float = 3.5, @@ -826,23 +799,38 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]: warnings.warn("Equivalent sites could not be found for some indices. Surface unchanged.") -def center_slab(slab: Slab) -> Slab: +def center_slab(slab: Structure) -> Structure: """Relocate the Slab region to the center. Deprecated, please use the center_slab method of Slab directly. Args: - slab (Slab): The Slab to center. + slab (Structure): The Slab to center. Returns: - Slab: The centered Slab. + Structure: The centered Slab. """ - warnings.warn("Please use the center_slab method of Slab directly.", DeprecationWarning) + # Calculate shift vector + center_of_mass = np.average( + slab.frac_coords, + weights=[site.species.weight for site in slab], + axis=0, + ) + shift = 0.5 - center_of_mass[2] + + # Move the slab to the center + slab.translate_sites( + indices=list(range(len(slab))), # all sites + vector=[0, 0, shift], + ) - return slab.center_slab() + return slab -def get_slab_regions(slab: Slab, blength: float = 3.5) -> list[tuple[float, float]]: +def get_slab_regions( + slab: Slab, + blength: float = 3.5, +) -> list[tuple[float, float]]: """Find the z-ranges for the slab region. Deprecated, please use the get_slab_regions method of Slab directly. From 50bf70a1de5e35e047230ca2c410c256186e4659 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sun, 12 May 2024 19:03:52 +0800 Subject: [PATCH 15/38] revert "overlap position check with isclose" --- pymatgen/core/surface.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 93f6d139245..7fef9b80ba9 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1260,9 +1260,7 @@ def gen_possible_shifts(ftol: float) -> list[float]: return sorted(shifts) - def get_z_ranges( - bonds: dict[tuple[Species | Element, Species | Element], float], tol: float - ) -> list[tuple[float, float]]: + def get_z_ranges(bonds: dict[tuple[Species | Element, Species | Element], float]) -> list[tuple[float, float]]: """Collect occupied z ranges where each z_range is a (lower_z, upper_z) tuple. This method examines all sites in the oriented unit cell (OUC) @@ -1272,7 +1270,6 @@ def get_z_ranges( Args: bonds (dict): A {(species1, species2): max_bond_dist} dict. - tol (float): Fractional tolerance for determine overlapping positions. """ # Sanitize species in dict keys bonds = {(get_el_sp(s1), get_el_sp(s2)): dist for (s1, s2), dist in bonds.items()} @@ -1295,13 +1292,15 @@ def get_z_ranges( z_ranges.extend([(0, z_range[1]), (z_range[0] + 1, 1)]) # Neglect overlapping positions - elif not isclose(z_range[0], z_range[1], abs_tol=tol): + # TODO (@DanielYang59): use the following for equality check + # elif not isclose(z_range[0], z_range[1], abs_tol=tol): + elif z_range[0] != z_range[1]: z_ranges.append(z_range) return z_ranges # Get occupied z_ranges - z_ranges = [] if bonds is None else get_z_ranges(bonds, tol=tol) + z_ranges = [] if bonds is None else get_z_ranges(bonds) slabs = [] for shift in gen_possible_shifts(ftol=ftol): From 8541a8345ad92396c890f561a1ab486e24819c88 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sun, 12 May 2024 19:10:46 +0800 Subject: [PATCH 16/38] revert deprecation of center_slab and get_slab_regions --- pymatgen/core/surface.py | 155 ++++++++++++++++++--------------------- 1 file changed, 72 insertions(+), 83 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 7fef9b80ba9..21d5ab27385 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -317,76 +317,6 @@ def is_polar(self, tol_dipole_per_unit_area: float = 1e-3) -> bool: dip_per_unit_area = self.dipole / self.surface_area return np.linalg.norm(dip_per_unit_area) > tol_dipole_per_unit_area - def get_slab_regions( - self, - blength: float = 3.5, - ) -> list[tuple[float, float]]: - """Find the z-ranges for the slab region. - - Useful for discerning where the slab ends and vacuum begins - if the slab is not fully within the cell. - - TODO (@DanielYang59): project all z coordinates to 1D? - - Args: - blength (float): The bond length between atoms in Angstrom. - You generally want this value to be larger than the actual - bond length in order to find atoms that are part of the slab. - """ - frac_coords: list = [] # TODO (@DanielYang59): zip site and coords? - indices: list = [] - - all_indices: list = [] - - for site in self: - neighbors = self.get_neighbors(site, blength) - for nn in neighbors: - # TODO (@DanielYang59): use z coordinate (z<0) to check - # whether a Slab being continuous is suspicious (Slab could locate - # entirely below z=0) - - # Find sites with z < 0 (sites non-continuous within cell) - if nn[0].frac_coords[2] < 0: - frac_coords.append(nn[0].frac_coords[2]) - indices.append(nn[-2]) - - if nn[-2] not in all_indices: - all_indices.append(nn[-2]) - - # If slab is non-continuous - if frac_coords: - # Locate the lowest site within the upper Slab - last_fcoords = [] - last_indices = [] - while frac_coords: - last_fcoords = copy.copy(frac_coords) - last_indices = copy.copy(indices) - - site = self[indices[frac_coords.index(min(frac_coords))]] - neighbors = self.get_neighbors(site, blength, include_index=True, include_image=True) - frac_coords, indices = [], [] - for nn in neighbors: - if 1 > nn[0].frac_coords[2] > 0 and nn[0].frac_coords[2] < site.frac_coords[2]: - # Sites are non-continuous within cell - frac_coords.append(nn[0].frac_coords[2]) - indices.append(nn[-2]) - if nn[-2] not in all_indices: - all_indices.append(nn[-2]) - - # Locate the highest site within the lower Slab - upper_fcoords: list = [] - for site in self: - if all(nn.index not in all_indices for nn in self.get_neighbors(site, blength)): - upper_fcoords.append(site.frac_coords[2]) - coords: list = copy.copy(frac_coords) if frac_coords else copy.copy(last_fcoords) - min_top = self[last_indices[coords.index(min(coords))]].frac_coords[2] - return [(0, max(upper_fcoords)), (min_top, 1)] - - # If the entire slab region is within the cell, just - # set the range as the highest and lowest site in the Slab - sorted_sites = sorted(self, key=lambda site: site.frac_coords[2]) - return [(sorted_sites[0].frac_coords[2], sorted_sites[-1].frac_coords[2])] - def get_surface_sites(self, tag: bool = False) -> dict[str, list]: """Get the surface sites and their indices in a dictionary. Useful for analysis involving broken bonds and for finding adsorption sites. @@ -800,15 +730,23 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]: def center_slab(slab: Structure) -> Structure: - """Relocate the Slab region to the center. - - Deprecated, please use the center_slab method of Slab directly. - - Args: - slab (Structure): The Slab to center. - - Returns: - Structure: The centered Slab. + """Relocate the Slab to the center such that its center + (the slab region) is close to z=0.5. + + This makes it easier to find surface sites and apply + operations like doping. + + There are two possible cases: + 1. When the slab region is completely positioned between + two vacuum layers in the cell but is not centered, we simply + shift the Slab to the center along z-axis. + 2. If the Slab completely resides outside the cell either + from the bottom or the top, we iterate through all sites that + spill over and shift all sites such that it is now + on the other side. An edge case being, either the top + of the Slab is at z = 0 or the bottom is at z = 1. + + TODO (@DanielYang59): this should be a method for `Slab`? """ # Calculate shift vector center_of_mass = np.average( @@ -833,8 +771,6 @@ def get_slab_regions( ) -> list[tuple[float, float]]: """Find the z-ranges for the slab region. - Deprecated, please use the get_slab_regions method of Slab directly. - Useful for discerning where the slab ends and vacuum begins if the slab is not fully within the cell. @@ -843,10 +779,63 @@ def get_slab_regions( blength (float): The bond length between atoms in Angstrom. You generally want this value to be larger than the actual bond length in order to find atoms that are part of the slab. + + TODO (@DanielYang59): this should be a method for `Slab`? + TODO (@DanielYang59): maybe project all z coordinates to 1D? """ - warnings.warn("Please use the get_slab_regions method of Slab directly.", DeprecationWarning) + frac_coords: list = [] # TODO (@DanielYang59): zip site and coords? + indices: list = [] + + all_indices: list = [] + + for site in slab: + neighbors = slab.get_neighbors(site, blength) + for nn in neighbors: + # TODO (@DanielYang59): use z coordinate (z<0) to check + # if a Slab is contiguous is suspicious (Slab could locate + # entirely below z=0) + + # Find sites with z < 0 (sites noncontiguous within cell) + if nn[0].frac_coords[2] < 0: + frac_coords.append(nn[0].frac_coords[2]) + indices.append(nn[-2]) + + if nn[-2] not in all_indices: + all_indices.append(nn[-2]) + + # If slab is noncontiguous + if frac_coords: + # Locate the lowest site within the upper Slab + last_fcoords = [] + last_indices = [] + while frac_coords: + last_fcoords = copy.copy(frac_coords) + last_indices = copy.copy(indices) + + site = slab[indices[frac_coords.index(min(frac_coords))]] + neighbors = slab.get_neighbors(site, blength, include_index=True, include_image=True) + frac_coords, indices = [], [] + for nn in neighbors: + if 1 > nn[0].frac_coords[2] > 0 and nn[0].frac_coords[2] < site.frac_coords[2]: + # Sites are noncontiguous within cell + frac_coords.append(nn[0].frac_coords[2]) + indices.append(nn[-2]) + if nn[-2] not in all_indices: + all_indices.append(nn[-2]) - return slab.get_slab_regions(blength=blength) + # Locate the highest site within the lower Slab + upper_fcoords: list = [] + for site in slab: + if all(nn.index not in all_indices for nn in slab.get_neighbors(site, blength)): + upper_fcoords.append(site.frac_coords[2]) + coords: list = copy.copy(frac_coords) if frac_coords else copy.copy(last_fcoords) + min_top = slab[last_indices[coords.index(min(coords))]].frac_coords[2] + return [(0, max(upper_fcoords)), (min_top, 1)] + + # If the entire slab region is within the cell, just + # set the range as the highest and lowest site in the Slab + sorted_sites = sorted(slab, key=lambda site: site.frac_coords[2]) + return [(sorted_sites[0].frac_coords[2], sorted_sites[-1].frac_coords[2])] class SlabGenerator: From 60271c3ed1b055014015320a7480a0dac12c5145 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sun, 12 May 2024 19:17:59 +0800 Subject: [PATCH 17/38] revert changes for `center_slab` function --- pymatgen/core/surface.py | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 21d5ab27385..945f36f4e44 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -748,20 +748,34 @@ def center_slab(slab: Structure) -> Structure: TODO (@DanielYang59): this should be a method for `Slab`? """ - # Calculate shift vector - center_of_mass = np.average( - slab.frac_coords, - weights=[site.species.weight for site in slab], - axis=0, - ) + # Get all site indices + all_indices = list(range(len(slab))) + + # Get a reasonable cutoff radius to sample neighbors + bond_dists = sorted(nn[1] for nn in slab.get_neighbors(slab[0], 10) if nn[1] > 0) + # TODO (@DanielYang59): magic number for cutoff radius (would 3 be too large?) + cutoff_radius = bond_dists[0] * 3 + + # TODO (@DanielYang59): do we need the following complex method? + # Why don't we just calculate the center of the Slab and move it to z=0.5? + # Before moving we need to ensure there is only one Slab layer though + + # If structure is case 2, shift all the sites + # to the other side until it is case 1 + for site in slab: # DEBUG (@DanielYang59): Slab position changes during loop? + # DEBUG (@DanielYang59): sites below z=0 is not considered (only check coord > c) + if any(nn[1] >= slab.lattice.c for nn in slab.get_neighbors(site, cutoff_radius)): + # TODO (@DanielYang59): the magic offset "0.05" seems unnecessary, + # as the Slab would be centered later anyway + shift = 1 - site.frac_coords[2] + 0.05 + slab.translate_sites(all_indices, [0, 0, shift]) + + # Now the slab is case 1, move it to the center + weights = [site.species.weight for site in slab] + center_of_mass = np.average(slab.frac_coords, weights=weights, axis=0) shift = 0.5 - center_of_mass[2] - # Move the slab to the center - slab.translate_sites( - indices=list(range(len(slab))), # all sites - vector=[0, 0, shift], - ) - + slab.translate_sites(all_indices, [0, 0, shift]) return slab From b138d180e7e3f1057cdabd48cdc610598300445d Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sun, 12 May 2024 19:24:26 +0800 Subject: [PATCH 18/38] recover docstring --- pymatgen/core/surface.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 945f36f4e44..18d8b8bf79d 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -747,6 +747,12 @@ def center_slab(slab: Structure) -> Structure: of the Slab is at z = 0 or the bottom is at z = 1. TODO (@DanielYang59): this should be a method for `Slab`? + + Args: + slab (Structure): The Slab to center. + + Returns: + Slab: The centered Slab. """ # Get all site indices all_indices = list(range(len(slab))) @@ -776,6 +782,7 @@ def center_slab(slab: Structure) -> Structure: shift = 0.5 - center_of_mass[2] slab.translate_sites(all_indices, [0, 0, shift]) + return slab @@ -1134,7 +1141,10 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # Center the slab layer around the vacuum if self.center_slab: - slab = center_slab(slab) + # TODO: (DanielYang59): use following center_slab + # slab = center_slab(slab) + c_center = np.mean([coord[2] for coord in slab.frac_coords]) + slab.translate_sites(list(range(len(slab))), [0, 0, 0.5 - c_center]) # Reduce to primitive cell if self.primitive: @@ -1263,7 +1273,9 @@ def gen_possible_shifts(ftol: float) -> list[float]: return sorted(shifts) - def get_z_ranges(bonds: dict[tuple[Species | Element, Species | Element], float]) -> list[tuple[float, float]]: + def get_z_ranges( + bonds: dict[tuple[Species | Element, Species | Element], float], + ) -> list[tuple[float, float]]: """Collect occupied z ranges where each z_range is a (lower_z, upper_z) tuple. This method examines all sites in the oriented unit cell (OUC) From 20aca0bc5d87f120c66d1f6c5060db380fc9947e Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:04:59 +0800 Subject: [PATCH 19/38] rename slab to struct (type is Structure) --- pymatgen/core/surface.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 18d8b8bf79d..dcb3c491b18 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1123,18 +1123,17 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No props = {k: v * n_layers_slab for k, v in props.items()} # Generate Slab - # DEBUG (@DanielYang): Structure created instead of a Slab - slab = Structure(new_lattice, species * n_layers_slab, all_coords, site_properties=props) + struct: Structure = Structure(new_lattice, species * n_layers_slab, all_coords, site_properties=props) # (Optionally) Post-process the Slab # Orthogonalize the structure (through LLL lattice basis reduction) if self.lll_reduce: # Sanitize Slab (LLL reduction + site sorting + map frac_coords) - lll_slab = slab.copy(sanitize=True) - slab = lll_slab + lll_slab = struct.copy(sanitize=True) + struct = lll_slab # Apply reduction on the scaling factor - mapping = lll_slab.lattice.find_mapping(slab.lattice) + mapping = lll_slab.lattice.find_mapping(struct.lattice) if mapping is None: raise RuntimeError("LLL reduction has failed") scale_factor = np.dot(mapping[2], scale_factor) @@ -1143,22 +1142,22 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No if self.center_slab: # TODO: (DanielYang59): use following center_slab # slab = center_slab(slab) - c_center = np.mean([coord[2] for coord in slab.frac_coords]) - slab.translate_sites(list(range(len(slab))), [0, 0, 0.5 - c_center]) + c_center = np.mean([coord[2] for coord in struct.frac_coords]) + struct.translate_sites(list(range(len(struct))), [0, 0, 0.5 - c_center]) # Reduce to primitive cell if self.primitive: - prim_slab = slab.get_primitive_structure(tolerance=tol) - slab = prim_slab + prim_slab = struct.get_primitive_structure(tolerance=tol) + struct = prim_slab if energy is not None: - energy *= prim_slab.volume / slab.volume + energy *= prim_slab.volume / struct.volume # Reorient the lattice to get the correctly reduced cell ouc = self.oriented_unit_cell.copy() if self.primitive: # Find a reduced OUC - slab_l = slab.lattice + slab_l = struct.lattice ouc = ouc.get_primitive_structure( constrain_latt={ "a": slab_l.a, @@ -1173,15 +1172,15 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No ouc = ouc if (slab_l.a == ouc.lattice.a and slab_l.b == ouc.lattice.b) else self.oriented_unit_cell return Slab( - slab.lattice, - slab.species_and_occu, - slab.frac_coords, + struct.lattice, + struct.species_and_occu, + struct.frac_coords, self.miller_index, ouc, shift, scale_factor, reorient_lattice=self.reorient_lattice, - site_properties=slab.site_properties, + site_properties=struct.site_properties, energy=energy, ) From 6c2803ca4a1f1eaed6bc6a399fbba58150dbfd7b Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:06:28 +0800 Subject: [PATCH 20/38] revert adding minus sign --- pymatgen/core/surface.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index dcb3c491b18..4abe2eed878 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1231,7 +1231,8 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Clustering does not work when there is only one atom if n_atoms == 1: - shift = -frac_coords[0][2] + 0.5 # shift to center + # DEBUG (DanielYang59): missing minus sign below + shift = frac_coords[0][2] + 0.5 # shift to center return [shift - math.floor(shift)] # Compute a Cartesian z-coordinate distance matrix From df7f0f9b2314e12b07bd9b9dd57631713b385e26 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:07:34 +0800 Subject: [PATCH 21/38] fix slab type --- pymatgen/core/surface.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 4abe2eed878..ea751a31d18 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -730,7 +730,7 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]: def center_slab(slab: Structure) -> Structure: - """Relocate the Slab to the center such that its center + """Relocate the slab to the center such that its center (the slab region) is close to z=0.5. This makes it easier to find surface sites and apply @@ -739,20 +739,20 @@ def center_slab(slab: Structure) -> Structure: There are two possible cases: 1. When the slab region is completely positioned between two vacuum layers in the cell but is not centered, we simply - shift the Slab to the center along z-axis. - 2. If the Slab completely resides outside the cell either + shift the slab to the center along z-axis. + 2. If the slab completely resides outside the cell either from the bottom or the top, we iterate through all sites that spill over and shift all sites such that it is now on the other side. An edge case being, either the top - of the Slab is at z = 0 or the bottom is at z = 1. + of the slab is at z = 0 or the bottom is at z = 1. TODO (@DanielYang59): this should be a method for `Slab`? Args: - slab (Structure): The Slab to center. + slab (Structure): The slab to center. Returns: - Slab: The centered Slab. + Structure: The centered slab. """ # Get all site indices all_indices = list(range(len(slab))) From 06b2410aad5553f08d5be39ca282cfe5be8ab227 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:42:46 +0800 Subject: [PATCH 22/38] keep shift def consistent (need clarify) --- pymatgen/core/surface.py | 40 +++++++++++++++++++--------------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index deb46efd0cd..fa7061fd502 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -40,7 +40,7 @@ from collections.abc import Sequence from typing import Any - from numpy.typing import ArrayLike + from numpy.typing import ArrayLike, NDArray from typing_extensions import Self from pymatgen.core.composition import Element, Species @@ -1107,7 +1107,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # Shift all atoms frac_coords = self.oriented_unit_cell.frac_coords frac_coords = np.array(frac_coords) + np.array([0, 0, shift])[None, :] - frac_coords -= np.floor(frac_coords) # wrap frac_coords to the [0, 1) range + frac_coords -= np.floor(frac_coords) # wrap to the [0, 1) range # Scale down z-coordinate by the number of layers frac_coords[:, 2] = frac_coords[:, 2] / n_layers @@ -1195,7 +1195,7 @@ def get_slabs( repair: bool = False, ) -> list[Slab]: """Generate slabs with shift values calculated from the internal - calculate_possible_shifts method. If the user decide to avoid breaking + gen_possible_shifts method. If the user decide to avoid breaking any polyhedral bond (by setting `bonds`), any shift value that do so would be filtered out. @@ -1228,7 +1228,7 @@ def gen_possible_shifts(ftol: float) -> list[float]: two atoms are on the same plane. """ frac_coords = self.oriented_unit_cell.frac_coords - n_atoms = len(frac_coords) + n_atoms: int = len(frac_coords) # Clustering does not work when there is only one atom if n_atoms == 1: @@ -1238,7 +1238,7 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Compute a Cartesian z-coordinate distance matrix # TODO (@DanielYang59): account for periodic boundary condition - dist_matrix = np.zeros((n_atoms, n_atoms)) + dist_matrix: NDArray = np.zeros((n_atoms, n_atoms)) for i, j in itertools.combinations(list(range(n_atoms)), 2): if i != j: z_dist = frac_coords[i][2] - frac_coords[j][2] @@ -1250,34 +1250,32 @@ def gen_possible_shifts(ftol: float) -> list[float]: z_matrix = linkage(squareform(dist_matrix)) clusters = fcluster(z_matrix, ftol, criterion="distance") - # Generate a cluster to z coordinate mapping - clst_loc = {c: frac_coords[i][2] for i, c in enumerate(clusters)} + # Generate cluster to z-coordinate mapping + clst_loc: dict[Any, float] = {clst: frac_coords[idx][2] for idx, clst in enumerate(clusters)} # Wrap all clusters into the unit cell ([0, 1) range) - possible_clst = [coord - math.floor(coord) for coord in sorted(clst_loc.values())] + possible_clst: list[float] = [coord - math.floor(coord) for coord in sorted(clst_loc.values())] # Calculate shifts - n_shifts = len(possible_clst) - shifts = [] - for i in range(n_shifts): - # Handle the special case for the first-last - # z coordinate (because of periodic boundary condition) - if i == n_shifts - 1: - # TODO (@DanielYang59): Why calculate the "center" of the - # two clusters, which is not actually the shift? - shift = (possible_clst[0] + 1 + possible_clst[i]) * 0.5 + n_shifts: int = len(possible_clst) + shifts: list[float] = [] + for idx in range(n_shifts): + # Handle the special case for the first-last pair of + # z coordinates (because of periodic boundary condition) + if idx == n_shifts - 1: + shift = (possible_clst[0] + 1 + possible_clst[idx]) * 0.5 else: - shift = (possible_clst[i] + possible_clst[i + 1]) * 0.5 + shift = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 - shifts.append(shift - math.floor(shift)) + shifts.append(-(shift - math.floor(shift))) return sorted(shifts) def get_z_ranges( bonds: dict[tuple[Species | Element, Species | Element], float], ) -> list[tuple[float, float]]: - """Collect occupied z ranges where each z_range is a (lower_z, upper_z) tuple. + """Collect occupied z ranges where each range is a (lower_z, upper_z) tuple. This method examines all sites in the oriented unit cell (OUC) and considers all neighboring sites within the specified bond distance @@ -1324,7 +1322,7 @@ def get_z_ranges( # position fall within the z_range occupied by a bond) bonds_broken = 0 for z_range in z_ranges: - if z_range[0] <= shift <= z_range[1]: + if z_range[0] <= -shift <= z_range[1]: bonds_broken += 1 # DEBUG(@DanielYang59): number of bonds broken passed to energy From 7f552c8dce07d6b53fe601b82bb0ac832c480092 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:46:12 +0800 Subject: [PATCH 23/38] unify usage of shift outside Slab --- pymatgen/analysis/interfaces/coherent_interfaces.py | 4 ++-- tests/core/test_surface.py | 2 +- tests/transformations/test_advanced_transformations.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index eb23a2f0ea2..055a759d65e 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -194,8 +194,8 @@ def get_interfaces( film_shift, sub_shift = self._terminations[termination] - film_slab = film_sg.get_slab(shift=film_shift) - sub_slab = sub_sg.get_slab(shift=sub_shift) + film_slab = film_sg.get_slab(shift=-film_shift) + sub_slab = sub_sg.get_slab(shift=-sub_shift) for match in self.zsl_matches: # Build film superlattice diff --git a/tests/core/test_surface.py b/tests/core/test_surface.py index ae7345c5349..013377c32ef 100644 --- a/tests/core/test_surface.py +++ b/tests/core/test_surface.py @@ -363,7 +363,7 @@ def setUp(self): def test_get_slab(self): struct = self.get_structure("LiFePO4") gen = SlabGenerator(struct, [0, 0, 1], 10, 10) - struct = gen.get_slab(0.25) + struct = gen.get_slab(-0.25) assert struct.lattice.abc[2] == approx(20.820740000000001) fcc = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3), ["Fe"], [[0, 0, 0]]) diff --git a/tests/transformations/test_advanced_transformations.py b/tests/transformations/test_advanced_transformations.py index 71ef890053d..99b849851f7 100644 --- a/tests/transformations/test_advanced_transformations.py +++ b/tests/transformations/test_advanced_transformations.py @@ -552,7 +552,7 @@ def test_apply_transformation(self): struct = self.get_structure("LiFePO4") trans = SlabTransformation([0, 0, 1], 10, 10, shift=0.25) gen = SlabGenerator(struct, [0, 0, 1], 10, 10) - slab_from_gen = gen.get_slab(0.25) + slab_from_gen = gen.get_slab(-0.25) slab_from_trans = trans.apply_transformation(struct) assert_allclose(slab_from_gen.lattice.matrix, slab_from_trans.lattice.matrix) assert_allclose(slab_from_gen.cart_coords, slab_from_trans.cart_coords) From 01a94580d3b4db26aba1bf5dce1309bc631e20dd Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:48:57 +0800 Subject: [PATCH 24/38] fix unit test --- tests/transformations/test_advanced_transformations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/transformations/test_advanced_transformations.py b/tests/transformations/test_advanced_transformations.py index 99b849851f7..5add679ca22 100644 --- a/tests/transformations/test_advanced_transformations.py +++ b/tests/transformations/test_advanced_transformations.py @@ -550,7 +550,7 @@ def test_find_codopant(self): class TestSlabTransformation(PymatgenTest): def test_apply_transformation(self): struct = self.get_structure("LiFePO4") - trans = SlabTransformation([0, 0, 1], 10, 10, shift=0.25) + trans = SlabTransformation([0, 0, 1], 10, 10, shift=-0.25) gen = SlabGenerator(struct, [0, 0, 1], 10, 10) slab_from_gen = gen.get_slab(-0.25) slab_from_trans = trans.apply_transformation(struct) From 0a33adac0eda9431a6c20525795c50b8a1303917 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 10:56:38 +0800 Subject: [PATCH 25/38] add missing minus sign for single atom shift --- pymatgen/core/surface.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index fa7061fd502..6138c957a57 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1232,9 +1232,10 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Clustering does not work when there is only one atom if n_atoms == 1: - # DEBUG (DanielYang59): missing minus sign below - shift = frac_coords[0][2] + 0.5 # shift to center - return [shift - math.floor(shift)] + shift = -frac_coords[0][2] + 0.5 # shift to center + return [ + shift, + ] # Compute a Cartesian z-coordinate distance matrix # TODO (@DanielYang59): account for periodic boundary condition @@ -1264,7 +1265,6 @@ def gen_possible_shifts(ftol: float) -> list[float]: # z coordinates (because of periodic boundary condition) if idx == n_shifts - 1: shift = (possible_clst[0] + 1 + possible_clst[idx]) * 0.5 - else: shift = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 From d75ee6117002883d46cdc59bc74aa7c3b56cb6a9 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 11:07:39 +0800 Subject: [PATCH 26/38] [need confirm] fix test for coh-interface --- pymatgen/analysis/interfaces/coherent_interfaces.py | 8 ++++---- tests/analysis/interfaces/test_coherent_interface.py | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index 055a759d65e..3b1234bc379 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -133,11 +133,11 @@ def _find_terminations(self): film_slabs = film_sg.get_slabs() sub_slabs = sub_sg.get_slabs() - film_shifts = [s.shift for s in film_slabs] - film_terminations = [label_termination(s) for s in film_slabs] + film_shifts = [slab.shift for slab in film_slabs] + film_terminations = [label_termination(slab) for slab in film_slabs] - sub_shifts = [s.shift for s in sub_slabs] - sub_terminations = [label_termination(s) for s in sub_slabs] + sub_shifts = [slab.shift for slab in sub_slabs] + sub_terminations = [label_termination(slab) for slab in sub_slabs] self._terminations = { (film_label, sub_label): (film_shift, sub_shift) diff --git a/tests/analysis/interfaces/test_coherent_interface.py b/tests/analysis/interfaces/test_coherent_interface.py index 081cc3c4b62..5c706de4c12 100644 --- a/tests/analysis/interfaces/test_coherent_interface.py +++ b/tests/analysis/interfaces/test_coherent_interface.py @@ -40,7 +40,7 @@ def test_coherent_interface_builder(self): substrate_miller=(1, 1, 1), ) - assert len(builder.terminations) == 2 + assert len(builder.terminations) == 3 # SP: this test is super fragile and the result fluctuates between 6, 30 and 42 for # no apparent reason. The author should fix this. assert len(list(builder.get_interfaces(termination=("O2_Pmmm_1", "Si_R-3m_1")))) >= 6 From a77d6bd097f4c9bcd20093859faa694335416bfa Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 11:18:24 +0800 Subject: [PATCH 27/38] add TODO tag --- pymatgen/core/surface.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 6138c957a57..69b25cd934a 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1268,6 +1268,7 @@ def gen_possible_shifts(ftol: float) -> list[float]: else: shift = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 + # TODO (@DanielYang59): clarify the need for the minus sign below shifts.append(-(shift - math.floor(shift))) return sorted(shifts) @@ -1322,6 +1323,7 @@ def get_z_ranges( # position fall within the z_range occupied by a bond) bonds_broken = 0 for z_range in z_ranges: + # TODO (@DanielYang59): clarify the need for the minus sign below if z_range[0] <= -shift <= z_range[1]: bonds_broken += 1 From 71e102f9327be6f2bddcfd6f7f3fb62ae4b8f897 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 11:25:36 +0800 Subject: [PATCH 28/38] clarify shift definition --- pymatgen/core/surface.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 69b25cd934a..df5bae2a713 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1230,12 +1230,10 @@ def gen_possible_shifts(ftol: float) -> list[float]: frac_coords = self.oriented_unit_cell.frac_coords n_atoms: int = len(frac_coords) - # Clustering does not work when there is only one atom + # Skip clusterring when there is only one atom if n_atoms == 1: shift = -frac_coords[0][2] + 0.5 # shift to center - return [ - shift, - ] + return [shift] # Compute a Cartesian z-coordinate distance matrix # TODO (@DanielYang59): account for periodic boundary condition @@ -1268,8 +1266,8 @@ def gen_possible_shifts(ftol: float) -> list[float]: else: shift = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 - # TODO (@DanielYang59): clarify the need for the minus sign below - shifts.append(-(shift - math.floor(shift))) + # Wrap shift to [0, 1) range + shifts.append(shift - math.floor(shift)) return sorted(shifts) @@ -1323,8 +1321,7 @@ def get_z_ranges( # position fall within the z_range occupied by a bond) bonds_broken = 0 for z_range in z_ranges: - # TODO (@DanielYang59): clarify the need for the minus sign below - if z_range[0] <= -shift <= z_range[1]: + if z_range[0] <= shift <= z_range[1]: bonds_broken += 1 # DEBUG(@DanielYang59): number of bonds broken passed to energy From 3d664d25d8b3f9f820018de360709bdc9150dcd6 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 11:40:58 +0800 Subject: [PATCH 29/38] revert ALL change related to the direction of shift --- pymatgen/analysis/interfaces/coherent_interfaces.py | 4 ++-- pymatgen/core/surface.py | 2 +- tests/analysis/interfaces/test_coherent_interface.py | 2 +- tests/core/test_surface.py | 2 +- tests/transformations/test_advanced_transformations.py | 4 ++-- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index 3b1234bc379..c7fb9b08736 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -194,8 +194,8 @@ def get_interfaces( film_shift, sub_shift = self._terminations[termination] - film_slab = film_sg.get_slab(shift=-film_shift) - sub_slab = sub_sg.get_slab(shift=-sub_shift) + film_slab = film_sg.get_slab(shift=film_shift) + sub_slab = sub_sg.get_slab(shift=sub_shift) for match in self.zsl_matches: # Build film superlattice diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index df5bae2a713..efd28ad85a5 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1106,7 +1106,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # Shift all atoms frac_coords = self.oriented_unit_cell.frac_coords - frac_coords = np.array(frac_coords) + np.array([0, 0, shift])[None, :] + frac_coords = np.array(frac_coords) + np.array([0, 0, -shift])[None, :] frac_coords -= np.floor(frac_coords) # wrap to the [0, 1) range # Scale down z-coordinate by the number of layers diff --git a/tests/analysis/interfaces/test_coherent_interface.py b/tests/analysis/interfaces/test_coherent_interface.py index 5c706de4c12..081cc3c4b62 100644 --- a/tests/analysis/interfaces/test_coherent_interface.py +++ b/tests/analysis/interfaces/test_coherent_interface.py @@ -40,7 +40,7 @@ def test_coherent_interface_builder(self): substrate_miller=(1, 1, 1), ) - assert len(builder.terminations) == 3 + assert len(builder.terminations) == 2 # SP: this test is super fragile and the result fluctuates between 6, 30 and 42 for # no apparent reason. The author should fix this. assert len(list(builder.get_interfaces(termination=("O2_Pmmm_1", "Si_R-3m_1")))) >= 6 diff --git a/tests/core/test_surface.py b/tests/core/test_surface.py index 013377c32ef..ae7345c5349 100644 --- a/tests/core/test_surface.py +++ b/tests/core/test_surface.py @@ -363,7 +363,7 @@ def setUp(self): def test_get_slab(self): struct = self.get_structure("LiFePO4") gen = SlabGenerator(struct, [0, 0, 1], 10, 10) - struct = gen.get_slab(-0.25) + struct = gen.get_slab(0.25) assert struct.lattice.abc[2] == approx(20.820740000000001) fcc = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3), ["Fe"], [[0, 0, 0]]) diff --git a/tests/transformations/test_advanced_transformations.py b/tests/transformations/test_advanced_transformations.py index 5add679ca22..71ef890053d 100644 --- a/tests/transformations/test_advanced_transformations.py +++ b/tests/transformations/test_advanced_transformations.py @@ -550,9 +550,9 @@ def test_find_codopant(self): class TestSlabTransformation(PymatgenTest): def test_apply_transformation(self): struct = self.get_structure("LiFePO4") - trans = SlabTransformation([0, 0, 1], 10, 10, shift=-0.25) + trans = SlabTransformation([0, 0, 1], 10, 10, shift=0.25) gen = SlabGenerator(struct, [0, 0, 1], 10, 10) - slab_from_gen = gen.get_slab(-0.25) + slab_from_gen = gen.get_slab(0.25) slab_from_trans = trans.apply_transformation(struct) assert_allclose(slab_from_gen.lattice.matrix, slab_from_trans.lattice.matrix) assert_allclose(slab_from_gen.cart_coords, slab_from_trans.cart_coords) From adc6e603c40f48e4d2aa564e3bddabf0fed2dddc Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 11:58:45 +0800 Subject: [PATCH 30/38] rename arg shift to termination --- .../interfaces/coherent_interfaces.py | 8 +-- pymatgen/core/surface.py | 69 +++++++++++-------- 2 files changed, 44 insertions(+), 33 deletions(-) diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index c7fb9b08736..d916db14f8a 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -78,8 +78,8 @@ def _find_matches(self) -> None: reorient_lattice=False, # This is necessary to not screw up the lattice ) - film_slab = film_sg.get_slab(shift=0) - sub_slab = sub_sg.get_slab(shift=0) + film_slab = film_sg.get_slab(termination=0) + sub_slab = sub_sg.get_slab(termination=0) film_vectors = film_slab.lattice.matrix substrate_vectors = sub_slab.lattice.matrix @@ -194,8 +194,8 @@ def get_interfaces( film_shift, sub_shift = self._terminations[termination] - film_slab = film_sg.get_slab(shift=film_shift) - sub_slab = sub_sg.get_slab(shift=sub_shift) + film_slab = film_sg.get_slab(termination=film_shift) + sub_slab = sub_sg.get_slab(termination=sub_shift) for match in self.zsl_matches: # Build film superlattice diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index efd28ad85a5..90af3efe2a6 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1067,23 +1067,33 @@ def calculate_scaling_factor() -> np.ndarray: _a, _b, c = self.oriented_unit_cell.lattice.matrix self._proj_height = abs(np.dot(normal, c)) - def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = None) -> Slab: - """Generate a slab based on a given shift value along the lattice c direction. + def get_slab( + self, + termination: float = 0, + tol: float = 0.1, + energy: float | None = None, + shift: float | None = None, # TODO: remove after 2025-05-13 + ) -> Slab: + """[Private method] Generate a slab based on a given termination coordinate + along the lattice c direction. - Note: - You should rarely use this (private) method directly, which is - intended for other generation methods. + You should NOT use this method directly. Args: - shift (float): The shift value along the lattice c direction - in fractional coordinates. + termination (float): The termination coordinate + along the lattice c direction in fractional coordinates. tol (float): Tolerance to determine primitive cell. energy (float): The energy to assign to the slab. + shift (float): Deprecated, confusing. The termination coordinate + along the lattice c direction in fractional coordinates. Returns: Slab: from a shifted oriented unit cell. """ - scale_factor = self.slab_scale_factor + # Check for usage of deprecated arg "shift" + if shift is not None: + termination = shift + warnings.warn("shift is deprecated in get_slab, please use termination.") # Calculate total number of layers height = self._proj_height @@ -1106,7 +1116,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # Shift all atoms frac_coords = self.oriented_unit_cell.frac_coords - frac_coords = np.array(frac_coords) + np.array([0, 0, -shift])[None, :] + frac_coords = np.array(frac_coords) + np.array([0, 0, -termination])[None, :] frac_coords -= np.floor(frac_coords) # wrap to the [0, 1) range # Scale down z-coordinate by the number of layers @@ -1128,6 +1138,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No # (Optionally) Post-process the Slab # Orthogonalize the structure (through LLL lattice basis reduction) + scale_factor = self.slab_scale_factor if self.lll_reduce: # Sanitize Slab (LLL reduction + site sorting + map frac_coords) lll_slab = struct.copy(sanitize=True) @@ -1178,7 +1189,7 @@ def get_slab(self, shift: float = 0, tol: float = 0.1, energy: float | None = No struct.frac_coords, self.miller_index, ouc, - shift, + termination, scale_factor, reorient_lattice=self.reorient_lattice, site_properties=struct.site_properties, @@ -1220,8 +1231,8 @@ def get_slabs( sorted by the number of bonds broken. """ - def gen_possible_shifts(ftol: float) -> list[float]: - """Generate possible shifts by clustering z coordinates. + def gen_possible_terminations(ftol: float) -> list[float]: + """Generate possible terminations by clustering z coordinates. Args: ftol (float): Threshold for fcluster to check if @@ -1232,8 +1243,8 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Skip clusterring when there is only one atom if n_atoms == 1: - shift = -frac_coords[0][2] + 0.5 # shift to center - return [shift] + # TODO (DanielYang59): why the magic number 0.5? + return [frac_coords[0][2] + 0.5] # Compute a Cartesian z-coordinate distance matrix # TODO (@DanielYang59): account for periodic boundary condition @@ -1255,21 +1266,21 @@ def gen_possible_shifts(ftol: float) -> list[float]: # Wrap all clusters into the unit cell ([0, 1) range) possible_clst: list[float] = [coord - math.floor(coord) for coord in sorted(clst_loc.values())] - # Calculate shifts - n_shifts: int = len(possible_clst) - shifts: list[float] = [] - for idx in range(n_shifts): + # Calculate terminations + n_terms: int = len(possible_clst) + terminations: list[float] = [] + for idx in range(n_terms): # Handle the special case for the first-last pair of # z coordinates (because of periodic boundary condition) - if idx == n_shifts - 1: - shift = (possible_clst[0] + 1 + possible_clst[idx]) * 0.5 + if idx == n_terms - 1: + termination = (possible_clst[0] + 1 + possible_clst[idx]) * 0.5 else: - shift = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 + termination = (possible_clst[idx] + possible_clst[idx + 1]) * 0.5 - # Wrap shift to [0, 1) range - shifts.append(shift - math.floor(shift)) + # Wrap termination to [0, 1) range + terminations.append(termination - math.floor(termination)) - return sorted(shifts) + return sorted(terminations) def get_z_ranges( bonds: dict[tuple[Species | Element, Species | Element], float], @@ -1316,19 +1327,19 @@ def get_z_ranges( z_ranges = [] if bonds is None else get_z_ranges(bonds) slabs = [] - for shift in gen_possible_shifts(ftol=ftol): - # Calculate total number of bonds broken (how often the shift - # position fall within the z_range occupied by a bond) + for termination in gen_possible_terminations(ftol=ftol): + # Calculate total number of bonds broken (how often the + # termination fall within the z_range occupied by a bond) bonds_broken = 0 for z_range in z_ranges: - if z_range[0] <= shift <= z_range[1]: + if z_range[0] <= termination <= z_range[1]: bonds_broken += 1 # DEBUG(@DanielYang59): number of bonds broken passed to energy # As per the docstring this is to sort final Slabs by number # of bonds broken, but this may very likely lead to errors # if the "energy" is used literally (Maybe reset energy to None?) - slab = self.get_slab(shift=shift, tol=tol, energy=bonds_broken) + slab = self.get_slab(termination=termination, tol=tol, energy=bonds_broken) if bonds_broken <= max_broken_bonds: slabs.append(slab) From 6703c951d13e4e0f7c93e9097ea3574b95843640 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 12:08:18 +0800 Subject: [PATCH 31/38] tweak comments --- pymatgen/core/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 90af3efe2a6..0950a4580a5 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1114,7 +1114,7 @@ def get_slab( species = self.oriented_unit_cell.species_and_occu - # Shift all atoms + # Shift all atoms to the termination frac_coords = self.oriented_unit_cell.frac_coords frac_coords = np.array(frac_coords) + np.array([0, 0, -termination])[None, :] frac_coords -= np.floor(frac_coords) # wrap to the [0, 1) range @@ -1206,7 +1206,7 @@ def get_slabs( repair: bool = False, ) -> list[Slab]: """Generate slabs with shift values calculated from the internal - gen_possible_shifts method. If the user decide to avoid breaking + gen_possible_terminations method. If the user decide to avoid breaking any polyhedral bond (by setting `bonds`), any shift value that do so would be filtered out. From e0fc5f3e1928c4f2922eabec635fe2abd8ffd722 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 12:20:29 +0800 Subject: [PATCH 32/38] clarify `shift` attrib in `Slab` --- pymatgen/core/surface.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 0950a4580a5..643a91fa012 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -115,8 +115,8 @@ def __init__( you should supply the conventional structure. oriented_unit_cell (Structure): The oriented_unit_cell from which this Slab is created (by scaling in the c-direction). - shift (float): The shift in the c-direction applied to get the - termination. + shift (float): The NEGATIVE of shift in the c-direction applied + to get the termination. scale_factor (np.ndarray): scale_factor Final computed scale factor that brings the parent cell to the surface cell. reorient_lattice (bool): reorients the lattice parameters such that @@ -125,7 +125,8 @@ def __init__( that are less than 0.01 Ang apart. Defaults to False. reconstruction (str): Type of reconstruction. Defaults to None if the slab is not reconstructed. - to_unit_cell (bool): Translates fractional coordinates into the unit cell. Defaults to False. + to_unit_cell (bool): Translates fractional coordinates 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 From 9cc5a0434604f9d76ad5906e43ea15774332a40d Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 12:24:01 +0800 Subject: [PATCH 33/38] clarify magic number --- pymatgen/core/surface.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 643a91fa012..4d528f2fff5 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1244,8 +1244,9 @@ def gen_possible_terminations(ftol: float) -> list[float]: # Skip clusterring when there is only one atom if n_atoms == 1: - # TODO (DanielYang59): why the magic number 0.5? - return [frac_coords[0][2] + 0.5] + # Put the atom to the center + termination = frac_coords[0][2] + 0.5 + return [termination - math.floor(termination)] # Compute a Cartesian z-coordinate distance matrix # TODO (@DanielYang59): account for periodic boundary condition From 79430a6199cda59019665daf3212dac815c38da8 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Mon, 13 May 2024 12:27:21 +0800 Subject: [PATCH 34/38] simplify doc --- pymatgen/core/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 4d528f2fff5..955a89861d1 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1081,8 +1081,8 @@ def get_slab( You should NOT use this method directly. Args: - termination (float): The termination coordinate - along the lattice c direction in fractional coordinates. + termination (float): The coordinate along the lattice c + direction in fractional coordinates. tol (float): Tolerance to determine primitive cell. energy (float): The energy to assign to the slab. shift (float): Deprecated, confusing. The termination coordinate From b6075367feae91fff8292d08c4e7659f0e390dd9 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Thu, 16 May 2024 12:24:12 +0800 Subject: [PATCH 35/38] reuse `center_slab` func --- pymatgen/core/surface.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 955a89861d1..33e3fca0fc7 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1153,10 +1153,7 @@ def get_slab( # Center the slab layer around the vacuum if self.center_slab: - # TODO: (DanielYang59): use following center_slab - # slab = center_slab(slab) - c_center = np.mean([coord[2] for coord in struct.frac_coords]) - struct.translate_sites(list(range(len(struct))), [0, 0, 0.5 - c_center]) + struct = center_slab(struct) # Reduce to primitive cell if self.primitive: @@ -1207,7 +1204,7 @@ def get_slabs( repair: bool = False, ) -> list[Slab]: """Generate slabs with shift values calculated from the internal - gen_possible_terminations method. If the user decide to avoid breaking + gen_possible_terminations func. If the user decide to avoid breaking any polyhedral bond (by setting `bonds`), any shift value that do so would be filtered out. @@ -1242,7 +1239,7 @@ def gen_possible_terminations(ftol: float) -> list[float]: frac_coords = self.oriented_unit_cell.frac_coords n_atoms: int = len(frac_coords) - # Skip clusterring when there is only one atom + # Skip clustering when there is only one atom if n_atoms == 1: # Put the atom to the center termination = frac_coords[0][2] + 0.5 From 191b959ae73f94193bd1432facbbd638f4c5c85c Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Thu, 16 May 2024 12:29:23 +0800 Subject: [PATCH 36/38] remove unneeded TODO tag --- pymatgen/core/surface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 33e3fca0fc7..204da1339b7 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -748,8 +748,6 @@ def center_slab(slab: Structure) -> Structure: on the other side. An edge case being, either the top of the slab is at z = 0 or the bottom is at z = 1. - TODO (@DanielYang59): this should be a method for `Slab`? - Args: slab (Structure): The slab to center. From ba4f83ac4330b62b45c2c8522e20af400f670332 Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sat, 18 May 2024 10:21:15 +0800 Subject: [PATCH 37/38] revert shift to termination rename --- .../interfaces/coherent_interfaces.py | 8 +++---- pymatgen/core/surface.py | 24 +++++++------------ 2 files changed, 12 insertions(+), 20 deletions(-) diff --git a/pymatgen/analysis/interfaces/coherent_interfaces.py b/pymatgen/analysis/interfaces/coherent_interfaces.py index d916db14f8a..c7fb9b08736 100644 --- a/pymatgen/analysis/interfaces/coherent_interfaces.py +++ b/pymatgen/analysis/interfaces/coherent_interfaces.py @@ -78,8 +78,8 @@ def _find_matches(self) -> None: reorient_lattice=False, # This is necessary to not screw up the lattice ) - film_slab = film_sg.get_slab(termination=0) - sub_slab = sub_sg.get_slab(termination=0) + film_slab = film_sg.get_slab(shift=0) + sub_slab = sub_sg.get_slab(shift=0) film_vectors = film_slab.lattice.matrix substrate_vectors = sub_slab.lattice.matrix @@ -194,8 +194,8 @@ def get_interfaces( film_shift, sub_shift = self._terminations[termination] - film_slab = film_sg.get_slab(termination=film_shift) - sub_slab = sub_sg.get_slab(termination=sub_shift) + film_slab = film_sg.get_slab(shift=film_shift) + sub_slab = sub_sg.get_slab(shift=sub_shift) for match in self.zsl_matches: # Build film superlattice diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index 204da1339b7..d6b900c4063 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1068,32 +1068,24 @@ def calculate_scaling_factor() -> np.ndarray: def get_slab( self, - termination: float = 0, + shift: float = 0, tol: float = 0.1, energy: float | None = None, - shift: float | None = None, # TODO: remove after 2025-05-13 ) -> Slab: - """[Private method] Generate a slab based on a given termination coordinate - along the lattice c direction. + """[Private method] Generate a slab based on a given termination + coordinate along the lattice c direction. - You should NOT use this method directly. + You should RARELY use this method directly. Args: - termination (float): The coordinate along the lattice c + shift (float): The termination coordinate along the lattice c direction in fractional coordinates. tol (float): Tolerance to determine primitive cell. energy (float): The energy to assign to the slab. - shift (float): Deprecated, confusing. The termination coordinate - along the lattice c direction in fractional coordinates. Returns: Slab: from a shifted oriented unit cell. """ - # Check for usage of deprecated arg "shift" - if shift is not None: - termination = shift - warnings.warn("shift is deprecated in get_slab, please use termination.") - # Calculate total number of layers height = self._proj_height height_per_layer = round(height / self.parent.lattice.d_hkl(self.miller_index), 8) @@ -1115,7 +1107,7 @@ def get_slab( # Shift all atoms to the termination frac_coords = self.oriented_unit_cell.frac_coords - frac_coords = np.array(frac_coords) + np.array([0, 0, -termination])[None, :] + frac_coords = np.array(frac_coords) + np.array([0, 0, -shift])[None, :] frac_coords -= np.floor(frac_coords) # wrap to the [0, 1) range # Scale down z-coordinate by the number of layers @@ -1185,7 +1177,7 @@ def get_slab( struct.frac_coords, self.miller_index, ouc, - termination, + shift, scale_factor, reorient_lattice=self.reorient_lattice, site_properties=struct.site_properties, @@ -1336,7 +1328,7 @@ def get_z_ranges( # As per the docstring this is to sort final Slabs by number # of bonds broken, but this may very likely lead to errors # if the "energy" is used literally (Maybe reset energy to None?) - slab = self.get_slab(termination=termination, tol=tol, energy=bonds_broken) + slab = self.get_slab(shift=termination, tol=tol, energy=bonds_broken) if bonds_broken <= max_broken_bonds: slabs.append(slab) From c4f871bbf598d2a12a68cc6cce67bc245b7a0bbc Mon Sep 17 00:00:00 2001 From: "Haoyu (Daniel)" Date: Sat, 18 May 2024 17:56:11 +0800 Subject: [PATCH 38/38] remove finished TODO tag --- pymatgen/core/surface.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py index d6b900c4063..bf36eef56b9 100644 --- a/pymatgen/core/surface.py +++ b/pymatgen/core/surface.py @@ -1305,8 +1305,6 @@ def get_z_ranges( z_ranges.extend([(0, z_range[1]), (z_range[0] + 1, 1)]) # Neglect overlapping positions - # TODO (@DanielYang59): use the following for equality check - # elif not isclose(z_range[0], z_range[1], abs_tol=tol): elif z_range[0] != z_range[1]: z_ranges.append(z_range)