Pripps — SAXS Intensity from Molecular Structures

Pripps computes small-angle X-ray scattering (SAXS) intensity profiles from PDB, XYZ, or GRO structures, with optional solvent corrections (excluded volume + hydration shell). It is fast enough for iterative fitting against experimental data.

Quick start (Python)

Installation

Install straight from GitHub with uv. The extension is built from source, so a Rust toolchain and Python ≥ 3.9 are required:

uv pip install "git+https://github.com/mlund/pripps.git#subdirectory=pripps-py"

To also get the example structures and form-factor assets used below, clone the repo and install in editable mode instead:

git clone https://github.com/mlund/pripps.git
cd pripps/pripps-py
uv pip install -e .          # with uv (recommended)
# pip install maturin && maturin develop   # without uv

Computing a SAXS profile

import pripps

pr = pripps.Pripps(atomfile="assets/foxs_formfactors.yaml")
model = pr.multipole_model(
    structure_file="tests/lysozyme/6lyz.pdb",
    qmin=0.0,
    qmax=0.5,
    nq=500,
    excluded="foxs",          # or pripps.Excluded.Foxs — 16π Fraser excluded volume
    hydration="sasa",         # or pripps.Hydration.Sasa
    probe=1.8,  # SASA probe radius (Å)
    bulk_electron_density=0.334,  # e/ų, pure water
)

# Fit knobs (volume_scale, contrast_density) are applied here, not at build time.
res = model.intensity(volume_scale=1.0, contrast_density=0.03)   # -> pripps.Intensity
print(res.q[:3], res.total[:3])

intensity() returns a typed Intensity with NumPy-array attributes:

Attribute	Description
.q	Scattering vector magnitudes (Å⁻¹)
.total	Total intensity I(q)
.atoms	Protein-only (vacuum) contribution
.excluded	Excluded-volume contribution
.hydration	Hydration-shell contribution

Cached models always populate .atoms / .excluded / .hydration (zero where a contribution is disabled). They are None only when there is no per-contribution breakdown — currently trajectory averages (the one-shot Pripps.intensity(..., xtc=...)).

Fitting against experimental data

multipole_model does the expensive per-structure setup once and caches the intermediate amplitudes; model.intensity(volume_scale, contrast_density) then only applies the fit knobs — typically under a millisecond. There are two ways to fit:

import numpy as np

# experimental data: q (Å⁻¹), I(q), σ(q)
q_exp, i_exp, sigma = np.loadtxt("experimental.dat", comments="#", unpack=True)

# (a) built-in fit — returns a typed FitResult:
fit = model.fit(q_exp, i_exp, sigma)
print(fit.volume_scale, fit.contrast_density, fit.chi2, fit.scale)

# (b) manual fit with scipy for full control:
from scipy.optimize import minimize
def chi2(params):
    out = model.intensity(*params)
    theory = np.interp(q_exp, out.q, out.total)
    scale = np.sum(theory * i_exp / sigma**2) / np.sum((theory / sigma) ** 2)
    return np.sum(((scale * theory - i_exp) / sigma) ** 2) / len(q_exp)
result = minimize(chi2, x0=[1.0, 0.03],
                  bounds=[(0.95, 1.12), (-2.0, 4.0)], method="L-BFGS-B")

debye_model (exact, for small structures) and direct_model (explicit discrete-q) accept the same solvent arguments and expose the same intensity / fit / q_values interface.

A full working example is in pripps-py/examples/fit_lysozyme.py, which fits lysozyme against the shipped tests/lysozyme/lyzexp.dat and reproduces FoXS's χ²/N ≈ 0.20 (demonstrating both fitting styles):

python pripps-py/examples/fit_lysozyme.py

Trajectory averaging

To average I(q) over an MD trajectory, use the one-shot Pripps.intensity (the cached models are single-structure). The model= argument selects the calculator; trajectory averages carry no per-contribution breakdown:

avg = pr.intensity(
    structure_file="protein.pdb",   # topology / template frame
    model="multipole",              # or "debye" / "direct", or a pripps.Scheme
    xtc="trajectory.xtc",
    stride=10,                      # use every 10th frame
    weights="weights.dat",          # optional per-frame weights
    excluded="foxs",
    hydration="sasa",
    contrast_density=0.03,
)
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

$$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:

$$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

$$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

$$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

$$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,

$$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.

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$.

pripps -a atoms.yaml -i box.xyz -x traj.xtc --num-molecules 100 direct -p 100

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
excluded	"disabled", "fraser", "foxs", "grid", "voronoi", "per-kind"	Excluded-volume model
hydration	"disabled", "sasa", "grid", "voronoi"	Hydration-shell model
volume_scale	float (default 1.0)	Scales excluded volume (PepsiSAXS r₀). In Python this is a fit knob passed to model.intensity() / model.fit(), not the model builder
excluded_spacing	float (default 1.0)	Grid spacing (Å) for excluded="grid"
contrast_density	float (default 0.03)	Fractional hydration-layer excess density (c₂ in A − c₁·C + c₂·B) — pepsi's d_rho / ρ_water, ~0.1 for a 10 % excess layer. In Python a fit knob on model.intensity() / model.fit()
probe	float (default: 1.8 for SASA/grid, 6.0 for voronoi)	Water probe radius in Å
spacing	float (default 3.10516)	Grid spacing (Å) for hydration="grid"; the bulk-water linear spacing, so each cell holds ~1 water molecule
thickness	float (default 3.0)	Shell thickness (Å) for hydration="grid"
bulk_electron_density	float (default 0.334)	Bulk solvent density in e/ų

