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.

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.

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.

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.

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:

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:

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.

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:

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.

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:

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.