diff --git a/mdpath/mdpath.py b/mdpath/mdpath.py
index 74a5a07..1909877 100644
--- a/mdpath/mdpath.py
+++ b/mdpath/mdpath.py
@@ -21,7 +21,7 @@
 import pickle
 
 from mdpath.src.structure import StructureCalculations, DihedralAngles
-from mdpath.src.mutual_information import NMICalculator
+from mdpath.src.mutual_information import NMICalculator, binning_diagnostic
 from mdpath.src.graph import GraphBuilder
 from mdpath.src.cluster import PatwayClustering
 from mdpath.src.visualization import MDPathVisualize
@@ -216,6 +216,26 @@ def main():
     parser.add_argument("-null-seed", type=int, default=0)
     parser.add_argument("-ddof", type=int, choices=(0, 1), default=0)
     parser.add_argument("-null-outdir", default="null_analysis")
+    parser.add_argument(
+        "-hist-range",
+        dest="hist_range",
+        nargs="+",
+        default=["-180", "180"],
+        help="Fixed histogram support for the dihedral-movement histograms, in "
+        "degrees, as two floats (default: -180 180). Fixing the bin edges is what "
+        "makes NMI comparable across replicas, blocks and bootstrap samples. Pass "
+        '"none" to restore the pre-fix behaviour, where each input derived its own '
+        "support from its own extrema.",
+    )
+    parser.add_argument(
+        "-diagnose-bins",
+        dest="diagnose_bins",
+        action="store_true",
+        help="Write binning_diagnostic.csv: per-residue histogram occupancy under "
+        "the chosen support and branch-cut exposure (frames with |dphi| > 180 "
+        "before wrapping). Use it to check whether the fixed support is too wide "
+        "for the observed dihedral movements.",
+    )
     parser.add_argument(
         "-experimental-c-corr",
         dest="experimental_c_corr",
@@ -225,6 +245,15 @@ def main():
     )
 
     args = parser.parse_args()
+    # -hist-range: two floats, or the literal "none" to restore data-derived binning.
+    if len(args.hist_range) == 1 and str(args.hist_range[0]).lower() == "none":
+        hist_range = None
+    elif len(args.hist_range) == 2:
+        hist_range = (float(args.hist_range[0]), float(args.hist_range[1]))
+        if hist_range[0] >= hist_range[1]:
+            raise ValueError("-hist-range needs an increasing (low, high) pair")
+    else:
+        raise ValueError('-hist-range takes two floats or the literal "none"')
     if not args.topology or not args.trajectory:
         print("Both trajectory and topology files are required!")
         exit()
@@ -333,8 +362,16 @@ def main():
             )
     print("\033[1mTrajectory is processed and ready for analysis.\033[0m")
 
+    if args.diagnose_bins:
+        binning_diagnostic(
+            df_all_residues,
+            num_bins=35,
+            hist_range=hist_range,
+            csv_path="binning_diagnostic.csv",
+        )
+
     # Calculate the mutual information and build the graph
-    nmi_calc = NMICalculator(df_all_residues, invert=invert)
+    nmi_calc = NMICalculator(df_all_residues, invert=invert, hist_range=hist_range)
     nmi_calc.entropy_df.to_csv("entropy_df.csv", index=False)
     nmi_calc.nmi_df.to_csv("nmi_df.csv", index=False)
     graph_builder = GraphBuilder(
@@ -354,6 +391,7 @@ def main():
             per_replica_dfs,
             pairs=graph_pairs,
             ddof=args.ddof,
+            hist_range=hist_range,
         )
         edge_conf_df = edge_conf_calc.build_edge_confidence_df()
         edge_conf_df_graph = EdgeConfidenceCalculator.filter_to_graph_edges(
@@ -401,6 +439,7 @@ def main():
                 seed=args.null_seed,
                 ddof=args.ddof,
                 n_jobs=num_parallel_processes,
+                hist_range=hist_range,
             )
             null_outputs = null_calc.run(args.null_outdir, args.null_block_ns)
             status = null_outputs["edge_null_confidence"].corr_status.value_counts(
@@ -471,6 +510,7 @@ def main():
             topology,
             structure_calc.last_res_num,
             graphdist,
+            hist_range=hist_range,
         )
         file_name = "path_confidence_intervals.txt"
         bootstrap_analysis.bootstrap_write(file_name)
