diff --git a/structuretoolkit/analyse/pymatgen.py b/structuretoolkit/analyse/pymatgen.py new file mode 100644 index 000000000..41655485b --- /dev/null +++ b/structuretoolkit/analyse/pymatgen.py @@ -0,0 +1,75 @@ +from pymatgen.io.ase import AseAtomsAdaptor +from pymatgen.analysis.local_env import VoronoiNN +from pymatgen.core import Structure, Element +import numpy as np +import pandas as pd + + +def get_stats(property_list, property_str): + """ + Calculate statistical properties of a list of values. + Parameters: + property_list (list): A list of numerical values for which statistics are calculated. + property_str (str): A string prefix to be used in the resulting statistical property names. + Returns: + dict: A dictionary containing statistical properties with keys in the format: + "{property_str}_{statistic}" where statistic can be "std" (standard deviation), + "mean" (mean), "min" (minimum), and "max" (maximum). + Example: + >>> values = [1, 2, 3, 4, 5] + >>> get_stats(values, "example") + {'example_std': 1.4142135623730951, + 'example_mean': 3.0, + 'example_min': 1, + 'example_max': 5} + """ + return { + f"{property_str}_std": np.std(property_list), + f"{property_str}_mean": np.mean(property_list), + f"{property_str}_min": np.min(property_list), + f"{property_str}_max": np.max(property_list), + } + + +def VoronoiSiteFeaturiser(structure, site): + """ + Calculate various Voronoi-related features for a specific site in a crystal structure. + Parameters: + structure (ase.Atoms or pymatgen.Structure): The crystal structure. + site (int): The index of the site in the crystal structure. + Returns: + pandas.DataFrame: A DataFrame containing computed Voronoi features for the specified site. + Columns include VorNN_CoordNo, VorNN_tot_vol, VorNN_tot_area, as well as + statistics for volumes, vertices, areas, and distances. + Example: + >>> from pymatgen import Structure + >>> structure = Structure.from_file("example.cif") + >>> VoronoiSiteFeaturiser(structure, 0) + VorNN_CoordNo VorNN_tot_vol VorNN_tot_area volumes_std volumes_mean ... + 0 7.0 34.315831 61.556747 10.172586 34.315831 ... + """ + structure = AseAtomsAdaptor().get_structure(structure) + coord_no = VoronoiNN().get_cn(structure=structure, n=site) + site_info_dict = VoronoiNN().get_voronoi_polyhedra(structure, site) + volumes = [site_info_dict[polyhedra]["volume"] for polyhedra in site_info_dict] + vertices = [site_info_dict[polyhedra]["n_verts"] for polyhedra in site_info_dict] + distances = [site_info_dict[polyhedra]["face_dist"] for polyhedra in site_info_dict] + areas = [site_info_dict[polyhedra]["area"] for polyhedra in site_info_dict] + + total_area = np.sum(areas) + total_volume = np.sum(volumes) + + data = { + "VorNN_CoordNo": coord_no, + "VorNN_tot_vol": total_volume, + "VorNN_tot_area": total_area, + } + + data_str_list = ["volumes", "vertices", "areas", "distances"] + + for i, value_list in enumerate([volumes, vertices, areas, distances]): + stats = get_stats(value_list, f"VorNN_{data_str_list[i]}") + data.update(stats) + + df = pd.DataFrame(data, index=[site]) + return df diff --git a/structuretoolkit/analyse/strain.py b/structuretoolkit/analyse/strain.py index 3632c9dda..f9ebb6aaa 100644 --- a/structuretoolkit/analyse/strain.py +++ b/structuretoolkit/analyse/strain.py @@ -6,7 +6,6 @@ class Strain: - """ Calculate local strain of each atom following the Lagrangian strain tensor: diff --git a/structuretoolkit/analyse/symmetry.py b/structuretoolkit/analyse/symmetry.py index bebeeae37..2da660fee 100644 --- a/structuretoolkit/analyse/symmetry.py +++ b/structuretoolkit/analyse/symmetry.py @@ -24,7 +24,6 @@ class Symmetry(dict): - """ Return a class for operations related to box symmetries. Main attributes: diff --git a/tests/test_pymatgen.py b/tests/test_pymatgen.py index 1fc887975..746308d68 100644 --- a/tests/test_pymatgen.py +++ b/tests/test_pymatgen.py @@ -4,12 +4,11 @@ from ase.constraints import FixAtoms from structuretoolkit.common import pymatgen_to_ase, ase_to_pymatgen - try: from pymatgen.core import Structure, Lattice - + from structuretoolkit.analyse.pymatgen import VoronoiSiteFeaturiser skip_pymatgen_test = False -except ImportError: +except (ImportError, ModuleNotFoundError): skip_pymatgen_test = True @@ -196,3 +195,56 @@ def test_pyiron_to_pymatgen_conversion(self): ), "Failed to produce equivalent sel_dyn when both magmom + sel_dyn are present!", ) + +@unittest.skipIf( + skip_pymatgen_test, "pymatgen is not installed, so the pymatgen tests are skipped." +) +class TestVoronoiSiteFeaturiser(unittest.TestCase): + def setUp(self): + self.example_structure = bulk("Fe") + + def assertListsAlmostEqual(self, list1, list2, decimal=4): + """ + Check if two lists are approximately equal up to a specified number of decimal places. + + Parameters: + list1 (list): The first list for comparison. + list2 (list): The second list for comparison. + decimal (int): The number of decimal places to consider for comparison. + + Raises: + AssertionError: Raised if the lists are not approximately equal. + """ + self.assertEqual(len(list1), len(list2), "Lists have different lengths") + + for i in range(len(list1)): + self.assertAlmostEqual(list1[i], list2[i], places=decimal, + msg=f"Lists differ at index {i}: {list1[i]} != {list2[i]}") + + def test_VoronoiSiteFeaturiser(self): + # Calculate the expected output manually + expected_output = { + "VorNN_CoordNo": 14, + "VorNN_tot_vol": 11.819951, + "VorNN_tot_area": 27.577769, + "VorNN_volumes_std": 0.304654, + "VorNN_volumes_mean": 0.844282, + "VorNN_volumes_min": 0.492498, + "VorNN_volumes_max": 1.10812, + "VorNN_vertices_std": 0.989743, + "VorNN_vertices_mean": 5.142857, + "VorNN_vertices_min": 4, + "VorNN_vertices_max": 6, + "VorNN_areas_std": 0.814261, + "VorNN_areas_mean": 1.969841, + "VorNN_areas_min": 1.029612, + "VorNN_areas_max": 2.675012, + "VorNN_distances_std": 0.095141, + "VorNN_distances_mean": 1.325141, + "VorNN_distances_min": 1.242746, + "VorNN_distances_max": 1.435 + } + + # Call the function with the example structure + df = VoronoiSiteFeaturiser(self.example_structure, 0) + self.assertListsAlmostEqual(df.values.tolist()[0], list(expected_output.values()), decimal=4)