Metadata-Version: 2.4
Name: SincCity
Version: 0.2.0
Summary: Generic Scattering Calculator: I(q) of a real-space structure (PDB, GRO, SLD or VTK)
Author: SincCity contributors
License-Expression: BSD-3-Clause
Project-URL: Documentation, https://easyscience.github.io/SincCity/
Project-URL: Source, https://github.com/easyscience/SincCity
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=2.0
Requires-Dist: periodictable>=1.6
Provides-Extra: plot
Requires-Dist: plotly>=5.0; extra == "plot"
Provides-Extra: sasdata
Requires-Dist: sasdata>=0.8; extra == "sasdata"
Provides-Extra: docs
Requires-Dist: mkdocs<2,>=1.6; extra == "docs"
Requires-Dist: mkdocs-material>=9.6; extra == "docs"
Requires-Dist: mkdocstrings[python]>=1.0; extra == "docs"
Provides-Extra: dev
Requires-Dist: sinccity[plot,sasdata]; extra == "dev"
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: ruff>=0.14; extra == "dev"
Dynamic: license-file

# SincCity: Generic Scattering Calculator

The orientationally averaged I(q) of a real-space structure, by the Debye
formula: atoms from a **PDB** or GROMACS **GRO** file, pixels from a
SasView **SLD** file, or the cells of a legacy ASCII **VTK** mesh. It is the
nuclear, 1D part of SasView's *Generic Scattering Calculator*, as a plain
Python package for scripts and notebooks. There are no 2D, magnetic or
polarisation options.

The name comes from the sinc function, sinc x = sin x / x, at the heart of
the Debye formula: every pair of points a distance r apart adds its weight
times sinc(q r) to the intensity.

