diff --git a/cmeutils/structure.py b/cmeutils/structure.py index 4ee10f7..b63f598 100644 --- a/cmeutils/structure.py +++ b/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 diff --git a/cmeutils/tests/assets/AB-traj.gsd b/cmeutils/tests/assets/AB-traj.gsd index 58755f5..1774fb5 100644 Binary files a/cmeutils/tests/assets/AB-traj.gsd and b/cmeutils/tests/assets/AB-traj.gsd differ diff --git a/cmeutils/tests/test_structure.py b/cmeutils/tests/test_structure.py index 9c11abb..34edc6a 100644 --- a/cmeutils/tests/test_structure.py +++ b/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(