feat: add beta(q) decoupling parameter to Multipole and Direct CSV output by pinedaps · Pull Request #4 · mlund/pripps

pinedaps · 2026-05-05T07:23:28Z

Summary

Adds β(q) = ⟨F⟩²_Ω / P(q) as a new beta column in the Multipole and Direct CSV outputs
Covers vacuum and solvent-corrected paths for both schemes
Trajectory averaging correctly accumulates ⟨F⟩ and P separately, then computes β = ⟨F⟩²/⟨P⟩
β output is suppressed for multi-protein boxes (physically meaningless there)
Debye path unaffected

Physics

β(q) 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)

Multipole path: ⟨F(q)⟩_Ω = √(4π)·A₀⁰(q), so β(q) = 4π|A₀⁰|²/P(q)

Direct path: ⟨F(q)⟩_Ω = Σ_i f_i·sinc(q|r_i − centroid|), so β(q) = ⟨F⟩²/P(q)

Both routes give the same physical quantity. β(q) → 1 at q = 0 and decreases with q for asymmetric proteins.

Trajectory averaging (Direct + XTC)

For a trajectory the correct formula is β = (⟨F⟩_conf)² / ⟨P⟩_conf. Averaging per-frame β_t = ⟨F⟩_t² / P_t is biased when P_t fluctuates (e.g. a rotating protein sampled through a small number of Miller directions). Instead:

Per frame: accumulate mean sinc amplitude A_t = Σ_i f_i sinc(q|r_i − centroid|)
After loop: β = (⟨A⟩)² / ⟨P⟩

Single-protein detection: β is only written when the XTC frame atom count equals the topology atom count. For multi-protein boxes P(q) includes inter-protein terms, making β meaningless.

PBC fix

complex_amplitude used for I(q) is immune to PBC wrapping (exp(iq·r) is periodic for integer Miller indices). But mean_sinc_amplitude uses physical distances |r_i − centroid|, which are corrupted when the protein crosses a box boundary. Fix: unwrap the molecule using the minimum-image convention before computing the centroid and distances. Per-frame β now stays ≈ 0.9999 at low q across all 100 frames instead of dropping to 0.03–0.24 in frames where the protein straddled the boundary.

CSV output

Vacuum direct/multipole (single protein, XTC or static structure):

q,total,beta
0.027318,2965661.190826,0.999856
0.038634,2817950.481171,0.999562
...

Solvent-corrected path:

q,total,atoms,excluded,hydration,beta

Test plan

Clean build, no warnings
Static PDB/XYZ: β(q→0) ≈ 1.000, decreases with q — matches multipole
Trajectory (4LZT, 100 frames, rotating protein): β(q=0.027) ≈ 0.9999 after PBC fix
Direct and Multipole β values agree at low q on same protein
Multi-protein trajectory: beta column absent (suppressed correctly)
Unit test: direct_beta_near_one_at_low_q_and_decays
Debye CSV output unchanged

🤖 Generated with Claude Code

…CSS alias - `--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 <noreply@anthropic.com>

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<Vec<f64>>` 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 <noreply@anthropic.com>

pinedaps · 2026-05-05T07:30:46Z

Follow-up if approved: can we also expose mean_amplitude_sq in the Python bindings (see original todo_mean_amplitude_sq.txt for the implementation sketch).

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 <noreply@anthropic.com>

…method 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 <noreply@anthropic.com>

Copilot

Pull request overview

Adds the decoupling parameter β(q) as an optional beta column to Pripps’ Multipole and Direct outputs, including both vacuum and solvent-corrected paths, and introduces PBC-unwrapping support for the Direct mean-sinc amplitude used in β.

Changes:

Extend Intensity with beta: Option<Vec<f64>> and emit beta for Multipole and Direct computations.
Update CSV writing to conditionally include beta (with/without solvent contribution columns).
Add PBC-aware mean-sinc amplitude helpers in the Direct implementation and adjust trajectory averaging logic for β.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 5 comments.

Show a summary per file

File	Description
`src/trajectory.rs`	Adds trajectory-time accumulation logic intended to compute β via averaged amplitudes.
`src/multipole/mod.rs`	Computes β from the monopole coefficient (vacuum + solvent-corrected multipole).
`src/lib.rs`	Adds `beta` to `Intensity` and applies a default box for non-trajectory runs when none is provided.
`src/explicit.rs`	Implements mean-sinc amplitude + PBC unwrapping and computes β for Direct (vacuum + solvent). Adds a unit test.
`src/cli.rs`	Makes `--box` optional for Direct and extends CSV output to include the `beta` column when present.
`assets/single_bead.yaml`	Adds a `CSS` alias entry to the form-factor YAML.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+        // Accumulate mean sinc amplitude (= sqrt(β·P)) rather than β itself.
+        // After the loop, β = (Σ w·<F>)² / (Σ 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();
+            }


@@ -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()]);
+            }


+/// 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<Vector3> {
+    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
+}


+fn mean_sinc_amplitude(
+    q_norm: f64,
+    positions: &[Vector3],
+    form_factors: &[f64],
+    box_sides: Option<&Vector3>,
+) -> f64 {
+    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
+    }).sum()
+}


+        let form_factors = vec![1.0, 1.0];
+        let center = centroid(&positions);
+


Sebastian Pineda and others added 2 commits May 4, 2026 14:08

pinedaps changed the title ~~feat: add beta(q) decoupling parameter to multipole CSV output~~ feat: add beta(q) decoupling parameter to Multipole and Direct CSV output May 6, 2026

mlund requested a review from Copilot May 29, 2026 09:54

Copilot started reviewing on behalf of mlund May 29, 2026 09:54 View session

Copilot AI reviewed May 29, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

feat: add beta(q) decoupling parameter to Multipole and Direct CSV output#4

feat: add beta(q) decoupling parameter to Multipole and Direct CSV output#4
pinedaps wants to merge 4 commits into
masterfrom
feat/beta-q-output

pinedaps commented May 5, 2026 •

edited

Loading

Uh oh!

pinedaps commented May 5, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

		let form_factors = vec![1.0, 1.0];
		let center = centroid(&positions);

Conversation

pinedaps commented May 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Physics

Trajectory averaging (Direct + XTC)

PBC fix

CSV output

Test plan

Uh oh!

pinedaps commented May 5, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

pinedaps commented May 5, 2026 •

edited

Loading