Peaks, indexing and a Le Bail check

Written by Claude Code, Anthropic’s coding agent, for the rietx project.

When the phase is unknown there is no CIF to start from. The cell has to come from the peak positions alone: find the peaks, then search for a cell that puts a reflection under each one. That search is called indexing. We index a synthetic pattern whose answer we know, read why rietx declines to name a single cell, and check the best candidate with a Le Bail fit.

You need rietx 1.7 or later, which the next cell installs. Notebook 01 introduces Le Bail fitting.

The data is synthetic aragonite, CaCO₃, made from its published structure: Dal Negro and Ungaretti (1971), American Mineralogist 56, 768, as entry 9000229 of the Crystallography Open Database (public domain). The displacement parameters are the isotropic equivalents of that entry’s anisotropic ones.

Runtime is under a minute on a laptop, most of it the indexing search.

%pip install rietx

A synthetic pattern

import tempfile
from pathlib import Path

import numpy as np

import rietx as rx
from rietx.indexing import SearchSpec, structure_from_candidate

work = Path(tempfile.mkdtemp())
cif = work / "aragonite.cif"
cif.write_text(
    "data_aragonite\n_space_group_name_H-M_alt 'P m c n'\n"
    "_cell_length_a 4.9616\n_cell_length_b 7.9705\n_cell_length_c 5.7394\n"
    "_cell_angle_alpha 90\n_cell_angle_beta 90\n_cell_angle_gamma 90\n"
    "loop_\n_atom_site_label\n_atom_site_type_symbol\n_atom_site_fract_x\n"
    "_atom_site_fract_y\n_atom_site_fract_z\n_atom_site_B_iso_or_equiv\n"
    "Ca1 Ca 0.25 0.4151 0.2403 0.67\n"
    "C1  C  0.25 0.7627 0.0850 0.80\n"
    "O1  O  0.25 0.9231 0.0952 1.25\n"
    "O2  O  0.4729 0.6801 0.0870 1.19\n",
    encoding="utf-8")
truth = rx.Structure.from_cif(cif)

true_instrument = rx.Instrument.bragg_brentano(radiation="CuKa")
true_instrument.profile.w.value = 0.004
true_instrument.profile.x.value = 0.05
true_instrument.background = rx.BackgroundChebyshev.with_terms(1)
true_instrument.background.coefficients[0].value = 200.0

two_theta = np.arange(10, 80, 0.02)
y = rx.Refinement(truth, true_instrument).predict(two_theta) - 200.0
truth.phases[0].scale.value *= 20000 / y.max()   # strongest peak about 20 000 counts
y = rx.Refinement(truth, true_instrument).predict(two_theta)
data = rx.PatternData(two_theta=two_theta,
                      intensity=np.random.default_rng(0).poisson(y).astype(float))
data.plot()

From here on we pretend not to know the structure. The true cell is a = 4.9616, b = 7.9705, c = 5.7394 Å, and the indexer should find it in some order of axes.

Find the peaks

pick_peaks finds every resolvable peak and fits each one for its position, with an esd. It needs an instrument for the wavelength and a first guess at the peak width. We give it a default copper instrument, as you would for an unknown sample.

lab = rx.Instrument.bragg_brentano(radiation="CuKa")
peaks = rx.pick_peaks(data, lab)
peaks

PeakList: 105 peaks, 59 usable, 10–79.98° 2θ, λ = 1.54059 Å, source fitted

#	two_theta	two_theta_esd	intensity	fwhm	origin	flags
0	21.0739	0.0024	31.4	0.101	fitted
1	22.2946	0.0056	9.6	0.084	fitted
2	26.2221	0.0003	2175.4	0.096	fitted	asymmetry_unmodelled
3	26.6192	0.0106	18.1	0.095	fitted	asymmetry_unmodelled
4	27.2200	0.0004	1241.7	0.095	fitted	asymmetry_unmodelled
5	31.1394	0.0016	92.8	0.110	fitted
6	32.7501	0.0010	170.8	0.096	fitted	asymmetry_unmodelled
7	32.7600	0.0019	168.7	0.097	fitted	position_at_bound, asymmetry_unmodelled, duplicate_line
8	33.1523	0.0003	1066.9	0.096	fitted	asymmetry_unmodelled
9	33.1524	0.0006	1065.8	0.097	fitted	asymmetry_unmodelled, duplicate_line
10	36.1488	0.0006	946.2	0.125	fitted
11	37.2578	0.0025	324.2	0.098	fitted	duplicate_line
12	37.2651	0.0011	298.2	0.100	fitted	asymmetry_unmodelled
13	37.8987	0.0335	953.5	0.117	fitted	duplicate_line
14	37.9064	0.0043	849.3	0.100	fitted	asymmetry_unmodelled, duplicate_line
15	37.9071	0.0004	870.7	0.097	fitted
16	37.9074	0.0005	867.3	0.098	fitted	asymmetry_unmodelled, duplicate_line
17	38.3949	0.0069	554.1	0.117	fitted	duplicate_line
18	38.4024	0.0005	620.7	0.097	fitted
19	38.4033	0.0007	644.1	0.098	fitted	asymmetry_unmodelled, duplicate_line
20	38.5784	0.0098	378.0	0.117	fitted	position_at_bound, duplicate_line
21	38.6312	0.0007	401.9	0.097	fitted
22	41.1949	0.0008	270.7	0.094	fitted	asymmetry_unmodelled
23	42.9070	0.0008	426.4	0.100	fitted
24	45.4833	0.0021	71.2	0.098	fitted
25	45.5000	0.0038	70.0	0.099	fitted	position_at_bound, asymmetry_unmodelled, duplicate_line
26	45.8568	0.0015	1462.2	0.101	fitted	asymmetry_unmodelled, duplicate_line
27	45.8588	0.0003	1531.2	0.098	fitted
28	45.8589	0.0005	1528.6	0.099	fitted	asymmetry_unmodelled, duplicate_line
29	46.0134	0.0137	39.2	0.101	fitted	asymmetry_unmodelled, kalpha2_residual
30	46.3251	0.1170	1.3	0.101	fitted	position_at_bound, no_intensity, asymmetry_unmodelled
31	48.3096	0.0007	639.6	0.100	fitted	fit_failed
32	48.3097	0.0007	632.9	0.100	fitted
33	48.4600	0.0011	419.8	0.100	fitted	position_at_bound
34	48.4639	405.0436	130.1	0.100	fitted	fit_failed, unresolved_shoulder, no_intensity, position_unmeasured, asymmetry_unmodelled, kalpha2_residual
35	48.4639	186.4488	282.6	0.100	fitted	fit_failed, unresolved_shoulder, no_intensity, position_unmeasured, asymmetry_unmodelled
36	49.8839	0.0062	31.2	0.099	fitted	asymmetry_unmodelled, duplicate_line
37	49.8849	0.0060	32.7	0.100	fitted	asymmetry_unmodelled
38	50.2397	0.0007	621.6	0.100	fitted	duplicate_line
39	50.2398	0.0007	628.0	0.099	fitted	asymmetry_unmodelled
40	51.9129	0.0038	57.9	0.101	fitted	asymmetry_unmodelled
41	52.4908	0.0013	657.3	0.101	fitted	asymmetry_unmodelled, duplicate_line
42	52.4910	0.0019	671.8	0.098	fitted	asymmetry_unmodelled, duplicate_line
43	52.4911	0.0004	662.8	0.100	fitted
44	52.9106	0.0012	143.9	0.100	fitted
45	52.9111	0.0064	141.6	0.098	fitted	asymmetry_unmodelled, duplicate_line
46	52.9918	0.0297	490.7	0.240	fitted	fit_failed, asymmetry_unmodelled
47	53.0543	0.0007	298.4	0.100	fitted
48	53.0549	0.0042	307.2	0.098	fitted	asymmetry_unmodelled, duplicate_line
49	53.3713	0.0168	9.2	0.110	fitted	asymmetry_unmodelled
50	53.9571	0.0025	45.8	0.110	fitted	asymmetry_unmodelled
51	56.1494	0.0017	55.1	0.102	fitted	duplicate_line
52	56.1494	0.0015	55.9	0.101	fitted
53	56.2020	0.0125	19.7	0.083	fitted
54	56.2973	0.0056	27.3	0.083	fitted
55	56.4387	0.0040	18.1	0.101	fitted
56	56.4391	0.0043	18.4	0.102	fitted	duplicate_line
57	56.7956	0.0208	23.9	0.102	fitted	duplicate_line
58	56.8086	0.0033	39.3	0.083	fitted	asymmetry_unmodelled, duplicate_line
59	56.8089	0.0021	35.6	0.101	fitted
60	59.2447	0.0017	152.5	0.101	fitted	asymmetry_unmodelled
61	59.6004	0.0158	7.7	0.107	fitted
62	60.2030	0.0029	35.9	0.107	fitted
63	61.8198	0.0012	120.7	0.107	fitted
64	62.9101	0.0019	58.1	0.106	fitted	asymmetry_unmodelled
65	62.9102	0.0020	58.1	0.107	fitted	asymmetry_unmodelled, duplicate_line
66	63.3338	0.0011	149.6	0.106	fitted	asymmetry_unmodelled, duplicate_line
67	63.3340	0.0010	148.8	0.107	fitted
68	65.0381	0.0091	8.2	0.090	fitted
69	66.1143	0.0024	77.9	0.097	fitted
70	66.1486	0.0042	124.4	0.154	fitted	duplicate_line
71	66.2126	0.0031	61.7	0.097	fitted
72	66.5357	0.0064	72.1	0.154	fitted	asymmetry_unmodelled, duplicate_line
73	66.5444	0.0020	72.2	0.097	fitted
74	68.6427	0.0015	74.8	0.098	fitted
75	68.6430	0.0015	73.8	0.100	fitted	asymmetry_unmodelled, duplicate_line
76	69.1001	0.0012	104.4	0.099	fitted	duplicate_line
77	69.1006	0.0012	102.8	0.098	fitted
78	69.1011	0.0013	102.4	0.100	fitted	asymmetry_unmodelled, duplicate_line
79	69.5874	0.0026	39.5	0.099	fitted
80	69.5875	0.0030	39.8	0.109	fitted	asymmetry_unmodelled, duplicate_line
81	69.5885	0.0027	40.8	0.098	fitted	duplicate_line
82	70.0819	0.0024	51.0	0.109	fitted	asymmetry_unmodelled, duplicate_line
83	70.0830	0.0023	54.3	0.099	fitted
84	70.2870	0.0165	12.4	0.114	fitted
85	70.8747	0.0050	35.2	0.114	fitted
86	72.2882	0.0039	15.7	0.058	fitted	asymmetry_unmodelled
87	72.2900	0.0055	11.8	0.088	fitted	duplicate_line
88	72.4957	0.0126	5.1	0.088	fitted
89	75.2903	0.0010	148.5	0.108	fitted
90	76.6601	0.0031	63.8	0.106	fitted	asymmetry_unmodelled
91	76.6603	0.0034	63.5	0.106	fitted	asymmetry_unmodelled, duplicate_line
92	76.6604	0.0033	63.5	0.106	fitted	asymmetry_unmodelled, duplicate_line
93	76.7802	0.0021	115.9	0.106	fitted
94	76.7802	0.0022	116.0	0.106	fitted	asymmetry_unmodelled, duplicate_line
95	76.7803	0.0022	116.2	0.106	fitted	asymmetry_unmodelled, duplicate_line
96	77.0997	0.0011	191.0	0.106	fitted
97	77.0998	0.0012	192.1	0.106	fitted	asymmetry_unmodelled, duplicate_line
98	77.1007	0.0013	193.8	0.106	fitted	asymmetry_unmodelled, duplicate_line
99	78.0360	0.0021	123.7	0.115	fitted	duplicate_line
100	78.0376	0.0015	122.0	0.116	fitted
101	78.6949	0.0052	20.8	0.115	fitted
102	78.7385	0.0473	64.8	0.116	fitted	position_at_bound, no_intensity
103	79.3343	0.0163	35.7	0.115	fitted
104	79.4033	0.0012	208.5	0.123	fitted	asymmetry_unmodelled

  [ info  ] PEAK_UNRESOLVED_SHOULDER: 2 line(s) never separated from a neighbour by half a FWHM, so their positions are correlated with it and their σ says so
           at: 48.4639°, 48.4639°
  [ info  ] PEAK_KALPHA2_RESIDUAL: 2 weak component(s) sit at a strong group-mate's predicted Kα2 maximum at a few per cent of its area — the residual of a *modelled* doublet, re-created by a re-seed pass after detection dropped the alias candidate (WP-1043)
           at: 46.0134°, 48.4639°
  [warning] PEAK_ASYMMETRY_UNMODELLED: 49 line(s) leave a significant odd residual after fitting, i.e. an asymmetry the declared profile does not carry; an unmodelled one-sided aberration biases the centroid in one direction, which σ cannot see
           at: 26.2221°, 26.6192°, 27.2200°, 32.7501°, 32.7600°, 33.1523°, 33.1524°, 37.2651°, 37.9064°, 37.9074°, 38.4033°, 41.1949°, 45.5000°, 45.8568°, 45.8589°, 46.0134°, 46.3251°, 48.4639°, 48.4639°, 49.8839°, 49.8849°, 50.2398°, 51.9129°, 52.4908°, 52.4910°, 52.9111°, 52.9918°, 53.0549°, 53.3713°, 53.9571°, 56.8086°, 59.2447°, 62.9101°, 62.9102°, 63.3338°, 66.5357°, 68.6430°, 69.1011°, 69.5875°, 70.0819°, 72.2882°, 76.6601°, 76.6603°, 76.6604°, 76.7802°, 76.7803°, 77.0998°, 77.1007°, 79.4033°
  [warning] PEAK_WIDTH_LAW_MISMATCH: measured FWHM 0.1565° against 0.0322° from the instrument's U,V,W,X,Y — a factor of 4.9; seeds and windows were scaled by 4.86 to match the data
           at: median of the 12 most prominent lines
  [ info  ] PEAK_KALPHA2_ALIAS: 26 candidate(s) dropped as the Kα2 maximum of a stronger line rather than lines of their own; the doublet is fitted as a constrained pair, so each reported position is a Kα1 position
           at: 33.2400°, 38.0200°, 38.5000°, 41.3000°, 43.0200°, 45.9800°, 48.6000°, 50.3800°, 52.6200°, 53.2000°, 56.3000°, 59.4000°, 60.3600°, 62.0000°, 63.0800°, 63.5200°, 66.7200°, 68.8400°, 69.2800°, 69.7600°, 70.2800°, 71.0800°, 75.5000°, 77.3200°, 78.2600°, 79.6400°
  [ info  ] PEAK_SHOULDER_SEEDED: 1 curvature seed(s) offered for components with no local maximum; survival was decided by ΔBIC, not by detection
           at: 51 groups

Each row is one peak: its position, the esd of that position, its area and its width. flags records anything the fit noticed. A flagged peak can still be usable, and the count line says how many are.

The diagnostics under the table explain the list. PEAK_WIDTH_LAW_MISMATCH says the peaks are about five times wider than the default instrument predicts, so the picker scaled its windows to match. PEAK_KALPHA2_ALIAS counts candidates it dropped as the Kα2 partner of a stronger line: each reported position is a Kα1 position.

Measure lines you name

fit_peaks fits exactly the positions you give it, and nothing else. Use it to measure lines you have chosen: a few reflections for a size analysis, or one line as a quick check.

chosen = rx.fit_peaks(data, lab, [26.22, 27.22, 33.13])
chosen

PeakList: 3 peaks, 3 usable, 10–79.98° 2θ, λ = 1.54059 Å, source fitted

#	two_theta	two_theta_esd	intensity	fwhm	origin	flags
0	26.2221	0.0003	2175.4	0.096	manual	asymmetry_unmodelled
1	27.2199	0.0006	1246.2	0.094	manual	asymmetry_unmodelled, unnamed_neighbour
2	33.1514	0.0031	1113.5	0.092	manual	asymmetry_unmodelled, unnamed_neighbour

  [warning] PEAK_ASYMMETRY_UNMODELLED: 3 line(s) leave a significant odd residual after fitting, i.e. an asymmetry the declared profile does not carry; an unmodelled one-sided aberration biases the centroid in one direction, which σ cannot see
           at: 26.2221°, 27.2199°, 33.1514°
  [warning] PEAK_WIDTH_LAW_MISMATCH: measured FWHM 0.1565° against 0.0322° from the instrument's U,V,W,X,Y — a factor of 4.9; seeds and windows were scaled by 4.86 to match the data
           at: median of the 12 most prominent lines
  [ info  ] PEAK_KALPHA2_ALIAS: 26 candidate(s) dropped as the Kα2 maximum of a stronger line rather than lines of their own; the doublet is fitted as a constrained pair, so each reported position is a Kα1 position
           at: 33.2400°, 38.0200°, 38.5000°, 41.3000°, 43.0200°, 45.9800°, 48.6000°, 50.3800°, 52.6200°, 53.2000°, 56.3000°, 59.4000°, 60.3600°, 62.0000°, 63.0800°, 63.5200°, 66.7200°, 68.8400°, 69.2800°, 69.7600°, 70.2800°, 71.0800°, 75.5000°, 77.3200°, 78.2600°, 79.6400°
  [ info  ] PEAK_SHOULDER_SEEDED: 1 curvature seed(s) offered for components with no local maximum; survival was decided by ΔBIC, not by detection
           at: 51 groups

Rebuilding a whole list from edited positions in code works less well: a real line you leave out still has intensity, and its neighbours absorb it.

The same work in the GUI

Editing a peak list is easier by hand, and rietx has a graphical mode for it. Run rietx gui in a terminal, open or create a project with your pattern, and choose the Peaks tab. The plot then becomes an editing surface:

Gesture	Does
click empty space	add a line, measured by fit_peaks
drag a marker	move it
shift-click	exclude a line from indexing
right-click	remove it

The same tab runs the indexing search, draws a candidate’s predicted lines over the pattern, and adopts a candidate as the starting point of a Le Bail fit. Every gesture also has a typed equivalent in the panel. The GUI guide covers the tab in full.

Index

index_pattern searches for cells that index the usable peaks. A search over every crystal system can take minutes, so we narrow it. We search only orthorhombic cells, with a volume between 200 and 250 ų and no axis longer than 8.5 Å. Those limits come from knowing the answer. With a real unknown, start from the formula and a wide window, and let the search run longer.

spec = SearchSpec(systems=("orthorhombic",), min_volume=200, max_volume=250, max_d_axis=8.5)
indexing = rx.index_pattern(peaks, data=data, instrument=lab, spec=spec)
indexing

engines: dichotomy, svd, trial_error   systems: orthorhombic   lines: 59   validated: yes
search INCOMPLETE in: orthorhombic — a negative result there is not evidence

#	confidence	lattice	a b c (Å)	α β γ (°)	V (ų)	indexed	chi2_red	Le Bail Rwp	confidence_caveats
1	low	orthorhombic P	4.96150 5.73930 7.97038	90.000 90.000 90.000	226.96	50/59	0.22	0.2570	predicted_but_absent, indexed_fraction_low, search_incomplete

NO CELL: the result abstains — see the diagnostics below.
  [ info  ] PEAK_POSITION_PRECISION: median σ(2θ) 0.0019°, worst 0.0168°; median σ(Q)/Q = 5.76e-05 and σ(Q) is 0.002 of the mean line spacing
           - 2θ 10.00-79.98°
  [warning] INDEX_SEARCH_INCOMPLETE: the svd search did not finish orthorhombic within 30 s per system, so finding no cell there is not evidence that none exists
           - orthorhombic
  [ info  ] INDEX_SHIFT_FROM_PAIRS: a systematic shift of 0.0052° 2θ was measured from 8 harmonic reflection pairs agreeing to 0.0049° (of 17 admitted, z = 6.6 against 200 structureless replicates, p = 0.005) — no reference positions were needed; the cause is not named: no template was refuted
           - allowance = 0.0052° 2θ (measured)
           - best template = cos_theta
  [warning] INDEX_ABSTAINED: no cell reached the confidence gate; the best candidate (orthorhombic P, V = 227.0 ų) is low because of: predicted_but_absent, indexed_fraction_low, search_incomplete
           - searched: orthorhombic
  [ info  ] INDEX_SYSTEMS_NOT_COVERED: 6 crystal system(s) were not searched, so this result says nothing about them
           - cubic
           - hexagonal
           - trigonal
           - tetragonal
           - monoclinic
           - triclinic
  [ info  ] INDEX_CELL_SYSTEMATIC_UNQUANTIFIED: a Bragg-Brentano cell carries ≈ ±85 ppm that no esd reports: sweeping the goniometer radius over 180-320 mm moves Rwp by only 0.029 points, so the data cannot identify it, while specimen displacement absorbs the change 4.6× and the rest lands on the lattice parameters
           - R = 217.5 mm (declared, not measured by the fit)
  [warning] INDEX_PREDICTED_BUT_ABSENT (candidate orthorhombic P V=227.0): 4 of 105 reflections this lattice predicts have no net intensity above the fitted background — the oversized-cell signature, and the one M20 cannot see (its N_poss denominator penalises it only weakly, Oishi-Tomiyasu 2013)
  [warning] INDEX_IMPURITY_LINES (candidate orthorhombic P V=227.0): 96 observed peak(s) have no calculated reflection nearby after a whole-pattern Le Bail fit of this cell

The top candidate is the true cell, with its axes in a different order. All three search engines found it.

Yet the result abstains: best_or_none() returns None, and the verdict line says NO CELL. rietx names a single cell only when exactly one candidate reaches high confidence, and this one is low for the reasons in its row.

• predicted_but_absent: the indexer checks a candidate in the lattice’s highest-symmetry space group, which has no systematic absences. Aragonite’s real space group, Pmcn, has glide planes that remove some reflections. Those reflections are predicted but absent.

• indexed_fraction_low: some usable lines are not indexed. A few of them are weak false lines the picker found beside strong peaks.

The candidate’s INDEX_IMPURITY_LINES warning comes from the indexer’s own quick Le Bail fit, whose Rwp in the table is poor. A poor fit leaves many peaks unmatched, so that count is not evidence of a second phase. Our own Le Bail fit below settles it.

print(indexing.best_or_none())
candidate = indexing.candidates[0]
print(candidate.cell, candidate.confidence_caveats)

None
(4.961504341308777, 5.7392991240231845, 7.9703799717669055, 90.0, 90.0, 90.0) ['predicted_but_absent', 'indexed_fraction_low', 'search_incomplete']

Check the candidate with Le Bail

structure_from_candidate turns a candidate into a structure for Le Bail fitting. It uses the lattice’s highest-symmetry space group (here P m m m), so the fit tests the lattice and nothing else. As in notebook 01, we start the background at the data’s floor and run the plan twice.

structure = structure_from_candidate(candidate)
print(structure.phases[0].space_group)

instrument = rx.Instrument.bragg_brentano(radiation="CuKa")
instrument.background = rx.BackgroundChebyshev.with_terms(4)
instrument.background.coefficients[0].value = float(np.percentile(data.intensity, 5))
ref = rx.Refinement(structure, instrument)
passes = [ref.fit(data, mode="lebail", plan="profile_only") for _ in range(2)]
best = min(passes, key=lambda r: r.statistics.rwp)
print(f"Rwp by pass: {[round(r.statistics.rwp, 4) for r in passes]}")
for edge in "abc":
    p = best.parameter(f"phases.0.cell.{edge}")
    print(f"{edge} = {p.value:.5f} +/- {p.stderr:.5f} Å")

P m m m

Rwp by pass: [0.0707, 0.0607]
a = 4.96156 +/- 0.00008 Å
b = 5.73953 +/- 0.00009 Å
c = 7.97054 +/- 0.00013 Å

best.plot()

The fit puts a reflection under every peak and leaves only small residuals, so the lattice is right. The refined edges match the true cell, in the candidate’s order of axes.

Some ticks stand over no peak, such as those near 11° and 15.5°. They are reflections P m m m allows and aragonite’s glide planes remove: the indexer’s predicted-but-absent reflections.

What this does not settle is the space group. That is a question about which reflections are absent, and it comes after the cell is refined.

Checking an agent’s work

When an agent reports a cell it indexed, check these. The rules are in the rietx agent skill (SKILL.md §6, and references/diagnostics-indexing.md).

• Did best_or_none() return a cell? If not, the agent chose one. Ask why, and read the chosen candidate’s confidence_caveats (rule 19).

• What was searched? A search narrowed by system, volume or axis length says nothing about the cells outside it. INDEX_SYSTEMS_NOT_COVERED and INDEX_SEARCH_INCOMPLETE say what was left out.

• Was the cell checked with a whole-pattern fit? A Le Bail fit of the candidate should leave no peak unexplained.

• Is the space group a separate claim? A cell from indexing comes with the lattice’s highest-symmetry group. The real space group needs its own evidence.