Skip to content

Commit

Permalink
Update displacements tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@kavanase
kavanase committed Dec 7, 2024
1 parent fbedff5 commit fa928da
Show file tree
Hide file tree
Showing 3 changed files with 140 additions and 93 deletions.
13 changes: 9 additions & 4 deletions doped/utils/displacements.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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),
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand Down
4 changes: 2 additions & 2 deletions doped/utils/parsing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
216 changes: 129 additions & 87 deletions tests/test_displacements.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down

0 comments on commit fa928da

Please sign in to comment.