diff --git a/src/py/mat3ra/made/basis.py b/src/py/mat3ra/made/basis.py index da34c13c..cf08c295 100644 --- a/src/py/mat3ra/made/basis.py +++ b/src/py/mat3ra/made/basis.py @@ -6,14 +6,14 @@ from pydantic import BaseModel from .cell import Cell -from .utils import ArrayWithIds +from .utils import ArrayWithIds, get_overlapping_coordinates class Basis(RoundNumericValuesMixin, BaseModel): elements: ArrayWithIds = ArrayWithIds(values=["Si"]) coordinates: ArrayWithIds = ArrayWithIds(values=[0, 0, 0]) units: str = AtomicCoordinateUnits.crystal - cell: Optional[Cell] = None + cell: Cell = Cell() labels: Optional[ArrayWithIds] = ArrayWithIds(values=[]) constraints: Optional[ArrayWithIds] = ArrayWithIds(values=[]) @@ -76,7 +76,42 @@ def to_crystal(self): self.coordinates.map_array_in_place(self.cell.convert_point_to_crystal) self.units = AtomicCoordinateUnits.crystal - def add_atom(self, element="Si", coordinate=[0.5, 0.5, 0.5]): + def add_atom( + self, + element="Si", + coordinate: Optional[List[float]] = None, + use_cartesian_coordinates: bool = False, + force: bool = False, + ): + """ + Add an atom to the basis. + + Before adding the atom at the specified coordinate, checks that no other atom is overlapping within a threshold. + + Args: + element (str): Element symbol of the atom to be added. + coordinate (List[float]): Coordinate of the atom to be added. + use_cartesian_coordinates (bool): Whether the coordinate is in Cartesian units (or crystal by default). + force (bool): Whether to force adding the atom even if it overlaps with another atom. + """ + if coordinate is None: + coordinate = [0, 0, 0] + if use_cartesian_coordinates and self.is_in_crystal_units: + coordinate = self.cell.convert_point_to_crystal(coordinate) + if not use_cartesian_coordinates and self.is_in_cartesian_units: + coordinate = self.cell.convert_point_to_cartesian(coordinate) + cartesian_coordinates_for_overlap_check = [ + self.cell.convert_point_to_cartesian(coord) for coord in self.coordinates.values + ] + cartesian_coordinate_for_overlap_check = self.cell.convert_point_to_cartesian(coordinate) + if get_overlapping_coordinates( + cartesian_coordinate_for_overlap_check, cartesian_coordinates_for_overlap_check, threshold=0.01 + ): + if force: + print(f"Warning: Overlapping coordinates found for {coordinate}. Adding atom anyway.") + else: + print(f"Warning: Overlapping coordinates found for {coordinate}. Not adding atom.") + return self.elements.add_item(element) self.coordinates.add_item(coordinate) diff --git a/src/py/mat3ra/made/material.py b/src/py/mat3ra/made/material.py index e7931088..2b46f31c 100644 --- a/src/py/mat3ra/made/material.py +++ b/src/py/mat3ra/made/material.py @@ -108,3 +108,8 @@ def set_new_lattice_vectors( self.basis = new_basis lattice = Lattice.from_vectors_array([lattice_vector1, lattice_vector2, lattice_vector3]) self.lattice = lattice + + def add_atom(self, element: str, coordinate: List[float]) -> None: + new_basis = self.basis.copy() + new_basis.add_atom(element, coordinate) + self.basis = new_basis diff --git a/src/py/mat3ra/made/tools/analyze.py b/src/py/mat3ra/made/tools/analyze.py index 47dcc109..7bf67b1f 100644 --- a/src/py/mat3ra/made/tools/analyze.py +++ b/src/py/mat3ra/made/tools/analyze.py @@ -1,10 +1,13 @@ from typing import Callable, List, Literal, Optional import numpy as np +from scipy.spatial import cKDTree from ..material import Material from .convert import decorator_convert_material_args_kwargs_to_atoms, to_pymatgen +from .enums import SurfaceTypes from .third_party import ASEAtoms, PymatgenIStructure, PymatgenVoronoiNN +from .utils import decorator_handle_periodic_boundary_conditions @decorator_convert_material_args_kwargs_to_atoms @@ -265,6 +268,8 @@ def get_atom_indices_with_condition_on_coordinates( def get_nearest_neighbors_atom_indices( material: Material, coordinate: Optional[List[float]] = None, + tolerance: float = 0.1, + cutoff: float = 13.0, ) -> Optional[List[int]]: """ Returns the indices of direct neighboring atoms to a specified position in the material using Voronoi tessellation. @@ -272,6 +277,10 @@ def get_nearest_neighbors_atom_indices( Args: material (Material): The material object to find neighbors in. coordinate (List[float]): The position to find neighbors for. + tolerance (float): tolerance parameter for near-neighbor finding. Faces that are smaller than tol fraction + of the largest face are not included in the tessellation. (default: 0.1). + as per: https://pymatgen.org/pymatgen.analysis.html#pymatgen.analysis.local_env.VoronoiNN + cutoff (float): The cutoff radius for identifying neighbors, in angstroms. Returns: List[int]: A list of indices of neighboring atoms, or an empty list if no neighbors are found. @@ -280,17 +289,26 @@ def get_nearest_neighbors_atom_indices( coordinate = [0, 0, 0] structure = to_pymatgen(material) voronoi_nn = PymatgenVoronoiNN( - tol=0.5, - cutoff=15.0, - allow_pathological=False, + tol=tolerance, + cutoff=cutoff, weight="solid_angle", - extra_nn_info=True, + extra_nn_info=False, compute_adj_neighbors=True, ) - structure.append("X", coordinate, validate_proximity=False) - neighbors = voronoi_nn.get_nn_info(structure, len(structure.sites) - 1) + coordinates = material.basis.coordinates + site_index = coordinates.get_element_id_by_value(coordinate) + remove_dummy_atom = False + if site_index is None: + structure.append("X", coordinate, validate_proximity=False) + site_index = len(structure.sites) - 1 + remove_dummy_atom = True + try: + neighbors = voronoi_nn.get_nn_info(structure, site_index) + except ValueError: + return None neighboring_atoms_pymatgen_ids = [n["site_index"] for n in neighbors] - structure.remove_sites([-1]) + if remove_dummy_atom: + structure.remove_sites([-1]) all_coordinates = material.basis.coordinates all_coordinates.filter_by_indices(neighboring_atoms_pymatgen_ids) @@ -324,3 +342,130 @@ def get_atomic_coordinates_extremum( coordinates = new_material.basis.coordinates.to_array_of_values_with_ids() values = [coord.value[{"x": 0, "y": 1, "z": 2}[axis]] for coord in coordinates] return getattr(np, extremum)(values) + + +def is_height_within_limits(z: float, z_extremum: float, depth: float, surface: SurfaceTypes) -> bool: + """ + Check if the height of an atom is within the specified limits. + + Args: + z (float): The z-coordinate of the atom. + z_extremum (float): The extremum z-coordinate of the surface. + depth (float): The depth from the surface to look for exposed atoms. + surface (SurfaceTypes): The surface type (top or bottom). + + Returns: + bool: True if the height is within the limits, False otherwise. + """ + return (z >= z_extremum - depth) if surface == SurfaceTypes.TOP else (z <= z_extremum + depth) + + +def is_shadowed_by_neighbors_from_surface( + z: float, neighbors_indices: List[int], surface: SurfaceTypes, coordinates: np.ndarray +) -> bool: + """ + Check if any one of the neighboring atoms shadow the atom from the surface by being closer to the specified surface. + + Args: + z (float): The z-coordinate of the atom. + neighbors_indices (List[int]): List of indices of neighboring atoms. + surface (SurfaceTypes): The surface type (top or bottom). + coordinates (np.ndarray): The coordinates of the atoms. + + Returns: + bool: True if the atom is not shadowed, False otherwise. + """ + return not any( + (coordinates[n][2] > z if surface == SurfaceTypes.TOP else coordinates[n][2] < z) for n in neighbors_indices + ) + + +@decorator_handle_periodic_boundary_conditions(cutoff=0.1) +def get_surface_atom_indices( + material: Material, surface: SurfaceTypes = SurfaceTypes.TOP, shadowing_radius: float = 2.5, depth: float = 5 +) -> List[int]: + """ + Identify exposed atoms on the top or bottom surface of the material. + + Args: + material (Material): Material object to get surface atoms from. + surface (SurfaceTypes): Specify "top" or "bottom" to detect the respective surface atoms. + shadowing_radius (float): Radius for atoms shadowing underlying from detecting as exposed. + depth (float): Depth from the surface to look for exposed atoms. + + Returns: + List[int]: List of indices of exposed surface atoms. + """ + new_material = material.clone() + new_material.to_cartesian() + coordinates = np.array(new_material.basis.coordinates.values) + ids = new_material.basis.coordinates.ids + kd_tree = cKDTree(coordinates) + + z_extremum = np.max(coordinates[:, 2]) if surface == SurfaceTypes.TOP else np.min(coordinates[:, 2]) + + exposed_atoms_indices = [] + for idx, (x, y, z) in enumerate(coordinates): + if is_height_within_limits(z, z_extremum, depth, surface): + neighbors_indices = kd_tree.query_ball_point([x, y, z], r=shadowing_radius) + if is_shadowed_by_neighbors_from_surface(z, neighbors_indices, surface, coordinates): + exposed_atoms_indices.append(ids[idx]) + + return exposed_atoms_indices + + +def get_coordination_numbers( + material: Material, + indices: Optional[List[int]] = None, + cutoff: float = 3.0, +) -> List[int]: + """ + Calculate the coordination numbers of atoms in the material. + + Args: + material (Material): Material object to calculate coordination numbers for. + indices (List[int]): List of atom indices to calculate coordination numbers for. + cutoff (float): The cutoff radius for identifying neighbors. + + Returns: + List[int]: List of coordination numbers for each atom in the material. + """ + new_material = material.clone() + new_material.to_cartesian() + if indices is not None: + new_material.basis.coordinates.filter_by_indices(indices) + coordinates = np.array(new_material.basis.coordinates.values) + kd_tree = cKDTree(coordinates) + + coordination_numbers = [] + for idx, (x, y, z) in enumerate(coordinates): + neighbors = kd_tree.query_ball_point([x, y, z], r=cutoff) + # Explicitly remove the atom itself from the list of neighbors + neighbors = [n for n in neighbors if n != idx] + coordination_numbers.append(len(neighbors)) + + return coordination_numbers + + +@decorator_handle_periodic_boundary_conditions(cutoff=0.1) +def get_undercoordinated_atom_indices( + material: Material, + indices: List[int], + cutoff: float = 3.0, + coordination_threshold: int = 3, +) -> List[int]: + """ + Identify undercoordinated atoms among the specified indices in the material. + + Args: + material (Material): Material object to identify undercoordinated atoms in. + indices (List[int]): List of atom indices to check for undercoordination. + cutoff (float): The cutoff radius for identifying neighbors. + coordination_threshold (int): The coordination number threshold for undercoordination. + + Returns: + List[int]: List of indices of undercoordinated atoms. + """ + coordination_numbers = get_coordination_numbers(material, indices, cutoff) + undercoordinated_atoms_indices = [i for i, cn in enumerate(coordination_numbers) if cn <= coordination_threshold] + return undercoordinated_atoms_indices diff --git a/src/py/mat3ra/made/tools/build/__init__.py b/src/py/mat3ra/made/tools/build/__init__.py index a6341f0a..e58f6617 100644 --- a/src/py/mat3ra/made/tools/build/__init__.py +++ b/src/py/mat3ra/made/tools/build/__init__.py @@ -117,11 +117,11 @@ def get_material( ) -> Material: return self.get_materials(configuration, selector_parameters, post_process_parameters)[0] - def _update_material_name(self, material, configuration): + def _update_material_name(self, material, configuration) -> Material: # Do nothing by default return material - def _update_material_metadata(self, material, configuration): + def _update_material_metadata(self, material, configuration) -> Material: if "build" not in material.metadata: material.metadata["build"] = {} material.metadata["build"]["configuration"] = configuration.to_json() diff --git a/src/py/mat3ra/made/tools/build/defect/builders.py b/src/py/mat3ra/made/tools/build/defect/builders.py index a9c25abd..637562bf 100644 --- a/src/py/mat3ra/made/tools/build/defect/builders.py +++ b/src/py/mat3ra/made/tools/build/defect/builders.py @@ -30,8 +30,7 @@ get_closest_site_id_from_coordinate_and_element, ) from ....utils import get_center_of_coordinates -from ...utils import transform_coordinate_to_supercell -from ...utils import coordinate as CoordinateCondition +from ...utils import transform_coordinate_to_supercell, coordinate as CoordinateCondition from ..utils import merge_materials from ..slab import SlabConfiguration, create_slab, Termination from ..supercell import create_supercell @@ -246,7 +245,7 @@ def get_equidistant_position( ) neighboring_atoms_ids_in_supercell = get_nearest_neighbors_atom_indices( - supercell_material, adatom_coordinate_in_supercell + material=supercell_material, coordinate=adatom_coordinate_in_supercell ) if neighboring_atoms_ids_in_supercell is None: raise ValueError("No neighboring atoms found. Try reducing the distance_z.") @@ -392,6 +391,10 @@ class IslandSlabDefectBuilder(SlabDefectBuilder): _ConfigurationType: type(IslandSlabDefectConfiguration) = IslandSlabDefectConfiguration # type: ignore _GeneratedItemType: Material = Material + @staticmethod + def _default_condition(coordinate: List[float]): + return True + def create_island( self, material: Material, @@ -410,28 +413,26 @@ def create_island( Returns: The material with the island added. """ - new_material = material.clone() - original_max_z = get_atomic_coordinates_extremum(new_material, use_cartesian_coordinates=False) + original_max_z = get_atomic_coordinates_extremum(new_material, use_cartesian_coordinates=True) material_with_additional_layers = self.create_material_with_additional_layers(new_material, thickness) - added_layers_max_z = get_atomic_coordinates_extremum(material_with_additional_layers) - + added_layers_max_z = get_atomic_coordinates_extremum( + material_with_additional_layers, use_cartesian_coordinates=True + ) if condition is None: - - def condition(coordinate: List[float]): - return True + condition = self._default_condition atoms_within_island = filter_by_condition_on_coordinates( material=material_with_additional_layers, condition=condition, use_cartesian_coordinates=use_cartesian_coordinates, ) - - # Filter atoms in the added layers + # Filter atoms in the added layers between the original and added layers island_material = filter_by_box( material=atoms_within_island, min_coordinate=[0, 0, original_max_z], - max_coordinate=[1, 1, added_layers_max_z], + max_coordinate=[material.lattice.a, material.lattice.b, added_layers_max_z], + use_cartesian_coordinates=True, ) return self.merge_slab_and_defect(island_material, new_material) diff --git a/src/py/mat3ra/made/tools/build/passivation/__init__.py b/src/py/mat3ra/made/tools/build/passivation/__init__.py new file mode 100644 index 00000000..32ee79e4 --- /dev/null +++ b/src/py/mat3ra/made/tools/build/passivation/__init__.py @@ -0,0 +1,18 @@ +from typing import Union + +from mat3ra.made.material import Material +from .configuration import PassivationConfiguration +from .builders import ( + SurfacePassivationBuilder, + CoordinationBasedPassivationBuilder, + SurfacePassivationBuilderParameters, +) + + +def create_passivation( + configuration: PassivationConfiguration, + builder: Union[SurfacePassivationBuilder, CoordinationBasedPassivationBuilder, None] = None, +) -> Material: + if builder is None: + builder = SurfacePassivationBuilder(build_parameters=SurfacePassivationBuilderParameters()) + return builder.get_material(configuration) diff --git a/src/py/mat3ra/made/tools/build/passivation/builders.py b/src/py/mat3ra/made/tools/build/passivation/builders.py new file mode 100644 index 00000000..dd383d04 --- /dev/null +++ b/src/py/mat3ra/made/tools/build/passivation/builders.py @@ -0,0 +1,214 @@ +from typing import List + +import numpy as np +from mat3ra.made.material import Material +from pydantic import BaseModel + +from ...enums import SurfaceTypes +from ...analyze import ( + get_surface_atom_indices, + get_undercoordinated_atom_indices, + get_nearest_neighbors_atom_indices, + get_coordination_numbers, +) +from ...modify import translate_to_z_level +from ...build import BaseBuilder +from .configuration import ( + PassivationConfiguration, +) + + +class PassivationBuilder(BaseBuilder): + """ + Base class for passivation builders. + """ + + _GeneratedItemType = Material + _ConfigurationType = PassivationConfiguration + + def _generate(self, configuration: BaseBuilder._ConfigurationType) -> List[Material]: + return [self.create_passivated_material(configuration)] + + def _update_material_name( + self, material: BaseBuilder._GeneratedItemType, configuration: BaseBuilder._ConfigurationType + ) -> BaseBuilder._GeneratedItemType: + material = super()._update_material_name(material, configuration) + material.name += f" {configuration.passivant}-passivated" + return material + + def create_passivated_material(self, configuration: BaseBuilder._ConfigurationType) -> Material: + material = translate_to_z_level(configuration.slab, "center") + return material + + def _add_passivant_atoms(self, material: Material, coordinates: list, passivant: str) -> Material: + """ + Add passivant atoms to the provided coordinates in the material. + + Args: + material (Material): The material object to add passivant atoms to. + coordinates (list): The coordinates to add passivant atoms to. + passivant (str): The chemical symbol of the passivating atom (e.g., 'H'). + + Returns: + Material: The material object with passivation atoms added. + """ + for coord in coordinates: + material.add_atom(passivant, coord) + return material + + +class SurfacePassivationBuilderParameters(BaseModel): + """ + Parameters for the SurfacePassivationBuilder, defining how atoms near the surface are + detected and passivated. + + Args: + shadowing_radius (float): Radius around each surface atom to exclude underlying atoms from passivation. + depth (float): Depth from the topmost (or bottommost) atom into the material to consider for passivation, + accounting for features like islands, adatoms, and terraces. + """ + + shadowing_radius: float = 2.5 + depth: float = 5.0 + + +class SurfacePassivationBuilder(PassivationBuilder): + """ + Builder for passivating a surface. + + Detects surface atoms looking along Z axis and passivates either the top or bottom surface or both. + """ + + build_parameters: SurfacePassivationBuilderParameters + _DefaultBuildParameters = SurfacePassivationBuilderParameters() + _ConfigurationType = PassivationConfiguration + + def create_passivated_material(self, configuration: PassivationConfiguration) -> Material: + material = super().create_passivated_material(configuration) + passivant_coordinates_values = [] + + if configuration.surface in (SurfaceTypes.BOTTOM, SurfaceTypes.BOTH): + passivant_coordinates_values.extend( + self._get_passivant_coordinates(material, SurfaceTypes.BOTTOM, configuration) + ) + if configuration.surface in (SurfaceTypes.TOP, SurfaceTypes.BOTH): + passivant_coordinates_values.extend( + self._get_passivant_coordinates(material, SurfaceTypes.TOP, configuration) + ) + + return self._add_passivant_atoms(material, passivant_coordinates_values, configuration.passivant) + + def _get_passivant_coordinates( + self, material: Material, surface: SurfaceTypes, configuration: PassivationConfiguration + ): + """ + Calculate the coordinates for placing passivants based on the specified surface type. + + Args: + material (Material): Material to passivate. + surface (SurfaceTypes): Surface type (TOP or BOTTOM). + configuration (SurfacePassivationConfiguration): Configuration for passivation. + + Returns: + list: Coordinates where passivants should be added. + """ + surface_atoms_indices = get_surface_atom_indices( + material, surface, self.build_parameters.shadowing_radius, self.build_parameters.depth + ) + surface_atoms_coordinates = [ + material.basis.coordinates.get_element_value_by_index(i) for i in surface_atoms_indices + ] + bond_vector = ( + [0, 0, configuration.bond_length] if surface == SurfaceTypes.TOP else [0, 0, -configuration.bond_length] + ) + passivant_bond_vector_crystal = material.basis.cell.convert_point_to_crystal(bond_vector) + return (np.array(surface_atoms_coordinates) + np.array(passivant_bond_vector_crystal)).tolist() + + +class CoordinationBasedPassivationBuilderParameters(SurfacePassivationBuilderParameters): + """ + Parameters for the CoordinationPassivationBuilder. + Args: + coordination_threshold (int): The coordination number threshold for atom to be considered undercoordinated. + """ + + coordination_threshold: int = 3 + + +class CoordinationBasedPassivationBuilder(PassivationBuilder): + """ + Builder for passivating material based on coordination number of each atom. + + Detects atoms with coordination number below a threshold and passivates them. + """ + + _BuildParametersType = CoordinationBasedPassivationBuilderParameters + _DefaultBuildParameters = CoordinationBasedPassivationBuilderParameters() + + def create_passivated_material(self, configuration: PassivationConfiguration) -> Material: + material = super().create_passivated_material(configuration) + surface_atoms_indices = get_surface_atom_indices( + material=material, + surface=SurfaceTypes.TOP, + shadowing_radius=self.build_parameters.shadowing_radius, + depth=self.build_parameters.depth, + ) + undercoordinated_atoms_indices = get_undercoordinated_atom_indices( + material=material, + indices=surface_atoms_indices, + cutoff=self.build_parameters.shadowing_radius, + coordination_threshold=self.build_parameters.coordination_threshold, + ) + passivant_coordinates_values = self._get_passivant_coordinates( + material, configuration, undercoordinated_atoms_indices + ) + return self._add_passivant_atoms(material, passivant_coordinates_values, configuration.passivant) + + def _get_passivant_coordinates( + self, material: Material, configuration: PassivationConfiguration, undercoordinated_atoms_indices: list + ): + """ + Calculate the coordinates for placing passivating atoms based on the specified edge type. + + Args: + material (Material): Material to passivate. + configuration (SurfacePassivationConfiguration): Configuration for passivation. + undercoordinated_atoms_indices (list): Indices of undercoordinated atoms. + """ + passivant_coordinates = [] + for idx in undercoordinated_atoms_indices: + nearest_neighbors = get_nearest_neighbors_atom_indices( + material=material, + coordinate=material.basis.coordinates.get_element_value_by_index(idx), + cutoff=self.build_parameters.shadowing_radius, + ) + if nearest_neighbors is None: + continue + average_coordinate = np.mean( + [material.basis.coordinates.get_element_value_by_index(i) for i in nearest_neighbors], axis=0 + ) + bond_vector = material.basis.coordinates.get_element_value_by_index(idx) - average_coordinate + bond_vector = bond_vector / np.linalg.norm(bond_vector) * configuration.bond_length + passivant_bond_vector_crystal = material.basis.cell.convert_point_to_crystal(bond_vector) + + passivant_coordinates.append( + material.basis.coordinates.get_element_value_by_index(idx) + passivant_bond_vector_crystal + ) + + return passivant_coordinates + + def get_unique_coordination_numbers(self, material: Material): + """ + Get unique coordination numbers for all atoms in the material for current builder parameters. + + Args: + material (Material): The material object. + + Returns: + set: The coordination numbers for all atoms in the material. + """ + + coordination_numbers = set( + get_coordination_numbers(material=material, cutoff=self.build_parameters.shadowing_radius) + ) + return coordination_numbers diff --git a/src/py/mat3ra/made/tools/build/passivation/configuration.py b/src/py/mat3ra/made/tools/build/passivation/configuration.py new file mode 100644 index 00000000..a106f3ec --- /dev/null +++ b/src/py/mat3ra/made/tools/build/passivation/configuration.py @@ -0,0 +1,30 @@ +from mat3ra.made.material import Material + +from ...enums import SurfaceTypes +from ...build import BaseConfiguration + + +class PassivationConfiguration(BaseConfiguration): + """ + Configuration for a passivation. + + Args: + slab (Material): The Material object. + passivant (str): The passivating element. + bond_length (float): The bond length. + """ + + slab: Material + passivant: str = "H" + bond_length: float = 1.0 + surface: SurfaceTypes = SurfaceTypes.BOTH + + @property + def _json(self): + return { + "type": self.get_cls_name(), + "slab": self.slab.to_json(), + "passivant": self.passivant, + "bond_length": self.bond_length, + "surface": self.surface.value, + } diff --git a/src/py/mat3ra/made/tools/enums.py b/src/py/mat3ra/made/tools/enums.py new file mode 100644 index 00000000..c525714e --- /dev/null +++ b/src/py/mat3ra/made/tools/enums.py @@ -0,0 +1,7 @@ +from enum import Enum + + +class SurfaceTypes(str, Enum): + TOP = "top" + BOTTOM = "bottom" + BOTH = "both" diff --git a/src/py/mat3ra/made/tools/utils/__init__.py b/src/py/mat3ra/made/tools/utils/__init__.py index 95c961eb..6a2b6a16 100644 --- a/src/py/mat3ra/made/tools/utils/__init__.py +++ b/src/py/mat3ra/made/tools/utils/__init__.py @@ -2,6 +2,8 @@ from typing import Callable, List, Optional import numpy as np +from mat3ra.made.material import Material +from mat3ra.made.utils import ArrayWithIds from mat3ra.utils.matrix import convert_2x2_to_3x3 from ..third_party import PymatgenStructure @@ -106,3 +108,70 @@ def transform_coordinate_to_supercell( if reverse: converted_array = (np_coordinate - np_translation_vector) * np_scaling_factor return converted_array.tolist() + + +def decorator_handle_periodic_boundary_conditions(cutoff): + """ + Decorator to handle periodic boundary conditions. + + Copies atoms near boundaries within the cutoff distance beyond the opposite side of the cell + creating the effect of periodic boundary conditions for edge atoms. + + Results of the function are filtered to remove atoms or coordinates outside the original cell. + + Args: + cutoff (float): The cutoff distance for a border slice in crystal coordinates. + + Returns: + Callable: The decorated function. + """ + + def decorator(func): + @wraps(func) + def wrapper(material, *args, **kwargs): + augmented_material, last_id = augment_material_with_periodic_images(material, cutoff) + result = func(augmented_material, *args, **kwargs) + + if isinstance(result, list): + if all(isinstance(x, int) for x in result): + result = [id for id in result if id <= last_id] + elif all(isinstance(x, list) and len(x) == 3 for x in result): + result = [coord for coord in result if all(0 <= c < 1 for c in coord)] + return result + + return wrapper + + return decorator + + +def augment_material_with_periodic_images(material: Material, cutoff: float = 0.1): + """ + Augment the material's dataset by adding atoms from periodic images within a cutoff distance from the boundaries by + copying them to the opposite side of the cell, translated by the cell vector beyond the boundary. + + Args: + material (Material): The material to augment. + cutoff (float): The cutoff value for filtering atoms near boundaries, in crystal coordinates. + + Returns: + Tuple[Material, int]: The augmented material and the original count of atoms. + """ + last_id = material.basis.coordinates.ids[-1] + coordinates = np.array(material.basis.coordinates.values) + elements = np.array(material.basis.elements.values) + augmented_material = material.clone() + new_basis = augmented_material.basis.copy() + + for axis in range(3): + for direction in [-1, 1]: + mask = (coordinates[:, axis] < cutoff) if direction == 1 else (coordinates[:, axis] > (1 - cutoff)) + filtered_coordinates = coordinates[mask] + filtered_elements = elements[mask] + translation_vector = np.zeros(3) + translation_vector[axis] = direction + translated_coordinates = filtered_coordinates + translation_vector + for coord, elem in zip(translated_coordinates, filtered_elements): + new_basis.add_atom(elem, coord) + + augmented_material.basis = new_basis + return augmented_material, last_id diff --git a/src/py/mat3ra/made/tools/utils/coordinate.py b/src/py/mat3ra/made/tools/utils/coordinate.py index 683c54b3..f4533dc2 100644 --- a/src/py/mat3ra/made/tools/utils/coordinate.py +++ b/src/py/mat3ra/made/tools/utils/coordinate.py @@ -20,9 +20,10 @@ def is_coordinate_in_cylinder( Returns: bool: True if the coordinate is inside the cylinder, False otherwise. """ - return (coordinate[0] - center_position[0]) ** 2 + (coordinate[1] - center_position[1]) ** 2 <= radius**2 and ( - min_z <= coordinate[2] <= max_z - ) + np_coordinate = np.array(coordinate) + np_center_position = np.array(center_position) + distance_squared = np.sum((np_coordinate[:2] - np_center_position[:2]) ** 2) + return distance_squared <= radius**2 and min_z <= np_coordinate[2] <= max_z def is_coordinate_in_sphere(coordinate: List[float], center_position: List[float], radius: float = 0.25) -> bool: @@ -42,9 +43,7 @@ def is_coordinate_in_sphere(coordinate: List[float], center_position: List[float return distance_squared <= radius**2 -def is_coordinate_in_box( - coordinate: List[float], min_coordinate: List[float] = [0, 0, 0], max_coordinate: List[float] = [1, 1, 1] -) -> bool: +def is_coordinate_in_box(coordinate: List[float], min_coordinate=None, max_coordinate=None) -> bool: """ Check if a coordinate is inside a box. Args: @@ -54,9 +53,14 @@ def is_coordinate_in_box( Returns: bool: True if the coordinate is inside the box, False otherwise. """ - x_min, y_min, z_min = min_coordinate - x_max, y_max, z_max = max_coordinate - return x_min <= coordinate[0] <= x_max and y_min <= coordinate[1] <= y_max and z_min <= coordinate[2] <= z_max + if max_coordinate is None: + max_coordinate = [1, 1, 1] + if min_coordinate is None: + min_coordinate = [0, 0, 0] + np_coordinate = np.array(coordinate) + np_min_coordinate = np.array(min_coordinate) + np_max_coordinate = np.array(max_coordinate) + return bool(np.all(np_min_coordinate <= np_coordinate) and np.all(np_coordinate <= np_max_coordinate)) def is_coordinate_within_layer( diff --git a/src/py/mat3ra/made/utils.py b/src/py/mat3ra/made/utils.py index 21d180a2..d0acf855 100644 --- a/src/py/mat3ra/made/utils.py +++ b/src/py/mat3ra/made/utils.py @@ -58,6 +58,25 @@ def get_center_of_coordinates(coordinates: List[List[float]]) -> List[float]: return np.mean(np.array(coordinates), axis=0).tolist() +def get_overlapping_coordinates( + coordinate: List[float], + coordinates: List[List[float]], + threshold: float = 0.01, +) -> List[List[float]]: + """ + Find coordinates that are within a certain threshold of a given coordinate. + + Args: + coordinate (List[float]): The coordinate. + coordinates (List[List[float]]): The list of coordinates. + threshold (float): The threshold for the distance, in the units of the coordinates. + + Returns: + List[List[float]]: The list of overlapping coordinates. + """ + return [c for c in coordinates if np.linalg.norm(np.array(c) - np.array(coordinate)) < threshold] + + class ValueWithId(RoundNumericValuesMixin, BaseModel): id: int = 0 value: Any = None @@ -94,9 +113,18 @@ def to_array_of_values_with_ids(self) -> List[ValueWithId]: def get_element_value_by_index(self, index: int) -> Any: return self.values[index] if index < len(self.values) else None + def get_element_id_by_value(self, value: Any) -> Optional[int]: + try: + return self.ids[self.values.index(value)] + except ValueError: + return None + def filter_by_values(self, values: Union[List[Any], Any]): - values_to_keep = set(values) if isinstance(values, list) else {values} - filtered_items = [(v, i) for v, i in zip(self.values, self.ids) if v in values_to_keep] + def make_hashable(value): + return tuple(value) if isinstance(value, list) else value + + values_to_keep = set(make_hashable(v) for v in values) if isinstance(values, list) else {make_hashable(values)} + filtered_items = [(v, i) for v, i in zip(self.values, self.ids) if make_hashable(v) in values_to_keep] if filtered_items: values_unpacked, ids_unpacked = zip(*filtered_items) self.values = list(values_unpacked) diff --git a/tests/py/unit/fixtures.py b/tests/py/unit/fixtures.py index d857147f..d6cdc874 100644 --- a/tests/py/unit/fixtures.py +++ b/tests/py/unit/fixtures.py @@ -226,6 +226,66 @@ "isUpdated": True, } + +SI_SLAB_PASSIVATED = { + "name": "Si8(001), termination Si_P4/mmm_1, Slab H-passivated", + "basis": { + "elements": [ + {"id": 0, "value": "Si"}, + {"id": 1, "value": "Si"}, + {"id": 2, "value": "Si"}, + {"id": 3, "value": "Si"}, + {"id": 4, "value": "H"}, + {"id": 5, "value": "H"}, + ], + "coordinates": [ + {"id": 0, "value": [0.5, 0.5, 0.312499993]}, + {"id": 1, "value": [0.5, 0.0, 0.437499993]}, + {"id": 2, "value": [0.0, 0.0, 0.562499993]}, + {"id": 3, "value": [0.0, 0.5, 0.687499993]}, + {"id": 4, "value": [0.5, 0.5, 0.177186054]}, + {"id": 5, "value": [0.0, 0.5, 0.822813932]}, + ], + "units": "crystal", + "cell": [[3.867, 0.0, 0.0], [-0.0, 3.867, 0.0], [0.0, 0.0, 10.937528]], + "labels": [], + }, + "lattice": { + "a": 3.867, + "b": 3.867, + "c": 10.937527692, + "alpha": 90.0, + "beta": 90.0, + "gamma": 90.0, + "units": {"length": "angstrom", "angle": "degree"}, + "type": "TRI", + "vectors": { + "a": [3.867, 0.0, 0.0], + "b": [-0.0, 3.867, 0.0], + "c": [0.0, 0.0, 10.937527692], + "alat": 1, + "units": "angstrom", + }, + }, + "isNonPeriodic": False, + "_id": "", + "metadata": { + "boundaryConditions": {"type": "pbc", "offset": 0}, + "build": { + "configuration": { + "type": "PassivationConfiguration", + "slab": SI_SLAB, + "passivant": "H", + "bond_length": 1.48, + "surface": "both", + }, + "termination": "Si_P4/mmm_1", + }, + }, + "isUpdated": True, +} + + SI_SLAB_VACUUM = copy.deepcopy(SI_SLAB) SI_SLAB_VACUUM["basis"]["coordinates"] = [ {"id": 0, "value": [0.5, 0.5, 0.386029718]}, diff --git a/tests/py/unit/test_tools_build_defect.py b/tests/py/unit/test_tools_build_defect.py index 1ee349fa..b8153cab 100644 --- a/tests/py/unit/test_tools_build_defect.py +++ b/tests/py/unit/test_tools_build_defect.py @@ -115,19 +115,20 @@ def test_create_crystal_site_adatom(): def test_create_island(): condition = CoordinateCondition.CylinderCoordinateCondition( - center_position=[0.625, 0.5], radius=0.25, min_z=0, max_z=1 + center_position=[0.5, 0.5], radius=0.15, min_z=0, max_z=1 ) island_config = IslandSlabDefectConfiguration( - crystal=SLAB_111, + crystal=SLAB_001, defect_type="island", condition=condition, - thickness=1, + number_of_added_layers=1, ) defect = create_slab_defect(configuration=island_config, builder=IslandSlabDefectBuilder()) - # Only one atom is in the island for this configuration - assert len(defect.basis.elements.values) == len(SLAB_111.basis.elements.values) + 1 + # Only 1 atoms in the island were added for this configuration with 001 slab orientation + NUMBER_OF_ATOMS_IN_ISLAND = 1 + assert len(defect.basis.elements.values) == len(SLAB_001.basis.elements.values) + NUMBER_OF_ATOMS_IN_ISLAND assert defect.basis.elements.values[-1] == "Si" diff --git a/tests/py/unit/test_tools_build_passivation.py b/tests/py/unit/test_tools_build_passivation.py new file mode 100644 index 00000000..469b73cc --- /dev/null +++ b/tests/py/unit/test_tools_build_passivation.py @@ -0,0 +1,15 @@ +from mat3ra.made.material import Material +from mat3ra.made.tools.build.passivation.builders import SurfacePassivationBuilder, SurfacePassivationBuilderParameters +from mat3ra.made.tools.build.passivation.configuration import PassivationConfiguration +from mat3ra.utils import assertion as assertion_utils + +from .fixtures import SI_SLAB, SI_SLAB_PASSIVATED + + +def test_passivate_surface(): + config = PassivationConfiguration(slab=Material(SI_SLAB), passivant="H", bond_length=1.48, surface="both") + builder = SurfacePassivationBuilder( + build_parameters=SurfacePassivationBuilderParameters(shadowing_radius=2.5, depth=2.0) + ) + passivated_material = builder.get_material(config) + assertion_utils.assert_deep_almost_equal(SI_SLAB_PASSIVATED, passivated_material.to_json())