diff --git a/src/py/mat3ra/made/tools/analyze.py b/src/py/mat3ra/made/tools/analyze.py index d2578baa..74dd5c4a 100644 --- a/src/py/mat3ra/made/tools/analyze.py +++ b/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", diff --git a/src/py/mat3ra/made/tools/build/__init__.py b/src/py/mat3ra/made/tools/build/__init__.py index 99e10373..c1e5f2c9 100644 --- a/src/py/mat3ra/made/tools/build/__init__.py +++ b/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 diff --git a/src/py/mat3ra/made/tools/build/defect/__init__.py b/src/py/mat3ra/made/tools/build/defect/__init__.py index d447e72c..3cc599fd 100644 --- a/src/py/mat3ra/made/tools/build/defect/__init__.py +++ b/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) diff --git a/src/py/mat3ra/made/tools/build/defect/builders.py b/src/py/mat3ra/made/tools/build/defect/builders.py index e1582f28..3c9988c2 100644 --- a/src/py/mat3ra/made/tools/build/defect/builders.py +++ b/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 diff --git a/src/py/mat3ra/made/tools/build/defect/configuration.py b/src/py/mat3ra/made/tools/build/defect/configuration.py index c0897a92..2deff393 100644 --- a/src/py/mat3ra/made/tools/build/defect/configuration.py +++ b/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, + } diff --git a/src/py/mat3ra/made/tools/build/defect/enums.py b/src/py/mat3ra/made/tools/build/defect/enums.py index 8e9c092f..da3b92fd 100644 --- a/src/py/mat3ra/made/tools/build/defect/enums.py +++ b/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" diff --git a/src/py/mat3ra/made/tools/third_party.py b/src/py/mat3ra/made/tools/third_party.py index f1f86695..9013c7de 100644 --- a/src/py/mat3ra/made/tools/third_party.py +++ b/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", ] diff --git a/src/py/mat3ra/made/utils.py b/src/py/mat3ra/made/utils.py index f5d5c1fa..9e002125 100644 --- a/src/py/mat3ra/made/utils.py +++ b/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 diff --git a/tests/py/unit/test_tools_build_defect.py b/tests/py/unit/test_tools_build_defect.py index bb8423e2..050c9896 100644 --- a/tests/py/unit/test_tools_build_defect.py +++ b/tests/py/unit/test_tools_build_defect.py @@ -1,8 +1,28 @@ from mat3ra.made.material import Material -from mat3ra.made.tools.build.defect import PointDefectBuilderParameters, PointDefectConfiguration, create_defect +from mat3ra.made.tools.build.defect import ( + AdatomSlabDefectConfiguration, + EquidistantAdatomSlabDefectBuilder, + PointDefectBuilderParameters, + PointDefectConfiguration, + create_defect, + create_slab_defect, +) +from mat3ra.made.tools.build.slab import SlabConfiguration, create_slab, get_terminations +from mat3ra.utils import assertion as assertion_utils clean_material = Material.create(Material.default_config) +slab_config = SlabConfiguration( + clean_material, + (1, 1, 1), + thickness=3, + vacuum=6, + use_orthogonal_z=True, + xy_supercell_matrix=[[2, 0, 0], [0, 2, 0], [0, 0, 1]], +) +t = get_terminations(slab_config)[0] +slab = create_slab(slab_config, t) + def test_create_vacancy(): # vacancy in place of 0 element @@ -49,3 +69,26 @@ def test_create_defect_from_site_id(): {"id": 0, "value": "Si"}, {"id": 1, "value": "Ge"}, ] + + +def test_create_adatom(): + # Adatom of Si at 0.5, 0.5 position + configuration = AdatomSlabDefectConfiguration( + crystal=slab, position_on_surface=[0.5, 0.5], distance_z=2, chemical_element="Si" + ) + defect = create_slab_defect(configuration=configuration, builder=None) + + assert defect.basis.elements.values[-1] == "Si" + assertion_utils.assert_deep_almost_equal([0.5, 0.5, 0.389826], defect.basis.coordinates.values[-1]) + + +def test_create_adatom_equidistant(): + # Adatom of Si at approximate 0.5, 0.5 position + configuration = AdatomSlabDefectConfiguration( + crystal=slab, position_on_surface=[0.5, 0.5], distance_z=2, chemical_element="Si" + ) + defect = create_slab_defect(configuration=configuration, builder=EquidistantAdatomSlabDefectBuilder()) + + assert defect.basis.elements.values[-1] == "Si" + # We expect adatom to shift from provided position + assertion_utils.assert_deep_almost_equal([0.4583333333, 0.541666667, 0.389826], defect.basis.coordinates.values[-1])