diff --git a/src/py/mat3ra/made/cell.py b/src/py/mat3ra/made/cell.py index 01dd95581..15bff1a1a 100644 --- a/src/py/mat3ra/made/cell.py +++ b/src/py/mat3ra/made/cell.py @@ -51,3 +51,7 @@ def convert_point_to_crystal(self, point: List[float]) -> List[float]: 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() + + @property + def volume(self) -> float: + return np.linalg.det(np.array(self.vectors_as_array)) diff --git a/src/py/mat3ra/made/material.py b/src/py/mat3ra/made/material.py index 2b46f31c5..be0dafbba 100644 --- a/src/py/mat3ra/made/material.py +++ b/src/py/mat3ra/made/material.py @@ -109,7 +109,7 @@ def set_new_lattice_vectors( lattice = Lattice.from_vectors_array([lattice_vector1, lattice_vector2, lattice_vector3]) self.lattice = lattice - def add_atom(self, element: str, coordinate: List[float]) -> None: + def add_atom(self, element: str, coordinate: List[float], use_cartesian_coordinates=False) -> None: new_basis = self.basis.copy() - new_basis.add_atom(element, coordinate) + new_basis.add_atom(element, coordinate, use_cartesian_coordinates) self.basis = new_basis diff --git a/src/py/mat3ra/made/tools/analyze/__init__.py b/src/py/mat3ra/made/tools/analyze/__init__.py new file mode 100644 index 000000000..76cbb7740 --- /dev/null +++ b/src/py/mat3ra/made/tools/analyze/__init__.py @@ -0,0 +1,21 @@ +import numpy as np +from mat3ra.made.material import Material +from scipy.spatial.distance import pdist + + +class BaseMaterialAnalyzer: + def __init__(self, material: Material): + self.material = material.clone() + self.material.to_cartesian() + + @property + def volume(self): + return self.material.basis.cell.volume + + @property + def atomic_density(self): + return len(self.material.coordinates_array) / self.volume + + @property + def pairwise_distances(self): + return pdist(np.array(self.material.coordinates_array)) diff --git a/src/py/mat3ra/made/tools/analyze/coordination.py b/src/py/mat3ra/made/tools/analyze/coordination.py new file mode 100644 index 000000000..7814a1e62 --- /dev/null +++ b/src/py/mat3ra/made/tools/analyze/coordination.py @@ -0,0 +1,56 @@ +from typing import List, Optional + +from mat3ra.made.material import Material + +from ..convert import to_pymatgen +from ..third_party import PymatgenVoronoiNN + + +def get_voronoi_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. + """ + if coordinate is None: + coordinate = [0, 0, 0] + structure = to_pymatgen(material) + voronoi_nn = PymatgenVoronoiNN( + tol=tolerance, + cutoff=cutoff, + weight="solid_angle", + extra_nn_info=False, + compute_adj_neighbors=True, + ) + 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] + if remove_dummy_atom: + structure.remove_sites([-1]) + + all_coordinates = material.basis.coordinates + all_coordinates.filter_by_indices(neighboring_atoms_pymatgen_ids) + return all_coordinates.ids diff --git a/src/py/mat3ra/made/tools/analyze/material.py b/src/py/mat3ra/made/tools/analyze/material.py new file mode 100644 index 000000000..6cdddf3d3 --- /dev/null +++ b/src/py/mat3ra/made/tools/analyze/material.py @@ -0,0 +1,275 @@ +from typing import List, Optional, Tuple + +import numpy as np +from mat3ra.made.material import Material +from mat3ra.made.utils import ArrayWithIds +from scipy.spatial._ckdtree import cKDTree + +from ..bonds import BondDirections, BondDirectionsTemplatesForElement +from ..site import CrystalSite, CrystalSiteList +from .rdf import RadialDistributionFunction +from .utils import decorator_handle_periodic_boundary_conditions + + +class MaterialWithCrystalSites(Material): + crystal_sites: CrystalSiteList = CrystalSiteList(values=[]) + + @classmethod + def from_material(cls, material: Material): + new_material = material.clone() + new_material.to_cartesian() + config = new_material.to_json() + return cls(config) + + def analyze(self): + """ + Analyze the material in place and generate the crystal sites with properties. + + Generated properties: + - nearest_neighbor_vectors: The nearest neighbor vectors for all sites. + """ + self.nearest_neighbor_vectors = self.get_neighbors_vectors_for_all_sites(cutoff=3.0) + self.crystal_sites = CrystalSiteList( + values=[CrystalSite(nearest_neighbor_vectors=item) for item in self.nearest_neighbor_vectors.values], + ids=self.nearest_neighbor_vectors.ids, + ) + + @property + def coordinates_as_kdtree(self): + return cKDTree(self.basis.coordinates.values) + + def get_neighbors_vectors_for_site( + self, site_index: int, cutoff: float = 3.0, max_number_of_neighbors: Optional[int] = None + ) -> List[np.ndarray]: + """ + Get the nearest neighbor vectors for a specific site. + + Args: + site_index (int): The index of the site to + cutoff (float): The maximum cutoff radius for identifying neighbors. + max_number_of_neighbors (int): The max number of neighbors possible. + + Returns: + List[np.ndarray]: List of vectors to the nearest neighbors. + """ + coordinates = self.basis.coordinates.values + neighbors, distances = self.get_neighbors_for_site(site_index, cutoff, max_number_of_neighbors) + vectors = [np.array(coordinates[neighbor]) - np.array(coordinates[site_index]) for neighbor in neighbors] + return vectors + + @decorator_handle_periodic_boundary_conditions(cutoff=0.25) + def get_neighbors_vectors_for_all_sites( + self, cutoff: float = 3.0, max_number_of_neighbors: Optional[int] = None + ) -> ArrayWithIds: + """ + Get the nearest neighbor vectors for all sites. + + Args: + cutoff (float): The maximum cutoff radius for identifying neighbors. + max_number_of_neighbors (int): The max number of neighbors possible. + + Returns: + ArrayWithIds: Array with the nearest neighbor vectors for all sites. + """ + nearest_neighbors = ArrayWithIds() + for site_index in range(len(self.basis.coordinates.values)): + vectors = self.get_neighbors_vectors_for_site(site_index, cutoff, max_number_of_neighbors) + nearest_neighbors.add_item(vectors, site_index) + return nearest_neighbors + + def get_neighbors_for_site( + self, + site_index: int, + cutoff: float = 3.0, + max_number_of_neighbors: Optional[int] = None, + nearest_only: bool = True, + ) -> Tuple[np.ndarray, np.ndarray]: + """ + Get the nearest neighbors for a specific site. + + Args: + site_index (int): The index of the site to + cutoff (float): The maximum cutoff radius for identifying neighbors. + max_number_of_neighbors (int): The max number of neighbors possible. + nearest_only (bool): If True, only consider the first shell of neighbors. + + Returns: + Tuple[np.ndarray, np.ndarray]: Tuple of neighbors and distances. + """ + coordinates = self.basis.coordinates + max_number_of_neighbors = ( + len(coordinates.values) if max_number_of_neighbors is None else max_number_of_neighbors + ) + coordinate = coordinates.get_element_value_by_index(site_index) + distances, neighbors = self.coordinates_as_kdtree.query( + coordinate, k=max_number_of_neighbors, distance_upper_bound=cutoff + ) + valid_indices = (distances != np.inf) & (distances != 0) + distances = distances[valid_indices] + neighbors = neighbors[valid_indices] + + if nearest_only: + rdf = RadialDistributionFunction.from_material(self) + valid_indices = np.where([rdf.is_within_first_peak(distance) for distance in distances]) + distances = distances[valid_indices] + neighbors = neighbors[valid_indices] + return neighbors, distances + + def get_nearest_neighbors_for_all_sites( + self, cutoff: float = 3.0, max_number_of_neighbors: Optional[int] = None + ) -> ArrayWithIds: + """ + Get the nearest neighbors for all sites. + + Args: + cutoff (float): The maximum cutoff radius for identifying neighbors. + max_number_of_neighbors (int): The max number of neighbors possible. + + Returns: + ArrayWithIds: Array with the nearest neighbors for all sites. + """ + nearest_neighbors = ArrayWithIds() + for site_index in self.basis.coordinates.ids: + neighbors, distances = self.get_neighbors_for_site(site_index, cutoff, max_number_of_neighbors) + nearest_neighbors.add_item(neighbors, site_index) + return nearest_neighbors + + @decorator_handle_periodic_boundary_conditions(cutoff=0.25) + def get_coordination_numbers(self, cutoff: float = 3.0) -> ArrayWithIds: + """ + Calculate the coordination numbers for all atoms in the material. + + Args: + cutoff (float): The cutoff radius for identifying neighbors. + + Returns: + ArrayWithIds: Array with the coordination numbers for all sites. + """ + nearest_neighbors = self.get_nearest_neighbors_for_all_sites(cutoff) + coordination_numbers = [len(neighbors) for neighbors in nearest_neighbors.values] + return ArrayWithIds(values=coordination_numbers, ids=nearest_neighbors.ids) + + def get_undercoordinated_atom_indices(self, cutoff: float, coordination_threshold: int) -> List[int]: + """ + Identify undercoordinated atoms based on the coordination threshold (inclusive). + + Args: + cutoff (float): The cutoff radius for identifying neighbors. + coordination_threshold (int): The coordination number threshold for an atom + to be considered undercoordinated, inclusive. + + Returns: + List[int]: List of indices of undercoordinated atoms. + """ + coordination_numbers = self.get_coordination_numbers(cutoff) + return [ + item.id + for item in coordination_numbers.to_array_of_values_with_ids() + if item.value <= coordination_threshold + ] + + def get_unique_coordination_numbers(self, cutoff: float) -> List[int]: + """ + Get the unique coordination numbers for all atoms in the material. + + Args: + self (MaterialWithCrystalSites): The material object with crystal sites. + cutoff (float): The cutoff radius for identifying neighbors. + + Returns: + List[int]: A list of unique coordination numbers present in the material. + """ + coordination_numbers = self.get_coordination_numbers(cutoff) + return sorted(list(set(coordination_numbers.values))) + + @decorator_handle_periodic_boundary_conditions(cutoff=0.25) + def find_missing_bonds_for_all_sites( + self, + bond_templates_list: List[BondDirectionsTemplatesForElement], + max_bonds_to_add: int = 1, + angle_tolerance: float = 0.1, + ) -> List[BondDirections]: + """ + Find missing bonds for all sites in the material. + + Args: + self (MaterialWithCrystalSites): The material object with crystal sites. + bond_templates_list (List[BondDirectionsTemplatesForElement]): List of bond templates for each element. + + Returns: + List[BondDirections]: List of missing bonds for each site. + """ + missing_bonds: List[BondDirections] = [] + for coordinate in self.basis.coordinates.to_array_of_values_with_ids(): + existing_bond_directions = BondDirections(self.get_neighbors_vectors_for_site(coordinate.id, cutoff=3.0)) + element = self.basis.elements.get_element_value_by_index(coordinate.id) + bond_templates = [ + bond_templates_for_element.to_ndarray() + for bond_templates_for_element in bond_templates_list + if bond_templates_for_element.element == element + ] + if bond_templates is None: + continue + missing_bonds_for_site = existing_bond_directions.find_missing_directions( + bond_templates=bond_templates, + angle_tolerance=angle_tolerance, + max_bonds_to_add=max_bonds_to_add, + ) + missing_bonds.append(BondDirections(missing_bonds_for_site)) + + return missing_bonds + + def find_unique_bond_directions(self) -> List[BondDirectionsTemplatesForElement]: + """ + Find unique bond templates for each element type. + + Args: + self (Material): The material object. + + Returns: + List[BondDirectionsTemplatesForElement]: List of unique bond templates for each element. + + Example schematic structure for Graphene with 2 unique bond directions: + {"C": [[0, 1, 0], [-2, -1, 0], [2, -1, 0]], [[0,-1,0], [-2, 1, 0], [2, 1, 0]]} + """ + element_templates: List[BondDirectionsTemplatesForElement] = [] + + for element in set(self.basis.elements.values): + element_templates.append(self.find_bond_directions_for_element(element)) + + return element_templates + + def find_bond_directions_for_element(self, element: str) -> BondDirectionsTemplatesForElement: + """ + Find bond directions templates for each element. + + Args: + self (MaterialWithCrystalSites): The material object with crystal sites. + element (str): Chemical element to find templates for. + + Returns: + List[BondDirectionsTemplatesForElement]: List of unique bond templates for the element. + """ + atom_vectors = self.nearest_neighbor_vectors + atom_elements = self.basis.elements.values + + element_indices = [i for i, e in enumerate(atom_elements) if e == element] + element_vector_lists = [np.array(atom_vectors.get_element_value_by_index(i)) for i in element_indices] + + if not element_vector_lists: + return BondDirectionsTemplatesForElement(bond_directions_templates=[], element=element) + + max_coordination_number = max(len(vectors) for vectors in element_vector_lists) + max_coordination_number_vectors = [v for v in element_vector_lists if len(v) == max_coordination_number] + + unique_templates: BondDirectionsTemplatesForElement = BondDirectionsTemplatesForElement( + bond_directions_templates=[], element=element + ) + for template in max_coordination_number_vectors: + bond_direction = BondDirections(template) + if not any( + np.array_equal(bond_direction, existing) for existing in unique_templates.bond_directions_templates + ): + unique_templates.bond_directions_templates.append(bond_direction) + + return unique_templates diff --git a/src/py/mat3ra/made/tools/analyze.py b/src/py/mat3ra/made/tools/analyze/other.py similarity index 78% rename from src/py/mat3ra/made/tools/analyze.py rename to src/py/mat3ra/made/tools/analyze/other.py index 4c5962380..33c257898 100644 --- a/src/py/mat3ra/made/tools/analyze.py +++ b/src/py/mat3ra/made/tools/analyze/other.py @@ -4,10 +4,11 @@ from mat3ra.made.material import Material from scipy.spatial import cKDTree -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_convert_position_to_coordinate, decorator_handle_periodic_boundary_conditions +from ..convert import decorator_convert_material_args_kwargs_to_atoms, to_pymatgen +from ..enums import SurfaceTypes +from ..third_party import ASEAtoms, PymatgenIStructure +from ..utils import decorator_convert_position_to_coordinate +from .utils import decorator_handle_periodic_boundary_conditions @decorator_convert_material_args_kwargs_to_atoms @@ -278,56 +279,6 @@ def get_atom_indices_with_condition_on_coordinates( return selected_indices -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. - """ - if coordinate is None: - coordinate = [0, 0, 0] - structure = to_pymatgen(material) - voronoi_nn = PymatgenVoronoiNN( - tol=tolerance, - cutoff=cutoff, - weight="solid_angle", - extra_nn_info=False, - compute_adj_neighbors=True, - ) - 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] - if remove_dummy_atom: - 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", @@ -393,7 +344,7 @@ def is_shadowed_by_neighbors_from_surface( ) -@decorator_handle_periodic_boundary_conditions(cutoff=0.1) +@decorator_handle_periodic_boundary_conditions(cutoff=0.25) def get_surface_atom_indices( material: Material, surface: SurfaceTypes = SurfaceTypes.TOP, shadowing_radius: float = 2.5, depth: float = 5 ) -> List[int]: @@ -427,63 +378,6 @@ def get_surface_atom_indices( 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 - - @decorator_convert_position_to_coordinate def get_local_extremum_atom_index( material: Material, diff --git a/src/py/mat3ra/made/tools/analyze/rdf.py b/src/py/mat3ra/made/tools/analyze/rdf.py new file mode 100644 index 000000000..3c676935d --- /dev/null +++ b/src/py/mat3ra/made/tools/analyze/rdf.py @@ -0,0 +1,74 @@ +import numpy as np +from mat3ra.made.material import Material +from pydantic import BaseModel + +from . import BaseMaterialAnalyzer + + +class RadialDistributionFunction(BaseModel): + rdf: np.ndarray + bin_centers: np.ndarray + + class Config: + arbitrary_types_allowed = True + + @classmethod + def from_material(cls, material: Material, cutoff: float = 10.0, bin_size: float = 0.1): + analyzer = BaseMaterialAnalyzer(material) + distances = analyzer.pairwise_distances + density = analyzer.atomic_density + distances = distances[distances <= cutoff] + + # Bin distances into a histogram + bins = np.arange(0, cutoff + bin_size, bin_size) # Bin edges + hist, bin_edges = np.histogram(distances, bins=bins, density=False) + + # Convert to radial distribution function + bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + shell_volumes = ( + (4 / 3) * np.pi * (np.power(bin_edges[1:], 3) - np.power(bin_edges[:-1], 3)) + ) # Volume of spherical shells + + rdf = hist / (shell_volumes * density) + + return cls(rdf=rdf, bin_centers=bin_centers) + + @property + def first_peak_index(self): + return np.argmax(self.rdf[1:]) + 1 + + @property + def first_peak_value(self): + return self.rdf[self.first_peak_index] + + @property + def first_peak_width(self): + half_max = 0.5 * self.first_peak_value + + # Find left boundary + left_index = self.first_peak_index + while left_index > 0 and self.rdf[left_index] > half_max: + left_index -= 1 + + # Find right boundary + right_index = self.first_peak_index + while right_index < len(self.rdf) - 1 and self.rdf[right_index] > half_max: + right_index += 1 + + # Compute the width + left_boundary = self.bin_centers[left_index] + right_boundary = self.bin_centers[right_index] + first_peak_width = right_boundary - left_boundary + + return float(first_peak_width) + + @property + def first_peak_distance(self): + return self.bin_centers[self.first_peak_index] + + def is_within_first_peak(self, distance: float, tolerance: float = 0.1) -> bool: + return ( + self.first_peak_distance - 0.5 * self.first_peak_width - tolerance + < distance + < self.first_peak_distance + 0.5 * self.first_peak_width + tolerance + ) diff --git a/src/py/mat3ra/made/tools/analyze/utils.py b/src/py/mat3ra/made/tools/analyze/utils.py new file mode 100644 index 000000000..09a7abee0 --- /dev/null +++ b/src/py/mat3ra/made/tools/analyze/utils.py @@ -0,0 +1,87 @@ +from functools import wraps +from typing import TYPE_CHECKING, Union + +import numpy as np +from mat3ra.made.material import Material +from mat3ra.made.utils import ArrayWithIds + +if TYPE_CHECKING: + from .material import MaterialWithCrystalSites + + +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): + original_basis_is_in_cartesian = material.basis.is_in_cartesian_units + if original_basis_is_in_cartesian: + material.to_crystal() + augmented_material, last_id = augment_material_with_periodic_images(material, cutoff) + if original_basis_is_in_cartesian: + augmented_material.to_cartesian() + result = func(augmented_material, *args, **kwargs) + if original_basis_is_in_cartesian: + material.to_cartesian() + original_ids = material.basis.coordinates.ids + if isinstance(result, ArrayWithIds): + result.filter_by_ids(original_ids) + 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: Union[Material, "MaterialWithCrystalSites"], cutoff: float = 1.0 # type: ignore +): + """ + 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/bonds.py b/src/py/mat3ra/made/tools/bonds.py new file mode 100644 index 000000000..cad7cad3e --- /dev/null +++ b/src/py/mat3ra/made/tools/bonds.py @@ -0,0 +1,173 @@ +from enum import Enum +from typing import List + +import numpy as np +from pydantic import BaseModel, ConfigDict + + +class BondDirectionsTemplatesEnum(list, Enum): + OCTAHEDRAL = [[1, 0, 0], [-1, 0, 0], [0, 1, 0], [0, -1, 0], [0, 0, 1], [0, 0, -1]] + TETRAHEDRAL = [[1, 1, 1], [-1, -1, 1], [1, -1, -1], [-1, 1, -1]] + PLANAR = [[1, 0, 0], [-1, 2, 0], [-1, -2, 0]] + LINEAR = [[1, 0, 0], [-1, 0, 0]] + + +class BondDirections(np.ndarray): + def __new__(cls, input_array, *args, **kwargs): + """ + Create a new instance of BondDirections from an input array. + + Args: + input_array (array-like): Input data to initialize the array. + + Returns: + BondDirections: A new BondDirections instance. + """ + obj = np.asarray(input_array).view(cls) + return obj + + def __eq__(self, other): + if not isinstance(other, BondDirections): + return False + if self.shape != other.shape: + return False # Prevents shape mismatches + return self.are_bond_directions_similar(self, other) + + @staticmethod + def are_bond_directions_similar(directions1: np.ndarray, directions2: np.ndarray, tolerance: float = 0.1) -> bool: + """ + Check if two bond direction templates are similar. + + Args: + directions1 (np.ndarray): First template of bond vectors. + directions2 (np.ndarray): Second template of bond vectors. + tolerance (float): Angle tolerance for comparison. + + Returns: + bool: True if the templates are similar, False otherwise. + """ + if len(directions1) != len(directions2): + return False + + dot_matrix = np.dot(directions1, directions2.T) + norms1 = np.linalg.norm(directions1, axis=1) + norms2 = np.linalg.norm(directions2, axis=1) + cosine_matrix = dot_matrix / np.outer(norms1, norms2) + angles_matrix = np.arccos(np.clip(cosine_matrix, -1.0, 1.0)) + + unmatched = list(range(len(directions2))) + for angle_row in angles_matrix: + matches = np.where(angle_row < tolerance)[0] + if len(matches) == 0: + return False + unmatched.remove(matches.tolist()[0]) + + return True + + def find_missing_directions( + self, + bond_templates: List[np.ndarray], + angle_tolerance: float = 0.1, + max_bonds_to_add: int = 1, + ) -> "BondDirections": + """ + Reconstruct missing bonds for the atom based on the bond template that matches most of the existing bonds. + + Args: + bond_templates (List[np.ndarray]): List of bond templates to match against. + angle_tolerance (float): Tolerance for comparing angles between bond vectors. + max_bonds_to_add (int): Maximum number of bonds to add. + + Returns: + BondDirections: List of reconstructed bond vectors. + """ + max_coordination_number = max(len(template) for template in bond_templates) + + if len(self) >= max_coordination_number: + return BondDirections([]) + + best_missing_directions = None + highest_match_count = -1 + + bond_templates = self._flatten_bond_templates(bond_templates) + + for bond_template in bond_templates: + # Check which existing bonds match the template bonds + is_matched = [ + any(self._are_vectors_similar(existing_bond, candidate_bond, angle_tolerance) for existing_bond in self) + for candidate_bond in bond_template + ] + + current_match_count = sum(is_matched) + missing_directions = [ + candidate_bond for candidate_bond, matched in zip(bond_template, is_matched) if not matched + ] + + # Select the template with the highest match count + if current_match_count > highest_match_count: + highest_match_count = current_match_count + best_missing_directions = missing_directions + + if best_missing_directions is not None: + num_bonds_to_add = min( + len(best_missing_directions), + max_bonds_to_add, + max_coordination_number - len(self), + ) + return BondDirections(best_missing_directions[:num_bonds_to_add]) + + return BondDirections([]) + + @staticmethod + def _are_vectors_similar(vector1: np.ndarray, vector2: np.ndarray, tolerance: float = 0.1) -> bool: + """ + Check if two bond vectors are similar. + + Args: + vector1 (np.ndarray): First bond vector. + vector2 (np.ndarray): Second bond vector. + tolerance (float): Angle tolerance for comparison. + + Returns: + bool: True if the bond vectors are similar, False otherwise. + """ + dot_product = np.dot(vector1, vector2) + norms = np.linalg.norm(vector1) * np.linalg.norm(vector2) + cosine_angle = dot_product / norms if norms != 0 else 0.0 + angle = np.arccos(np.clip(cosine_angle, -1.0, 1.0)) + return angle < tolerance + + @staticmethod + def _flatten_bond_templates(bond_templates: List[np.ndarray]) -> List[np.ndarray]: + """ + Ensure that bond_templates is a list of np.ndarray with no extra nesting. + + Args: + bond_templates (list): List of bond templates to preprocess. + + Returns: + list: Flattened list of bond templates (np.ndarray). + """ + flattened_templates = [] + for template in bond_templates: + if isinstance(template, list) and all(isinstance(sub_template, np.ndarray) for sub_template in template): + # If the template is a list of np.ndarray, extend the flattened list + flattened_templates.extend(template) + else: + # Otherwise, add the template directly + flattened_templates.append(template) + return flattened_templates + + +class BondDirectionsTemplatesForElement(BaseModel): + bond_directions_templates: List[BondDirections] + element: str + + model_config = ConfigDict(arbitrary_types_allowed=True) + + def to_ndarray(self): + return [np.array(template) for template in self.bond_directions_templates] + + +class BondDirectionsForElementList(List[BondDirectionsTemplatesForElement]): + pass diff --git a/src/py/mat3ra/made/tools/build/defect/builders.py b/src/py/mat3ra/made/tools/build/defect/builders.py index 24b738012..d9a6524f7 100644 --- a/src/py/mat3ra/made/tools/build/defect/builders.py +++ b/src/py/mat3ra/made/tools/build/defect/builders.py @@ -24,13 +24,13 @@ ) from ...build import BaseBuilder from ...convert import to_pymatgen -from ...analyze import ( - get_nearest_neighbors_atom_indices, +from ...analyze.other import ( get_atomic_coordinates_extremum, get_closest_site_id_from_coordinate, get_closest_site_id_from_coordinate_and_element, get_local_extremum_atom_index, ) +from ...analyze.coordination import get_voronoi_nearest_neighbors_atom_indices from ...utils import transform_coordinate_to_supercell, coordinate as CoordinateCondition from ..utils import merge_materials from ..slab import SlabConfiguration, create_slab, Termination @@ -346,7 +346,7 @@ def get_equidistant_position( coordinate=adatom_coordinate, scaling_factor=scaling_factor, translation_vector=translation_vector ) - neighboring_atoms_ids_in_supercell = get_nearest_neighbors_atom_indices( + neighboring_atoms_ids_in_supercell = get_voronoi_nearest_neighbors_atom_indices( material=supercell_material, coordinate=adatom_coordinate_in_supercell ) if neighboring_atoms_ids_in_supercell is None: diff --git a/src/py/mat3ra/made/tools/build/defect/configuration.py b/src/py/mat3ra/made/tools/build/defect/configuration.py index 7eda27cdc..b4a297733 100644 --- a/src/py/mat3ra/made/tools/build/defect/configuration.py +++ b/src/py/mat3ra/made/tools/build/defect/configuration.py @@ -4,7 +4,7 @@ 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 ...analyze.other import get_closest_site_id_from_coordinate, get_atomic_coordinates_extremum from ...utils.coordinate import ( CylinderCoordinateCondition, SphereCoordinateCondition, diff --git a/src/py/mat3ra/made/tools/build/grain_boundary/builders.py b/src/py/mat3ra/made/tools/build/grain_boundary/builders.py index d3621bf50..94ea6d35a 100644 --- a/src/py/mat3ra/made/tools/build/grain_boundary/builders.py +++ b/src/py/mat3ra/made/tools/build/grain_boundary/builders.py @@ -4,7 +4,7 @@ from mat3ra.made.material import Material from ...third_party import PymatgenInterface -from ...analyze import get_chemical_formula +from ...analyze.other import get_chemical_formula from ..slab import SlabConfiguration, get_terminations, create_slab from ..interface import ZSLStrainMatchingInterfaceBuilderParameters, InterfaceConfiguration diff --git a/src/py/mat3ra/made/tools/build/interface/builders.py b/src/py/mat3ra/made/tools/build/interface/builders.py index 2c72c1c6a..ddb9461a3 100644 --- a/src/py/mat3ra/made/tools/build/interface/builders.py +++ b/src/py/mat3ra/made/tools/build/interface/builders.py @@ -17,7 +17,7 @@ add_vacuum_sides, add_vacuum, ) -from ...analyze import get_chemical_formula +from ...analyze.other import get_chemical_formula from ...convert import to_ase, from_ase, to_pymatgen, PymatgenInterface, ASEAtoms from ...build import BaseBuilder, BaseBuilderParameters from ..nanoribbon import NanoribbonConfiguration, create_nanoribbon diff --git a/src/py/mat3ra/made/tools/build/passivation/__init__.py b/src/py/mat3ra/made/tools/build/passivation/__init__.py index c26e63b04..9d467b1e4 100644 --- a/src/py/mat3ra/made/tools/build/passivation/__init__.py +++ b/src/py/mat3ra/made/tools/build/passivation/__init__.py @@ -1,4 +1,4 @@ -from typing import Union +from typing import Union, List from mat3ra.made.material import Material from .configuration import PassivationConfiguration @@ -6,7 +6,9 @@ SurfacePassivationBuilder, CoordinationBasedPassivationBuilder, SurfacePassivationBuilderParameters, + CoordinationBasedPassivationBuilderParameters, ) +from ...analyze.material import MaterialWithCrystalSites def create_passivation( @@ -18,8 +20,19 @@ def create_passivation( return builder.get_material(configuration) -def get_unique_coordination_numbers(configuration: PassivationConfiguration) -> set: +def get_unique_coordination_numbers( + configuration: PassivationConfiguration, + cutoff: float = 3.0, +) -> List[int]: """ - Get the unique coordination numbers for the provided configuration as a set type. + Get the unique coordination numbers for the provided passivation configuration and cutoff radius. + + Args: + configuration (PassivationConfiguration): The configuration object. + cutoff (float): The cutoff radius for defining neighbors. + Returns: + set: The unique coordination numbers. """ - return CoordinationBasedPassivationBuilder().get_unique_coordination_numbers(material=configuration.slab) + material_with_crystal_sites = MaterialWithCrystalSites.from_material(configuration.slab) + material_with_crystal_sites.analyze() + return material_with_crystal_sites.get_unique_coordination_numbers(cutoff=cutoff) diff --git a/src/py/mat3ra/made/tools/build/passivation/builders.py b/src/py/mat3ra/made/tools/build/passivation/builders.py index 19aa6ff27..a5f61a766 100644 --- a/src/py/mat3ra/made/tools/build/passivation/builders.py +++ b/src/py/mat3ra/made/tools/build/passivation/builders.py @@ -1,21 +1,15 @@ from typing import List - -import numpy as np from mat3ra.made.material import Material +from mat3ra.made.tools.bonds import BondDirections from pydantic import BaseModel +import numpy as np +from ...analyze.material import MaterialWithCrystalSites from ...enums import SurfaceTypes -from ...analyze import ( - get_surface_atom_indices, - get_undercoordinated_atom_indices, - get_nearest_neighbors_atom_indices, - get_coordination_numbers, -) +from ...analyze.other import get_surface_atom_indices from ...modify import translate_to_z_level from ...build import BaseBuilder -from .configuration import ( - PassivationConfiguration, -) +from .configuration import PassivationConfiguration class PassivationBuilder(BaseBuilder): @@ -40,7 +34,9 @@ def create_passivated_material(self, configuration: BaseBuilder._ConfigurationTy material = translate_to_z_level(configuration.slab, "center") return material - def _add_passivant_atoms(self, material: Material, coordinates: list, passivant: str) -> Material: + def _add_passivant_atoms( + self, material: Material, coordinates: list, passivant: str, use_cartesian_coordinates=False + ) -> Material: """ Add passivant atoms to the provided coordinates in the material. @@ -48,12 +44,13 @@ def _add_passivant_atoms(self, material: Material, coordinates: list, passivant: 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'). + use_cartesian_coordinates (bool): Whether the coordinates are in Cartesian units (or crystal by default). Returns: Material: The material object with passivation atoms added. """ for coord in coordinates: - material.add_atom(passivant, coord) + material.add_atom(passivant, coord, use_cartesian_coordinates) return material @@ -127,88 +124,98 @@ def _get_passivant_coordinates( class CoordinationBasedPassivationBuilderParameters(SurfacePassivationBuilderParameters): """ - Parameters for the CoordinationPassivationBuilder. + Parameters for the CoordinationPassivationBuilder. + Args: coordination_threshold (int): The coordination number threshold for atom to be considered undercoordinated. + bonds_to_passivate (int): The maximum number of bonds to passivate for each undercoordinated atom. + symmetry_tolerance (float): The tolerance for symmetry comparison of vectors for bonds. """ - coordination_threshold: int = 3 + coordination_threshold: int = 3 # The coordination number threshold for an atom to be considered undercoordinated. + bonds_to_passivate: int = 1 # The maximum number of bonds to passivate for each undercoordinated atom. + symmetry_tolerance: float = 0.1 # The tolerance for symmetry comparison of vectors for bonds. 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. + Builder for passivating material surfaces based on coordination number analysis. + + This builder analyzes atomic coordination environments and reconstructs missing bonds + using templates derived from fully coordinated atoms. It works by: + + 1. Finding undercoordinated atoms that need passivation + 2. Creating bond vector templates for each element type by: + - Collecting vectors from atoms with the highest valid coordination + - Grouping similar vector arrangements into unique templates + 3. Reconstructing missing bonds by: + - Matching existing bonds against templates + - Finding the best-matching template for each atom + - Adding missing bond vectors from the template + 4. Placing passivant atoms (typically H) along the reconstructed bond vectors """ - _BuildParametersType = CoordinationBasedPassivationBuilderParameters + build_parameters: CoordinationBasedPassivationBuilderParameters _DefaultBuildParameters = CoordinationBasedPassivationBuilderParameters() + _ConfigurationType = PassivationConfiguration def create_passivated_material(self, configuration: PassivationConfiguration) -> Material: + """ + Create a passivated material based on the configuration. + """ 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, + material_with_crystal_sites = MaterialWithCrystalSites.from_material(material) + material_with_crystal_sites.analyze() + + undercoordinated_atoms_indices = material_with_crystal_sites.get_undercoordinated_atom_indices( cutoff=self.build_parameters.shadowing_radius, coordination_threshold=self.build_parameters.coordination_threshold, ) + # TODO: bonds_templates will be passed from the configuration in the "controlled" version of this class + bonds_templates = material_with_crystal_sites.find_unique_bond_directions() + reconstructed_bonds = material_with_crystal_sites.find_missing_bonds_for_all_sites( + bond_templates_list=bonds_templates, + max_bonds_to_add=self.build_parameters.bonds_to_passivate, + angle_tolerance=self.build_parameters.symmetry_tolerance, + ) passivant_coordinates_values = self._get_passivant_coordinates( - material, configuration, undercoordinated_atoms_indices + material, + configuration, + undercoordinated_atoms_indices, + reconstructed_bonds, ) - return self._add_passivant_atoms(material, passivant_coordinates_values, configuration.passivant) + return self._add_passivant_atoms(material, passivant_coordinates_values, configuration.passivant, False) def _get_passivant_coordinates( - self, material: Material, configuration: PassivationConfiguration, undercoordinated_atoms_indices: list + self, + material: Material, + configuration: PassivationConfiguration, + undercoordinated_atoms_indices: list, + reconstructed_bonds: List[BondDirections], ): """ - Calculate the coordinates for placing passivating atoms based on the specified edge type. + Calculate the coordinates for placing passivating atoms based on reconstructed bonds. Args: material (Material): Material to passivate. - configuration (SurfacePassivationConfiguration): Configuration for passivation. + configuration (PassivationConfiguration): Configuration for passivation. undercoordinated_atoms_indices (list): Indices of undercoordinated atoms. + reconstructed_bonds (List[BondDirections]): Reconstructed bond directions. + + Returns: + list: Coordinates where passivants should be added. """ 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 - ) + for idx in undercoordinated_atoms_indices: + for bond_vector in reconstructed_bonds[idx]: + bond_vector_np = np.array(bond_vector) + if np.linalg.norm(bond_vector_np) == 0: + continue # Avoid division by zero + normalized_bond = bond_vector_np / np.linalg.norm(bond_vector_np) * configuration.bond_length + normalized_bond_crystal = material.basis.cell.convert_point_to_crystal(normalized_bond) + passivant_coordinates.append( + material.basis.coordinates.get_element_value_by_index(idx) + normalized_bond_crystal + ) return passivant_coordinates - - def get_unique_coordination_numbers(self, material: Material) -> set: - """ - 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/slab/builders.py b/src/py/mat3ra/made/tools/build/slab/builders.py index 378c80918..282ec27e8 100644 --- a/src/py/mat3ra/made/tools/build/slab/builders.py +++ b/src/py/mat3ra/made/tools/build/slab/builders.py @@ -5,7 +5,7 @@ from ...modify import add_vacuum from ...third_party import PymatgenSlab, PymatgenSlabGenerator, label_pymatgen_slab_termination -from ...analyze import get_chemical_formula +from ...analyze.other import get_chemical_formula from ...convert import to_pymatgen from ...build import BaseBuilder from ...build.mixins import ConvertGeneratedItemsPymatgenStructureMixin diff --git a/src/py/mat3ra/made/tools/calculate/__init__.py b/src/py/mat3ra/made/tools/calculate/__init__.py index 88b728f53..b1686439f 100644 --- a/src/py/mat3ra/made/tools/calculate/__init__.py +++ b/src/py/mat3ra/made/tools/calculate/__init__.py @@ -1,7 +1,7 @@ from typing import List, Optional, Union from ...material import Material -from ..analyze import get_surface_area +from ..analyze.other import get_surface_area from ..build.interface.utils import get_slab from ..convert import decorator_convert_material_args_kwargs_to_atoms, from_ase from ..third_party import ASEAtoms, ASECalculator, ASECalculatorEMT diff --git a/src/py/mat3ra/made/tools/calculate/calculators.py b/src/py/mat3ra/made/tools/calculate/calculators.py index b80749766..67fe2b123 100644 --- a/src/py/mat3ra/made/tools/calculate/calculators.py +++ b/src/py/mat3ra/made/tools/calculate/calculators.py @@ -4,7 +4,7 @@ from mat3ra.made.material import Material from pydantic import BaseModel -from ..analyze import get_surface_atom_indices +from ..analyze.other import get_surface_atom_indices from ..convert.utils import InterfacePartsEnum from ..enums import SurfaceTypes from ..modify import interface_get_part diff --git a/src/py/mat3ra/made/tools/material.py b/src/py/mat3ra/made/tools/material.py deleted file mode 100644 index cd5376c39..000000000 --- a/src/py/mat3ra/made/tools/material.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -from mat3ra.code.constants import AtomicCoordinateUnits - - -def scale_one_lattice_vector(material, key="a", factor=1.0): - """ - Scales one lattice vector for the given material in place. - - Args: - material: The material to scale. - key: The key of the lattice vector to scale. - factor: The factor to scale the lattice vector by - """ - material.to_cartesian() - - lattice = getattr(material, "lattice") - setattr(lattice, key, getattr(lattice, key) * factor) - - setattr(material, "lattice", lattice) - - material.to_crystal() - - -def get_basis_config_translated_to_center(material): - """ - Updates the basis of a material by translating the coordinates - so that the center of the material and lattice are aligned. - Args: - material: The material to update. - """ - original_units = material.basis.units - material.to_cartesian() - updated_basis = material.basis - center_of_coordinates = updated_basis.center_of_coordinates_point() - center_of_lattice = 0.5 * np.sum(material.lattice.vector_arrays, axis=0) - translation_vector = center_of_lattice - center_of_coordinates - updated_basis.translate_by_vector(translation_vector) - material.set_basis(updated_basis.to_json()) - if original_units != AtomicCoordinateUnits.cartesian: - material.to_crystal() diff --git a/src/py/mat3ra/made/tools/modify.py b/src/py/mat3ra/made/tools/modify.py index ee39a8ac4..e128aab1d 100644 --- a/src/py/mat3ra/made/tools/modify.py +++ b/src/py/mat3ra/made/tools/modify.py @@ -3,7 +3,7 @@ import numpy as np from mat3ra.made.material import Material -from .analyze import ( +from .analyze.other import ( get_atom_indices_with_condition_on_coordinates, get_atom_indices_within_radius_pbc, get_atomic_coordinates_extremum, diff --git a/src/py/mat3ra/made/tools/site.py b/src/py/mat3ra/made/tools/site.py new file mode 100644 index 000000000..36a0d93e9 --- /dev/null +++ b/src/py/mat3ra/made/tools/site.py @@ -0,0 +1,25 @@ +from typing import List, Optional + +import numpy as np +from mat3ra.made.utils import ArrayWithIds +from pydantic import BaseModel + + +class CrystalSite(BaseModel): + # element: str + # coordinate: List[float] + nearest_neighbor_vectors: List[np.ndarray] = [] + # coordination_number: int = 0 + # see https://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-wp-list for an example + wyckoff_letter: Optional[str] = None + + class Config: + arbitrary_types_allowed = True + + @property + def coordination_number(self): + return len(self.nearest_neighbor_vectors) + + +class CrystalSiteList(ArrayWithIds): + values: List[CrystalSite] diff --git a/src/py/mat3ra/made/tools/utils/__init__.py b/src/py/mat3ra/made/tools/utils/__init__.py index c3d67f283..7240dab73 100644 --- a/src/py/mat3ra/made/tools/utils/__init__.py +++ b/src/py/mat3ra/made/tools/utils/__init__.py @@ -1,8 +1,7 @@ from functools import wraps -from typing import Any, Callable, List, Optional, Tuple +from typing import Callable, List, Optional import numpy as np -from mat3ra.made.material import Material from mat3ra.utils.matrix import convert_2x2_to_3x3 from ..third_party import PymatgenStructure @@ -128,70 +127,3 @@ 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 f4533dc28..56faec13b 100644 --- a/src/py/mat3ra/made/tools/utils/coordinate.py +++ b/src/py/mat3ra/made/tools/utils/coordinate.py @@ -160,6 +160,9 @@ def is_coordinate_behind_plane( class CoordinateCondition(BaseModel): + class Config: + arbitrary_types_allowed = True + def condition(self, coordinate: List[float]) -> bool: raise NotImplementedError diff --git a/tests/py/unit/fixtures.py b/tests/py/unit/fixtures.py index 34883958c..aa7ba5b37 100644 --- a/tests/py/unit/fixtures.py +++ b/tests/py/unit/fixtures.py @@ -492,6 +492,99 @@ "isUpdated": True, } +GRAPHENE_ZIGZAG_NANORIBBON_PASSIVATED = { + "name": "Graphene (Zigzag nanoribbon) H-passivated", + "basis": { + "elements": [ + {"id": 0, "value": "C"}, + {"id": 1, "value": "C"}, + {"id": 2, "value": "C"}, + {"id": 3, "value": "C"}, + {"id": 4, "value": "C"}, + {"id": 5, "value": "C"}, + {"id": 6, "value": "C"}, + {"id": 7, "value": "C"}, + {"id": 8, "value": "C"}, + {"id": 9, "value": "C"}, + {"id": 10, "value": "C"}, + {"id": 11, "value": "C"}, + {"id": 12, "value": "C"}, + {"id": 13, "value": "C"}, + {"id": 14, "value": "C"}, + {"id": 15, "value": "C"}, + {"id": 16, "value": "H"}, + {"id": 17, "value": "H"}, + {"id": 18, "value": "H"}, + {"id": 19, "value": "H"}, + {"id": 20, "value": "H"}, + {"id": 21, "value": "H"}, + {"id": 22, "value": "H"}, + {"id": 23, "value": "H"}, + ], + "coordinates": [ + {"id": 0, "value": [0.062499937, 0.366666681, 0.5]}, + {"id": 1, "value": [0.312499937, 0.366666681, 0.5]}, + {"id": 2, "value": [0.187500062, 0.433333286, 0.5]}, + {"id": 3, "value": [0.187499937, 0.566666695, 0.5]}, + {"id": 4, "value": [0.062500062, 0.6333333, 0.5]}, + {"id": 5, "value": [0.562499937, 0.366666681, 0.5]}, + {"id": 6, "value": [0.437500062, 0.433333286, 0.5]}, + {"id": 7, "value": [0.437499937, 0.566666695, 0.5]}, + {"id": 8, "value": [0.312500062, 0.6333333, 0.5]}, + {"id": 9, "value": [0.812499938, 0.366666681, 0.5]}, + {"id": 10, "value": [0.687500062, 0.433333286, 0.5]}, + {"id": 11, "value": [0.687499937, 0.566666695, 0.5]}, + {"id": 12, "value": [0.562500062, 0.6333333, 0.5]}, + {"id": 13, "value": [0.937500063, 0.433333286, 0.5]}, + {"id": 14, "value": [0.937499937, 0.566666695, 0.5]}, + {"id": 15, "value": [0.812500063, 0.6333333, 0.5]}, + {"id": 16, "value": [0.062500067, 0.228137674, 0.5]}, + {"id": 17, "value": [0.312500067, 0.228137674, 0.5]}, + {"id": 18, "value": [0.062499932, 0.771862307, 0.5]}, + {"id": 19, "value": [0.562500067, 0.228137674, 0.5]}, + {"id": 20, "value": [0.312499932, 0.771862307, 0.5]}, + {"id": 21, "value": [0.812500068, 0.228137674, 0.5]}, + {"id": 22, "value": [0.562499932, 0.771862307, 0.5]}, + {"id": 23, "value": [0.812499933, 0.771862307, 0.5]}, + ], + "units": "crystal", + "cell": [[9.869164, 0.0, 0.0], [-0.0, 10.683683, 0.0], [0.0, 0.0, 20.0]], + "labels": [], + }, + "lattice": { + "a": 9.869164, + "b": 10.683683422, + "c": 20.0, + "alpha": 90.0, + "beta": 90.0, + "gamma": 90.0, + "units": {"length": "angstrom", "angle": "degree"}, + "type": "TRI", + "vectors": { + "a": [9.869164, 0.0, 0.0], + "b": [-0.0, 10.683683422, 0.0], + "c": [0.0, 0.0, 20.0], + "alat": 1, + "units": "angstrom", + }, + }, + "isNonPeriodic": False, + "_id": "", + "metadata": { + "boundaryConditions": {"type": "pbc", "offset": 0}, + "build": { + "configuration": { + "type": "PassivationConfiguration", + "slab": GRAPHENE_ZIGZAG_NANORIBBON, + "passivant": "H", + "bond_length": 1.48, + "surface": "both", + } + }, + }, + "isUpdated": True, +} + GRAPHENE_NICKEL_INTERFACE = { "name": "C2(001)-Ni4(111), Interface, Strain 0.105pct", diff --git a/tests/py/unit/test_tools_analyze.py b/tests/py/unit/test_tools_analyze.py index 86d366668..c795301d8 100644 --- a/tests/py/unit/test_tools_analyze.py +++ b/tests/py/unit/test_tools_analyze.py @@ -1,8 +1,10 @@ import numpy as np from ase.build import bulk -from mat3ra.made.tools.analyze import get_average_interlayer_distance, get_surface_area +from mat3ra.made.material import Material +from mat3ra.made.tools.analyze.other import get_average_interlayer_distance, get_surface_area +from mat3ra.made.tools.analyze.rdf import RadialDistributionFunction -from .fixtures import INTERFACE_ATOMS +from .fixtures import GRAPHENE_ZIGZAG_NANORIBBON, INTERFACE_ATOMS def test_calculate_average_interlayer_distance(): @@ -14,3 +16,29 @@ def test_calculate_surface_area(): atoms = bulk("Si", cubic=False) area = get_surface_area(atoms) assert np.isclose(area, 12.7673) + + +def test_radial_distribution_function(): + material = Material(GRAPHENE_ZIGZAG_NANORIBBON) + + rdf = RadialDistributionFunction.from_material(material, cutoff=10.0, bin_size=0.1) + + # Test that RDF is non-negative + assert np.all(rdf.rdf >= 0), "RDF contains negative values." + + # Test the first peak properties + assert rdf.first_peak_index > 0, "First peak index should be greater than 0." + assert rdf.first_peak_value > 0, "First peak value should be greater than 0." + assert rdf.first_peak_width > 0, "First peak width should be greater than 0." + assert rdf.first_peak_distance > 0, "First peak distance should be greater than 0." + + # Test that `is_within_first_peak` works as expected + assert rdf.is_within_first_peak(rdf.first_peak_distance), "First peak distance should be within the first peak." + assert not rdf.is_within_first_peak(0), "Zero distance should not be within the first peak." + assert not rdf.is_within_first_peak(10.0), "Cutoff distance should not be within the first peak." + + # Test that RDF drops to zero near the cutoff + assert np.isclose(rdf.rdf[-1], 0, atol=1e-2), "RDF should approach zero near the cutoff." + + # Specific material related tests + assert np.isclose(rdf.first_peak_distance, 1.42, atol=0.1), "First peak distance should be close to 1.42." diff --git a/tests/py/unit/test_tools_build_passivation.py b/tests/py/unit/test_tools_build_passivation.py index 469b73cc2..7724a2fa5 100644 --- a/tests/py/unit/test_tools_build_passivation.py +++ b/tests/py/unit/test_tools_build_passivation.py @@ -1,9 +1,15 @@ from mat3ra.made.material import Material -from mat3ra.made.tools.build.passivation.builders import SurfacePassivationBuilder, SurfacePassivationBuilderParameters +from mat3ra.made.tools.build.passivation import get_unique_coordination_numbers +from mat3ra.made.tools.build.passivation.builders import ( + CoordinationBasedPassivationBuilder, + CoordinationBasedPassivationBuilderParameters, + 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 +from .fixtures import GRAPHENE_ZIGZAG_NANORIBBON, GRAPHENE_ZIGZAG_NANORIBBON_PASSIVATED, SI_SLAB, SI_SLAB_PASSIVATED def test_passivate_surface(): @@ -13,3 +19,21 @@ def test_passivate_surface(): ) passivated_material = builder.get_material(config) assertion_utils.assert_deep_almost_equal(SI_SLAB_PASSIVATED, passivated_material.to_json()) + + +def test_get_unique_coordination_numbers(): + config = PassivationConfiguration( + slab=Material(GRAPHENE_ZIGZAG_NANORIBBON), passivant="H", bond_length=1.48, surface="both" + ) + unique_coordination_numbers = get_unique_coordination_numbers(config, cutoff=3.0) + assert unique_coordination_numbers == [2, 3] + + +def test_passivate_coordination_based(): + config = PassivationConfiguration(slab=Material(GRAPHENE_ZIGZAG_NANORIBBON), passivant="H", bond_length=1.48) + params = CoordinationBasedPassivationBuilderParameters( + shadowing_radius=2.5, coordination_threshold=2, bonds_to_passivate=1 + ) + builder = CoordinationBasedPassivationBuilder(build_parameters=params) + passivated_material = builder.get_material(config) + assertion_utils.assert_deep_almost_equal(GRAPHENE_ZIGZAG_NANORIBBON_PASSIVATED, passivated_material.to_json())