Fixed failing tests, added link to the issue that describes problem with old semicircle function, and type-hinted functions

Author: Nils Edvin Richard Zimmermann

src/pymatgen/analysis/local_env.py

@@ -19,13 +19,12 @@
 import numpy as np
 from monty.dev import deprecated, requires
 from monty.serialization import loadfn
-from ruamel.yaml import YAML
-from scipy.spatial import Voronoi
-
 from pymatgen.analysis.bond_valence import BV_PARAMS, BVAnalyzer
 from pymatgen.analysis.graphs import MoleculeGraph, StructureGraph
 from pymatgen.analysis.molecule_structure_comparator import CovalentRadius
 from pymatgen.core import Element, IStructure, PeriodicNeighbor, PeriodicSite, Site, Species, Structure
+from ruamel.yaml import YAML
+from scipy.spatial import Voronoi
 
 try:
     from openbabel import openbabel
@@ -35,10 +34,9 @@
 if TYPE_CHECKING:
     from typing import Any
 
-    from typing_extensions import Self
-
     from pymatgen.core.composition import SpeciesLike
     from pymatgen.util.typing import Tuple3Ints
+    from typing_extensions import Self
 
 
 __author__ = "Shyue Ping Ong, Geoffroy Hautier, Sai Jayaraman, "
@@ -1161,7 +1159,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)
@@ -1218,7 +1216,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:
@@ -1984,7 +1982,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'.
@@ -2755,7 +2753,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
@@ -2764,10 +2762,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!")
@@ -3990,7 +3990,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())
@@ -4047,10 +4047,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
@@ -4075,6 +4078,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):
         """

tests/analysis/test_local_env.py

@@ -7,8 +7,6 @@
 import numpy as np
 import pytest
 from numpy.testing import assert_allclose
-from pytest import approx
-
 from pymatgen.analysis.graphs import MoleculeGraph, StructureGraph
 from pymatgen.analysis.local_env import (
     BrunnerNNReal,
@@ -38,6 +36,7 @@
 )
 from pymatgen.core import Element, Lattice, Molecule, Structure
 from pymatgen.util.testing import TEST_FILES_DIR, PymatgenTest
+from pytest import approx
 
 TEST_DIR = f"{TEST_FILES_DIR}/analysis/local_env/fragmenter_files"
 
@@ -445,6 +444,10 @@ 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
@@ -1178,7 +1181,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))
 
@@ -1202,6 +1205,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,
@@ -1210,7 +1214,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)
 
@@ -1222,11 +1227,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()