From 58f9a53cdbcee94b395ff73a3691c3aff8690734 Mon Sep 17 00:00:00 2001 From: srikar jy Date: Sat, 15 Aug 2026 16:59:39 -0400 Subject: [PATCH 1/2] Fix PersistenceLength normalization for sliced run (Issue #5453) PersistenceLength._conclude() normalized bond_autocorrelation using the total number of frames in the trajectory instead of the number of frames actually analyzed. This gave wrong results whenever run() was called with start, stop, or step. Now uses self.n_frames, which AnalysisBase already tracks for this purpose. Added a regression test comparing full trajectory and sliced runs against manually computed normalization values. --- package/AUTHORS | 1 + package/CHANGELOG | 37 +++++++++++++- package/MDAnalysis/analysis/polymer.py | 2 +- .../analysis/test_persistencelength.py | 48 +++++++++++++++++++ 4 files changed, 86 insertions(+), 2 deletions(-) diff --git a/package/AUTHORS b/package/AUTHORS index ed82321a9e8..3ed92cc27ba 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -284,6 +284,7 @@ Chronological list of authors - Sai Udayagiri - Apoorva Verma - Aryaman Chaudhri + - Srikar JY External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index 3e7c3d309d3..d0fb038f75b 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -18,11 +18,12 @@ The rules for this file: spyke7, talagayev, tanii1125, BradyAJohnston, hejamu, jeremyleung521, harshitgajjela-droid, kunjsinha, aygarwal, jauy123, Dreamstick9, ollyfutur, Amarendra22, charity-g, ParthUppal523, apoorva-01, RMeli, - raulloiscuns, Aryaman-Chaudhri + raulloiscuns, Aryaman-Chaudhri, srikarjy * 2.11.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * `AtomGroup.rotate()` and the `rotateby` trajectory transformation now also rotate velocities and forces besides positions. This also affects `MDAnalysis.analysis.align.alignto()` and `AlignTraj`, since they @@ -120,6 +121,7 @@ Deprecations * 2.10.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix incorrect conversion factors for speed units A/fs, A/us, A/ms in MDAnalysis.units (Issue #5051, PR #5053) * `analysis.polymer.sort_backbone` is now working for discontinuous polymers @@ -212,6 +214,7 @@ Deprecations * 2.9.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Add support for TPR files produced by GROMACS 2024.4 and GROMACS 2025.0 * Fixes invalid default unit from Angstrom to Angstrom^{-3} for convert_density() function. (Issue #4829) @@ -254,6 +257,7 @@ Changes * 2.8.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix `test_imd.py` test failures after imdclient 0.2.4 increased default timeout by ensuring EOF sent (Issue #5442, PR #5443) * Allows bond/angle/dihedral connectivity to be guessed additively with @@ -389,6 +393,7 @@ Deprecations * 2.7.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix doctest errors of lib/transformations.py (Issue #3925, PR #4370) * Fix documentation building errors due to `html_static_path` and unindented block of code (Issue #4362, PR #4365) @@ -458,6 +463,7 @@ Deprecations * 2.6.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix Atom type guessing error (PR #4168, Issue #4167) Enhancements @@ -482,6 +488,7 @@ Deprecations * 2.6.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix AtomGroup.center_of_charge(..., unwrap=True) giving inconsistent (but scientifically correct) results on Intel macOS (Issue #4211) @@ -541,6 +548,7 @@ Deprecations * 2.5.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed using 'rb' mode inside GSDReader to avoid FutureWarning from GSDFile (Issue #4152) * Fixed ChainReader performing transformations correctly across boundaries in the individual trajectories (Issue #3657 #4008 PR #3906) @@ -631,6 +639,7 @@ Changes * 2.4.3 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed DCD reading for large (>2Gb) files (Issue #4039). This was broken for versions 2.4.0, 2.4.1 and 2.4.2 * Fix element parsing from PSF files tests read via Parmed (Issue #4015) @@ -641,6 +650,7 @@ Fixes * 2.4.2 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * np.histogramdd calls in :class:`DensityAnalysis` now pass the `density` argument rather than the NumPy 1.24 removed `normed` (PR #3976) * visualization.streamlines_3D and visualization.streamlines no longer @@ -653,6 +663,7 @@ Fixes * 2.4.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixes pyproject.toml minimum pin for biopython @@ -662,6 +673,7 @@ Fixes * 2.4.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Update hbond analysis doc string to use exclusive bound language (Issue #3847) * XTC and TRR readers now fail with IOError when a status except EOK (=0) is @@ -755,6 +767,7 @@ Deprecations * 2.3.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixes awk call in deploy.yaml tests for macos runners (Issue #3693) * add a 0.5 for correct midpoints in hole analysis (Issue #3715) * Fix reading error when PDB CONECT records are corrupt. (Issue #988) @@ -798,6 +811,7 @@ Deprecations * 2.2.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed issue where Reader.copy() did not preserve optional arguments and kwargs passed to the original class (Issue #3664, PR #3685) * Iterating a SingleFrameReaderBase trajectory now rewinds the @@ -895,6 +909,7 @@ Deprecations * 2.1.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Use uint64_t loop counters in C level distance functions to avoid overflow for large arrays (Issue #3512, PR #3513). * Prevents attempts to close an already closed NamedStream (Issue #3386) @@ -946,6 +961,7 @@ Deprecations * 2.0.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Added an `rdkit_kwargs` parameter to pass arguments to the RDKitConverter from the `select_atoms` method (Issue #3319, PR #3324) * Fixes AtomGroup.unique, ResidueGroup.unique, SegmentGroup.unique not @@ -1309,6 +1325,7 @@ Deprecations * 1.1.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Remove absolute paths from package upload to pypi. @@ -1318,6 +1335,7 @@ Fixes * 1.1.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Removes use of absolute paths in setup.py to avoid Windows installation failures (Issue #3129) * Adds test for crashes caused by small box NSGrid searches (Issue #2670) @@ -1361,6 +1379,7 @@ Enhancements * 1.0.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Due to issues with the reliability/accuracy of `nsgrid`, this method is currently not recommended for use. It has also been removed as an option from lib.capped_distance and lib.self_capped_distance. Please use PKDTree @@ -1423,6 +1442,7 @@ Deprecations * 1.0.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * MOL2Writer now accepts both Universes and AtomgGroups (Issue #2717) * Use user-provided `remark` in `XYZWriter` (Issue #2692) * Added more informative error messages about topology attributes @@ -1625,6 +1645,7 @@ Deprecations * 0.20.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * The NetCDF writer now writes `cell_angle` units as `degree` instead of `degrees` in accordance with the AMBER NetCDF convention (Issue #2327). * Fixed installation without Cython (Issue #2337) @@ -1698,6 +1719,7 @@ Changes * changed the water bridge analysis output format (PR #2087) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * fixed ChainReader setting format with format keyword (Issue #2334) * fixed lack of check for scaling of NCDFReader velocities (Issue #2323) * fixed PDBReader and PDBWriter newlines for PDB header (Issue #2324) @@ -1769,6 +1791,7 @@ Changes * 0.19.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * limit output of Chainreader __repr__ (#2109) * added missing docs for lib.pkdtree (#2104) * Added sphinx markup for FrameIterator (#2106) @@ -1857,6 +1880,7 @@ Enhancements (where possible) (Issue #1041 PR #2080) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Rewind in the SingleFrameReader now reads the frame from the file (Issue #1929) * Fixed order of indices in Angle/Dihedral/Improper repr * coordinates.memory.MemoryReader now takes np.ndarray only (Issue #1685) @@ -1956,6 +1980,7 @@ Enhancements * Added AtomGroup, ResidueGroup and SegmentGroup to the top module namespace Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed MPI fork() warning when importing MDAnalysis in an Infiniband-enabled MPI environment (PR #1794) * Fixed waterdynamics SurvivalProbability ignoring the t0 start time @@ -2034,6 +2059,7 @@ Deprecations * timeseries keyword format is replaced by order. The format keyword is marked for deprecation in 1.0 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed analysis.psa.dist_mat_to_vec not returning int values (Issue #1507) * Fixed triclinic PBC transform for a- and b- axes (Issue #1697) * Fixed nuclinfo.tors() not converting delta (Issue #1572) @@ -2098,6 +2124,7 @@ Deprecations * deprecated the core.flag registry for 1.0 (Issue #782) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * fixed GROWriter truncating long resids from the wrong end (Issue #1395) * Fixed dtype of numpy arrays to accomodate 32 bit architectures (Issue #1362) * Groups are hashable on python 3 (Issue #1397) @@ -2129,6 +2156,7 @@ Enhancements * made online docs responsive with the Alabaster Sphinx theme (#378) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * In Universe.transfer_to_memory(): dt is now adjusted with step (Issue #1310) * Various documentation sphinx errors (PR #1312) * Bugfix in confdistmatrix.get_distance_matrix; now works on all trajectory types. @@ -2222,6 +2250,7 @@ Enhancements weights="size") by about 5x (partially Issue #1191) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Trajectory slicing made completely Pythonic (Issue #918 PR #1195) * Argument validation of dist_mat_to_vec is fixed (#597 PR #1183) * Give correct error when the topology file format is not recognized (Issue #982) @@ -2357,6 +2386,7 @@ Enhancements * Iteration and seeking in PDB files made faster (Issue #848) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed TypeError in PSAnalysis heatmap-dendrogram plotting (Issue #1018) * ENT file format added to PDB Readers/Writers/Parsers (Issue #834) * rmsd now returns proper value when given array of weights (Issue #814) @@ -2451,6 +2481,7 @@ Changes before trying to deduce the format from file extension. (Issue #712) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Syntax error corrected in psa.py (Issue #738) * XDR file seeking and telling working again for large files (Issue #677). @@ -2530,6 +2561,7 @@ Changes * MOL2Writer.write now only writes a single frame (Issue #521) Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed select_atoms requiring a trajectory be loaded (Issue #270) * AtomGroup timesteps no longer cached (Issue #606) @@ -2577,6 +2609,7 @@ Enhancements Changes Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed OpenMP detection on Linux/OSX #459 * Fixed reading of LAMMPS trajectory times: default unit ought to be fs and not ps @@ -2631,6 +2664,7 @@ Changes * To build with debug-symbols use 'MDA_DEBUG_CFLAGS' instead of 'MDA_DEBUG_CFLAGS' Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed minor issue in lib.mdamath.make_whole where if all bonds were correctly sized, it wouldn't notice that multiple fragments had been given. (Issue #445) @@ -3227,6 +3261,7 @@ Testsuite * 0.7.5.1 Fixes + * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * added: missing files (Issue 95) * removed: unused delaunay-related files diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 7a0510eab03..0d832d1198c 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -317,7 +317,7 @@ def fit(self): def _conclude(self): norm = np.linspace(self.chainlength - 1, 1, self.chainlength - 1) - norm *= len(self._atomgroups) * self._trajectory.n_frames + norm *= len(self._atomgroups) * self.n_frames self.results.bond_autocorrelation = ( self.results.raw_bond_autocorr / norm ) diff --git a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py index 54005c8a667..a0c228c5ab2 100644 --- a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py +++ b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py @@ -177,3 +177,51 @@ def test_circular(self): u.add_TopologyAttr(Bonds(bondlist)) with pytest.raises(ValueError, match="Cyclical"): polymer.sort_backbone(u.atoms) + + +class TestPersistenceLengthSlicedNormalization(object): + # Regression test for Issue #5453. + @staticmethod + @pytest.fixture(scope="class") + def u(): + return mda.Universe(TRZ_psf, TRZ) + + @staticmethod + @pytest.fixture(scope="class") + def chains(u): + backbones = [ + chain.select_atoms("not name O* H*") for chain in u.atoms.fragments + ] + return [polymer.sort_backbone(bb) for bb in backbones] + + def test_full_trajectory_unchanged(self, chains): + p = polymer.PersistenceLength(chains).run() + + chainlength = len(chains[0]) + norm = np.linspace(chainlength - 1, 1, chainlength - 1) + norm *= len(chains) * p.n_frames + + assert p.n_frames == p._trajectory.n_frames + expected = p.results.raw_bond_autocorr / norm + assert_almost_equal(p.results.bond_autocorrelation, expected, decimal=6) + + def test_sliced_run_normalized_by_frames_analyzed(self, chains, u): + n_sliced_frames = 3 + assert n_sliced_frames < u.trajectory.n_frames + + p = polymer.PersistenceLength(chains).run(stop=n_sliced_frames) + + assert p.n_frames == n_sliced_frames + assert p.n_frames != p._trajectory.n_frames + + chainlength = len(chains[0]) + norm = np.linspace(chainlength - 1, 1, chainlength - 1) + norm *= len(chains) * n_sliced_frames + expected = p.results.raw_bond_autocorr / norm + + assert_almost_equal(p.results.bond_autocorrelation, expected, decimal=6) + + wrong_norm = np.linspace(chainlength - 1, 1, chainlength - 1) + wrong_norm *= len(chains) * u.trajectory.n_frames + wrong = p.results.raw_bond_autocorr / wrong_norm + assert not np.allclose(p.results.bond_autocorrelation, wrong) From c25337c2a07befa4b241db2139e931b797294b1a Mon Sep 17 00:00:00 2001 From: srikar jy Date: Mon, 17 Aug 2026 19:58:17 -0400 Subject: [PATCH 2/2] Address review feedback on PersistenceLength normalization fix * Remove the CHANGELOG entry that was accidentally duplicated across about 25 historical release sections. It now appears only once, under the current unreleased 2.11.0 header. * Replace the tautological regression test, where the expected value was derived from the class under test, with an independent reference implementation that recomputes the normalized bond autocorrelation directly from atom positions, without touching PersistenceLength or its results. * Extend test coverage to step=, start/stop/step, frames=[...], and the multiprocessing backend, beyond the original stop= only case. * Add a guard in _conclude() that raises a clear ValueError when a run selects zero frames, instead of silently returning a fabricated all zero curve before this fix, or crashing with an opaque scipy NaN error after this fix, found while testing the above. * Document, without fixing since it is out of scope for this issue, two existing unrelated edge cases found during testing: a one atom chain crashing identically before and after this fix, and a numpy boolean frames= mask being misread by AnalysisBase computation grouping, pinned as xfail for visibility. --- package/CHANGELOG | 34 --- package/MDAnalysis/analysis/polymer.py | 5 + .../analysis/test_persistencelength.py | 203 ++++++++++++++++-- 3 files changed, 190 insertions(+), 52 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index d0fb038f75b..ab1803bb484 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -121,7 +121,6 @@ Deprecations * 2.10.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix incorrect conversion factors for speed units A/fs, A/us, A/ms in MDAnalysis.units (Issue #5051, PR #5053) * `analysis.polymer.sort_backbone` is now working for discontinuous polymers @@ -214,7 +213,6 @@ Deprecations * 2.9.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Add support for TPR files produced by GROMACS 2024.4 and GROMACS 2025.0 * Fixes invalid default unit from Angstrom to Angstrom^{-3} for convert_density() function. (Issue #4829) @@ -257,7 +255,6 @@ Changes * 2.8.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix `test_imd.py` test failures after imdclient 0.2.4 increased default timeout by ensuring EOF sent (Issue #5442, PR #5443) * Allows bond/angle/dihedral connectivity to be guessed additively with @@ -393,7 +390,6 @@ Deprecations * 2.7.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix doctest errors of lib/transformations.py (Issue #3925, PR #4370) * Fix documentation building errors due to `html_static_path` and unindented block of code (Issue #4362, PR #4365) @@ -463,7 +459,6 @@ Deprecations * 2.6.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix Atom type guessing error (PR #4168, Issue #4167) Enhancements @@ -488,7 +483,6 @@ Deprecations * 2.6.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fix AtomGroup.center_of_charge(..., unwrap=True) giving inconsistent (but scientifically correct) results on Intel macOS (Issue #4211) @@ -548,7 +542,6 @@ Deprecations * 2.5.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed using 'rb' mode inside GSDReader to avoid FutureWarning from GSDFile (Issue #4152) * Fixed ChainReader performing transformations correctly across boundaries in the individual trajectories (Issue #3657 #4008 PR #3906) @@ -639,7 +632,6 @@ Changes * 2.4.3 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed DCD reading for large (>2Gb) files (Issue #4039). This was broken for versions 2.4.0, 2.4.1 and 2.4.2 * Fix element parsing from PSF files tests read via Parmed (Issue #4015) @@ -650,7 +642,6 @@ Fixes * 2.4.2 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * np.histogramdd calls in :class:`DensityAnalysis` now pass the `density` argument rather than the NumPy 1.24 removed `normed` (PR #3976) * visualization.streamlines_3D and visualization.streamlines no longer @@ -663,7 +654,6 @@ Fixes * 2.4.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixes pyproject.toml minimum pin for biopython @@ -673,7 +663,6 @@ Fixes * 2.4.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Update hbond analysis doc string to use exclusive bound language (Issue #3847) * XTC and TRR readers now fail with IOError when a status except EOK (=0) is @@ -767,7 +756,6 @@ Deprecations * 2.3.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixes awk call in deploy.yaml tests for macos runners (Issue #3693) * add a 0.5 for correct midpoints in hole analysis (Issue #3715) * Fix reading error when PDB CONECT records are corrupt. (Issue #988) @@ -811,7 +799,6 @@ Deprecations * 2.2.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed issue where Reader.copy() did not preserve optional arguments and kwargs passed to the original class (Issue #3664, PR #3685) * Iterating a SingleFrameReaderBase trajectory now rewinds the @@ -909,7 +896,6 @@ Deprecations * 2.1.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Use uint64_t loop counters in C level distance functions to avoid overflow for large arrays (Issue #3512, PR #3513). * Prevents attempts to close an already closed NamedStream (Issue #3386) @@ -961,7 +947,6 @@ Deprecations * 2.0.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Added an `rdkit_kwargs` parameter to pass arguments to the RDKitConverter from the `select_atoms` method (Issue #3319, PR #3324) * Fixes AtomGroup.unique, ResidueGroup.unique, SegmentGroup.unique not @@ -1325,7 +1310,6 @@ Deprecations * 1.1.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Remove absolute paths from package upload to pypi. @@ -1335,7 +1319,6 @@ Fixes * 1.1.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Removes use of absolute paths in setup.py to avoid Windows installation failures (Issue #3129) * Adds test for crashes caused by small box NSGrid searches (Issue #2670) @@ -1379,7 +1362,6 @@ Enhancements * 1.0.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Due to issues with the reliability/accuracy of `nsgrid`, this method is currently not recommended for use. It has also been removed as an option from lib.capped_distance and lib.self_capped_distance. Please use PKDTree @@ -1442,7 +1424,6 @@ Deprecations * 1.0.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * MOL2Writer now accepts both Universes and AtomgGroups (Issue #2717) * Use user-provided `remark` in `XYZWriter` (Issue #2692) * Added more informative error messages about topology attributes @@ -1645,7 +1626,6 @@ Deprecations * 0.20.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * The NetCDF writer now writes `cell_angle` units as `degree` instead of `degrees` in accordance with the AMBER NetCDF convention (Issue #2327). * Fixed installation without Cython (Issue #2337) @@ -1719,7 +1699,6 @@ Changes * changed the water bridge analysis output format (PR #2087) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * fixed ChainReader setting format with format keyword (Issue #2334) * fixed lack of check for scaling of NCDFReader velocities (Issue #2323) * fixed PDBReader and PDBWriter newlines for PDB header (Issue #2324) @@ -1791,7 +1770,6 @@ Changes * 0.19.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * limit output of Chainreader __repr__ (#2109) * added missing docs for lib.pkdtree (#2104) * Added sphinx markup for FrameIterator (#2106) @@ -1880,7 +1858,6 @@ Enhancements (where possible) (Issue #1041 PR #2080) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Rewind in the SingleFrameReader now reads the frame from the file (Issue #1929) * Fixed order of indices in Angle/Dihedral/Improper repr * coordinates.memory.MemoryReader now takes np.ndarray only (Issue #1685) @@ -1980,7 +1957,6 @@ Enhancements * Added AtomGroup, ResidueGroup and SegmentGroup to the top module namespace Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed MPI fork() warning when importing MDAnalysis in an Infiniband-enabled MPI environment (PR #1794) * Fixed waterdynamics SurvivalProbability ignoring the t0 start time @@ -2059,7 +2035,6 @@ Deprecations * timeseries keyword format is replaced by order. The format keyword is marked for deprecation in 1.0 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed analysis.psa.dist_mat_to_vec not returning int values (Issue #1507) * Fixed triclinic PBC transform for a- and b- axes (Issue #1697) * Fixed nuclinfo.tors() not converting delta (Issue #1572) @@ -2124,7 +2099,6 @@ Deprecations * deprecated the core.flag registry for 1.0 (Issue #782) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * fixed GROWriter truncating long resids from the wrong end (Issue #1395) * Fixed dtype of numpy arrays to accomodate 32 bit architectures (Issue #1362) * Groups are hashable on python 3 (Issue #1397) @@ -2156,7 +2130,6 @@ Enhancements * made online docs responsive with the Alabaster Sphinx theme (#378) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * In Universe.transfer_to_memory(): dt is now adjusted with step (Issue #1310) * Various documentation sphinx errors (PR #1312) * Bugfix in confdistmatrix.get_distance_matrix; now works on all trajectory types. @@ -2250,7 +2223,6 @@ Enhancements weights="size") by about 5x (partially Issue #1191) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Trajectory slicing made completely Pythonic (Issue #918 PR #1195) * Argument validation of dist_mat_to_vec is fixed (#597 PR #1183) * Give correct error when the topology file format is not recognized (Issue #982) @@ -2386,7 +2358,6 @@ Enhancements * Iteration and seeking in PDB files made faster (Issue #848) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed TypeError in PSAnalysis heatmap-dendrogram plotting (Issue #1018) * ENT file format added to PDB Readers/Writers/Parsers (Issue #834) * rmsd now returns proper value when given array of weights (Issue #814) @@ -2481,7 +2452,6 @@ Changes before trying to deduce the format from file extension. (Issue #712) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Syntax error corrected in psa.py (Issue #738) * XDR file seeking and telling working again for large files (Issue #677). @@ -2561,7 +2531,6 @@ Changes * MOL2Writer.write now only writes a single frame (Issue #521) Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed select_atoms requiring a trajectory be loaded (Issue #270) * AtomGroup timesteps no longer cached (Issue #606) @@ -2609,7 +2578,6 @@ Enhancements Changes Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed OpenMP detection on Linux/OSX #459 * Fixed reading of LAMMPS trajectory times: default unit ought to be fs and not ps @@ -2664,7 +2632,6 @@ Changes * To build with debug-symbols use 'MDA_DEBUG_CFLAGS' instead of 'MDA_DEBUG_CFLAGS' Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * Fixed minor issue in lib.mdamath.make_whole where if all bonds were correctly sized, it wouldn't notice that multiple fragments had been given. (Issue #445) @@ -3261,7 +3228,6 @@ Testsuite * 0.7.5.1 Fixes - * `PersistenceLength._conclude()` now normalizes `results.bond_autocorrelation` using the number of frames actually analyzed instead of the total trajectory frame count, fixing incorrect results when `run()` is called with `start`/`stop`/`step` (Issue #5453). * added: missing files (Issue 95) * removed: unused delaunay-related files diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 0d832d1198c..70fea80fa01 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/package/MDAnalysis/analysis/polymer.py @@ -316,6 +316,11 @@ def fit(self): return self.results.fit def _conclude(self): + if self.n_frames == 0: + raise ValueError( + "PersistenceLength.run() analyzed zero frames; check the " + "start/stop/step or frames selection passed to run()." + ) norm = np.linspace(self.chainlength - 1, 1, self.chainlength - 1) norm *= len(self._atomgroups) * self.n_frames self.results.bond_autocorrelation = ( diff --git a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py index a0c228c5ab2..360d4c66513 100644 --- a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py +++ b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py @@ -31,7 +31,7 @@ import matplotlib import matplotlib.pyplot as plt -from numpy.testing import assert_almost_equal, assert_equal +from numpy.testing import assert_allclose, assert_almost_equal, assert_equal from MDAnalysisTests.datafiles import Plength, TRZ_psf, TRZ @@ -179,6 +179,39 @@ def test_circular(self): polymer.sort_backbone(u.atoms) +def _reference_bond_autocorrelation(chains, frame_indices): + """Independently compute the expected normalized bond autocorrelation + for exactly the given sequence of trajectory frame indices. + + Mirrors the mathematical definition directly from atom positions, + without calling into :class:`~MDAnalysis.analysis.polymer. + PersistenceLength` or reusing any of its results, so this stays a + meaningful check even if the class under test is itself broken. + Frame indices are visited explicitly (in the order given, duplicates + allowed) rather than via a slice, so this can independently mirror + any of `run`'s `start`/`stop`/`step`/`frames` selections. See Issue + #5453. + """ + chainlength = len(chains[0]) + raw = np.zeros(chainlength - 1, dtype=np.float64) + + universe = chains[0].universe + n_frames = 0 + for frame in frame_indices: + universe.trajectory[frame] + n_frames += 1 + for chain in chains: + vecs = chain.positions[1:] - chain.positions[:-1] + vecs = vecs / np.sqrt((vecs * vecs).sum(axis=1))[:, None] + inner_pr = np.inner(vecs, vecs) + for i in range(chainlength - 1): + raw[: (chainlength - 1) - i] += inner_pr[i, i:] + + norm = np.linspace(chainlength - 1, 1, chainlength - 1) + norm *= len(chains) * n_frames + return raw / norm, n_frames + + class TestPersistenceLengthSlicedNormalization(object): # Regression test for Issue #5453. @staticmethod @@ -194,16 +227,23 @@ def chains(u): ] return [polymer.sort_backbone(bb) for bb in backbones] - def test_full_trajectory_unchanged(self, chains): + def test_full_trajectory_unchanged(self, chains, u): p = polymer.PersistenceLength(chains).run() - chainlength = len(chains[0]) - norm = np.linspace(chainlength - 1, 1, chainlength - 1) - norm *= len(chains) * p.n_frames + assert p.n_frames == u.trajectory.n_frames - assert p.n_frames == p._trajectory.n_frames - expected = p.results.raw_bond_autocorr / norm - assert_almost_equal(p.results.bond_autocorrelation, expected, decimal=6) + expected, n_frames = _reference_bond_autocorrelation( + chains, range(u.trajectory.n_frames) + ) + assert n_frames == u.trajectory.n_frames + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) def test_sliced_run_normalized_by_frames_analyzed(self, chains, u): n_sliced_frames = 3 @@ -212,16 +252,143 @@ def test_sliced_run_normalized_by_frames_analyzed(self, chains, u): p = polymer.PersistenceLength(chains).run(stop=n_sliced_frames) assert p.n_frames == n_sliced_frames - assert p.n_frames != p._trajectory.n_frames - - chainlength = len(chains[0]) - norm = np.linspace(chainlength - 1, 1, chainlength - 1) - norm *= len(chains) * n_sliced_frames - expected = p.results.raw_bond_autocorr / norm + assert p.n_frames != u.trajectory.n_frames - assert_almost_equal(p.results.bond_autocorrelation, expected, decimal=6) + expected, n_frames = _reference_bond_autocorrelation( + chains, range(n_sliced_frames) + ) + assert n_frames == n_sliced_frames + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) - wrong_norm = np.linspace(chainlength - 1, 1, chainlength - 1) - wrong_norm *= len(chains) * u.trajectory.n_frames - wrong = p.results.raw_bond_autocorr / wrong_norm + # Sanity check: reproduce the old (buggy) behavior by normalizing + # the *sliced* raw sum with the *total* trajectory frame count + # instead of n_sliced_frames, and confirm that no longer matches + # the class's output -- i.e. this test is sensitive to the bug + # described in Issue #5453. + wrong = expected * n_sliced_frames / u.trajectory.n_frames assert not np.allclose(p.results.bond_autocorrelation, wrong) + + @pytest.mark.parametrize( + "run_kwargs, frame_indices", + [ + pytest.param({"step": 2}, [0, 2, 4], id="step_only"), + pytest.param( + {"start": 1, "stop": 6, "step": 2}, + [1, 3, 5], + id="start_stop_step", + ), + pytest.param({"frames": [0, 2, 5]}, [0, 2, 5], id="frames_list"), + pytest.param( + {"frames": [4, 1, 3]}, + [4, 1, 3], + id="frames_list_unordered", + ), + ], + ) + def test_various_frame_selections_normalized_correctly( + self, chains, u, run_kwargs, frame_indices + ): + # Regression test for Issue #5453: the normalization fix relies + # on self.n_frames, which AnalysisBase derives the same way + # regardless of whether the run was sliced via start/stop/step + # or via an explicit frames= list. Exercise more than just a + # plain stop= slice to guard against that assumption breaking. + p = polymer.PersistenceLength(chains).run(**run_kwargs) + + assert p.n_frames == len(frame_indices) + assert p.n_frames != u.trajectory.n_frames + + expected, n_frames = _reference_bond_autocorrelation( + chains, frame_indices + ) + assert n_frames == len(frame_indices) + # atol/rtol (looser than a plain decimal=6 comparison) account for + # PersistenceLength accumulating raw_bond_autocorr in float32 while + # this reference accumulates in float64; their rounding drift is + # a couple of float32 ULPs at this trajectory's frame/chain count, + # not a sign of disagreement. + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + ) + + def test_sliced_run_parallel_matches_serial(self, chains, u): + # Regression test for Issue #5453: the bug was introduced by the + # multi-worker refactor in PR #5074, but neither existing test + # actually exercised backend="multiprocessing" together with a + # sliced run. _conclude() only ever runs once, on the main + # process, after worker results are merged -- so self.n_frames + # at normalization time should be the total sliced frame count + # regardless of how work was chunked across workers. Confirm + # that holds, rather than only checking it by hand. + n_sliced_frames = 5 + assert n_sliced_frames < u.trajectory.n_frames + + serial = polymer.PersistenceLength(chains).run(stop=n_sliced_frames) + parallel = polymer.PersistenceLength(chains).run( + stop=n_sliced_frames, backend="multiprocessing", n_workers=2 + ) + + assert serial.n_frames == parallel.n_frames == n_sliced_frames + assert_allclose( + serial.results.bond_autocorrelation, + parallel.results.bond_autocorrelation, + rtol=1e-4, + atol=5e-6, + ) + + def test_empty_frame_selection_raises_clear_error(self, chains): + # Edge case found while verifying the fix for Issue #5453: an + # empty frame selection drives self.n_frames to 0, which used + # to silently divide by the (nonzero) total trajectory frame + # count and return a fabricated all-zero curve with no warning. + # After the fix, the zero-frame divisor would instead produce + # NaNs that surface as an opaque scipy error deep inside + # curve_fit. Raise a clear, immediate error instead. + with pytest.raises(ValueError, match="analyzed zero frames"): + polymer.PersistenceLength(chains).run(start=2, stop=2) + + def test_single_atom_chain_unaffected_by_normalization_fix(self, u): + # Not a regression from Issue #5453's fix: a chain with only one + # atom has no bonds to correlate (chainlength - 1 == 0), so + # raw_bond_autocorr and norm are both length-0 arrays regardless + # of the run() call's start/stop/step. This has always raised + # from inside _perform_fit()'s curve_fit call, independent of + # which frame count _conclude() normalizes by. Documented here + # so it doesn't get mistaken for something this PR should fix. + one_atom_chain = [u.atoms.fragments[0][:1]] + with pytest.raises(ValueError): + polymer.PersistenceLength(one_atom_chain).run() + + @pytest.mark.xfail( + reason=( + "Pre-existing AnalysisBase bug, independent of Issue #5453: " + "_setup_computation_groups checks `isinstance(obj, bool)`, " + "which is False for numpy bool scalars, so a numpy boolean " + "frames= mask is treated as an integer index array (True/" + "False read as 1/0) instead of a mask. A plain Python list " + "of bools is unaffected. Left xfail to document the " + "framework-level issue without attempting to fix " + "AnalysisBase in this PR." + ), + strict=True, + ) + def test_numpy_bool_mask_frames_selection(self, chains, u): + mask = np.zeros(u.trajectory.n_frames, dtype=bool) + mask[[0, 2, 4]] = True + + p = polymer.PersistenceLength(chains).run(frames=mask) + + assert p.n_frames == int(mask.sum()) + expected, n_frames = _reference_bond_autocorrelation( + chains, np.flatnonzero(mask) + ) + assert_allclose( + p.results.bond_autocorrelation, expected, rtol=1e-4, atol=5e-6 + )