diff --git a/doped/utils/displacements.py b/doped/utils/displacements.py index 9214d036..d45c2c5a 100644 --- a/doped/utils/displacements.py +++ b/doped/utils/displacements.py @@ -40,6 +40,7 @@ def calc_site_displacements( relaxed_distances: bool = False, relative_to_defect: bool = False, vector_to_project_on: Optional[list] = None, + threshold: float = 2.0, ) -> pd.DataFrame: """ Calculates the site displacements in the defect supercell, relative to the @@ -65,7 +66,10 @@ def calc_site_displacements( the `Displacement wrt defect` key of the returned dictionary. vector_to_project_on (list): Direction to project the site displacements along - (e.g. [0, 0, 1]). Defaults to None. + (e.g. [0, 0, 1]). Defaults to ``None``. + threshold (float): + If the distance between a pair of matched sites is larger than this, + then a warning will be thrown. Default is 2.0 Å. Returns: ``pandas`` ``DataFrame`` with site displacements (compared to pristine supercell), @@ -75,7 +79,7 @@ def calc_site_displacements( defect_site = defect_sc_with_site[defect_site_index] # Map sites in defect supercell to bulk supercell: - mappings = get_site_mapping_indices(defect_sc_with_site, bulk_sc) + mappings = get_site_mapping_indices(defect_sc_with_site, bulk_sc, threshold=threshold) mappings_dict = {i[1]: i[2] for i in mappings} # {defect_sc_index: bulk_sc_index} disp_dict: dict[str, list] = { # mapping defect site index (in defect sc) to displacement @@ -138,10 +142,11 @@ def calc_site_displacements( # sort each list in disp dict by index of species in bulk element list, then by distance to defect: element_list = _get_element_list(defect_entry) disp_df = pd.DataFrame(disp_dict) - # Sort by species, then distance to defect, then index: + # Sort by species, then distance to defect, then displacement magnitude (reversed), then index: disp_df = disp_df.sort_values( - by=["Species", "Distance to defect", "Index (defect supercell)"], + by=["Species", "Distance to defect", "Displacement", "Index (defect supercell)"], key=lambda col: col.map(element_list.index) if col.name == "Species" else col, + ascending=[True, True, False, True], ) disp_df = _round_floats(disp_df, 4) # round numerical values to 4 dp # reorder columns as species, distance, displacement, etc, then index last: diff --git a/doped/utils/parsing.py b/doped/utils/parsing.py index f87dbd80..96aa679e 100644 --- a/doped/utils/parsing.py +++ b/doped/utils/parsing.py @@ -894,8 +894,8 @@ def get_site_mapping_indices( then a warning will be thrown. Default is 2.0 Å. dists_only (bool): Whether to return only the distances between matched sites, rather - than a list of lists containing the distance, index in struct1 - and index in struct2. Default is False. + than a list of lists containing the distance, index in ``struct1`` + and index in ``struct2``. Default is ``False``. Returns: list: diff --git a/tests/test_displacements.py b/tests/test_displacements.py index 116e837d..de7030e5 100644 --- a/tests/test_displacements.py +++ b/tests/test_displacements.py @@ -45,82 +45,115 @@ def setUp(self): def test_calc_site_displacements(self): """ - Test calc_site_displacements() function. + Test ``calc_site_displacements()`` function. """ - defect_entry = self.v_Cd_0_defect_entry # Neutral Cd vacancy - disp_df = calc_site_displacements(defect_entry) - for i, disp in [ - (0, [0.0572041, 0.00036486, -0.01794981]), - (15, [0.11715445, -0.03659073, 0.01312027]), - ]: - np.allclose(disp_df["Displacement"].iloc[i], np.array(disp)) - # Check distance - for i, dist in [ - (0, 0.0), - (15, 6.543643778883931), - ]: - assert np.isclose(disp_df["Distance to defect"].iloc[i], dist) - # Test displacements added to defect_entry: - disp_metadata = defect_entry.calculation_metadata["site_displacements"] - for i, disp in [ - (0, [0.0572041, 0.00036486, -0.01794981]), - (15, [0.11715445, -0.03659073, 0.01312027]), - ]: - np.allclose(disp_metadata["displacements"][i], np.array(disp)) - # Test displacement of vacancy removed before adding to calculation_metadata - assert len(disp_metadata["distances"]) == 63 # Cd Vacancy so 63 sites - # Test relative displacements from defect - disp_df = calc_site_displacements(defect_entry, relative_to_defect=True) - for i, disp in [ - (0, 0.0), - (1, -0.1166), - ]: - assert np.isclose(disp_df["Displacement wrt defect"].iloc[i], disp, atol=1e-3) - # Test projection along Te dimer direction (1,1,0) - disp_df = calc_site_displacements(defect_entry, vector_to_project_on=[1, 1, 0]) - for i, dist, disp in [ - (32, 2.177851642, 0.980779911), # index, distance, displacement - (33, 2.234858368, -0.892888144), - ]: - assert np.isclose(disp_df["Displacement projected along vector"].iloc[i], disp, atol=1e-3) - assert np.isclose(disp_df["Distance to defect"].iloc[i], dist, atol=1e-3) - # Test projection along (-1,-1,-1) for V_Cd^-1 - defect_entry = self.v_Cd_m1_defect_entry - disp_df = calc_site_displacements(defect_entry, vector_to_project_on=[-1, -1, -1]) - indexes = (32, 33, 34, 35) # Defect NNs - distances = (2.5850237041739614, 2.5867590623267396, 2.5867621810347914, 3.0464198655727284) - disp_parallel = (0.10857369429255698, 0.10824441910793342, 0.10824525022932621, 0.2130514712472405) - disp_perpendicular = ( - 0.22517568241969493, - 0.22345433515763177, - 0.22345075557153446, - 0.0002337424502264542, - ) - for index, dist, disp_paral, disp_perp in zip( - indexes, distances, disp_parallel, disp_perpendicular - ): - assert np.isclose(disp_df["Distance to defect"].iloc[index], dist, atol=1e-3) - assert np.isclose( - disp_df["Displacement projected along vector"].iloc[index], disp_paral, atol=1e-3 + for relaxed_distances in [False, True]: + print(f"Testing calc_site_displacements with relaxed_distances={relaxed_distances}") + defect_entry = self.v_Cd_0_defect_entry # Neutral Cd vacancy + disp_df = calc_site_displacements(defect_entry, relaxed_distances=relaxed_distances) + for i, disp in [ + (0, [0.0572041, 0.00036486, -0.01794981]), + (15, [0.11715445, -0.03659073, 0.01312027]), + ]: + np.allclose(disp_df["Displacement"].iloc[i], np.array(disp)) + # Check distance + for i, dist in [ + (0, 0.0), + (15, 6.54), + ]: + assert np.isclose(disp_df["Distance to defect"].iloc[i], dist, atol=2e-2) + # Test displacements added to defect_entry: + disp_metadata = defect_entry.calculation_metadata["site_displacements"] + for i, disp in [ + (0, [0.0572041, 0.00036486, -0.01794981]), + (15, [0.11715445, -0.03659073, 0.01312027]), + ]: + np.allclose(disp_metadata["displacements"][i], np.array(disp)) + # Test displacement of vacancy removed before adding to calculation_metadata + assert len(disp_metadata["distances"]) == 63 # Cd Vacancy so 63 sites + # Test relative displacements from defect + disp_df = calc_site_displacements( + defect_entry, relaxed_distances=relaxed_distances, relative_to_defect=True ) - assert np.isclose( - disp_df["Displacement perpendicular to vector"].iloc[index], disp_perp, atol=1e-3 + disp_tuples = [ + (0, 0.0), + ] + ( + [ + (1, -0.1166), + ] + if relaxed_distances + else [ + (4, -0.0796), + ] ) + for i, disp in disp_tuples: + assert np.isclose(disp_df["Displacement wrt defect"].iloc[i], disp, atol=1e-3) - # Substitution: - disp_df = calc_site_displacements(self.Te_Cd_1_defect_entry) - for i, disp in [ - (0, [0.00820645, 0.00821417, -0.00815738]), - (15, [-0.00639524, 0.00639969, -0.01407927]), - ]: - np.allclose(disp_df["Displacement"].iloc[i], np.array(disp), atol=1e-3) - # Interstitial: - disp_df = calc_site_displacements(self.Te_i_1_defect_entry) - for i, disp in [ - (0, [-0.03931121, 0.01800569, 0.04547194]), - (15, [-0.04850126, -0.01378455, 0.05439607]), - ]: - np.allclose(disp_df["Displacement"].iloc[i], np.array(disp), atol=1e-3) + # Test projection along Te dimer direction (1,1,0) + disp_df = calc_site_displacements( + defect_entry, relaxed_distances=relaxed_distances, vector_to_project_on=[1, 1, 0] + ) + if relaxed_distances: + disp_tuples = [ + (32, 2.177851642, 0.980779911), # index, distance, displacement + (33, 2.234858368, -0.892888144), + ] + else: + disp_tuples = [ + (32, 2.83337, 0.980779911), # index, distance, displacement + (33, 2.83337, -0.892888144), + ] + for i, dist, disp in disp_tuples: + assert np.isclose(disp_df["Displacement projected along vector"].iloc[i], disp, atol=1e-3) + assert np.isclose(disp_df["Distance to defect"].iloc[i], dist, atol=1e-3) + + # Test projection along (-1,-1,-1) for V_Cd^-1 + defect_entry = self.v_Cd_m1_defect_entry + disp_df = calc_site_displacements( + defect_entry, relaxed_distances=relaxed_distances, vector_to_project_on=[-1, -1, -1] + ) + indexes = (32, 33, 34, 35) # Defect NNs + distances = (2.5850237041739614, 2.5867590623267396, 2.5867621810347914, 3.0464198655727284) + disp_parallel = ( + 0.10857369429255698, + 0.10824441910793342, + 0.10824525022932621, + 0.2130514712472405, + ) + disp_perpendicular = ( + 0.22517568241969493, + 0.22345433515763177, + 0.22345075557153446, + 0.0002337424502264542, + ) + for index, dist, disp_paral, disp_perp in zip( + indexes, distances, disp_parallel, disp_perpendicular + ): + if relaxed_distances: + assert np.isclose(disp_df["Distance to defect"].iloc[index], dist, atol=1e-3) + else: # all the same NN distance: + assert np.isclose(disp_df["Distance to defect"].iloc[index], 2.83337, atol=1e-3) + assert np.isclose( + disp_df["Displacement projected along vector"].iloc[index], disp_paral, atol=1e-3 + ) + assert np.isclose( + disp_df["Displacement perpendicular to vector"].iloc[index], disp_perp, atol=1e-2 + ) + + # Substitution: + disp_df = calc_site_displacements(self.Te_Cd_1_defect_entry) + for i, disp in [ + (0, [0.00820645, 0.00821417, -0.00815738]), + (15, [-0.00639524, 0.00639969, -0.01407927]), + ]: + np.allclose(disp_df["Displacement"].iloc[i], np.array(disp), atol=1e-3) + # Interstitial: + disp_df = calc_site_displacements(self.Te_i_1_defect_entry) + for i, disp in [ + (0, [-0.03931121, 0.01800569, 0.04547194]), + (15, [-0.04850126, -0.01378455, 0.05439607]), + ]: + np.allclose(disp_df["Displacement"].iloc[i], np.array(disp), atol=1e-3) def test_plot_site_displacements_error(self): # Check ValueError raised if user sets both separated_by_direction and vector_to_project_on @@ -188,21 +221,30 @@ def test_calc_displacements_ellipsoid(self): ), ]: print("Testing displacement ellipsoid for", entry.name) - ellipsoid_center, ellipsoid_radii, ellipsoid_rotation, anisotropy_df = ( - calc_displacements_ellipsoid(entry, quantile=0.8) - ) - assert np.allclose(ellipsoid_center, np.array(ellipsoid_center_benchmark), atol=1e-3) - assert np.allclose(ellipsoid_radii, np.array(ellipsoid_radii_benchmark), atol=1e-3) - assert np.allclose(ellipsoid_rotation, np.array(ellipsoid_rotation_benchmark), atol=1e-3) - assert anisotropy_df["Longest Radius"].to_numpy()[0] == ellipsoid_radii[2] - assert ( - anisotropy_df["2nd_Longest/Longest"].to_numpy()[0] - == ellipsoid_radii[1] / ellipsoid_radii[2] - ) - assert ( - anisotropy_df["3rd_Longest/Longest"].to_numpy()[0] - == ellipsoid_radii[0] / ellipsoid_radii[2] - ) + for relaxed_distances in [False, True]: + ellipsoid_center, ellipsoid_radii, ellipsoid_rotation, anisotropy_df = ( + calc_displacements_ellipsoid(entry, quantile=0.8, relaxed_distances=relaxed_distances) + ) + if relaxed_distances: + assert np.allclose(ellipsoid_center, np.array(ellipsoid_center_benchmark), atol=1e-3) + assert np.allclose(ellipsoid_radii, np.array(ellipsoid_radii_benchmark), atol=1e-3) + assert np.allclose( + ellipsoid_rotation, np.array(ellipsoid_rotation_benchmark), atol=1e-3 + ) + + else: + assert np.allclose(ellipsoid_center, np.array(ellipsoid_center_benchmark), atol=2.0) + assert np.allclose(ellipsoid_radii, np.array(ellipsoid_radii_benchmark), atol=2.0) + + assert anisotropy_df["Longest Radius"].to_numpy()[0] == ellipsoid_radii[2] + assert ( + anisotropy_df["2nd_Longest/Longest"].to_numpy()[0] + == ellipsoid_radii[1] / ellipsoid_radii[2] + ) + assert ( + anisotropy_df["3rd_Longest/Longest"].to_numpy()[0] + == ellipsoid_radii[0] / ellipsoid_radii[2] + ) @custom_mpl_image_compare(filename="v_Cd_0_disp_proj_plot.png") def test_plot_site_displacements_proj(self):