Added _quadrant_integral function in CrystalNN and switched over to use that function instead of _semicircle_integral which did not compute circular segments of unit quadrants correctly

Author: Nils Edvin Richard Zimmermann

src/pymatgen/analysis/local_env.py

@@ -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())
@@ -4075,6 +4075,35 @@ def _semicircle_integral(dist_bins, idx):
 
         return (area1 - area2) / (0.25 * math.pi * radius**2)
 
+    @staticmethod
+    def _quadrant_integral(dist_bins, idx):
+        """
+        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

@@ -444,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