From e0b0ecd1c3e72fb25a86fc332a8091f12efc335b Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Mon, 4 May 2026 14:08:18 +0200 Subject: [PATCH 1/4] fix: allow direct scheme with XTC by making box length optional; add CSS alias MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - `--box` for `direct` is now optional (no default_value in clap); when omitted without an XTC the 200 Å fallback is applied in lib.rs, preserving existing single-frame behaviour - Removes the hard error that fired even when -b was not supplied, allowing trajectory-averaged SAXS via `direct` without specifying a box - Add CSS (disulfide-bonded cysteine) as an alias for CYS in single_bead.yaml Co-Authored-By: Claude Sonnet 4.6 --- assets/single_bead.yaml | 1 + src/cli.rs | 6 +++--- src/lib.rs | 2 ++ 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/assets/single_bead.yaml b/assets/single_bead.yaml index d417599..1ae15b3 100644 --- a/assets/single_bead.yaml +++ b/assets/single_bead.yaml @@ -824,6 +824,7 @@ HLYS: { ff: !Alias LYS } HTYR: { ff: !Alias TYR } HCYS: { ff: !Alias CYS } CYX: { ff: !Alias CYS } +CSS: { radius: 3.6, formfactor: !Alias CYS } # PDB atom name aliases — map multi-character PDB names to their # element form factors so atomic-resolution structures work with diff --git a/src/cli.rs b/src/cli.rs index e577088..56b4f47 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -47,8 +47,8 @@ pub enum Commands { #[clap(short = 'p', long, default_value = "6")] pmax: usize, /// Cubic box side length (Å) - #[clap(short = 'b', long = "box", default_value = "200.0")] - side_length: f64, + #[clap(short = 'b', long = "box")] + side_length: Option, #[clap(flatten)] solvent: SolventArgs, }, @@ -291,7 +291,7 @@ fn command_to_scheme( solvent, } => ( IntensityScheme::Direct { pmax }, - Some([side_length; 3]), + side_length.map(|s| [s; 3]), solvent.to_config(), ), Commands::Multipole { diff --git a/src/lib.rs b/src/lib.rs index 01b2aad..806c0bc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -212,6 +212,8 @@ impl Pripps { )); } structure.box_sides = Some(Vector3::from(box_sides)); + } else if xtcfile.is_none() { + structure.box_sides = Some(Vector3::from([200.0; 3])); } let active_solvent = solvent.filter(|cfg| cfg.is_active()); From 358a13273d2b63658963e96b03df53f3b164cce4 Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Tue, 5 May 2026 09:22:47 +0200 Subject: [PATCH 2/4] feat: add beta(q) decoupling parameter to multipole CSV output MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Exposes β(q) = ⟨F⟩²_Ω / P(q) = 4π|A₀⁰|² / P(q) for all multipole paths (vacuum and solvent-corrected). This is the decoupling parameter needed to convert a true S(q) into the experimentally measured effective structure factor: S_eff(q) = 1 + β(q)·(S(q) − 1), required for concentrated-solution and mixture modelling. - Add `beta: Option>` to `Intensity` struct - Compute beta inline in `vacuum_intensity` (atoms-only path) - Compute beta inline in `MultipoleModel::intensity` (solvent path), using the combined A₀⁰ − v·C₀⁰ + ρc·B₀⁰ monopole coefficient - Update `write_csv` to emit a `beta` column when present; Debye and Direct paths are unaffected (beta stays None) Co-Authored-By: Claude Sonnet 4.6 --- src/cli.rs | 28 ++++++++++++++++++++++++---- src/explicit.rs | 1 + src/lib.rs | 6 ++++++ src/multipole/mod.rs | 21 ++++++++++++++++----- 4 files changed, 47 insertions(+), 9 deletions(-) diff --git a/src/cli.rs b/src/cli.rs index 56b4f47..cca6974 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -327,14 +327,20 @@ fn write_csv(path: &Path, curve: &Intensity) -> Result<()> { let mut out = BufWriter::new(file); log::info!("Writing {}", path.display()); - match &curve.contributions { - None => { + match (&curve.contributions, &curve.beta) { + (None, None) => { writeln!(out, "q,total")?; for (q, total) in curve.q.iter().zip(&curve.total) { writeln!(out, "{q:.6},{total:.6}")?; } } - Some(c) => { + (None, Some(beta)) => { + writeln!(out, "q,total,beta")?; + for ((q, total), beta) in curve.q.iter().zip(&curve.total).zip(beta) { + writeln!(out, "{q:.6},{total:.6},{beta:.6}")?; + } + } + (Some(c), None) => { writeln!(out, "q,total,atoms,excluded,hydration")?; for ((((q, total), atoms), excluded), hydration) in curve .q @@ -343,10 +349,24 @@ fn write_csv(path: &Path, curve: &Intensity) -> Result<()> { .zip(&c.atoms) .zip(&c.excluded) .zip(&c.hydration) + { + writeln!(out, "{q:.6},{total:.6},{atoms:.6},{excluded:.6},{hydration:.6}")?; + } + } + (Some(c), Some(beta)) => { + writeln!(out, "q,total,atoms,excluded,hydration,beta")?; + for (((((q, total), atoms), excluded), hydration), beta) in curve + .q + .iter() + .zip(&curve.total) + .zip(&c.atoms) + .zip(&c.excluded) + .zip(&c.hydration) + .zip(beta) { writeln!( out, - "{q:.6},{total:.6},{atoms:.6},{excluded:.6},{hydration:.6}" + "{q:.6},{total:.6},{atoms:.6},{excluded:.6},{hydration:.6},{beta:.6}" )?; } } diff --git a/src/explicit.rs b/src/explicit.rs index 6cfecbb..78ffc48 100644 --- a/src/explicit.rs +++ b/src/explicit.rs @@ -194,6 +194,7 @@ impl DirectTransform { excluded, hydration, }), + beta: None, }) } } diff --git a/src/lib.rs b/src/lib.rs index 806c0bc..7f9a83e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -48,11 +48,17 @@ pub(crate) type Vector3 = nalgebra::Vector3; /// `contributions` is populated when Multipole runs with an active /// [`SolventConfig`], carrying the per-contribution diagnostic /// columns. It is `None` for vacuum Multipole, Debye, and Direct. +/// +/// `beta` is populated for all Multipole paths (vacuum and solvent). +/// It is the decoupling parameter β(q) = ⟨F⟩²/P(q) = 4π|A₀⁰|²/P(q), +/// needed to convert a true S(q) into the experimentally measured +/// effective structure factor: S_eff(q) = 1 + β(q)·(S(q) − 1). #[derive(Debug, Clone, Default, PartialEq)] pub struct Intensity { pub q: Vec, pub total: Vec, pub contributions: Option, + pub beta: Option>, } /// Per-contribution intensity curves emitted by the three-contribution diff --git a/src/multipole/mod.rs b/src/multipole/mod.rs index af66c2c..3820515 100644 --- a/src/multipole/mod.rs +++ b/src/multipole/mod.rs @@ -340,16 +340,22 @@ fn vacuum_intensity( l_max_per_q, )?; - let total: Vec = coeffs_per_q + let (total, beta): (Vec, Vec) = coeffs_per_q .par_iter() .zip(l_max_per_q) - .map(|(c, &l_q)| intensity_from_multipoles(c, l_q)) - .collect(); + .map(|(c, &l_q)| { + let p = intensity_from_multipoles(c, l_q); + let mean_amp_sq = 4.0 * std::f64::consts::PI * c[(0, 0)].norm_sqr(); + let beta = if p > 0.0 { mean_amp_sq / p } else { 1.0 }; + (p, beta) + }) + .unzip(); Ok(Intensity { q: sampler.q_values.clone(), total, contributions: None, + beta: Some(beta), }) } @@ -467,7 +473,7 @@ impl MultipoleModel { // overhead would dominate. The expensive O(N·Q·L²) coefficient // build in `build_model` correctly uses `par_iter`. pub fn intensity(&self, volume_scale: f64, contrast_density: f64) -> Intensity { - let per_q: Vec<(f64, f64, f64, f64)> = izip!( + let per_q: Vec<(f64, f64, f64, f64, f64)> = izip!( &self.atoms_coeffs, &self.excluded_coeffs, &self.hydration_coeffs, @@ -475,16 +481,20 @@ impl MultipoleModel { ) .map(|(a, e, h, &l_q)| { let total = combine_contributions(a, e, h, l_q, volume_scale, contrast_density); + let mono = a[(0, 0)] - volume_scale * e[(0, 0)] + contrast_density * h[(0, 0)]; + let mean_amp_sq = 4.0 * std::f64::consts::PI * mono.norm_sqr(); + let beta = if total > 0.0 { mean_amp_sq / total } else { 1.0 }; ( total, intensity_from_multipoles(a, l_q), intensity_from_multipoles(e, l_q), intensity_from_multipoles(h, l_q), + beta, ) }) .collect(); - let (total, atoms, excluded, hydration): (Vec<_>, Vec<_>, Vec<_>, Vec<_>) = + let (total, atoms, excluded, hydration, beta): (Vec<_>, Vec<_>, Vec<_>, Vec<_>, Vec<_>) = per_q.into_iter().multiunzip(); Intensity { @@ -495,6 +505,7 @@ impl MultipoleModel { excluded, hydration, }), + beta: Some(beta), } } From a65ef21fd0a63840f754160658fff28e3e5d6967 Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Wed, 6 May 2026 14:17:36 +0200 Subject: [PATCH 3/4] feat: add beta(q) to direct Fourier transform output MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Extends the decoupling parameter β(q) = ⟨F⟩²_Ω / P(q) to the direct (Miller-index) scheme, for both vacuum and solvent-corrected paths. ⟨F(q)⟩_Ω is computed analytically as Σ_i f_i·sinc(q|r_i − center|), the isotropic (l=0) component — equivalent to the multipole A₀⁰ route. All contributions (atoms, excluded volume, hydration) are centred on the protein geometric centroid so that β is translation-invariant. β → 1 at q = 0 and decreases with q, matching the multipole values. - Add `mean_sinc_amplitude(q, pos, ff, center)` and `centroid()` helpers - `vacuum_intensities`: compute beta after Miller-direction averaging - `solvent_intensities`: compute combined mean amplitude using protein centroid for atoms, excluded-volume, and hydration dummies - Update `lib.rs` doc to note Direct also populates `beta` - Add unit test verifying β ≈ 1 at low q and β < 0.5 at high q Co-Authored-By: Claude Sonnet 4.6 --- src/explicit.rs | 100 ++++++++++++++++++++++++++++++++++++++++++++---- src/lib.rs | 5 ++- 2 files changed, 96 insertions(+), 9 deletions(-) diff --git a/src/explicit.rs b/src/explicit.rs index 78ffc48..9c969df 100644 --- a/src/explicit.rs +++ b/src/explicit.rs @@ -44,6 +44,33 @@ const MILLER_INDEX: [Vector3; 13] = [ Vector3::new(1.0, 1.0, -1.0), ]; +/// Orientationally averaged amplitude ⟨F(q)⟩_Ω = Σ_i f_i · sin(q|r_i − center|)/(q|r_i − center|). +/// +/// `center` must be the particle centroid so that |r_i − center| are +/// distances from the scattering centre. This is the isotropic component, +/// equivalent to the l=0 multipole coefficient: 4π|A₀⁰|² = mean_sinc_amplitude². +/// Used to compute β(q) = ⟨F⟩²_Ω / P(q) for the decoupling approximation. +fn mean_sinc_amplitude( + q_norm: f64, + positions: &[Vector3], + form_factors: &[f64], + center: &Vector3, +) -> f64 { + positions.iter().zip(form_factors).map(|(r, &f)| { + let qr = q_norm * (r - center).norm(); + let sinc = if qr < 1e-10 { 1.0 } else { qr.sin() / qr }; + f * sinc + }).sum() +} + +/// Geometric centroid of a set of positions. +fn centroid(positions: &[Vector3]) -> Vector3 { + if positions.is_empty() { + return Vector3::zeros(); + } + positions.iter().sum::() / positions.len() as f64 +} + /// Calculate I(q) using discrete q vectors /// /// References: @@ -123,9 +150,22 @@ impl DirectTransform { .collect::>>()?; let (q, total): (Vec, Vec) = average_duplicates(table).into_iter().unzip(); + + let center = centroid(&structure.pos); + let beta = q + .iter() + .zip(&total) + .map(|(&qi, &p)| { + let ff = eval_atom_form_factors(qi, &structure.ids, &self.formfactors)?; + let mean_amp = mean_sinc_amplitude(qi, &structure.pos, &ff, ¢er); + Ok(if p > 0.0 { mean_amp * mean_amp / p } else { 1.0 }) + }) + .collect::>>()?; + Ok(Intensity { q, total, + beta: Some(beta), ..Default::default() }) } @@ -178,12 +218,34 @@ impl DirectTransform { let mut atoms = Vec::with_capacity(grouped.len()); let mut excluded = Vec::with_capacity(grouped.len()); let mut hydration = Vec::with_capacity(grouped.len()); - for (qi, t, a, e, h) in grouped { - q.push(qi); - total.push(t); - atoms.push(a); - excluded.push(e); - hydration.push(h); + let mut beta = Vec::with_capacity(grouped.len()); + + for (qi, t, a, e, h) in &grouped { + q.push(*qi); + total.push(*t); + atoms.push(*a); + excluded.push(*e); + hydration.push(*h); + } + + // All three contributions are centered on the same protein centroid so + // that β(q) is frame-of-reference independent (P(q) is already invariant + // under translation because |A|² gains only a common phase factor). + let center = centroid(&structure.pos); + for (&qi, &p) in q.iter().zip(&total) { + let atom_ff = eval_atom_form_factors(qi, &structure.ids, &self.formfactors)?; + let mean_a = mean_sinc_amplitude(qi, &structure.pos, &atom_ff, ¢er); + + let mut ev_ff = Vec::new(); + excluded_dummies.form_factors_at(qi, &mut ev_ff)?; + let mean_c = mean_sinc_amplitude(qi, excluded_dummies.positions(), &ev_ff, ¢er); + + let mut hyd_ff = Vec::new(); + hydration_dummies.form_factors_at(qi, &mut hyd_ff)?; + let mean_b = mean_sinc_amplitude(qi, hydration_dummies.positions(), &hyd_ff, ¢er); + + let mean_amp = mean_a - volume_scale * mean_c + contrast_density * mean_b; + beta.push(if p > 0.0 { mean_amp * mean_amp / p } else { 1.0 }); } Ok(Intensity { @@ -194,7 +256,7 @@ impl DirectTransform { excluded, hydration, }), - beta: None, + beta: Some(beta), }) } } @@ -258,4 +320,28 @@ mod tests { let result = average_duplicates(v); assert_eq!(result, vec![(1.0, 1.5), (2.0, 3.5), (3.0, 0.1)]); } + + #[test] + fn direct_beta_near_one_at_low_q_and_decays() { + // Two atoms offset from origin: centroid correction must be applied. + let positions = vec![ + Vector3::new(10.0, 10.0, 10.0), + Vector3::new(10.0, 10.0, 14.0), + ]; + let form_factors = vec![1.0, 1.0]; + let center = centroid(&positions); + + let beta_low = { + let amp = mean_sinc_amplitude(0.001, &positions, &form_factors, ¢er); + amp * amp / (2.0 * 2.0) // P(q→0) = (Σ f_i)² + }; + let beta_high = { + let amp = mean_sinc_amplitude(5.0, &positions, &form_factors, ¢er); + let p = 4.0; // |f1|² + |f2|² (incoherent limit) + amp * amp / p + }; + + assert!((beta_low - 1.0).abs() < 1e-4, "β at q→0 should be ≈1, got {beta_low}"); + assert!(beta_high < 0.5, "β should decrease at high q, got {beta_high}"); + } } diff --git a/src/lib.rs b/src/lib.rs index 7f9a83e..cfc76d9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -49,10 +49,11 @@ pub(crate) type Vector3 = nalgebra::Vector3; /// [`SolventConfig`], carrying the per-contribution diagnostic /// columns. It is `None` for vacuum Multipole, Debye, and Direct. /// -/// `beta` is populated for all Multipole paths (vacuum and solvent). -/// It is the decoupling parameter β(q) = ⟨F⟩²/P(q) = 4π|A₀⁰|²/P(q), +/// `beta` is populated for all Multipole and Direct paths (vacuum and solvent). +/// It is the decoupling parameter β(q) = ⟨F⟩²/P(q), /// needed to convert a true S(q) into the experimentally measured /// effective structure factor: S_eff(q) = 1 + β(q)·(S(q) − 1). +/// For Multipole, ⟨F⟩ = √(4π)|A₀⁰|. For Direct, ⟨F⟩ = Σ_i f_i sinc(q|r_i|). #[derive(Debug, Clone, Default, PartialEq)] pub struct Intensity { pub q: Vec, From 53fa3a57afd7ac65c7a06bb68e37aacafc51a923 Mon Sep 17 00:00:00 2001 From: Sebastian Pineda Date: Wed, 6 May 2026 17:08:56 +0200 Subject: [PATCH 4/4] fix: correct beta(q) trajectory averaging and PBC wrapping in direct method MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two bugs in the beta(q) calculation for XTC trajectories: 1. Averaging formula: was accumulating per-frame β_t = ⟨F⟩²/P_t and averaging those, which is biased when P_t fluctuates (rotating protein sampled through few Miller directions). Correct formula is β = (⟨⟨F⟩_Ω⟩_conf)² / ⟨P⟩_conf. Fix: accumulate mean sinc amplitude A_t = sqrt(β_t·P_t) per frame, then compute β = ⟨A⟩²/⟨P⟩ after the loop. 2. PBC wrapping: mean_sinc_amplitude uses physical distances |r_i − centroid|, which are corrupted when the protein crosses a periodic box boundary. complex_amplitude is immune (exp(iq·r) is periodic for integer Miller indices), but the sinc-based ⟨F⟩_Ω is not. Fix: unwrap the molecule using minimum-image convention before computing the centroid and distances. Per-frame β now stays ≈ 0.9999 at low q instead of dropping to 0.03–0.24 in frames where the protein straddled the boundary. Also: β is suppressed for multi-protein boxes (frame atom count ≠ topology atom count), where it is physically meaningless. Co-Authored-By: Claude Sonnet 4.6 --- src/explicit.rs | 60 +++++++++++++++++++++++++++++++++-------------- src/trajectory.rs | 20 ++++++++++++++++ 2 files changed, 63 insertions(+), 17 deletions(-) diff --git a/src/explicit.rs b/src/explicit.rs index 9c969df..7a00f82 100644 --- a/src/explicit.rs +++ b/src/explicit.rs @@ -44,19 +44,50 @@ const MILLER_INDEX: [Vector3; 13] = [ Vector3::new(1.0, 1.0, -1.0), ]; +/// Apply minimum-image convention to a displacement vector given orthorhombic box sides. +fn minimum_image(delta: &Vector3, box_sides: &Vector3) -> Vector3 { + delta.zip_map(box_sides, |d, l| d - l * (d / l).round()) +} + +/// Unwrap PBC-wrapped positions so the whole molecule lies in one image. +/// +/// Anchors on the first atom and folds each subsequent atom to the nearest +/// image of the previous one. Returns the atom positions with intra-molecular +/// distances preserved. +fn unwrap_positions(positions: &[Vector3], box_sides: &Vector3) -> Vec { + if positions.is_empty() { + return vec![]; + } + let mut out = Vec::with_capacity(positions.len()); + out.push(positions[0]); + for r in &positions[1..] { + let prev = *out.last().unwrap(); + let delta = minimum_image(&(r - prev), box_sides); + out.push(prev + delta); + } + out +} + /// Orientationally averaged amplitude ⟨F(q)⟩_Ω = Σ_i f_i · sin(q|r_i − center|)/(q|r_i − center|). /// -/// `center` must be the particle centroid so that |r_i − center| are -/// distances from the scattering centre. This is the isotropic component, -/// equivalent to the l=0 multipole coefficient: 4π|A₀⁰|² = mean_sinc_amplitude². -/// Used to compute β(q) = ⟨F⟩²_Ω / P(q) for the decoupling approximation. +/// When `box_sides` is provided positions are first unwrapped so distances +/// from the centroid are physically correct even if the molecule straddles +/// a periodic boundary. fn mean_sinc_amplitude( q_norm: f64, positions: &[Vector3], form_factors: &[f64], - center: &Vector3, + box_sides: Option<&Vector3>, ) -> f64 { - positions.iter().zip(form_factors).map(|(r, &f)| { + let (pos, center) = match box_sides { + Some(b) => { + let unwrapped = unwrap_positions(positions, b); + let c = centroid(&unwrapped); + (unwrapped, c) + } + None => (positions.to_vec(), centroid(positions)), + }; + pos.iter().zip(form_factors).map(|(r, &f)| { let qr = q_norm * (r - center).norm(); let sinc = if qr < 1e-10 { 1.0 } else { qr.sin() / qr }; f * sinc @@ -151,13 +182,12 @@ impl DirectTransform { let (q, total): (Vec, Vec) = average_duplicates(table).into_iter().unzip(); - let center = centroid(&structure.pos); let beta = q .iter() .zip(&total) .map(|(&qi, &p)| { let ff = eval_atom_form_factors(qi, &structure.ids, &self.formfactors)?; - let mean_amp = mean_sinc_amplitude(qi, &structure.pos, &ff, ¢er); + let mean_amp = mean_sinc_amplitude(qi, &structure.pos, &ff, Some(box_sides)); Ok(if p > 0.0 { mean_amp * mean_amp / p } else { 1.0 }) }) .collect::>>()?; @@ -228,21 +258,17 @@ impl DirectTransform { hydration.push(*h); } - // All three contributions are centered on the same protein centroid so - // that β(q) is frame-of-reference independent (P(q) is already invariant - // under translation because |A|² gains only a common phase factor). - let center = centroid(&structure.pos); for (&qi, &p) in q.iter().zip(&total) { let atom_ff = eval_atom_form_factors(qi, &structure.ids, &self.formfactors)?; - let mean_a = mean_sinc_amplitude(qi, &structure.pos, &atom_ff, ¢er); + let mean_a = mean_sinc_amplitude(qi, &structure.pos, &atom_ff, Some(box_sides)); let mut ev_ff = Vec::new(); excluded_dummies.form_factors_at(qi, &mut ev_ff)?; - let mean_c = mean_sinc_amplitude(qi, excluded_dummies.positions(), &ev_ff, ¢er); + let mean_c = mean_sinc_amplitude(qi, excluded_dummies.positions(), &ev_ff, Some(box_sides)); let mut hyd_ff = Vec::new(); hydration_dummies.form_factors_at(qi, &mut hyd_ff)?; - let mean_b = mean_sinc_amplitude(qi, hydration_dummies.positions(), &hyd_ff, ¢er); + let mean_b = mean_sinc_amplitude(qi, hydration_dummies.positions(), &hyd_ff, Some(box_sides)); let mean_amp = mean_a - volume_scale * mean_c + contrast_density * mean_b; beta.push(if p > 0.0 { mean_amp * mean_amp / p } else { 1.0 }); @@ -332,11 +358,11 @@ mod tests { let center = centroid(&positions); let beta_low = { - let amp = mean_sinc_amplitude(0.001, &positions, &form_factors, ¢er); + let amp = mean_sinc_amplitude(0.001, &positions, &form_factors, None); amp * amp / (2.0 * 2.0) // P(q→0) = (Σ f_i)² }; let beta_high = { - let amp = mean_sinc_amplitude(5.0, &positions, &form_factors, ¢er); + let amp = mean_sinc_amplitude(5.0, &positions, &form_factors, None); let p = 4.0; // |f1|² + |f2|² (incoherent limit) amp * amp / p }; diff --git a/src/trajectory.rs b/src/trajectory.rs index 655cbad..9f9971e 100644 --- a/src/trajectory.rs +++ b/src/trajectory.rs @@ -158,6 +158,7 @@ pub fn trajectory_average( let mut mean = Intensity::default(); let mut structure = structure.clone(); + let single_protein = structure.pos.len(); for (frame, weight) in reader { structure.set_positions(&frame)?; @@ -166,15 +167,34 @@ pub fn trajectory_average( if mean.q.is_empty() { mean.q = curve.q.clone(); mean.total = vec![0.0; curve.total.len()]; + // β(q) is only meaningful for a single-protein box + if frame.natoms() == single_protein && curve.beta.is_some() { + mean.beta = Some(vec![0.0; curve.total.len()]); + } } for (acc, v) in zip(&mut mean.total, &curve.total) { *acc += weight * v; } + // Accumulate mean sinc amplitude (= sqrt(β·P)) rather than β itself. + // After the loop, β = (Σ w·)² / (Σ w·P) is computed from these. + if let (Some(acc), Some(src_beta)) = (mean.beta.as_mut(), curve.beta.as_ref()) { + for ((a, &beta_t), &p_t) in acc.iter_mut().zip(src_beta).zip(&curve.total) { + *a += weight * (beta_t * p_t).sqrt(); + } + } + #[cfg(feature = "cli")] progress.inc(1); } + // Convert accumulated mean amplitudes to β = ² /

+ if let Some(mean_amps) = mean.beta.as_mut() { + for (amp, &p) in mean_amps.iter_mut().zip(&mean.total) { + *amp = if p > 0.0 { (*amp * *amp) / p } else { 1.0 }; + } + } + Ok(mean) }