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.

Sign up for GitHub

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

Jump to bottom

Added _quadrant_integral function in CrystalNN #3974

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Added _quadrant_integral function in CrystalNN #3974

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fee3843
Fixed failing tests, added link to the issue that describes problem w…
Aug 7, 2024
91df1d0
pre-commit auto-fixes
pre-commit-ci[bot] Aug 7, 2024
304d42b
Merge branch 'master' into master
janosh Aug 24, 2024
4e61a26
Merge branch 'master' into master
janosh Sep 2, 2024
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
50 changes: 42 additions & 8 deletions src/pymatgen/analysis/local_env.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1162,7 +1162,7 @@ def _is_in_targets(site, targets):
targets ([Element]) List of elements

Returns:
bool: Whether this site contains a certain list of elements
boolean: Whether this site contains a certain list of elements
Copy link
Member

Choose a reason for hiding this comment

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

this should not be reverted

"""
elems = _get_elements(site)
return all(elem in targets for elem in elems)
Expand Down Expand Up @@ -1219,7 +1219,7 @@ def __init__(

# Update any user preference elemental radii
if el_radius_updates:
self.el_radius |= el_radius_updates
self.el_radius.update(el_radius_updates)
Copy link
Member

Choose a reason for hiding this comment

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

same here


@property
def structures_allowed(self) -> bool:
Expand Down Expand Up @@ -1985,7 +1985,7 @@ def get_okeeffe_distance_prediction(el1, el2):
"""Get an estimate of the bond valence parameter (bond length) using
the derived parameters from 'Atoms Sizes and Bond Lengths in Molecules
and Crystals' (O'Keeffe & Brese, 1991). The estimate is based on two
experimental parameters: r and c. The value for r is based off radius,
experimental parameters: r and c. The value for r is based off radius,
Copy link
Member

Choose a reason for hiding this comment

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

and here

while c is (usually) the Allred-Rochow electronegativity. Values used
are *not* generated from pymatgen, and are found in
'okeeffe_params.json'.
Expand Down Expand Up @@ -2765,7 +2765,7 @@ def get_type(self, index):
raise ValueError("Index for getting order parameter type out-of-bounds!")
return self._types[index]

def get_parameters(self, index: int) -> list[float]:
def get_parameters(self, index):
"""Get list of floats that represents
the parameters associated
with calculation of the order
Expand All @@ -2774,10 +2774,12 @@ def get_parameters(self, index: int) -> list[float]:
inputted because of processing out of efficiency reasons.

Args:
index (int): index of order parameter for which to return associated params.
index (int):
index of order parameter for which associated parameters
are to be returned.

Returns:
list[float]: parameters of a given OP.
[float]: parameters of a given OP.
"""
if index < 0 or index >= len(self._types):
raise ValueError("Index for getting parameters associated with order parameter calculation out-of-bounds!")
Expand Down Expand Up @@ -4000,7 +4002,7 @@ def get_nn_data(self, structure: Structure, n: int, length=None):
nn_info.append(entry)
cn = len(nn_info)
cn_nninfo[cn] = nn_info
cn_weights[cn] = self._semicircle_integral(dist_bins, idx)
cn_weights[cn] = self._quadrant_integral(dist_bins, idx)

