Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Fix consideration of translation symmetry for some (extremely rare) edge cases in LobsterEnv #4148

Draft
wants to merge 7 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions src/pymatgen/electronic_structure/cohp.py
Original file line number Diff line number Diff line change
Expand Up @@ -1042,6 +1042,17 @@ def is_spin_polarized(self) -> bool:
"""
return self._is_spin_polarized

@property
def translation(self) -> list[int, int, int]:
"""
Returns the translation vector with respect to the origin cell
as defined in LOBSTER.

Returns:
list[int, int, int]
"""
return self._translation

def icohpvalue(self, spin: Spin = Spin.up) -> float:
"""
Args:
Expand Down
63 changes: 46 additions & 17 deletions src/pymatgen/io/lobster/lobsterenv.py
Original file line number Diff line number Diff line change
Expand Up @@ -790,14 +790,9 @@ def _evaluate_ce(
raise ValueError("Please give two limits or leave them both at None")

# Find environments based on ICOHP values
(
list_icohps,
list_keys,
list_lengths,
list_neighisite,
list_neighsite,
list_coords,
) = self._find_environments(additional_condition, lowerlimit, upperlimit, only_bonds_to)
(list_icohps, list_keys, list_lengths, list_neighisite, list_neighsite, list_coords) = self._find_environments(
additional_condition, lowerlimit, upperlimit, only_bonds_to
)

self.list_icohps = list_icohps
self.list_lengths = list_lengths
Expand Down Expand Up @@ -939,11 +934,13 @@ def _find_environments(
lengths_from_ICOHPs,
neighbors_from_ICOHPs,
selected_ICOHPs,
translations_ICOHPs,
) = additional_conds

if len(neighbors_from_ICOHPs) > 0:
centralsite = site

copysite = copy.copy(centralsite)
centralsite.frac_coords - copysite.to_unit_cell().frac_coords
neighbors_by_distance_start = self.structure.get_sites_in_sphere(
pt=centralsite.coords,
r=np.max(lengths_from_ICOHPs) + 0.5,
Expand All @@ -955,48 +952,63 @@ def _find_environments(
list_distances = []
index_here_list = []
coords = []
translations_by_distance = []
for neigh_new in sorted(neighbors_by_distance_start, key=lambda x: x[1]):
site_here = neigh_new[0].to_unit_cell()
index_here = neigh_new[2]
index_here_list.append(index_here)
cell_here = neigh_new[3]

new_coords = [
site_here.frac_coords[0] + float(cell_here[0]),
site_here.frac_coords[1] + float(cell_here[1]),
site_here.frac_coords[2] + float(cell_here[2]),
]
coords.append(site_here.lattice.get_cartesian_coords(new_coords))

# new_site = PeriodicSite(
# species=site_here.species_string,
# coords=site_here.lattice.get_cartesian_coords(new_coords),
# lattice=site_here.lattice,
# to_unit_cell=False,
# coords_are_cartesian=True,
# )
neighbors_by_distance.append(neigh_new[0])
list_distances.append(neigh_new[1])
translations_by_distance.append(cell_here)
_list_neighsite = []
_list_neighisite = []
copied_neighbors_from_ICOHPs = copy.copy(neighbors_from_ICOHPs)
copied_distances_from_ICOHPs = copy.copy(lengths_from_ICOHPs)
copied_translations_from_ICOHPs = copy.copy(translations_ICOHPs)
_neigh_coords = []
_neigh_frac_coords = []

for neigh_idx, neigh in enumerate(neighbors_by_distance):
index_here2 = index_here_list[neigh_idx]

for dist_idx, dist in enumerate(copied_distances_from_ICOHPs):
if (
np.isclose(dist, list_distances[neigh_idx], rtol=1e-4)
and copied_neighbors_from_ICOHPs[dist_idx] == index_here2
and (
(
copied_translations_from_ICOHPs[dist_idx][0]
== -translations_by_distance[neigh_idx][0]
and copied_translations_from_ICOHPs[dist_idx][1]
== -translations_by_distance[neigh_idx][1]
and copied_translations_from_ICOHPs[dist_idx][2]
== -translations_by_distance[neigh_idx][2]
)
or (
copied_translations_from_ICOHPs[dist_idx][0]
== translations_by_distance[neigh_idx][0]
and copied_translations_from_ICOHPs[dist_idx][1]
== translations_by_distance[neigh_idx][1]
and copied_translations_from_ICOHPs[dist_idx][2]
== translations_by_distance[neigh_idx][2]
)
Comment on lines +995 to +1002
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could raise a warning if this fall-back option is taken. I think with updating to reading the CONTCAR, most issues should be fixed, however.

)
):
_list_neighsite.append(neigh)
_list_neighisite.append(index_here2)
_neigh_coords.append(coords[neigh_idx])
_neigh_frac_coords.append(neigh.frac_coords)
del copied_distances_from_ICOHPs[dist_idx]
del copied_neighbors_from_ICOHPs[dist_idx]
del copied_translations_from_ICOHPs[dist_idx]
break

list_neighisite.append(_list_neighisite)
Expand Down Expand Up @@ -1042,6 +1054,7 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs: list[float] = []
neighbors_from_ICOHPs: list[int] = []
icohps_from_ICOHPs: list[IcohpValue] = []
translation_from_ICOHPs: list[list[int, int, int]] = []

for key, icohp in icohps.items():
atomnr1 = self._get_atomnumber(icohp._atom1)
Expand All @@ -1062,11 +1075,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
# add translation to icohp value
translation_from_ICOHPs.append(icohp.translation)
elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# ONLY_ANION_CATION_BONDS
elif additional_condition == 1:
Expand All @@ -1076,12 +1092,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# NO_ELEMENT_TO_SAME_ELEMENT_BONDS
elif additional_condition == 2:
Expand All @@ -1091,12 +1109,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# ONLY_ANION_CATION_BONDS_AND_NO_ELEMENT_TO_SAME_ELEMENT_BONDS
elif additional_condition == 3:
Expand All @@ -1108,12 +1128,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# ONLY_ELEMENT_TO_OXYGEN_BONDS
elif additional_condition == 4:
Expand All @@ -1123,12 +1145,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# DO_NOT_CONSIDER_ANION_CATION_BONDS
elif additional_condition == 5:
Expand All @@ -1138,12 +1162,14 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

# ONLY_CATION_CATION_BONDS
elif additional_condition == 6 and val1 > 0.0 and val2 > 0.0: # type: ignore[operator]
Expand All @@ -1152,18 +1178,21 @@ def _find_relevant_atoms_additional_condition(
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

elif atomnr2 == site_idx:
neighbors_from_ICOHPs.append(atomnr1)
lengths_from_ICOHPs.append(icohp._length)
icohps_from_ICOHPs.append(icohp.summed_icohp)
keys_from_ICOHPs.append(key)
translation_from_ICOHPs.append(icohp.translation)

return (
keys_from_ICOHPs,
lengths_from_ICOHPs,
neighbors_from_ICOHPs,
icohps_from_ICOHPs,
translation_from_ICOHPs,
)

@staticmethod
Expand Down
Loading