Excluded-volume models

The excluded volume accounts for bulk solvent displaced by the protein. Each model places dummy scatterers whose contribution is subtracted from the atomic scattering via |A − volume_scale·C + contrast_density·B|².

Model	Description	Requires
fraser	One Gaussian dummy per atom, sized by displaced_volume from the YAML. Dummy form factor decays faster with q → smaller excluded-volume subtraction at high q. Matches PepsiSAXS / CRYSOL.	displaced_volume per atom kind
foxs	Same formula as fraser but with the FoXS 16π denominator: dummy form factor decays more slowly with q → larger excluded-volume subtraction at high q. Matches FoXS.	displaced_volume per atom kind
grid	Fills the molecular envelope with uniform grid cells, each carrying a Fraser Gaussian sized by cell volume (spacing³). Independent of per-atom displaced volumes.	radius per atom kind
voronoi	Per-atom Fraser Gaussians sized by individual Voronoi cell volumes. Derives displaced volumes geometrically — no tabulated volumes needed. Achieves χ²≈0.20 on lysozyme, comparable to fraser.	radius per atom kind
per-kind	Uses per-kind excluded-solvent form factors stored as excluded_solvent beside each atom/residue formfactor in the same YAML file.	excluded_solvent for referenced atom kinds

The fraser and foxs models use per-atom displaced volumes from the form-factor YAML (derived from radius as (4/3)πr³ when displaced_volume is not explicitly set). The grid and voronoi models need only atomic radii.

Hydration-shell models

The hydration shell accounts for the first water layer around the protein, which has a higher electron density than bulk solvent.

All three schemes weight their dummies so weight × F_water(q) represents the scattering from bulk-density water in the dummy's associated volume — so contrast_density has a consistent physical meaning (fractional excess density above bulk) across schemes. Pepsi's d_rho in e/ų maps as contrast_density ≈ d_rho / ρ_water.

Model	Description
sasa	FoXS-style: one dummy per atom, weight = SASA area × 3 Å / 30 ų (accessible-area × shell-thickness → bulk-water count). Best with atomic-resolution YAMLs.
grid	PepsiSAXS-style: water dummies on a Cartesian grid in a shell around the molecular envelope, weight = waters per cell at bulk density.
voronoi	Explicit solvent spheres on the molecular surface via voronota-ltr, weight = cell volume / 30 ų. Dummies sit at actual water-accessible positions (not atom centers). Needs only radius. Larger probe (default 6 Å) smooths the granular shell into an envelope-like layer that avoids spurious first-minimum interference.

Form factor definitions

Atom definitions (atomic or coarse-grained) are specified in a YAML file. Each entry maps a species name to a form factor model and optional physical properties:

H: { formfactor: !Constant 1.0, radius: 1.0 }
N: { ff: !Table [ [0.0, 8.9], [0.1, 7.8], ...] }
C: { ff: !Gaussian { qrange: [0.0, 1.0], ab: [[0.1, 2.0]], c: 0.4 },
     radius: 1.577 }
# United atom with explicit displaced_volume (includes implicit H):
NH: { radius: 1.009, displaced_volume: 7.606,
      ff: !United { center: N, orbit: [[H, 1.009, 1]] } }
CA: { ff: !Alias C }
# Linear combination of other entries (no geometry), F(q) = Σ wᵢ Fᵢ(q):
H2O: { ff: !Linear [[O, 1.0], [H, 2.0]] }
ALA: { ff: !Polynomial { qrange: [0.0, 0.75], coeff: [8.9, 0.0, 5.8, ...]},
       excluded_solvent: !Polynomial { qrange: [0.0, 0.75], coeff: [7.1, 0.0, 4.2, ...]} }

Supported form factor models:

Model	Keyword	Formula
Constant	!Constant	F(q) = c
Polynomial	!Polynomial	F(q) = a + bq + cq² + …
Gaussian	!Gaussian	F(q) = c + Σ aᵢ exp(−bᵢs²), s = q/(4π)
United atom	!United	F(q) = Fc(q) + Σ wᵢ Fᵢ(q) sinc(q rᵢ)
Linear combination	!Linear	F(q) = Σ wᵢ Fᵢ(q)
Tabulated	!Table	Array of [q, F] pairs
Alias	!Alias	Reference to another entry

Property	Unit	Used by	Auto-derived?
radius	Å	SASA hydration, Fraser excluded volume	No — must be specified
displaced_volume	ų	Fraser excluded volume	Yes — defaults to (4/3)πr³ from radius if omitted
excluded_solvent	form factor	per-kind excluded volume	No — optional; required per species when using --excluded per-kind

