Fix issue in counting number of excluded pairs in RDF.

cmeutils/structure.py

@@ -502,13 +502,16 @@ def gsd_rdf(
             # These 2 *_indices variables store global particle indices (Before filtering by type)
             # If excluding by bond depth, these need to be passed into filter_nlist
             type_A_indices = np.where(type_A)[0]
+            type_A_set = set(type_A_indices)
             type_B_indices = np.where(type_B)[0]
+            type_B_set = set(type_B_indices)
             exclude_ii = A_name == B_name
         else:  # Use all particles for this RDF
             type_A = type_B = np.ones(
                 snap.particles.N, dtype=bool
             )  # Array of True at all indices
             type_A_indices = type_B_indices = np.arange(snap.particles.N)
+            type_A_set = type_B_set = type_A_indices
             exclude_ii = True
 
         # Build up pair exclusions if exclude_bond_depth or exclude_all_bonded
@@ -528,7 +531,14 @@ def gsd_rdf(
                 [i * max_idx + j for i, j in excluded_pairs]
             )
 
-            n_excluded = len(excluded_pairs)
+            # Only count excluded pairs if they are being used in the RDF calculation
+            n_excluded = sum(
+                1
+                for i, j in excluded_pairs
+                if (i in type_A_set and j in type_B_set)
+                or (j in type_A_set and i in type_B_set)
+            )
+
             if A_name == B_name or not any(
                 [A_name, B_name]
             ):  # Using same type or all particles

cmeutils/tests/test_structure.py

@@ -365,6 +365,32 @@ def test_gsd_rdf_exclusions(self, AB_chain_gsd):
         rdf_values = rdf.rdf[zero_indices]
         assert np.allclose(rdf_values, zero_array)
 
+    def test_gsd_rdf_normalization(self, AB_chain_gsd):
+        # Chain has 10 particles total
+        # All pairs are being used, N total = (10 * 9) / 2 = 45
+        # Exclude bond depth of 1 should exclude all bonded pairs = 9
+        # 45 / (45 - 9) = 1.25
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_bond_depth=1,
+            exclude_all_bonded=False,
+        )
+        assert np.isclose(scale_factor, 1.25)
+
+        # Excl bond depth of 2 should result in 9 (bonds) + 8 (angles) exclusions
+        # 45 / (45 - 17) = 1.607
+
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_bond_depth=2,
+            exclude_all_bonded=False,
+        )
+        assert np.isclose(scale_factor, 1.6071428571428572)
+
     def test_gsd_rdf_r_max(self, LJ_gsd):
         """Test 2 RDFs with different r_cuts. The values of the shared r_cut region should be very close"""
         rdf, scale_factor = gsd_rdf(