Merge pull request #154 from Exabyte-io/feature/SOF-7413

feature/SOF 7413

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)
 

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

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,13 +268,19 @@ 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.
 
     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

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

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)

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)