`sinccity` is developed alongside [IQuana](https://github.com/ai4se1dk/IQuana)
but does not depend on it. The common interface is `profile.to_data1d()`, which
returns a sasdata `Data1D` that `IQuana.set_data` can ingest.

## Installation

```bash
pip install -e .                  # numpy and periodictable only
pip install -e ".[plot,sasdata]"  # plus plotly (plots) and sasdata (Data1D output)
```

The required dependencies are `numpy` and `periodictable` (scattering lengths
and densities of atoms). Two extras are optional: `plot` installs `plotly`,
which `plot_points` and `plot_profile` need, and `sasdata` installs `sasdata`,
which `to_data1d` needs. Without them those calls raise an `ImportError` that
names the extra.

## Usage

```python
from sinccity import GenericScatteringCalculator

gsc = GenericScatteringCalculator()
gsc.load_structure('particle.sld')  # .pdb, .gro, .sld or .vtk
gsc.plot_points()  # SasView's "Draw Points", in 3D
gsc.set_q_grid(3e-4, 0.3, 30)  # optional: SasView's defaults
gsc.set_params(sld_solvent=6.34)  # optional: 10^-6 Å^-2, D2O here
profile = gsc.calculate()
profile.describe()
gsc.plot_profile()  # or gsc.plot_profile(data=some_data1d)
data = profile.to_data1d()  # a sasdata Data1D, ready for a fit
```

An aggregate from a molecular-dynamics run, in several contrasts:

```python
gsc.load_structure(
    'CTA178-eq5.gro',
    residues=['CTA'],  # the micelle, not its solvent
    masses={'CTA': cta_masses},  # from the topology, in atom order
)
gsc.set_params(scale=0.01, size_pd=0.1)
profiles = gsc.calculate_series([
    {'label': 'h-CTA in D2O', 'params': {'sld_solvent': 6.36}},
    {'label': 'h-CTA in H2O', 'params': {'sld_solvent': -0.56}},
    {'label': 'd-CTA in H2O', 'params': {'sld_solvent': -0.56, 'deuterated': ['CTA']}},
])
```

The calculation is the Debye sum over all points:

```text
w_j  = (sld_j - sld_solvent * 1e-6) * volume_j        SLD in Å^-2, volume in Å^3
D(q) = Σ_j Σ_k w_j w_k sin(q r_jk) / (q r_jk)
I(q) = scale * 1e8 * D(q) / V_norm + background       1/cm
```

`sld_solvent` is in 10⁻⁶ Å⁻², like every sasmodels SLD parameter. SasView's
calculator takes Å⁻² instead; a value that looks like Å⁻² triggers a warning.
The structure's own SLDs stay in Å⁻², as the files store them.

**Normalisation.** `V_norm` is the volume of the points that scatter, those
whose SLD differs from the solvent's. For a particle of pixels or cells in a
box of solvent, the curve at `scale=1` then equals sasmodels' form factor at
`scale=1`, so `scale` is a volume fraction. When the solvent exactly matches a
population of points, those points stop counting and the curve jumps by the
volume ratio; for contrast-variation series set
`set_params(normalisation_volume=...)` explicitly (`'auto'` restores the
default). SasView calls the override `total_volume` and, when some points
have zero contrast, also multiplies by the ratio of all points' volume to the
scattering volume. This calculator does not.

**Size polydispersity.** `set_params(size_pd=0.1)` sums blown-up and shrunk
copies of the whole structure, as sasmodels does for `radius_pd`: `size_pd_n`
size factors (35) evenly over ±`size_pd_nsigma` (3) standard deviations,
factors at or below 0 dropped, each weighted by the Gaussian's height there.
A factor f stretches every distance by f and every amplitude by f³, so

```text
I(q)   = scale * 1e8 * Σ n_f f⁶ D(f q) / V_norm + background
V_norm = Σ n_f f³ V                     V: the volume of the points that scatter
```

with the shares n_f adding up to 1, so `scale` stays a volume fraction. The
lattice, histogram and average kernels evaluate their one distance spectrum
at f·q for every size, which costs only the sums over q; `'pairs'` repeats
the whole pair sum per size and warns. The automatic histogram bin and the
cell-size warning use the largest factor, `1 + size_pd_nsigma * size_pd`.

**Deuteration.** `set_params(deuterated=['CTA'])` gives the hydrogens of
these residue names the SLD of D at `calculate()`; volumes are unchanged,
since D and H have the same natural atomic volume. It needs a structure with
residue names (PDB, GRO or `Structure.from_atoms` with `resnames=`); a name
that is not in the structure, or whose atoms include no hydrogen, is an error.
`structure.deuterated(['CTA'])` returns the same change as a new structure.

## File formats

Positions are in Å and SLDs in Å⁻² in all four (a `.gro` file's nm are converted).

| Format | What is read |
|---|---|
| PDB | `ATOM` records of the first model. `include_hetatm=True` adds `HETATM` records, never waters (`HOH`, `DOD`, `WAT` and the MD names `SOL`, `TIP3`, `TIP4`, `SPC`, `SPCE`, `OPC`). One alternate-location ID for the whole structure, occupancy not applied. Residue names from columns 18–21, where MD programs write four-letter names. Elements as below. Each element gets SasView's SLD and *natural atomic volume* from periodictable; an unknown element is an error. `CONECT` bonds are drawn. A `CRYST1` box is stored unless it is unreadable or a placeholder (a side of 1 Å or less). Residue numbers that are not decimal integers (hybrid-36) are renumbered in order. |
| GRO | All atoms, cut at GROMACS's fixed columns (names may touch numbers; `-ndec` files are read). Velocities are ignored; the box comes from the last line and is used to unwrap the aggregate. Elements as below. |
| SLD | One header line, then `x y z sld`, `x y z sld mx my mz` or `x y z sld mx my mz volume`. Magnetic columns are ignored; six-column magnetic-only files are refused. Without a volume column the points must form a lattice and each pixel's volume is the product of the steps. A sparse subset that skips planes looks like a coarser lattice, so give such files the volume column. |
| VTK | Legacy ASCII `UNSTRUCTURED_GRID` (classic or 5.1 layout) with tetrahedra, voxels, hexahedra, wedges and pyramids. The nuclear SLD is the one single-component `SCALARS` field, per cell or per point (averaged to cells as SasView does). Each cell becomes one point at its volume centroid. Binary files are refused with an export recipe. |

**Selecting residues.** `load_structure(path, residues=['CTA'])` reads only
the atoms of these residue names, from PDB and GRO files. When given, it is
the whole selection: the PDB water and `HETATM` rules do not apply. A `.gro`
file of a solvated system needs it; without it the reader notes how many
residue names it read. Options a format does not have are ignored with a
note. The structure keeps the options it was read with
(`structure.reader_options`), so `get_settings()` records them however the
structure was installed.

**Elements.** For PDB and GRO files each residue name gets its elements from
the first source that has them:

1. `masses={'CTA': [14.007, 12.011, ...]}`: the atom masses in u in topology
   order. Each mass gives the element of nearest mass, H to U and D, within
   0.01 u. Hydrogen mass repartitioning (H at 3.024 u) or united atoms (CH₂
   at 14.027 u) are no element's and an error: pass `elements=` for them.
2. `elements={'CTA': ['N', 'C', ...]}`: the element symbols in order.
3. The PDB element column, when it is filled for every atom of the residue.
4. A guess from the atom names, with a note. PDB files use SasView's column
   rule (` CA ` is carbon, `CA` starting in column 13 calcium). GRO files have
   no column convention: leading digits are dropped, then the first letter is
   the element, except that a one-atom residue whose name is a two-letter
   element (`NA`, `CL`) is that element. So `CA` in a protein is carbon and a
   `CA` ion calcium. The note lists every two-letter guess, since those are
   the likeliest to be wrong.

Masses and element lists are matched by the atom's position in its residue
(a run of atoms with one residue number and name), so the atom names of the structure
and the topology may differ. `structure.element_sources` records the source
per residue name. `structure.with_symbols(symbols)` returns a copy with
corrected elements, and `Structure.from_atoms(xyz, symbols, resids=,
resnames=, names=)` builds atoms that were read elsewhere.

**Periodic boxes.** GROMACS writes every position back into the box, which
can split molecules and the aggregate. With `unwrap='auto'` (the default) a
`.gro` file is unwrapped: each molecule is made whole around its first atom,
then moved by whole box lengths next to the aggregate's centre, the circular
mean of the molecule centres. A note counts the atoms moved. It works for one
aggregate smaller than half the box. A `.gro` file whose box is zero (no
periodic box) is read as it is, with a note. PDB files are unwrapped only with
`unwrap=True`, because many programs write a placeholder `CRYST1`; a missing
or placeholder box is then an error. Triclinic boxes cannot be unwrapped
(`unwrap=False` reads such files without the box).

The readers follow SasView's (`sas.sascalc.calculator.sas_gen`, BSD-3-Clause,
Copyright (c) 2009-2026, SasView Developers) with one deviation: the VTK cell
geometry is recomputed, because SasView's volumes for hexahedra, wedges and
pyramids are wrong. The GRO reader and the residue options come from
coords2sans.

## Things to know

- **PDB structures in a solvent are only qualitative.** The natural atomic
  volumes use bulk elemental densities, liquid H₂, N₂ and O₂ for H, N and O.
  A typical residue gets about 310 Å³ instead of about 133 Å³, so in D₂O a
  protein comes out about 8× too intense at q → 0. There is no hydration
  shell and no H/D exchange either. `calculate()` warns with the actual
  volume ratio. In vacuum (`sld_solvent=0`) the result is right. For atoms,
  `scale` is not a volume fraction.
- **Points stand in for cells.** Pixels and VTK cells are treated as points,
  which degrades as q × cell size approaches 1. `calculate()` warns above
  that, with the cell size taken as the lattice step or the cell diameter
  divided by √3.
- **Kernels.** `set_params(kernel=...)` picks how D(q) is summed; the
  profile's `method` records it.

  | `kernel` | How | Exact? | Time |
  |---|---|---|---|
  | `'lattice'` | FFT autocorrelation of the pixel grid | yes | well under a second for 10⁵ pixels |
  | `'histogram'` | pairs counted by distance, then one sum per q over the occupied distances | to about 1e-4 | about 2 s for 11 000 atoms, N² |
  | `'average'` | ⟨A⟩²: each point's distance from the bounding-box centre | only for spherically symmetric structures | linear in N |
  | `'pairs'` | every pair at every q | to the last digit | 2–3 s for 4 000 points, 25 s for 8 000 |

  `'auto'` (the default) is `'lattice'` for pixels on a lattice and
  `'histogram'` for everything else. The histogram bin is `'auto'`:
  `min(0.01 Å, 1e-3 / (q_max · f_max))`, with f_max the largest size factor
  (1 without polydispersity), which keeps the rounding phase error below
  5e-4 rad; `histogram_bin=` sets it, and `profile.histogram_bin` records it.
  A cached spectrum with a finer bin is reused.
  `'average'` reproduces SasView's `is_avg=True` and suits round aggregates
  above about 10⁵ atoms. `calculate(progress=...)` can report progress or
  stop the pair passes of `'histogram'` and `'pairs'`.
- **Contrast series are cheap.** The lattice, histogram and average kernels
  keep the scattering lengths and the volumes apart in their distance
  spectrum, so the calculator keeps the last one and a new `sld_solvent`,
  `scale`, `background` or `normalisation_volume` costs only the sums over
  q. A new structure, deuteration, kernel, a finer bin or a change of which
  points match the solvent repeats the pair pass.
- **The curve is theory.** `profile.to_data1d()` returns a `Data1D` with no
  dI and no dQ, with `qmin`, `qmax` and an empty `mask` set. A fit of it
  should state its weighting and resolution explicitly; pass
  `to_data1d(dI=...)` to attach uncertainties.
- **Reproducibility.** `gsc.get_settings()` records the file path and
  SHA-256, the reader options, the q grid, every parameter and the package
  versions. `GenericScatteringCalculator.from_settings(settings)` rebuilds
  the calculator and warns if the file or the periodictable version changed;
  records of version 0.1 replay with the new parameters at their defaults,
  except the kernel: `'auto'` for pixels on a lattice and otherwise
  `'pairs'`, the pair sum 0.1 used, so a record gives the curve it was
  made with.
  A structure no file and options can rebuild (`from_arrays`, `from_atoms`,
  or a `deuterated` / `with_symbols` copy) is recorded as `'in-memory'`, which
  `from_settings` refuses; a derived copy also records `derived_from`.
- **Series and worker processes.** `gsc.calculate_series(series)` computes
  one profile per entry, each a `label` and `set_params` keys that apply to
  that entry only. All entries are checked before anything runs, entries
  that share a deuteration state and kernel share one pair pass, and the
  profiles come back in the order given, each with its `label`. The
  parameters are left as they were. A settings record plus a series is plain
  JSON, so another process can run it:
  `GenericScatteringCalculator.from_settings(settings).calculate_series(series)`.
  Call `sinccity.set_verbosity('quiet')` there; the logger is per process.
- **Console output.** Progress lines go through the `sinccity` logger;
  `sinccity.set_verbosity('quiet')` silences them.

## Coming from coords2sans

`sinccity` takes over coords2sans's features; its names and units differ in
places.

| coords2sans | sinccity |
|---|---|
| `Structure.read(path, residues, masses, elements)` | `sinccity.load(path, residues=, masses=, elements=)` or `gsc.load_structure(...)` |
| `Structure(resids, resnames, names, positions, elements)` | `Structure.from_atoms(xyz, symbols, resids=, resnames=, names=)` |
| `structure.positions`, `structure.elements` (editable) | `structure.xyz`, `structure.symbols`; `structure.with_symbols(new)` |
| `structure.sld(deuterated=[...])` | `structure.deuterated([...]).sld`, or `set_params(deuterated=[...])` |
| `SansCalculator(solvent_sld=6.36e-6)` | `set_params(sld_solvent=6.36)`, in 10⁻⁶ Å⁻² |
| `q_min`, `q_max`, `q_points`, `q_log` | `set_q_grid(qmin, qmax, npoints, spacing='log' or 'linear')`; the default grid is SasView's, 3e-4 to 0.3 Å⁻¹ in 30 points |
| `polydispersity`, `size_points`, `SIZE_RANGE` | `size_pd`, `size_pd_n`, `size_pd_nsigma` |
| `average_amplitude=True` | `kernel='average'` |
| `BIN_WIDTH` (0.01 Å) | `histogram_bin`, by default at most 0.01 Å and smaller at high q or with a size spread |
| `calculate(structure, label)` → `Curve` | `set_structure(structure)`, then `calculate(label=)` → `ScatteringProfile` |
| `Curve.save`, `Curve.plot(ax)` | `profile.save_csv(path)` (the same two columns); `ax.plot(profile.q, profile.intensity)` or `profile.plot()` (plotly) |
| `SansTask` | `get_settings()` plus a series, run by `from_settings(settings).calculate_series(series)` |

Two rules differ on purpose. PDB files without an element column use
SasView's name rule (` CA ` is carbon), and in `.gro` files a two-letter
element name counts only in a one-atom residue (so `CA` in a protein is
carbon). V_norm is the volume of the atoms that scatter, which equals
coords2sans's volume of all atoms unless an atom matches the solvent.

## Tests

```bash
pytest tests
```

## Licence

BSD 3-Clause, like IQuana. The test structures in `tests/data/` are
SasView test data, BSD 3-Clause, Copyright (c) 2009-2026, SasView Developers.
