Steering a refinement#
We take the fluorapatite refinement from notebook 02 and change what it refines. First we tie three displacement parameters together. Then we run a plan that frees two parameters the data cannot separate, go back in the history, and hold one of them. By the end you can tie, go back and hold, and read the verdict rietx gives for a structure.
Written by Claude Code, Anthropic’s coding agent, for the rietx project.
You need rietx 1.10 or later, which the next cell installs. Notebook 02 introduces the first fit and how to read its parameters.
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.
Runtime is under 30 seconds on a laptop.
%pip install rietx
The starting fit#
This cell repeats notebook 02’s first refinement, which that notebook reads line by line.
import rietx as rx
from rietx.examples import examples_dir
data = rx.read_pattern(examples_dir() / "FAP.XRA")
structure = rx.Structure.from_cif(examples_dir() / "fluorapatite.cif")
instrument = rx.Instrument.bragg_brentano(radiation="CuKa")
instrument.geometry.axial_sl.value = 0.02
instrument.geometry.axial_hl.value = 0.02
instrument.background = rx.BackgroundChebyshev.with_terms(6)
ref = rx.Refinement(structure, instrument)
result = ref.fit(data, plan="mccusker_structural", two_theta_limits=(15, 130))
result.status
'converged'
A constraint between three atoms#
Fluorapatite has three oxygen sites in its phosphate group: O5, O6 and O7. Their Biso values 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",
]
[result.parameter(path) for path in OXYGENS]
[phases.0.atoms.4.biso 0.24(15) vary,
phases.0.atoms.5.biso 0.47(16) vary,
phases.0.atoms.6.biso 0.36(11) vary]
The three values agree within their esds.
No two differ by much more than the esd of their difference.
The data cannot tell them apart, so one shared value loses nothing and refines two fewer parameters.
tie_equal makes the second and third follow the first.
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))
[tied.parameter(path) for path in OXYGENS]
[phases.0.atoms.4.biso 0.360(79) vary,
phases.0.atoms.5.biso 0.360(79) tied,
phases.0.atoms.6.biso 0.360(79) tied]
The three rows now carry one value, and two of them read tied.
The shared esd is smaller than any of the three separate ones.
Rwp says little about a constraint.
Judge one by its premise and its esds.
How to tie parameters together shows the check in full.
Two parameters the data cannot separate#
A flat-plate diffractometer has two corrections that shift every peak.
The zero shift is constant in 2θ, and the specimen displacement goes as cos θ.
The lab_bragg_brentano plan frees both.
Run it from where the tied fit ended, and expect a warning about the pair.
SHIFTS = [
"instrument.zero_shift",
"instrument.geometry.sample_displacement",
"phases.0.cell.a",
]
both = ref.fit(data, plan="lab_bragg_brentano", two_theta_limits=(15, 130))
[both.parameter(path) for path in SHIFTS]
[instrument.zero_shift -0.166(38) vary,
instrument.geometry.sample_displacement -0.234(68) vary,
phases.0.cell.a 9.3690(10) vary]
Each correction refined to several esds from zero. Ask the diagnostics which pairs they found correlated with the zero shift.
[
d
for d in both.diagnostics
if d.code == "HIGH_CORRELATION" and "instrument.zero_shift" in d.where
]
[warning[HIGH_CORRELATION]: Correlated pair: instrument.zero_shift ~ instrument.geometry.sample_displacement (ρ=+1.000). Flagged in stages: zero_disp, cell, profile_w, profile, lines_axial, roughness, closing.
at: instrument.zero_shift, instrument.geometry.sample_displacement
help: Fix one of the two at a known value, or free them in separate stages.,
warning[HIGH_CORRELATION]: Correlated pair: phases.0.cell.a ~ instrument.zero_shift (ρ=+0.997). Flagged in stages: cell, profile_w, profile, lines_axial, roughness, closing.
at: phases.0.cell.a, instrument.zero_shift
help: Fix one of the two at a known value, or free them in separate stages.,
warning[HIGH_CORRELATION]: Correlated pair: phases.0.cell.c ~ instrument.zero_shift (ρ=+0.995). Flagged in stages: cell, profile_w, profile, lines_axial, roughness, closing.
at: phases.0.cell.c, instrument.zero_shift
help: Fix one of the two at a known value, or free them in separate stages.]
The first HIGH_CORRELATION warning says the zero shift and the displacement move together completely.
The fit could trade one against the other without changing the pattern, so neither value is a measurement.
The cell’s a and c are correlated with the zero shift too.
About correlations and esds lists the ways out.
Go back, and hold one of the pair#
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.
Every fit committed its stages to the refinement’s history.
checkout returns the working state to the tied fit’s last node, tie included.
set_values puts the zero shift at 0, and hold keeps it there whatever the plan frees.
ref.checkout(tied.node_id)
ref.set_values({"instrument.zero_shift": 0.0})
ref.hold(["instrument.zero_shift"])
['instrument.zero_shift']
hold returns the paths it now holds.
Run the same plan again.
held = ref.fit(data, plan="lab_bragg_brentano", two_theta_limits=(15, 130))
[held.parameter(path) for path in SHIFTS[1:]]
[instrument.geometry.sample_displacement 0.0630(15) vary,
phases.0.cell.a 9.373322(81) vary]
The list starts at the displacement, because a held parameter is not refined and the result has no row for it. The displacement’s esd and the cell’s are both far smaller than they were with the pair free.
The plan still tried to free the zero shift, and the result says so in HOLD_BLOCKED_PLAN.
[d for d in held.diagnostics if d.code == "HOLD_BLOCKED_PLAN"]
[info[HOLD_BLOCKED_PLAN]: 1 held parameter was matched by the turn_on glob of stages 'closing', 'zero_disp' and stayed fixed: instrument.zero_shift.
at: instrument.zero_shift
help: Call unhold() on them first if you meant the plan to refine them.]
Use ref.hold for this, and leave the parameter’s vary flag alone.
A plan sets which parameters vary at each stage, so it would free a parameter marked vary=False.
A hold outranks the plan.
How to keep a parameter fixed whatever the plan frees covers releasing one too.
The final structural fit#
lab_bragg_brentano refines the profile and the instrument but leaves the atoms alone.
Refine the structure once more, with the zero shift held, the displacement at the value the last fit measured, and the tie in place.
final = ref.fit(data, plan="mccusker_structural", two_theta_limits=(15, 130))
final.plot()
The difference curve looks as it did after notebook 02’s first fit. What changed is which numbers are measurements. The displacement now has a usable esd, and the zero shift is a value we stated.
The history#
Every stage of every fit is a node you can return to.
ref.history
tcb8d2285 56 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 stage:closing Rwp 0.0898
n0013 set_tie —
n0014 stage:scale_bkg Rwp 0.0895
n0015 stage:zero Rwp 0.0895
n0016 stage:cell Rwp 0.0895
n0017 stage:profile_w Rwp 0.0895
n0018 stage:profile Rwp 0.0895
n0019 stage:extra_components Rwp 0.0895
n0020 stage:coordinates Rwp 0.0895
n0021 stage:biso Rwp 0.0894
n0022 stage:preferred_orientation Rwp 0.0894
n0023 stage:extinction Rwp 0.0894
n0024 stage:roughness Rwp 0.0894
n0025 stage:closing Rwp 0.0898
├─ n0026 stage:scale_bkg Rwp 0.0894
│ n0027 stage:zero_disp Rwp 0.0894
│ n0028 stage:cell Rwp 0.0889
│ n0029 stage:profile_w Rwp 0.0889
│ n0030 stage:profile Rwp 0.0889
│ n0031 stage:lines_axial Rwp 0.0883
│ n0032 stage:roughness Rwp 0.0883
│ n0033 stage:closing Rwp 0.0887
└─ n0034 set_value —
n0035 set_hold —
n0036 stage:scale_bkg Rwp 0.3860
n0037 stage:zero_disp Rwp 0.0986
n0038 stage:cell Rwp 0.0897
n0039 stage:profile_w Rwp 0.0897
n0040 stage:profile Rwp 0.0897
n0041 stage:lines_axial Rwp 0.0890
n0042 stage:roughness Rwp 0.0890
n0043 stage:closing Rwp 0.0894
n0044 stage:scale_bkg Rwp 0.0890
n0045 stage:zero Rwp 0.0890
n0046 stage:cell Rwp 0.0890
n0047 stage:profile_w Rwp 0.0890
n0048 stage:profile Rwp 0.0890
n0049 stage:extra_components Rwp 0.0890
n0050 stage:coordinates Rwp 0.0890
n0051 stage:biso Rwp 0.0890
n0052 stage:preferred_orientation Rwp 0.0890
n0053 stage:extinction Rwp 0.0890
n0054 stage:roughness Rwp 0.0890
*n0055 stage:closing Rwp 0.0894
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.
The set_tie, set_value and set_hold nodes are the three edits.
* 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 (3 evals, ftol=1e-06)
stage profile_w: converged (4 evals, ftol=1e-06)
stage profile: converged (6 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 (4 evals, ftol=1e-06)
stage extinction: converged (3 evals, ftol=1e-06)
stage roughness: converged (1 evals, ftol=1e-06)
stage closing: converged (9 evals, ftol=solver default), max|Δθ|/esd=0.000
diagnostics under structure: 2 warnings, 0 errors, 3 info
warning[FLOOR_HIT]: phases.0.extinction refined to its floor (ρ=+0.064).
at: phases.0.extinction
help: Quote an upper limit instead of a value, or remove the parameter if its phase or term is unsupported. For a displacement parameter, check absorption, species and occupancy.
warning[PATTERN_UNDERSAMPLED]: The median fitted peak spans 4.6 steps across its FWHM (76 peaks), below the 5 the guidelines ask for.
help: Re-collect at a step near FWHM/5 if the intensities must be quotable.
info: RESOLUTION_UNCONSTRAINED at instrument.profile.u, .v, .w (warning under microstructure, calibration)
info: CONSTRAINT_ACTIVE
info: HOLD_BLOCKED_PLAN at instrument.zero_shift
report:
62 misfit regions. The top 15 carry 47% of χ² and the other 47 carry 30%.
5 observed peak(s) match no reflection.
Layer 1 abstained (unreadable). Misfitting regions carry 45% of χ², but only 20% of that sits in regions the local gates accept. Layer 1 needs 40%.
region 31.5-32.5°: 9% of χ², max|Δ|/σ=12.5
region 103.2-111.3°: 4% of χ², max|Δ|/σ=5.0
region 32.6-33.4°: 4% of χ², max|Δ|/σ=9.5
serial correlation (DW): 1.01, esds inflated ×2.20
next: free phases.0.atoms.2.occ, predicted ΔBIC +30.4 at N_eff 1191 (Δχ² 452.5)
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, backend=numpy, solver=lm
Rwp 0.0894 / Rexp 0.0559 (GoF 1.60), Rp 0.0695, χ² 2.6, DW 1.01
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.
reportreads the difference curve region by region. It counts the misfit regions and the observed peaks with no reflection under them, and lists the regions carrying most of χ². Find those peaks before quoting a structure, because a weak second phase would explain them.When too little of the misfit sits in regions it can explain, Layer 1 abstains from naming a cause. An abstention is an answer: no single correction is indicated.
nextnames the parameter whose release rietx predicts would help most, and the ΔBIC that release would have to pass.exchangeabilitylists held parameters whose effect a free one could imitate. Fitting each rival alone, with the other held, settles which one the data choose.protocolis what you report alongside the numbers: the plan, what was held, and the observation count.
About the report explains each layer, and How to choose the next parameter to free acts on next.