diff --git a/src/pymatgen/analysis/local_env.py b/src/pymatgen/analysis/local_env.py index 18fb7941dba..e93fef3d685 100644 --- a/src/pymatgen/analysis/local_env.py +++ b/src/pymatgen/analysis/local_env.py @@ -1185,7 +1185,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 """ elems = _get_elements(site) return all(elem in targets for elem in elems) @@ -1240,7 +1240,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) @property def structures_allowed(self) -> bool: @@ -2016,7 +2016,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, while c is (usually) the Allred-Rochow electronegativity. Values used are *not* generated from pymatgen, and are found in 'okeeffe_params.json'. @@ -2818,7 +2818,9 @@ def get_parameters(self, index: int) -> dict[str | int, float] | None: 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: dict[str, float]: parameters of a given OP. @@ -4075,7 +4077,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()) @@ -4132,10 +4134,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 @@ -4160,6 +4165,35 @@ def _semicircle_integral(dist_bins, idx): return (area1 - area2) / (0.25 * math.pi * radius**2) + @staticmethod + 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): """ diff --git a/tests/analysis/test_local_env.py b/tests/analysis/test_local_env.py index 2c5c4d6039b..b84701a713f 100644 --- a/tests/analysis/test_local_env.py +++ b/tests/analysis/test_local_env.py @@ -470,6 +470,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 @@ -1211,7 +1213,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)) @@ -1235,6 +1237,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, @@ -1243,7 +1246,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) @@ -1255,11 +1259,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()