Skip to content

PersistenceLength.run(start=, stop=) rescales entire results.bond_autocorrelation array by n_frames_processed / n_frames_total_trajectory #5453

Description

@lslumass

When run() is given a partial start/stop range, the entire results.bond_autocorrelation curve — not just the x=0 point — comes out uniformly rescaled by n_frames_processed / n_frames_total_trajectory, rather than correctly normalized by the number of frames actually processed.

Reproduction (30-atom backbone, 50,000-frame trajectory), plotting the measured C(x) and fitted curve for two ranges on the same chain:

import MDAnalysis as mda
from MDAnalysis.analysis import polymer
import matplotlib.pyplot as plt

u = mda.Universe("conf.psf", "system.xtc")
bb = u.select_atoms("segid RA0 and name P")

for start, stop in [(1, 25000), (1, 50000)]:
    plen = polymer.PersistenceLength([bb])
    plen.run(start=start, stop=stop, step=1)
    plt.figure()
    plt.plot(plen.results.x, plen.results.bond_autocorrelation, "ro", label="Result")
    plt.plot(plen.results.x, plen.results.fit, "b-", label="Fit")
    plt.title(f"start={start}, stop={stop}")
    plt.xlabel("x"); plt.ylabel("C(x)"); plt.legend()
    plt.savefig(f"lp_{start}_{stop}.png")

Comparing the two resulting curves point-by-point, every value in the stop=25000 curve is almost exactly half the corresponding value in the stop=50000 curve (e.g. C(0)=0.50 vs 1.00, C(x≈6)=0.27 vs 0.55, C(x≈12)=0.18 vs 0.37) — matching 25000/50000 = 0.5 exactly. This shows the mis-normalization is a constant multiplicative rescaling of the whole array, not a single-point artifact, and it corrupts the resulting exponential fit and reported lp for any partial-range analysis. Only start=0/stop=full trajectory length happens to self-correct (numerator == denominator). No warning or error is raised. MDAnalysis version: 2.10.0

Do you have any misunderstanding?

Image Image

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions