A simple Rietveld refinement¶

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

We refine the fluorapatite structure against the pattern notebook 01 fitted with Le Bail. This time the peak intensities come from the atoms, so the fit can move them. Along the way we look at the plan’s stages, the parameter table, two parameters the data cannot tell apart, and a constraint between three atoms.

You need rietx 1.7 or later, which the next cell installs. Notebook 01 introduces the pattern and how to read a fit’s summary.

The data is FAP.XRA from the GSAS-II LabData tutorial, with a CIF transcribed from the same tutorial’s experiment file. Both ship with rietx. Fluorapatite, Ca₅(PO₄)₃F, has seven atomic sites in space group P6₃/m.

Runtime is under a minute on a laptop.

%pip install rietx

The plan and its stages¶

A plan frees parameters in groups, in a fixed order, and runs each group to convergence before freeing the next. The order follows McCusker et al. (1999): background and scale, then zero shift, cell and peak widths, then atomic coordinates and displacement parameters. Freeing everything at once from a poor start tends to land in a wrong minimum.

print(rx.PLAN_INFO["mccusker_structural"].description)
for stage in result.stages:
    print(f"{stage.name:22} {stage.status:10} {len(stage.freed):3} freed")

The standard order continued into the structure: coordinates as site-symmetry DOFs, then displacement parameters (biso or the anisotropic patterns, whichever each site declares), then preferred orientation, extinction and surface roughness where declared.
scale_bkg              converged    7 freed
zero                   converged    1 freed
cell                   converged    2 freed
profile_w              converged    1 freed
profile                converged    4 freed
extra_components       converged    0 freed
coordinates            converged   12 freed
biso                   converged    7 freed
preferred_orientation  converged    0 freed
extinction             converged    1 freed
roughness              converged    0 freed

A stage that frees nothing belongs to a correction this model does not declare, such as preferred orientation. It converges at once and changes nothing.

The parameter table¶

ref.parameters() lists every parameter the model has, whether it varies or not. Here are the rows for one oxygen site, O7, which sits on a general position.

for row in ref.parameters():
    if row.path.startswith("phases.0.atoms.6."):
        state = "varies" if row.vary else (row.held_because or "fixed")
        esd = f"{row.esd:.5f}" if row.esd is not None else "-"
        print(f"{row.path:28} {row.value:10.5f}  esd {esd:8}  {state}")

phases.0.atoms.6.x              0.33986  esd 0.00063   tied: = 0.339862 + 1·phases.0.atoms.6.dof.0
phases.0.atoms.6.y              0.25772  esd 0.00062   tied: = 0.25772 + 1·phases.0.atoms.6.dof.1
phases.0.atoms.6.z              0.07048  esd 0.00061   tied: = 0.0704785 + 1·phases.0.atoms.6.dof.2
phases.0.atoms.6.dof.0          0.00000  esd 0.00063   varies
phases.0.atoms.6.dof.1          0.00000  esd 0.00062   varies
phases.0.atoms.6.dof.2          0.00000  esd 0.00061   varies
phases.0.atoms.6.occ            1.00000  esd -         fixed
phases.0.atoms.6.biso           0.38548  esd 0.11281   varies

The atom’s x, y and z are not refined directly. rietx refines one dof per direction the site symmetry allows. O7 is on a general position, so it has three. An atom on a mirror plane has two, and one on a three-fold axis has one. The coordinates follow from the dof values, so the symmetry always holds.

A constraint between three atoms¶

Fluorapatite has three oxygen sites in its phosphate group: O5, O6 and O7. Their displacement parameters (Biso, the mean-square vibration amplitude times 8π²) are expected to be similar. Check that against the data before constraining anything.

OXYGENS = ["phases.0.atoms.4.biso", "phases.0.atoms.5.biso", "phases.0.atoms.6.biso"]
for path in OXYGENS:
    p = result.parameter(path)
    print(f"{path:24} {p.value:.3f} +/- {p.stderr:.3f} Å²")

phases.0.atoms.4.biso    0.264 +/- 0.149 Å²
phases.0.atoms.5.biso    0.467 +/- 0.158 Å²
phases.0.atoms.6.biso    0.385 +/- 0.113 Å²

Each value lies within about one esd of the others. The data cannot tell them apart, so refining one shared value loses nothing and spends two fewer parameters. tie_equal makes the second and third follow the first. We refit with the same plan, so the tie is the only difference from the first fit.

ref.tie_equal(OXYGENS)
tied = ref.fit(data, plan="mccusker_structural", two_theta_limits=(15, 130))

