Technical Overview & Methodology

This document is a technical reference for the Brutus powder indexing software. It explains the main algorithms, search parameters, and methodology, and is intended for users familiar with powder X-ray diffraction and crystallography. As the name suggests, this is a brute-force method, the search is exhaustive rather than heuristic, and is not necessarily fast.

Core Goal

The aim of ab initio powder indexing is to determine the unit-cell parameters ($a, b, c, \alpha, \beta, \gamma$) from a list of observed diffraction peak positions ($2\theta$). Brutus performs this task using a symmetry-specific, (almost) exhaustive search algorithm.

The central assumption is that a small subset of the most intense, low-angle reflections corresponds to simple crystal planes with low-integer Miller indices $(hkl)$. For a given crystal system, the program chooses some observed peaks and solves the resulting system of linear equations.

Q-Space Formulation

All peak positions are first converted from $2\theta$ to Q-space, where $Q = 1/d^2$. The general quadratic relationship between $Q$, the Miller indices, and the reciprocal cell parameters ($A, B, C, D, E, F$) is:

$$ Q_{hkl} = Ah^2 + Bk^2 + Cl^2 + Dkl + Ehl + Fhk $$

where $A = a^{*2}$, $B = b^{*2}$, $C = c^{*2}$, and the cross-terms involve the reciprocal angles. Working in Q-space linearises the problem: for a given assignment of Miller indices to observed peaks, the unknown parameters $\{A, B, C, \ldots\}$ appear as the solution to a system of linear equations.

Brutus solves for these reciprocal parameters (or the subset relevant to the current crystal system), then converts them back to real-space cell parameters. Each candidate cell is then refined with a weighted least-squares fit and scored against the full observed peak list via figures of merit.

System-Specific Parameterisation

The quadratic form simplifies for each crystal system, reducing the number of free parameters. The design matrix column vector for the least-squares system is:

SystemFree parameters ($k$)Design row getLSDesignRow(hkl)
Cubic1$[h^2+k^2+l^2]$
Tetragonal2$[h^2+k^2,\; l^2]$
Hexagonal2$[\tfrac{4}{3}(h^2+hk+k^2),\; l^2]$
Orthorhombic3$[h^2,\; k^2,\; l^2]$
Monoclinic4$[h^2,\; k^2,\; l^2,\; hl]$
Triclinic6$[h^2,\; k^2,\; l^2,\; kl,\; hl,\; hk]$

The solved parameters relate to reciprocal lattice constants. For example, in orthorhombic, the three fitted parameters are $A=1/a^2$, $B=1/b^2$, $C=1/c^2$, giving $a=1/\sqrt{A}$ etc. For monoclinic, the four parameters are $A=1/(a^2\sin^2\beta)$, $B=1/b^2$, $C=1/(c^2\sin^2\beta)$, $D=-2\cos\beta/(ac\sin^2\beta)$, from which $\beta$ is extracted as $\cos\beta = -D/(2\sqrt{AC})$. For the CPU indexing, more than one system can be selected at a time, the work will be sequantially done for each system. In the case of the GPU, only one system can be selected at a time, this is mainly because each system has specific parameters to set (hkl pools, etc).

Architecture

Brutus runs entirely in the browser, requiring no installation, and it can be completely run locally in a server. The program is split across several files:

Quick Start Guide

Use the following workflow for a typical single-phase powder pattern.

  1. Load the data file. Click Load File. Supported formats include .xy, .xrdml, .xml (Panalytical XRDML or Bruker RawDataFile), .ras, .xra, .uxd, .udf, GSAS .esd, and plain text two-column data.
  2. Detect peaks. On the Peaks tab, adjust the Min peak (%), Radius (pts), and Points sliders until the automatically detected peaks match the visual pattern.
  3. Curate the peak list. Carefully review all peaks:
    • Edit $2\theta$ positions for accuracy.
    • Delete spurious peaks (noise, impurities, Kα2 if stripping is off).
    • Add important missing peaks, e.g. low-angle reflections, using Ctrl + Click on the chart.
    A clean list of about 15–20 peaks, free of impurities at low angle, is ideal. The result of the indexing will depend only on the quality of the data.
  4. Set parameters. On the Parameters tab:
    • Select the correct X-ray Radiation Preset (e.g. Cu Kα).
    • Choose whether to enable Strip K-alpha2. When enabled, the Ka1 wavelength is used; when disabled, the average Kα wavelength is used. The default is OFF (average Kα wavelength).
    • Set a chemically reasonable Max Volume (ų) to limit the search space.
    • Set 2θ Error (°) according to your data quality (e.g. ≈0.02° for synchrotron, ≈0.05° for a typical lab diffractometer).
    • Leave Refine Zero-Point Error enabled unless you have a specific reason to turn it off.
    • Select the crystal systems to search. Enabling Orthorhombic, Monoclinic, or Triclinic activates GPU-accelerated searches. If these systems are disabled, it is likely that the GPU acceleration is not activated in the browser. See the Troubleshooting section.
  5. (Optional) Tune GPU Parameters. If searching Orthorhombic, Monoclinic, or Triclinic, a new field GPU Search Parameters appears.
    • HKL Basis Size: Number of simple HKLs to use for the search (Defaults: Ortho 300, Mono 100, Tri 40).
    • Peaks to Combine: Number of observed peaks used to form combinations (Defaults: Ortho 7, Mono 7, Tri 9).
    • Leave these at their defaults unless you have trouble finding a solution.
  6. Start indexing. Click Start Indexing. Progress is shown on the main progress bar. If the number of GPU buffers is too large, a warning will appear.
  7. Inspect solutions. On the Solutions tab:
    • Sort solutions by M(20).
    • Click a row to display calculated (blue) and observed (red) tick marks on the chart.
    • A plausible solution will show excellent alignment and reasonable space-group suggestions.
    • A right click on a solution will open a menu with several additional functions, see below.

The User Interface

The application window is divided into a Controls Panel (left) and a Results Area (right), separated by a draggable divider that allows you to resize both panels.

Controls Panel

The controls are organized into three main tabs.

1. Peaks Tab

2. Parameters Tab

GPU Search Parameters

When Orthorhombic, Monoclinic, or Triclinic is selected, a dedicated panel appears with four parameters controlling the scope, speed, and memory usage of the GPU search.

Note: The status text at the bottom of the panel updates in real-time as you adjust these values, showing the total estimated number of cells to be tested. Ensure this number remains reasonable (e.g. under 500 Billion for a result within a few minutes). For a regular GPU card (e.g. Nvidia T1000, 8 GB) the effective search speed is typically in the range of tens of millions of cells per second.

3. Solutions Tab

Results Area

Saving & Converting Data

Brutus can export the pattern it is displaying, which makes it a convenient format converter as well as an indexer. The Save as button sits next to Load File and becomes active once a data file has been loaded. Clicking it opens a small dialog where you choose an output format; pressing Save then hands off to the browser's own save dialog, so you pick the folder and file name in the normal way.

What is exported

The exported pattern is exactly what is plotted on the chart, so two view choices carry through to the file:

The abscissa is always written in 2θ regardless of whether the chart is currently showing θ, d or Q — those are display transforms only, and the underlying pattern is stored in 2θ.

Formats

The plain-text column formats store an explicit 2θ value for every point, so they preserve a non-uniform scan faithfully. The instrument formats store a single start position, a step and a list of intensities, so they assume a constant step; if the loaded scan is not evenly spaced the dialog warns that positions may shift slightly when written to one of these.

Format Extension Stores λ Grid
XY (2θ intensity).xynoper-point 2θ
CSV (2θ,intensity).csvnoper-point 2θ
Two-column whitespace.datnoper-point 2θ
Panalytical XRDML.xrdmlyesconstant step
Bruker RawDataFile XML.xmlyesconstant step
Bruker UXD.uxdyesconstant step
GSAS ESD.esdnoconstant step

For the formats that record a wavelength, the value written is the one in the Ka1 Wavelength box, and the dialog shows it before you save. The default file name is taken from the loaded file; change it in the browser's save dialog if you want.

A note on the informational messages. The lines in the dialog are labelled Info, not warnings — they describe what is being written (Kα2 state, visible range, wavelength, constant-step assumption) so you can confirm the output matches your intent. None of them blocks the save.

Peak Finding in Detail

Accurate peak positions are the single most important input for successful indexing. Brutus uses a multi-step procedure to detect peaks from raw intensity data.

Algorithm Steps

  1. 2 stripping (optional): If Strip K-alpha2 is enabled, the Zhang (2026) dual-wavelength $q$-space differencing method is applied to the raw intensities. By operating in reciprocal space, it cleanly decouples overlapping doublet reflections without needing prior background subtraction or empirical peak-shape parameters.
  2. Background subtraction: A rolling-ball style algorithm estimates and removes the background. The Radius slider controls the ball radius in data points. A larger radius smooths over broader background humps.
  3. Data smoothing: A Savitzky–Golay filter is applied to the background-subtracted signal to reduce noise while preserving peak shapes. The Points slider sets the window width.
  4. Initial peak detection: Local maxima above the Min peak (%) threshold are identified.
  5. Position refinement: For each detected peak, a five-point least-squares quadratic fit (Savitzky–Golay coefficients) is carried out around the local maximum to obtain a sub-step precision position. If the peak is too close to the data edge, the algorithm falls back to a three-point fit.

Suspected Kα2 Peaks & Wavelength Treatment

If your data has not been stripped of Kα2, Brutus automatically scans the peak list to flag suspected Kα2 ghosts. A candidate must satisfy two independent tests:

When no Kα2 exists, none of this applies. The entire Kα2 tagger is inactive for radiation with no doublet component — that is, Custom wavelengths, any 1-only preset, and data for which Kα2 stripping is enabled. In those cases no peak is ever marked Kα2-suspect, so no reflection is withheld from the indexing search or from the M(20) calculation, and no extinction violation is demoted to "soft" on Kα2 grounds. Switching to such a radiation also clears any flags raised under a previous doublet preset, so nothing stale survives the change.

Flagged peaks are then treated as follows:

Practical Recommendations

Indexing Algorithm and Search Strategy

Brutus uses an exhaustive, symmetry-specific trial-and-refine indexing algorithm. For each crystal system, the program generates trial solutions from combinations of low-angle observed peaks and low-index Miller indices, solves the corresponding linear system in reciprocal space, rejects unphysical cells, and then performs a full weighted least-squares refinement on surviving candidates.

Linear System Formulation

The search is formulated as a system of linear equations based on the quadratic form. Given a trial assignment of Miller indices $\{(hkl)_i\}$ to observed peaks $\{Q_{obs,i}\}$:

$$ Q_{obs,i} = \sum_j P_j \cdot H_{j,i} $$

where $H_{j,i}$ is the $j$-th element of the design row for $(hkl)_i$ (see the table in the Overview), and $P_j$ are the unknown reciprocal-lattice parameters. Stacking $k$ such equations gives a $k \times k$ square system that is solved exactly (or in the over-determined case, by weighted least-squares).

CPU Searches (High Symmetry) CPU

High-symmetry systems have fewer degrees of freedom ($k \leq 2$) and are solved on the CPU in a background Web Worker. Recent updates have significantly expanded the search scope to ensure large unit cells are not missed.

GPU-Accelerated Searches (Low Symmetry) GPU

For Orthorhombic, Monoclinic, and Triclinic systems, the combinatorial explosion makes CPU searching impractical. Brutus offloads these tasks to the WebGPU API, allowing billions of trial cells to be evaluated in minutes.

The Combinatorial Search Space

The total number of cells tested depends on two user-configurable parameters.

Total Trials = Peak Combinations × HKL Combinations × Permutations $$ N_{total} = C(N_{peaks},\, k) \;\times\; C(N_{hkl},\, k) \;\times\; k! $$

System Specifics

System$k$Default $N_{hkl}$Default $N_{peaks}$Estimated trials (defaults)
Orthorhombic33007 $C(7,3) \times C(300,3) \times 6 \approx 936\;\text{Million}$
Monoclinic41007 $C(7,4) \times C(100,4) \times 24 \approx 3.3\;\text{Billion}$
Triclinic6409 $C(9,6) \times C(40,6) \times 720 \approx 232\;\text{Billion}$

These numbers can grow very rapidly. For example, increasing the Monoclinic $N_{hkl}$ from 100 to 150 multiplies the HKL combinations by roughly $(150/100)^4 \approx 5\times$, and increasing $N_{peaks}$ from 7 to 9 multiplies the peak combinations by $C(9,4)/C(7,4) = 126/35 \approx 3.6\times$.

Axial Front-Loading (finding "strange" axis ratios)

The HKL basis is sorted by index magnitude $h^2+k^2+l^2$, but that magnitude is only a proxy for a reflection's actual $Q$ when the cell is roughly isotropic. For a cell with one very long (or very short) axis, the reflection that indexes a given low-angle peak can be a high-order axial reflection — e.g. a long $c$ axis makes $(0\,0\,l)$ with large $l$ appear at low angle. Under a pure magnitude sort such a reflection sits far down the list and is discarded when the basis is truncated to $N_{hkl}$, so the cell cannot be found.

To prevent this, for orthorhombic and monoclinic Brutus partitions the sorted list into axial reflections (two of $h,k,l$ equal to zero) and everything else, and places all axials at the front before truncation. This guarantees that the axis-defining reflections survive regardless of $N_{hkl}$, giving the brute-force solver a chance to assign a high-order axial to a low-angle peak and so recover an unusual axis ratio. Because $N_{hkl}$ (and therefore $C(N_{hkl},k)$) is fixed, front-loading changes only which reflections are kept — it swaps the highest-magnitude tail reflections for the axials — and does not increase the search cost.

The depth of the axial set matters. Orthorhombic generates reflections to index 12, so its axials already reach $(0\,0\,12)$. Monoclinic's general grid only reaches index 6, which was too shallow for genuinely long axes (a $c \approx 30\;\text{Å}$ axis needs $(0\,0\,l)$ well beyond $l=6$ to index its observable peaks), so Brutus additionally injects deep axial reflections $(h00)/(0k0)/(00l)$ out to index 12 for monoclinic. These are cheap (three reflections per order) and are front-loaded like the rest; the only overhead is that an all-axial 4-tuple is rank-deficient in the $D=hl$ (i.e. $\beta$) column and is discarded by the near-zero-determinant filter — roughly $1.5\%$ of combinations at $N_{hkl}=100$.

Triclinic is deliberately excluded from front-loading. With $k=6$, the reciprocal angles ($D=kl$, $E=hl$, $F=hk$) are carried only by mixed reflections, and the basis is small (default $N_{hkl}=40$). Promoting axials there would crowd out the mixed reflections needed to determine the angles and drive most trials singular, so triclinic keeps the plain magnitude ordering. The practical consequence is that triclinic cells with one very long axis are intrinsically harder for this tool; raising $N_{hkl}$ (up to its ceiling of 123) is the available lever.

GPU Memory Management & Safe Batching

Testing hundreds of billions of cells generates massive amounts of data. To prevent crashing your browser or freezing your computer, Brutus uses a Safe Batching architecture:

  1. Dispatch chunking (thread budget): The HKL-combination dimension is split into chunks so that no single GPU dispatch launches more than a fixed per-system thread budget: 500,000 trials/dispatch for Orthorhombic and Monoclinic, and a tighter 50,000 for Triclinic (whose 6×6 solve and 720-permutation inner loop make each thread far heavier). The engine computes hklsPerChunk = safeWorkgroupsY × workgroupSizeY, iterates the full search as ceil(totalHklCombos / hklsPerChunk) sequential dispatches, and yields to the event loop between chunks so the browser tab stays responsive. The yield is a MessageChannel post rather than setTimeout(0): browsers clamp nested timers to a ~4 ms floor, so a per-chunk setTimeout(0) in a long run would waste several seconds sleeping in the clamp. The MessageChannel task is unclamped, and the engine only yields when the last yield was more than one display frame (~16 ms) ago, so keeping the tab responsive costs almost nothing.
  2. WebGPU dimension ceiling: A compute dispatch may request at most 65,535 workgroups per dimension; Brutus additionally clamps the Y (HKL) dimension to 16,383 workgroups per dispatch as a conservative cross-driver limit. This bound, combined with the thread budget above, is what fixes the chunk size and it allows indexing on old GPU cards.
  3. Time Limit (TDR Prevention): On Windows, if a single GPU dispatch takes longer than ~2 seconds, the OS assumes the driver has frozen and resets it (DXGI_ERROR_DEVICE_HUNG). The smaller Triclinic thread budget exists specifically to keep each dispatch well under this watchdog timeout.
  4. Buffer sizing: At device-request time Brutus asks the adapter for its maximum maxBufferSize and maxStorageBufferBindingSize. The results buffer is then sized to Buffer/Candidate (kCells) × cell-struct-size (32 bytes/cell for Triclinic, 16 for the others). All per-run GPU buffers are explicitly destroyed after each run (in a finally block) to avoid VRAM leaks across repeated searches. The GPU device itself is created once and reused by every search for the lifetime of the page, rather than re-created per run, and is released on page unload. If the browser or driver loses the device mid-search (a driver reset or TDR, or the tab being suspended), Brutus reports it and disables GPU searches until the page is reloaded, since a lost device cannot be revived. Note that an inactive panel may be throttled by the OS, so for full speed you may need to keep the panel active.
  5. Incremental read-back: After each dispatch Brutus copies back only the candidate cells written since the previous dispatch, not the whole results buffer. Reading the full buffer every chunk would re-copy every earlier cell over and over, turning read-back into an $O(n^2)$ cost that on a full buffer dominates the run; the incremental copy keeps it linear.
What this means for you: You do not need to configure memory settings. If you see the progress bar moving in many small steps, the safety system is working correctly.

If the search is too slow: Adjusting internal GPU memory or chunk parameters is not useful. Instead, curate the peak list, lower the HKL Basis Size, or reduce Peaks to Combine.

Two-Stage GPU Filtering

Most candidate cells are invalid and are discarded on the GPU to save CPU time. Two filters are applied:

  1. Geometric Filter: Cells are rejected immediately if the reciprocal metric is singular (near-zero determinant), if any edge length falls outside $2\;\text{Å} \le a,b,c \le 50\;\text{Å}$, or if the volume is below $20\;\text{Å}^3$ or above the user's Max Volume. Angle bounds are applied on the GPU directly and are system-specific:
    • Monoclinic: $\beta$ is folded to its obtuse representation and required to lie in $89° \le \beta \le 150°$.
    • Triclinic: each of $\alpha, \beta, \gamma$ must lie in $60° \le \angle \le 150°$.
    (These GPU pre-filter windows are deliberately tighter than the CPU refinement stage, which accepts a broader $10°$–$170°$ range; the GPU stage only needs to reject obviously degenerate trial cells cheaply.)
  2. Figure of Merit Pre-filter: For surviving cells, the GPU calculates a fast internal score: the mean of the per-peak normalised residual $\langle |\Delta Q_i| / \Delta Q_{tol,i}\rangle$, evaluated over the first $\min(\max(N_{peaks}, 10),\, 32)$ observed peaks. The lower bound of 10 guarantees a meaningful sample even for very small $N_{peaks}$; the upper bound is the shader constant MAX_FOM_PEAKS = 32, even for very small $N_{peaks}$; the upper bound is the shader constant MAX_FOM_PEAKS = 32, which is the same for all three GPU solvers, so a cell is judged on the same amount of evidence whichever system found it. The matching tolerance $\Delta Q_{tol,i}$ is the per-peak Q tolerance derived from the user's 2θ Error (see Q-Space Tolerance). If the mean normalised residual exceeds FoM Tolerance, the cell is discarded. A running-sum "fail-fast" test aborts the peak loop as soon as the partial sum can no longer pass, so most rejects touch only the first 2–3 peaks. Only cells passing this filter are written to the GPU results buffer for CPU refinement.

    Impurity tolerance. If Impurity Peaks is non-zero, the shader takes a second path: it scores every peak, then averages only the best $(N_{fom} - N_{impurity})$ residuals via a partial selection sort, so the worst-fitting peaks are discarded rather than penalising the cell. Note this path cannot use the fail-fast abort — every peak must be scored before the worst can be identified — so a non-zero Impurity Peaks makes the GPU stage noticeably slower. The impurity count is clamped inside the shader to at most $N_{fom} - 1$, so at least one peak always contributes to the average regardless of what is entered in the UI.

Least-Squares Refinement

Every candidate cell passing the GPU (or CPU) filters is subjected to a full weighted least-squares refinement before being scored and reported. This section details the weighting scheme, zero-point correction, and error propagation implemented in worker-logic.js.

Q-Space Tolerance

The matching tolerance used during peak assignment is not a fixed ΔQ but a peak-angle-dependent quantity derived from the user's stated $\sigma_{2\theta}$. From $Q = (4/\lambda^2)\sin^2\theta$, differentiating:

$$ \Delta Q_i = \frac{8\,\sin\theta_i\,\cos\theta_i}{\lambda^2}\,\frac{\sigma_{2\theta}\,\pi}{360} = \frac{2\,\sin(2\theta_i)}{\lambda^2}\,\sigma_{2\theta}^{(\text{rad})} $$

This means the Q tolerance is largest at intermediate angles (where $\sin 2\theta$ is largest) and converges to zero at $2\theta = 0°$ and $2\theta = 180°$. The implementation adds a small constant ($10^{-9}$) to prevent division by zero.

Least-Squares Weighting Scheme

For a measurement with constant $2\theta$ uncertainty $\sigma_{2\theta}$, the uncertainty in $Q_{obs}$ is:

$$ \sigma_{Q,i} = \frac{2\,\sin(2\theta_i)}{\lambda^2}\,\sigma_{2\theta} $$

The statistically optimal weight for each peak in the least-squares system is therefore:

$$ w_i = \frac{1}{\sigma_{Q,i}^2} \propto \frac{1}{\sin^2(2\theta_i)} $$

This scheme gives high weight to low-angle peaks (which have small $\sigma_Q$) and low weight to high-angle peaks (which have large $\sigma_Q$). This is the opposite of using $w_i = Q_{obs,i}$ (which would over-weight high-angle peaks), and it is particularly important for the zero-point correction: the historical un-weighted approach significantly biased both the cell parameters and the zero shift. A floor of $w_\text{min} = 1/\sin^2(178°)$ prevents weights from blowing up for hypothetical peaks near $0°$ or $180°$.

The least-squares normal equations solved are therefore:

$$ (M^T W M)\,\mathbf{x} = M^T W\,\mathbf{q}_{obs} $$

where $M$ is the design matrix (rows = indexed peaks, columns = reciprocal-lattice parameters), $W = \text{diag}(w_1, \ldots, w_n)$, and $\mathbf{x}$ is the vector of reciprocal parameters. This system is solved numerically via Cholesky decomposition for symmetric positive-definite matrices (implemented in choleskyDecomposition / choleskySolve).

Zero-Point Error Refinement (Two-Round Strategy)

When Refine Zero-Point Error is enabled, a zero-point parameter $\Delta_{zero}$ (in degrees) is added as an extra column to the design matrix. The physical model is that the true $2\theta$ of a peak is $2\theta_{true} = 2\theta_{obs} - \Delta_{zero}$, which shifts the observed Q as:

$$ Q_{true} = \frac{4\sin^2[(\theta_{obs} - \Delta_{zero}/2)]}{\lambda^2} $$

Linearising gives the zero-error design-matrix column element:

$$ \frac{\partial Q}{\partial \Delta_{zero}} \approx \frac{2}{\lambda^2}\sin(2\theta_{obs}) $$

A two-round strategy is used to handle non-trivial zero shifts robustly:

  1. Round 1: Peak assignment is done assuming $\Delta_{zero}=0$. The full weighted LS system (cell + zero) is solved to obtain a first estimate $\Delta_{zero}^{(1)}$.
  2. Round 2: Peak assignment is repeated, this time correcting the observed $Q$ values by $\Delta_{zero}^{(1)}$. The LS system is solved again on the new peak set. This dramatically improves robustness when the zero shift is large enough to cause mis-assignments in Round 1.

The same two-round strategy is applied to every candidate regardless of its origin: GPU-searched systems (orthorhombic, monoclinic, triclinic) and CPU-searched systems (cubic, tetragonal, hexagonal). GPU candidates are refined in parallel by a pool of dedicated Web Workers (see Architecture).

How the zero shift appears on the plot. The measured pattern is drawn exactly as recorded, including its zero error, because that is the only thing actually observed. The refined cell, in contrast, describes the sample once the zero error has been removed, so its calculated lines live in the corrected frame $2\theta_{true} = 2\theta_{obs} - \Delta_{zero}$. To draw a calculated line beside the observed peak it explains, Brutus therefore maps it back the other way, $2\theta_{plot} = 2\theta_{calc} + \Delta_{zero}$. Without that step every blue tick would sit a distance $\Delta_{zero}$ from its red partner — a small, perfectly uniform offset across the whole pattern, easy to misread as a slightly wrong cell. If you see a uniform offset between red and blue ticks, suspect the zero error rather than the cell parameters; the refined value is shown as $Z$ in the solutions table.

Minimum Peaks Required for Acceptance

A solution is accepted only if a minimum number of observed peaks can be indexed. These thresholds are hard-coded and slightly exceed the number of free parameters to ensure the fit is over-determined:

SystemFree parametersMin indexed peaks
Cubic14
Tetragonal / Hexagonal25
Orthorhombic36
Monoclinic47
Triclinic67

Standard Deviation Propagation

The covariance matrix $V = \sigma^2 (M^T W M)^{-1}$ (with $\sigma^2 = SSR/\nu$, $\nu$ = degrees of freedom) is propagated to standard deviations of the cell parameters using analytical derivatives:

Computational Limits & The Scale of Search

Because Brutus performs a brute-force combinatorial search, the number of trial cells can become astronomical. While modern GPUs are incredibly fast, there are hard limits defined by both hardware (WebGPU specifications) and practical time constraints.

The Hierarchy of Limits

The "Total Trials" displayed in the application is:

$$ \text{Trials} = \underbrace{C(N_{hkl}, k)}_{\text{Limit: } \approx 4.29 \times 10^9} \times \underbrace{C(N_{peaks}, k)}_{\text{User dependent}} \times \underbrace{k!}_{\text{Constant per system}} $$

The HKL Combos count $C(N_{hkl}, k)$ is the critical hardware constraint, and the reason is specific to how the shaders work. Rather than materialising an enormous list of HKL tuples in memory, each GPU thread receives a single linear index and reconstructs its own unique $k$-combination on the fly using the combinatorial number system (get_combinadic_indices), which reads a precomputed table of binomial coefficients. Both that linear index and every entry of the binomial table are 32-bit unsigned integers (u32). Consequently the entire enumeration must fit in $C(N_{hkl}, k) < 2^{32}-1 \approx 4.29 \times 10^9$: past that point the thread index and the unranking arithmetic would overflow, silently producing wrong HKL assignments and truncating the search. Widening the host-side table would not help — the limit is intrinsic to the u32 shader.

This translates into a hard ceiling on the HKL Basis Size, distinct for each system:

System$k$Max $N_{hkl}$$C(\text{Max}, k)$
Orthorhombic32,954$4.292 \times 10^9$
Monoclinic4568$4.291 \times 10^9$
Triclinic6123$4.249 \times 10^9$

Each value is the largest $N_{hkl}$ for which $C(N_{hkl}, k)$ still fits in u32 (e.g. $C(124, 6) \approx 4.47 \times 10^9$ already overflows, so Triclinic stops at 123). Brutus enforces this on two levels: a pre-flight check recomputes $C(N_{hkl}, k)$ from your inputs and, if it exceeds the ceiling, disables the Start Indexing button (turning it red with a "Too Large" label) and shows an explanatory status message; and a backstop clamp inside the GPU task caps $N_{hkl}$ at the per-system maximum and warns if it ever had to intervene. In normal use you will hit the button block first — the clamp only guards against a stale or programmatically-set value slipping through.

The u32 ceiling is not the only pre-flight limit. A search can stay well under it and still be impractically slow, because the number of GPU dispatches grows as you raise Peaks to Combine: the per-dispatch thread budget is fixed, so more peak combinations mean fewer HKL combinations per dispatch and therefore more dispatches to cover the same space. At the extreme — triclinic with the maximum peak count — the chunk collapses to a handful of HKLs and the run would need on the order of a billion dispatches: it never errors, it just never finishes. The pre-flight check therefore also estimates the dispatch count and treats it like the u32 ceiling: a warning in the status line above a soft threshold, and a red Start Indexing button above a hard one. If you see a “too slow” block, lower Peaks to Combine first — it has a much larger effect on dispatch count than HKL Basis Size.

System-Specific Breakdown

1. Triclinic ($k=6$)

2. Monoclinic ($k=4$)

3. Orthorhombic ($k=3$)

Practical Implications

Is the "Trials" number real?
Yes. It represents the total number of mathematical systems the GPU would solve. However, most invalid cells are rejected after checking only the first 2 or 3 peaks ("Fail-Fast" optimisation), meaning the GPU rarely performs the full calculation for every trial. The effective throughput is substantially higher than the raw theoretical number implies.

Can you actually run 120 Quadrillion trials?
Technically yes — the Safe Batching system prevents your computer from crashing. However, even at 200 Million trials/sec, 120 Quadrillion trials would take ~19 years to finish.

Practical Limit: For a search to finish in a "coffee break" time frame (5–10 minutes), keep the Total Trials count under 100–500 Billion. For modern desktop GPUs (e.g. Nvidia RTX series), speeds of 100 Mcells/sec are typical for Orthorhombic; Triclinic is slower due to the larger $k$ and the more complex 6×6 linear solve. The size of the GPU buffer, or the accepted solutions, will greatly affect the indexing time because each cell will pass a CPU Least Square fit.

If you obtain no solutions, or too many, consider the following adjustments.

If no solutions (or only poor ones) are found

If you get too many solutions (GPU buffer fills)

The GPU buffer is limited (default 50,000). If it fills up, the search stops early ("Buffer Full"). In that case, if there is no solution:

  1. Tighten the FoM Tolerance (e.g. from 3 to 1.5). This discards more candidates on the GPU, so fewer reach the CPU buffer.
  2. Lower the Max Volume if you can constrain the cell size further.
  3. Increase the Buffer size (e.g. from 50 to 200 kCells) if you suspect the correct solution is present but ranked below the buffer cut-off. Note that larger buffers mean more CPU work during post-processing.
  4. Tighten 2θ Error if your data quality warrants it — this also tightens the Q-tolerance used in matching.

Evaluating Solutions

The indexing search usually produces several candidate cells. To keep the live ledger bounded, whenever the in-memory list grows past 500 entries it is sorted by M(20) and pruned back to the best 400. This reservoir is deliberately large: it exists only to cap memory, not to select results, and pruning too early can discard a strong cell that happens to arrive after the list has briefly filled with near-duplicates. Only the leading few dozen rows are rendered in the Solutions table at any moment, so the larger reservoir costs memory but not redraw time. At the end of the run a final de-duplication pass (applyFinalSieve) removes symmetry-equivalent and transformed duplicates, so the displayed set is typically a few tens of distinct cells ranked by M(20). Selecting the correct one requires interpreting figures of merit and checking the refined fit visually.

de Wolff Figure of Merit: M(20) / M(N)

The primary ranking indicator is the de Wolff Figure of Merit, M(N), calculated from the first $N$ observed reflections (typically $N = 20$ or the full peak list if fewer than 20 peaks are available). It quantifies both positional accuracy and completeness. The implementation uses:

$$ M(N) = \frac{Q_N}{2\,\langle|\Delta Q|\rangle\,N_{calc}} $$

A high M(N) therefore requires both a small average positional error (large 1/$\langle|\Delta Q|\rangle$) and high completeness (small $N_{calc}$ relative to $N$). Both M(20) (from the first 20 peaks) and M(all) (from the full peak list) are computed and stored. The M(20) value is used for ranking; M(all) appears in the PDF report.

M(20) valueInterpretation
> 20Very likely correct. A cell with both high M(20) and plausible chemistry is almost certainly the true solution.
> 10Likely correct, provided the cell volume is chemically plausible.
5–10Plausible; requires further inspection. Check all peaks visually.
< 5Probably spurious; treat with great caution.

A minimum M(20) of 2.0 is required for a solution to be retained at all. Solutions below this threshold are silently discarded.

F(N) Figure of Merit

As a complementary metric, Brutus computes the F(N) figure of merit (Smith & Snyder, 1979). While M(N) is expressed in Q-space, F(N) uses $2\theta$ directly:

$$ F(N) = \frac{N}{\langle|\Delta(2\theta)|\rangle \cdot N_{calc}} $$

A high F(20) indicates a precise fit with low average angular error. A solution with both high M(20) and high F(20) is highly reliable. Both F(20) and F(all) are reported in the PDF.

Least-Squares Refinement and Standard Deviations

For each promising candidate, Brutus performs the two-round weighted least-squares refinement described in the Refinement section. The resulting cell parameters are reported with standard deviations propagated from the covariance matrix.

When Refine Zero-Point Error is enabled, the reported zero correction ($\Delta_{zero}$) represents a systematic misalignment of the diffractometer. A large value (e.g. >0.05° for a lab instrument) may indicate a calibration issue and should be investigated independently.

Duplicate Suppression

To avoid reporting multiple equivalent cells, Brutus generates a canonical "key" for each refined cell: the system name plus the three axis lengths rounded to two decimal places, and (for monoclinic/triclinic) the angles. If a new solution matches an existing key, it replaces the existing entry only if its M(20) is higher.

Space Group Analysis

After a high-quality unit cell is obtained, Brutus can suggest likely space groups based on systematic absences. This serves as input for subsequent structure solution or Rietveld refinement.

Note: For some space groups that have different origin choices but identical extinction rules (e.g. Pmmn), both settings appear in the list. This is intentional — they are not distinguishable by powder diffraction alone, and the reminder is useful for subsequent Rietveld refinement.

Method

  1. Generate unique reflections. Using the refined cell, Brutus computes a complete list of theoretical reflections up to the maximum observed $2\theta$, applying Friedel's law to keep only symmetry-unique $(hkl)$ indices (e.g. $l > 0$, or $l=0,\,k > 0$, or $l=k=0,\,h > 0$ for triclinic).
  2. Index observed peaks (1.5× Tolerance Relaxation). All peaks in the curated list are matched against the theoretical pattern. During space group analysis—and when generating the final PDF report table—Brutus applies a 1.5× safety multiplier to your stated 2θ Error window. For example, an input tolerance of 0.04° opens an effective 0.06° matching window. This relaxation prevents slight sample zero-shifts or high-angle refinement drift from rejecting valid reflections, which would otherwise trigger false systematic absences and disqualify the correct space group.
  3. Retain all indexed reflections. Every uniquely indexed peak is passed to the extinction analysis, including peaks that overlap another theoretical $(hkl)$ within the 2θ Error window. Filtering overlapped peaks out would be self-defeating on exactly the cells that need the most care: when a lattice of high symmetry is described in a lower-symmetry setting, the reflections that overlap are the symmetry-related ones, so such a filter removes most of the pattern (typically 70–80% of peaks) and keeps the biased remnant least able to falsify a centering. Overlap is therefore handled one reflection at a time, through the Ambiguous demotion described below, rather than by excluding peaks up front.
  4. Determine centering and extinctions. The indexed peaks are compared against the extinction rules (centering conditions and glide/screw symmetry conditions) for all candidate space groups. To prevent noise or overlap from unfairly disqualifying a valid space group, Brutus classifies violations into two categories:
    • Hard violations: A strong, real (non-Kα2), unambiguous reflection clearly breaks a rule. These heavily penalise the space group and lower its ranking.
    • Soft violations: A reflection breaks a rule, but the peak is considered unreliable. These do not penalise the space group's ranking, but are listed in the report for the user to review. A violation is flagged as "soft" if the peak is:
      • 2-suspect: Flagged as a likely Kα2 companion (position and intensity ratio — see below). Never applies to Custom or Kα1-only radiation, or to Kα2-stripped data — see the note below on stripping residuals.
      • Weak: The peak height is below 5% of a local maximum, measured over a ±15° 2θ window centered on the peak. The window is local because diffracted intensity falls off steeply with angle (Lorentz–polarisation, Debye–Waller, absorption). Measured against a single global maximum, the same 5% threshold would act as an angle cutoff in disguise: above roughly 90° 2θ nearly every genuine reflection would fall below it, and the whole back-reflection region would silently lose its power to falsify an absence rule. Comparing each peak against peaks at comparable 2θ preserves the crystallographic intent — a truly forbidden reflection should have no intensity — without the angular bias.
      • Ambiguous: The assigned $(hkl)$ violates a rule, but an allowed theoretical reflection explains the observed peak about as well. The alternative must be competitive, not merely present: it qualifies only if it lies within your stated 2θ Error of the observed position, or fits no more than worse than the assigned reflection does. Merely falling somewhere inside the overlap window is not enough — on a crowded pattern almost every reflection has some allowed neighbour, and accepting those would excuse every violation in turn. See the note on wavelength below.

    A note on Kα2-stripped data. Stripping is never perfect. Rachinger-type algorithms subtract a scaled, shifted copy of the pattern, and the error accumulates towards high angle where the doublet separation grows, leaving ripples or “ringing” beside strong peaks. It is reasonable to ask whether disabling the Kα2-suspect test on stripped data lets those artefacts through as hard evidence.

    Enabling the test would not catch them, and would cause harm. The Kα2-suspect test looks for a peak at a specific predicted offset from a lower-angle parent, with an intensity ratio near 0.497. A stripping residual has neither property: the companion it would have matched has already been removed, and residuals sit immediately beside the parent on either side, not at the Kα2 offset. Running the test on stripped data would therefore flag mostly genuine reflections that happen to fall near the old Kα2 position of a strong neighbour — silently removing real evidence.

    What actually covers stripping residuals is the Weak demotion. Ringing is typically a few percent of its parent peak, so it falls below the 5%-of-local-maximum floor and is demoted to soft automatically — and because that floor is local rather than global, it keeps working at high angle where residuals are worst. The Ambiguous demotion catches a further share, since residuals sit close to a real reflection by construction. The gap worth knowing about is a residual beside a very intense peak that clears the 5% floor: that one will count as a hard violation. If a space group is being rejected by one or two violations sitting right next to your strongest reflections, over-stripping is a realistic explanation — mark those peaks as excluded, or re-strip with a lower ratio, before trusting the rejection.

    Consensus override. Each demotion above asserts that one reflection is poor evidence; none of them licenses ignoring a systematic trend. Noise, peak tails and Kα2 ghosts are not correlated with reflection parity, so if many independent reflections all break the same rule, the rule is genuinely broken and the demotions are concealing a real result. Brutus therefore promotes an entire soft pile back to hard when it passes both an absolute and a proportional floor: at least 5 violations and at least 15% of the reflections examined. A couple of weak stragglers will not disqualify a genuine centering; a systematic breach will.

    This matters most for pseudo-symmetric cells. A hexagonal lattice indexed in a monoclinic setting, for example, produces a false body-centering signal: roughly half the reflections appear absent, and without the override every surviving counter-example is demoted as weak or overlapped, so the cell is reported as I-centered with zero hard violations. That false centering then propagates into the Niggli reduction (see Advanced Topics).

    Why Primitive (P) is not listed as an alternative

    A primitive lattice imposes no centering condition: its test function forbids nothing, so P can never accumulate a single violation. It is therefore always among the zero-violation candidates, for every pattern, by construction. Reporting “A-centered (A) or Primitive (P)” would state something structurally guaranteed rather than something measured, and it was never done for F or I, which are collapsed to a single symbol. The Lattice Centering line therefore names only the genuine centering found — P appears there only when no centering survives at all.

    P is still retained internally and passed to the space-group ranking as an allowed centering; removing it there would delete every primitive space group (P212121, Pbcn, Pca21 …) from the candidate list. The question “is this lattice actually primitive?” is answered by the Possible Extinctions list and the space-group ranking, not by the centering line.

    How to read this. A genuine centering produces a condition on general $hkl$ — hkl: $k+l=2n$ for A, $h+k+l=2n$ for I, and so on — and that condition appears in the extinction list. If the extinction list contains only zonal conditions ($h0l$, $0kl$, $hk0$: glide planes) and axial conditions ($h00$, $0k0$, $00l$: screw axes) with no general-$hkl$ rule, the lattice is primitive and the centering shown survived only because no measured reflection happened to fall in its forbidden parity class. Treat such an assignment as unconfirmed and prefer the primitive space groups in the ranking.

    Beyond this, P remains meaningful in one other place: cell reduction, where it simply means “no centering transform applied” before the Křivý–Gruber algorithm runs (see Advanced Topics).

    Redundant rules are collapsed (4n beats 2n)

    A rule is confirmed when nothing observed contradicts it. That test is one-sided, so a zone whose only observed reflection is 004 confirms 00l: $l=2n$ and 00l: $l=4n$ — the second implies the first, so both survive. Reporting both is redundant and actively understates the result: $l=4n$ is the stronger claim, because it additionally forbids 002 and 006, and it is the only one carrying information. Brutus therefore drops, within each zone, any confirmed rule that is implied by another confirmed rule in the same zone, and reports only 00l: $l=4n$.

    Implication is tested by enumeration over the zone rather than algebraically, so compound shorthand (“h, l=2n”) and any modulus are handled without special cases. The collapse is evidence-driven, not cosmetic: if 002 or 006 were actually observed, $l=4n$ would be contradicted and never confirmed in the first place, so the reported rule correctly falls back to $l=2n$.

    Redundancy also arises without any single rule implying another. Zektzerite confirms hk0: $h=2n$, hk0: $k=2n$ and hk0: $h+k=2n$ — the third follows from the first two (even + even is even) but from neither alone, so a pairwise test cannot see it. A rule is therefore dropped when every reflection allowed by the conjunction of all the others already satisfies it, applied iteratively so a mutually-redundant set collapses to a minimal equivalent subset instead of vanishing. Rules from other zones count too, so a general hkl condition can subsume a zonal restatement of itself. Zektzerite's eight reported conditions reduce to four independent ones, with the discrimination unchanged: Abma and Ac2a still explain 4/4, Abmm 3/4, Pccn 2/4.

    The extinction test and the ranking use the same demotions

    A rule is confirmed when every reliable reflection in its zone satisfies it. Reliability excludes Kα2-suspect peaks and peaks below the local intensity floor — and now also applies the same ambiguity demotion the ranking uses: a reflection that formally breaks the rule is forgiven when another calculated hkl within the same peak's tolerance window is not forbidden by it, because the peak could equally well be that reflection.

    Without this the two analyses contradict each other. In zektzerite the four soft violations all have $k+l$ odd, so they blocked confirmation of hkl: $k+l=2n$ even though the centering test had already forgiven the same four peaks and reported A:0(+4 soft). The A-centering condition was consequently missing from the detected list, and P-lattice candidates such as Pccn were never penalised for failing to explain it.

    The same consensus guard as the ranking applies: forgiving one overlap is reasonable, forgiving a systematic trend is not. If five or more reflections — and at least 15% of those examined in the zone — all break the same rule, the rule is treated as genuinely broken however excusable each individual case looks.

    Acentric settings are flagged

    Systematic absences are blind to an inversion centre. Abma (#64) and Ac2a (#41) have byte-identical reflection conditions, so no powder extinction analysis can separate them — only $|E|$ statistics or the solved structure can. Settings lacking a centre of symmetry are therefore marked [acentric] in the list, so an apparent tie is read as the genuine ambiguity it is rather than as indecision by the ranking.

    The same logic applies to ranking. A confirming reflection is weighted by how much of its zone the condition forbids, normalised so an ordinary $2n$ rule keeps weight 1.0 and $l=4n$ scores 1.5. A stricter condition is harder evidence: a reflection such as 004 is consistent with $l=2n$ and with $l=4n$ alike, but surviving the $l=4n$ test excludes far more of the zone. Weighting the two equally would leave P41, P43, P4122 and P4322 tied with P42 and P4222 on identical evidence, and the ranking would fall through to the space-group-number tiebreak — which happens to favour the higher-numbered, less constrained groups. Because the weight is normalised rather than raw, an ordinary $2n$-based score is untouched: the scheme promotes genuinely stronger rules instead of shifting the whole ranking against the rule-free baseline.

    A reflection belongs to every zone it fits

    Zone labels are nested, not exclusive. In ITA notation hk0 means simply $l=0$ with h and k unrestricted, so 330 is an hk0 reflection — and simultaneously an hh0/hhl one, since $h=k$. Likewise 004 is 00l and also hhl (with $h=k=0$), and 300 is h00 and also hk0 and h0l. A reflection must satisfy the condition of every family it belongs to.

    Brutus therefore tests each reflection against all applicable zones rather than against a single classifying label. The membership is deliberately inclusive, which is self-consistent with the printed tables: I-4̄2d's hhl: $2h+l=4n$ specialises at $h=k=0$ to $l=4n$, exactly its separately-listed 00l condition, and Ia-3̄'s hk0: $h,k=2n$ specialises at $k=0$ to $h=2n$, exactly its listed h00 condition.

    Why it matters — the database omits axial conditions on purpose. The bundled cctbx_space_groups_all_settings_v4.json lists only the conditions not already implied by a more general class. Pna21 (33), for example, is stored as just {0kl: k+l=2n, h0l: h=2n}. The axial conditions that International Tables prints for it — 00l: $l=2n$, 0k0: $k=2n$, h00: $h=2n$ — are absent, because each follows by restriction: $k+l=2n$ at $k=0$ is $l=2n$. Inheritance is therefore not an optional refinement; without it those extinctions are never tested at all.

    Checked against structure factors computed directly from the space-group operations, over all reflections to $h,k,l\le 6$:

    Space groupMissed by single-zoneMissed by inclusive
    Pna21 (33)9  (001 003 005 010 030 050 100 300 500)0
    Pbcn (60)90
    Fddd (70)6  (002 006 020 060 200 600)0
    I-4̄2d (122)5  (002 006 110 330 550)0
    P43212 (96)00
    Ia-3̄ (206)00

    The single-zone test carries real weight: across the whole database it decides the verdict for 264 of 527 settings, almost all of them axial reflections that a naive test would treat as allowed. A cross-check over all 527 settings finds no case in which a general-class condition contradicts a listed special-class condition, so inheritance can never manufacture a false absence — it only ever restores a missing one.

    Because inheritance only ever adds applicable rules, each space group is held to every condition that genuinely applies to it. Groups that would otherwise show zero hard violations therefore accumulate them, and fall in the ranking or drop out of the candidate list. That is the test working as intended: those violations were always present in the data, they were simply not being checked against.

    Rewarding settings that explain the observed absences

    Violation counting is one-sided: it only ever looks at reflections that are present, so it can punish a setting for predicting an absence that did not happen, but it can never reward one for predicting an absence that did. The consequence is structural — if one setting's condition set is a strict subset of another's, it can never accumulate more violations, so the least-constrained candidate wins by default.

    Zektzerite is the worked example. The true structure is Cmca (#64); in a cyclically permuted axis setting that is Abma, whose conditions are {hkl: k+l=2n, 0kl: k=2n & l=2n, hk0: h=2n} — exactly Abmm (#67) plus hk0: h=2n. The observed absences include h00: $h=2n$, which Abmm cannot explain: A-centering gives $k+l=0$ for every h00 reflection, which is even for all h, so Abmm predicts odd-h h00 reflections that were never seen. Only the hk0: h=2n of Abma accounts for them.

    Each rule from the extinction analysis is therefore treated as evidence in its own right. A setting is tested for whether it entails each detected condition — is every reflection the setting allows in that zone already required to satisfy it? — and is rewarded or penalised accordingly, weighted by the condition's selectivity and by how many observed reflections back it. The result is reported per setting as explains n/N absences, with any it fails to account for listed explicitly.

    Ranking order is now: hard violations ascending, then unexplained absences ascending, then the confirmation figure of merit, then soft violations, then space-group number. The absence pattern is deliberately promoted above the figure of merit, because the figure of merit counts reflections that are present and is dominated by the well-populated zones that every candidate explains equally; systematic absences are the discriminating evidence in space-group determination. For zektzerite, Abma and Acam explain 4/4 detected conditions while Abmm, Acmm, Abm2, Ac2m and Ammm explain 3/4.

    Hard violations still come first. A genuinely observed reflection in a forbidden class is strong evidence against a setting, so no amount of absence agreement will lift a candidate above one with fewer hard violations. If an expected setting is missing from the list entirely, it has exceeded MAX_VIOLATIONS (2) and the reflections responsible are worth inspecting individually — they are frequently overlaps that the ambiguity demotion failed to rescue under the stricter rule set.

    Swap hkl — correcting an assignment by hand

    Indexing assigns each observed peak to the nearest calculated line. When two lines straddle a peak that choice can be decided by a few thousandths of a degree, and it can be wrong without producing any violation at all: in a permissive space group both candidates are allowed, so nothing flags it. No rule-driven search can find those cases, because such a search only ever sees assignments that some space group forbids. The decision is therefore yours.

    Right-click a solution and choose Swap hkl. The dialog shows the assignment the refit will actually use, including the one-line-per-peak rule: where two peaks compete for the same calculated line the closer one keeps it and the other is shown unindexed, rather than both being fitted to the same reflection. The dialog lists the first twelve peaks with their current indices, observed d, calculated 2θ and the residual — the low-angle end, where a mis-assignment does most damage to the refined cell. Edit any indices you disagree with; edited fields are highlighted. Apply re-fits the cell (and the zero error, if that refinement is enabled) with your assignments imposed and adds the result to the solutions list as an ordinary solution, tagged swap×n. Cancel discards everything.

    Because the new solution is ordinary in every respect — its $M$ and $F$ are computed against the full calculated line list, exactly like any other entry — you can compare it directly with its parent. A correct swap improves the figures; a wrong one degrades them visibly, which is the point. Applying Swap hkl to a swapped solution again lets you build up several changes step by step, with the effect of each one visible.

    The swapped assignments are carried on the new solution and honoured everywhere it is used: the space-group analysis re-runs against them, and the hkl table in the report shows them marked (manual) rather than silently reverting to nearest-line matching — which would otherwise reproduce the very mistake you corrected, since the whole reason for swapping is that the nearest line was the wrong one.

    Several mistakes are caught before refinement: an incomplete row (h, k and l must all be given); assigning a reflection that another peak already holds, which would otherwise be quietly averaged into a distorted cell; two edits aimed at the same peak; and an edit whose peak is no longer inside the 2θ range. The duplicate check is made against every peak the refit will use, not only the twelve rows on screen — a collision with a reflection held by peak thirteen or beyond counts just the same. A manual assignment is honoured even when it is not the nearest line and even when it falls outside the indexing window — overruling the nearest-line rule is exactly what the command is for.

    Worked example. In PbSO4 (anglesite) the peak at 16.423° indexes as (1,0,0) at 0.015° rather than (0,1,1) at −0.002° in the reduced setting. Swapping it and re-refining recovers a −0.012° zero shift that the original assignment had been masking, and returns a cell within 0.002 Å of the published one.

    Nothing is re-indexed automatically

    Indexing is extinction-blind: each observed peak is assigned to the nearest calculated line, whatever that line's parity, because the absence rules do not exist yet at that stage. When two lines straddle a peak, the nearest one can easily be systematically absent. PbSO4 (anglesite, Pnma) is the classic case — the peak at 16.423° takes (0,1,0), 0.019° off, over (1,0,1) at 0.024°, though 010 breaks 0kl: $k+l=2n$ and cannot exist in that space group.

    Brutus does not correct this by itself, and that is a deliberate choice. Re-assigning such peaks automatically once the rules are known, then re-fitting the cell to match, would alter a solution the user never asked to change; with several peaks involved the resulting corrections are hard to interpret; and it would still miss the cases that matter most — an assignment that is wrong while remaining allowed, which happens whenever the candidate space group is permissive enough that neither choice violates anything. An automatic fix would give the appearance of thoroughness while leaving the hardest errors untouched.

    The report therefore shows exactly what the indexer decided, with no silent edits. Where a decision looks wrong, correct it deliberately with Swap hkl, which leaves the original solution untouched and adds the corrected one alongside it for comparison.

    This is not in tension with the automatic combinatorial swap search, which does relabel peaks without being asked. The distinction is what happens to your solution: the automatic search never edits an existing entry, it refits the alternative labelling into a separate cell and offers it alongside as a new candidate to be judged on its own figures. What Brutus declines to do is silently rewrite the assignments of a solution already in your table, or re-derive them from absence rules the assignments themselves were used to establish.

    One consequence is worth knowing when reading figures of merit. The de Wolff figure is $M_{20} = Q_{20} / (2\,\langle|\Delta Q|\rangle\,N_{20})$, and $N_{20}$ counts possible calculated lines. Being extinction-blind inflates that count — for anglesite, 13 distinct lines to $d=3.0$ Å against 9 once Pnma is imposed — so $M_{20}$ is systematically conservative for a centred or glide-rich lattice. A correct assignment is never penalised into being discarded by the figure of merit; it is simply never considered, because assignment happens before the space group is known.

    Why the Ambiguous test is proximity-weighted

    Peak separation follows $\Delta(2\theta) = 2\tan\theta \cdot (\Delta d / d)$. A fixed 2θ window therefore corresponds to a lattice resolution of $\Delta d/d = \Delta(2\theta) / (2\tan\theta)$, which grows without bound as θ falls. Short-wavelength anodes push the whole pattern to low 2θ, so the same window swallows far more reflections — around 20 neighbours per reflection for Ag against 6 for Cu. Under a mere presence test, essentially every violation on Mo or Ag data would be demoted and the extinction analysis would stop working.

    Note that the remedy is not a wavelength-dependent tolerance. Instrumental 2θ uncertainty is set by alignment, sample displacement and detector resolution, and is very nearly independent of the anode: a 0.04° tolerance means 0.04° whichever tube is fitted, so widening or narrowing the window per wavelength would misrepresent the instrument. The competitiveness requirement addresses the problem from the other end, and adapts on its own: on a crowded short-wavelength pattern the extra neighbours are mostly far from the observed position and so do not excuse the violation, while a genuine near-coincidence still does.

    The zone assignment for each $(hkl)$ (e.g. h00, 0k0, hk0, hkl etc.) is used to apply the correct subset of conditions for that zone.
  5. Rank space groups. Based on the crystal system and detected centering, the internal space-group database is filtered. Candidate space groups are evaluated using a Weighted Extinction Confidence score, which actively rewards groups whose systematic absence rules are confirmed by experimental reflections. Higher-order centerings (A, B, C, I, F, R) receive a statistical weighting bonus over primitive (P) lattices because they constrain more reciprocal space. Candidates are sorted primarily by minimizing hard violations, then ranked by their confidence score (with minor penalties applied for soft violations such as low-intensity peaks or suspected Kα2 overlaps). To prevent interface clutter from combinatorial settings, the output is strictly limited to the top 20 highest-ranking candidates, displaying the conventional space-group name consistent with the orientation found by the program.

Interpreting the Output

Space Group MC — refining under each hypothesis

The analysis above scores every candidate against a single cell, and that cell was refined without knowing about any of them: its calculated line list contains every reflection the lattice allows, forbidden ones included. On real data the two weaknesses compound. The least-squares fit is pulled toward reflections that should not be there; peaks near them are then paired with the wrong line; and those mispairings are precisely the violations that condemn the correct group. The scoring is asymmetric as well — a candidate can be punished for predicting an absence that did not happen, but never rewarded for predicting one that did — so a rule set that is a strict subset of another can never lose on violations alone, and the least constrained candidate tends to win by default.

Space Group MC (right-click a solution) asks the opposite question: if this space group were true, how well would the pattern index? For each candidate the forbidden reflections are deleted from the calculated line list, the cell is re-refined against what remains, and the evidence is recounted. Every number in the resulting row (the peak-to-line pairing, the cell, the zero error, $M(20)$, the line count) then belongs to one and the same hypothesis, so the rows compare hypotheses against each other rather than bookkeeping tallies against a cell fitted to forbidden lines.

What a row is: extinction classes

Candidates are grouped before anything is refined. Powder absences reveal only which reflections are missing, so two space groups that forbid exactly the same reflections are indistinguishable in principle — no amount of computation separates them. Brutus therefore collects the settings of the crystal system into extinction classes and refines once per class, listing every group that shares the row.

The grouping happens at two levels, and the difference between them matters:

The second step exists because presenting such rule sets separately — each with its own $M(20)$, as though the numbers meant something — invents a discrimination the experiment cannot support. That is the single most common way a space-group table misleads. Merging on observable behaviour rather than arithmetic is what makes the row count honest: a short scan on a small cell may collapse a dozen candidates into two or three rows, and that collapse is the result.

Calculated lines too close together to be told apart are merged the same way, into what the code calls resolution groups. Two reflections that no observed peak could separate count as one position, so a rule set gets no credit for forbidding a line buried under an allowed neighbour. A group is capped at one matching tolerance wide: it is not enough for each line to sit within tolerance of the one before it, because in a crowded high-angle region a chain of such lines would otherwise merge into a single group many tolerances across, whose two ends are perfectly resolvable from each other. That kind of over-merging costs evidence twice, once by counting several distinct forbidden positions as one and again by letting a peak anywhere in the chain excuse an absence at the far end.

Row labels are extinction symbols

Each row is named by its extinction symbol (Pbca, Pn-a, P2₁/c, Fd-, I4₁/a--), derived by probing what the class forbids along each symmetry direction of the crystal system rather than being read from a table. A dash means no condition in that position; a “?” means the probe found absences it could not name, which is reported rather than guessed at.

Labelling a row with the lowest-numbered group it contains — the obvious alternative — reads like a determination. A user who sees P2₁2₁2₁ at the head of the table takes it for the answer, when what the absences establish is only the extinction symbol and the row may hold half a dozen groups. The Space groups column beside it always carries the full membership.

How the ranking is made

This is the part worth reading carefully, because the criterion is not the one most powder software uses, and the reason it is not is worth stating plainly.

Why $M(20)$ cannot do this job

The de Wolff figure of merit is $M(20) = Q_{20} / (2\langle|\Delta Q|\rangle N_{20})$, where $N_{20}$ counts the possible lines below the twentieth observed one. Deleting any line from the calculated list lowers $N_{20}$ and therefore raises $M(20)$. A more restrictive rule set is rewarded automatically, whether or not the lines it deleted had anything to do with the data.

The usual defence is that an over-restrictive rule set pays for its extra restrictions in violations — it deletes lines that were observed, and gets caught. That defence only holds if every allowed line is actually observed. Real patterns are full of reflections that are simply too weak to see, lost in an overlap, or below the noise. A rule set that happens to delete a handful of unobserved lines gains $M(20)$ for free and incurs no violation at all. Rank on $M(20)$ and you systematically favour the most restrictive class the data failed to contradict, which is not the same thing as the class the data support.

$M(20)$ remains in the table because it is a good figure of merit for a cell, most readers have a feel for its scale, and it can be sorted on like any other column. It is simply not a criterion for choosing between line lists, and so it is not what the table is ordered by. The Score column is, and the table leads with it.

Sorting by $M(20)$ is still occasionally useful — it answers “which of these cells refined best?” — but read it knowing that it puts the most restrictive surviving class on top, because deleting lines is what raises it. That is a different question from “which hypothesis do the absences support?”, and only the second one has an answer the data can defend.

What a nat is

The Score column is a log-odds, measured in nats. A nat is the natural-logarithm unit of evidence: if one hypothesis is $R$ times more probable than another given the data, that is $\ln R$ nats of evidence in its favour. It is the same quantity as a “bit” of evidence, only using base $e$ instead of base 2 (1 nat = 1.443 bits).

The useful property is that nats add. Two independent pieces of evidence worth 1.2 and 0.8 nats together give 2.0 nats, which is why every term below can simply be summed. Convert back to odds by exponentiating:

$$ \text{odds} = e^{\Delta} \qquad\text{so}\qquad \Delta = 1 \Rightarrow 2.7\!:\!1,\quad \Delta = 2.3 \Rightarrow 10\!:\!1,\quad \Delta = 4.6 \Rightarrow 100\!:\!1 $$

Only differences between rows are meaningful; the absolute value of a score carries no information on its own. The Score column therefore shows each row's difference from the best row of its own tier (see falsification below), so the leader reads 0.00 and everything else in that tier is negative. A row at $-0.4$ is barely distinguishable from the leader; a row at $-15$ is out of the running.

The 2.3-nat rule of thumb. Brutus treats a lead of less than about 2.3 nats (10:1 odds) as no lead at all, and says so in words beneath the table: “No decisive winner: N classes lie within 2.3 nats of each other.” If you see that message, the absences in your pattern genuinely do not choose between those classes, and no amount of staring at the $M(20)$ column will change that.
The model

Let $p$ be the probability that a symmetry-allowed line inside the measured range actually produces a detectable peak — at its simplest, the fraction of allowed line positions at which a peak was in fact observed. A clean, strongly diffracting sample measured over a wide range gives $p$ near 1; a sparse or noisy pattern gives $p$ near 0.3. The value is shown in the status line.

$p$ is estimated once and shared by every row, so that no class can inflate the value of its own evidence by discarding lines. Estimating it needs a reference class, and the choice is not innocent: measured on the most permissive class it comes out too low whenever the true lattice is centred, because every systematically absent position then sits in the denominator and can never be observed. Brutus therefore makes two passes — a first estimate from the most permissive class, and a second from whichever class that pass favoured, with every row rescored against the final value. One iteration is enough in practice, and the estimate is clamped at both ends either way.

When peak heights are available $p$ stops being a single number and becomes a per-position quantity, described under intensity weighting below. The formulae that follow are unchanged; only the value of $p$ at each position differs, and the same two-pass calibration then applies to the intensity scale as well.

For a candidate rule set $H$ the score adds up four contributions:

$$ \Delta \;=\; \underbrace{S\!\left(\sum_{j\,\in\,\text{clean}} \ln\frac{1-\varepsilon_0}{1-p_j}\right)}_{\text{absences that held}} \;+\; \underbrace{\sum_{i\,\in\,\text{viol}} \min\!\left(0,\; \ln\frac{\varepsilon_i}{p_i}\right)}_{\text{absences that failed}} \;+\; \underbrace{n_u \ln \varepsilon_u}_{\text{peaks nothing explains}} \;-\; \underbrace{\left[\tfrac{1}{2} k \ln n_{\text{obs}} \;+\; \frac{k(k+1)}{n_{\text{con}}-k-1}\right]}_{\text{parameter charge}} $$
  1. Absences that held. Every forbidden position that is informative — meaning no allowed line shares it, so its emptiness is attributable — and that is in fact empty contributes $\ln[(1-\varepsilon_0)/(1-p_j)]$. Read it as: under $H$ that position should be empty (bar a small spurious-peak rate $\varepsilon_0$); under “no extinction” it would have been empty only with probability $1-p_j$. This is the term that makes the criterion self-limiting. When $p$ is small — most allowed lines unobserved anyway — each clean absence is worth almost nothing, which is correct: in a sparse pattern a missing line proves very little. When $p$ approaches 1, each clean absence is worth a couple of nats. Adding a restriction that deletes an unobserved line therefore buys $\ln[1/(1-p)]$ and no more, instead of the free $M(20)$ increase it used to buy.

    Where no per-position weight is available the contribution is $\bar{p}\,\ln[(1-\varepsilon_0)/(1-p)]$ — the same term scaled by the mean detectability $\bar{p}$. The reason is that conditioning on “this position is empty” preferentially selects the positions where nothing would have shown up anyway, whose true worth is near zero; paying every empty position the population average over-credits the restrictive class systematically. The factor is the same for every row, so it rescales the term without disturbing the order it induces.
  2. Absences that failed. Every violation — an observed peak that only a forbidden line can explain — contributes $\ln(\varepsilon_i/p_i)$, a penalty. The rate $\varepsilon_i$ is set per peak by how believable that particular peak is as real evidence, which is where intensity enters:
    Peak$\varepsilon_i$Rationale
    Strong, clean, unambiguous$\approx N_{stray}/N_{pos}$Near-fatal. Derived from the pattern, not fixed — see below.
    Weak (<5% of the local maximum)0.25More likely a tail, a wing, an overlap or noise than a genuine reflection.
    2 ghost0.50An artefact of the source, not the structure.
    Ambiguous (an allowed line within 1.5× tolerance)0.50The assignment is a judgement call, not a fact.
    Between the weak threshold and full strength the rate is interpolated, so a 6%-of-local peak is not treated as identical to a 100% one. Intensity is measured against a local baseline — the strongest peak within 15° — because diffracted intensity falls off with angle, and a fixed fraction of the global maximum would quietly strip the whole back-reflection region of any power to falsify a rule.

    The three soft rates are genuine probabilities for that peak: each names an identified alternative explanation (a Kα1 parent at the right offset, an allowed line inside the window), so they do not depend on how large the pattern is. The rate for a strong, clean, unexplained peak is different in kind — it is a rate per forbidden position — and is therefore derived as roughly one stray peak across the whole pattern, $N_{stray}/N_{pos}$. It used to be pinned at 0.02, which on a pattern with 700 resolvable positions asserts that fourteen spurious peaks are expected; the derived value there is 0.0014, making each hard violation cost 5.3 nats instead of 2.6. On a small pattern of about 50 positions it returns 0.02 and reproduces the old behaviour exactly. $N_{pos}$ is counted once on the parent lattice and shared, so $\varepsilon_0$ is a property of the pattern and identical for every row.

    The $\min(0,\cdot)$ is not decoration. The ratio $\varepsilon_i/p_i$ exceeds one whenever the detectability model says a reflection at that position would have been invisible, and the peak that contradicts the hypothesis then scores in its favour — up to +6 nats apiece, in a measured case where a deliberately mis-set cell collected seventeen of them and beat the cell that indexed the pattern exactly. When $p_i < \varepsilon_i$ the honest reading is not “this peak supports $H$” but “the intensity model is wrong about this position”: a peak is sitting where the model says nothing could be seen, so it is the model the observation bears on, not the hypothesis. Clamping at zero says exactly that. A violation is at best uninformative, and never support.
  3. Peaks nothing explains. Each unindexed peak contributes $\ln \varepsilon_u$ with $\varepsilon_u = 0.15$. These used to be excluded from the ranking on the grounds that every class carries them equally — true of the old fixed-cell scoring, but false as soon as each class gets its own refined cell.
  4. Parameter charge. $k$ is the number of refined parameters (cell parameters, plus the zero error if you are refining it), $n_{\text{obs}}$ the number of observed peaks, and $n_{\text{con}}$ the number of lines actually constraining the fit. The BIC part, $-\tfrac{1}{2}k\ln n_{\text{obs}}$, is the same for every row and cancels out of every comparison, which is as it should be — charging it on the indexed count instead, as an earlier version did, handed a discount to exactly the restrictive classes that match fewest peaks. What does not cancel is the small-sample correction $k(k+1)/(n_{\text{con}}-k-1)$, which blows up as the number of constraining lines approaches the number of free parameters. That is the “absorb the error into my own zero-point” case: a class with very few allowed lines can fit almost anything, and this is what it pays for the privilege.

The impurity peaks setting is spent here too, on the least believable unexplained peaks first — weakest relative intensity before anything else, and violations before unindexed peaks, since a foreign phase that merely adds an unindexable line falsifies nothing and needs no forgiving. The previous behaviour subtracted the allowance from the raw violation count, which let a class buy forgiveness for the strongest peak in the pattern.

Why the absence term saturates. $S$ in the formula above is a ceiling: $S(x) = C\tanh(x/C) + 10^{-3}x$ with $C = 46$ nats, twenty times the decisive threshold. Absences are the one term that grows with how much a hypothesis forbids rather than with what the pattern actually shows, so a class that deletes half of reciprocal space can accumulate hundreds of nats from positions that were never going to show a peak. A hard cap would flatten every row above it into a tie; $\tanh$ instead squeezes smoothly, staying within a percent of $x$ while $x$ is small compared with $C$ and never exceeding $C$ however much a class forbids. The tiny linear term is there only to keep the ordering strict, since $\tanh$ reaches 1.0 exactly in floating-point arithmetic and would otherwise make the largest rows indistinguishable.

$C$ has to sit above the largest separation absences can legitimately produce, or the ceiling eats real evidence instead of runaway evidence. Two clean classes differing by thirteen informative absences at $p = 0.9$ are 26.9 nats apart; at $C = 23$ that reads as 0.95 nats — an overwhelming result reported as a tie. At $C = 46$ it reads 8.0, while the PbSO4 failure below still cannot buy its way out.
Intensity weighting — a limited borrowing from ExtSym

Everything above uses peak positions. Peak heights carry more, and the statistical method of Markvardsen, David, Johnston & Shankland (ExtSym) shows what can be done with intensities. Brutus borrows part of it. It is worth being exact about which part, because the difference matters.

What ExtSym does that Brutus cannot. ExtSym works from integrated intensities extracted by Pawley refinement, together with their full covariance matrix. That is what lets it ask the sharp question: given all the overlaps, what is the intensity at this forbidden position? Brutus has a peak list, not a profile. Where no peak was picked there is no measurement at all — only absence. The sharp question is therefore out of reach here, and any claim to answer it from a peak list would be pretending.

What is reachable, and why it is worth having. The weakest part of the scoring is that single global $p$. It says a strong low-angle reflection of multiplicity 24 and a weak high-angle one of multiplicity 2 were equally likely to be seen, which is plainly false. Wilson statistics turn that one number into a per-reflection probability, using exactly the information a peak list does carry:

  1. Each observed height is corrected to a quantity proportional to $|F|^2$ by dividing out the Lorentz-polarisation factor and the symmetry multiplicity that the lattice implies.
  2. The corrected intensities are averaged in shells of $s^2 = (\sin\theta/\lambda)^2$. Binning is not a refinement but the definition: the Wilson relation $\langle I \rangle = K \exp(-2Bs^2)$ holds for the shell mean, while individual reflections scatter enormously about it — for an acentric structure $|F|^2$ is exponentially distributed, so its standard deviation equals its mean.
  3. The shell means are corrected for truncation. A peak list contains only what rose above the detection limit, and that limit bites hardest where intensities are weakest, so the surviving means are biased upward more and more with angle. Since $E[y \mid y > t] = t + \mu$ for an exponential variate, subtracting each reflection's own threshold before averaging recovers the true mean without needing to know how many reflections were lost — which is fortunate, because that is unknowable from a peak list.
  4. The weakest peak actually picked sets the detection limit. At each position it converts into the $|E|^2$ a reflection there would have needed in order to be seen, and $p$ at that position is the Wilson tail probability of exceeding it.
Wilson supplies the shape; your data supply the level. The tail probability computed as above turned out to be systematically too generous — measured against synthetic patterns in which 62% of lines were detectable, it averaged about 80%, inflating every absence from 0.94 nats to 1.59. Multiplied over a hundred forbidden positions that was enough for a class with 28 hard violations to outscore one with none: the over-restriction failure this scoring exists to prevent, arriving by a different route. So the shape of $p$ across the pattern is taken from Wilson, while its overall level is pinned to the observed detection rate by a single calibration factor. What the weighting contributes is the relative judgement — this absence is worth more than that one — which is exactly what a global $p$ cannot express.

That calibration follows the same two passes as $p$ itself, and against the same reference class each time. The two are the same physical quantity measured two ways, and pinning them to different hypotheses makes the clean and violation terms disagree about how detectable the pattern is: absences weighted against a no-extinction lattice, compared with a $p$ measured on the winner's own allowed positions. On a synthetic test the level moved by more than a factor of two between the two passes, so this is not a refinement of a detail.

The same scale gives every violating peak an $|E|^2$: its intensity in units of the mean at that angle, corrected for Lp and multiplicity. This is a better variable than raw height for deciding how damning a violation is, because it already accounts for angle and multiplicity, and it appears in the Viol tooltip. A value near or above 1 is unmistakably a real reflection; well below 0.1 is the sort of thing a tail or a little noise produces.

Honest limitations. Peak heights are not integrated intensities: the ratio between them is the peak width, which varies with angle. That matters less than it first appears, because a smooth monotonic error is absorbed almost entirely into the fitted $B$ — which is consequently not a physically meaningful temperature factor and is not reported as one. Overlapped reflections remain a genuine error, since their shared height cannot be divided between them without a profile. And the whole mechanism is switched off, with a note in the status line, whenever the pattern has no heights, too few peaks or too few usable shells — in which case the scoring falls back to the single global $p$ and behaves exactly as described above.

What it buys, measured. On synthetic patterns built with Wilson-distributed intensities, a real Debye–Waller falloff and a detection cut-off, the weighting sharpened a correct Pbca identification from a 3.5-nat margin to 12.1. More importantly, in the degraded regimes where neither method can identify the class, the weighting recognises that the surviving absences are the ones that prove least and reports a tie — where the unweighted score confidently announced the wrong class at 5.9 nats. Refusing to answer is the more valuable half of the improvement.

Falsification comes before the score

A systematic absence is not a statistical tendency. If a space group has an a-glide then $|F|$ is exactly zero for $h0l$ with $h$ odd, and a single genuine reflection there rules the group out however many other absences hold. The likelihood score cannot express that on its own, because it multiplies evidence across reflections: a class that forbids a great deal of reciprocal space accrues absence credit in proportion to how much it forbids, and with enough of it any number of violations can be outweighed.

This was a real failure, not a hypothetical. On a PbSO4 pattern — 192 peaks, about 700 resolvable positions, so only 27% of the possible lines observed — the I-centred class collected 219 empty forbidden positions worth $+65.8$ nats while its nineteen hard violations cost only $-46.6$. It therefore outscored the correct $P2_1/a$, which explained all 48 of its absences and violated nothing.

Rows are therefore ranked in two tiers. A row carrying hard violations that the impurity allowance does not cover is falsified and cannot outrank an unfalsified one, whatever its score. Within each tier the score decides, so the ordering among survivors is still the full likelihood comparison, and the falsified rows are still ordered least-bad first — knowing which groups the data exclude, and by how much, is half the answer. Only hard violations falsify: Kα2 ghosts, peaks below the local weak threshold, and peaks with an allowed line inside the matching window are graded soft, because none of them is a reliable reflection.

Tiering and the saturation ceiling attack the same failure from opposite ends. Tiering stops a pile of absences from burying a reflection that is demonstrably there; the ceiling stops it from burying everything else — soft violations, unindexed peaks, the parameter charge — which no tier protects. Neither alone is enough.

The margin quoted beneath the table is measured within the surviving tier, and among rows refined to the same depth. Comparing an unfalsified row against a falsified one would report a lead that means nothing, since the two are not competing on the same question; comparing a Monte-Carlo row against a stage-1 row measures partly who got more compute, which is not a fact about the crystal. When the restriction bites, the status line says how many rows the comparison covered. And if every class is contradicted, it says so rather than crowning the least bad one — that verdict usually means the cell is wrong, the range too short, or the pattern genuinely contains foreign lines and the impurity allowance needs raising.

A worked reading

Suppose a scan reports $p = 90\%$ on a pattern of about 300 resolvable positions ($\varepsilon_0 \approx 0.0033$), the leading row shows Clean 34/34 with no violations, and the runner-up shows Clean 21/21, also with none. Each clean absence is worth $0.9 \times \ln(0.997/0.10) = 2.07$ nats, so the leader's absences come to 70.4 nats raw against the runner-up's 43.5 — but both are well into the ceiling, which brings them to 41.9 and 34.0. The lead is therefore about 8 nats: decisive, roughly 3000:1, and the extra restrictions are all supported. The raw difference of 27 nats would have claimed $10^{11}\!:\!1$ from thirteen missing lines, which no powder pattern has ever earned.

Now suppose the same pattern gives $p = 35\%$. Each absence drops to $0.35 \times \ln(0.997/0.65) = 0.15$ nats, nothing comes near the ceiling, and the same thirteen-absence difference is worth 1.9 nats — below the threshold, reported as a tie. A single strong violation ($\ln(0.0033/0.35) = -4.7$) would more than wipe it out. Same absences, very different confidence, and that is the whole point: what an absence proves depends entirely on whether a line would have been visible had it been there.

Using the dialog

The scan runs in two stages. Stage 1 gives every class one constrained least-squares refit against its own restricted line list, and keeps the result. That cell is the hypothesis — the cell fitted to the lines the hypothesis allows — so it is accepted whether or not $M(20)$ happened to rise; judging a hypothesis on a cell fitted to reflections it forbids is the exact bias this tool exists to remove. (Earlier versions scored every class on the parent cell as it stood, which re-imported that bias wholesale: the correct class could be demoted at stage 1 and never reach the Monte Carlo pass at all.)

Stage 2 spends the annealing budget on the shortlist. It offers three candidates per class — the parent cell, the least-squares cell, and the annealed one — and keeps whichever scores best, so stage 2 can never leave a class worse off than stage 1 did. Where the score cannot separate two candidates, and it often cannot, $M(20)$ breaks the tie: within a single class the line list is fixed, which makes $M(20)$ a legitimate comparison there and the only one of the two that can see the difference between a good cell and a mediocre one. Three fields control the stage:

FieldMeaning
IterationsAnnealing steps per run, as in Refine MC. Default 600. Below ~200 the walk has not cooled.
RestartsIndependent runs per class. Default 4. Every class walks the same pseudo-random sequence, so a difference between two rows is a difference between hypotheses and not between two draws.
Refine topHow many classes reach stage 2 on score. Up to four more are added for coverage: the most permissive class (the null hypothesis the score is measured against), the class with the most clean absences, and one representative per centering type. The rest keep their stage-1 result, marked with a small dot beside the score.

Progress is reported per class and the run can be stopped at any point with the same button (rows already scored are kept). Kα2-suspect peaks are excluded from the scan exactly as they are from indexing — a ghost sitting a few hundredths of a degree off its parent is precisely the sort of peak that lands on a forbidden line and manufactures a violation against the correct rule set.

The results table has one row per distinguishable class, sorted by Score so that it leads with the hypothesis the evidence favours. Click any header to sort by that column instead, and click it again to reverse; errored rows always sink to the bottom. The row leading on Score stays bold wherever the sort places it, so re-sorting never hides which row actually won.

ColumnMeaning
ClassThe extinction symbol. Rows merged because your data cannot separate them show every symbol they cover, joined by “≡”.
Space groupsEvery group and setting sharing these absences. Hover for the full list; hover the row for the reflection conditions.
LinesDistinct line positions the restricted list produces within the measured range. Lines the scan never reached are not counted, so this agrees with every other tally in the row.
M20de Wolff figure of merit of the refined cell under this class. Describes the cell, not the hypothesis — see above.
ScoreLog-odds against the best row of the same tier, in nats. 0.00 marks the leader; more negative is worse. Within about 2.3 nats of the leader means “not distinguishable”. Falsified rows are measured against the least-bad falsified row, not against the winner, since they are answering a different question.
CleanInformative forbidden positions that are empty, out of the number that are informative at all. 34/34 means every absence this class predicts is borne out; 3/19 means most are not.
ViolHard violations, with soft ones (Kα2, weak, ambiguous) in brackets. Hover for the offending 2θ, its relative intensity and how far it sits into the tolerance window.
UnidxObserved peaks no line explains at all.
ZeroThe zero-point error this class refined to. A row whose zero is far from its neighbours’ is fitting the shift rather than the structure.

Selecting a row and pressing Add as solution inserts its refined cell into the solutions list immediately above its parent, tagged with an SG badge. The badge is green when that row led decisively and amber (SG?) when it did not; its tooltip carries the class, the groups it covers, the reflection conditions, the evidence tally, the estimated $p$ and the margin over the runner-up. That cell was fitted against a restricted line list, so it only means anything alongside the hypothesis it assumes — which is what the badge is for.

Interpreting the table

Limits worth knowing

Trigonal and rhombohedral groups are only ever tested in their hexagonal setting. Brutus classifies cells by their metric, and a trigonal or rhombohedral lattice described in hexagonal axes has $a = b$, $\gamma = 120^\circ$ — which the program calls hexagonal, correctly. There is no separate trigonal system to search and none is needed: every rhombohedral lattice can be indexed on the hexagonal triple cell, and that is what you will get. The consequences are worth stating explicitly:

Monoclinic cells are b-unique throughout. The refinement solves for $\beta$ with $\alpha = \gamma = 90^\circ$, so only the b-unique settings of the monoclinic groups are considered; the a-unique and c-unique descriptions of the same groups carry conditions written for indices Brutus does not produce and are excluded. If your reference data uses one of those settings, transform it before comparing.

What the answer actually is. Systematic absences determine an extinction symbol, not a space group. A class containing five groups means the pattern cannot choose between those five, and the table says so rather than picking one. Going further requires intensities — the statistical test of Markvardsen et al. (ExtSym) works from integrated intensities and can separate candidates this scan must leave tied — or structure solution itself. The scoring above is in the same spirit as that work, reduced to what peak positions alone can support.

Advanced Topics: Enhanced Search and Sieving

Beyond the core indexing routine, Brutus applies several "fishing" strategies and reduction steps to improve robustness and simplify the set of final solutions. These are applied in findTransformedSolutions() after the main search completes.

1. Niggli Cell Reduction & Symmetry Upgrade

For each solution, Brutus computes the Niggli reduced cell — the standardised, most compact cell for the lattice — by applying the Křivý–Gruber reduction directly to the solved cell.

The reduction deliberately does not use the centering inferred from systematic absences. The Niggli cell and the symmetry upgrade are questions about the metric of the lattice actually solved, whereas the detected centering is a space-group guess, and that guess is least reliable exactly when reduction matters most. If a pseudo-symmetric cell is wrongly assigned a centering, applying the corresponding primitive transform extracts a half-volume sublattice bearing no relation to the true symmetry — the reduction then reports a spurious triclinic-looking cell and the symmetry upgrade never fires. Reducing the solved lattice on its own terms makes this step robust to an incorrect absence analysis.

As a worked example: a hexagonal cell forced into a monoclinic setting ($a \approx c \approx 9.40$ Å, $b \approx 6.90$ Å, $\beta \approx 120°$, $V = 527$ ų) is reduced correctly to a cell with two near-equal edges and a 120° angle, revealing the hexagonal metric. Had the false I-centering been applied first, the result would instead have been a 264 ų sublattice with no recognisable symmetry.

After reduction, the code calls getSymmetry() with a deliberately loose tolerance (0.25 Å / 0.25°) to detect pseudo-symmetry: a cell found as monoclinic might, after Niggli reduction, reveal all angles equal to 90° and two equal axes — i.e. it is actually tetragonal. If the detected symmetry is higher than the search symmetry (and the user has selected that higher-symmetry system), a new trial cell is constructed with the correct symmetry and re-refined and re-scored.

A true-primitive reduction of a genuinely centered lattice is available to callers that want it, by passing an explicit centering to the reduction routine. It is not taken from the absence analysis automatically: the centering derived there is a ranked hypothesis, not an established fact, and feeding a hypothesis silently into the cell reduction would make the reported cell depend on a guess the user never sees.

Two cells in the PDF report, and what they are called

Because the internal reduction ignores centering, reporting only that cell would be misleading: for a centered lattice it is not primitive, yet a heading of “Niggli cell” implies that it is. This is a definitional point, not a stylistic one. In International Tables for Crystallography Vol. A, §9.2, the reduced cell is defined as a primitive cell, and Niggli reduction is specified on the primitive basis. Both “reduced cell” and “Niggli” therefore already assert primitivity, so a heading such as “Reduced Cell (Niggli, centering not applied)” asserts and then retracts the same property, and printing it above a line reading “Lattice Centering: Body-centered (I)” is self-contradictory.

The report therefore prints both cells and reserves the Niggli name for the one that has earned it:

For the same reason, the pseudo-symmetry alert reads “HIGHER SYMMETRY DETECTED IN CONVENTIONAL CELL METRIC” and its entries are labelled Conventional-cell metric: that check runs on the first block, not on the primitive cell.

Worked example. A body-centered cubic solution with a = 18.8800 Å, V = 6729.89 Å3 reports a first block identical to the input (the conventional cell is already reduced), and a second block giving the rhombohedral primitive cell: a = b = c = 16.3506 Å, α = β = γ = 109.471°, V = 3364.93 Å3 — exactly half, as expected for I-centering. Seeing both makes it immediately clear that the full-volume cell is the conventional setting and not a failed reduction.

Niggli cells are also used for:

2. Matrix-Based Cell Transformations

A set of seven crystallographic transformation matrices $P$ is applied to each candidate cell. Each transforms the real-space metric tensor $G \to P^T G P$, generating a new lattice description. The matrices include:

The resulting cell is classified by symmetry and, if the system is among those selected by the user, sent to refinement and scoring.

3. HKL Divisor Analysis (Sub-cell Detection)

The list of indexed $(hkl)$ values for each solution is inspected for common divisors. The GCD of all $|h|$ values, all $|k|$ values, and all $|l|$ values are computed separately. If any GCD exceeds 1 (e.g. all $h$ values are even), the program tests a shrunken cell with the corresponding axis divided by that GCD. This detects cases where the indexing algorithm found a super-cell of the true cell.

4. Orthorhombic–Hexagonal Relationship

A hexagonal lattice can sometimes be described as C-centered orthorhombic with $b/a \approx \sqrt{3}$ (or permutations). All orthorhombic solutions are checked against this condition (within 3% tolerance), and potential hexagonal equivalents are generated and evaluated.

5. Combinatorial Swap Search (Swap Fishing)

Mis-assignment is the commonest way a correct cell ends up with a poor M(20): two calculated lines straddle one observed peak and the nearest-line rule picks the wrong one by a few thousandths of a degree. Nothing flags it — in a permissive space group both labels are allowed — and, crucially, the wrong labelling always has the smaller residual at the cell you already have, because that is how it was chosen. The only way to tell is to relabel, refit the cell, and rescore.

A manual swap is available: move one peak to a neighbouring reflection, or exchange the hkl of the two closest peaks. An automatic swap fishing is implemented on all cells that are retained in the final list.

What is enumerated

For every observed peak, every calculated line inside the search window is a candidate label. The search then walks the Cartesian product of those per-peak candidate sets, so any number of peaks may be relabelled simultaneously.

That product is astronomically large, so it is visited cheapest-first: each assignment carries a penalty $\sum_i (\Delta q_i / \mathrm{tol}_i)^2$ measuring how far it departs from the nearest-line labelling, and a priority queue expands states in order of increasing penalty starting from all-nearest. Truncating at any budget therefore leaves the plausible region exhaustively covered rather than sampled. Because cost is set by the fit budget and not by the size of the space, candidate generation is deliberately unrestricted — capping it would only delete reachable answers without saving time.

Two-stage evaluation

A full rescore means regenerating the whole calculated line list, which is too expensive to do hundreds of times per solution. Each candidate labelling is therefore first fitted cheaply: the least-squares design rows are built straight from the labels against a line list generated once for the parent cell, and the normal equations solved. No line regeneration, no nearest-line rescan. These are ranked by weighted RMS residual, and only the best few distinct cells go on to a full M and F recomputation.

Note: the residual ranks candidates; it never accepts one. At a fixed cell a swapped labelling always shows a larger $|\Delta q|$ by construction, so no residual comparison can detect a crossing. It is only after the cell has been refitted to the new labels that the comparison becomes meaningful, and the accepted figure of merit is what decides.

Rounds — why one pass is not enough

Fixing a crossing moves the cell, and moving the cell changes which calculated lines sit near which observed peaks. The second round therefore enumerates a genuinely different candidate set — one the first round could not see at any budget, because those candidates did not exist at the starting geometry. The search runs up to four rounds, re-centring on the improved cell each time and stopping as soon as a round stops gaining.

On the PbSO4 test pattern available in this repository: a single round lifts the best cell from M(20) 43.7 to 49.4 and stops there; four rounds reach 93.6, with the number of indexed lines rising from 44/47 to 46/47.

What counts as better

M(20) is computed on the first twenty peaks and is blind to everything beyond them. A relabelling can inflate it while quietly wrecking the fit to the rest of the pattern, and the combinatorial search is very good at finding exactly such labellings. Candidates are therefore ranked on how many of all N observed lines they index first, and on M(20) only as a tie-break..

Pruning and constraints

Window, budget and what is kept

Alternatives are collected in a window wider than the indexing tolerance (2.5×). A competing reflection just outside tolerance is precisely the interesting case: if it is the right assignment, the cell will shift to meet it. At the strict tolerance the PbSO4 example is invisible — (1,0,0) and (0,1,1) sit 0.048° apart against a 0.040° window, so no alternative would ever be found.

The search runs on the de-duplicated solution list, not the raw one. The candidates arriving from refinement are overwhelmingly redundant — a real PbSO4 run produced 49 solutions that the sieve collapsed to 2, roughly two dozen copies of each distinct lattice — so relabelling all of them repeats the same work over and over to reach the same answers. Sieving first cut that pass from 2.3 s to 0.13 s on the same data, and the saving is spent on a deeper search per surviving parent rather than pocketed: with only a handful of genuinely distinct lattices the budget can be several times the default and still cost a fraction of the old pass. Per round it is bounded by a fixed number of cheap fits and a smaller number of full rescores, and at most a few children are posted per parent so that one solution cannot swamp the table. A child is only posted if its M(20) is at least its parent's, so the search can only ever widen the solution set — it never modifies or replaces a solution you already have.

The scan needs at least as many indexed peaks as the crystal system has parameters to refine, which means six for triclinic and one more again when the zero error is being refined. This is a requirement rather than a tuning choice: below it there is nothing to fit and the search declines to run.

Tuning. All limits live in the SWAP_CFG object in worker-logic.js: window width, alternatives kept per peak, peaks allowed to vary at once, fits per round, number of rounds, and how many parent solutions are searched. Setting DEBUG: true prints one line per parent to the worker console — rounds run, fits spent, lines indexed before and after, and M(20) before and after — which is the quickest way to see whether the search is reaching a given solution at all.

6. Monte Carlo Cell Polish (Refine MC)

Least squares converges to the nearest minimum of the residual, but the M(20) score is not the residual: a cell a few tenths of a percent away can index more lines and score considerably higher. Brutus includes a Refine MC tool (accessible via the right-click context menu on the Solutions tab) to run a symmetry-constrained stochastic local hunt around a refined solution.

Instead of explicitly mutating or manually swapping HKL labels, this algorithm takes continuous random-walk steps in the reciprocal-space parameters (and the zero-error shift, if enabled). As the trial cell geometry shifts, the evaluation functions dynamically pair each observed peak to whichever theoretical line is now closest.

If this dynamic re-pairing captures a better-fitting theoretical reflection and improves the overall M(20) score, the step is accepted. This allows the solution to smoothly escape local minima and naturally resolve its own HKL assignments without requiring manual intervention.

The same walk is reused, under a restricted line list, by Space Group MC, which refines the cell once per space group of its crystal system and ranks the results.

7. Final De-Duplication (Sieving)

After all transformations and re-refinements, Brutus applies a final de-duplication step. If two solutions have volumes within 1% of each other:

The final retained set contains at most 50 solutions, ranked by M(20).

Troubleshooting & FAQ

Error: "WebGPU not supported in this browser"

This error usually indicates a software or configuration issue rather than a lack of hardware capability. Even integrated graphics (e.g. Intel UHD) support WebGPU if configured correctly.

To verify: Type chrome://gpu in your address bar and look for "WebGPU". If it says "Software only" or "Disabled", the GPU search will not function.

Why were no solutions found?

Why is M(20) low although the fit looks good?

The progress bar is stuck / search seems frozen

Test Files

Found a Bug?

If you suspect you have found a bug (e.g. the interface freezes unexpectedly or valid data crashes the GPU search), please report it.

Report a Bug on GitHub

Note: You will need a free GitHub account. The form will ask for your browser version and GPU model to help reproduce the issue.

References

Brutus was developed by Nita Dragoe at Université Paris-Saclay (2024–2025) as a successor to the earlier program Powder (1999–2000). If you use Brutus in your work, please cite:
https://doi.org/10.13140/rg.2.2.18182.84806

For further background on the methodology, the following references are recommended:

  1. M(N) figure of merit
    de Wolff, P. M. (1968). "A Simplified Criterion for the Reliability of a Powder Pattern Indexing." Journal of Applied Crystallography 1, 108–113.
  2. F(N) figure of merit
    Smith, G. S. & Snyder, R. L. (1979). "F(N): A Criterion for Rating Powder Diffraction Patterns and Evaluating the Reliability of Powder-Pattern Indexing." Journal of Applied Crystallography 12, 60–65.
  3. General powder diffraction text
    Klug, H. P. & Alexander, L. E. (1974). X-Ray Diffraction Procedures for Polycrystalline and Amorphous Materials, 2nd ed. New York: Wiley-Interscience.
  4. 2 Elimination Algorithm (q-space differencing)
    Zhang, Z., Lv, B., Zhao, G., Sun, X., Chang, S., Zhang, J., Rodriguez-Carvajal, J., Chen, X. & Feng, Z. (2026). "Parameter-free Kα2 elimination for X-ray diffraction via dual-wavelength q-space differencing." Measurement 276, 121448.
  5. Savitzky–Golay smoothing
    Savitzky, A. & Golay, M. J. E. (1964). "Smoothing and Differentiation of Data by Simplified Least Squares Procedures." Analytical Chemistry 36(8), 1627–1639.
  6. Alternative indexing approaches
    Ito, T. (1949). "A General Powder X-ray Photography." Nature 164, 755–756.
    Werner, P.-E., Eriksson, L. & Westdahl, M. (1985). "TREOR, a Semi-exhaustive Trial-and-Error Powder Indexing Program for All Symmetries." Journal of Applied Crystallography 18, 367–370.
    Visser, J. W. (1969). "A Fully Automatic Program for Finding the Unit Cell from Powder Data." Journal of Applied Crystallography 2, 89–95.
    Le Bail, A. (2004). "Monte Carlo Indexing with McMaille." Powder Diffraction 19(3), 249–254.
    Boultif, A. & Louër, D. (2004). "Powder Pattern Indexing with the Dichotomy Method." Journal of Applied Crystallography 37, 724–731.
  7. Space-group determination from powder absences
    Markvardsen, A. J., David, W. I. F., Johnson, J. C. & Shankland, K. (2001). “A Probabilistic Approach to Space-Group Determination from Powder Diffraction Data.” Acta Crystallographica A57, 47–54.
  8. Previous software
    Dragoe, N. (2001). "PowderV2: A Suite of Applications for Powder X-Ray Diffraction Calculations." Journal of Applied Crystallography 34, 535.