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

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

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<f64>,
         #[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 {
@@ -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}"
                 )?;
             }
         }

src/explicit.rs

@@ -44,6 +44,64 @@ 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<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
+}
+
+/// Orientationally averaged amplitude ⟨F(q)⟩_Ω = Σ_i f_i · sin(q|r_i − center|)/(q|r_i − center|).
+///
+/// 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],
+    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()
+}
+
+/// Geometric centroid of a set of positions.
+fn centroid(positions: &[Vector3]) -> Vector3 {
+    if positions.is_empty() {
+        return Vector3::zeros();
+    }
+    positions.iter().sum::<Vector3>() / positions.len() as f64
+}
+
 /// Calculate I(q) using discrete q vectors
 ///
 /// References:
@@ -123,9 +181,21 @@ impl DirectTransform {
             .collect::<Result<Vec<(f64, f64)>>>()?;
 
         let (q, total): (Vec<f64>, Vec<f64>) = average_duplicates(table).into_iter().unzip();
+
+        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, Some(box_sides));
+                Ok(if p > 0.0 { mean_amp * mean_amp / p } else { 1.0 })
+            })
+            .collect::<Result<Vec<f64>>>()?;
+
         Ok(Intensity {
             q,
             total,
+            beta: Some(beta),
             ..Default::default()
         })
     }
@@ -178,12 +248,30 @@ 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);
+        }
+
+        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, 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, 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, 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 });
         }
 
         Ok(Intensity {
@@ -194,6 +282,7 @@ impl DirectTransform {
                 excluded,
                 hydration,
             }),
+            beta: Some(beta),
         })
     }
 }
@@ -257,4 +346,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, 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, None);
+            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}");
+    }
 }

src/lib.rs

@@ -48,11 +48,18 @@ pub(crate) type Vector3 = nalgebra::Vector3<f64>;
 /// `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 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<f64>,
     pub total: Vec<f64>,
     pub contributions: Option<Contributions>,
+    pub beta: Option<Vec<f64>>,
 }
 
 /// Per-contribution intensity curves emitted by the three-contribution
@@ -212,6 +219,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());

src/multipole/mod.rs

@@ -340,16 +340,22 @@ fn vacuum_intensity(
         l_max_per_q,
     )?;
 
-    let total: Vec<f64> = coeffs_per_q
+    let (total, beta): (Vec<f64>, Vec<f64>) = 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,24 +473,28 @@ 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,
             &self.l_max_per_q,
         )
         .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),
         }
     }
 

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·<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();
+            }
+        }
+
         #[cfg(feature = "cli")]
         progress.inc(1);
     }
 
+    // Convert accumulated mean amplitudes to β = <F>² / <P>
+    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)
 }