Merge pull request #143 from Exabyte-io/feature/SOF-7389

Feature/SOF-7389 feat: adatom defect

src/py/mat3ra/made/tools/analyze.py

@@ -1,10 +1,10 @@
-from typing import Callable, List, Optional, Literal
+from typing import Callable, List, Literal, Optional
 
 import numpy as np
 
 from ..material import Material
 from .convert import decorator_convert_material_args_kwargs_to_atoms, to_pymatgen
-from .third_party import ASEAtoms, PymatgenIStructure
+from .third_party import ASEAtoms, PymatgenIStructure, PymatgenVoronoiNN
 
 
 @decorator_convert_material_args_kwargs_to_atoms
@@ -231,6 +231,41 @@ def get_atom_indices_with_condition_on_coordinates(
     return selected_indices
 
 
+def get_nearest_neighbors_atom_indices(
+    material: Material,
+    position: Optional[List[float]] = None,
+) -> 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.
+        position (List[float]): The position to find neighbors for.
+
+    Returns:
+        List[int]: A list of indices of neighboring atoms, or an empty list if no neighbors are found.
+    """
+    if position is None:
+        position = [0, 0, 0]
+    structure = to_pymatgen(material)
+    voronoi_nn = PymatgenVoronoiNN(
+        tol=0.5,
+        cutoff=13.0,
+        allow_pathological=False,
+        weight="solid_angle",
+        extra_nn_info=True,
+        compute_adj_neighbors=True,
+    )
+    structure.append("X", position, validate_proximity=False)
+    neighbors = voronoi_nn.get_nn_info(structure, len(structure.sites) - 1)
+    neighboring_atoms_pymatgen_ids = [n["site_index"] for n in neighbors]
+    structure.remove_sites([-1])
+
+    all_coordinates = material.basis.coordinates
+    all_coordinates.filter_by_indices(neighboring_atoms_pymatgen_ids)
+    return all_coordinates.ids
+
+
 def get_atomic_coordinates_extremum(
     material: Material,
     extremum: Literal["max", "min"] = "max",

src/py/mat3ra/made/tools/build/__init__.py

@@ -65,6 +65,8 @@ def _select(
     def _post_process(
         self, items: List[_GeneratedItemType], post_process_parameters: Optional[_PostProcessParametersType]
     ) -> List[Material]:
+        if self._GeneratedItemType == Material:
+            return items
         return [Material(self._convert_generated_item(item)) for item in items]
 
     @staticmethod

src/py/mat3ra/made/tools/build/defect/__init__.py

@@ -1,11 +1,16 @@
-from typing import Optional
+from typing import Optional, Union
 
 from mat3ra.utils.factory import BaseFactory
 from mat3ra.made.material import Material
 
-from .builders import PointDefectBuilderParameters
-from .configuration import PointDefectConfiguration
-from .enums import PointDefectTypeEnum
+from .builders import (
+    PointDefectBuilderParameters,
+    SlabDefectBuilderParameters,
+    AdatomSlabDefectBuilder,
+    EquidistantAdatomSlabDefectBuilder,
+)
+from .configuration import PointDefectConfiguration, AdatomSlabDefectConfiguration
+from .enums import PointDefectTypeEnum, SlabDefectTypeEnum
 
 
 class DefectBuilderFactory(BaseFactory):
@@ -17,8 +22,8 @@ class DefectBuilderFactory(BaseFactory):
 
 
 def create_defect(
-    configuration: PointDefectConfiguration,
-    builder_parameters: Optional[PointDefectBuilderParameters] = None,
+    configuration: Union[PointDefectConfiguration, AdatomSlabDefectConfiguration],
+    builder_parameters: Union[PointDefectBuilderParameters, SlabDefectBuilderParameters, None] = None,
 ) -> Material:
     """
     Return a material with a selected defect added.
@@ -34,3 +39,22 @@ def create_defect(
     builder = BuilderClass(builder_parameters)
 
     return builder.get_material(configuration) if builder else configuration.crystal
+
+
+def create_slab_defect(
+    configuration: Union[AdatomSlabDefectConfiguration],
+    builder: Optional[Union[AdatomSlabDefectBuilder, EquidistantAdatomSlabDefectBuilder]] = None,
+) -> Material:
+    """
+    Return a material with a selected slab defect added.
+
+    Args:
+        configuration: The configuration of the defect to be added.
+        builder: The builder to be used to create the defect.
+
+    Returns:
+        The material with the defect added.
+    """
+    if builder is None:
+        builder = AdatomSlabDefectBuilder(build_parameters=SlabDefectBuilderParameters())
+    return builder.get_material(configuration)

src/py/mat3ra/made/tools/build/defect/builders.py

@@ -1,7 +1,10 @@
-from typing import List, Callable
+from typing import List, Callable, Optional
 
-from mat3ra.made.material import Material
+from mat3ra.made.tools.build.supercell import create_supercell
+from mat3ra.made.tools.modify import add_vacuum
 from pydantic import BaseModel
+from mat3ra.made.material import Material
+
 
 from ...third_party import (
     PymatgenStructure,
@@ -12,8 +15,10 @@
 )
 from ...build import BaseBuilder
 from ...convert import to_pymatgen
+from ...analyze import get_nearest_neighbors_atom_indices, get_atomic_coordinates_extremum
+from ....utils import get_center_of_coordinates
 from ..mixins import ConvertGeneratedItemsPymatgenStructureMixin
-from .configuration import PointDefectConfiguration
+from .configuration import PointDefectConfiguration, AdatomSlabDefectConfiguration
 
 
 class PointDefectBuilderParameters(BaseModel):
@@ -67,3 +72,124 @@ class SubstitutionPointDefectBuilder(PointDefectBuilder):
 
 class InterstitialPointDefectBuilder(PointDefectBuilder):
     _generator: PymatgenInterstitial = PymatgenInterstitial
+
+
+class SlabDefectBuilderParameters(BaseModel):
+    auto_add_vacuum: bool = True
+    vacuum_thickness: float = 5.0
+
+
+class SlabDefectBuilder(BaseBuilder):
+    _BuildParametersType = SlabDefectBuilderParameters
+    _DefaultBuildParameters = SlabDefectBuilderParameters()
+
+
+class AdatomSlabDefectBuilder(SlabDefectBuilder):
+    _ConfigurationType: type(AdatomSlabDefectConfiguration) = AdatomSlabDefectConfiguration  # type: ignore
+    _GeneratedItemType: Material = Material
+
+    def create_adatom(
+        self,
+        material: Material,
+        chemical_element: str = "Si",
+        position_on_surface: Optional[List[float]] = None,
+        distance_z: float = 2.0,
+    ) -> List[Material]:
+        """
+        Create an adatom at the specified position on the surface of the material.
+
+        Args:
+            material: The material to add the adatom to.
+            chemical_element: The chemical element of the adatom.
+            position_on_surface: The position on the surface of the material.
+            distance_z: The distance of the adatom from the surface.
+
+        Returns:
+            The material with the adatom added.
+        """
+        if position_on_surface is None:
+            position_on_surface = [0.5, 0.5]
+        new_material = material.clone()
+        new_basis = new_material.basis
+        adatom_position = self._calculate_position_from_2d(material, position_on_surface, distance_z)
+        new_basis.add_atom(chemical_element, adatom_position)
+        new_material.basis = new_basis
+        return [new_material]
+
+    def _calculate_position_from_2d(
+        self, material: Material, position_on_surface: List[float], distance_z: float
+    ) -> List[float]:
+        max_z = get_atomic_coordinates_extremum(material)
+        distance_z = distance_z
+        distance_in_crystal_units = distance_z / material.lattice.c
+        position = position_on_surface.copy()
+        position = position[:2]
+        position.append(max_z + distance_in_crystal_units)
+        return position
+
+    def _generate(self, configuration: _ConfigurationType) -> List[_GeneratedItemType]:
+        return self.create_adatom(
+            material=configuration.crystal,
+            chemical_element=configuration.chemical_element,
+            position_on_surface=configuration.position_on_surface,
+            distance_z=configuration.distance_z,
+        )
+
+
+class EquidistantAdatomSlabDefectBuilder(AdatomSlabDefectBuilder):
+    def create_adatom(
+        self,
+        material: Material,
+        chemical_element: str = "Si",
+        position_on_surface: Optional[List[float]] = None,
+        distance_z: float = 2.0,
+    ) -> List[Material]:
+        """
+        Create an adatom with an equidistant XY position among the nearest neighbors
+        at the given distance from the surface.
+
+        Args:
+            material: The material to add the adatom to.
+            chemical_element: The chemical element of the adatom.
+            position_on_surface: The position on the surface of the material.
+            distance_z: The distance of the adatom from the surface.
+
+        Returns:
+            The material with the adatom added.
+        """
+        if position_on_surface is None:
+            position_on_surface = [0.5, 0.5]
+        equidistant_position = self.get_equidistant_position(material, position_on_surface, distance_z)
+        new_material = material.clone()
+        if equidistant_position[2] > 1:
+            if self.build_parameters.auto_add_vacuum:
+                new_material = add_vacuum(material, self.build_parameters.vacuum_thickness)
+                equidistant_position = self.get_equidistant_position(new_material, position_on_surface, distance_z)
+            else:
+                raise ValueError("Not enough vacuum space to place the adatom.")
+
+        return super().create_adatom(new_material, chemical_element, equidistant_position, distance_z)
+
+    def get_equidistant_position(
+        self, material: Material, position_on_surface: List[float], distance_z: float = 2.0
+    ) -> List[float]:
+        new_basis = material.basis
+        adatom_position = self._calculate_position_from_2d(material, position_on_surface, distance_z)
+        neighboring_atoms_ids = get_nearest_neighbors_atom_indices(material, adatom_position)
+        # We need to check if neighboring atoms number is the same in pbc
+        supercell_material = create_supercell(material, [[3, 0, 0], [0, 3, 0], [0, 0, 1]])
+        # Move the coordinate to the central unit cell of the supercell (crystal coordinates)
+        supercell_adatom_position = [1 / 3 + adatom_position[0] / 3, 1 / 3 + adatom_position[1] / 3, adatom_position[2]]
+        supercell_neighboring_atoms_ids = get_nearest_neighbors_atom_indices(
+            supercell_material, supercell_adatom_position
+        )
+        if neighboring_atoms_ids is None or supercell_neighboring_atoms_ids is None:
+            raise ValueError("No neighboring atoms found. Try reducing the distance_z.")
+        if len(supercell_neighboring_atoms_ids) != len(neighboring_atoms_ids):
+            raise ValueError("Number of neighboring atoms is not the same in PBC. Try increasing the supercell size.")
+        neighboring_atoms_coordinates = [new_basis.coordinates.values[atom_id] for atom_id in neighboring_atoms_ids]
+
+        equidistant_position = get_center_of_coordinates(neighboring_atoms_coordinates)
+        equidistant_position[2] = adatom_position[2]
+
+        return equidistant_position

src/py/mat3ra/made/tools/build/defect/configuration.py

@@ -5,7 +5,7 @@
 from mat3ra.made.material import Material
 
 from ...analyze import get_closest_site_id_from_position
-from .enums import PointDefectTypeEnum
+from .enums import PointDefectTypeEnum, SlabDefectTypeEnum
 
 
 class BaseDefectConfiguration(BaseModel):
@@ -23,7 +23,7 @@ def from_site_id(
         cls, crystal: Material, defect_type: PointDefectTypeEnum, site_id: int, chemical_element: Optional[str] = None
     ):
         if not crystal:
-            RuntimeError("Crystal is not defined")
+            raise RuntimeError("Crystal is not defined")
         position = crystal.coordinates_array[site_id]
         return cls(crystal=crystal, defect_type=defect_type, position=position, chemical_element=chemical_element)
 
@@ -36,7 +36,7 @@ def from_approximate_position(
         chemical_element: Optional[str] = None,
     ):
         if not crystal:
-            RuntimeError("Crystal is not defined")
+            raise RuntimeError("Crystal is not defined")
         closest_site_id = get_closest_site_id_from_position(crystal, approximate_position)
         return cls.from_site_id(
             crystal=crystal, defect_type=defect_type, site_id=closest_site_id, chemical_element=chemical_element
@@ -50,3 +50,24 @@ def _json(self):
             "position": self.position,
             "chemical_element": self.chemical_element,
         }
+
+
+class SlabDefectConfiguration(BaseDefectConfiguration, InMemoryEntity):
+    pass
+
+
+class AdatomSlabDefectConfiguration(SlabDefectConfiguration):
+    defect_type: SlabDefectTypeEnum = SlabDefectTypeEnum.ADATOM
+    position_on_surface: List[float] = [0.5, 0.5]
+    distance_z: float = 2.0
+    chemical_element: str = "Si"
+
+    @property
+    def _json(self):
+        return {
+            "type": "AdatomSlabDefectConfiguration",
+            "defect_type": self.defect_type.name,
+            "position_on_surface": self.position_on_surface,
+            "distance_z": self.distance_z,
+            "chemical_element": self.chemical_element,
+        }

src/py/mat3ra/made/tools/build/defect/enums.py

@@ -5,3 +5,7 @@ class PointDefectTypeEnum(str, Enum):
     VACANCY = "vacancy"
     SUBSTITUTION = "substitution"
     INTERSTITIAL = "interstitial"
+
+
+class SlabDefectTypeEnum(str, Enum):
+    ADATOM = "adatom"

src/py/mat3ra/made/tools/third_party.py

@@ -6,6 +6,7 @@
 from pymatgen.analysis.defects.core import Interstitial as PymatgenInterstitial
 from pymatgen.analysis.defects.core import Substitution as PymatgenSubstitution
 from pymatgen.analysis.defects.core import Vacancy as PymatgenVacancy
+from pymatgen.analysis.local_env import VoronoiNN as PymatgenVoronoiNN
 from pymatgen.core import IStructure as PymatgenIStructure
 from pymatgen.core import PeriodicSite as PymatgenPeriodicSite
 from pymatgen.core.interface import Interface as PymatgenInterface
@@ -40,4 +41,5 @@
     "ase_add_vacuum",
     "PymatgenAseAtomsAdaptor",
     "PymatgenPoscar",
+    "PymatgenVoronoiNN",
 ]

src/py/mat3ra/made/utils.py

@@ -1,6 +1,7 @@
 import json
 from typing import Any, Callable, Dict, List, Optional, Union
 
+import numpy as np
 from mat3ra.utils.array import convert_to_array_if_not
 from mat3ra.utils.mixins import RoundNumericValuesMixin
 from pydantic import BaseModel
@@ -44,6 +45,19 @@ def are_arrays_equal_by_id_value(array1: List[Dict[str, Any]], array2: List[Dict
     return map_array_with_id_value_to_array(array1) == map_array_with_id_value_to_array(array2)
 
 
+def get_center_of_coordinates(coordinates: List[List[float]]) -> List[float]:
+    """
+    Calculate the center of the coordinates.
+
+    Args:
+        coordinates (List[List[float]]): The list of coordinates.
+
+    Returns:
+        List[float]: The center of the coordinates.
+    """
+    return list(np.mean(np.array(coordinates), axis=0))
+
+
 class ValueWithId(RoundNumericValuesMixin, BaseModel):
     id: int = 0
     value: Any = None