# Superspace: incommensurate crystallography in (3+d) dimensions

`pyextal.superspace` handles modulated structures — crystals whose diffraction
pattern is sharp but does not sit on a three-dimensional lattice, so three
integers cannot index it.

Every reflection is labelled by $3+d$ integers:

$$\mathbf{S}(H) = h\,\mathbf a^* + k\,\mathbf b^* + l\,\mathbf c^*
  + m_1 \mathbf q_1 + \cdots + m_d \mathbf q_d$$

Main reflections have every $m_j = 0$; the rest are satellites, of order
$\sum_j |m_j|$. `d = 1, 2, 3` are all supported.

The worked example is 1T-TaS$_2$ — see
`example/superspace/TaS2_superspace_SAED.ipynb`.

```python
from pyextal.superspace import tas2, saed_pattern, plot_pattern

lat, atoms = tas2.structure()                    # the published refinement
pattern = saed_pattern(lat, atoms, zone=(0, 0, 1),
                       condition=tas2.reflection_condition)
plot_pattern(pattern)
```

`tas2.structure()` returns the model refined by Spijkerman *et al.* — 21
harmonics on Ta and 12 on S, in superspace group $P\bar3(\alpha\beta\tfrac13)$.
`model='standin'` gives a small invented triple-q wave instead, which is faster
and useful where the real amplitudes are not the point.

## Modules

| Module | Role |
| :--- | :--- |
| `lattice` | `SuperspaceLattice`: the metric, the q-vectors, index arithmetic, rank and commensurability |
| `indexing` | Enumerating and labelling reflections; arrays are `(n, 3+d)` throughout |
| `modulation` | Displacive and occupational modulation as Fourier series in internal space; `ModulatedAtom` |
| `structurefactor` | $F(H)$ by quadrature over internal space, plus the Bessel closed form |
| `symmetry` | Superspace operations, the $\varepsilon$ solver, reflection transformation |
| `scattering` | Form factors read from the engine's own `sftabl.sct` |
| `saed` | Kinematic zone-axis patterns with relrod weighting |
| `approximant` | Commensurate supercell → `.dat` for the unmodified Bloch engine |
| `engine` | Drives the Fortran engine with a modulated structure (Stage B) |
| `tas2` | The 1T-TaS$_2$ example structure |

## Rank: how many dimensions the structure actually needs

`d` q-vectors give a $(3{+}d)$D description **only if they are independent
modulo the main lattice**. If some integer combination $\sum_j n_j \mathbf q_j$
lands on a reciprocal lattice point, satellites coincide and a refinement in
that many dimensions fits parameters the data cannot separate.

1T-TaS$_2$ is the case in point. Its charge-density wave is triple-q, three arms
related by the three-fold axis, which suggests $(3{+}3)$D. But with the refined
wave vectors the third arm is **exactly** $\mathbf q_2 - \mathbf q_1$ — no
lattice vector needed to make it work — so the three span a rank-2 module and
the correct description is $(3{+}2)$D, which is what the published refinement
used. The paper's own satellite grouping says the same: it lists $(1,0)$,
$(0,1)$ and $(-1,1)$ as one symmetry-equivalent set.

```python
q1, q2 = tas2.published_q('tripled')
lat = SuperspaceLattice(tas2.CELL_TRIPLED, np.vstack([q1, q2, q2 - q1]))
# UserWarning: ... not independent modulo the main lattice: q1 - q2 + q3 = 0
lat.rank()                 # 2, not 3
lat.integer_relations()    # [(array([1, -1, 1]), array([0, 0, 0]))]
```

Commensurate (rational) q-vectors are rank-deficient too, but that is a
legitimate thing to build, so it gets a different message pointing at the
supercell rather than reporting a mistake.

## Satellite order is not `sum |m|`

The obvious measure of how far out a satellite is, $\sum_j |m_j|$, is **not
invariant under the crystal's symmetry** once $d > 1$, and using it as a cutoff
silently destroys the symmetry of everything downstream.

