Merge pull request #158 from Exabyte-io/feature/SOF-7426

feature/SOF-7426 feat: displace film minimizing norm of distances

src/py/mat3ra/made/basis.py

@@ -1,5 +1,5 @@
 import json
-from typing import Dict, List, Optional
+from typing import Dict, List, Optional, Union
 
 from mat3ra.code.constants import AtomicCoordinateUnits
 from mat3ra.utils.mixins import RoundNumericValuesMixin
@@ -115,13 +115,24 @@ def add_atom(
         self.elements.add_item(element)
         self.coordinates.add_item(coordinate)
 
-    def remove_atom_by_id(self, id=None):
+    def remove_atom_by_id(self, id: int):
         self.elements.remove_item(id)
         self.coordinates.remove_item(id)
-        self.labels.remove_item(id)
+        if self.labels is not None:
+            self.labels.remove_item(id)
 
-    def filter_atoms_by_ids(self, ids):
+    def filter_atoms_by_ids(self, ids: Union[List[int], int]) -> "Basis":
         self.elements.filter_by_ids(ids)
         self.coordinates.filter_by_ids(ids)
         if self.labels is not None:
             self.labels.filter_by_ids(ids)
+        return self
+
+    def filter_atoms_by_labels(self, labels: Union[List[str], str]) -> "Basis":
+        if self.labels is None:
+            return self
+        self.labels.filter_by_values(labels)
+        ids = self.labels.ids
+        self.elements.filter_by_ids(ids)
+        self.coordinates.filter_by_ids(ids)
+        return self

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

@@ -1,9 +1,9 @@
 from typing import Callable, List, Literal, Optional
 
 import numpy as np
+from mat3ra.made.material import Material
 from scipy.spatial import cKDTree
 
-from ..material import Material
 from .convert import decorator_convert_material_args_kwargs_to_atoms, to_pymatgen
 from .enums import SurfaceTypes
 from .third_party import ASEAtoms, PymatgenIStructure, PymatgenVoronoiNN

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

@@ -1,6 +1,11 @@
-from typing import Union, List, Optional
+from typing import Union, List, Optional, Tuple
+
+import numpy as np
 
 from mat3ra.made.material import Material
+from ...calculate.calculators import InterfaceMaterialCalculator
+from ...modify import displace_interface_part
+from ...optimize import evaluate_calculator_on_xy_grid
 from .builders import (
     SimpleInterfaceBuilder,
     SimpleInterfaceBuilderParameters,
@@ -24,3 +29,52 @@ def create_interface(
     if builder is None:
         builder = SimpleInterfaceBuilder(build_parameters=SimpleInterfaceBuilderParameters())
     return builder.get_material(configuration)
+
+
+def get_optimal_film_displacement(
+    material: Material,
+    grid_size_xy: Tuple[int, int] = (10, 10),
+    grid_offset_position: List[float] = [0, 0],
+    grid_range_x=(-0.5, 0.5),
+    grid_range_y=(-0.5, 0.5),
+    use_cartesian_coordinates=False,
+    calculator: InterfaceMaterialCalculator = InterfaceMaterialCalculator(),
+):
+    """
+    Calculate the optimal displacement in of the film to minimize the interaction energy
+        between the film and the substrate. The displacement is calculated on a grid.
+
+    This function evaluates the interaction energy between the film and substrate
+    over a specified grid of (x,y) displacements. It returns the displacement vector that
+    results in the minimum interaction energy.
+
+    Args:
+        material (Material): The interface Material object.
+        grid_size_xy (Tuple[int, int]): The size of the grid to search for the optimal displacement.
+        grid_offset_position (List[float]): The offset position of the grid.
+        grid_range_x (Tuple[float, float]): The range of the grid in x.
+        grid_range_y (Tuple[float, float]): The range of the grid in y.
+        use_cartesian_coordinates (bool): Whether to use Cartesian coordinates.
+        calculator (InterfaceMaterialCalculator): The calculator to use for the calculation of the interaction energy.
+
+    Returns:
+        List[float]: The optimal displacement vector.
+
+    """
+    xy_matrix, results_matrix = evaluate_calculator_on_xy_grid(
+        material=material,
+        calculator_function=calculator.get_energy,
+        modifier=displace_interface_part,
+        modifier_parameters={},
+        grid_size_xy=grid_size_xy,
+        grid_offset_position=grid_offset_position,
+        grid_range_x=grid_range_x,
+        grid_range_y=grid_range_y,
+        use_cartesian_coordinates=use_cartesian_coordinates,
+    )
+    min_index = np.unravel_index(np.argmin(results_matrix), results_matrix.shape)
+
+    optimal_x = xy_matrix[0][min_index[0]]
+    optimal_y = xy_matrix[1][min_index[1]]
+
+    return [optimal_x, optimal_y, 0]

src/py/mat3ra/made/tools/calculate.py → src/py/mat3ra/made/tools/calculate/__init__.py

@@ -1,10 +1,11 @@
 from typing import Optional
 
-from ..material import Material
-from .analyze import get_surface_area
-from .build.interface.utils import get_slab
-from .convert import decorator_convert_material_args_kwargs_to_atoms
-from .third_party import ASEAtoms, ASECalculator, ASECalculatorEMT
+from ...material import Material
+from ..analyze import get_surface_area
+from ..build.interface.utils import get_slab
+from ..convert import decorator_convert_material_args_kwargs_to_atoms
+from ..third_party import ASEAtoms, ASECalculator, ASECalculatorEMT
+from .interaction_functions import sum_of_inverse_distances_squared
 
 
 @decorator_convert_material_args_kwargs_to_atoms

src/py/mat3ra/made/tools/calculate/calculators.py

@@ -0,0 +1,122 @@
+from typing import Callable
+
+import numpy as np
+from mat3ra.made.material import Material
+from pydantic import BaseModel
+
+from ..analyze import get_surface_atom_indices
+from ..convert.utils import InterfacePartsEnum
+from ..enums import SurfaceTypes
+from ..modify import get_interface_part
+from .interaction_functions import sum_of_inverse_distances_squared
+
+
+class MaterialCalculatorParameters(BaseModel):
+    """
+    Defines the parameters for a material calculator.
+
+    Args:
+        interaction_function (Callable): A function used to calculate the interaction metric between
+        sets of coordinates. The default function is sum_of_inverse_distances_squared.
+    """
+
+    interaction_function: Callable = sum_of_inverse_distances_squared
+
+
+class InterfaceMaterialCalculatorParameters(MaterialCalculatorParameters):
+    """
+    Parameters specific to the calculation of interaction energies between
+    an interface material's film and substrate.
+
+    Args:
+        shadowing_radius (float): Radius used to determine the surface atoms of the film or substrate
+                                  for interaction calculations. Default is 2.5 Å.
+    """
+
+    shadowing_radius: float = 2.5
+
+
+class MaterialCalculator(BaseModel):
+    """
+    A base class for performing calculations on materials.
+
+    This class uses the parameters defined in MaterialCalculatorParameters to calculate
+    interaction metrics between atoms or sets of coordinates within the material.
+
+    Args:
+        calculator_parameters (MaterialCalculatorParameters): Parameters controlling the calculator,
+        including the interaction function.
+    """
+
+    calculator_parameters: MaterialCalculatorParameters = MaterialCalculatorParameters()
+
+    def get_energy(self, material: Material):
+        """
+        Calculate the energy (or other metric) for a material.
+
+        Args:
+            material (Material): The material to calculate the interaction energy for.
+
+        Returns:
+            float: The interaction energy between the coordinates of the material,
+            calculated using the specified interaction function.
+        """
+        return self.calculator_parameters.interaction_function(material.coordinates, material.coordinates)
+
+
+class InterfaceMaterialCalculator(MaterialCalculator):
+    """
+    A specialized calculator for computing the interaction energy between a film and a substrate
+    in an interface material.
+
+    This class extends MaterialCalculator and uses additional parameters specific to interface materials,
+    such as the shadowing radius to detect surface atoms.
+
+    Args:
+        calculator_parameters (InterfaceMaterialCalculatorParameters): Parameters that include the
+        shadowing radius and interaction function.
+    """
+
+    calculator_parameters: InterfaceMaterialCalculatorParameters = InterfaceMaterialCalculatorParameters()
+
+    def get_energy(
+        self,
+        material: Material,
+        shadowing_radius: float = 2.5,
+        interaction_function: Callable = sum_of_inverse_distances_squared,
+    ) -> float:
+        """
+        Calculate the interaction energy between the film and substrate in an interface material.
+
+        This method uses the shadowing radius to detect surface atoms and applies the given
+        interaction function to calculate the interaction between the film and substrate.
+
+        Args:
+            material (Material): The interface Material object consisting of both the film and substrate.
+            shadowing_radius (float): The radius used to detect surface atoms for the interaction
+                                      calculation. Defaults to 2.5 Å.
+            interaction_function (Callable): A function to compute the interaction between the film and
+                                             substrate. Defaults to sum_of_inverse_distances_squared.
+
+        Returns:
+            float: The calculated interaction energy between the film and substrate.
+        """
+        film_material = get_interface_part(material, part=InterfacePartsEnum.FILM)
+        substrate_material = get_interface_part(material, part=InterfacePartsEnum.SUBSTRATE)
+
+        film_surface_atom_indices = get_surface_atom_indices(
+            film_material, SurfaceTypes.BOTTOM, shadowing_radius=shadowing_radius
+        )
+        substrate_surface_atom_indices = get_surface_atom_indices(
+            substrate_material, SurfaceTypes.TOP, shadowing_radius=shadowing_radius
+        )
+
+        film_surface_atom_coordinates = film_material.basis.coordinates
+        film_surface_atom_coordinates.filter_by_ids(film_surface_atom_indices)
+        substrate_surface_atom_coordinates = substrate_material.basis.coordinates
+        substrate_surface_atom_coordinates.filter_by_ids(substrate_surface_atom_indices)
+
+        film_coordinates_values = np.array(film_surface_atom_coordinates.values)
+        substrate_coordinates_values = np.array(substrate_surface_atom_coordinates.values)
+
+        return interaction_function(film_coordinates_values, substrate_coordinates_values)

src/py/mat3ra/made/tools/calculate/interaction_functions.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+
+def sum_of_inverse_distances_squared(
+    coordinates_1: np.ndarray, coordinates_2: np.ndarray, epsilon: float = 1e-12
+) -> float:
+    """
+    Calculate the sum of inverse squares of distances between two sets of coordinates.
+
+    Args:
+        coordinates_1 (np.ndarray): The first set of coordinates, shape (N1, 3).
+        coordinates_2 (np.ndarray): The second set of coordinates, shape (N2, 3).
+        epsilon (float): Small value to prevent division by zero.
+
+    Returns:
+        float: The calculated sum.
+    """
+    differences = coordinates_1[:, np.newaxis, :] - coordinates_2[np.newaxis, :, :]  # Shape: (N1, N2, 3)
+    distances_squared = np.sum(differences**2, axis=2)  # Shape: (N1, N2)
+    distances_squared = np.where(distances_squared == 0, epsilon, distances_squared)
+    inv_distances_squared = -1 / distances_squared
+    total = np.sum(inv_distances_squared)
+    return float(total)

src/py/mat3ra/made/tools/modify.py

@@ -466,7 +466,7 @@ def rotate_material(material: Material, axis: List[int], angle: float) -> Materi
     return Material(from_ase(atoms))
 
 
-def displace_interface(
+def displace_interface_part(
     interface: Material,
     displacement: List[float],
     label: InterfacePartsEnum = InterfacePartsEnum.FILM,
@@ -500,3 +500,15 @@ def displace_interface(
     new_material.to_crystal()
     new_material = wrap_to_unit_cell(new_material)
     return new_material
+
+
+def get_interface_part(
+    interface: Material,
+    part: InterfacePartsEnum = InterfacePartsEnum.FILM,
+) -> Material:
+    if interface.metadata["build"]["configuration"]["type"] != "InterfaceConfiguration":
+        raise ValueError("The material is not an interface.")
+    interface_part_material = interface.clone()
+    film_atoms_basis = interface_part_material.basis.filter_atoms_by_labels([int(part)])
+    interface_part_material.basis = film_atoms_basis
+    return interface_part_material

src/py/mat3ra/made/tools/optimize.py

@@ -0,0 +1,55 @@
+from typing import Any, Callable, Dict, List, Optional, Tuple
+
+import numpy as np
+from mat3ra.made.material import Material
+
+
+def evaluate_calculator_on_xy_grid(
+    material: Material,
+    calculator_function: Callable[[Material], Any],
+    modifier: Optional[Callable] = None,
+    modifier_parameters: Dict[str, Any] = {},
+    grid_size_xy: Tuple[int, int] = (10, 10),
+    grid_offset_position: List[float] = [0, 0],
+    grid_range_x: Tuple[float, float] = (-0.5, 0.5),
+    grid_range_y: Tuple[float, float] = (-0.5, 0.5),
+    use_cartesian_coordinates: bool = False,
+) -> Tuple[List[np.ndarray], np.ndarray]:
+    """
+    Calculate a property on a grid of x-y positions.
+
+    Args:
+        material (Material): The material object.
+        modifier (Callable): The modifier function to apply to the material.
+        modifier_parameters (Dict[str, Any]): The parameters to pass to the modifier.
+        calculator_function (Callable): The calculator function to apply to the modified material.
+        grid_size_xy (Tuple[int, int]): The size of the grid in x and y directions.
+        grid_offset_position (List[float]): The offset position of the grid, in Angstroms or crystal coordinates.
+        grid_range_x (Tuple[float, float]): The range to search in x direction, in Angstroms or crystal coordinates.
+        grid_range_y (Tuple[float, float]): The range to search in y direction, in Angstroms or crystal coordinates.
+        use_cartesian_coordinates (bool): Whether to use Cartesian coordinates.
+
+    Returns:
+        Tuple[List[np.ndarray[float]], np.ndarray[float]]: The x-y positions and the calculated property values.
+    """
+    x_values = np.linspace(grid_range_x[0], grid_range_x[1], grid_size_xy[0]) + grid_offset_position[0]
+    y_values = np.linspace(grid_range_y[0], grid_range_y[1], grid_size_xy[1]) + grid_offset_position[1]
+
+    xy_matrix = [x_values, y_values]
+    results_matrix = np.zeros(grid_size_xy)
+
+    for i, x in enumerate(x_values):
+        for j, y in enumerate(y_values):
+            if modifier is None:
+                modified_material = material
+            else:
+                modified_material = modifier(
+                    material,
+                    displacement=[x, y, 0],
+                    use_cartesian_coordinates=use_cartesian_coordinates,
+                    **modifier_parameters,
+                )
+            result = calculator_function(modified_material)
+            results_matrix[i, j] = result
+
+    return xy_matrix, results_matrix

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

@@ -3,18 +3,9 @@
 
 import numpy as np
 from mat3ra.made.material import Material
-from mat3ra.made.utils import ArrayWithIds
 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
 
 DEFAULT_SCALING_FACTOR = np.array([3, 3, 3])
 DEFAULT_TRANSLATION_VECTOR = 1 / DEFAULT_SCALING_FACTOR