XRD CSV to phase insight: a minimal honest pipeline

An XRD pattern is one of the most honest measurements in a lab and one of the easiest to oversell. The numbers going in are simple — intensity versus 2θ — but between that array and "the sample is phase-pure γ-alumina" there are five decisions, and each one can turn into a silent error.

This is the smallest pipeline that produces defensible numbers, using nothing but numpy and scipy. It does not try to be a Rietveld engine. It ends where identification begins: with d-spacings and a reference comparison, clearly labelled as a hypothesis.

Step 0 — load with the metadata, not just the numbers

The community .xye format is three columns: 2θ, intensity, estimated standard deviation. numpy.loadtxt is genuinely enough — the format is its own beauty:

import numpy as np

# .xye: 2theta_deg, intensity, esd
two_theta, intensity, esd = np.loadtxt("sample.xye", unpack=True)

step = np.median(np.diff(two_theta))
print(f"{len(two_theta)} points, step ~{step:.4f} deg, "
      f"range {two_theta[0]:.2f}-{two_theta[-1]:.2f}")

The step size and the anode wavelength (Cu Kα1 = 1.5406 Å for most lab diffractometers) are not optional trivia — every later calculation depends on them. If the export dropped them, find them in the instrument log before continuing. Vendor binaries (.raw, .brml, .xrdml) carry this metadata; tools like pymatgen and xrayutilities read several of them, and /docs/instruments/xrd covers the format zoo.

Step 1 — background: choose a method and admit it

The background is air scatter, fluorescence, and any amorphous or nanocrystalline fraction. It is not a constant, and different labs subtract it differently. Two defensible options:

• SNIP — statistics-sensitive iterative clipping. Robust for powder patterns, deterministic, ~15 lines.
• Rolling minimum — the simplest defensible method: estimate the background as a moving minimum, then smooth it. Fast; adequate for flat backgrounds.

from scipy.ndimage import minimum_filter1d, uniform_filter1d

def rolling_background(y, window=101, smooth=51):
    """Moving-minimum background with a smoothing pass. window in points."""
    bg = minimum_filter1d(y, size=window)
    bg = uniform_filter1d(bg, size=smooth)
    return bg

bg = rolling_background(intensity)
corrected = intensity - bg

# Sanity check: the background should never be above a genuine peak.
assert (bg <= intensity + 1e-9).all(), "background exceeds data — widen the window"

Record the window sizes. A weak reflection that appears or disappears with the window size was never a reflection.

Step 2 — strip Kα2 (or decide not to)

Cu Kα radiation is a doublet: Kα1 at 1.5406 Å and Kα2 at 1.5444 Å, the second at half intensity. At low angles the doublet merges; at high angles it visibly splits every peak. For lattice-parameter work you want it stripped; for a quick "is there a peak here" check, you can skip it and say so.

The Rachinger method is the classical approach: iterate from low angle, subtracting half of each peak shifted by Δ(2θ) = 2·Δλ/λ·tanθ. Most people should use a tested implementation (for example the one in xrayutilities) rather than write their own — the bookkeeping is unforgiving. What matters here is the check afterwards:

# After stripping, peaks shift slightly. Quantify the shift per peak.
# If positions move by more than ~a step size, the method was too aggressive.
assert np.max(np.abs(positions_after - positions_before)) < 3 * step, \
    "Kalpha2 stripping moved peaks further than expected — inspect it"

Step 3 — find peaks, then fit them properly

scipy.signal.find_peaks is fine for detection. For positions and widths you want a fit, because the measured peak position is the fit centre, not the highest channel:

from scipy.signal import find_peaks
from scipy.optimize import curve_fit

def pseudo_voigt(x, center, width, area, eta):
    g = np.exp(-((x - center) ** 2) / (2 * width ** 2))
    l = 1.0 / (1.0 + ((x - center) / width) ** 2)
    return area * ((1 - eta) * g + eta * l)

peaks, props = find_peaks(corrected, prominence=0.05 * corrected.max(),
                          distance=int(np.ceil(0.2 / step)))
