Merge pull request #17 from Materials-Data-Science-and-Informatics/ad…

…d_polyhedral_template_maching

Feature: add polyhedral matching

src/atomid/annotate.py

@@ -4,9 +4,13 @@
 
 import atomrdf as ardf
 from ase.io import read as ase_read
+from ovito.data import DataCollection
+from ovito.io import import_file
+from ovito.modifiers import PolyhedralTemplateMatchingModifier
 
 from atomid.crystal.structure_identification import (
-    find_lattice_parameter,
+    analyse_polyhedral_template_matching_data,
+    find_lattice_parameter_2,
     get_crystal_structure_using_cna,
 )
 from atomid.point_defect_analysis.wigner_seitz_method import analyze_defects
@@ -32,23 +36,55 @@ def read_crystal_structure_file(self, data_file: str, format: str) -> None:
         """
         self.ase_crystal = ase_read(data_file, format=format)
         kg = ardf.KnowledgeGraph()
+
         crystal_structure = ardf.System.read.file(
             filename=self.ase_crystal, format="ase", graph=kg
         )
 
+        self.ovito_pipeline = import_file(data_file)
+
         self.kg = kg
         self.system = crystal_structure
 
+    def get_polyhedral_template_matching_data(self) -> DataCollection:
+        """Get the polyhedral template matching data from the ovito pipeline.
+
+        Parameters
+        ----------
+        ovito_pipeline : OvitoPipeline
+            The ovito pipeline object.
+
+        Returns
+        -------
+        dict
+            The polyhedral template matching data.
+        """
+        self.ovito_pipeline.modifiers.append(
+            PolyhedralTemplateMatchingModifier(output_interatomic_distance=True)
+        )
+        data = self.ovito_pipeline.compute()
+        return data
+
     def annotate_crystal_structure(self) -> None:
         """Identify and annotate the crystal structure.
 
         This method identifies the crystal structure using Common Neighbour Analysis
         and lattice constant using radial distribution function.
         """
         crystal_type = get_crystal_structure_using_cna(self.system)
+        # get crystal structure from polyhedral template matching
+        structure_data = self.get_polyhedral_template_matching_data()
+        structure_type_atoms = structure_data.particles["Structure Type"][...]  # noqa
+        structure_id, crystal_type = analyse_polyhedral_template_matching_data(
+            structure_type_atoms
+        )
 
         if crystal_type != "others":
-            lattice_constants = find_lattice_parameter(self.system, crystal_type)
+            interatomic_distance = structure_data.particles["Interatomic Distance"][...]  # noqa
+            lattice_constants = find_lattice_parameter_2(
+                interatomic_distance, structure_type_atoms, int(structure_id)
+            )
+            # lattice_constants = find_lattice_parameter(self.system, crystal_type)
             self.system = None
             self.kg = ardf.KnowledgeGraph()
             self.system = ardf.System.read.file(

src/atomid/crystal/structure_identification.py

@@ -2,11 +2,48 @@
 
 import logging
 from math import sqrt
+from typing import Tuple
 
+import numpy as np
 from atomrdf import System
 from scipy.signal import find_peaks
 
 
+def get_crsytal_structure_from_id(id: int) -> str:
+    """
+    Get the crystal structure type from the given ID.
+
+    Parameters
+    ----------
+    id : int
+        The ID of the crystal structure.
+
+    Returns
+    -------
+    str
+        The crystal structure type.
+    """
+    structure_type = {
+        0: "other",
+        1: "fcc",
+        2: "hcp",
+        3: "bcc",
+        4: "ico",
+        5: "sc",
+        6: "cubic diamond",
+        7: "hex diamond",
+        8: "graphene",
+    }
+    return structure_type.get(id, "other")
+
+
+def analyse_polyhedral_template_matching_data(atoms_structure_type: np.ndarry) -> Tuple:
+    """Analyse polyhedral template data to identify crystal structure."""
+    unique, counts = np.unique(atoms_structure_type, return_counts=True)
+    structure_id = unique[np.argmax(counts)]
+    return structure_id, get_crsytal_structure_from_id(structure_id)
+
+
 def get_crystal_structure_using_cna(pyscal_system: System) -> str:
     """
     Get the crystal structure using adaptive common neighbour analysis.
@@ -78,6 +115,39 @@ def analyse_diamond_structures(pyscal_system: System) -> str:
     return str(diamond_type)
 
 
+def find_lattice_parameter_2(
+    interatomic_distance: np.ndarray,
+    structure_type_atoms: np.ndarray,
+    crystal_type_id: int,
+) -> float:
+    # create mask for the structure type matching the crystal type id
+    mask = structure_type_atoms == crystal_type_id
+
+    # get the interatomic distances for the given structure type
+    filtered_interatomic_distance = interatomic_distance[mask]
+
+    # calculate mean of the interatomic distances
+    mean_interactomic_distance = np.mean(filtered_interatomic_distance)
+
+    # calculate multiplier for the lattice parameter
+
+    multiplier = {
+        1: sqrt(2),
+        2: sqrt(2),
+        3: sqrt(4 / 3),
+        4: 4 / sqrt(3),
+        5: 1,
+        6: sqrt(16 / 3),
+        7: sqrt(16 / 3),
+        8: sqrt(2),
+    }
+
+    lattice_parameter: float = round(
+        mean_interactomic_distance * multiplier[crystal_type_id], 3
+    )
+    return lattice_parameter
+
+
 def find_lattice_parameter(crystal_system: System, lattice_type: str) -> float:
     """
     Calculate the lattice parameter of a crystal structure.

src/atomid/point_defect_analysis/wigner_seitz_method.py

@@ -136,6 +136,7 @@ def calculate_defects(
         },
         "Substitutions": {
             "count": len(substitutions),
-            "fraction": len(substitutions) / len(reference_array),
+            "fraction": (len(substitutions) - len(interstitials))
+            / len(reference_array),
         },
     }