# add zero coord
cn0_weight = 1 - sum(cn_weights.values())
Expand Down Expand Up @@ -4057,10 +4059,13 @@ def get_cn_dict(self, structure: Structure, n: int, use_weights: bool = False, *
return super().get_cn_dict(structure, n, use_weights)

@staticmethod
def _semicircle_integral(dist_bins, idx):
def _semicircle_integral(dist_bins: list, idx: int) -> float:
"""
An internal method to get an integral between two bounds of a unit
semicircle. Used in algorithm to determine bond probabilities.
This function has an issue, which is detailed at
https://github.com/materialsproject/pymatgen/issues/3973.
Therefore, _quadrant_integral is now the method of choice.

Args:
dist_bins: (float) list of all possible bond weights
Expand All @@ -4085,6 +4090,35 @@ def _semicircle_integral(dist_bins, idx):

return (area1 - area2) / (0.25 * math.pi * radius**2)

@staticmethod
Copy link
Member

Choose a reason for hiding this comment

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

can you add a link to your issue to the _semicircle_integral doc string for later reference?

also, would be good to type-hint both the old and your new method

def _quadrant_integral(dist_bins: list, idx: int) -> float:
"""
An internal method to get an integral between two bounds of a unit
quadrant. Used in algorithm to determine bond probabilities.

Args:
dist_bins: (float) list of all possible bond weights
idx: (float) index of starting bond weight

Returns:
float: integral of portion of unit quadrant
"""
radius = 1

x1 = dist_bins[idx]
x2 = dist_bins[idx + 1]

areaquarter = 0.25 * math.pi * radius**2

area1 = areaquarter - 0.5 * (
radius**2 * math.acos(1 - x1 / radius) - (radius - x1) * math.sqrt(2 * radius * x1 - x1**2)
)
area2 = areaquarter - 0.5 * (
radius**2 * math.acos(1 - x2 / radius) - (radius - x2) * math.sqrt(2 * radius * x2 - x2**2)
)

return (area2 - area1) / areaquarter

@staticmethod
def transform_to_length(nn_data, length):
"""
Expand Down
12 changes: 8 additions & 4 deletions tests/analysis/test_local_env.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -445,6 +445,8 @@ def test_all_nn_classes(self):
assert voronoi_nn.get_cn(self.cscl, 0) == 8
assert voronoi_nn.get_cn(self.lifepo4, 0) == 6

assert CrystalNN._quadrant_integral([1, 0.36], 0) == approx(0.7551954297486029)
assert CrystalNN._quadrant_integral([1, 0.36, 0], 1) == approx(1 - 0.7551954297486029)
crystal_nn = CrystalNN()
assert crystal_nn.get_cn(self.diamond, 0) == 4
assert crystal_nn.get_cn(self.nacl, 0) == 6
Expand Down Expand Up @@ -1178,7 +1180,7 @@ def test_sanity(self):
def test_discrete_cn(self):
cnn = CrystalNN()
cn_array = []
expected_array = 8 * [6] + 20 * [4]
expected_array = 8 * [6] + 6 * [4] + [1] + 2 * [4] + [1] + 4 * [4] + [1] + +2 * [4] + [1] + 2 * [4]
for idx, _ in enumerate(self.lifepo4):
cn_array.append(cnn.get_cn(self.lifepo4, idx))

Expand All @@ -1202,6 +1204,7 @@ def test_weighted_cn(self):

def test_weighted_cn_no_oxid(self):
cnn = CrystalNN(weighted_cn=True)
cn_array = []
# fmt: off
expected_array = [
5.8962, 5.8996, 5.8962, 5.8996, 5.7195, 5.7195, 5.7202, 5.7194, 4.0012, 4.0012,
Expand All @@ -1210,7 +1213,8 @@ def test_weighted_cn_no_oxid(self):
]
# fmt: on
struct = self.lifepo4.copy().remove_oxidation_states()
cn_array = [cnn.get_cn(struct, idx, use_weights=True) for idx in range(len(struct))]
for idx in range(len(struct)):
cn_array.append(cnn.get_cn(struct, idx, use_weights=True))

assert_allclose(expected_array, cn_array, 2)

Expand All @@ -1222,11 +1226,11 @@ def test_fixed_length(self):

def test_cation_anion(self):
cnn = CrystalNN(weighted_cn=True, cation_anion=True)
assert cnn.get_cn(self.lifepo4, 0, use_weights=True) == approx(5.8630, abs=1e-2)
assert cnn.get_cn(self.lifepo4, 0, use_weights=True) == approx(5.5426, abs=1e-2)

def test_x_diff_weight(self):
cnn = CrystalNN(weighted_cn=True, x_diff_weight=0)
assert cnn.get_cn(self.lifepo4, 0, use_weights=True) == approx(5.8630, abs=1e-2)
assert cnn.get_cn(self.lifepo4, 0, use_weights=True) == approx(5.5426, abs=1e-2)

def test_noble_gas_material(self):
cnn = CrystalNN()
Expand Down