print(f"free parameters: {result.statistics.n_free_parameters} before, "
      f"{tied.statistics.n_free_parameters} after")
p = tied.parameter(OXYGENS[0])
print(f"shared oxygen Biso: {p.value:.3f} +/- {p.stderr:.3f} Å²")
print(f"Rwp: {result.statistics.rwp:.4f} before, {tied.statistics.rwp:.4f} after")

free parameters: 35 before, 33 after
shared oxygen Biso: 0.377 +/- 0.078 Å²
Rwp: 0.0893 before, 0.0894 after

The shared esd is smaller than any of the three separate ones. That is what the constraint bought. Rwp barely moved, and it could not have shown the gain: a constraint is judged by its premise and its esds.

Two parameters the data cannot separate¶

A flat-plate diffractometer has two corrections that shift every peak: the zero shift (constant in 2θ) and the specimen displacement (proportional to cos θ). The lab_bragg_brentano plan frees both. Run it from where the tied fit ended and see what happens.

both = ref.fit(data, plan="lab_bragg_brentano", two_theta_limits=(15, 130))
for path in ("instrument.zero_shift", "instrument.geometry.sample_displacement", "phases.0.cell.a"):
    p = both.parameter(path)
    print(f"{path:40} {p.value:9.5f} +/- {p.stderr:.5f}")
for d in both.diagnostics:
    if d.code == "HIGH_CORRELATION" and "instrument.zero_shift" in d.where:
        print(f"[{d.level}] {d.code}: {d.message.split(' — ')[0]}")

instrument.zero_shift                     -0.16551 +/- 0.03728
instrument.geometry.sample_displacement   -0.23273 +/- 0.06668
phases.0.cell.a                            9.36897 +/- 0.00098
[warning] HIGH_CORRELATION: instrument.zero_shift ~ instrument.geometry.sample_displacement (ρ=+1.000)
[warning] HIGH_CORRELATION: phases.0.cell.a ~ instrument.zero_shift (ρ=+0.997)
[warning] HIGH_CORRELATION: phases.0.cell.c ~ instrument.zero_shift (ρ=+0.995)

The zero shift and the displacement are correlated at ρ = 1.000. The fit could trade one against the other without changing the pattern, so neither value is a measurement. The cell’s esd grew too, because the cell shifts peaks in a similar way.

The cure is to hold one of the pair at a value known from outside this fit. GSAS’s own refinement of this file held the zero shift at 0 and refined the displacement. We do the same. First go back to the tied fit: every fit is a node in the refinement’s history, and checkout restores one, tie included.

ref.checkout(tied.node_id)
ref.set_values({"instrument.zero_shift": 0.0})
ref.hold(["instrument.zero_shift"])

held = ref.fit(data, plan="lab_bragg_brentano", two_theta_limits=(15, 130))
for path in ("instrument.geometry.sample_displacement", "phases.0.cell.a"):
    p = held.parameter(path)
    print(f"{path:40} {p.value:9.5f} +/- {p.stderr:.5f}")
for d in held.diagnostics:
    if d.code == "HOLD_BLOCKED_PLAN":
        print(f"[{d.level}] {d.code}: {d.message}")

instrument.geometry.sample_displacement    0.06322 +/- 0.00147
phases.0.cell.a                            9.37331 +/- 0.00008
[info] HOLD_BLOCKED_PLAN: 1 held parameter was matched by the turn_on glob of stage 'zero_disp' and stayed fixed: instrument.zero_shift. The hold you declared outranks the plan's glob

The displacement now has a small esd, and so does the cell. HOLD_BLOCKED_PLAN reports that the plan tried to free the zero shift and the hold stopped it.

Use ref.hold for this, not vary=False on the parameter. A plan sets which parameters vary at each stage, so it would free a parameter marked vary=False. A hold outranks the plan.

lab_bragg_brentano refines the profile and the instrument, not the atoms. Refine the structure once more, with the zero shift held, the displacement as just measured, and the tie in place.

final = ref.fit(data, plan="mccusker_structural", two_theta_limits=(15, 130))
print(f"{final.status}  Rwp={final.statistics.rwp:.4f}")

converged  Rwp=0.0890

Is the structure chemically sensible?¶

A fit can match the pattern with atoms in impossible places. result.geometry lists the bonds the refined structure implies, with esds. A phosphate’s P–O bonds are close to 1.54 Å.

for bond in final.geometry.bonds:
    if bond.atom_1 == "P3" and bond.atom_2.startswith("O"):
        print(f"{bond.atom_1}-{bond.atom_2}  {bond.distance:.3f} +/- {bond.stderr:.3f} Å")

P3-O7  1.526 +/- 0.005 Å
P3-O7  1.526 +/- 0.005 Å
P3-O5  1.553 +/- 0.007 Å
P3-O6  1.574 +/- 0.006 Å

The list has one row per bond, and O7 sits on both sides of the mirror plane, so it appears twice. All four bonds are within 0.04 Å of 1.54 Å. That is chemically sensible, but some differ from it by several esds. An esd counts only the noise in the data, and not errors in the model, so a real structure can sit further from a reference than its esd suggests.

The refined structure¶

rx.viz.render_structure draws a structure the way the GUI’s 3D viewer does, from the refined values the model holds now. view says where it is seen from: an axis ("a", "b", "c"), a direction [u, v, w], the normal to a plane {"hkl": (h, k, l)}, or "auto". Looking down c shows what makes this an apatite: fluorine channels at the cell corners, each ringed by Ca2, with the PO₄ tetrahedra and the Ca1 columns between them.

rx.viz.render_structure(ref.structure, view="c")

Atoms outside the cell, drawn to complete a polyhedron or a bond, are shaded darker, so one element shows in two shades. The orange tetrahedra are the phosphate groups.

Down an axis the picture is flat. turn adds small rotations after the view, in degrees about the screen’s x and y axes. Here the default view is tipped 15° about the screen’s x axis, so the cell reads as a box and the tetrahedra show their stacking along c.

rx.viz.render_structure(ref.structure, turn="-15x")

Each call takes a fraction of a second, so try other views. mode="ellipsoid" draws the displacement parameters, atom_labels=True names each site, and path="apatite.png" writes the picture to a file.

The history¶

Every stage of every fit is a node you can return to.

print(ref.history.summary())

tcb8d2285  51 nodes  data=FAP.XRA
 n0000  root                   —
 n0001  stage:scale_bkg        Rwp 0.5704
 n0002  stage:zero             Rwp 0.4556
 n0003  stage:cell             Rwp 0.4545
 n0004  stage:profile_w        Rwp 0.1941
 n0005  stage:profile          Rwp 0.0939
 n0006  stage:extra_components Rwp 0.0897
 n0007  stage:coordinates      Rwp 0.0895
 n0008  stage:biso             Rwp 0.0893
 n0009  stage:preferred_orientation Rwp 0.0893
 n0010  stage:extinction       Rwp 0.0893
 n0011  stage:roughness        Rwp 0.0893
 n0012  set_tie                —
 n0013  stage:scale_bkg        Rwp 0.0895
 n0014  stage:zero             Rwp 0.0895
 n0015  stage:cell             Rwp 0.0895
 n0016  stage:profile_w        Rwp 0.0895
 n0017  stage:profile          Rwp 0.0895
 n0018  stage:extra_components Rwp 0.0895
 n0019  stage:coordinates      Rwp 0.0895
 n0020  stage:biso             Rwp 0.0894
 n0021  stage:preferred_orientation Rwp 0.0894
 n0022  stage:extinction       Rwp 0.0894
 n0023  stage:roughness        Rwp 0.0894
├─  n0024  stage:scale_bkg        Rwp 0.0894
│   n0025  stage:zero_disp        Rwp 0.0894
│   n0026  stage:cell             Rwp 0.0889
│   n0027  stage:profile_w        Rwp 0.0889
│   n0028  stage:profile          Rwp 0.0889
│   n0029  stage:lines_axial      Rwp 0.0883
│   n0030  stage:roughness        Rwp 0.0883
└─  n0031  set_value              —
    n0032  set_hold               —
    n0033  stage:scale_bkg        Rwp 0.3859
    n0034  stage:zero_disp        Rwp 0.0986
    n0035  stage:cell             Rwp 0.0897
    n0036  stage:profile_w        Rwp 0.0897
    n0037  stage:profile          Rwp 0.0897
    n0038  stage:lines_axial      Rwp 0.0890
    n0039  stage:roughness        Rwp 0.0890
    n0040  stage:scale_bkg        Rwp 0.0890
    n0041  stage:zero             Rwp 0.0890
    n0042  stage:cell             Rwp 0.0890
    n0043  stage:profile_w        Rwp 0.0890
    n0044  stage:profile          Rwp 0.0890
    n0045  stage:extra_components Rwp 0.0890
    n0046  stage:coordinates      Rwp 0.0890
    n0047  stage:biso             Rwp 0.0890
    n0048  stage:preferred_orientation Rwp 0.0890
    n0049  stage:extinction       Rwp 0.0890
   *n0050  stage:roughness        Rwp 0.0890

Read it from the top. A straight run of stages stays in one column. Where we checked out an earlier node and fitted again, the history forks, and each branch is indented one level: ├─ for the first and └─ for the last. * marks the node the model stands at now.

The verdict for a declared purpose¶

What counts as finished depends on what the fit is for. ref.summary prints the rows that decide a given deliverable. Here the deliverable is the structure.

print(ref.summary(deliverable="structure"))

Refinement.summary: converged (rietveld)
  stage scale_bkg: converged (2 evals, ftol=1e-06)
  stage zero: converged (1 evals, ftol=1e-06)
  stage cell: converged (2 evals, ftol=1e-06)
  stage profile_w: converged (2 evals, ftol=1e-06)
  stage profile: converged (4 evals, ftol=1e-06)
  stage extra_components: converged (4 evals, ftol=1e-06)
  stage coordinates: converged (5 evals, ftol=1e-06)
  stage biso: converged (5 evals, ftol=1e-06)
  stage preferred_orientation: converged (3 evals, ftol=1e-06)
  stage extinction: converged (1 evals, ftol=1e-06)
  stage roughness: converged (1 evals, ftol=solver default)
  diagnostics: 5 unresolved
    WARNING SITE_SNAPPED_TO_SPECIAL_POSITION: 1 site(s) in phase 'fluorapatite' sit within 0.0001 of a special position of 'P 63/m' without being on it, and were expanded at it: Ca1 (3.3e-11 → multiplicity 4). Largest shift 3.33e-11 in fractional coordinates — the stored coordinates are unchanged and the fit is unaffected, but the multiplicity is the special one — it is what decides how many atoms the site puts in the cell, hence ZMV and every weight fraction. Check these sites against the source's own multiplicities: if they agree, nothing is wrong and the coordinates are simply rounded; if they disagree, the coordinates and the space group are telling you different things
    WARNING FLOOR_HIT: phases.0.extinction refined to its floor (ρ=+0.064) — the data push this parameter below zero, where it cannot go. The usual reasons are a phase or term the data do not support, a width the instrument already supplies, or, for a displacement parameter, absorption or a wrong species or occupancy. Fix it at 0
    WARNING RESOLUTION_UNCONSTRAINED: U, V, W refined on a predominantly Lorentzian pattern (Γ_L exceeds Γ_G at 4951 of 5751 fitted points; the Gaussian carries 41% of the width) — on a pattern of this character the Gaussian resolution terms are not determined by the data, so what they converged to is a fitted artefact rather than a measurement of the instrument — measure the instrument instead of refining it here: lab_calibrate on a standard with its certified cell held fixed, save_instrument_profile, then load_instrument_profile before this fit, which holds U, V and W at the instrumental values. Do not quote the widths, and treat any crystallite size or strain read off this profile as unmeasured
    INFO HOLD_BLOCKED_PLAN: 1 held parameter was matched by the turn_on glob of stage 'zero' and stayed fixed: instrument.zero_shift. The hold you declared outranks the plan's glob — this is the hold working. If you meant the plan to refine these, unhold() them first; if you meant them held, the plan is doing less than its name suggests and a narrower one would say so
    WARNING PATTERN_UNDERSAMPLED: 4.6 steps across the FWHM (median over 76 fitted peaks) — below the 5 the guidelines ask for — a data-collection limit, not a refinement one: the integrated intensities were not measured finely enough and no choice made here recovers them. Re-collect at a step size near FWHM/5 if the intensities have to be quotable
  layer0/1: Rwp=0.0890 GoF=1.59; 63 regions, top 15 shown (47% of χ²); 4 unmatched observed peak(s); remaining 48 regions carry 30% of χ²; Layer 1 abstained — regions carrying 40% of χ² misfit, but only 20% of that sits in regions the local gates accept (need 40%); Layer 1 abstains
    abstained: unreadable — regions carrying 40% of χ² misfit, but only 20% of that sits in regions the local gates accept (need 40%); Layer 1 abstains
    region 31.5-32.5°: 9% of χ², max|Δ|/σ=12.7
    region 103.2-111.3°: 4% of χ², max|Δ|/σ=5.0
    region 32.6-33.4°: 4% of χ², max|Δ|/σ=9.6
    unmatched obs peaks: 4
    serial correlation (DW): 1.03, esds inflated ×2.17
  next: free phases.0.atoms.2.occ, predicted ΔBIC +31.4 at N_eff 1220 (Δχ² 449)
  deliverable: structure
    exchangeability: 2 row(s)
      instrument.geometry.sample_displacement ~ phases.0.cell.a, phases.0.cell.c, instrument.profile.v, instrument.profile.w (R²=0.68)
      instrument.geometry.sample_transparency ~ phases.0.cell.a, phases.0.cell.c (R²=0.90)
  protocol:
    plan: scale_bkg, zero, cell, profile_w, profile, extra_components, coordinates, biso, preferred_orientation, extinction, roughness
    held: 16 symmetry-locked, 17 user-fixed, 2 wavelength-blocked
    N points=5751, N reflections=325, effective obs=143
    σ source: Poisson √max(y,1)
  provenance: rietx 1.10.0.dev0, backend=numpy, solver=lm
  Rwp 0.0890 / Rexp 0.0559 (GoF 1.59), Rp 0.0693, χ² 2.5, DW 1.03
  visual check: call summary(plot='fit.png') or plot_for_vlm(result, report) to draw these regions

The summary repeats the stages and diagnostics, then adds what needs the whole model.

layer0/1 reads the difference curve region by region and ranks the regions by their share of χ². When too little of the misfit sits in regions it can explain, it abstains from attributing a cause. An abstention is an answer: no single correction is indicated.
unmatched obs peaks counts observed peaks with no reflection under them. Find them before quoting a structure, because a weak second phase would explain them.
next names the parameter whose release the package predicts would help most, and the ΔBIC that release would have to pass.
exchangeability lists held parameters whose effect a free one could imitate. The agent skill’s rule 14 settles each with two short fits, one per rival.
protocol is what you report alongside the numbers: the plan, what was held, and the observation count.

The final fit¶

final.plot()

Saving the work as a project¶

Everything so far lives in this Python session and goes when it ends. A project keeps a refinement on disk as a directory whose name ends .rex. The GUI opens the same directory (rietx gui fluorapatite.rex), so a project started in code can be continued by hand, and the reverse.

Project.create takes the pattern file, the models and the protocol, and starts a history of its own. project.fit() refines and records each stage there.

work = Path(tempfile.mkdtemp())
project = rx.Project.create(work / "fluorapatite.rex", pattern=examples_dir() / "FAP.XRA",
                            structure=ref.structure, instrument=ref.instrument,
                            plan="mccusker_structural", two_theta_limits=(15, 130))
saved = project.fit()
project.save()
print(saved.status)
for item in sorted((work / "fluorapatite.rex").iterdir()):
    print(f"{item.name + ('/' if item.is_dir() else ''):16} {item.stat().st_size if item.is_file() else ''}")
print(sorted(json.loads((work / "fluorapatite.rex" / "project.json").read_text(encoding="utf-8"))))

converged
FAP.XRA          47396
exports/         
history.jsonl    163395
live/            
project.json     4967
['created_utc', 'excluded_regions', 'format_version', 'history_file', 'indexing', 'mode', 'package_version', 'patterns', 'plan', 'two_theta_limits', 'ui', 'updated_utc']

Each entry holds one kind of fact, and no fact is stored twice.

FAP.XRA is the pattern, copied byte for byte. project.json records its checksum, which reader opened it, and a fingerprint of the numbers that reader produced. A changed file and a changed reader can then be told apart.
project.json holds the settings: the plan, the fitted range, excluded regions, and the GUI’s own choices. It holds no refined values.
history.jsonl holds every node of the history, one JSON line each, appended as each stage finishes. The node at its head is the working state, so the refined values live here.
live/ holds each run’s record while it happens, which rietx watch and the GUI draw from.
exports/ is where exported files go, such as a CIF of the refined structure.

rx.Project.open reopens a project in code. Every verb writes into the directory, so to look at someone’s project without changing it, open a copy (rietx gui --scratch makes one).

Checking an agent’s work¶

When an agent reports a structure refinement, check these in order. The numbers are the rules in the rietx agent skill (SKILL.md §4).

Status and every diagnostic first (rule 9). Ask which diagnostics limit what it quotes.
Correlations near 1 (§3, rule 7). Two parameters at ρ ≈ 1 are one measurement. Ask which was held, and what value it was held at.
Chemistry (rule 12). Bond lengths and Biso values should be physically possible.
Constraints (§3, rule 8). A tie needs its premise checked in an unconstrained fit first.
Rwp last (rule 16), and only against a fit of the same channels.