When radius is present but displaced_volume is omitted, the displaced volume is automatically derived as the sphere volume (4/3)πr³. Specify displaced_volume explicitly only for united-atom groups (e.g. NH, SH) where the displaced volume includes implicit hydrogens that the effective radius does not account for.

The optional excluded_solvent field stores a second form factor (any of the models above) for the bulk solvent each species displaces, sitting beside its in-vacuo formfactor in the same entry. It is read only by the per-kind excluded-volume model (--excluded per-kind, or excluded="per-kind" from Python), which subtracts this form factor — 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/pepsi_atoms.yaml	Atomic Waasmaier-Kirfel five-Gaussians + PepsiSAXS 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	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

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 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.

For atomic-resolution PDB files, use foxs_formfactors.yaml with --excluded foxs — this reproduces FoXS's χ²/N ≈ 0.29 on lysozyme (unfitted) and ≈ 0.20 (fitted with hydration).

PDB atom mapping: for PDB files, pripps uses residue context to map each heavy atom to its united-atom form-factor type (e.g. ALA:CB → CH₃, LYS:NZ → NH₃), folding in implicit hydrogen scattering — same convention as FoXS.

Performance

The multipole expansion is linear in atom count, so it beats the O(N²) Debye schedule above a few hundred atoms and stays practical for atomic-resolution PDBs up into the hundreds of thousands of atoms.

• Auto-selected truncation order per q-point. Low q uses only a handful of multipoles, high q uses a convergent upper bound. The default is conservative on non-globular shapes (elongated complexes, stacked rings, multi-domain assemblies) where radius-of-gyration-based heuristics in other tools can silently under-converge at high q.
• Auto-selected q-grid density at the Nyquist rate (2·L + 1 points) when --nq is omitted. The theory I(q) is band-limited by the multipole truncation, so sampling any denser just computes redundant points — experimental q values are reached via natural cubic-spline interpolation on the theory grid. Shannon-Nyquist sampling applied to SAXS is discussed by Moore, J. Appl. Cryst. 13, 168-175 (1980).
• Parallel across CPU cores and SIMD-vectorised. No tuning needed; the workload scales with hardware.
• Cached fit setup. multipole_model amortises the per-structure work so fit loops spend their time in the cheap recombination, not the expensive scattering calculation.

Illustrative wall-clock on a 90 672-atom structure, 100 q-points, full atomic-resolution form factors (vacuum):

Hardware	wall
10-core Apple Silicon laptop	~4 s
64-thread AMD Threadripper 3970X (-C target-cpu=native)	~9 s

Tuning knobs

• --lmax <N> (CLI) or l_max=N (Python) overrides the auto-selected truncation order. Useful for reproducing a fixed reference, timing sweeps, or trading a little high-q accuracy for extra speed on very elongated structures — the default is correct but conservative for those.
• --nq <N> (CLI) overrides the auto-selected q-grid density. The default 2·L_max + 1 is enough to represent the band- limited theory without aliasing; pass a larger value to oversample the output curve for visualisation.

Command-line interface (Rust)

Installing

Requires a Rust toolchain:

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

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

# With solvent correction — CSV output with 5 columns
pripps -a atoms.yaml -i structure.pdb -o intensity.csv \
    multipole --excluded fraser --hydration sasa
pripps -a atoms.yaml -i structure.pdb -o intensity.csv \
    debye --excluded fraser --hydration sasa

# Fit against experimental data (box-constrained Nelder-Mead over
# c₁ ∈ [0.5, 2.0] and c₂ ∈ [0.0, 0.1]; warns if the fit hits a bound)
pripps -a assets/foxs_formfactors.yaml -i structure.pdb \
    -f experimental.dat -o fit.csv \
    multipole --excluded foxs --hydration sasa
pripps -a assets/foxs_formfactors.yaml -i structure.pdb \
    -f experimental.dat -o fit.csv \
    debye --excluded foxs --hydration sasa
pripps -a assets/foxs_formfactors.yaml -i structure.pdb \
    -f experimental.dat -o fit.csv \
    direct --excluded foxs --hydration sasa

# Trajectory averaging — use `direct` so each frame's periodic box
# is honoured (multipole ignores PBC).
pripps -a atoms.yaml -i structure.pdb -x trajectory.xtc direct

Choosing a scheme

Scheme	Best for	Notes
multipole	Single structures, fitting loops	Fastest for large atomic PDBs; no periodic boundaries.
direct	MD trajectories with periodic box	Evaluates I(q) on discrete Miller-index q-vectors; supports the same solvent fit knobs for single-frame fits.
debye	Small clusters, reference checks	Exact pair-distance sum with the same solvent and fit knobs as multipole; quadratic in atoms, fine under a few hundred.

Output is CSV: q,total for vacuum, or q,total,atoms,excluded,hydration when solvent is active. 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
FF_TABLE	If set, export resolved form factors as tabulated_formfactors.yaml
RUST_LOG	Log level: Debug, Trace, Warn, Info