In 1T-TaS$_2$ the three-fold sends $(1,0) \to (0,1) \to (-1,1)$. The first two
have $\sum|m| = 1$ and the third has 2, so a cutoff at 1 keeps two thirds of a
symmetry orbit and the computed pattern comes out visibly **not** three-fold.
The reference makes the same point in section II: *"'order of a satellite'
cannot be represented by a single number anymore"*.

What is invariant is the length of the modulation part,
$|\sum_j m_j \mathbf q_j|$, since any symmetry operation permutes the
$\mathbf q$-vectors. `enumerate_reflections` cuts on that, in units of the
shortest $\mathbf q$, and `ReflectionSet.order` ranks it, so
symmetry-equivalent satellites always share an order. For $d = 1$ the two
measures coincide. `ReflectionSet.index_sum` still gives $\sum|m|$ for
labelling — just never filter on it.

## The structure factor

Each atom is smeared over the internal coordinate, so the sum over atoms becomes
a sum of integrals:

$$F(H) = \sum_\mu f_\mu(s)\,T_\mu(s) \int_{[0,1]^d} \mathrm d\bar x\;
  P_\mu(\bar x)\;
  e^{2\pi i [\, \mathbf h\cdot\mathbf r_\mu + \mathbf m\cdot\bar x
  + \mathbf S(H)\cdot \mathbf u_\mu(\bar x)\,]}$$

`structure_factor` evaluates this by quadrature. For a single harmonic the
integral is $J_{-m}(2\pi\,\mathbf S\!\cdot\!\mathbf A)$, and the quadrature
agrees with it to ~1e-22 relative — the rectangle rule converges exponentially
on a smooth periodic integrand, so nothing is given up by taking the general
route that also handles many harmonics and occupational modulation.

The one thing that fails silently is an under-resolved internal grid, whose
required size grows with satellite order and modulation amplitude.
`auto_n_internal` sizes it, and `structure_factor` uses that by default.

Form factors come from `sftabl.sct`, the same Cromer-Mann table the Fortran
engine reads, with the engine's own Mott-Bethe conversion — so the kinematic and
dynamical halves use identical scattering.

## Dynamical calculations: the commensurate approximant

When the q-vectors are rational the structure is not incommensurate at all — it
is an ordinary periodic crystal with a bigger cell. That supercell is plain 3D,
so **the existing Bloch-wave engine computes it with no modification**.

```python
lat, atoms = tas2.structure(commensurate=True)     # published, approximated
res = dynamical_pattern(lat, atoms, thickness=200.0, gmax=0.7)
res['H']          # superspace indices of every output beam
res['intensity']  # dynamical intensities
```

`supercell_matrix()` derives the cell from the q-vectors alone — for
1T-TaS$_2$ it returns the $\sqrt{13}a \times \sqrt{13}a \times 3c$ "Star of
David" cell, 117 atoms. `assign_superspace_index` maps the engine's supercell
reflections back to $(3{+}d)$ labels, so the dynamical result lines up with the
kinematic one reflection by reflection.

The supercell reciprocal lattice is $|\det M|$ times denser than the parent's,
so keep `gmax` modest: a value that is unremarkable for the basic cell can
select thousands of beams.

## The 1T-TaS2 model

`tas2` carries Tables III–V of the reference verbatim. Two settings are available:

* `setting='tripled'` (default) — $a_0 \times a_0 \times 3c_0$, the cell the
  refinement uses. The wave vectors are purely in-plane and the three layers
  come from superspace centring translations that carry an internal phase
  shift. Nine atoms.
* `setting='basic'` — $a_0 \times a_0 \times c_0$, with the wave vectors picking
  up components of $\pm\tfrac13$ along $\mathbf c_0^*$. `published_q`,
  `tripled_to_basic` and `basic_to_tripled` support it, and the *indexing* is
  exact: both settings place every reflection identically to 2e-16, and give the
  same structure factors to 8e-16 for the **unmodulated** structure.

  **The modulated structure is only available in the tripled setting.** In the
  basic cell $\mathbf q$ has a $\mathbf c^*$ component, so an atom's own $z$
  enters its internal coordinate — $\bar x_{basic} = \bar x_{tripled} + q_z z$ —
  while the published modulation functions are written against the tripled
  convention, where it does not. The compensating phase is atom-dependent;
  neither a constant internal shift nor a per-atom shift by $q_z z$ reproduces
  the tripled model (both leave 2–4% residuals). `published_atoms` raises rather
  than return a structure that is quietly wrong.