diff --git a/mdpath/src/bootstrap.py b/mdpath/src/bootstrap.py
index 3e30ac0..c7a652b 100644
--- a/mdpath/src/bootstrap.py
+++ b/mdpath/src/bootstrap.py
@@ -37,6 +37,8 @@ class BootstrapAnalysis:
 
         num_bins (int): Number of bins to group dihedral angle movements into for NMI calculation. Defaults to 35.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support forwarded to the NMI calculation. Each resample otherwise derives its own support, so path-confidence values inherit that inconsistency. Defaults to ``(-180.0, 180.0)``.
+
         common_counts (np.ndarray): Array with the counts of common paths between the original sample and bootstrap samples.
 
         path_confidence_intervals (dict): Dictionary with the confidence intervals for each path.
@@ -53,6 +55,7 @@ def __init__(
         last_residue: int,
         graphdist: int,
         num_bins: int = 35,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         self.df_all_residues = df_all_residues
         self.df_distant_residues = df_distant_residues
@@ -63,6 +66,7 @@ def __init__(
         self.last_residue = last_residue
         self.graphdist = graphdist
         self.num_bins = num_bins
+        self.hist_range = hist_range
         self.common_counts, self.path_confidence_intervals = self.bootstrap_analysis()
 
     def create_bootstrap_sample(self, df_dihedral: pd.DataFrame) -> tuple:
@@ -100,7 +104,9 @@ def process_bootstrap_sample(
             bootstrap_pathways (list): List of paths within the bootstrap sample.
         """
         bootstrap_sample = self.create_bootstrap_sample(self.df_all_residues)
-        nmi_calculator = NMICalculator(bootstrap_sample, num_bins=self.num_bins)
+        nmi_calculator = NMICalculator(
+            bootstrap_sample, num_bins=self.num_bins, hist_range=self.hist_range
+        )
         bootstrap_mi_diff = nmi_calculator.nmi_df
         graph = GraphBuilder(
             self.pdb, self.last_residue, bootstrap_mi_diff, self.graphdist
diff --git a/mdpath/src/confidence.py b/mdpath/src/confidence.py
index a94c4b4..f8584d0 100644
--- a/mdpath/src/confidence.py
+++ b/mdpath/src/confidence.py
@@ -11,8 +11,18 @@
 These scores measure *consistency of the estimated edge magnitude* across
 replicas; they are not calibrated probabilities, significance thresholds, or
 statements about effect size (an edge with consistently negligible mutual
-information can still score highly). The finite-sampling block reference and its
-experimental corrected score live in :mod:`mdpath.src.null_model`.
+information can still score highly).
+
+.. important::
+
+   Per-replica values are only comparable to each other under a **fixed histogram
+   support**. With ``hist_range=None`` every replica derives its own bin edges from
+   its own extrema, so a single wider excursion in one replica changes that
+   replica's binning and the difference propagates directly into ``MI_std`` ->
+   ``MI_cv`` -> ``confidence``. The default fixed support exists for this reason.
+
+The finite-sampling block reference and its experimental corrected score live in
+:mod:`mdpath.src.null_model`.
 
 Classes
 --------
@@ -43,6 +53,12 @@ class EdgeConfidenceCalculator:
 
         num_bins (int): Number of histogram bins forwarded to the mutual information calculator.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support forwarded to the
+            mutual information calculator. Per-replica values are only comparable under a
+            fixed support: with ``None`` each replica derives its own bin edges from its own
+            extrema, and that difference flows straight into ``MI_std`` -> ``MI_cv`` ->
+            ``confidence``. ``None`` is retained only to reproduce pre-fix numbers.
+
         pairs: Optional restriction of residue pairs forwarded to the mutual information calculator; ``None`` evaluates all pairs.
 
         ddof (int): Delta degrees of freedom (0 or 1) used for the standard deviation and the finite-replica floor.
@@ -56,6 +72,7 @@ def __init__(
         num_bins: int = 35,
         pairs=None,
         ddof: int = 0,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         if len(per_replica_dfs) < 2:
             raise ValueError("EdgeConfidenceCalculator requires at least two replicas.")
@@ -64,6 +81,7 @@ def __init__(
         self.per_replica_dfs = per_replica_dfs
         self.num_bins = num_bins
         self.pairs, self.ddof = pairs, ddof
+        self.hist_range = hist_range
         self.per_replica_nmi_dfs = self._compute_per_replica_nmi()
 
     def _compute_per_replica_nmi(self) -> List[pd.DataFrame]:
@@ -78,6 +96,7 @@ def _compute_per_replica_nmi(self) -> List[pd.DataFrame]:
                 num_bins=self.num_bins,
                 pairs=self.pairs,
                 show_progress=False,
+                hist_range=self.hist_range,
             ).nmi_df
             for frame in self.per_replica_dfs
         ]
diff --git a/mdpath/src/mutual_information.py b/mdpath/src/mutual_information.py
index 97d2861..07430ce 100644
--- a/mdpath/src/mutual_information.py
+++ b/mdpath/src/mutual_information.py
@@ -9,6 +9,11 @@
 --------
 
 :class:`NMICalculator`
+
+Functions
+---------
+
+:func:`binning_diagnostic`
 """
 
 import pandas as pd
@@ -29,6 +34,14 @@ class NMICalculator:
 
         num_bins (int): Number of bins to use for histogram calculations. Default is 35.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` support for the
+            dihedral-movement histograms, in degrees. Default ``(-180, 180)``
+            fixes the bin edges so that values computed from different
+            replicas, blocks or bootstrap samples share one binning and are
+            directly comparable. Passing ``None`` restores the previous
+            behaviour, in which each input derived its own support from its
+            own extrema.
+
         pairs (list | None): Optional list of residue pairs restricting which pairs
             are evaluated (each pair is any two residue labels, e.g. ``42`` or
             ``"Res 42"``). When ``None`` every pair is computed as before. This only
@@ -48,12 +61,14 @@ def __init__(
         invert=False,
         pairs=None,
         show_progress: bool = True,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         self.df_all_residues = df_all_residues
         self.num_bins = num_bins
         self.invert = invert
         self.pairs = pairs
         self.show_progress = show_progress
+        self.hist_range = hist_range
         self.nmi_df, self.entropy_df = self.NMI_calcs()
 
     def _column_pairs(self, columns):
@@ -88,7 +103,11 @@ def NMI_calcs(self):
         histograms = {}
         entropys = {}
         for col in columns:
-            hist, _ = np.histogram(self.df_all_residues[col], bins=self.num_bins)
+            hist, _ = np.histogram(
+                self.df_all_residues[col],
+                bins=self.num_bins,
+                range=self.hist_range,
+            )
             histograms[col] = hist
             entropys[col] = entropy(hist)
 
@@ -103,6 +122,14 @@ def NMI_calcs(self):
                 self.df_all_residues[col1],
                 self.df_all_residues[col2],
                 bins=self.num_bins,
+                # np.histogram wants (lo, hi); np.histogram2d wants
+                # [(lo, hi), (lo, hi)] and raises otherwise. Both accept None,
+                # so the None branch is a true no-op.
+                range=(
+                    None
+                    if self.hist_range is None
+                    else [self.hist_range, self.hist_range]
+                ),
             )
             mi = mutual_info_score(
                 histograms[col1], histograms[col2], contingency=hist_joint
@@ -121,3 +148,77 @@ def NMI_calcs(self):
             max_nmi_diff = nmi_df["MI Difference"].max()
             nmi_df["MI Difference"] = max_nmi_diff - nmi_df["MI Difference"]
         return nmi_df, entropy_df
+
+
+def binning_diagnostic(
+    df_all_residues: pd.DataFrame,
+    num_bins: int = 35,
+    hist_range=(-180.0, 180.0),
+    csv_path: str = "binning_diagnostic.csv",
+) -> pd.DataFrame:
+    """Report histogram occupancy and branch-cut exposure per residue.
+
+    A fixed support makes mutual information comparable across replicas, blocks and
+    bootstrap samples, but it costs resolution when the dihedral movements are narrow:
+    bins outside the sampled range sit empty. This quantifies that trade-off so the
+    width of the support can be chosen from data rather than guessed.
+
+    Args:
+        df_all_residues (pd.DataFrame): Dihedral angle movements, one column per residue.
+
+        num_bins (int): Number of histogram bins. Default is 35.
+
+        hist_range (tuple | None): Fixed ``(low, high)`` support to evaluate. ``None``
+            reports the data-derived binning instead.
+
+        csv_path (str): Where the per-residue table is written.
+
+    Returns:
+        pd.DataFrame: One row per residue with ``residue``, ``occupied_bins``,
+        ``occupied_fraction``, ``movement_min``, ``movement_max`` and, when the
+        upstream calculation recorded it, ``branch_cut_crossings`` (the number of
+        frames whose *unwrapped* difference exceeded 180 degrees).
+    """
+    crossings = df_all_residues.attrs.get("branch_cut_crossings", {})
+    records = []
+    for column in df_all_residues.columns:
+        values = np.asarray(df_all_residues[column], dtype=float)
+        values = values[np.isfinite(values)]
+        counts, _ = np.histogram(values, bins=num_bins, range=hist_range)
+        record = {
+            "residue": column,
+            "occupied_bins": int(np.count_nonzero(counts)),
+            "occupied_fraction": float(np.count_nonzero(counts) / num_bins),
+            "movement_min": float(values.min()) if values.size else np.nan,
+            "movement_max": float(values.max()) if values.size else np.nan,
+        }
+        try:
+            record["branch_cut_crossings"] = int(crossings[residue_id(column)])
+        except (KeyError, ValueError):
+            record["branch_cut_crossings"] = np.nan
+        records.append(record)
+
+    frame = pd.DataFrame(records)
+    if csv_path:
+        frame.to_csv(csv_path, index=False)
+
+    occupied = frame["occupied_bins"].to_numpy()
+    low, high = np.percentile(occupied, [5, 95]) if occupied.size else (np.nan, np.nan)
+    sparse = int(np.count_nonzero(occupied < 10))
+    print(
+        f"Binning diagnostic ({num_bins} bins, support "
+        f"{'data-derived' if hist_range is None else f'{hist_range[0]:g}..{hist_range[1]:g}'}): "
+        f"occupied bins median {np.median(occupied):.0f}, 5th/95th pct {low:.0f}/{high:.0f}; "
+        f"{sparse}/{len(occupied)} residues ({sparse / max(1, len(occupied)):.1%}) occupy <10 bins."
+    )
+    exposure = frame["branch_cut_crossings"]
+    if exposure.notna().any():
+        affected = int((exposure.fillna(0) > 0).sum())
+        print(
+            f"  branch-cut exposure: {affected}/{len(exposure)} residues "
+            f"({affected / max(1, len(exposure)):.1%}) had at least one frame with "
+            f"|dphi| > 180 before wrapping; {int(exposure.fillna(0).sum())} frames total."
+        )
+    if csv_path:
+        print(f"  per-residue values written to {csv_path}")
+    return frame
diff --git a/mdpath/src/null_model.py b/mdpath/src/null_model.py
index d789c65..f3c0ef4 100644
--- a/mdpath/src/null_model.py
+++ b/mdpath/src/null_model.py
@@ -14,6 +14,11 @@
 score, so its outputs are best read as a sampling-stability reference rather than
 a formal null distribution.
 
+Per-block and matched-observed values are likewise only comparable under a **fixed
+histogram support** (``hist_range``); with ``None`` each block derives its own bin
+edges and the block-to-block difference in binning is absorbed into the reference as
+if it were sampling noise.
+
 The design sweeps block length ``L`` (ns) and replica count ``R`` so every edge
 is compared at the *same* ``(L, R)`` as the matched observed estimate. For each
 source and ``(L, R)`` the reference is estimated by **resampling** ``draws``
@@ -71,6 +76,7 @@ def _edge_vector(
     frame: pd.DataFrame,
     pairs: Sequence[Pair],
     num_bins: int,
+    hist_range=(-180.0, 180.0),
 ) -> np.ndarray:
     """Computes the per-edge mutual-information difference for a single trajectory block.
 
@@ -85,6 +91,10 @@ def _edge_vector(
 
         num_bins (int): Number of histogram bins for the mutual information calculator.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support. Blocks are
+            only comparable to each other, and to the matched observed windows, under a
+            shared support; ``None`` restores per-block data-derived binning.
+
     Returns:
         np.ndarray: One MI-difference value per requested pair, ordered like ``pairs``.
     """
@@ -93,6 +103,7 @@ def _edge_vector(
         num_bins=num_bins,
         pairs=pairs,
         show_progress=False,
+        hist_range=hist_range,
     ).nmi_df
 
     lookup = {
@@ -196,8 +207,8 @@ def _source_task(payload: tuple) -> Tuple[pd.DataFrame, List[dict]]:
 
     Args:
         payload (tuple): Packed arguments ``(source_id, source, lengths, r_values,
-            pairs, num_bins, ddof, draws, seed)`` so the task can run inside a
-            worker process.
+            pairs, num_bins, ddof, draws, seed, hist_range)`` so the task can run
+            inside a worker process.
 
     Returns:
         Tuple[pd.DataFrame, List[dict]]: The per-edge reference table (one row per
@@ -216,6 +227,7 @@ def _source_task(payload: tuple) -> Tuple[pd.DataFrame, List[dict]]:
         ddof,
         draws,
         seed,
+        hist_range,
     ) = payload
     n_edges = len(pairs)
     res_u = np.fromiter((pair[0] for pair in pairs), dtype=int, count=n_edges)
@@ -230,7 +242,9 @@ def _source_task(payload: tuple) -> Tuple[pd.DataFrame, List[dict]]:
             )
             for index in range(n_available)
         ]
-        weights = np.vstack([_edge_vector(block, pairs, num_bins) for block in blocks])
+        weights = np.vstack(
+            [_edge_vector(block, pairs, num_bins, hist_range) for block in blocks]
+        )
         for r_null in r_values[block_ns]:
             groups = _block_groups(
                 n_available, r_null, draws, seed, source_id, block_frames
@@ -348,14 +362,14 @@ def _observed_task(payload: tuple) -> Dict[float, np.ndarray]:
 
     Args:
         payload (tuple): Packed arguments ``(frame, lengths, pairs, num_bins,
-            observed_window)`` so the task can run inside a worker process.
+            observed_window, hist_range)`` so the task can run inside a worker process.
 
     Returns:
         Dict[float, np.ndarray]: Mapping from block length in nanoseconds to a
         ``(n_windows, n_edges)`` weight matrix for this replica (a single row in
         ``"leading"`` mode).
     """
-    frame, lengths, pairs, num_bins, observed_window = payload
+    frame, lengths, pairs, num_bins, observed_window, hist_range = payload
     result: Dict[float, np.ndarray] = {}
     for block_ns, block_frames in lengths:
         n_windows = 1 if observed_window == "leading" else len(frame) // block_frames
@@ -426,6 +440,11 @@ class BlockSplitNull:
 
         num_bins (int): Number of histogram bins passed to the MI calculator.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support passed to the
+            MI calculator. Per-block and matched-observed values are only comparable under
+            a shared support; ``None`` restores per-input data-derived binning and is
+            retained only to reproduce pre-fix numbers.
+
         draws (int): Number of random resampling draws for both the block reference and the matched observed estimate.
 
         seed (int): Seed for the resampling random generators.
@@ -455,6 +474,7 @@ def __init__(
         ddof: int = 0,
         n_jobs: int = 1,
         observed_window: str = "random",
+        hist_range=(-180.0, 180.0),
     ) -> None:
         if len(source_dfs) < 1 or len(observed_dfs) < 2:
             raise ValueError(
@@ -493,6 +513,7 @@ def __init__(
         self.draws, self.seed, self.ddof = int(draws), int(seed), ddof
         self.n_jobs = max(1, int(n_jobs))
         self.observed_window = observed_window
+        self.hist_range = hist_range
 
         # Feasibility is evaluated per length: a length is usable only if every
         # source yields at least two non-overlapping blocks and every observed
@@ -693,6 +714,7 @@ def run(
                 self.ddof,
                 self.draws,
                 self.seed,
+                self.hist_range,
             )
             for source_id, source in enumerate(self.sources)
         ]
@@ -720,7 +742,14 @@ def run(
         ).reset_index()
 
         observed_payloads = [
-            (frame, self.lengths, self.pairs, self.num_bins, self.observed_window)
+            (
+                frame,
+                self.lengths,
+                self.pairs,
+                self.num_bins,
+                self.observed_window,
+                self.hist_range,
+            )
             for frame in self.observed
         ]
         observed_results = self._map(_observed_task, observed_payloads)
@@ -1024,6 +1053,7 @@ def _write_metadata(self, outdir, grid, native, fixed_r, primary_l, primary_r):
             "ddof": self.ddof,
             "n_jobs": self.n_jobs,
             "observed_window": self.observed_window,
+            "hist_range": list(self.hist_range) if self.hist_range else None,
             "block_lengths": [
                 {
                     "L_ns": block_ns,
diff --git a/mdpath/src/structure.py b/mdpath/src/structure.py
index 6955c51..6c13463 100644
--- a/mdpath/src/structure.py
+++ b/mdpath/src/structure.py
@@ -225,7 +225,13 @@ def calculate_dihedral_movements_multi_traj(
                 f"dihedral data across replicas.\033[0m"
             )
 
-        all_dfs = [df[common_cols] for df in all_dfs]
+        aligned = []
+        for df in all_dfs:
+            attrs = df.attrs.copy()
+            df = df[common_cols].copy()
+            df.attrs = attrs
+            aligned.append(df)
+        all_dfs = aligned
         concatenated = pd.concat(all_dfs, axis=0, ignore_index=True)
         if return_per_replica:
             return concatenated, all_dfs
@@ -258,14 +264,17 @@ def __init__(
         self.last_res_num = last_res_num
         self.num_residues = num_residues
 
-    def calc_dihedral_angle_movement(self, res_id: int) -> tuple:
-        """Calculates dihedral angle movement for a residue over the course of the MD trajectory.
+    def _calc_dihedral_movement_with_crossings(self, res_id: int) -> tuple:
+        """Calculate wrapped dihedral movement and count branch-cut crossings.
 
         Args:
             res_id (int): Residue number.
 
         Returns:
-            tuple[int, np.ndarray] | None: Tuple of (residue_id, dihedral_angles) if successful, None if failed.
+            tuple[int, np.ndarray, int] | None: ``(residue_id, movement, crossings)``
+            where ``crossings`` is the number of frames whose *unwrapped* difference
+            exceeded 180 degrees, i.e. the residue's exposure to the periodic branch
+            cut. ``None`` if the dihedral could not be calculated.
         """
         try:
             res = self.traj.residues[res_id - self.first_res_num]
@@ -274,14 +283,34 @@ def calc_dihedral_angle_movement(self, res_id: int) -> tuple:
                 return None
             R = Dihedral(ags).run()
             dihedrals = R.results.angles
-            dihedral_angle_movement = np.diff(dihedrals, axis=0)
-            return res_id, dihedral_angle_movement
+            raw = np.diff(dihedrals, axis=0)
+            # phi is periodic on [-180, 180); an unwrapped difference turns a
+            # small step across the branch cut into a spurious ~360 deg jump,
+            # which then dominates the histogram support for that residue.
+            dihedral_angle_movement = (raw + 180.0) % 360.0 - 180.0
+            crossings = int(np.count_nonzero(np.abs(raw) > 180.0))
+            return res_id, dihedral_angle_movement, crossings
         except (TypeError, AttributeError, IndexError) as e:
             logging.debug(
                 f"Failed to calculate dihedral for residue {res_id}: {str(e)}"
             )
             return None
 
+    def calc_dihedral_angle_movement(self, res_id: int) -> tuple:
+        """Calculates dihedral angle movement for a residue over the course of the MD trajectory.
+
+        The difference is wrapped onto ``[-180, 180)`` because phi is periodic; see
+        :meth:`_calc_dihedral_movement_with_crossings`.
+
+        Args:
+            res_id (int): Residue number.
+
+        Returns:
+            tuple[int, np.ndarray] | None: Tuple of (residue_id, dihedral_angles) if successful, None if failed.
+        """
+        result = self._calc_dihedral_movement_with_crossings(res_id)
+        return None if result is None else result[:2]
+
     def calculate_dihedral_movement_parallel(
         self,
         num_parallel_processes: int,
@@ -295,6 +324,7 @@ def calculate_dihedral_movement_parallel(
             pd.DataFrame: DataFrame with all residue dihedral angle movements.
         """
         collected = []
+        branch_cut_crossings = {}
 
         try:
             with Pool(processes=num_parallel_processes) as pool:
@@ -304,7 +334,7 @@ def calculate_dihedral_movement_parallel(
                     desc="\033[1mProcessing residue dihedral movements\033[0m",
                 ) as pbar:
                     results = pool.imap_unordered(
-                        self.calc_dihedral_angle_movement,
+                        self._calc_dihedral_movement_with_crossings,
                         range(self.first_res_num, self.last_res_num + 1),
                     )
 
@@ -313,7 +343,8 @@ def calculate_dihedral_movement_parallel(
                             pbar.update(1)
                             continue
 
-                        res_id, dihedral_data = result
+                        res_id, dihedral_data, crossings = result
+                        branch_cut_crossings[res_id] = crossings
                         try:
                             collected.append(
                                 pd.DataFrame(dihedral_data, columns=[f"Res {res_id}"])
@@ -329,5 +360,8 @@ def calculate_dihedral_movement_parallel(
             logging.error(f"Parallel processing failed: {str(e)}")
 
         if collected:
-            return pd.concat(collected, axis=1)
+            result = pd.concat(collected, axis=1)
+            # branch-cut exposure per residue, for the binning diagnostic
+            result.attrs["branch_cut_crossings"] = branch_cut_crossings
+            return result
         return pd.DataFrame()
diff --git a/mdpath/tests/test_mutual_information.py b/mdpath/tests/test_mutual_information.py
index 56bf668..de25bdc 100644
--- a/mdpath/tests/test_mutual_information.py
+++ b/mdpath/tests/test_mutual_information.py
@@ -41,11 +41,13 @@ def test_nmi_calculator_invert():
         result["MI Difference"] >= 0
     ).all(), "MI Difference values should be non-negative"
 
-    # Add your expected values here
+    # Recomputed under the fixed (-180, 180) histogram support introduced with the
+    # circular-coordinate fix; the pre-fix values are still reproduced exactly by
+    # hist_range=None (see test_hist_range_none_is_exact_no_op).
     expected_values = {
-        ("Residue1", "Residue2"): 1.1102230246251565e-16,
-        ("Residue1", "Residue3"): 0.006253825923384082,
-        ("Residue2", "Residue3"): 0.006287366618772272,
+        ("Residue1", "Residue2"): 0.019723970154380432,
+        ("Residue1", "Residue3"): 0.0,
+        ("Residue2", "Residue3"): 0.0005904872893245372,
     }
     for pair, expected in expected_values.items():
         actual = result.loc[result["Residue Pair"] == pair, "MI Difference"].values[0]
@@ -89,14 +91,89 @@ def test_nmi_calculator():
         result["MI Difference"] >= 0
     ).all(), "MI Difference values should be non-negative"
 
-    # Add your expected values here
+    # Recomputed under the fixed (-180, 180) histogram support introduced with the
+    # circular-coordinate fix; the pre-fix values are still reproduced exactly by
+    # hist_range=None (see test_hist_range_none_is_exact_no_op).
     expected_values = {
-        ("Residue1", "Residue2"): 0.6729976399990001,
-        ("Residue1", "Residue3"): 0.6667438140756161,
-        ("Residue2", "Residue3"): 0.6667102733802279,
+        ("Residue1", "Residue2"): 0.6510110199334083,
+        ("Residue1", "Residue3"): 0.6707349900877887,
+        ("Residue2", "Residue3"): 0.6701445027984642,
     }
     for pair, expected in expected_values.items():
         actual = result.loc[result["Residue Pair"] == pair, "MI Difference"].values[0]
         assert np.isclose(
             actual, expected, atol=1e-5
         ), f"For pair {pair}, expected {expected} but got {actual}"
+
+
+# --------------------------------------------------------------------------- #
+# circular-coordinate / fixed-support regression tests                         #
+# --------------------------------------------------------------------------- #
+def test_dihedral_difference_is_wrapped():
+    """A step across the +-180 branch cut must stay a small step."""
+    from mdpath.src.structure import DihedralAngles
+
+    angles = np.array([[179.0], [-179.0], [178.0]])
+    raw = np.diff(angles, axis=0)
+    wrapped = (raw + 180.0) % 360.0 - 180.0
+
+    # unwrapped, a 2 deg step becomes a ~358 deg jump
+    assert np.allclose(raw.ravel(), [-358.0, 357.0])
+    # 179 -> -179 is +2 deg (179 + 2 = 181 == -179); -179 -> 178 is -3 deg
+    assert np.allclose(wrapped.ravel(), [2.0, -3.0])
+    assert np.all(np.abs(wrapped) <= 180.0)
+    # and this is the code path the package uses
+    assert hasattr(DihedralAngles, "_calc_dihedral_movement_with_crossings")
+
+
+def test_fixed_support_gives_identical_bin_edges_across_inputs():
+    """Inputs with different extrema must share bin edges under a fixed support."""
+    rng = np.random.default_rng(0)
+    narrow = pd.DataFrame({"Res 1": rng.uniform(-20, 20, 200)})
+    wide = pd.DataFrame({"Res 1": rng.uniform(-150, 150, 200)})
+
+    edges_fixed = [
+        np.histogram(f["Res 1"], bins=35, range=(-180.0, 180.0))[1]
+        for f in (narrow, wide)
+    ]
+    assert np.array_equal(edges_fixed[0], edges_fixed[1])
+
+    edges_legacy = [np.histogram(f["Res 1"], bins=35)[1] for f in (narrow, wide)]
+    assert not np.array_equal(edges_legacy[0], edges_legacy[1])
+
+
+def test_hist_range_none_is_exact_no_op():
+    """hist_range=None must reproduce the pre-fix values bit-for-bit."""
+    np.random.seed(0)
+    df = pd.DataFrame(
+        {f"Residue{i + 1}": np.random.uniform(-180, 180, size=100) for i in range(3)}
+    )
+    legacy = NMICalculator(df, show_progress=False, hist_range=None).nmi_df
+    pre_fix = {
+        ("Residue1", "Residue2"): 0.6729976399990001,
+        ("Residue1", "Residue3"): 0.6667438140756161,
+        ("Residue2", "Residue3"): 0.6667102733802279,
+    }
+    for pair, expected in pre_fix.items():
+        actual = legacy.loc[legacy["Residue Pair"] == pair, "MI Difference"].values[0]
+        assert actual == expected, f"{pair}: {actual!r} != {expected!r}"
+
+
+def test_one_contaminated_frame_does_not_collapse_the_histogram():
+    """A single branch-cut frame must not distort a residue once wrapped."""
+    rng = np.random.default_rng(1)
+    clean = rng.uniform(-15, 15, 1000)
+
+    contaminated = clean.copy()
+    contaminated[0] = -358.0  # what an unwrapped cut looks like
+    wrapped = (contaminated + 180.0) % 360.0 - 180.0
+
+    def occupied(values, **kwargs):
+        return np.count_nonzero(np.histogram(values, bins=35, **kwargs)[0])
+
+    # fixed support: wrapping restores the clean occupancy
+    assert occupied(wrapped, range=(-180.0, 180.0)) == occupied(
+        clean, range=(-180.0, 180.0)
+    )
+    # unwrapped + data-derived support is what collapses it
+    assert occupied(contaminated) < occupied(clean) / 3
diff --git a/mdpath/tests/test_null_model.py b/mdpath/tests/test_null_model.py
index 36a1c8a..8e12471 100644
--- a/mdpath/tests/test_null_model.py
+++ b/mdpath/tests/test_null_model.py
@@ -181,8 +181,8 @@ def test_random_observed_windows_use_the_whole_replica(tmp_path):
     # "leading" evaluates one window per replica; "random" evaluates every window
     reps = _windowed_replicas()
     lengths = [(30.0, 30)]
-    lead = _observed_task((reps[0], lengths, [(1, 2)], 8, "leading"))
-    rand = _observed_task((reps[0], lengths, [(1, 2)], 8, "random"))
+    lead = _observed_task((reps[0], lengths, [(1, 2)], 8, "leading", (-180.0, 180.0)))
+    rand = _observed_task((reps[0], lengths, [(1, 2)], 8, "random", (-180.0, 180.0)))
     assert lead[30.0].shape[0] == 1
     assert rand[30.0].shape[0] == len(reps[0]) // 30 == 4
 
diff --git a/mdpath/tests/test_structure.py b/mdpath/tests/test_structure.py
index 8c4bc38..dbf060f 100644
--- a/mdpath/tests/test_structure.py
+++ b/mdpath/tests/test_structure.py
@@ -96,12 +96,14 @@ def test_calc_dihedral_angle_movement(mock_dihedral_class):
 @patch("mdpath.src.structure.Pool")  # Mocking multiprocessing Pool
 def test_calculate_dihedral_movement_parallel(mock_pool, mock_traj):
     angles = DihedralAngles(mock_traj, 1, 2, 1)
+    # workers now yield (res_id, movement, branch_cut_crossings)
     mock_pool.return_value.__enter__.return_value.imap_unordered.return_value = [
-        (1, np.array([1, 2]))
+        (1, np.array([1, 2]), 0)
     ]
     df = angles.calculate_dihedral_movement_parallel(num_parallel_processes=1)
     assert isinstance(df, pd.DataFrame)
     assert not df.empty
+    assert df.attrs["branch_cut_crossings"] == {1: 0}
 
 
 class TestCalculateDihedralMovementsMultiTraj: