Joint refinement#
This page gives MultiHistogramRefinement and refine_multi, which refine one structure against several patterns at once: their arguments, which parameters the patterns share, two radiations in one fit, a refinable wavelength, the joint result and its per-pattern half, the events and cancellation of a joint fit, and a report per pattern.
The patterns are stacked into one residual [Von Dreele, 1997].
GSAS-II calls this a combined refinement, and calls each pattern a histogram, as rietx does.
Shared structural parameters draw information from every pattern, while each pattern keeps its own Instrument for its wavelength, geometry, resolution and background.
A series of patterns of a changing specimen is a chain instead (Sequential refinement), and a parameter that follows an equation of each pattern’s coordinates is declared on a joint fit (Parametric refinement).
How to fit an X-ray and a neutron pattern together runs a joint fit from start to finish.
Entry points#
import rietx as rx
joint = rx.MultiHistogramRefinement(structure, [instrument_a, instrument_b])
result = joint.fit([pattern_a, pattern_b], plan="mccusker_default")
result = rx.refine_multi([pattern_a, pattern_b], structure, [instrument_a, instrument_b])
refine_multi is the same run as one call, with the patterns first.
Both return an ordinary RefinementResult, with the per-pattern half on RefinementResult.histograms (Joint result).
MultiHistogramRefinement.fit takes the patterns in the order of the instruments given to the constructor, and the two lists must be the same length.
Argument |
Taken by |
Default |
Does |
|---|---|---|---|
|
constructor, |
required |
one |
|
constructor, |
required |
one |
|
constructor, |
|
which parameters the patterns share (Shared parameters) |
|
constructor, |
|
a parameter declared as an equation of each pattern’s coordinates (Parametric refinement) |
|
constructor, |
|
as for a single fit (Entry points) |
|
|
|
a preset name or a |
|
|
|
any other mode raises |
|
|
|
one |
|
|
|
one positive relative weight per pattern’s residual block |
|
|
|
as on |
Attribute |
Holds |
|---|---|
|
how many patterns there are |
|
the multi-histogram parameter table: one ordinary table per pattern, threaded by a column map that folds the shared columns onto one |
|
the parametric declarations as resolved at construction (Reference and scale) |
|
after a fit, each pattern’s refined models; the shared parameters are identical across the structures |
|
after a fit, the result |
|
the |
A joint fit runs in Rietveld mode only. Le Bail and Pawley intensities are per-pattern extractions rather than shared quantities, so a joint fit of them would be independent single fits.
Two radiations#
A joint fit may mix radiations.
An X-ray pattern and a constant-wavelength neutron pattern of one specimen refine one structure, the combined refinement [Von Dreele, 1997] was written for.
There is nothing extra to declare: each Instrument carries its own source, and every correction that depends on the radiation reads the source of the pattern it is computing.
import rietx as rx
xray = rx.Instrument.debye_scherrer(wavelength=0.4132950)
neutron = rx.Instrument.constant_wavelength_neutron(1.54040, fwhm_deg=0.30)
assert (xray.source.kind, neutron.source.kind) == ("xray_cw", "neutron_cw")
assert xray.source.dispersion is not None and neutron.source.dispersion is None
The split between shared and per-pattern parameters has nothing radiation-specific in it. The specimen’s parameters are shared: the cell, the coordinates, the occupancies, the displacement parameters and a magnetic moment. Each pattern keeps its own scale, background, zero shift, peak width and wavelength.
What the second radiation adds is a different weighting of the same structure factor.
X-rays scatter from electrons, so the heavier atoms dominate the form factor f, while a nucleus’s scattering length b follows no such order.
Corundum ranks its two sites oppositely: f(0) is 13 electrons for Al and 8 for O, while b is 3.449 fm for Al and 5.803 fm for O (Sears’s table, which b_Sears.dat holds).
On synthetic corundum patterns with a known answer, the X-ray pattern alone determined z(Al) about 3 times better than the neutron pattern alone, and x(O) about as well.
Refined jointly, the esd of z(Al) stayed close to the X-ray pattern’s own, and the esd of x(O) fell to about two thirds of it: 0.62-0.68 across the builds tested.
The neutron pattern added the oxygen and left the aluminium to the X-rays.
In the same fits the two scales, 3.5 times apart, and two zero shifts of opposite sign each came back on their own pattern’s value.
Correction |
On an X-ray pattern ( |
On a neutron pattern ( |
|---|---|---|
Scattering amplitude |
the form factor f₀ of the species, falling with s = sin θ/λ, plus f′ + i·f″ |
b of the species, one number at every s |
Anomalous dispersion |
|
none: |
Polarisation, |
refinable |
pinned at K = 1 and locked |
Magnetic structure factor |
never built |
built for a phase that declares a moment |
µR estimated from composition |
when the geometry declares a capillary radius, from X-ray attenuation |
when the geometry declares a capillary radius, from neutron cross-sections (Sears 1992) |
A magnetic moment is shared like any other structural parameter, so the X-ray pattern carries its column and adds no gradient to it.
The magnetic term is never built there: rietx.model.forward.magnetic_wanted decides, from rietx.crystallography.magnetic.scattering.MAGNETIC_SOURCE_KINDS.
The moment is measured by the neutron pattern alone, against a structure both patterns constrain.
Three shared size coefficients, lor_size, gauss_size and magnetic_lor_size, are a different number per pattern.
That is a matter of wavelength rather than radiation (Size in a joint fit, rietx.params.multi.SIZE_LAMBDA_POWER).
The three diagnostics that key on the radiation, DISPERSION_NEGLECTED, NEUTRON_RESONANT_ABSORBER and SPECIES_FALLBACK_NEUTRAL, land on the pattern whose radiation raised them, never at the top level (Diagnostic placement).
HistogramResult carries no instrument, so the result does not say which radiation each pattern was.
Read MultiHistogramRefinement.fitted_instruments instead, whose entry for each pattern carries its source.kind.
The real-data joint fit this package is tested against is one Nd₂Ru₂O₇ specimen measured on the APS 11-BM synchrotron and the NCNR BT-1 constant-wavelength neutron diffractometer, with the neutron wavelength refined.
Wavelength refinement#
In a single fit a wavelength and the cell trade exactly. Bragg’s law (2.4) gives the pattern access to λ/(2 sin θ), so only the product of λ with the reciprocal cell is measurable. A free λ beside a free cell is therefore refused, and a free λ beside a held cell, a certified standard’s, is how a wavelength is calibrated (Wavelength refinement).
Several patterns of one specimen share one cell, and that breaks the degeneracy with the cell free. Holding one wavelength pins the cell’s scale. The remaining N − 1 are then over-determined by the shared cell and measurable, each against the same lattice. Holding all of them forces every monochromator’s calibration error into the shared cell, and freeing all of them puts the flat direction back. So a free wavelength requires its pattern’s cell to be shared with at least one pattern whose wavelength is held.
“Exactly one held, at most N − 1 free” is that rule when every pattern shares one cell, which is the SharingMap default.
Un-share the cell and the single-pattern degeneracy is back, per pattern, inside a joint fit where it looks solved, so that combination is refused too.
Each refusal names its cause and states the count: “2 of 2 wavelengths are free; hold one”.
Hold the wavelength of the pattern that determines the cell. On a synchrotron-plus-neutron pair that is the synchrotron. Its wavelength calibration is the better known, and its angular resolution makes its cell the better determined. What the neutron data supplies is the structural content, meaning oxygen positions, site occupancies and displacement parameters, through scattering-length contrast X-rays do not have. The neutron wavelengths then refine against a cell the X-ray pattern has pinned, which is the only way their monochromator calibrations can be measured with the cell free.
A wavelength is per-pattern under the default sharing rule, so a scoped glob frees it. The bare path would free every pattern’s copy and be refused:
import rietx as rx
from rietx.strategy.staged import PLAN_PRESETS, Stage
plan = PLAN_PRESETS["mccusker_structural"]()
# pattern 0 is the synchrotron and keeps its declared λ; pattern 1's neutron λ
# refines against the cell pattern 0 has pinned. The cell and the coordinates
# refine with it, because λ trades against both.
plan.stages.append(Stage("wavelength",
["hist.1.instrument.source.lines.0.wavelength",
"phases.*.cell.*", "phases.*.atoms.*.dof.*"]))
result = rx.refine_multi([xray, neutron], structure, [ins_x, ins_n], plan=plan)
The number to read is the WAVELENGTH_CALIBRATION diagnostic on that pattern.
It reports how far λ moved from its declared value in ppm, the unit a calibration error is quoted in, and how that compares with its own esd.
A refined λ inside its own esd measured nothing, and the diagnostic says so.
Only line 0’s wavelength refines, because the ratio between lines of one source is atomic physics.
A known ratio such as a λ/2 harmonic is a Refinement.tie with scale=0.5 in a single fit, and a neutron Harmonic follows its fundamental.
The reflection list, evaluation windows and quadrature node counts are frozen from the declared λ at each stage compile.
So a free λ moves peaks inside their frozen windows, as a free cell does.
That holds while the motion is small against the window, and a few hundred ppm is: 250 ppm of λ moves a peak at 2θ = 150° by 0.11°, against a 0.30° FWHM.
All of this is constant wavelength. The same rule carries to a time-of-flight multi-bank fit, where each bank carries its own DIFC calibration and one of them has to be held to pin the cell. Nothing here implements it.
Profile scan#
A joint fit often has one parameter that is weakly determined for a physical reason: an antisite fraction between two elements whose scattering lengths barely differ, a site occupancy, a small strain. Freeing it with everything else has two failure modes. The solver can land in a local minimum, and its esd comes from a linearised curvature the data does not have.
A profile scan holds the parameter on a grid, refines everything else at each point, and reads the value off the minimum of χ² and the interval off its width.
It measures the curvature from which rietx.optimize.statistics.normal_covariance builds an esd, and ParameterTable.unmeasured_rows names what a direction the data does not move reached.
There is no scan verb.
On a single Refinement the scan composes verbs that exist:
Step |
Verb |
|---|---|
fork a fresh line of work per grid point |
|
pin the scanned parameter at the grid value |
|
refine the rest |
|
read the curve |
|
set_values and hold commit history nodes, so the whole scan is recoverable.
Branching per grid point keeps one point from seeding the next.
MultiHistogramRefinement keeps no history and has no set_values or hold verb.
Scan a joint fit with one fit per grid point, each from a structure carrying the grid value, under a plan whose stages do not free the scanned path, and read each result’s RefinementResult.statistics.
Set the value on the shared model.
SharingMap decides whether the scanned path is shared, and scanning a per-pattern path pins it in one pattern only, a different experiment.
The minimum gives the value, and the interval where χ² rises by 1 above the minimum is the 1σ range for one parameter. A scan measures only the surface it searched. A coarse grid can step over a narrow minimum, and the interval is only as good as the grid spacing, so refine the grid near the minimum before quoting a number.
Joint result#
RefinementResult.histograms is a list of HistogramResult, one per pattern.
An empty list means an ordinary single-pattern fit.
Field |
Holds |
|---|---|
|
that pattern’s name |
|
the relative weight applied to its residual block |
|
its own agreement indices |
|
its curves |
|
its per-point σ |
|
its reflection positions, by phase |
|
which reflection each of those is, paired by index |
|
its phase quantities |
|
always |
|
its RBragg and RF, per phase |
|
what that pattern reported |
A pooled Rwp is never quoted alone. A single pooled number can hide a badly fitting pattern, so each pattern reports its own. On two synthetic LaB6 patterns of one crystal, at λ = 0.41390 Å (4200 points) and λ = 0.71070 Å (8000 points), the pooled Rwp was 0.0514 while the two patterns were at 0.0411 and 0.0613. The pooled figure describes neither.
HistogramResult.weight is 1.0 at unit weight, where each point’s own esd governs.
Provenance.notes records the weighting under histogram_weights, as unit (each point's esd governs) or as each pattern’s weight, so a non-unit weighting is never silent.
RefinementResult.data_support is None for a joint fit (Diagnostic placement).
RefinementResult.for_histogram returns a single-pattern view of one pattern.
Its curves, statistics, ticks, phase quantities and diagnostics move to the top level and histograms is cleared, so result.for_histogram(0).plot() and a report built from it work per pattern.
The view’s parameters keep the shared rows and that pattern’s own rows.
Its own rows lose their hist.<h>. prefix, so they are spelled as a single fit spells them.
Every other pattern’s rows are dropped.
The view carries no tick_hkl and no phase_agreement, so read those on the HistogramResult.
RefinementResult.statistics carries the joint solve’s Statistics.max_shift_over_esd, one shift per column of the joint parameter vector against its esd, and every view carries the same figure.
On the two LaB6 patterns a converged fit reads between 5 × 10⁻⁷ and 2.2 × 10⁻³ across both solvers, and a stage starved to one iteration reads 4.0 under LM and 17.9 under TRF.
Diagnostic placement#
A joint fit reports each diagnostic in one of two places, and the place follows from what the diagnostic is about.
A finding about one pattern goes on that pattern’s HistogramResult.diagnostics, once per pattern: its radiation, its instrument, its counts, or how well the model fits it.
A finding about the shared structure, the joint solve or the plan goes once on RefinementResult.diagnostics.
So a mixed X-ray and neutron fit whose X-ray source declines anomalous dispersion raises DISPERSION_NEGLECTED on the X-ray pattern only.
The same holds for SPECIES_FALLBACK_NEUTRAL, and NEUTRON_RESONANT_ABSORBER goes on the neutron pattern only.
STAGE_MAX_ITER and SPACE_GROUP_SETTING_ASSUMED appear at the top level, once each.
BACKGROUND_ABSORPTION is the one finding on both lists.
It is raised on its pattern’s list and repeated once on RefinementResult.diagnostics, where it keeps its hist.<h>. path.
Each copy is the same finding about one pattern, so count it once.
The rule is written as data.
rietx.multi.DIAGNOSTIC_SCOPES places the output of each function that writes diagnostics, and rietx.multi.GUARD_SCOPES each field of a GuardReport, each row with its reason.
Six rows are marked absent, for diagnostics a joint fit does not compute: the two Pawley ones (a joint fit is Rietveld-only), HOLD_BLOCKED_PLAN (a joint fit has no hold verb), RESTRAINT_TENSION (restraints are refused), and the two magnetic-width checks, STAGE_FREES_MAGNETIC_WIDTH_WITH_MOMENT and MAGNETIC_WIDTH_MOVED_MOMENT.
The magnetic-width checks read the single-pattern plan, so a joint fit’s width-before-moment ordering is unchecked.
DATA_SUPPORT_LOW is withheld on purpose.
Its ratio of observations to parameters has no agreed definition across patterns, and one pattern’s reflections counted against the shared columns would warn about a fit the other patterns support.
Its two sibling checks, PATTERN_UNDERSAMPLED and PATTERN_DEAD_CHANNELS, are about one measurement, so each pattern gets them.
Events and cancellation#
MultiHistogramRefinement.fit takes events=, cancel=, telemetry=, label= and progress=, and each means what it means on Refinement.fit (Events, Progress lines).
The fit records itself into one run directory however many patterns it holds.
joint = rx.MultiHistogramRefinement(structure, [instrument_a, instrument_b])
result = joint.fit([pattern_a, pattern_b], events="joint.jsonl")
The joint solve is one solve, so each event kind is emitted once per fit, stage or evaluation.
On the two LaB6 patterns the log held one fit_start, five stage_start, 40 eval, five stage_end and one fit_end.
An eval measures the stacked residual.
Its values line up with the free_paths of the latest stage_start, where a row that belongs to one pattern is spelled hist.<h>. and a shared row is bare.
A fact about one pattern rides in a histograms list on fit_start, stage_start, stage_end and fit_end, one row per pattern, each stamped with its histogram_index.
Every one of those events also carries n_histograms.
A fit_start row holds n_points and n_fitted, and a stage_start row holds n_points.
A stage_end row holds that pattern’s Rwp at the end of the stage, beside the pooled rwp: on the last stage of that fit, 0.0411 and 0.0613 against 0.0514, the numbers the result reports.
A stage_start row holds unreached, the globs that matched rows of another pattern and none of this one.
A cancelled stage is abandoned.
The parameter table and every pattern’s models go back to where the stage found them, and the stage still ends on a stage_end with status="cancelled".
RefinementCancelled.completed_stages names the stages that finished, and its node_id is None, because a joint fit keeps no history (Cancellation).
A joint run writes no snapshot.json, so rietx watch lists it and follows its status with no picture.
The live view draws one pattern, and a joint fit has several.
MultiHistogramRefinement.fit takes no stage_reports=.
A report at each stage would be one report per pattern per stage, and a report build takes about 0.36 s however fast the fit is.
Until that changes, a joint fit’s stage trajectory is its event log, with each pattern’s Rwp on every stage_end.
Report per pattern#
MultiHistogramRefinement.report builds the FitReport for one pattern of the last fit, at Layers 0 to 2 (The fit report).
Its other keywords go to build_report.
It raises RuntimeError before a fit, and IndexError for a pattern that does not exist.
report = joint.report(1, plan="mccusker_default")
Each report reads that pattern’s curves and statistics through RefinementResult.for_histogram, and that pattern’s compiled model for Layers 1 and 2.
The background section reads that pattern’s own absorption screen, which the joint result keeps in Identifiability.background_absorption under hist.<h>. keys.
The joint path runs no correlation, soft-mode or extra-peak screen and no screen for exchangeable parameters, so those tables are empty.
The specimen has one crystallite size, so every pattern’s report carries the joint result’s microstructure block, read off the shared column at pattern 0’s wavelength.
On a two-wavelength pair in tests/test_joint_report.py both reports quote 408.82(4.86) Å, and the two patterns refined alone give 408.71(8.05) Å and 409.20(6.43) Å.
The separability caveat beside the size is each report’s own, from that pattern’s 2θ range.
A joint result’s curves come from the last stage’s compile, sized at the values the stage started from. A single fit compiles the model again at its answer, and a joint fit does not. If the last stage broadened the peaks, their tails stop at the old windows and Layer 0 reports the cut as unmatched peaks. On that size pair the reports find 44 and 27 unmatched peaks where a compile at the answer finds none, so read Layer 0 there against the plot.
Precision of a joint fit#
On the two LaB6 patterns, the shared cell came back identical in both entries of MultiHistogramRefinement.fitted_structures: a = 4.156597 Å, against the 4.15660 Å the patterns were built from, 0.6 ppm out.
The per-pattern zero shifts separated, at 0.005997° and −0.010008° against the 0.006° and −0.010° that went in.
The cell’s esd was 2.43 × 10⁻⁶ Å, against 4.70 × 10⁻⁶ Å and 2.83 × 10⁻⁶ Å from the two patterns refined alone: 1.94 and 1.17 times smaller than either.
Those numbers come from synthetic patterns with known answers. Real data has no truth to check against, so watch the per-pattern Rwp spread instead.