The centring imposes $l + 2m_1 + m_2 \equiv 0 \pmod 3$
(`tas2.reflection_condition`), which has a consequence worth knowing: the [001]
zero-order Laue zone is $l = 0$, so it needs $2m_1 + m_2 \equiv 0$, and **no
first-order satellite satisfies that**. A [001] pattern shows the main
reflections and the second-order $(1,1)$-type satellites; the first-order ones
live in $l_0 = n \pm \tfrac13$.

Nothing filters the forbidden reflections — the three centred layers cancel them
to fourteen orders of magnitude, which is the check that the internal phase
shifts are right.

### What was verified

Tables IV and V quote amplitudes "in angstrom" for the tripled cell, which is
ambiguous, so the convention was checked against the paper's own figures rather
than assumed:

| quantity | recomputed | paper |
| :--- | ---: | ---: |
| Ta max $\|u_{xy}\|$ (Fig. 6) | 0.3080 Å | 0.308 Å |
| Ta–Ta(100) shortest (Figs. 3–5) | 3.0865 Å | 3.086 Å |
| Ta–Ta(100) longest (Figs. 3–5) | 3.8079 Å | 3.806 Å |
| S $u_z$ max (Fig. 10) | +0.116 Å | +0.106 Å |
| S $u_z$ min (Fig. 10) | −0.114 Å | −0.124 Å |