fitted = []
for i, p in enumerate(peaks):
    lo, hi = p - int(0.5 / step), p + int(0.5 / step)
    x, y = two_theta[lo:hi], corrected[lo:hi]
    guess = [two_theta[p], 0.05, y.max(), 0.5]
    popt, _ = curve_fit(pseudo_voigt, x, y, p0=guess,
                        bounds=([x[0], 0.005, 0, 0], [x[-1], 0.5, 1e6, 1]))
    fitted.append(popt)  # center, width, area, eta

Report the fitted centres with their errors, not the raw channel positions. A peak reported to four decimals without a fit is a number pretending to be more precise than the data.

Step 4 — d-spacings, then a reference check (not an identification)

Bragg's law converts fitted positions to d-spacings:

LAMBDA = 1.5406  # Cu Kalpha1, Angstrom

def d_spacing(two_theta_deg, wavelength=LAMBDA):
    theta = np.radians(np.asarray(two_theta_deg) / 2.0)
    return wavelength / (2.0 * np.sin(theta))

ds = d_spacing([f[0] for f in fitted])
for tt, d in zip((f[0] for f in fitted), ds):
    print(f"2theta {tt:7.3f}  ->  d = {d:.4f} A")

Now the honest part. A phase match is: does the measured set of d-spacings (with intensities and their ratios) agree with a reference pattern from COD/ICDD — within the instrument's resolution and the sample's state? The following illustrates the mechanics with a cubic example; do not use the placeholder reference below on real data — pull the actual reference from a database:

# Illustrative reference: cubic phase, lattice parameter a (Angstrom).
# Replace with a real reference pattern (hkl + d + relative intensity).
a = 5.431  # example only
def d_cubic(hkl, a):
    h, k, l = hkl
    return a / np.sqrt(h * h + k * k + l * l)

reference = [(hkl, d_cubic(hkl, a)) for hkl in
             [(1,1,1), (2,0,0), (2,2,0), (3,1,1), (2,2,2), (4,0,0)]]

TOL = 0.02  # Angstrom; tighten/loosen per your 2theta range and broadening
for hkl, d_ref in reference:
    hit = min(ds, key=lambda d: abs(d - d_ref))
    marker = "match" if abs(hit - d_ref) < TOL else "no match"
    print(f"{hkl}: ref d={d_ref:.4f}, closest measured d={hit:.4f}  [{marker}]")

What this produces is *evidence for a hypothesis*. Calling it an identification requires more: relative intensities consistent with the reference, no unexplained lines, and — for anything you will publish or act on — a check against the known polymorphs of your system. Overlapping reflections from a second phase hide exactly this way.

The five ways this pipeline still lies to you

1. Preferred orientation. Plate-like or needle-like crystallites boost or suppress specific reflections. Intensities stop matching any powder reference; positions still look fine. If intensity ratios are off but d-spacings fit, suspect texture, not a new phase.
2. Nanocrystalline broadening. Crystallite size and microstrain both widen peaks, and ordinary fitting cannot separate them without a Williamson–Hall analysis across many reflections.
3. Amorphous content. A broad hump under the pattern is invisible to peak fitting and can be tens of percent of the sample.
4. An unrecorded background method. Two analysts subtract backgrounds differently and report different intensities for the same reflection. The figure is not reproducible without the parameters.
5. Database mismatch. Reference patterns are measured on other instruments; zero-offset, sample displacement and specimen transparency shift your entire pattern. Fit the offset before rejecting a phase.

Free path, provenance path

If you want the parsed numbers without the plumbing, drop an .xye or CSV pattern at /tools/xrd — the free parser returns the series, detected peaks and derived analytics in memory, nothing stored. For a full dataset — many patterns, a reference library, review before modelling — ingestion in a free account keeps the original files as provenance and versions the dataset, which is what makes the background-method question answerable a year later. See /docs/instruments/xrd for the format-specific details.

For the modelling step that usually follows — small-n regression with calibrated uncertainty — start with the conformal prediction primer. It is the difference between a predicted lattice parameter and a predicted lattice parameter with an interval that means something.