Grid fitting

Cargo.toml

@@ -26,6 +26,11 @@ serde = { version = "1.0", features = ["derive"] }
 yaml_serde = "0.10"
 voronota-ltr = { git = "https://github.com/mlund/voronota-ltr" }
 
+# GPU (wgpu) Direct Fourier transform — see `gpu` feature below.
+wgpu = { version = "29", optional = true }
+bytemuck = { version = "1", features = ["derive"], optional = true }
+pollster = { version = "0.4", optional = true }
+
 [dev-dependencies]
 approx = "0.5"
 criterion = "0.5"
@@ -34,6 +39,8 @@ tempfile = "3"
 [features]
 default = ["cli"]
 cli = ["dep:clap", "dep:pretty_env_logger", "dep:indicatif"]
+# GPU-accelerated Direct Fourier transform via wgpu/WGSL.
+gpu = ["dep:wgpu", "dep:bytemuck", "dep:pollster"]
 
 [[bin]]
 name = "pripps"

README.md

@@ -129,6 +129,107 @@ avg = pr.intensity(
 print(avg.q, avg.total)
 ```
 
+### Effective structure factor (`--num-molecules`)
+
+For a periodic box of $`N`$ identical molecules — a concentrated solution
+simulated with the `direct` scheme, optionally averaged over an XTC trajectory —
+`--num-molecules N` (CLI) or `num_molecules=N` (Python) separates the
+single-molecule **form factor** $`P(q)`$ from the inter-molecular **structure
+factor**.
+
+Each molecule $`m`$ scatters with amplitude
+
+```math
+F_m(\mathbf{q}) = \sum_{i \in m} f_i(q)\, e^{\,i\,\mathbf{q}\cdot\mathbf{r}_i},
+```
+
+where $`f_i(q)`$ is the atomic form factor and $`\mathbf{r}_i`$ the atom position.
+The box intensity is the squared modulus of the total amplitude, which splits
+into intra-molecular (self) and inter-molecular (cross) terms:
+
+```math
+I(\mathbf{q}) = \left|\sum_{m=1}^{N} F_m(\mathbf{q})\right|^2
+            = \underbrace{\sum_{m} \left|F_m\right|^2}_{\text{self}}
+            + \underbrace{\sum_{m \neq m'} F_m F_{m'}^{*}}_{\text{cross}} .
+```
+
+Writing $`\langle\,\cdot\,\rangle`$ for the average over scattering-vector
+orientations at fixed $`q = |\mathbf{q}|`$ (and over trajectory frames), the form
+factor and the **effective** (apparent) structure factor are
+
+```math
+P(q) = \frac{1}{N}\left\langle \sum_{m} \left|F_m(\mathbf{q})\right|^2 \right\rangle
+     = \bigl\langle \left|F(\mathbf{q})\right|^2 \bigr\rangle_{\Omega},
+\qquad
+S_\mathrm{eff}(q) = \frac{\langle I(\mathbf{q})\rangle}{N\,P(q)}
+                  = \frac{\bigl\langle \left|\sum_m F_m\right|^2 \bigr\rangle}
+                         {\bigl\langle \sum_m \left|F_m\right|^2 \bigr\rangle} .
+```
+
+pripps emits exactly these two columns, `p` and `s_eff`. $`S_\mathrm{eff}(q)`$ is
+the quantity a scattering experiment measures — the box intensity normalised by
+$`N`$ independent molecules — and is obtained directly from the self/cross split,
+so it is **exact**: no assumption about particle shape or correlations is made.
+
+**Relation to the true structure factor.** The thermodynamic centre-to-centre
+structure factor, the Fourier transform of the centre–centre pair correlation
+$`g(r)`$, is
+
+```math
+S(q) = 1 + \frac{1}{N}\left\langle \sum_{m \neq m'}
+        e^{\,i\,\mathbf{q}\cdot(\mathbf{R}_m - \mathbf{R}_{m'})} \right\rangle,
+```
+
+with $`\mathbf{R}_m`$ the molecular centre. For non-spherical molecules
+$`S_\mathrm{eff} \neq S`$. Under the *decoupling approximation* — orientation
+uncorrelated with position — they are related by
+
+```math
+S_\mathrm{eff}(q) = 1 + \beta(q)\,\bigl[S(q) - 1\bigr],
+\qquad
+\beta(q) = \frac{\bigl|\langle F(\mathbf{q})\rangle_{\Omega}\bigr|^2}
+                {\bigl\langle |F(\mathbf{q})|^2 \bigr\rangle_{\Omega}},
+```
+
+where the decoupling parameter obeys $`0 \le \beta(q) \le 1`$, with $`\beta = 1`$ for
+centrosymmetric (e.g. spherical) particles, giving $`S_\mathrm{eff} = S`$. The true
+$`S(q)`$ follows by inversion,
+
+```math
+S(q) = 1 + \frac{S_\mathrm{eff}(q) - 1}{\beta(q)},
+```
+
+which is ill-conditioned near the nodes of $`\beta(q)`$; pripps therefore reports
+$`P`$ and $`S_\mathrm{eff}`$ and leaves $`\beta(q)`$ and the (judgement-dependent)
+inversion to the user. See Kotlarchyk & Chen, *J. Chem. Phys.* **79**, 2461
+(1983), [doi:10.1063/1.446055](https://doi.org/10.1063/1.446055).
+
+**Notes.**
+- The amplitudes use the box-commensurate vectors $`\mathbf{q} = 2\pi p\,(h,k,l)/L`$,
+  for which $`e^{\,i\,\mathbf{q}\cdot\mathbf{r}}`$ is unchanged under
+  $`\mathbf{r}\to\mathbf{r}+\mathbf{L}`$. Both $`I`$ and the self term are thus
+  PBC-invariant and need no coordinate unwrapping; the orientational average is
+  built from the sampled $`(h,k,l)`$ directions, and the trajectory average is a
+  ratio of frame-averaged numerator and denominator.
+- $`P(q)`$ is the **at-concentration** form factor — the molecules as they actually
+  are in the crowded box. For a flexible molecule it depends on concentration,
+  unlike the dilute-limit $`P`$ an experiment divides out, so comparing the two
+  reveals conformational change. It is on an absolute scale,
+  $`P(0) = \bigl(\sum_i f_i\bigr)^2`$, not normalised to $`P(0)=1`$.
+
+```sh
+pripps -a atoms.yaml -i box.xyz -x traj.xtc --num-molecules 100 direct -p 100
+```
+
+```python
+out = pr.intensity("box.xyz", model="direct", box_sides=200.0, num_molecules=100)
+print(out.s_eff, out.p)
+```
+
+The box atom count must be a multiple of $`N`$ (contiguous, identical molecules),
+and the run must be vacuum. The `gpu` feature accelerates this path too; XTC
+averaging assumes a fixed box (constant volume).
+
 ### Solvent models
 
 | Parameter | Values | Description |
@@ -234,27 +335,39 @@ scaled by `volume_scale` — instead of building a Gaussian dummy from
 `radius`/`displaced_volume`. Every species present in the structure
 must define it when `--excluded per-kind` is selected.
 
+#### Merging multiple files
+
+Several form-factor files can be supplied and are merged in order — repeat
+`--atomfile` on the CLI, or pass a list to the Python constructor — letting
+you compose, e.g., a coarse-grained residue file with an atomic one. If the
+same species appears in more than one file the **last file listed wins** and
+a warning is logged. Merging happens before symbolic resolution, so an
+`!Alias` / `!United` / `!Linear` entry in one file may reference a species
+defined in another.
+
 ### Shipped form-factor files
 
 | File | Description | Best for |
 |------|-------------|----------|
 | `assets/foxs_formfactors.yaml` | Cromer-Mann (Int. Tables Vol. C) atomic form factors, FoXS-compatible volumes and SASA radii | Atomic PDBs with `--excluded foxs` |
-| `assets/gaussian_atoms.yaml` | Mixed table: residue single-bead polynomials plus atomic/united-atom Gaussians (CH/CH₂/NH/…) | Atomic PDBs with PepsiSAXS-style solvent models |
-| `assets/stovgaard2010.yaml` | Per-residue form-factor centroids (one-body model, q ∈ [0, 0.75] Å⁻¹) from Stovgaard et al. 2010, <https://doi.org/10.1186/1471-2105-11-429> | Residue / Cα single-bead models |
+| `assets/pepsi_atoms.yaml` | Atomic [Waasmaier-Kirfel](https://doi.org/10.1107/S0108767394013292) five-Gaussians + [PepsiSAXS](https://doi.org/10.1107/S2059798317005745) united-atom groups (CH/CH₂/NH/…), K = 4π | Atomic PDBs with PepsiSAXS-style solvent models |
+| `assets/stovgaard2010.yaml` | Per-residue form-factor centroids (one-body model, q ∈ [0, 0.75] Å⁻¹) from [Stovgaard et al. 2010](https://doi.org/10.1186/1471-2105-11-429) | Residue / Cα single-bead models |
 | `assets/poly_amino_acid.yaml` | Amino-acid bead polynomials with a separate `excluded_solvent` polynomial per bead | Cα / single-bead with `--excluded per-kind` |
 | `assets/poly_amino_acid_combined.yaml` | Amino-acid bead polynomials with the excluded solvent already folded into the form factor | Cα / single-bead with `--excluded disabled` |
 | `assets/poly_martini.yaml` | Martini-3 bead polynomials with a separate `excluded_solvent` polynomial per bead | Martini-3 CG with `--excluded per-kind` |
 | `assets/poly_martini_combined.yaml` | Martini-3 bead polynomials with the excluded solvent already folded into the form factor | Martini-3 CG with `--excluded disabled` |
 
-Gaussian files use b values in the standard s² = (q/4π)² convention
-(International Tables / Cromer-Mann). The atomic coefficients differ
-(Cromer-Mann vs Stovgaard-2010) and the volume/radius tables match
-their respective tools.
+Files using `!Gaussian` form factors (`foxs_formfactors.yaml`,
+`pepsi_atoms.yaml`) give b values in the standard s² = (q/4π)² convention
+(International Tables / Cromer-Mann), so published exponents can be used
+directly. Their atomic coefficients come from different parametrizations
+(Cromer-Mann for FoXS vs Waasmaier-Kirfel for Pepsi) and the volume/radius
+tables match their respective tools.
 
 **Don't mix files with non-matching `--excluded` flags.** Each
 YAML + excluded-volume convention is a self-consistent bundle:
 `foxs_formfactors.yaml` was calibrated jointly with FoXS's
-`K = 16π` Gaussian, while `gaussian_atoms.yaml` was calibrated with
+`K = 16π` Gaussian, while `pepsi_atoms.yaml` was calibrated with
 PepsiSAXS's `K = 4π`. Swapping the denominator gives atomic and
 excluded contributions that were never fit together and yields
 subtly biased intensities.
@@ -322,6 +435,9 @@ Requires a [Rust toolchain](https://www.rust-lang.org/tools/install):
 
 ```sh
 cargo install --git https://github.com/mlund/pripps.git
+
+# …or with GPU acceleration for the `direct` scheme (see below)
+cargo install --git https://github.com/mlund/pripps.git --features gpu
 ```
 
 ### Usage
@@ -332,6 +448,9 @@ pripps --help
 # Debye formula (exact, O(N²))
 pripps -a atoms.yaml -i structure.pdb debye
 
+# Merge several form-factor files (last listed wins on duplicate species)
+pripps -a base.yaml -a overrides.yaml -i structure.pdb debye
+
 # Multipole expansion (linear in atoms, truncation auto-selected)
 pripps -a atoms.yaml -i structure.pdb multipole
 pripps -a atoms.yaml -i structure.pdb multipole --lmax 68  # fixed order
@@ -373,6 +492,14 @@ With `--fit`, output is `q,i_exp,i_fit,sigma` at experimental
 q points. The experimental file is auto-detected as CSV or
 whitespace-separated (3 columns: q, I, σ).
 
+### GPU acceleration (`gpu` feature)
+
+Building with `--features gpu` runs the **`direct`** scheme on the GPU,
+which speeds up trajectory averaging substantially (measured ~15× over
+the CPU on an Apple M4). It applies to the vacuum transform; runs with a
+solvent model, or without a usable GPU, fall back to the CPU
+automatically. Pass `--cpu` to opt out.
+
 ### Environment variables
 
 | Variable | Description |