Diffraction¶

An X-ray or neutron experiment never sees atomic positions. It shines a beam at a sample and records how much radiation scatters into each angle. The quantity that comes out is the static structure factor \(S(k)\), and if you want to compare a simulation against an experiment, this is the language you have to speak.

\(S(k)\) carries the same information as \(g(r)\). It is that information expressed in reciprocal space, where lengths become inverse lengths and "atoms 3.7 Å apart" becomes "a density wave with wavenumber \(k\)".

Where the formula comes from¶

Two atoms separated by \(\mathbf{r}_{ij}\) scatter waves that arrive at the detector with a phase difference \(\mathbf{k}\cdot\mathbf{r}_{ij}\). A liquid or powder has no preferred orientation, so every orientation of \(\mathbf{r}_{ij}\) occurs with equal probability. Averaging \(e^{i\mathbf{k}\cdot\mathbf{r}}\) over all directions of \(\mathbf{r}\) gives \(\sin(kr)/(kr)\), and summing over all pairs gives the Debye scattering equation:

\[ S(k) = \frac{1}{N}\left\langle \sum_{i}\sum_{j} \frac{\sin(k r_{ij})}{k r_{ij}} \right\rangle . \]

MolPy evaluates this directly from coordinates — no histogram and no Fourier transform of \(g(r)\).

One property of that evaluation matters more than it first appears: the sum runs over literal pair distances, not minimum-image distances. The estimator treats your coordinates as an isolated cluster in open space, even when the frame has a periodic box — which is why its compute takes positions and no box. This matches freud's Debye implementation, and it has two consequences. Pairs are counted out to the full diagonal of the box rather than being capped at \(L/2\), and the configuration has a surface, so the outermost atoms have fewer neighbours than bulk ones. For the 500-atom argon box used below, that surface is a real finite-size effect on peak heights — it is part of why the main peak comes out lower than the experimental value quoted later. Peak positions are robust; absolute intensities from a small box are not.

The equivalent statement in terms of \(g(r)\) makes the connection explicit:

\[ S(k) = 1 + 4\pi\rho\int_{0}^{\infty} r^{2}\,[g(r)-1]\, \frac{\sin(kr)}{kr}\,\mathrm{d}r . \]

Note what is being transformed: \(g(r) - 1\), the deviation from randomness. A structureless fluid has \(g = 1\) everywhere, the integral vanishes, and \(S(k) = 1\) at every \(k\). Structure in real space is what puts features in reciprocal space.

Reading a real curve¶

Figure 1. \(S(k)\) of liquid argon from the Debye equation, averaged over 50 configurations of the 500-atom reference run (seed 0, 88.3 K). The main peak sits at 1.98 Å⁻¹ and the oscillations damp out by about 6 Å⁻¹.

The main peak at \(k = 1.98\) Å⁻¹ is the dominant density wave of the liquid — the reciprocal-space fingerprint of the first coordination shell. Its position is the number a diffraction experiment reports; neutron diffraction on liquid argon at 85 K puts it close to 2.0 Å⁻¹. Its height here, 2.07, falls short of the experimental ≈ 2.7 for the finite-size reason given above; do not read that as the force field failing.

The damped oscillations past the main peak are the second and third shells. They die out around 6 Å⁻¹, and beyond that \(S(k) \to 1\): at short enough wavelengths the liquid looks structureless.

Two cautions come with this figure, and both catch people.

Do not convert the peak with \(2\pi/k\). That gives \(2\pi/1.98 = 3.18\) Å, while the actual first-neighbour distance from \(g(r)\) is 3.68 Å. The \(2\pi/k\) rule works for Bragg planes in a crystal, not for the main peak of a liquid, because that peak is a broad superposition of many pair distances rather than one repeat spacing. To go from \(S(k)\) to distances, transform the whole curve; do not read one point off it.

The figure starts at 1 Å⁻¹ on purpose. The Debye sum above includes the \(i = j\) terms, each contributing 1, so as \(k \to 0\) every term approaches 1 and \(S(k) \to N\). That runaway is bookkeeping, not the thermodynamic limit \(S(0) = \rho k_B T \chi_T\). Reading compressibility off small-\(k\) values from this estimator is wrong, and in a finite box \(k < 2\pi/L\) has no meaning anyway (0.22 Å⁻¹ for the 28.9 Å argon box).

Computing it¶

A crystal is the honest way to test the machinery, because you can predict the answer. FCC reflections appear at \(k = 2\pi\sqrt{h^2+k^2+l^2}/a\) with \(h,k,l\) all even or all odd.

import numpy as np
import molpy as mp

a = 5.26  # argon FCC lattice constant, Å
basis = np.array([[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5]])
xyz = np.array(
    [(np.array([i, j, k]) + b) * a
     for i in range(4) for j in range(4) for k in range(4) for b in basis]
)
frame = mp.Frame()
frame["atoms"] = {"x": xyz[:, 0], "y": xyz[:, 1], "z": xyz[:, 2]}
frame.box = mp.Box(mp.core.Cell.lengths([4 * a] * 3))

