diff --git a/src/py/mat3ra/made/basis.py b/src/py/mat3ra/made/basis.py index 298e9af14..da34c13ca 100644 --- a/src/py/mat3ra/made/basis.py +++ b/src/py/mat3ra/made/basis.py @@ -24,14 +24,14 @@ def from_dict( coordinates: List[Dict], units: str, labels: Optional[List[Dict]] = None, - cell: Optional[Dict] = None, + cell: Optional[List[List[float]]] = None, constraints: Optional[List[Dict]] = None, ) -> "Basis": return Basis( elements=ArrayWithIds.from_list_of_dicts(elements), coordinates=ArrayWithIds.from_list_of_dicts(coordinates), units=units, - cell=Cell.from_nested_array(cell), + cell=Cell.from_vectors_array(cell), labels=ArrayWithIds.from_list_of_dicts(labels) if labels else ArrayWithIds(values=[]), constraints=ArrayWithIds.from_list_of_dicts(constraints) if constraints else ArrayWithIds(values=[]), ) diff --git a/src/py/mat3ra/made/cell.py b/src/py/mat3ra/made/cell.py index 0f65bdbb7..01dd95581 100644 --- a/src/py/mat3ra/made/cell.py +++ b/src/py/mat3ra/made/cell.py @@ -1,4 +1,4 @@ -from typing import List +from typing import List, Optional import numpy as np from mat3ra.utils.mixins import RoundNumericValuesMixin @@ -13,13 +13,13 @@ class Cell(RoundNumericValuesMixin, BaseModel): __round_precision__ = 6 @classmethod - def from_nested_array(cls, nested_array): - if nested_array is None: - nested_array = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] - return cls(vector1=nested_array[0], vector2=nested_array[1], vector3=nested_array[2]) + def from_vectors_array(cls, vectors_array: Optional[List[List[float]]] = None) -> "Cell": + if vectors_array is None: + vectors_array = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] + return cls(vector1=vectors_array[0], vector2=vectors_array[1], vector3=vectors_array[2]) @property - def vectors_as_nested_array(self, skip_rounding=False) -> List[List[float]]: + def vectors_as_array(self, skip_rounding=False) -> List[List[float]]: if skip_rounding: return [self.vector1, self.vector2, self.vector3] return self.round_array_or_number([self.vector1, self.vector2, self.vector3]) @@ -32,22 +32,22 @@ def to_json(self, skip_rounding=False): self.vector3 if skip_rounding else _(self.vector3), ] - def clone(self): - return self.from_nested_array(self.vectors_as_nested_array) + def clone(self) -> "Cell": + return self.from_vectors_array(self.vectors_as_array) - def clone_and_scale_by_matrix(self, matrix): + def clone_and_scale_by_matrix(self, matrix: List[List[float]]) -> "Cell": new_cell = self.clone() new_cell.scale_by_matrix(matrix) return new_cell - def convert_point_to_cartesian(self, point): - np_vector = np.array(self.vectors_as_nested_array) + def convert_point_to_cartesian(self, point: List[float]) -> List[float]: + np_vector = np.array(self.vectors_as_array) return np.dot(point, np_vector) - def convert_point_to_crystal(self, point): - np_vector = np.array(self.vectors_as_nested_array) + def convert_point_to_crystal(self, point: List[float]) -> List[float]: + np_vector = np.array(self.vectors_as_array) return np.dot(point, np.linalg.inv(np_vector)) - def scale_by_matrix(self, matrix): - np_vector = np.array(self.vectors_as_nested_array) - self.vector1, self.vector2, self.vector3 = np.dot(matrix, np_vector).tolist() + def scale_by_matrix(self, matrix: List[List[float]]): + np_vector = np.array(self.vectors_as_array) + self.vector1, self.vector2, self.vector3 = np.dot(np.array(matrix), np_vector).tolist() diff --git a/src/py/mat3ra/made/lattice.py b/src/py/mat3ra/made/lattice.py index 3cf7ff579..e6fe59a46 100644 --- a/src/py/mat3ra/made/lattice.py +++ b/src/py/mat3ra/made/lattice.py @@ -8,6 +8,8 @@ from .cell import Cell HASH_TOLERANCE = 3 +DEFAULT_UNITS = {"length": "angstrom", "angle": "degree"} +DEFAULT_TYPE = "TRI" class Lattice(RoundNumericValuesMixin, BaseModel): @@ -17,11 +19,8 @@ class Lattice(RoundNumericValuesMixin, BaseModel): alpha: float = 90.0 beta: float = 90.0 gamma: float = 90.0 - units: Dict[str, str] = { - "length": "angstrom", - "angle": "degree", - } - type: str = "TRI" + units: Dict[str, str] = DEFAULT_UNITS + type: str = DEFAULT_TYPE @property def vectors(self) -> List[List[float]]: @@ -52,6 +51,29 @@ def vectors(self) -> List[List[float]]: [0.0, 0.0, c], ] + @classmethod + def from_vectors_array( + cls, vectors: List[List[float]], units: Optional[Dict[str, str]] = None, type: Optional[str] = None + ) -> "Lattice": + """ + Create a Lattice object from a nested array of vectors. + Args: + vectors (List[List[float]]): A nested array of vectors. + Returns: + Lattice: A Lattice object. + """ + a = np.linalg.norm(vectors[0]) + b = np.linalg.norm(vectors[1]) + c = np.linalg.norm(vectors[2]) + alpha = np.degrees(np.arccos(np.dot(vectors[1], vectors[2]) / (b * c))) + beta = np.degrees(np.arccos(np.dot(vectors[0], vectors[2]) / (a * c))) + gamma = np.degrees(np.arccos(np.dot(vectors[0], vectors[1]) / (a * b))) + if units is None: + units = DEFAULT_UNITS + if type is None: + type = DEFAULT_TYPE + return cls(a=float(a), b=float(b), c=float(c), alpha=alpha, beta=beta, gamma=gamma, units=units, type=type) + def to_json(self, skip_rounding: bool = False) -> Dict[str, Any]: __round__ = RoundNumericValuesMixin.round_array_or_number round_func = __round__ if not skip_rounding else lambda x: x @@ -78,7 +100,7 @@ def vector_arrays(self) -> List[List[float]]: @property def cell(self) -> Cell: - return Cell.from_nested_array(self.vector_arrays) + return Cell.from_vectors_array(self.vector_arrays) def volume(self) -> float: np_vector = np.array(self.vector_arrays) diff --git a/src/py/mat3ra/made/material.py b/src/py/mat3ra/made/material.py index 9b50ba47e..1a5d01625 100644 --- a/src/py/mat3ra/made/material.py +++ b/src/py/mat3ra/made/material.py @@ -90,3 +90,8 @@ def to_crystal(self) -> None: new_basis = self.basis.copy() new_basis.to_crystal() self.basis = new_basis + + def set_coordinates(self, coordinates: List[List[float]]) -> None: + new_basis = self.basis.copy() + new_basis.coordinates.values = coordinates + self.basis = new_basis diff --git a/src/py/mat3ra/made/tools/build/defect/builders.py b/src/py/mat3ra/made/tools/build/defect/builders.py index 7d8ee877b..a9c25abd1 100644 --- a/src/py/mat3ra/made/tools/build/defect/builders.py +++ b/src/py/mat3ra/made/tools/build/defect/builders.py @@ -30,7 +30,8 @@ get_closest_site_id_from_coordinate_and_element, ) from ....utils import get_center_of_coordinates -from ...utils import transform_coordinate_to_supercell, CoordinateConditionBuilder +from ...utils import transform_coordinate_to_supercell +from ...utils import coordinate as CoordinateCondition from ..utils import merge_materials from ..slab import SlabConfiguration, create_slab, Termination from ..supercell import create_supercell @@ -436,7 +437,7 @@ def condition(coordinate: List[float]): return self.merge_slab_and_defect(island_material, new_material) def _generate(self, configuration: _ConfigurationType) -> List[_GeneratedItemType]: - condition_callable, _ = configuration.condition + condition_callable = configuration.condition.condition return [ self.create_island( material=configuration.crystal, @@ -463,7 +464,7 @@ def _calculate_cut_direction_vector(self, material: Material, cut_direction: Lis The normalized cut direction vector in Cartesian coordinates. """ np_cut_direction = np.array(cut_direction) - direction_vector = np.dot(np.array(material.basis.cell.vectors_as_nested_array), np_cut_direction) + direction_vector = np.dot(np.array(material.basis.cell.vectors_as_array), np_cut_direction) normalized_direction_vector = direction_vector / np.linalg.norm(direction_vector) return normalized_direction_vector @@ -499,7 +500,7 @@ def _calculate_rotation_parameters( """ height_cartesian = self._calculate_height_cartesian(original_material, new_material) cut_direction_xy_proj_cart = np.linalg.norm( - np.dot(np.array(new_material.basis.cell.vectors_as_nested_array), normalized_direction_vector) + np.dot(np.array(new_material.basis.cell.vectors_as_array), normalized_direction_vector) ) # Slope of the terrace along the cut direction hypotenuse = np.linalg.norm([height_cartesian, cut_direction_xy_proj_cart]) @@ -592,10 +593,10 @@ def create_terrace( ) normalized_direction_vector = self._calculate_cut_direction_vector(material, cut_direction) - condition, _ = CoordinateConditionBuilder.plane( + condition = CoordinateCondition.PlaneCoordinateCondition( plane_normal=normalized_direction_vector, plane_point_coordinate=pivot_coordinate, - ) + ).condition atoms_within_terrace = filter_by_condition_on_coordinates( material=material_with_additional_layers, condition=condition, diff --git a/src/py/mat3ra/made/tools/build/defect/configuration.py b/src/py/mat3ra/made/tools/build/defect/configuration.py index dcd33e338..43e49f250 100644 --- a/src/py/mat3ra/made/tools/build/defect/configuration.py +++ b/src/py/mat3ra/made/tools/build/defect/configuration.py @@ -1,17 +1,25 @@ -from typing import Optional, List, Any, Callable, Dict, Tuple, Union +from typing import Optional, List, Union from pydantic import BaseModel from mat3ra.code.entity import InMemoryEntity from mat3ra.made.material import Material from ...analyze import get_closest_site_id_from_coordinate, get_atomic_coordinates_extremum -from ...utils import CoordinateConditionBuilder +from ...utils.coordinate import ( + CylinderCoordinateCondition, + SphereCoordinateCondition, + BoxCoordinateCondition, + TriangularPrismCoordinateCondition, + PlaneCoordinateCondition, +) from .enums import PointDefectTypeEnum, SlabDefectTypeEnum, AtomPlacementMethodEnum, ComplexDefectTypeEnum class BaseDefectConfiguration(BaseModel): - # TODO: fix arbitrary_types_allowed error and set Material class type - crystal: Any = None + crystal: Material = None + + class Config: + arbitrary_types_allowed = True @property def _json(self): @@ -169,16 +177,21 @@ class IslandSlabDefectConfiguration(SlabDefectConfiguration): """ defect_type: SlabDefectTypeEnum = SlabDefectTypeEnum.ISLAND - condition: Optional[Tuple[Callable[[List[float]], bool], Dict]] = CoordinateConditionBuilder().cylinder() + condition: Union[ + CylinderCoordinateCondition, + SphereCoordinateCondition, + BoxCoordinateCondition, + TriangularPrismCoordinateCondition, + PlaneCoordinateCondition, + ] = CylinderCoordinateCondition() @property def _json(self): - _, condition_json = self.condition return { **super()._json, "type": self.get_cls_name(), "defect_type": self.defect_type.name, - "condition": condition_json, + "condition": self.condition.get_json(), } diff --git a/src/py/mat3ra/made/tools/build/interface/builders.py b/src/py/mat3ra/made/tools/build/interface/builders.py index 8c2642ebd..a02abbdf3 100644 --- a/src/py/mat3ra/made/tools/build/interface/builders.py +++ b/src/py/mat3ra/made/tools/build/interface/builders.py @@ -1,6 +1,7 @@ from typing import Any, List, Optional import numpy as np +from mat3ra.made.tools.modify import translate_to_z_level from pydantic import BaseModel from ase.build.tools import niggli_reduce from pymatgen.analysis.interfaces.coherent_interfaces import ( @@ -144,9 +145,13 @@ class ZSLStrainMatchingInterfaceBuilder(ConvertGeneratedItemsPymatgenStructureMi def _generate(self, configuration: InterfaceConfiguration) -> List[PymatgenInterface]: generator = ZSLGenerator(**self.build_parameters.strain_matching_parameters.dict()) + substrate_with_atoms_translated_to_bottom = translate_to_z_level( + configuration.substrate_configuration.bulk, "bottom" + ) + film_with_atoms_translated_to_bottom = translate_to_z_level(configuration.film_configuration.bulk, "bottom") builder = CoherentInterfaceBuilder( - substrate_structure=to_pymatgen(configuration.substrate_configuration.bulk), - film_structure=to_pymatgen(configuration.film_configuration.bulk), + substrate_structure=to_pymatgen(substrate_with_atoms_translated_to_bottom), + film_structure=to_pymatgen(film_with_atoms_translated_to_bottom), substrate_miller=configuration.substrate_configuration.miller_indices, film_miller=configuration.film_configuration.miller_indices, zslgen=generator, diff --git a/src/py/mat3ra/made/tools/build/perturbation/__init__.py b/src/py/mat3ra/made/tools/build/perturbation/__init__.py new file mode 100644 index 000000000..4e799c31c --- /dev/null +++ b/src/py/mat3ra/made/tools/build/perturbation/__init__.py @@ -0,0 +1,37 @@ +from typing import Union, Optional + +from mat3ra.made.material import Material +from .builders import ( + SlabPerturbationBuilder, + DistancePreservingSlabPerturbationBuilder, + CellMatchingDistancePreservingSlabPerturbationBuilder, +) +from .configuration import PerturbationConfiguration + + +def create_perturbation( + configuration: PerturbationConfiguration, + preserve_distance: Optional[bool] = False, + builder: Union[ + SlabPerturbationBuilder, + DistancePreservingSlabPerturbationBuilder, + CellMatchingDistancePreservingSlabPerturbationBuilder, + None, + ] = None, +) -> Material: + """ + Return a material with a perturbation applied. + + Args: + configuration: The configuration of the perturbation to be applied. + preserve_distance: If True, the builder that preserves the distance between atoms is used. + builder: The builder to be used to create the perturbation. + + Returns: + The material with the perturbation applied. + """ + if builder is None: + builder = SlabPerturbationBuilder() + if preserve_distance: + builder = CellMatchingDistancePreservingSlabPerturbationBuilder() + return builder.get_material(configuration) diff --git a/src/py/mat3ra/made/tools/build/perturbation/builders.py b/src/py/mat3ra/made/tools/build/perturbation/builders.py new file mode 100644 index 000000000..f13ba6d33 --- /dev/null +++ b/src/py/mat3ra/made/tools/build/perturbation/builders.py @@ -0,0 +1,88 @@ +from typing import List, Optional, Any + +from mat3ra.made.material import Material +from mat3ra.made.tools.build import BaseBuilder + +from .configuration import PerturbationConfiguration +from ...modify import wrap_to_unit_cell, translate_to_z_level + + +class PerturbationBuilder(BaseBuilder): + _ConfigurationType: type(PerturbationConfiguration) = PerturbationConfiguration # type: ignore + _GeneratedItemType: Material = Material + _PostProcessParametersType: Any = None + + @staticmethod + def _prepare_material(configuration: _ConfigurationType) -> _GeneratedItemType: + new_material = configuration.material.clone() + new_material = translate_to_z_level(new_material, "center") + if configuration.use_cartesian_coordinates: + new_material.to_cartesian() + return new_material + + def _generate(self, configuration: _ConfigurationType) -> List[_GeneratedItemType]: + """Generate materials with applied continuous perturbation based on the given configuration.""" + new_material = self.create_perturbed_slab(configuration) + return [new_material] + + def _post_process( + self, + items: List[_GeneratedItemType], + post_process_parameters: Optional[_PostProcessParametersType], + ) -> List[Material]: + return [wrap_to_unit_cell(item) for item in items] + + def _update_material_name(self, material: Material, configuration: _ConfigurationType) -> Material: + perturbation_details = f"Perturbation: {configuration.perturbation_function_holder.get_json().get('type')}" + material.name = f"{material.name} ({perturbation_details})" + return material + + +class SlabPerturbationBuilder(PerturbationBuilder): + def create_perturbed_slab(self, configuration: PerturbationConfiguration) -> Material: + new_material = self._prepare_material(configuration) + new_coordinates = [ + configuration.perturbation_function_holder.apply_perturbation(coord) + for coord in new_material.basis.coordinates.values + ] + new_material.set_coordinates(new_coordinates) + return new_material + + +class DistancePreservingSlabPerturbationBuilder(SlabPerturbationBuilder): + def create_perturbed_slab(self, configuration: PerturbationConfiguration) -> Material: + new_material = self._prepare_material(configuration) + new_coordinates = [ + configuration.perturbation_function_holder.transform_coordinates(coord) + for coord in new_material.basis.coordinates.values + ] + new_coordinates = [ + configuration.perturbation_function_holder.apply_perturbation(coord) for coord in new_coordinates + ] + new_material.set_coordinates(new_coordinates) + return new_material + + +class CellMatchingDistancePreservingSlabPerturbationBuilder(DistancePreservingSlabPerturbationBuilder): + def _transform_lattice_vectors(self, configuration: PerturbationConfiguration) -> List[List[float]]: + cell_vectors = configuration.material.basis.cell.vectors_as_array + return [ + configuration.perturbation_function_holder.transform_coordinates( + configuration.perturbation_function_holder.transform_coordinates(coord) + ) + for coord in cell_vectors + ] + + def create_perturbed_slab(self, configuration: PerturbationConfiguration) -> Material: + new_material = super().create_perturbed_slab(configuration) + new_lattice_vectors = self._transform_lattice_vectors(configuration) + new_lattice = new_material.lattice.copy() + new_lattice = new_lattice.from_vectors_array(new_lattice_vectors) + new_material.lattice = new_lattice + + new_basis = new_material.basis.copy() + new_basis.to_cartesian() + new_basis.cell = new_basis.cell.from_vectors_array(new_lattice_vectors) + new_basis.to_crystal() + new_material.basis = new_basis + return new_material diff --git a/src/py/mat3ra/made/tools/build/perturbation/configuration.py b/src/py/mat3ra/made/tools/build/perturbation/configuration.py new file mode 100644 index 000000000..767ffa8a9 --- /dev/null +++ b/src/py/mat3ra/made/tools/build/perturbation/configuration.py @@ -0,0 +1,24 @@ +from mat3ra.code.entity import InMemoryEntity +from mat3ra.made.material import Material +from pydantic import BaseModel + +from ...utils.perturbation import SineWavePerturbationFunctionHolder + + +class PerturbationConfiguration(BaseModel, InMemoryEntity): + material: Material + perturbation_function_holder: SineWavePerturbationFunctionHolder = SineWavePerturbationFunctionHolder() + use_cartesian_coordinates: bool = True + + class Config: + arbitrary_types_allowed = True + + @property + def _json(self): + perturbation_function_json = self.perturbation_function_holder.get_json() + return { + "type": self.get_cls_name(), + "material": self.material.to_json(), + "perturbation_function": perturbation_function_json, + "use_cartesian_coordinates": self.use_cartesian_coordinates, + } diff --git a/src/py/mat3ra/made/tools/build/supercell.py b/src/py/mat3ra/made/tools/build/supercell.py index 52de72a5c..1ab80e1d1 100644 --- a/src/py/mat3ra/made/tools/build/supercell.py +++ b/src/py/mat3ra/made/tools/build/supercell.py @@ -2,28 +2,32 @@ import numpy as np from mat3ra.made.material import Material -from ..third_party import ASEAtoms, ase_make_supercell +from ..third_party import ase_make_supercell from ..utils import decorator_convert_2x2_to_3x3 -from ..convert import from_ase, decorator_convert_material_args_kwargs_to_atoms +from ..convert import from_ase, to_ase @decorator_convert_2x2_to_3x3 -@decorator_convert_material_args_kwargs_to_atoms def create_supercell( - atoms: ASEAtoms, supercell_matrix: Optional[List[List[int]]] = None, scaling_factor: Optional[List[int]] = None + material: Material, supercell_matrix: Optional[List[List[int]]] = None, scaling_factor: Optional[List[int]] = None ) -> Material: """ Create a supercell of the atoms. Args: - atoms (Material): The atoms to create a supercell of. + material (Material): The atoms to create a supercell of. supercell_matrix (List[List[int]]): The supercell matrix (e.g. [[3,0,0],[0,3,0],[0,0,1]]). scaling_factor (List[int], optional): The scaling factor instead of matrix (e.g. [3,3,1]). Returns: Material: The supercell of the atoms. """ + atoms = to_ase(material) if scaling_factor is not None: supercell_matrix = np.multiply(scaling_factor, np.eye(3)).tolist() supercell_atoms = ase_make_supercell(atoms, supercell_matrix) - return Material(from_ase(supercell_atoms)) + new_material = Material(from_ase(supercell_atoms)) + if material.metadata: + new_material.metadata = material.metadata + new_material.name = material.name + return new_material diff --git a/src/py/mat3ra/made/tools/modify.py b/src/py/mat3ra/made/tools/modify.py index 983640bd0..daefdd6b1 100644 --- a/src/py/mat3ra/made/tools/modify.py +++ b/src/py/mat3ra/made/tools/modify.py @@ -7,9 +7,9 @@ get_atom_indices_within_radius_pbc, get_atomic_coordinates_extremum, ) -from .convert import decorator_convert_material_args_kwargs_to_structure, from_ase, to_ase -from .third_party import PymatgenStructure, ase_add_vacuum -from .utils import ( +from .convert import from_ase, to_ase +from .third_party import ase_add_vacuum +from .utils.coordinate import ( is_coordinate_in_box, is_coordinate_in_cylinder, is_coordinate_in_triangular_prism, @@ -87,18 +87,18 @@ def translate_by_vector( return Material(from_ase(atoms)) -@decorator_convert_material_args_kwargs_to_structure -def wrap_to_unit_cell(structure: PymatgenStructure): +def wrap_to_unit_cell(material: Material) -> Material: """ - Wrap atoms to the cell + Wrap the material to the unit cell. Args: - structure (PymatgenStructure): The pymatgen PymatgenStructure object to normalize. + material (Material): The material to wrap. Returns: - PymatgenStructure: The wrapped pymatgen PymatgenStructure object. + Material: The wrapped material. """ - structure.make_supercell((1, 1, 1), to_unit_cell=True) - return structure + atoms = to_ase(material) + atoms.wrap() + return Material(from_ase(atoms)) def filter_material_by_ids(material: Material, ids: List[int], invert: bool = False) -> Material: diff --git a/src/py/mat3ra/made/tools/utils/__init__.py b/src/py/mat3ra/made/tools/utils/__init__.py new file mode 100644 index 000000000..6eb87397a --- /dev/null +++ b/src/py/mat3ra/made/tools/utils/__init__.py @@ -0,0 +1,109 @@ +from functools import wraps +from typing import Callable, List, Optional + +import numpy as np +from mat3ra.utils.matrix import convert_2x2_to_3x3 + +from ..third_party import PymatgenStructure +from .coordinate import ( + is_coordinate_behind_plane, + is_coordinate_in_box, + is_coordinate_in_cylinder, + is_coordinate_in_sphere, + is_coordinate_in_triangular_prism, +) +from .factories import PerturbationFunctionHolderFactory +from .functions import AXIS_TO_INDEX_MAP + +DEFAULT_SCALING_FACTOR = np.array([3, 3, 3]) +DEFAULT_TRANSLATION_VECTOR = 1 / DEFAULT_SCALING_FACTOR + + +# TODO: convert to accept ASE Atoms object +def translate_to_bottom_pymatgen_structure(structure: PymatgenStructure): + """ + Translate the structure to the bottom of the cell. + Args: + structure (PymatgenStructure): The pymatgen Structure object to translate. + + Returns: + PymatgenStructure: The translated pymatgen Structure object. + """ + min_c = min(site.c for site in structure) + translation_vector = [0, 0, -min_c] + translated_structure = structure.copy() + for site in translated_structure: + site.coords += translation_vector + return translated_structure + + +def decorator_convert_2x2_to_3x3(func: Callable) -> Callable: + """ + Decorator to convert a 2x2 matrix to a 3x3 matrix. + """ + + @wraps(func) + def wrapper(*args, **kwargs): + new_args = [convert_2x2_to_3x3(arg) if isinstance(arg, list) and len(arg) == 2 else arg for arg in args] + return func(*new_args, **kwargs) + + return wrapper + + +def get_distance_between_coordinates(coordinate1: List[float], coordinate2: List[float]) -> float: + """ + Get the distance between two coordinates. + Args: + coordinate1 (List[float]): The first coordinate. + coordinate2 (List[float]): The second coordinate. + + Returns: + float: The distance between the two coordinates. + """ + return float(np.linalg.norm(np.array(coordinate1) - np.array(coordinate2))) + + +def get_norm(vector: List[float]) -> float: + """ + Get the norm of a vector. + Args: + vector (List[float]): The vector. + + Returns: + float: The norm of the vector. + """ + return float(np.linalg.norm(vector)) + + +def transform_coordinate_to_supercell( + coordinate: List[float], + scaling_factor: Optional[List[int]] = None, + translation_vector: Optional[List[float]] = None, + reverse: bool = False, +) -> List[float]: + """ + Convert a crystal coordinate of unit cell to a coordinate in a supercell. + Args: + coordinate (List[float]): The coordinates to convert. + scaling_factor (List[int]): The scaling factor for the supercell. + translation_vector (List[float]): The translation vector for the supercell. + reverse (bool): Whether to convert in the reverse transformation. + + Returns: + List[float]: The converted coordinates. + """ + if scaling_factor is None: + np_scaling_factor = np.array([3, 3, 3]) + else: + np_scaling_factor = np.array(scaling_factor) + + if translation_vector is None: + np_translation_vector = np.array([0, 0, 0]) + else: + np_translation_vector = np.array(translation_vector) + + np_coordinate = np.array(coordinate) + converted_array = np_coordinate * (1 / np_scaling_factor) + np_translation_vector + if reverse: + converted_array = (np_coordinate - np_translation_vector) * np_scaling_factor + return converted_array.tolist() diff --git a/src/py/mat3ra/made/tools/utils.py b/src/py/mat3ra/made/tools/utils/coordinate.py similarity index 50% rename from src/py/mat3ra/made/tools/utils.py rename to src/py/mat3ra/made/tools/utils/coordinate.py index 88ef35c34..683c54b38 100644 --- a/src/py/mat3ra/made/tools/utils.py +++ b/src/py/mat3ra/made/tools/utils/coordinate.py @@ -1,72 +1,8 @@ -from functools import wraps -from typing import Callable, Dict, List, Optional, Tuple +# Place all functions acting on coordinates +from typing import Dict, List import numpy as np -from mat3ra.utils.matrix import convert_2x2_to_3x3 - -from .third_party import PymatgenStructure - -DEFAULT_SCALING_FACTOR = np.array([3, 3, 3]) -DEFAULT_TRANSLATION_VECTOR = 1 / DEFAULT_SCALING_FACTOR - - -# TODO: convert to accept ASE Atoms object -def translate_to_bottom_pymatgen_structure(structure: PymatgenStructure): - """ - Translate the structure to the bottom of the cell. - Args: - structure (PymatgenStructure): The pymatgen Structure object to translate. - - Returns: - PymatgenStructure: The translated pymatgen Structure object. - """ - min_c = min(site.c for site in structure) - translation_vector = [0, 0, -min_c] - translated_structure = structure.copy() - for site in translated_structure: - site.coords += translation_vector - return translated_structure - - -def decorator_convert_2x2_to_3x3(func: Callable) -> Callable: - """ - Decorator to convert a 2x2 matrix to a 3x3 matrix. - """ - - @wraps(func) - def wrapper(*args, **kwargs): - new_args = [convert_2x2_to_3x3(arg) if isinstance(arg, list) and len(arg) == 2 else arg for arg in args] - return func(*new_args, **kwargs) - - return wrapper - - -def get_distance_between_coordinates(coordinate1: List[float], coordinate2: List[float]) -> float: - """ - Get the distance between two coordinates. - Args: - coordinate1 (List[float]): The first coordinate. - coordinate2 (List[float]): The second coordinate. - - Returns: - float: The distance between the two coordinates. - """ - return float(np.linalg.norm(np.array(coordinate1) - np.array(coordinate2))) - - -def get_norm(vector: List[float]) -> float: - """ - Get the norm of a vector. - Args: - vector (List[float]): The vector. - - Returns: - float: The norm of the vector. - """ - return float(np.linalg.norm(vector)) - - -# Condition functions: +from pydantic import BaseModel, Field def is_coordinate_in_cylinder( @@ -219,106 +155,63 @@ def is_coordinate_behind_plane( return np.dot(np_plane_normal, np_coordinate - np_plane_point) < 0 -def transform_coordinate_to_supercell( - coordinate: List[float], - scaling_factor: Optional[List[int]] = None, - translation_vector: Optional[List[float]] = None, - reverse: bool = False, -) -> List[float]: - """ - Convert a crystal coordinate of unit cell to a coordinate in a supercell. - Args: - coordinate (List[float]): The coordinates to convert. - scaling_factor (List[int]): The scaling factor for the supercell. - translation_vector (List[float]): The translation vector for the supercell. - reverse (bool): Whether to convert in the reverse transformation. +class CoordinateCondition(BaseModel): + def condition(self, coordinate: List[float]) -> bool: + raise NotImplementedError - Returns: - List[float]: The converted coordinates. - """ - if scaling_factor is None: - np_scaling_factor = np.array([3, 3, 3]) - else: - np_scaling_factor = np.array(scaling_factor) + def get_json(self) -> Dict: + json = {"type": self.__class__.__name__} + json.update(self.dict()) + return json - if translation_vector is None: - np_translation_vector = np.array([0, 0, 0]) - else: - np_translation_vector = np.array(translation_vector) - np_coordinate = np.array(coordinate) - converted_array = np_coordinate * (1 / np_scaling_factor) + np_translation_vector - if reverse: - converted_array = (np_coordinate - np_translation_vector) * np_scaling_factor - return converted_array.tolist() - - -class CoordinateConditionBuilder: - @staticmethod - def create_condition(condition_type: str, evaluation_func: Callable, **kwargs) -> Tuple[Callable, Dict]: - condition_json = {"type": condition_type, **kwargs} - return lambda coordinate: evaluation_func(coordinate, **kwargs), condition_json - - @staticmethod - def cylinder(center_position=None, radius: float = 0.25, min_z: float = 0, max_z: float = 1): - if center_position is None: - center_position = [0.5, 0.5] - return CoordinateConditionBuilder.create_condition( - condition_type="cylinder", - evaluation_func=is_coordinate_in_cylinder, - center_position=center_position, - radius=radius, - min_z=min_z, - max_z=max_z, - ) +class CylinderCoordinateCondition(CoordinateCondition): + center_position: List[float] = Field(default_factory=lambda: [0.5, 0.5]) + radius: float = 0.25 + min_z: float = 0 + max_z: float = 1 - @staticmethod - def sphere(center_position=None, radius: float = 0.25): - if center_position is None: - center_position = [0.5, 0.5, 0.5] - return CoordinateConditionBuilder.create_condition( - condition_type="sphere", - evaluation_func=is_coordinate_in_sphere, - center_position=center_position, - radius=radius, - ) + def condition(self, coordinate: List[float]) -> bool: + return is_coordinate_in_cylinder(coordinate, self.center_position, self.radius, self.min_z, self.max_z) - @staticmethod - def triangular_prism( - position_on_surface_1: List[float] = [0, 0], - position_on_surface_2: List[float] = [1, 0], - position_on_surface_3: List[float] = [0, 1], - min_z: float = 0, - max_z: float = 1, - ): - return CoordinateConditionBuilder.create_condition( - condition_type="prism", - evaluation_func=is_coordinate_in_triangular_prism, - coordinate_1=position_on_surface_1, - coordinate_2=position_on_surface_2, - coordinate_3=position_on_surface_3, - min_z=min_z, - max_z=max_z, - ) - @staticmethod - def box(min_coordinate=None, max_coordinate=None): - if max_coordinate is None: - max_coordinate = [1, 1, 1] - if min_coordinate is None: - min_coordinate = [0, 0, 0] - return CoordinateConditionBuilder.create_condition( - condition_type="box", - evaluation_func=is_coordinate_in_box, - min_coordinate=min_coordinate, - max_coordinate=max_coordinate, - ) +class SphereCoordinateCondition(CoordinateCondition): + center_position: List[float] = Field(default_factory=lambda: [0.5, 0.5]) + radius: float = 0.25 + + def condition(self, coordinate: List[float]) -> bool: + return is_coordinate_in_sphere(coordinate, self.center_position, self.radius) + + +class BoxCoordinateCondition(CoordinateCondition): + min_coordinate: List[float] = Field(default_factory=lambda: [0, 0, 0]) + max_coordinate: List[float] = Field(default_factory=lambda: [1, 1, 1]) - @staticmethod - def plane(plane_normal: List[float], plane_point_coordinate: List[float]): - return CoordinateConditionBuilder.create_condition( - condition_type="plane", - evaluation_func=is_coordinate_behind_plane, - plane_normal=plane_normal, - plane_point_coordinate=plane_point_coordinate, + def condition(self, coordinate: List[float]) -> bool: + return is_coordinate_in_box(coordinate, self.min_coordinate, self.max_coordinate) + + +class TriangularPrismCoordinateCondition(CoordinateCondition): + position_on_surface_1: List[float] = [0, 0] + position_on_surface_2: List[float] = [1, 0] + position_on_surface_3: List[float] = [0, 1] + min_z: float = 0 + max_z: float = 1 + + def condition(self, coordinate: List[float]) -> bool: + return is_coordinate_in_triangular_prism( + coordinate, + self.position_on_surface_1, + self.position_on_surface_2, + self.position_on_surface_3, + self.min_z, + self.max_z, ) + + +class PlaneCoordinateCondition(CoordinateCondition): + plane_normal: List[float] + plane_point_coordinate: List[float] + + def condition(self, coordinate: List[float]) -> bool: + return is_coordinate_behind_plane(coordinate, self.plane_normal, self.plane_point_coordinate) diff --git a/src/py/mat3ra/made/tools/utils/factories.py b/src/py/mat3ra/made/tools/utils/factories.py new file mode 100644 index 000000000..22b6343f7 --- /dev/null +++ b/src/py/mat3ra/made/tools/utils/factories.py @@ -0,0 +1,7 @@ +from mat3ra.utils.factory import BaseFactory + + +class PerturbationFunctionHolderFactory(BaseFactory): + __class_registry__ = { + "sine_wave": "mat3ra.made.tools.utils.functions.SineWaveFunctionHolder", + } diff --git a/src/py/mat3ra/made/tools/utils/functions.py b/src/py/mat3ra/made/tools/utils/functions.py new file mode 100644 index 000000000..9d0ec4eaa --- /dev/null +++ b/src/py/mat3ra/made/tools/utils/functions.py @@ -0,0 +1,32 @@ +from typing import List + +from pydantic import BaseModel + +AXIS_TO_INDEX_MAP = {"x": 0, "y": 1, "z": 2} +EQUATION_RANGE_COEFFICIENT = 5 + + +class FunctionHolder(BaseModel): + def apply_function(self, coordinate: List[float]) -> float: + """ + Get the value of the function at the given coordinate. + """ + raise NotImplementedError + + def calculate_derivative(self, coordinate: List[float]) -> float: + """ + Get the derivative of the function at the given coordinate + """ + raise NotImplementedError + + def calculate_arc_length(self, a: float, b: float) -> float: + """ + Get the arc length of the function between a and b. + """ + raise NotImplementedError + + def get_json(self) -> dict: + """ + Get the json representation of the function holder. + """ + raise NotImplementedError diff --git a/src/py/mat3ra/made/tools/utils/perturbation.py b/src/py/mat3ra/made/tools/utils/perturbation.py new file mode 100644 index 000000000..8e5ff4cfd --- /dev/null +++ b/src/py/mat3ra/made/tools/utils/perturbation.py @@ -0,0 +1,77 @@ +from typing import List, Literal + +import numpy as np +from scipy.integrate import quad +from scipy.optimize import root_scalar + +from .functions import AXIS_TO_INDEX_MAP, EQUATION_RANGE_COEFFICIENT, FunctionHolder + + +class PerturbationFunctionHolder(FunctionHolder): + def calculate_arc_length_equation(self, w_prime: float, w: float) -> float: + """ + Get the arc length equation for the perturbation function. + """ + arc_length = quad( + lambda t: np.sqrt(1 + (self.calculate_derivative([t]) ** 2)), + a=0, + b=w_prime, + )[0] + return arc_length - w + + def transform_coordinates(self, coordinate: List[float]) -> List[float]: + """ + Transform coordinates to preserve the distance between points on a sine wave when perturbation is applied. + Achieved by calculating the integral of the length between [0,0,0] and given coordinate. + + Returns: + Callable[[List[float]], List[float]]: The coordinates transformation function. + """ + raise NotImplementedError + + def apply_perturbation(self, coordinate: List[float]) -> List[float]: + """ + Apply the perturbation to the given coordinate. + """ + raise NotImplementedError + + +class SineWavePerturbationFunctionHolder(PerturbationFunctionHolder): + amplitude: float = 0.05 + wavelength: float = 1 + phase: float = 0 + axis: Literal["x", "y"] = "x" + + def apply_function(self, coordinate: List[float]) -> float: + w = coordinate[AXIS_TO_INDEX_MAP[self.axis]] + return self.amplitude * np.sin(2 * np.pi * w / self.wavelength + self.phase) + + def calculate_derivative(self, coordinate: List[float]) -> float: + w = coordinate[AXIS_TO_INDEX_MAP[self.axis]] + return self.amplitude * 2 * np.pi / self.wavelength * np.cos(2 * np.pi * w / self.wavelength + self.phase) + + def apply_perturbation(self, coordinate: List[float]) -> List[float]: + return [coordinate[0], coordinate[1], coordinate[2] + self.apply_function(coordinate)] + + def transform_coordinates(self, coordinate: List[float]) -> List[float]: + index = AXIS_TO_INDEX_MAP[self.axis] + + w = coordinate[index] + # Find x' such that the integral from 0 to x' equals x + result = root_scalar( + self.calculate_arc_length_equation, + args=w, + bracket=[0, EQUATION_RANGE_COEFFICIENT * w], + method="brentq", + ) + coordinate[index] = result.root + return coordinate + + def get_json(self) -> dict: + return { + "type": self.__class__.__name__, + "amplitude": self.amplitude, + "wavelength": self.wavelength, + "phase": self.phase, + "axis": self.axis, + } diff --git a/tests/py/unit/fixtures.py b/tests/py/unit/fixtures.py index 384d21a14..193a9af55 100644 --- a/tests/py/unit/fixtures.py +++ b/tests/py/unit/fixtures.py @@ -114,7 +114,7 @@ } SI_SUPERCELL_2X2X1: Dict[str, Any] = { - "name": "Si8", + "name": "Silicon FCC", "basis": { "elements": [ {"id": 0, "value": "Si"}, @@ -260,3 +260,32 @@ ) t_001 = get_terminations(slab_001_config)[0] SLAB_001 = create_slab(slab_001_config, t_001) + +GRAPHENE = { + "name": "Graphene", + "basis": { + "elements": [{"id": 0, "value": "C"}, {"id": 1, "value": "C"}], + "coordinates": [{"id": 0, "value": [0, 0, 0]}, {"id": 1, "value": [0.333333, 0.666667, 0]}], + "units": "crystal", + "cell": [[2.467291, 0, 0], [-1.2336454999, 2.1367366845, 0], [0, 0, 20]], + "constraints": [], + }, + "lattice": { + "a": 2.467291, + "b": 2.467291, + "c": 20, + "alpha": 90, + "beta": 90, + "gamma": 120, + "units": {"length": "angstrom", "angle": "degree"}, + "type": "HEX", + "vectors": { + "a": [2.467291, 0, 0], + "b": [-1.233645, 2.136737, 0], + "c": [0, 0, 20], + "alat": 1, + "units": "angstrom", + }, + }, + "isNonPeriodic": False, +} diff --git a/tests/py/unit/test_tools_build_defect.py b/tests/py/unit/test_tools_build_defect.py index 647610b57..1ee349fae 100644 --- a/tests/py/unit/test_tools_build_defect.py +++ b/tests/py/unit/test_tools_build_defect.py @@ -19,7 +19,7 @@ PointDefectPairConfiguration, TerraceSlabDefectConfiguration, ) -from mat3ra.made.tools.utils import CoordinateConditionBuilder +from mat3ra.made.tools.utils import coordinate as CoordinateCondition from mat3ra.utils import assertion as assertion_utils from .fixtures import SLAB_001, SLAB_111 @@ -114,7 +114,9 @@ def test_create_crystal_site_adatom(): def test_create_island(): - condition = CoordinateConditionBuilder.cylinder(center_position=[0.625, 0.5], radius=0.25, min_z=0, max_z=1) + condition = CoordinateCondition.CylinderCoordinateCondition( + center_position=[0.625, 0.5], radius=0.25, min_z=0, max_z=1 + ) island_config = IslandSlabDefectConfiguration( crystal=SLAB_111, defect_type="island", diff --git a/tests/py/unit/test_tools_build_perturbation.py b/tests/py/unit/test_tools_build_perturbation.py new file mode 100644 index 000000000..86232a5dc --- /dev/null +++ b/tests/py/unit/test_tools_build_perturbation.py @@ -0,0 +1,50 @@ +from mat3ra.made.cell import Cell +from mat3ra.made.material import Material +from mat3ra.made.tools.build.perturbation import create_perturbation +from mat3ra.made.tools.build.perturbation.builders import SlabPerturbationBuilder +from mat3ra.made.tools.build.perturbation.configuration import PerturbationConfiguration +from mat3ra.made.tools.build.supercell import create_supercell +from mat3ra.made.tools.utils.perturbation import SineWavePerturbationFunctionHolder +from mat3ra.utils import assertion as assertion_utils + +from .fixtures import GRAPHENE + + +def test_sine_perturbation(): + material = Material(GRAPHENE) + slab = create_supercell(material, [[10, 0, 0], [0, 10, 0], [0, 0, 1]]) + + perturbation_config = PerturbationConfiguration( + material=slab, + perturbation_function=SineWavePerturbationFunctionHolder(amplitude=0.05, wavelength=1), + use_cartesian_coordinates=False, + ) + builder = SlabPerturbationBuilder() + perturbed_slab = builder.get_material(perturbation_config) + # Check selected atoms to avoid using 100+ atoms fixture + assertion_utils.assert_deep_almost_equal([0.0, 0.0, 0.5], perturbed_slab.basis.coordinates.values[0]) + assertion_utils.assert_deep_almost_equal([0.2, 0.1, 0.547552826], perturbed_slab.basis.coordinates.values[42]) + + +def test_distance_preserved_sine_perturbation(): + material = Material(GRAPHENE) + slab = create_supercell(material, [[10, 0, 0], [0, 10, 0], [0, 0, 1]]) + + perturbation_config = PerturbationConfiguration( + material=slab, + perturbation_function=SineWavePerturbationFunctionHolder(amplitude=0.05, wavelength=1, phase=0.25, axis="y"), + use_cartesian_coordinates=False, + ) + perturbed_slab = create_perturbation(configuration=perturbation_config, preserve_distance=True) + # Check selected atoms to avoid using 100+ atoms fixture + assertion_utils.assert_deep_almost_equal([0.0, 0.0, 0.5], perturbed_slab.basis.coordinates.values[0]) + assertion_utils.assert_deep_almost_equal( + [0.201132951, 0.1, 0.546942315], perturbed_slab.basis.coordinates.values[42] + ) + # Value taken from visually inspected notebook + expected_cell = Cell( + vector1=[23.517094, 0.0, 0.0], + vector2=[-11.758818, 21.367367, 0.0], + vector3=[0.0, 0.0, 20.0], + ) + assertion_utils.assert_deep_almost_equal(expected_cell.vectors_as_array, perturbed_slab.basis.cell.vectors_as_array)