The three Ta numbers agree to four significant figures, fixing the convention as
displacements in Å along $a$, $b$, $c$ of the tripled cell. The S extremes are
about 0.01 Å out; that is not a transcription error (flipping any harmonic
group's sign makes it worse, and every Table I symmetry relation closes), and
Fig. 10 plots a particular symmetry copy read off contours spaced 0.0177 Å
apart.

**One value is corrected against the source.** Table V's $B^y$ for the harmonic
$(-1,1)$ is entered as $-0.0107$, not $+0.0107$: Table I's own relation
$B^y_{-m_1-m_2,m_1} = -B^x_{-m_2,m_1+m_2}$ requires the minus, and with the plus
the sulfur modulation breaks its three-fold site symmetry by 0.043 Å — a third
of its amplitude. The PDF renders minus as the digit 2 and this one is lost on
extraction. Every dependent amplitude in both tables is checked against Table I
in the test suite.

Table III places S at $(1/3, 1/3, 0.0862)$ while Fig. 1 gives
$\pm(1/3, 2/3, z)$ and §V C refers to $(2/3, 1/3, 0.086)$. The standard CdI2
positions $\pm(1/3, 2/3, z)$ are used, as the only choice giving Ta octahedral
sulfur coordination — corroborated by the refinement reporting every Ta cosine
amplitude as zero, which requires Ta on the inversion centre.

## Modulated structures inside the Bloch engine

The engine carries superspace indices and a modulated structure factor of its
own:

* `blochwave.superspace` holds `dmod`, the wave vectors `qvec`, the satellite
  indices `hklm` (parallel to `difpar.hkl`), and the per-atom harmonics.
* `difpar.gvec` holds each beam's **effective** index $\mathbf{hkl} + \sum_j m_j
  \mathbf q_j$, which is not integral and so cannot live in `difpar.hkl`. Every
  geometry calculation reads `gvec`, so the ordinary and modulated cases follow
  one path.
* `sfact_mod` integrates the modulated structure factor over internal space, and
  `setugh` calls it instead of `sfact` whenever `dmod > 0`.

```python
from pyextal.superspace import engine, enumerate_reflections, plot_image, tas2

lat, atoms = tas2.structure(setting='tripled')
beams = enumerate_reflections(lat, gmax=0.55, max_index=3, max_order=2,
                              condition=tas2.reflection_condition)

# the whole sequence: write the cell, parse, load the modulation, solve
pattern = engine.dynamical_saed(lat, atoms, beams.H, thickness=200.0,
                                zone=(0, 0, 1), hv=200.0)
plot_image(pattern, sigma=3.0, gamma=1.5)
```

or step by step, to keep the engine loaded for something else:

```python
from pyextal.callBloch import bloch_parse

engine.write_dat('mod.dat', lat, atoms, beams.H, aper=2.0)
bloch_parse('mod.dat')
engine.load(lat, atoms, beams.H)       # U(g,h) is now modulated
engine.unload()                        # and back to the ordinary matrix
```

The output aperture is applied **during parsing**, before the modulation
exists, so it can only see the main index. `dynamical_saed` therefore keeps
every beam and cuts on the true effective $|g|$ afterwards; if you drive the
engine yourself, make `aper` generous and filter in Python.

### What is verified

| check | result |
| :--- | ---: |
| `dmod = 0` against the pre-change binary (`hkl`, U(g,h), LARBED) | bit-identical |
| `sfact_mod` with zero modulation vs `sfact` | 9e-8 |
| `gvec` vs `hkl + m.q` | 9e-8 |
| Fortran `sfact_mod` vs the Python `structure_factor`, 320 reflections | 2.3e-5 |

The first is the gate: `gvec` is seeded as `real(hkl)`, which is exact, so an
ordinary crystal is untouched. The last is the real validation — the Fortran
integrates over internal space in its own loop and the Python does it in numpy,
sharing only the Cromer-Mann table.

### Limits

* **P1 only.** Relating the modulation functions of symmetry-equivalent atoms
  needs a phase that depends on each atom's own position; `load` refuses
  anything else rather than putting the wrong modulation on generated atoms.
* **Explicit beam lists.** `sele` enumerates integer indices and cannot reach a
  satellite, so write the `hkl` block out. `engine.write_dat` does.
* **The `.dat` grammar is untouched.** The modulation is pushed in from Python
  after parsing rather than parsed from the file, which keeps the
  order-sensitive format out of the blast radius.
* **One crystal per engine** — the engine keeps the parsed crystal and the loaded
  modulation in module-global state, and `terminate()` ends it for the life of
  the process. `engine.run_isolated(func, ...)` puts each calculation in a fresh
  interpreter, which is how a notebook compares two crystals or two routes.

## The commensurate phase

Below 180 K the CDW locks in and the in-plane cell becomes
$\sqrt{13}a \times \sqrt{13}a$ — 13 basic cells, the Star of David.
`tas2.commensurate_structure()` builds that model, **in the paper's own
$(3\!+\!2)$ basis**: $\mathbf q_1 = (\sigma_1, \sigma_2, 0)$ and
$\mathbf q_2 = (-\sigma_2, \sigma_1{+}\sigma_2, 0)$ with
$(\sigma_1, \sigma_2) = (3/13,\, 1/13)$, which is `tas2.SIGMA_C`.

The basis is load-bearing. The harmonic labels $(n_1, n_2)$ of Tables IV and V
are written against that pair, so expressing the same lattice through the
three-fold triplet re-labels every published amplitude. It is silent and it is
large: the site residual goes from $4.7\times10^{-4}$ Å to $0.57$ Å, and the
diffraction pattern stops being three-fold.

The full C phase also stacks with a $13c$ repeat. That is **not** represented,
and cannot be in this basis: with arms $\{\mathbf q_1, \mathbf q_2,
\mathbf q_2 - \mathbf q_1\}$ the $c$ components are $q_z$, $q_z$ and $0$, so a
nonzero $q_z$ makes the third arm inequivalent to the other two and
`epsilon_for` reports the three-fold as not a symmetry. Representing the
stacking needs a 13-fold $c$ axis carrying the layer shifts as centring
translations, exactly as the paper puts the thirds in a tripled cell — and the
C-phase stacking sequence that would fix those translations is not published.
`commensurate_structure(cz=...)` accepts the parameter for exploring it and
warns that the result is not the C phase.

### Coincident reflections

Once $\mathbf q$ is rational, the Fourier module collapses onto a 3D lattice and
the map from a superspace index to a reciprocal-lattice point stops being
injective. If $\mathbf c \cdot \mathbf q$ is an integer vector $d\mathbf h$,
then $(\mathbf h, \mathbf m)$ and $(\mathbf h - d\mathbf h, \mathbf m +
\mathbf c)$ are the **same physical reflection**, and the crystal has one
Fourier coefficient there — the *coherent sum* of their amplitudes.

Two places need this, and they need it in different forms:

```python
from pyextal.superspace import merge_coincident, coincidence_coset

rep, group, ngroups = merge_coincident(lat, H)   # 437 indices -> 305 reflections
c, dh = coincidence_coset(lat)                   # 37 members within |dh| < 1.05
```

- **The reflection list.** `saed_pattern` calls `merge_coincident` by itself
  whenever `lat.is_commensurate()` (override with `merge=`). Without it the one
  intensity is split between the coincident indices in whatever way the
  enumeration happened to order them, so symmetry-equivalent spots get different
  shares: the worst three-fold mismatch is $0.63$ unmerged against
  $2.5\times10^{-3}$ merged.
- **Every $U(\mathbf g, \mathbf h)$ of the dynamical matrix.** Here merging the
  list is not enough — the index *differences* $\mathbf H_i - \mathbf H_j$ reach
  coincidences that are not in the beam list at all — so the sum happens inside
  the Fortran structure factor. `engine.load` switches it on whenever the
  lattice is commensurate (override with `coset=`, size it with `coset_dh=`).

Without the coset the superspace engine computes the **$t$-averaged** structure
rather than the actual crystal: integrating over a continuous internal
coordinate averages over the modulation phase, which for an incommensurate
structure is right and for a commensurate one is not. Its satellites come out
41% wrong.

**Truncate on the shift, never on the satellite order.** The coset is a lattice
and so infinite. $|\mathbf c \cdot \mathbf q|$ is a reciprocal-space length and
therefore invariant under the point group, so cutting there keeps whole shells
and converges monotonically. $\sum_j |c_j|$ is not invariant: cutting there
admits part of a symmetry orbit and drops the rest, which *adds* asymmetry.
Measured against the exact supercell:

| cut | members | worst relative difference |
| :--- | ---: | ---: |
| none | 1 | $4.1\times10^{-1}$ |
| $\lvert d\mathbf h\rvert < 0.40$ | 7 | $5.0\times10^{-2}$ |
| $\lvert d\mathbf h\rvert < 0.70$ | 19 | $9.7\times10^{-3}$ |
| $\lvert d\mathbf h\rvert < 1.05$ (default) | 37 | $1.2\times10^{-3}$ |
| order $\sum_j\lvert c_j\rvert \le 4$ | 3 | $1.7$ |

The coset multiplies the cost of every structure-factor evaluation by its
member count, so it is the expensive part of a commensurate dynamical run.
Lower `coset_dh` to trade accuracy for speed; the table above is the exchange
rate.

### Reduce the beam list first

A dynamical calculation must be given *distinct* reflections, and on a
commensurate lattice an enumerated list is not. At `gmax=0.40, max_index=3,
max_order=3` the C phase of 1T-TaS2 enumerates 581 superspace indices that are
only **209 distinct reflections**.

```python
from pyextal.superspace import distinct_beams

H = distinct_beams(lat, beams.H)      # 581 -> 209, and warns that it did
```

`dynamical_saed` and `supercell_saed` call this themselves, so **results come
back indexed by the reduced list** — `sup['H']` and `ss.H` — which is shorter
than what was passed in. Join on those. `engine.load` raises instead of
reducing, because by the time it runs the `.dat` has already been written.

Skipping the reduction does not produce a wrong number, it produces a crash a
long way from its cause: duplicate beams are identical rows in the structure
matrix, and the `.dat` written from such a list repeats indices in its `hkl`
block, so `bloch_parse` cannot resolve some `out` entries and appends them past
the array allocated from the declared `nbeams`. Expect a failed eigensolve, a
segmentation fault, or the child process killed outright.

### Two routes to the same answer

The commensurate structure can be computed two independent ways, and the
agreement is the strongest check on the superspace engine:

```python
from pyextal.superspace import engine, enumerate_reflections, supercell_saed, tas2

lat, atoms = tas2.commensurate_structure()
beams = enumerate_reflections(lat, gmax=0.40, max_index=2, max_order=1)

sup = engine.run_isolated(supercell_saed, lat, atoms, beams.H, thickness=200.0)
ss = engine.run_isolated(engine.dynamical_saed, lat, atoms, beams.H, thickness=200.0)
```

`supercell_saed` freezes the modulation into 39 atomic coordinates and hands an
ordinary crystal to the unmodified engine — it never mentions a satellite.
`dynamical_saed` keeps 3 atoms, integrates the structure factor over internal
space and sums the coincident amplitudes. They agree to **1.2e-3**, the
resolution of that quadrature. Both must be given the same beam list, or they
are solving different scattering problems.

Note also that the first-order satellites are visible in a $[001]$ pattern of
the commensurate phase and are not in the nearly commensurate one: $q_z = 1/3$
lifts them out of the zero-order Laue zone, and locking in puts $\mathbf q$ in
plane and brings them back.

## Diffraction patterns as images

`render_image` puts intensities on a pixel grid, giving a diffraction *image*
rather than a scatter of markers — something to threshold, line-scan or compare
against a plate. It convolves with the pseudo-Voigt the JAX refinement fits, in
the same functional form, so a refined `sigma` and `gamma` transfer directly.

```python
from pyextal.superspace import render_image, plot_image, saturation_level

image, extent = render_image(pattern, shape=(700, 700), sigma=3.0, gamma=1.5)
plot_image(pattern, sigma=3.0, gamma=1.5, colorbar=True)
```

`vmax` sets the intensity that maps to full scale, in the units of the returned
image; smaller brightens, larger dims, and the returned array is untouched
either way. Its main use is putting several patterns on **one** scale — left to
itself `plot_image` normalises each pattern to its own strongest diffracted
beam, so two panels look equally bright whatever the physics did:

```python
render = dict(shape=(650, 650), sigma=3.0, gamma=1.5)
level = saturation_level(kinematic, **render)          # the automatic choice
plot_image(kinematic, vmax=level, **render)
plot_image(dynamical, vmax=level, **render)            # now comparable
```

`saturation_level` computes that default without drawing anything. Give it the
same render arguments: the convolution spreads a fixed number of counts, so the
peak height depends on the grid and the kernel.

Two details are not cosmetic. The direct beam is some seven orders of magnitude
above the strongest diffracted beam, so `plot_image` references the display to
the *diffracted* maximum and lets the centre clip, as a plate does. And the
kernel is tapered to a disc: the Lorentzian decays only as $1/r^2$, so a square
truncation leaves a visible square halo around every strong reflection.

## Limits

- **Modulated ADPs are not implemented.** `B` is constant per atom.
- **Symmetry expansion of modulation functions** works only for operations that
  map the q-set with no reciprocal-lattice offset ($G_j = 0$). A nonzero $G_j$
  makes the phase relation between equivalent modulation functions depend on the
  atom's position; `expand_atoms` raises rather than guessing. List such
  structures in $P1$. `transform_reflection` is general and has no such limit.
- **Genuinely incommensurate structures have no dynamical path here.** They have
  no finite supercell. Reaching them means teaching the Fortran engine
  $(3{+}d)$ indices, which `difpar.hkl` — `integer(3, nbeams)` — cannot hold,
  since satellites are not at integral positions.
- **A dynamical calculation needs a commensurate approximant.** The refined NC
  structure is incommensurate, so `tas2.structure(commensurate=True)` evaluates
  the same refined modulation functions on
  $\mathbf q_c = (3\mathbf a_0^* + \mathbf b_0^*)/13$, giving the
  $\sqrt{13}a \times \sqrt{13}a \times 3c$ cell with 117 atoms.
- **Higher-order satellites carry the refinement's own limits.** Table II of the
  reference reports $R_F$ rising from 0.038 on main reflections to 0.15 at order
  $(2,0)$ and past 0.4 for the sparse high-order groups.

## Reference

A. Spijkerman, J. L. de Boer, A. Meetsma, G. A. Wiegers and S. van Smaalen,
"X-ray crystal-structure refinement of the nearly commensurate phase of
1T-TaS₂ in (3+2)-dimensional superspace", *Phys. Rev. B* **56**, 13757 (1997),
[doi:10.1103/PhysRevB.56.13757](https://doi.org/10.1103/PhysRevB.56.13757).