Where to evaluate \(S(k)\) is configuration, so the \(k\) grid goes to the constructor. The data is one snapshot's positions, read from the frame; the Debye sum runs over every pair as given and takes no box. The result carries k_values, sk and n_particles:

from molpy.compute import StaticStructureFactorDebye

k = np.linspace(1.0, 4.5, 140)
debye = StaticStructureFactorDebye(k)
sf = debye.compute(positions=frame.coords)

print(sf.n_particles)                              # -> 256
print(round(float(sf.k_values[np.argmax(sf.sk)]), 2))   # -> 2.06

The strongest peak lands at 2.06 Å⁻¹; the (111) reflection is predicted at \(2\pi\sqrt{3}/5.26 = 2.07\) Å⁻¹. Peaks for (200), (220), and (311) follow at 2.39, 3.38, and 3.96 Å⁻¹.

They are broad rather than sharp because the crystal is only four unit cells across. Peak width in a diffraction pattern is set by the size of the ordered domain — the same physics behind the Scherrer equation used to size nanocrystals. Enlarge the lattice and the peaks narrow.

S(k) is an ensemble average like \(g(r)\), and a trajectory is a loop: the one configured instrument applied to each frame's positions, the curves averaged afterwards:

frames = [frame, frame]                   # in practice, your trajectory
s_mean = np.mean([debye.compute(positions=f.coords).sk for f in frames], axis=0)
print(s_mean.shape)                       # -> (140,)

Cost. Every \(k\) point sums over every pair: \(\mathcal{O}(N^2 n_k)\) per frame. Doubling the atom count quadruples the work, and a 1000-point \(k\) grid costs ten times a 100-point one. Choose the coarsest grid that resolves your peaks.

The periodic estimator samples the reciprocal lattice of the box¶

The Debye sum treats the configuration as an isolated cluster. A periodic box admits a second estimator that respects the periodicity: evaluate the density's Fourier component \(\rho_{\mathbf k} = \sum_j e^{-i\mathbf k\cdot\mathbf r_j}\) only at wavevectors compatible with the box, \(\mathbf k = 2\pi\,\mathbf n / L\), and take \(S(\mathbf k) = |\rho_{\mathbf k}|^2 / N\). StaticStructureFactorDirect does that, either at explicit k_vecs or binned isotropically up to k_max in n_bins shells. At a reciprocal-lattice vector of a perfect crystal every atom scatters in phase, so \(S = N\) exactly:

from molpy.compute import StaticStructureFactorDirect

g111 = 2 * np.pi / a * np.array([1.0, 1.0, 1.0])
bragg = StaticStructureFactorDirect(k_vecs=[g111]).compute(
    positions=frame.coords, simbox=frame.box
)
print(bragg.sk, np.round(bragg.k_magnitudes, 2))   # -> [256.] [2.07]

The Debye peak at the same \(k\) is much lower, because it averages \(|\rho_{\mathbf k}|^2\) over every direction of \(\mathbf k\) — most of which miss the lattice. Use the Debye form to compare with a powder or liquid experiment, and the direct form when the periodicity of the box is part of the question.

For a picture rather than a curve, DiffractionPattern(n_grid, sigma) projects the snapshot along one axis (axis=2 by default), smears it onto an n_grid×n_grid image and returns that image beside its centred 2-D power spectrum diffraction — the pattern a beam along that axis would record.

When it goes wrong¶

\(S(k)\) shoots up at small \(k\). Expected, not a bug — the self terms, as above. Start your grid near \(2\pi/L\) at the very lowest, and do not interpret the rise.

The first point equals the number of atoms. Your grid includes \(k = 0\), where every \(\sin(kr)/(kr)\) is 1 and the Debye sum returns \(N\). That value is a count, not a measurement; start above zero.

\(S(k)\) does not approach 1 at large \(k\). Either far too few frames, or the coordinates are not what you think — check for a mangled box or duplicated atoms.

The peak position disagrees with experiment but \(g(r)\) looks fine. Check units. A GROMACS trajectory is in nm; a \(k\) grid in Å⁻¹ against coordinates in nm gives a peak off by a factor of 10.

The calculation takes forever. \(\mathcal{O}(N^2 n_k)\). Reduce the \(k\) grid first, then the number of frames, then the system size.

Check yourself¶

Scatter points at random in a box and compute \(S(k)\). You should get 1 at every \(k\) except the small-\(k\) self-term rise — the reciprocal-space statement of "no structure".
Compare the main peak position of argon (1.98 Å⁻¹) with \(2\pi/r_1\) using the first \(g(r)\) peak (3.68 Å). They disagree, and now you know why.
Double the FCC lattice to \(8^3\) cells and watch the Bragg peaks narrow while their positions stay put.

References¶

P. Debye, Ann. Phys. 351, 809 (1915) — the scattering equation.
J. L. Yarnell, M. J. Katz, R. G. Wenzel, S. H. Koenig, Phys. Rev. A 7, 2130 (1973) — neutron diffraction from liquid argon at 85 K.
J.-P. Hansen, I. R. McDonald, Theory of Simple Liquids, 4^th ed. (2013) — chapter 4 for \(S(k)\), its \(k\to 0\) limit, and the link to \(g(r)\).