Introduction

This document serves as a scientific and technical guide for the powder5 web application, a tool for the analysis of powder X-ray diffraction (PXRD) data via whole-pattern fitting. This technique, also known as pattern decomposition, is a crucial method in materials science and crystallography for refining structural and microstructural parameters when a complete structural model is either unknown or unnecessary.

This application facilitates the extraction of precise lattice parameters, peak profile information, and integrated intensities of Bragg reflections. It implements two decomposition algorithms: the iterative Le Bail method for rapid and stable convergence, and the Pawley method for unbiased intensity extraction. Peak profiles are modeled using phenomenological functions, including a versatile Simple pseudo-Voigt, an Asymmetric Split pseudo-Voigt, and Anisotropic model (TCH) based on the Thompson-Cox-Hastings formulation with a Stephens model for anisotropic line broadening. The background is modeled using a monotonic cubic spline interpolation between user-defined anchors, which are held fixed or fitted during refinement.

Two independent routes to a structure are provided, and either can be used without the other. Charge flipping is a dual-space method: it solves for a density map with no model at all, and is the first move when the composition is unknown. The Wyckoff search is a direct-space method: it fits atoms on Wyckoff positions of a known composition directly to the extracted intensities, needs no map, and copes better with light atoms in the presence of heavy ones - where charge-flipping ripples can be as large as some light atoms.

Getting Started: Data Input

Analysis commences with the loading of a powder diffraction data file. The application is designed to automatically parse numerous common ASCII-based file formats from major instrument manufacturers and standard crystallographic software.

• Supported Formats: Built-in parsers are included for Bruker (.brml, .uxd), PANalytical (.xrdml), Rigaku (.ras, .rasx), Philips (.udf), GSAS (.esd, .gsa, .std, .xra), and Jade (.mdi). Format detection is by content first and file extension second, so a correctly formatted file will usually be recognised even if it has been renamed. pdCIF files can also be read, at least those saved by this program or Brutus.
• Scientific notation (1.5E+03, 1.2e4) is accepted in the intensity or 2θ column.
• Ordering: data stored in descending 2θ is detected and sorted ascending on load. Duplicated 2θ values and non-finite points are removed, with a note in the browser console.
• Generic Data: Standard two-column ASCII files (.xy, .csv, .txt, .dat, .asc etc.) containing $2\theta$ and intensity values are also supported. The parser accommodates space, comma, or semicolon delimiters. Comment lines prefixed with #, !, ; or // are ignored, as are non-numeric header lines.
• Metadata Parsing: For many instrument-specific formats, instrument parameters such as the X-ray wavelength for Kα1 are read from the file's metadata and used to populate the relevant fields in the user interface. It is incumbent upon the user to verify the correctness of these automatically populated values.

Files & Offline Use

The application is a set of static files with no build step and no server-side component. Every file below must sit in the same directory:

The tree is:

Interactive Data Visualization

The diffraction pattern is rendered in an interactive plot to facilitate detailed inspection of the experimental data and the quality of the model fit. You can hide any of the plots by clicking on its legend at the top of the chart.

All navigation is implemented directly in the application, with no third-party charting plugin, so it behaves identically whether or not you have a network connection.

• Reflection Data: Hovering the cursor near a Bragg reflection marker (tick mark) displays a tooltip containing the corresponding Miller indices ($hkl$) and $d$-spacing.
• Hiding curves: Click an entry in the legend at the top of the chart to toggle that curve.

The status bar

One bar, at the foot of the results panel, reports whatever is running: a Le Bail or Pawley refinement, a charge-flipping run, or a Wyckoff search. It reads, for example

Wyckoff search: Z = 4 (2/3), wave 2/3, restart 1/4, gen 500/1000, wR2 13.0%
Charge flipping: start 2/4, cycle 1300/2000, R 0.2841
Refinement: 64%

A phase that cannot report a fraction - building the Wyckoff tables, compiling a GPU shader - shows a moving stripe rather than a bar stuck at zero. One-off messages (“Stopping…”, “complete - 3 solutions”) appear there too and clear themselves; warnings are also written to the Log tab and to Console, so anything important stays findable after the bar has moved on.

Pattern Decomposition Methodologies

Pattern decomposition enables the fitting of a powder diffraction pattern based on a unit cell and space group, without requiring a full structural model (atomic coordinates). This is essential for the precise determination of lattice parameters and the extraction of integrated intensities, which are requisite for ab initio structure determination.

The Le Bail Method

The Le Bail method is an iterative, sequential algorithm known for its computational efficiency and robust convergence. The process is as follows:

1. Initialization: A theoretical pattern is calculated from the user-supplied lattice, profile, and background parameters (defined by the spline points). Every reflection is given the same starting height. This flat start is deliberate: the method assumes no structural model, so the first decomposition partitions purely on profile overlap. Intensities are never carried over from a previous run, which would make the result depend on history rather than on the data.
2. Seeding passes: Before any parameter is allowed to move, the decomposition alone is iterated a few times so the intensities settle. The partitioning is a fixed-point iteration and converges in a handful of passes.
3. Intensity Extraction (Height Partitioning): The observed net intensity ($y_{i,obs} - y_{i,bkg}$) at each data point is partitioned among the calculated Bragg peaks contributing to that point. The contribution of each peak is proportional to its profile function value at that point. Summing these partitioned intensities for each reflection yields a new set of "observed" integrated intensities. These are then converted back to estimated peak heights using the current profile function's area.
4. Extract, then refine: The intensities are re-extracted at the top of every optimizer iteration and then held fixed for the whole of that iteration - for the baseline residual, for every column of the Jacobian, and for the trial step. The optimizer adjusts only the lattice and profile parameters; the background stays fixed at the spline points throughout.
5. Why not re-extract more often: Re-extracting inside the objective function (i.e. at every function evaluation) flattens the $\chi^2$ landscape: the peak widths can grow substantially while the extraction quietly compensates, so the residual barely responds and the width parameters become indeterminate. Holding the intensities fixed within an iteration keeps the minimum at the true width sharp. Alternating extraction with least-squares - rather than fusing them - is the actual Le Bail method.
6. Convergence: Because the alternation happens inside the iteration loop, a single press of Run Le Bail converges self-consistently; you do not need to click repeatedly for the intensities and parameters to agree. Re-running is useful only to restart from a different algorithm (see Recommended Refinement Strategy) or after changing a fixed input such as the background points.

The Pawley Method

The Pawley method treats the peak height ($I_{hkl}$) of every Bragg reflection as a free variable, fitted at the same time as the cell and the profile. The pattern is linear in every intensity: $$y(x) = \sum_j I_j\,\phi_j(x;\theta) + b(x)$$ so for any given set of non-linear parameters $\theta$ - cell, zero-point, profile widths - the optimal intensities are not something to iterate towards but the exact solution of a linear least-squares problem. powder5 solves that problem outright at every step and lets the minimiser search only $\theta$. This is variable projection (Golub & Pereyra, 1973).

Two consequences matter in practice. The intensities are exactly optimal at every iteration rather than partly converged at the end, so the residual reported is the residual of the model actually fitted. And the search space is the dozen or so genuinely non-linear parameters instead of several thousand, which is what makes parallel tempering usable on a Pawley fit at all.

The normal matrix of the linear step is banded: it has a non-zero entry for two reflections only where their profile windows overlap, and the reflection list is sorted by $2\theta$. The solve is therefore $O(n\,b^2)$ in the bandwidth $b$, not $O(n^3)$ - a few tens of milliseconds and about a megabyte for a pattern with 2500 reflections. The background is a spline through the user-defined anchors, whose heights are refined only if you ask; see Background Modeling for what that costs.

Intensities are bounded at zero using Non-Negative Least Squares (NNLS). An integrated intensity is physically non-negative, but unconstrained least squares can split exactly overlapping reflections into large positive and negative values that sum correctly. Imposing non-negativity via an active-set NNLS algorithm removes this physically impossible region. Crucially, the Estimated Standard Deviations (ESDs) are extracted from the unconstrained normal equations matrix. This ensures that a reflection pinned to zero still reports the true uncertainty of the measurement, rather than falsely claiming an ESD of zero, preserving the statistical honesty required for weak or systematically absent reflections.

• Advantages: The Pawley method is considered more rigorous as it avoids the iterative bias of the Le Bail method, particularly in cases of severe peak overlap. It can yield more accurate and statistically sound integrated intensities (reported as areas) and parameter uncertainties.
• Limitations: the cost is no longer the parameter count - the intensities are projected out rather than iterated - but the information content of the data is unchanged. Where reflections overlap severely their intensities remain strongly correlated, and where they are exactly coincident only the group sum is determined at all. See exactly coincident reflections below. A poor starting cell is still a poor starting cell: the separable form makes the surface in $\theta$ better behaved, not convex.

Minimization Algorithms

The goal of refinement is to minimize the sum-of-squares objective function, $\chi^2 = \sum w_i (y_{i,obs} - y_{i,calc})^2$, where $w_i$ is the statistical weight of each data point (typically $w_i = 1/y_{i,obs}$). This application provides several algorithms to navigate the complex parameter space and find the minimum of this function.

Levenberg-Marquardt (LM)

The LM algorithm is a standard gradient-based method for non-linear least-squares problems. It effectively interpolates between the Gauss-Newton algorithm and the method of gradient descent. By calculating the Jacobian matrix (the matrix of first partial derivatives of the calculated pattern with respect to the parameters), it determines the most efficient path toward the nearest local minimum.

• Characteristics: LM is a local minimizer, exhibiting rapid quadratic convergence when the initial parameters are close to the true minimum. It is the preferred method for final, high-precision refinement and is the only algorithm here that can calculate valid estimated standard deviations (ESDs) for the refined parameters from the covariance matrix.
• Limitations: It is susceptible to converging to a local minimum if the starting model is far from the global solution.
• Parameter bounds: several parameters are bounded - profile widths and Pawley intensities below at zero, the mixing parameter $\eta$ to $[0,1]$, cell angles to $(0^\circ, 180^\circ)$. A step that would cross a bound is clipped to it, and the gain ratio $\rho = \Delta\chi^2_{\mathrm{actual}} / \Delta\chi^2_{\mathrm{predicted}}$ is formed from the step that was actually applied. This matters in Pawley mode, where weak reflections routinely park at zero intensity: the damping parameter $\lambda$ is shared by the whole model, so measuring $\rho$ against a step that was never taken would let a handful of pinned intensities throttle the cell and profile terms as well.
• Stopping at a bound: if an entire step is blocked - the model sitting in a corner of the bound box with the gradient pointing out of it - the refinement stops and reports it, rather than spending its remaining iterations raising $\lambda$ against a wall.
• Pawley Mode: Generally the recommended algorithm for Pawley refinements due to stability and efficiency.

Parallel Tempering (Replica Exchange)

Parallel Tempering, also known as Replica Exchange MCMC, is an advanced stochastic algorithm designed to overcome the slow convergence of traditional search methods on complex landscapes. Instead of a single system, Parallel Tempering simulates multiple copies (or "replicas") of the system simultaneously, each at a different, fixed temperature in a predefined ladder ($T_1 < T_2 < ... < T_N$).

• Mechanism: Each replica evolves independently according to a standard Monte Carlo or Simulated Annealing-like algorithm at its respective temperature. The high-temperature replicas explore the parameter space broadly (high mobility, escaping local minima), while the low-temperature replicas perform a fine-grained search of local minima (high precision).
• The Swap Move: Periodically, the algorithm attempts to swap the entire set of parameters between adjacent replicas (e.g., between replica $i$ at temperature $T_i$ and replica $i+1$ at $T_{i+1}$). The swap is accepted with a Metropolis-like probability that depends on the energies (costs) and temperatures of the two replicas. This crucial step allows a good solution discovered by a high-temperature replica in a distant valley to "percolate down" to the low-temperature replicas, dramatically improving the efficiency of finding the global minimum compared to single-temperature methods.
• Advantages: Significantly more efficient at global exploration than simpler stochastic methods, making it a robust choice for complex problems or when the starting model is highly uncertain (primarily in Le Bail mode).
• Le Bail intensities under PT: Unlike the LM path, which re-extracts at the top of each iteration, the stochastic search extracts the intensities once before the run and holds them fixed throughout. A stochastic walker proposes many rejected moves, and re-extracting on each would make the cost function non-stationary and the Metropolis test meaningless.
• Pawley mode: the intensities are not searched. They are solved exactly at each cost evaluation (see the Pawley method), so a replica walks the same non-linear parameters it would in Le Bail mode and the cost it reports is the reduced objective - the residual after the intensities have been optimised. That is the function whose minimum in $\theta$ is the true minimum, so the Metropolis test and the replica swaps mean what they are supposed to mean.

Guide to Refinable Parameters

This section provides a detailed breakdown of the parameters you can control and refine.

Crystal System & Space Group

These parameters define the crystallographic symmetry of the material.

• The System selection imposes metrical constraints on the lattice parameters (e.g., for Cubic, $a=b=c$, $\alpha=\beta=\gamma=90^\circ$).
• The Space Group selection determines the systematic reflection conditions ($hkl$ absences) used to generate the list of Bragg peaks. The underlying logic for these conditions is consistent with established crystallographic libraries and was taken from Computational Crystallography Toolbox (cctbx).
• Space groups are chosen from a searchable modal listing every setting of all 230 groups. Where a group has several settings (for example No. 62: $Pnma$, $Pmnb$, $Pbnm$, $Pcmn$, $Pmcn$, $Pnam$) each is offered separately, since the reflection conditions differ between them. Rhombohedral groups are offered on hexagonal axes only, because the $d$-spacing formula is written for that indexing.
• Monoclinic settings. All three unique-axis choices are supported. The unique axis is read from the symmetry operators of the setting you select, and the lattice panel then offers the one angle that is free: $\beta$ for unique axis $b$ ($P12_1/c1$), $\gamma$ for unique axis $c$ ($P112_1/b$), $\alpha$ for unique axis $a$ ($P2_1/b11$). The other two are fixed at $90^\circ$ and are neither shown nor refined. The reciprocal metric, the reflection orbits and the multiplicities all follow that same axis, so the three settings of one lattice give identical peak positions and identical intensities. Switching between settings carries the angle across, so the cell is not lost.
• Startup default: the application opens on $Pnma$ (No. 62, standard setting) with an orthorhombic cell of $a = 8.478$, $b = 5.397$, $c = 6.958$ Å. This is only a starting point - change it to match your material before refining.

Reflection List & Multiplicities

The Bragg peak list is built from the rotation operators of the selected space-group setting. Every candidate $hkl$ inside the resolution sphere is mapped onto its orbit under the Laue group; the orbit is kept once, labelled by a canonical member, and its size becomes the multiplicity $m$ that scales the intensity, $I \propto m \cdot LP \cdot |F|^2$. Systematic absences are then tested on that representative, which decides the whole orbit, since an absence is a property of the orbit and not of any one index triple.

Working from the operators rather than from a fixed index range per crystal system matters wherever two Laue classes share a system, because the two classes have genuinely different sets of independent reflections:

In the hexagonal and trigonal systems the special forms are set by the Bravais–Miller index $i = -(h+k)$ as much as by $h$ and $k$: a reflection lies on a mirror whenever any two of $h$, $k$, $i$ are equal in magnitude, so $(1\bar{2}l)$ and $(2\bar{4}l)$ are special even though $|h| \neq |k|$. The orbit calculation accounts for this automatically.

Operators are taken from the sym_ops field of the space-group database (or the older rotations field, which serves equally well here - only the rotation parts affect a reflection orbit). The closed group is checked against the order_p recorded for the setting before it is used. If the database carries no operators, a built-in table of the Laue-class generators is used instead; if the Laue class itself cannot be identified, the generator falls back to the lowest Laue class of the crystal system. That last case is safe rather than wrong: orbits split into smaller ones at identical $d$, so peak positions and total intensities are unaffected, and only the number of listed entries changes. The browser console records which source was used whenever it is not the first.

Instrumental Parameters

Found under the "Sample" tab, these parameters model the diffractometer configuration.

• Radiation 1/2 (Å) & Ratio: Defines the X-ray source. For divergent-beam laboratory instruments, a Kα1/Kα2 doublet is typically used. For synchrotron radiation, the Ratio is set to 0.
• Polarisation: Selects the polarisation model used in the Lorentz–polarisation factor - Lab (default), Synchrotron, or None. A second field appears beside it: the monochromator angle $2\theta_M$ for a laboratory source, or the polarised fraction $f$ for a synchrotron. This setting does not affect the fit; it governs how an extracted intensity is converted into $|F|$. See Lorentz–Polarisation Factor.
• Zero: A refinable parameter that corrects for instrumental zero-point error in the $2\theta$ axis. It is highly correlated with lattice parameters and must be refined with caution.
• 2θ Min / Max: These sliders define the refinement range. It is standard practice to exclude regions of low signal-to-noise or known artifacts from the calculation.

Lorentz–Polarisation Factor (Lp)

The integrated intensity of a powder reflection is not $|F|^2$ alone. It carries two geometric weights - the multiplicity $m$ of the reflection and the Lorentz–polarisation factor $Lp$ - on top of an arbitrary overall scale $s$:

$Lp$ is therefore the last thing standing between an extracted Le Bail or Pawley intensity and a structure factor, and it is the reason the reflection table in the report prints $m$, $Lp$ and $|F_o|$ next to the intensity rather than the intensity on its own.

The Lorentz factor

For a conventional $\theta$–$2\theta$ powder scan, Powder 5 uses

This lumps together the rate at which a crystallite passes through the reflecting condition and the fraction of the Debye–Scherrer ring intercepted by the detector. It is often written $1/(4\sin^2\theta\cos\theta)$; the factor of 4 is a constant and is swallowed by $s$, so it makes no difference to anything reported here. $L$ has no adjustable parameters and is applied for every polarisation setting, including None.

The polarisation factor

All three models are the same expression with one constant $K$, normalised so that $P(0) = 1$ in every case:

With $2\theta_M = 0$ the Lab setting reduces to $P = \tfrac12(1 + \cos^2 2\theta)$, the classical unpolarised expression, and this is what Powder 5 uses by default. Note that a fully perpendicular-polarised synchrotron beam ($f = 1$) gives $K = 0$ and therefore exactly the same arithmetic as None; the two settings are kept distinct because they describe different geometry and are labelled differently in the report.

Choosing $2\theta_M$ for a laboratory monochromator

$2\theta_M$ is the take-off angle of the monochromator crystal at the wavelength in use, not a free parameter. It follows from Bragg's law applied to the monochromator itself,

where $d_M$ is the interplanar spacing of the reflection the crystal is cut for. The same crystal therefore gives a different angle on every anode: a graphite (002) monochromator sits at 26.6° on Cu but at 12.2° on Mo, and $K$ rises from 0.800 to 0.956 with it. Change the tube and this field must change too.

Each cell below gives $2\theta_M$ in degrees, with the resulting $K = \cos^2 2\theta_M$ in brackets. Enter the angle, not $K$ - Powder 5 squares the cosine for you. Angles are computed from the K$\alpha$ weighted-mean wavelengths ($\alpha_1{:}\alpha_2 = 2{:}1$); a monochromator that passes K$\alpha_1$ alone shifts $2\theta_M$ by less than 0.05°, which is far below anything $P$ can resolve.

$d_M$ is quoted so you can check an optic that is not listed: Ge $a = 5.6575$ Å and Si $a = 5.4310$ Å give $d_{hkl} = a/\sqrt{h^2+k^2+l^2}$, and graphite's $d_{002}$ is half its $c$ axis.

Choosing $f$ for a synchrotron

Bending-magnet and undulator radiation is polarised in the horizontal plane, with $f$ typically between 0.90 and 0.98 in the plane of the orbit. What matters is the orientation of the diffractometer relative to that:

• Vertical scattering plane (the usual arrangement for high-resolution powder diffraction): the electric field is perpendicular to the scattering plane, so $f \approx 0.95$–$1.0$ and $P$ is nearly flat. The correction is almost pure Lorentz, and this is precisely why the geometry is chosen.
• Horizontal scattering plane: the field lies in the scattering plane. Enter $f \approx 0.05$–$0.1$, which gives $P \to \cos^2 2\theta$ - a severe correction that goes to zero at $2\theta = 90^\circ$. Getting this backwards is a much worse error than misjudging $f$ by a few per cent.

How large is the effect?

$Lp$ falls steeply with angle whatever the model, and the polarisation term adds a further factor of two by $2\theta = 90^\circ$ for an unpolarised source. Values for Cu K$\alpha$:

$Lp$ spans a factor of about 90 between $2\theta = 10^\circ$ and $90^\circ$. Omitting it does not perturb a structure-solution calculation slightly - it weights the first few peaks roughly two orders of magnitude too heavily, and a charge-flipping map built on uncorrected intensities is dominated by them. The choice between models is a much smaller effect: at most a factor of two, and concentrated near $2\theta = 90^\circ$. Getting the correction on at all matters more than getting the model exactly right, but the models are cheap to select correctly.

Where it is applied

All of these read one shared calculation, evaluated at the corrected $2\theta$ - the same angle printed in the reflection table, including the zero shift and any displacement or transparency term. The charge-flipping worker takes the per-reflection value computed by the refinement rather than recomputing its own, so the map, the report and the prior cannot disagree about it. The charge-flipping summary panel reports the model the worker actually used, which is the check that they have not.

Background Modeling (Spline Interpolation)

The background contribution is modeled using a monotonic cubic Hermite spline interpolation between a series of user-defined points (spline points or knots). This approach provides flexibility and ensures a smooth, physically realistic background shape without introducing refinable background parameters into the least-squares minimization. The background shape is therefore considered fixed during the refinement process based on the current spline points.

Control of the background spline is located under the "Background" tab:

• Auto-generation: The application automatically estimates background points immediately upon loading a data file. You can adjust the density of these points using the Auto-points slider (10-40 points). Adjusting the slider automatically recalculates the points based on local intensity minima within intervals distributed across the current 2θ Min/Max slider range. The points at the exact 2θ Min and Max slider positions are always included and fixed to these $2\theta$ values.
• Manual Addition: Add individual points by holding Ctrl and clicking on the chart. The closest experimental point will be added to the list, provided it's within the current slider range and not an edge point.
• Editable List: The generated and manually added points appear in a list below the controls.
  • You can directly edit the $2\theta$ and Intensity values for any point, except for the $2\theta$ values of the first (Min) and last (Max) points, which are fixed by the sliders. Edits trigger recalculation of the spline.
  • Points can be deleted using the × button, except for the first and last points.
• Chart Display: The spline points can be toggled on/off on the chart using the "Show Points on Chart" checkbox. The calculated spline curve is always shown.

Whether you tick it changes what the numbers in the report mean.

The statistics block states which of the two applied on every run.

Simple pVoigt

This function models the peak shape as a linear combination of a Gaussian and a Lorentzian function: $pV(x) = \eta L(x) + (1-\eta)G(x)$. The angular dependence of the Full Width at Half Maximum (FWHM) for each component is modeled empirically. $$H_G^2 = GU \tan^2\theta + GV \tan\theta + GW + GP / \cos^2\theta$$ $$H_L = LX / \cos\theta$$

• GU, GV, GW, GP: Parameters describing the Gaussian FWHM ($H_G$). The terms are associated with strain ($GU$), instrumental factors ($GV, GW$), and particle size effects ($GP$). Note that $GU$, $GW$, and $GP$ should physically be non-negative.
• LX: Describes the Lorentzian FWHM ($H_L$), primarily associated with crystallite size broadening ($LX > 0$).
• eta: A simple linear mixing parameter ($0 \le \eta \le 1$; $\eta=0$ for pure Gaussian, $\eta=1$ for pure Lorentzian).
• shft & trns: Corrects for peak position shifts due to sample displacement and transparency, respectively.
  Unit and Scaling for the shft parameter:
  The refined shft parameter is a dimensionless, scaled coefficient, not a direct physical length. Its relationship to the physical specimen displacement ($s$) and the goniometer radius ($R$) is defined as follows:

  • The physical peak shift in radians is: $\Delta(2\theta)_{\text{rad}} = -2 \frac{s}{R} \cos(\theta)$
  • The program calculates this shift (in degrees) using the formula: $\Delta(2\theta)_{\text{deg}} = -(\text{shft} / 1000) \times \cos(\theta) \times (180 / \pi)$
  • Therefore, the relationship is: $\frac{\text{shft}}{1000} = \frac{2s}{R}$
  • To find the physical displacement $s$ from the refined parameter, use: $s = R \times (\text{shft} / 2000)$.
    Example: For a typical instrument with $R=240$ mm, a refined shft value of 1.0 corresponds to a physical displacement $s$ of $240 \times (1/2000) = 0.12$ mm.

TCH (Size/Strain/Aniso)

This is a physically rigorous profile function adapted from GSAS, convoluting a Thompson-Cox-Hastings (TCH) pseudo-Voigt with a Stephens model for anisotropic strain broadening.

Isotropic Broadening (TCH Model)

The TCH formulation models the FWHM of the Gaussian ($H_G$) and Lorentzian ($H_L$) components based on physical contributions to line broadening: $$H_G^2 = U \tan^2\theta + V \tan\theta + W$$ $$H_L = X \tan\theta + Y / \cos\theta$$

The total FWHM ($H$) and pseudo-Voigt mixing parameter ($\eta$) are then derived from these components using polynomial approximations. The final shape is $pV(x) = \eta L(x, H) + (1-\eta)G(x, H)$, where both functions share the same convoluted FWHM.

• U, V, W: Gaussian broadening parameters related to strain ($U, V$) and instrumental resolution ($W$). Physically, $U$ and $W$ should be non-negative.
• X, Y: Lorentzian broadening parameters related to strain ($X$) and crystallite size ($Y$). Physically, $X$ and $Y$ should be non-negative.

Peak Asymmetry

The S/L and H/L parameters introduce an angle-dependent asymmetry, primarily correcting for axial divergence effects at low $2\theta$.

Anisotropic Broadening (Stephens Model)

Anisotropic microstrain, where broadening varies with crystallographic direction, is modeled by adding terms to the Lorentzian component ($H_L$) that are dependent on the Miller indices ($hkl$). The model is a fourth-order polynomial in the reciprocal lattice vectors.

The refinable parameters (S400, S040, etc.) are the non-zero, symmetry-unique coefficients of this polynomial. The application automatically applies symmetry constraints based on the Laue class of the selected space group (e.g., for cubic, $S400=S040=S004$).

Split pVoigt (Asymmetric)

This profile function is a modification of the Simple pseudo-Voigt designed to model asymmetric peaks (e.g., from axial divergence or stacking faults). It achieves this by defining independent sets of profile width parameters for the left side (at $2\theta$ values less than the peak center) and the right side of the peak.

The shape is still a linear combination $pV(x) = \eta L(x) + (1-\eta)G(x)$, but the $H_G$ and $H_L$ parameters used in the calculation depend on whether $x$ is to the left or right of the peak center.

Gaussian Broadening (Left & Right)

The Gaussian FWHM ($H_G$) for each side is modeled as:

$$H_{G, \text{side}}^2 = GU_{\text{side}} \tan^2\theta + GV_{\text{side}} \tan\theta + GW_{\text{side}}$$

(Note: This model does not use the $GP$ term used in the Simple pVoigt profile.)

• GU-L, GV-L, GW-L: Parameters describing the Gaussian FWHM for the left side of the peak.
• GU-R, GV-R, GW-R: Parameters describing the Gaussian FWHM for the right side of the peak.

Lorentzian Broadening (Left & Right)

The Lorentzian FWHM ($H_L$) for each side is modeled as:

$$H_{L, \text{side}} = LX_{\text{side}} / \cos\theta$$

• LX-L: Describes the Lorentzian FWHM for the left side, primarily associated with size broadening.
• LX-R: Describes the Lorentzian FWHM for the right side, primarily associated with size broadening.

Peak Shape & Position

• eta (Mixing) (param: eta_split): A simple linear mixing parameter ($0 \le \eta \le 1$; $\eta=0$ for pure Gaussian, $\eta=1$ for pure Lorentzian). This single value is used for both sides of the peak.
• shft (Displ.) (param: shft_split): Corrects for peak position shifts due to sample displacement.
• trns (Transp.) (param: trns_split): Corrects for peak position shifts due to transparency.
  Unit and Scaling for the shft_split parameter:
  This parameter is a dimensionless, scaled coefficient, identical in function to the shft parameter in the Simple pVoigt profile.

  • The physical peak shift in radians is: $\Delta(2\theta)_{\text{rad}} = -2 \frac{s}{R} \cos(\theta)$
  • The program calculates this shift (in degrees) using the formula: $\Delta(2\theta)_{\text{deg}} = -(\text{shft\_split} / 1000) \times \cos(\theta) \times (180 / \pi)$
  • Therefore, the relationship is: $\frac{\text{shft\_split}}{1000} = \frac{2s}{R}$
  • To find the physical displacement $s$ from the refined parameter, use: $s = R \times (\text{shft\_split} / 2000)$.

Forward Simulation: Calculating a Pattern from a Structure

Every other mode in this program starts from a measured pattern and works backwards. Le Bail and Pawley extract intensities from it; charge flipping and the Wyckoff search turn those intensities into a structure. Simulation runs the same physics forwards — atoms in, pattern out — and needs no data file at all.

It is selected from the Method box on the control panel (the block previously headed “Refinement”), where it sits last, after Levenberg–Marquardt and Parallel Tempering. Those two are minimisers; this one is not, and the controls that only make sense during a refinement hide themselves when it is chosen.

What it shares with the fit

The simulation contains no profile function, no background, no Lorentz–polarisation factor, no hkl generator and no peak-area integral. Every one of those already exists in this program, is already used by the refinement, and is called here unchanged:

This is the whole point of the design. The header of js/profile.js records what happened the last time this program had two implementations of the pattern: the preview drew one and the fit minimised another, and every report was self-consistent about the wrong answer. A simulation with its own profile function would be a third such copy, and the failure would be worse than a drifting fit — a simulated pattern that cannot be compared with a measured one is not a simulation of anything. Because the two share the chain, a simulated pattern and a Rietveld calculation of the same model are computed through identical code.

The only new physics

The one thing the fit did not already have is the structure factor. For each reflection, summed over every atom in the unit cell:

$$F(hkl) = \sum_j n_j\, f_j(s)\, \exp\!\left(-B_j s^2/4\right)\, \exp\!\left(2\pi i\,(hx_j + ky_j + lz_j)\right)$$

with $n_j$ the site occupancy, $f_j$ the atomic scattering factor and $B_j$ the isotropic displacement parameter. $|F|^2$ is then handed to reflectionHeightsFromFsq(), which applies the multiplicity, the Lorentz–polarisation factor and the Kα₂ doublet exactly as the Pawley extraction does in reverse.

Occupancy and the Debye–Waller factor both belong to the site and depend on the reflection only through $s$, so they are folded into the site's scattering factor once per reflection rather than once per atom. For a high-symmetry structure that is the difference between 192 exponentials and one.

The Atoms tab

A new tab appears in the left panel while Simulation is active. It holds the asymmetric unit — the only thing none of the other tabs already carry. The cell, the space group, the wavelengths, the zero point and the peak shape all come from the Sample and Profile tabs as usual.

Each atom is one collapsed line: its label, its Wyckoff position, and its fractional coordinates. Clicking the line expands it into the editable fields; clicking again folds it away. Expand all / Collapse all appears once there is more than one atom. A newly added atom opens automatically, since you are almost certainly about to adjust it; a structure loaded from a CIF arrives collapsed, because expanding fifty sites at once buries the thing you actually want to see.

The collapsed line keeps the coordinates visible because that is what gets compared down a list. Occupancy and $B$ appear there only when they differ from 1.000 and 0.50, so an unusual value stands out instead of being lost in twenty identical repetitions.

Every field is live: x, y, z, occupancy, $B_{iso}$ and the element are edited in place and the pattern follows through a short debounce. Watching an intensity change as an atom moves is how a structure gets understood, and no dialog that has to be opened and dismissed per edit can give you that.

Physical ranges

Adding an atom

Add atom opens a picker. Type into the search box to filter the scattering table: T gives Ti, Ti²⁺, Ti³⁺, Ti⁴⁺, Tc, Te, Tb, Tm, Ta, Tl, Tl¹⁺, Th. Ions are listed immediately after the neutral atom they came from rather than in a separate block, because someone typing “T” is looking for an element and then deciding its charge. Matching is on the prefix, not a substring, so typing O gives O and O²⁻ and not Co, Mo, Ho and Po. Typing a complete label selects it, so Ti4+ followed by Enter adds the atom without touching the mouse.

The dialog shows, live, which Wyckoff position the coordinates you are entering will land on — you find out what the symmetry does to them before committing, rather than asserting a position and meeting the mismatch later in the pattern.

The Wyckoff position is derived, never entered

The letter beside each atom's name is computed from its coordinates and the current space group. It is not a label you assert. The pattern is calculated from the coordinates either way, so a typed letter could only ever be an assertion the coordinates might not support; showing the derived value instead makes the constraint visible. Move an atom onto a mirror and watch the multiplicity halve.

The multiplicity decides and the projection only breaks ties. The orbit of the point under the full operator set is the multiplicity — it is the number of atoms the structure factor will actually sum over — so the letter is looked up among the positions of that multiplicity and nowhere else. Several positions can share one multiplicity (Pnma has 4a, 4b and 4c), and those are separated by projecting onto each candidate's special_op. The test runs over every symmetry image of the point, not just the coordinates as typed, because special_op describes one representative: a point at $(0,0,\tfrac{1}{2})$ in $P112_1/a$ is on 2c, but its representative is $(\tfrac{1}{2},0,0)$ and only an image of the point ever reaches it.

Centring is handled by crossing the setting's sym_ops with its centring_translations and deduplicating, which is correct whether or not the database already folds centring into the operator list. Dropping the centring would not be a subtle error: an F-centred cell has four times the atoms of its primitive subgroup, so the pattern would come out with the right peak positions and uniformly wrong intensities — the hardest kind of error to see.

Load CIF

Load CIF reads an ordinary crystal-structure CIF and populates the atom list. It is a structure reader, not a pattern reader — parsePdCifFile() in js/data_io.js handles powder profiles, this handles the _atom_site_ loop and the cell.

• Standard uncertainties are stripped, both from tags and from inside loops: 0.30478(5) is a value, not a syntax error.
• Semicolon text fields are folded before parsing, so a publication title containing the words loop_ or data_ cannot derail the reader.
• Only the first data block is read. Multi-block CIFs are common, and silently merging their atom lists would produce a structure that exists in none of them. The status line says how many blocks were found.
• $U$ is converted to $B$ when the file carries _atom_site_U_iso_or_equiv, using $B = 8\pi^2 U$.
• The scatterer is resolved from _atom_site_type_symbol first, verbatim, since ionic spellings such as Ti4+ already match the table. Failing that the charge is dropped, because a table without Ti⁴⁺ still scatters perfectly well as Ti and a site contributing nothing is a worse answer than one with a slightly wrong form factor. The site label is a last resort — Ca1 carries the element only by convention, and that convention breaks often enough (M1, OW) that guessing silently would be worse than reporting the site as unidentified. Two-letter readings are tried before one-letter ones, so OW1 resolves to O.

The cell and the space group are applied with the atoms. Fractional coordinates are fractions of a particular cell under a particular set of operators and mean nothing apart from those two, so the whole triple loads together. The space group goes first, because it decides which lattice boxes exist — a cubic group leaves only a, and writing the angles before the system is set would put them into inputs that are about to be destroyed.

Loading a CIF replaces the atom list. A CIF is one whole structure, and merging it into whatever was already there would build a phase that exists in neither file.

Pattern controls

The 2θ range is the one on the Sample tab. Those sliders are normally disabled until a file is loaded, because outside Simulation they select a slice of the data; here they are enabled with no file open, since the range is the only thing that says how much pattern to compute. With data loaded they are left completely alone, which is what makes the simulated pattern come out on the same interval a fit would use.

Counting noise

The physics settles the design here. A diffractometer counts photons, and photon counting is a Poisson process: the variance of a measured intensity is the intensity, so the scatter goes as $\sqrt{I}$. That single fact is what makes a simulated pattern behave like a real one — strong peaks noisy in absolute terms and clean in relative terms, weak peaks the other way round, and a background whose own scatter sets the detection limit.

Gaussian noise of one fixed width reproduces none of that. It would swamp the weak reflections and be invisible on the strong ones, and any program weighting by $1/\sigma^2$ — which is what a Rietveld refinement does — would be handed a wrong weight at every point.

So the field is expressed in units of the counting statistics themselves:

Noise is applied to peaks plus background, because the detector counts both: a background of 10 carries its own scatter of about 3, and that is what sets how weak a peak can be and still be seen.

Exporting: the report button writes a pdCIF

On the Plot tab, while the Method is Simulation, the report button produces a pdCIF instead of a PDF. Every other tab, and every other Method, prints exactly what it printed before.

A PDF of a calculated curve is a picture of something the program can state exactly, and the numbers are the point of a simulation — they are what gets compared against a measurement, fed to another program, or checked by hand. The file carries the 2θ axis, the intensities, the reflection list with |F|², the cell, the space group with its Hall symbol, and the asymmetric unit.

A simulation with no measurement must not claim one. _pd_meas_2theta_scan and _pd_meas_intensity_total are measured quantities; writing a calculated axis under a _meas_ tag with a column of ? beside it would state that a scan was performed and its intensities lost, which is a different thing from a scan that never happened.

With noise switched on, though, the pattern is a synthetic observation — that is what noise means — so it is written as the observed column with the noise-free curve alongside as the calculated one. That is exactly the pair another Rietveld program needs: fit the observation, and check that it recovers the model that produced it. The block carries an explicit “these intensities are SIMULATED, not measured” comment, because tags alone would be indistinguishable from a real scan and a file that quietly claims to be a measurement is the kind of thing that ends up in a paper.

The atom-site loop uses standard core-CIF tags rather than pdCIF ones, which is correct: the structural model and the powder profile are different dictionaries describing the same block, exactly as in a Rietveld pdCIF. $B$ is written rather than $U$, because $B$ is what the panel holds and what the structure factor used; converting here would put a different number in the file from the one the pattern was computed with, distinguishable only by knowing which convention the writer chose. _atom_site_Wyckoff_symbol gets the letter and _atom_site_symmetry_multiplicity gets the number, per the dictionary.

Scattering factors

Tables live in scatters/, with index.json naming what is available and which is preferred. They are five-term Gaussian parameterisations, $f(\mathrm{stol}) = c + \sum_i a_i \exp(-b_i\,\mathrm{stol}^2)$, loaded lazily on first entry into Simulation so that a user who never opens it pays nothing.

The picker lists exactly what is tabulated, which is the quickest way to find out whether a particular ion is present. Not every oxidation state is: the ion lists in these tables are finite, and a highly charged cation such as S⁶⁺ is often absent. That is usually not a problem. A bare S⁶⁺ is a closed-shell Ne core, which describes a covalently bonded sulfate sulfur poorly, and the difference between the ionic and neutral form factors is confined to low $\sin\theta/\lambda$ — well below what a powder pattern distinguishes. Neutral S is standard practice for sulfates. If an ion you want is missing from the preferred table, the other table in scatters/ has a different ion list.

Only X-ray scattering is wired up. The neutron table loads but there is no radiation selector, so neutron simulation is not yet available.

Reading the plot

The simulation draws into the chart's existing Simulation dataset rather than adding a tenth curve — that dataset already means “a pattern computed from the model rather than fitted to the data”, and a second curve with the same meaning could only ever disagree with the first.

With no atoms defined the plot falls back to a stick pattern of reflection positions. That is the honest answer as well as a convenient one: with no atoms, the positions are the whole of what the model predicts.

Above 2000 reflections the Kα tick strip shows only the strongest, and the status line says how many were omitted. The pattern always uses all of them — the cap is on the tick marks only. Past that density the ticks are a solid band that tells the eye nothing, while handing the renderer one layout element per reflection is what turns a large cell over a wide range from slow into unusable.

Right-click resets the view, as everywhere else. With no data file open the simulation restores its own axes, since the usual machinery needs a loaded pattern to work out the vertical extent.

The session log

Entering and leaving Simulation is logged, along with the full model at the moment of entry; so is every atom added, edited, removed or loaded from a CIF, and every pdCIF written. Edits are logged on commit rather than per keystroke — a line for every character typed into a coordinate box is not a record of what you did, it is a record of you typing, and it would bury the entries that matter.

Limitations

The simulation is deliberately a simple one. It is a forward calculation of an ideal powder from an ideal specimen, and the following are not modelled:

What it is good for: checking that a candidate structure produces the pattern you measured, seeing which reflections a particular atom controls, generating clean or noisy test data for other software, and building intuition about how coordinates, occupancies and displacement parameters move intensity around.

Probabilistic Space-Group Determination

Once a Pawley refinement has completed, a small ? button becomes active on the space-group line of the control panel. It scores every space group compatible with the current Laue class against the extracted intensities and returns a ranked list of posterior probabilities. The method follows Markvardsen, David, Johnson & Shankland (2001).

The Principle

A candidate space group $S$ makes exactly one testable statement about a powder pattern: a particular subset $e$ of the reflections is systematically absent, so their true integrated intensities are identically zero. Everything else it leaves free. The question is therefore not “does this space group fit?” but “are the reflections it forbids consistent with zero, given how well the data actually determine them?”

A Pawley refinement is what makes this answerable. Because every intensity is a free least-squares parameter, the refinement returns not only the intensities $\hat{\mathbf{I}}$ but the normal matrix behind them, and therefore their full covariance. That matters enormously: in a powder pattern reflections overlap, and overlapping intensities are strongly correlated. A reflection that should be absent can carry a large apparent intensity purely by borrowing it from a neighbour it cannot be resolved from. Judging each reflection on its own $I/\sigma$ ignores this and can be badly misleading in either direction.

The Statistical Model

The true intensities are given a prior that encodes what the space group claims. Writing $\mathbf{P}_S$ for a diagonal matrix,

A variance of exactly zero forces that intensity to vanish; $\tau_j$ is the local Wilson-like intensity scale for reflections the group allows. The refinement itself contributes $\hat{\mathbf{I}} \mid \mathbf{I} \sim N(\mathbf{I}, \mathbf{Q}^{-1})$, where $\mathbf{Q}$ is the intensity precision matrix. Marginalising over the unknown true intensities collapses these into a single multivariate normal,

so the evidence for $S$ is one density evaluation. Writing $f$ for the set of reflections the group allows and dropping every term that does not depend on $S$:

The quadratic form is the evidence carried by the would-be-absent reflections, evaluated with the full correlation structure rather than reflection by reflection. The log-determinant is an Occam factor: it rewards a group for explaining more absences, and without it the group with no extinction conditions at all would win by construction, since it is never contradicted by anything.

Posterior probabilities follow by normalising over the candidate list.

Reading the Result

Each row of the table is a distinct extinction symbol, not a single space group. Powder diffraction cannot separate groups that predict identical absences - centrosymmetric and non-centrosymmetric pairs, for instance - so those share one probability and are listed together. Clicking any symbol adopts that setting. Note that the setting matters: $Pnma$ and $Pbnm$ are the same group but predict different absences relative to the fixed axes, so they appear separately.

• P(SG | data) - the posterior for that extinction symbol.
• absent - how many of the refined reflections the group forbids. A group forbidding nothing is always in the list as a baseline.
• ⟨I/σ⟩ and worst I/σ - diagnostics over the would-be-absent set, using diagonal errors only. These are the numbers a crystallographer looks at first, but they are not what drives the ranking; the probability uses the full correlation structure.

Requirements and Limitations

Two further points are worth keeping in mind.

• The Laue class is not under test. Only the absences are. The Laue class is fixed by the diffraction symmetry you have already chosen, since it determines the reflection orbits and the multiplicities. Running the test from $Pm\bar{3}m$ can never return a group belonging to Laue class $m\bar{3}$, such as $Pa\bar{3}$. If the Laue class itself is in doubt, run the test from each candidate in turn.
• Exact reflection coincidences are handled, but weaken the diagnostics. In cubic and rhombohedral cells reflections such as $(333)/(511)$ or $(300)/(221)$ fall at identical $2\theta$. Their individual intensities are then not separately determined and no covariance matrix exists. The test is formulated in precision space precisely so that it does not need one - $\mathbf{Q}$ may be singular while $\mathbf{B}_S$ never is - so the probabilities remain valid. The $I/\sigma$ columns fall back to conditional estimates in this case and are flagged as such.

Ab Initio Structure Solution: Charge Flipping

Once a Pawley refinement has extracted a set of individual integrated intensities, Powder 5 can attempt to solve the crystal structure ab initio - with no starting model - using the dual-space charge-flipping algorithm of Oszlányi and Sütő. The result is an electron-density map and a list of atomic positions, obtained purely from the diffraction intensities and the unit cell.

The Principle

A diffraction experiment measures the amplitudes \(|F_{hkl}|\) but not the phases \(\phi_{hkl}\); recovering those phases is the phase problem. Charge flipping solves it by iterating between real space (the electron density \(\rho\)) and reciprocal space (the structure factors \(F\)), enforcing one simple constraint in each:

1. Real space. A correct electron density is non-negative almost everywhere. Every voxel whose density falls below a small positive threshold \(\delta\) has its sign flipped: \(\rho \rightarrow -\rho\). This deliberately wrong perturbation is what drives the phases toward consistency.
2. Reciprocal space. Fourier-transform the flipped density, then reimpose everything that is actually known about the structure factors: the measured amplitudes, the point-group symmetry, and the systematic absences. The calculated phases are kept. Reflections that were never measured are left free - that freedom is what lets the map extend beyond the data.

Transforming back gives an improved density, and the cycle repeats. The threshold is expressed in units of the map's own RMS density \(\sigma\), so \(\delta\) is scale-independent; values around \(0.8\sigma\)–\(1.2\sigma\) work for most problems. The agreement between calculated and observed amplitudes,

\( R = \dfrac{\sum_{c} \left| \sqrt{I^{calc}_{c}} - \sqrt{I^{obs}_{c}} \right|}{\sum_{c} \sqrt{I^{obs}_{c}}} \)

is tracked each iteration. The sum runs over measured peaks \(c\), not over individual reflections: a powder pattern cannot distinguish reflections that overlap in \(2\theta\), so scoring them separately would be scoring noise. \(I^{calc}_{c}\) is the total calculated intensity of everything under peak \(c\). A run that has found the structure typically drops below \(R \approx 0.15\) and then oscillates. Powder 5 keeps the best iterate of the run, not the last.

Space-Group Constraints

Classical charge flipping deliberately works in \(P1\): it throws the symmetry away and lets the algorithm rediscover it, which is a useful safeguard when the space group is uncertain. Powder 5 keeps that mode available, but by default it uses the selected space group inside the iteration. Five separate constraints come out of it, and they are worth distinguishing because they behave very differently.

1. Orbits: symmetry equivalents share one intensity

A Pawley fit reports one intensity per symmetry-unique reflection. Its equivalents \(hR\) must be placed on the \(P1\) grid before a transform is possible, and they are generated from the actual symmetry operators of the selected group. The distinction matters: the Laue class is not the holohedry of the crystal system. \(P4\) has Laue class \(4/m\), so a general reflection has 8 equivalents, not the 16 that \(4/mmm\) would give; the same applies to \(\bar 3\) versus \(\bar 3 m\), \(6/m\) versus \(6/mmm\) and \(m\bar 3\) versus \(m\bar 3 m\). Spreading one measured intensity over twice as many positions as really exist puts half of it onto reflections that are not equivalent at all.

Within an orbit the amplitudes are equal by symmetry: \(|F| = \sqrt{I^{obs}/m}\), where \(m\) is the orbit size. This is a hard constraint, applied every cycle. Only the phases are free.

2. Overlap: distinct reflections may share one peak

The genuinely powder-specific ambiguity is different. Reflections that are not symmetry equivalent can still fall at the same \(2\theta\), either exactly (511/333 in a cubic cell) or accidentally. The pattern gives only their sum. Powder 5 groups such reflections into a cluster and re-partitions the cluster's total intensity between them each cycle in proportion to the current calculated values - the same idea as Le Bail extraction, applied inside the flipping loop. The split is free; the total is fixed.

3. Systematic absences are excluded and held at zero

Reflections forbidden by lattice centring, screw axes or glide planes are removed before the charge-flipping orbits are built, and their reciprocal-grid points are forced to \(F = 0\) every cycle. They are found directly from the operators: \(h\) is extinct if some operator \((R,\mathbf{t})\) satisfies \(hR = h\) with \(h \cdot \mathbf{t}\) non-integral. Leaving them free lets the density break the lattice centring outright, and for an \(F\) or \(I\) lattice that is most of the reciprocal grid.

The raw, signed Pawley area and its ESD are retained long enough to test the selected space group. A forbidden fitted HKL never becomes an ordinary CF reflection. If it overlaps an allowed reflection, however, its corrected intensity remains in the measured cluster total and is re-partitioned only among the allowed orbits; this preserves what the powder peak measured without assigning amplitude to a forbidden reciprocal point.

4. Symmetrisation of the phases

For a real-space operator \(\mathbf{x}' = \mathbf{x}R + \mathbf{t}\), the structure factors of a correctly positioned structure satisfy

\( F(hR) = F(h)\,\exp(-2\pi i\, h \cdot \mathbf{t}) \)

Each cycle the members of an orbit are reduced to a common representative, averaged, and pushed back out. The Symmetry in loop control sets how strongly: \(F \leftarrow (1-\lambda)F + \lambda F_{sym}\). At \(\lambda = 0\) nothing is imposed and the run is classical \(P1\) charge flipping. At \(\lambda = 1\) the projection is exact, which converges faster and has a useful side effect: the solution comes out on a standard origin instead of a random one, so the automatic origin search in the Structure block becomes a verification rather than a necessity.

5. Weak reflections are freed, not forced

The rules above impose \(|F_{\text{obs}}|\) on every ordinary reflection. The Weak Fraction control exempts the selected weakest units: their calculated moduli are kept and their phases are advanced by a constant \(\Delta\varphi = \pi/2\) each cycle instead (Oszlányi & Sütő, 2005). Their measured intensity is then used only for the knowledge that they are weak.

Weak classification chooses those units. The default, Unresolved clusters, ranks each measured peak cluster by its total observed intensity per reciprocal-grid member and assigns the same immutable weak flag to every orbit in a selected cluster. This avoids classifying an arbitrary, dynamically changing overlap partition as if it were an independent measurement. Individual reflections instead ranks each orbit by its own initial target intensity per orbit member. In both modes the classification is made once before iteration and cannot change as overlaps are re-partitioned.

The \(\pi/2\) rotation perturbs the iteration and can help it escape the flat, featureless maps the 2004 algorithm can settle into. The CPU and WebGPU paths apply the same immutable classification and phase operation, including the conjugate mate needed for a real density map.

Intensity scaling

The refined Pawley parameter is a peak height, because the profile functions Powder 5 uses are normalised to unit height, not unit area. What charge flipping needs is the integrated intensity, so each height is multiplied by the area of its own profile (summed over the Kα₁/Kα₂ doublet) before it is used. Skipping that step underweights the high-angle reflections by the full Caglioti broadening, which is roughly a factor of two across a typical scan.

The signed fitted height is preserved during that conversion. In particular, a negative Pawley estimate is not clipped before French–Wilson correction: both the integrated value and its ESD are transformed by the same profile-area factor, so weak and negative observations retain their statistical meaning. If French–Wilson cannot be applied, only then are non-positive amplitudes unavailable to the CF modulus constraint, and the result reports that fact.

The integrated intensity then carries \(I \propto m \cdot LP \cdot |F|^2\). The multiplicity \(m\) is absorbed by spreading the intensity across the \(m\) grid points of the orbit, but the Lorentz-polarisation factor \(LP = (1+\cos^2 2\theta)/(2\sin^2\theta\cos\theta)\) is divided out explicitly. Without it the low-angle reflections are weighted several times too heavily and the map is dominated by the first few peaks.

Workflow & Controls

The controls sit in the left panel; results appear in the Charge Flipping tab of the right panel.

Recommended Settings

The defaults are chosen to solve a typical small inorganic structure from good laboratory data. The table below is a starting point by problem type; the paragraphs after it explain what to change when a run fails.

Choosing the grid

The grid must be large enough to hold every reflection in the fit: an index \(h_{max}\) needs \(N \geq 2h_{max}+2\). If it is too small the summary reports how many reflections did not fit and what \(N\) they would need, and those reflections are simply lost. Beyond that requirement, aim for a voxel of roughly 0.25–0.35 Å (the summary prints the spacing). Finer than ~0.2 Å buys nothing when the data stop at 1 Å resolution and costs \(N^3\log N\): 64³ to 128³ is about an eight-fold increase in run time.

Choosing the threshold

\(\delta/\sigma\) is the single most sensitive parameter. Too low and almost nothing is flipped, so the iteration stalls with \(R\) flat and high; too high and the perturbation is so violent that the phases never settle, giving an \(R\) that jumps around without trend. The working range is 0.8–1.2 and 1.1 is a good default. If \(R\) plateaus above 0.3 with no oscillation, lower it in steps of 0.1; if \(R\) is erratic from the first cycles, raise it. Weak or noisy data generally want a slightly lower threshold.

Iterations and random starts

Charge flipping does not converge gradually - it stays flat and then drops abruptly when the phases lock in, often after several hundred cycles. Judge a run by whether that drop happened, not by the final \(R\). More starts is usually a better investment than more iterations: a run that has not found the solution by ~1000 cycles is unlikely to find it by 3000, whereas a fresh random start might find it immediately. For a difficult problem, 10–20 starts of 1000 cycles beats 2 starts of 10000.

Choosing λ

Use \(\lambda = 1\) when the space group is established (a clean space-group test, a known compound, an unambiguous extinction pattern): it converges in fewer cycles, needs fewer random starts, and delivers the map already on a standard origin. Use \(\lambda = 0\) when the group is in doubt, or as an independent check on a \(\lambda = 1\) solution - if the \(P1\) map symmetrises to a high correlation without having been forced to, the group is right. \(\lambda = 0.5\) is the compromise and the default: it accelerates convergence appreciably while still allowing the density to disagree with a slightly wrong group.

Choosing the overlap tolerance

Set it to roughly half the FWHM of a mid-angle peak. Too small and truly overlapped reflections are treated as independently measured, which imposes intensities the data never determined; too large and reflections that were resolved get to trade intensity freely, discarding real information. 0.05° suits sharp laboratory data; broad peaks or a high-symmetry cell with many coincidences want 0.10–0.15°. The summary reports how many reflections ended up sharing a peak - if that number is close to the total, the tolerance is almost certainly too large.

Choosing the weak fraction

Scan it; do not tune it. The published range is broad - the precursor experiment worked with anything from 10% to 50% of the data inside the resolution sphere - so the setting has a plateau rather than an optimum. Try 0, 0.10, 0.20 and 0.30 and compare across all random starts, not one. If 0.15 and 0.20 give visibly different answers, that difference is trial-to-trial scatter and not the parameter.

With the default cluster classification, the fraction is a fraction of unresolved peak clusters, not a fraction of individual HKLs. Exactly \(\lfloor fN\rfloor\) clusters are selected after ranking their total observed intensity per reciprocal-grid member; every orbit in those clusters receives the weak operation. Consequently the reported number of phase-shifted orbits can exceed the number of selected clusters. Individual mode instead selects exactly \(\lfloor fN\rfloor\) reflection orbits by initial target intensity per orbit member. Neither membership changes during the run.

Reading the Solution

Each run appends a timestamped entry to the Charge Flipping history, so earlier attempts remain available from the run selector. A solution contains:

• Best R. The amplitude agreement of the retained map, summed over measured peaks. Below ~0.15 is promising; above ~0.30 usually means the run has not converged - add random starts, adjust the threshold by ±0.1σ, or see Recommended Settings.
• Electron-density peaks. Fractional coordinates. At \(\lambda = 0\) the map is in \(P1\) with an arbitrary origin and hand, so absolute positions are meaningless - compare interatomic distances. At \(\lambda > 0\) the origin is pinned by the symmetrisation and the coordinates are directly comparable to the space group's Wyckoff positions. Read the Charge column, not Height, to tell a heavy atom from a light one - see the column descriptions for why height gets this backwards.
• Symmetry and absence diagnostics. The summary reports whether full operators or a Laue fallback expanded the data, the requested and effective in-loop symmetry, how many reciprocal points were forced to zero, how many forbidden fitted HKLs were excluded from the CF orbits, how many of those shared a cluster whose intensity was transferred to allowed orbits, and how many isolated ones had raw \(I/\sigma \geq 3\). A large systematic-zero count is normal for a centred lattice. A significant isolated forbidden observation is a direct inconsistency to investigate; a forbidden HKL inside an unresolved allowed peak may instead be an artefact of the Pawley partition.
• Density section. A slice through the map, with selectable axis and depth. The colour bar under the image is an absolute scale taken from the whole map, not from the slice on screen, so slices are directly comparable and an empty one looks empty. The caption gives the slice maximum in units of the map's \(\sigma\); anything under about 3σ is ripple, not an atom. The Stretch each slice to its own range checkbox restores per-slice auto-contrast, which is occasionally useful for tracing weak features but makes noise look like structure - the caption says so whenever it is on.
• 3D isosurface. The same map drawn as a surface of constant density, at a level you set in units of the map's \(\sigma\) — the same units as the section caption, so the two views can be read against each other. Drag to rotate, scroll to zoom. At 1–2σ you will see mostly noise; real atoms usually appear as isolated closed blobs somewhere above 3σ.
  • As solved is the raw charge-flipping density, exactly what the algorithm produced.
  • Symmetry-averaged is the same density after the Structure step has shifted it onto the space-group origin and averaged it over the symmetry operators. It is much cleaner, because averaging \(n\) equivalent copies of the cell suppresses the noise by \(\sqrt{n}\) — but it assumes the space group is right, so judge the solution on the raw map and use this one to read it. It only becomes selectable once the structure has been built.
  • None (atoms only) draws no density at all - just the cell edges and the atom spheres. Once the map has done its job and you trust the solution, the surface is mostly in the way: it hides the spheres it is there to justify, and judging bond geometry or spotting an impossible contact is far easier without it. Selecting it greys out the Level slider, since there is nothing for it to act on.

    This is also the fastest option, not merely the tidiest: the surface extraction is skipped outright rather than being computed and hidden, so switching to it is instant even at 128³. The atom coordinates are those of the symmetry-averaged frame, the same as for Symmetry-averaged - with no map on screen there is nothing for them to be aligned against, so the structure's own frame is the only meaningful choice.

  • Atom positions is a list under the viewer controls, one row per unique site from the Structure block: rank, colour, peak height, multiplicity, and a show/hide toggle. Every site starts hidden, so the density can be judged on its own first; click a row to bring one in, or use hide all / show all. Ten rows are shown at a time and the list scrolls beyond that. The header counts what is visible, and the caption under the viewer says how many positions are drawn and how many sites are hidden. Two buttons sit on the row beneath the header, each under the column it acts on: reset colours on the left, over the swatches, and hide all on the right, over the visibility markers.

    Each site is drawn expanded over the symmetry operators: a site of multiplicity 4 appears four times, once beside each of its four density lobes, and hiding its row removes the whole orbit at once. The Structure table lists only one representative per orbit, so drawing that list directly would leave most of the lobes bare. The sphere radius follows the peak height, so in PbSO₄ the four Pb lobes are obvious against the lighter atoms.

    Site colours can be changed. The swatch in each row is a colour picker: click it and the spheres recolour as you drag, so you can match the picture to what you believe the structure to be. The default palette holds sixteen distinct colours, enough that a solution with more than eight sites does not reuse one, but the defaults carry no chemical meaning - they are assigned in order of peak height and nothing more. That is exactly why recolouring is useful: charge flipping gives you density, not elements (see the note under Building a Structure), so the identification of which lobe is which atom is yours to make, and colour is the natural way to record it. Giving the two sites you read as oxygen the same red, or the heavy atom a colour that stands out from its neighbours, makes a crowded cell legible in a way the rank ordering cannot.

    The colour belongs to the site, not to the individual sphere, so every symmetry image of that atom follows at once - the same grouping the show/hide toggle uses. Clicking the swatch opens the picker without hiding the row, so recolouring and hiding stay independent. Reset colours puts every site back to its palette default and is available only once something has actually been changed.

    Nothing has been least-squares refined at this stage; these are located peaks. Seeing them sit centred inside the density lobes is the quickest visual confirmation that the solution hangs together, and bringing them in one at a time is the easiest way to tell which lobe belongs to which site when the cell is crowded.

  The three cell edges meeting at the origin are coloured and labelled a / b / c, so you can tell which axis you are looking down.

The two buttons beside the section viewer write the peak list (.csv) and the full density grid (.grd); the refined sites are exported separately, from the button under the Structure block.

3D Center of Mass Peak Picking

When dealing with heavy scatterers (such as Pb or U) and truncated high-angle diffraction data, Fourier series termination artifacts frequently cause the electron density map to form hollow "donuts" or rings rather than ideal 3D Gaussian spheres.

To prevent standard 1D parabolic interpolation from erroneously placing atomic coordinates on the rim of these donuts, Powder 5 evaluates candidate peak locations using a 3D Center of Mass (CoM) calculation over a $5 \times 5 \times 5$ voxel window ($\Delta h, \Delta k, \Delta l \in [-2, 2]$):

By weighting positive density across the entire $5 \times 5 \times 5$ neighborhood, the algorithm geometrically averages out the ring artifact and places the peak coordinate directly in the true center of the atomic site.

Building a Structure in the Space Group

This runs automatically at the end of every charge-flipping run, and again whenever a run is re-selected from the history (the result is cached, so re-selecting is instant). It takes the actual space-group symmetry operators, read from the bundled database, and:

1. finds the origin shift that makes the map obey those operators - computed in closed form from a single Fourier transform, then validated by symmetrising the map and measuring the correlation;
2. averages the density over the symmetry operators;
3. reduces the peaks to one representative per symmetry orbit, reporting the site multiplicity of each (a peak that maps onto itself sits on a special position).

The symmetry correlation reported alongside the sites is the number to trust: a high value means the density is genuinely consistent with the chosen space group; a low one means either the run has not converged or the space group is wrong. The resulting sites can be exported as a .cif file.

Result Fields Explained

The Charge Flipping results tab reports the quantities below. The Wyckoff search has its own results tab, with its own history selector and its own PDF and CIF exports; it reports an assignment, a free-coordinate count and \(wR(F^2)\), and it has no origin, hand or symmetry correlation to report because the Wyckoff operators fix those by construction.

Peak and site table columns

Requirements and Limitations

• Data quality. The solution is only as good as the extracted intensities. Severe peak overlap, a poor background, or an incorrect cell will all degrade or prevent a solution.
• Resolution. Atomic-resolution data (roughly \(d \lesssim 1.2\) Å) is normally needed to resolve individual atoms; lower resolution gives a blurred map.
• The space group is an assumption. With \(\lambda > 0\) it is imposed, not tested. A wrong group can still produce a low \(R\), because the constraint makes the calculated intensities agree with the observed ones by construction. Cross-check with a \(\lambda = 0\) run and with the symmetry correlation from the Structure block.
• Non-uniqueness. Random starts can land on different origins or the inverted structure. This is expected at \(\lambda = 0\); the structure builder resolves origin and hand against the space-group symmetry.
• Light atoms near heavy ones. As with any Fourier method, light atoms sitting close to much heavier ones may be lost in truncation ripples.
• Overlap is modelled, not solved. Re-partitioning the intensity of a cluster in proportion to the calculated values is a reasonable guess, not a measurement. Structures whose solution depends on resolving a specific pair of overlapped reflections may be beyond reach from powder data alone.

Direct-Space Structure Solution: Wyckoff Search

The Wyckoff search has its own tab in the left panel and its own results tab on the right. It is a peer of charge flipping, not a step inside it: it consumes a space group, a cell, a composition and the Pawley intensities, and not one of those comes from a density map. It does not require a charge-flipping run, and it cannot read one. The two methods share only the Pawley intensity list - the same measurement, seen independently.

This is direct-space (global-optimisation) structure solution, the family that includes FOX, EXPO, DASH and TOPAS's simulated annealing. It works where charge flipping struggles on powder data because of parameter counting: PbSO₄ in \(Pnma\) has nine free positional parameters against a few hundred observations, whereas charge flipping solves for thousands of grid amplitudes and phases from the same data. It has no phase problem, no origin ambiguity, no hand ambiguity and no series-termination ripples - and because the composition is enforced by construction, it cannot invent an atom where there is nothing to explain it.

The objective

A WebGPU swarm (replica-exchange MCMC) places atoms on trial Wyckoff assignments that satisfy the stoichiometry exactly. Each candidate is scored by $$f = (1 - wR2) - \text{pen}$$ where wR2 is the scale-only weighted residual between \(|F_{\text{calc}}|^2\) and the observed intensities, taken over groups of reflections the powder pattern cannot separate. An optimal multiplicative scale is eliminated analytically; no additive intensity intercept is fitted, because a constant shift in extracted reflection intensities is not a physical structure factor. The observed intensity of a group is \(\sum m |F_o|^2\) and the calculated one is accumulated the same way, so the comparison is made where the measurement actually exists. See the scoring section for the formula.

Scattering factors are the tabulated \(f(s)\) with an overall temperature factor folded in as \(f(s)\exp(-Bs^2/4)\), pre-computed per (reflection, element) on the host. \(B\) is estimated by a Wilson plot over resolution shells from the extracted intensities themselves. The Pawley intensities do not lack a temperature factor - they contain the real one, because they are a measurement - so the model needs one to reproduce their fall-off.

Penalties and hard constraints

Clash, bond-window and coordination penalties are soft and ramped during the search, which is correct while the swarm still needs to move through crowded configurations. They are not how the answer is decided. A candidate whose final geometry breaks any stated floor, upper bound, exact-count or at-least-count rule is rejected outright, and candidates are walked in score order until one passes; if every candidate fails, that is reported rather than swallowed, because it is a statement about the composition, \(Z\), or the constraints.

The same test is applied after the coordinate refinement. The refinement moves every free coordinate with no geometry term at all, so guarding only its input guards nothing: a refinement that breaks any rule is discarded and the search positions are kept, with the reason recorded.

Resolution ramp - and why the running wR2 is not the final one

Controls

Distance constraints

One rule per line, blank lines and lines starting with # ignored. Element symbols are matched case-insensitively and must appear in the target formula; anything else is inert and is reported as an error rather than skipped. The general form is

ElemA ElemB [count] dmin/dmax

and every field except the two symbols is optional, which is why the same intention can be written several ways. What matters is which of three different mechanisms a line ends up switching on:

Worked examples

Rule and slot costs are shown against the budget of 8 rules and 64 slots. The constraint box totals them live and refuses the run before it starts rather than dropping the surplus.

A complete example: PbSO₄

# anglesite, Pnma, Z = 4
S  O   4  1.40/1.80     # sulfate tetrahedron
Pb O   8  2.00/3.00     # lead coordination shell
O  O      2.40          # keep the tetrahedra from interpenetrating

Two rules and twelve slots, well inside both budgets. The third line costs nothing at all: it is a floor, and floors live in the clash matrix rather than in the rule table. Note what is not constrained - no Pb–S rule, because there is no Pb–S bond and the O–O floor plus the two coordination shells already fix the packing. Constraining a distance you have no chemical reason to constrain is the commonest way to make the right answer unreachable.

Direction is not decoration

A coordination number is tested from its first element outward. S O 4 1.40/1.80 asks every sulfur for four oxygens; O S 4 1.40/1.80 asks every oxygen for four sulfurs, which in a sulfate is false and will be paid for on every single evaluation. Writing both is not one bond counted once - they are two independent demands costing two rules and eight slots, and the box warns when it sees the pair. Put the count on the centre of the polyhedron.

Upper bounds without a count are the exception: they are symmetric statements about nearest neighbours, so they are applied both ways automatically. That is convenient and it is also why one typed line becomes two rules.

\(Z\), the density window, and the scan

\(Z\) is the number of formula units in the cell, and it is the one input to this search that is neither measured nor refined. Get it wrong and everything downstream is wrong in a way no figure of merit will name: the stoichiometry is enforced by construction, so the swarm will dutifully fit the wrong number of atoms and report a residual for it.

Leaving the box blank selects Auto, which does not estimate \(Z\) - it scans. Every plausible value is searched in full and the best structure for each is offered, because a rule of thumb about cell volume is a prior and the intensities are evidence, and the evidence has to be allowed to overrule the prior.

Where the candidate list comes from

The crystallographic density fixes the relationship between \(Z\) and the cell:

\( \rho = \dfrac{Z M}{N_A V} \qquad\Longleftrightarrow\qquad Z = \dfrac{\rho V}{1.66054\, M} \)

with \(M\) the formula mass in g/mol, \(V\) the cell volume in ų and \(\rho\) in g/cm³. Inverting it turns a density window into a \(Z\) range. The candidate list is then cut down twice more before any GPU time is spent:

1. The density window. Two sliders, whole g/cm³, 1 to 22, default 2 to 10. The bounds are rounded inward, so every candidate's density genuinely lies inside what the sliders say - a maximum of 10 will not admit a \(Z\) at 10.3. This is the only step that uses knowledge the program does not have, so it is the one you set.
2. What the space group permits. Pure arithmetic, described below, and in a centred or high-symmetry setting it removes most of the window at a stroke.
3. Whether a legal assignment exists. Each survivor is put through the same assignment enumeration the search itself uses, capped at one solution, so a \(Z\) that cannot be built at all is discarded in milliseconds instead of failing half a second into a run.

What the space group says about \(Z\)

Every atom of one element occupies some set of Wyckoff positions, so the number of them in the cell is a sum of multiplicities, repeats allowed. Let \(g\) be the greatest common divisor of every multiplicity the setting offers. A sum of multiples of \(g\) is a multiple of \(g\), so for each element with subscript \(k_e\) in the formula unit,

\( Z\,k_e \equiv 0 \pmod{g} \qquad\Longrightarrow\qquad Z \equiv 0 \pmod{\;\operatorname{lcm}_e \dfrac{g}{\gcd(g,\,k_e)}\;} \)

Worked through for a pyrochlore A₂B₂O₇ in \(Fd\bar{3}m\): the multiplicities are 8, 16, 32, 48, 96 and 192, so \(g = 8\). The oxygen requires \(7Z \equiv 0 \pmod 8\), and since 7 and 8 are coprime, \(Z\) must be a multiple of 8. A 2–10 g/cm³ window spans \(Z = 3\ldots11\) for a typical pyrochlore cell; the rule removes ten of those eleven before anything is enumerated, and \(Z = 8\) is the only candidate left.

One search per \(Z\), and a list of solutions

Every surviving \(Z\) gets its own complete search - all the restarts, the resolution ramp, the closing quench, the contact filter and the coordinate refinement - and the best structure it produced is kept. The plausibility ranking (volume per atom nearest 16 Å³, with a mild preference for a \(Z\) commensurate with the group order) decides only the order of the searches, so stopping part-way through means the most likely candidates have already been done.

The results tab then opens with a Solutions table: one row per \(Z\), with its density, atoms per cell, assignment, search wR2 and \(wR(F^2)\). Click a row and everything below it - the coordinate tables, the R factor, the CIF export, the contacts, the bond valences, the structure plot and the PDF report - switches to that solution. The run-history label carries both numbers: Z=8/3 means the best of three solutions was \(Z = 8\).

The ranking is by \(wR(F^2)\). A solution whose refinement did not run has no \(wR\), and is placed below every solution that has one and ranked among its peers on search wR2 - the two are different quantities measured against different things, and mixing them on a common scale would invent a comparison that does not exist. Values of \(Z\) that produced nothing are listed with their reasons rather than quietly omitted.

Typing a \(Z\) instead

A typed \(Z\) is honoured exactly as typed and the density window is not consulted. The panel still computes what that \(Z\) implies and shows it as you type, in either order - entering the formula after the \(Z\) works as well as before it:

Z = 8 gives a density of 6.73 g/cm³ (Gd2Ti2O7, 522.23 g/mol, cell 1030.3 ų).

How many sites a structure needs

Worth checking before a run rather than after, because the answer is one division and it decides whether the search can represent your structure at all:

\( \text{sites} \ \ge\ \left\lceil \dfrac{\text{atoms per cell}} {\text{largest multiplicity in the group}} \right\rceil \)

Atoms per cell is the formula unit times \(Z\), counting hydrogen. The largest multiplicity is the general position, which for a symmorphic group is the group order. Take cimetidine, C₁₀H₁₆N₆S with \(Z = 2\) - 66 atoms per cell:

So the same compound is solvable or not depending on the space group, and the cliff is steep: one step down in symmetry doubles the requirement. In \(P2_1/c\) the narrowest assignment cimetidine admits is 5×4e C, 8×4e H, 3×4e N and one sulfur on a 2-fold inversion site - 17 sites. The odd sulfur is what makes 66, which is not divisible by 4, work out at all.

Final coordinate refinement

The search finishes with a least-squares refinement of the free Wyckoff coordinates against \(|F_{\text{obs}}|^2\), with an overall intensity scale eliminated analytically and each atom re-projected onto its symmetry element every step. It uses the same reflection groups, scattering factors, \(B\) value and weighting policy as the search. If at least one complete powder group has a usable ESD, complete groups receive \(1/\sigma^2\) and incomplete groups receive zero weight; if none has a usable ESD, all groups receive unit weight. Positive weights are normalised internally without changing their ratios. If the search and refinement ever disagree on the scattering or weighting model, they are not fitting the same problem and their figures cannot be compared.

Both answers are reported: the search positions and the refined positions, in separate tables, with the second appearing only when the refinement actually moved something. They optimise different quantities, and when they disagree that disagreement is the interesting part. If you intend to run a Rietveld refinement afterwards, the search positions are usually the cleaner starting point.

Structure plot, contacts and bond valence

Everything above this point on the Wyckoff results tab answers one question: how well do these coordinates reproduce the measured intensities. That is not the same question as whether they are chemistry. A structure can reach a respectable \(wR(F^2)\) with an oxygen 1.3 Å from its sulfur, or with a cation whose coordination sphere is a physical impossibility, and no R factor will mention it. The block below the coordinate tables is the independent check: it uses no intensities at all, only the cell, the space group and the atom positions.

What is displayed

• Cell content, volume and calculated density. The formula per cell is built from the multiplicity counted from the symmetry expansion, not read off the Wyckoff label - a site that has drifted onto a special position during refinement generates fewer distinct images, and the counted value is the one the contacts and the mass are actually built from. A density far from what the compound should have usually means the wrong \(Z\), not a wrong coordinate - and if the search ran on \(Z\) = Auto, the other values of \(Z\) are already sitting in the Solutions table above.
• Structure plot. Drag to rotate, wheel to zoom, click an atom to list its contacts. Sites can be hidden individually, the cell frame and labels toggled, boundary images switched off, and the view saved as a PNG. Bonds are drawn between counter-ions closer than 1.4 × the sum of their covalent radii; they are a drawing convention, not part of the analysis.
• Contacts. Every neighbour of the selected site within the cutoff, sorted, tagged with the site it belongs to, with a rule drawn where the distance jumps by more than 0.45 Å. The gap after the first coordination sphere is what a reader is looking for, and a flat list hides it.
• Structure quality. Coordination number, mean bond length, the spread of those lengths, the bond-valence sum against the formal valence, and a verdict.

Coordination number

Taken over the counter-ions - anions around a cation, cations around an anion - and cut at the first clear break in the sorted distances rather than at the cutoff, separately for each neighbouring element. This info is not relevant for molecular structures.

Both of those choices are forced. A single cutoff cannot serve both a sulfate sulfur and a twelve-coordinate lead in the same structure: at 3.5 Å the sulfur is reported as CN 5 or 6, with a mean bond length halfway between the SO₄ tetrahedron and the next oxygen in the lattice and a spread that describes nothing. And an oxygen in that same sulfate sits 1.5 Å from its sulfur and 2.4 Å from the nearest lead, so a break taken over the merged list would stop at the S–O bond and report CN 1 - the sulfate group, not the oxygen's coordination. Per-element breaks give CN 12 for the lead, 4 for the sulfur and 4 for the oxygen, which is the structure as a chemist would describe it.

The spread is the standard deviation of the distances in that shell. On a cation surrounded by one element it is a distortion measure and a large value is worth a look. On an anion coordinated by two different elements it mixes two bond types and will be large for a perfectly good structure; the shell composition is given beside the CN for that reason.

Bond-valence sums

Each contact between a cation and an anion is assigned a valence from its length alone,

\( s_{ij} = \exp\!\left( \dfrac{R_0 - d_{ij}}{b} \right), \qquad V_i = \sum_j s_{ij} \)

with \(b = 0.37\) Å and \(R_0\) tabulated per cation–anion pair and oxidation state. In a correct structure the sum \(V_i\) over a site recovers that ion's formal valence. The sum runs over every counter-ion inside the cutoff rather than over the coordination shell defined above: the exponential already makes distant contacts negligible, so the sum does not need a shell boundary and is better off without one.

Reported as “V of n expected”, with a verdict that is nothing more than the relative deviation:

Which atoms are anions, and in what oxidation state

Decided per structure, not from a fixed list: the candidate anions present are compared by electronegativity, and those within 0.5 units of the most electronegative species in the cell act as the anions. Sulfur is therefore an anion in PbS and a cation in PbSO₄, which a fixed list gets wrong one way or the other every time. Oxidation states are taken from the common state of each element. They are no longer editable, and the net cell charge is no longer reported: both summed a default state over every site, which means something only for a simple ionic solid with one state per element. For a mixed-valence compound, a molecular or covalent structure, or anything with partial occupancy the total was arbitrary, and it marked sound structures as suspect. The per-site bond-valence sums below answer the same question from the geometry that was actually refined, and answer it honestly.

What the search actually maximises: (1 − wR2) − pen

Every candidate structure is scored by a single number, \\( \\text{score} = (1 - wR2) - \\mathrm{pen} \\), and the candidate list is ordered by it.

wR2 is the scale-only weighted residual over the observed group intensities. Eliminating the optimal multiplicative scale gives

\( wR2 = \sqrt{1- \dfrac{\left(\sum_i w_i I_{c,i}I_{o,i}\right)^2} {\left(\sum_i w_i I_{c,i}^2\right) \left(\sum_i w_i I_{o,i}^2\right)}} \)

This is invariant to overall intensity scale but deliberately not to an additive shift. Centering the intensities would silently fit a constant structure-factor intercept, which has no physical place in this model. The CPU check, WebGPU kernel and coordinate refinement use the same scale-only form.

Reflections the powder pattern cannot separate are summed into one group first, so the comparison is made where the measurement exists. Each group is weighted by \( 1/\sigma^2 \) on its total when usable Pawley ESDs exist. If at least one complete group has an ESD, incomplete groups receive zero weight rather than unit or \(1/I\) weight. Only when no group has a usable ESD is the whole problem unweighted. The same policy is used in final coordinate refinement.

pen is the geometry penalty, summed over every atom in the trial cell. It is what stops the diffraction score buying a chemically impossible structure. It has exactly three components:

The total is then multiplied by penScale, which ramps from 0.2 at the start of a run to 4.0 at the end. Geometry is nearly free early, so the swarm can move through crowded configurations that lie between two good basins, and expensive late, so the answer has to be physically legal. Archived bests are stored with the penalty unscaled and re-scored at the current weight whenever they are compared, so a candidate found early is never judged against a rival on a different scale.

When the ranking is withheld

A high fitness is not by itself evidence that a structure is right, and the search now says so when it cannot be. Two properties of the weighted residual make this necessary.

First, wR2 compares a DISTRIBUTION of intensity, and a model can reproduce the few highest-leverage points while every other reflection is wrong by any amount. A search figure and a refinement \( wR(F^2) \) can therefore disagree sharply: they are computed over different weightings of the same data, and when the effective observation count is small the search figure is constrained by very little.

Second, the number of reflections is not the number of observations. Because the residual is weighted, a group counts for its share of the observed denominator, \( w_i I_i^2 \), not merely for being in the list. Those shares are very uneven: \(1/\sigma^2\) rewards sharp low-angle reflections, whose ESDs are small, and all but erases high-angle overlaps a Pawley fit could not partition, whose ESDs can exceed their own values. The search reports the effective number of observations

\( N_{\text{eff}} = 1 / \sum_i p_i^2, \qquad p_i = w_i I_i^2 / \sum_j w_j I_j^2 \)

which equals \(N\) when all shares are equal and 1 when a single group is everything. A list of 35 groups is routinely 5 or 6 effective observations.

Distance constraints count as observations. A window is information about the structure that did not come from the pattern, and a fit that uses it rests on more than the reflections alone — the same reason Rietveld practice has counted restraints in the observation total for decades. A coordination-shell rule asking for \(n\) partners contributes \(n\) pseudo-observations per independent site, not per atom in the cell: the symmetry images of a site are the same distance seen again, and counting the orbit would inflate the total by the multiplicity. They are reported separately from \(N_{\text{eff}}\), because a structure held up by restraints is a different claim from one determined by the data.

To lift a withheld ranking, give the search more to work with: extract more reflections, particularly at high angle where severe overlap leaves \(\sigma/I\) near or above 1 and the weighting correctly discards them; or add distance constraints, which count. Reducing the free-parameter count — a smaller \(Z\), or sites of higher symmetry — works from the other side.

Recommended Refinement Strategy

A sequential and hierarchical refinement strategy is crucial for achieving a stable and physically meaningful solution. Attempting to refine all parameters simultaneously from a poor starting model will likely lead to divergence or convergence to a false minimum.

Phase 1: Initial Model Setup

1. Define the Model: Load data, select the crystal system and space group, and define the refinement range using the $2\theta$ sliders.
2. Set Background Points: The background is automatically estimated upon loading. Adjust the Auto-points slider under the "Background" tab to optimize the density of points so that they reasonably follow the experimental background. You can also edit the points manually or use Ctrl+Click on the chart. This background shape is fixed during refinement.
3. Peak Position Refinement: Using the Le Bail and LM algorithms, refine only the Lattice Parameter(s) and, if necessary, the instrumental Zero Shift. The goal is to align the calculated Bragg positions with the observed peak maxima.

Phase 2: Peak Profile Refinement

4. Isotropic Broadening: Once positions are correct, refine the primary isotropic peak shape parameters (e.g., W, Y, U, X in TCH, or GW, LX in Simple/Split pVoigt). This will account for the dominant size and strain contributions.
5. Asymmetry and Shape: Introduce asymmetry parameters (S/L in TCH, shft/trns in Simple/Split pVoigt) if there is a clear misfit, particularly at low angles. Refine the mixing parameter (eta) if needed.
6. Anisotropic Broadening (TCH only): If systematic misfits remain (e.g., some peaks are consistently broader than the model), introduce the anisotropic Stephens parameters (e.g., S400). Refine only the symmetry-independent terms.

Phase 3: Finalization and Intensity Extraction

7. Global Optimization (Optional): If the LM algorithm converges to a poor solution during the Le Bail steps, switch to Parallel Tempering (PT) for one or more Le Bail cycles to perform a global search. Afterwards, switch back to LM for a final, precise local minimization.
8. Pawley Refinement: With a stable and well-refined model from the Le Bail method, perform a final refinement using the Pawley method. It is generally recommended to use the Levenberg-Marquardt (LM) algorithm for stability, though Parallel Tempering (PT) can also be used (potentially requiring more iterations). This will provide the most statistically robust set of integrated intensities, suitable for subsequent structure solution.

Interpretation of Results

Assessing the quality of a refinement requires both statistical analysis and critical visual inspection of the fit.

Figures of Merit

Standard crystallographic R-factors are provided to quantify the quality of the fit.

• R-pattern ($R_p$): The unweighted residual error based on net intensities, sensitive primarily to the fit of high-intensity reflections. $$R_p = \frac{\sum |(y_{i,obs} - y_{i,bkg}) - (y_{i,calc} - y_{i,bkg})|}{\sum |y_{i,obs} - y_{i,bkg}|} \times 100\%$$
• Weighted R-pattern ($R_{wp}$): The primary figure of merit, weighted by the inverse of the observed gross intensities ($w_i = 1/y_{i,obs}$), which properly accounts for the counting statistics across the entire pattern. $$R_{wp} = \left[ \frac{\sum w_i (y_{i,obs} - y_{i,calc})^2}{\sum w_i y_{i,obs}^2} \right]^{1/2} \times 100\%$$
• Reduced Chi-squared ($\chi^2$, Goodness of Fit): The most statistically rigorous indicator. For a statistically perfect fit where the model correctly describes the data and the weights are accurate, $\chi^2$ should approach 1.0. $$\chi^2 = \frac{1}{N - P} \sum w_i (y_{i,obs} - y_{i,calc})^2 = \left(\frac{R_{wp}}{R_{exp}}\right)^2$$ where $N$ is the number of data points, $P$ is the number of refined parameters, and $R_{exp}$ is the statistically expected minimum $R_{wp}$.

Weighting, and what an ESD means here

Two different weights appear in the program and they are not interchangeable.

• Profile fitting weights each data point by counting statistics, \(w_i = 1/\max(y_i, \text{floor})\), the right model for a measured count.
• Structure fitting - the Wyckoff search and the final coordinate refinement - weights each observation by \(1/\sigma^2\) from the Pawley decomposition. Counting statistics are the wrong model there: the uncertainty on an extracted intensity comes from the decomposition, not from how many photons arrived, and for a heavily overlapped reflection it is far larger than \(\sqrt{I}\). A measured 200 with an ESD of 500 gets weight \(1/200\) under counting statistics and \(1/500^2\) from its real uncertainty - three orders of magnitude of difference in how much a number consistent with anything from −300 to 700 is allowed to say.

Mixed ESD availability has one rule in both Wyckoff stages. If at least one complete unresolved group has a usable ESD, only complete groups are weighted: they receive \(1/\sigma^2\), while incomplete groups receive zero weight. If no group has a usable ESD, all groups receive unit weight. Missing uncertainties are never replaced by \(1/I\), and positive weights are normalised only for numerical conditioning, without changing their relative influence.

The residual reported for a structure is therefore \(wR(F^2)\), the weighted residual actually minimised. An unweighted figure printed next to a weighted refinement describes a fit nobody performed.

How values and uncertainties are printed

A refined value is given in the crystallographic compact form, $a = 5.5139(3)$ Å: the parenthesis is the standard uncertainty in units of the last digit shown, so this reads $5.5139 \\pm 0.0003$.

The value is rounded to the digit its uncertainty reaches, and no further. A cell edge printed as 5.439771388246059 with the su in a separate column asks you to do that rounding yourself, and everything past the fifth decimal there is the binary representation of a double rather than a measurement. The raw su is still given in its own column for anyone who needs it arithmetically.

The su itself is quoted to one significant figure, except where its leading digit is 1, which gets two: rounding 0.00012 to 0.0001 misstates the uncertainty by 17%, and 0.00010 and 0.00019 would be indistinguishable. With a leading digit of 9 the same rounding is worth at most 5%.

Where the su reaches 10 or more the compact form has no last decimal to refer to, and the report falls back to $6 \\pm 2$. An su that large means the parameter is not determined, and the notation should not disguise it.

Parameters constrained by symmetry

In the parameter table, the Fitted column distinguishes three things, not two:

The middle row used to be printed as No with a blank ESD, which understated the fit: it read as though $b$ had been assumed when it was determined precisely as well as $a$. The distinction matters if you quote a cell edge with its uncertainty, and it is invisible in the fit flags, where a slaved parameter and a fixed one look identical.

The same applies to the Stephens anisotropic terms: $S_{040} = S_{400}$ in a hexagonal or tetragonal cell, and $S_{202} = S_{220}$ as well in a cubic one.

χ² well below 1 means the weights are wrong, not that the fit is good

\(\chi^2 \approx 1\) is the right target, but only when \(w_i\) is the real variance of \(y_i\). The profile weight above assumes \(\sigma^2 = y\), which is exact for raw counts and for nothing else. If the intensities on the vertical axis are not counts, the weights are wrong by whatever factor separates the two scales, and \(\chi^2\) is wrong by the same factor.

Three ways a real dataset stops being counts:

• Normalisation. Synchrotron data are routinely divided by a monitor, the ring current or the exposure time, and the result is written in arbitrary units. A pattern whose strongest peak reads 100 “a.u.” would, if those were counts, scatter by \(\pm 10\) at the peak and by about a quarter of its own height in the background. If the trace on screen is visibly smoother than that, the numbers are not counts.
• Averaging. Multi-analyser crystals, strip and area detectors, and repeated scans all average many measurements into each point. That cuts the true noise by \(\sqrt{n}\) while leaving the intensity scale alone, so \(\sqrt{y}\) overstates \(\sigma\) by the same \(\sqrt{n}\).
• Smoothing or rebinning. Both correlate neighbouring points, which removes point-to-point scatter without touching \(y\). This one cannot be undone after the fact - fit the unsmoothed data if you still have it.

What a wrong weighting scale does and does not change.

• Unaffected: the refined parameters. Multiplying every weight by a constant does not move the minimum, so the cell, the profile terms and the extracted intensities are exactly what they would have been.
• Unaffected: \(R_{wp}\) and \(R_p\). Both are ratios in which the scale cancels, which is why they remain usable while \(\chi^2\) does not.
• Unaffected: the ESDs. Powder 5 reports \(\sigma(p) = \sqrt{\chi^2 \, C_{pp}}\), and scaling the covariance by \(\chi^2\) cancels the weighting error exactly. Cell and profile ESDs come out independent of the intensity scale; intensity ESDs scale with the intensities, so relative ESDs are invariant. That is the correct behaviour, not a lucky accident, and it is why a low \(\chi^2\) is not a reason to distrust the numbers in the table.
• Broken: \(\chi^2\) as a figure of merit, and with it \(R_{exp} = R_{wp}/\sqrt{\chi^2}\) and any judgement of the form “\(\chi^2\) is near 1, so the model is adequate”. Until the scale is right, read \(R_{wp}\) and the difference curve instead.

Overlap clusters: which intensities mean anything

Two reflections closer together than their profile width have nearly parallel derivatives. Intensity can move from one to the other with almost no change to the calculated pattern, so the data cannot say how it divides between them: the normal matrix is near-singular in that direction and the diagonal of its inverse blows up. A large ESD there is the fit reporting, correctly, that this reflection's share is undetermined - not a failed fit, and the peak may look perfect on screen.

What is determined is the direction along which the reflections move together, so the report adds an Overlap Clusters section giving the cluster total and its ESD. That total uses the full covariance, \(\operatorname{var} = \sum_{ij}\mathrm{cov}[i][j]\), not the sum of the individual variances: the off-diagonal terms are large and negative, and it is that cancellation which makes the sum well determined. Adding variances in quadrature, as if the reflections were independent, gives a hugely pessimistic total and misses the point.

Kα₂ is included when deciding what overlaps - at high angle the satellite is displaced by several FWHM and lands on unrelated reflections, correlating pairs that look well separated in a list of Kα₁ positions. Members whose \(\sigma\) exceeds their intensity are marked with an asterisk. Use the cluster total, not the individual values, for anything downstream.

Exactly coincident reflections

An overlap cluster is a set of reflections the data separates badly. This is a set it cannot separate at all. In a cubic $F$ cell, $333$ and $511$ have the same $h^2+k^2+l^2$ and therefore the same $2\theta$, the same profile, and identical columns in the normal matrix. So do $700$ and $522$, and a dozen more pairs besides. Every lattice produces them, and accidental coincidences in low-symmetry cells do the same thing.

The linear system is then exactly singular in those directions. Powder 5 factors it with a modified Cholesky that repairs only the pivots that fail, so the several hundred intensities that are determined come back untouched - a single damping term large enough to get past every coincidence would have biased those too. The affected reflections are listed in the report under Exactly Coincident Reflections, with the group total.

Calculating Observed Intensities ($I_{obs}$) for Overlapping Peaks

A simple numerical integration over a fixed angular range is insufficient for accurately determining the observed integrated intensity ($I_{obs}$) of overlapping peaks. This tool employs a more robust intensity partitioning method.

At each point in the diffraction pattern, the net observed intensity ($y_{obs} - y_{bkg}$) is distributed among all contributing Bragg reflections. This distribution is proportional to the value of each peak's calculated profile function (including both Kα1 and Kα2 components, scaled by the refined peak height) at that specific point. By integrating these partitioned "slices" of intensity for each reflection across the entire pattern (using the trapezoidal rule), the method yields a reliable $I_{obs}$ value (reported as integrated area) that correctly deconvolutes contributions from neighboring peaks.

Visual Inspection

Numerical indicators can be misleading. Visual inspection of the difference plot (observed minus calculated) is the most critical step in evaluating the fit.

• The Difference Plot: A successful refinement should yield a difference plot that consists of random, uncorrelated noise centered on zero. The plot is scaled relative to the main pattern for visibility.
• Systematic Residuals: The presence of structured, non-random features in the difference plot (e.g., "M-shaped" residuals around peaks, broad humps where the spline is inadequate, or un-indexed peaks) is a clear indication of systematic errors in the model. These may arise from an incorrect peak shape, unmodeled anisotropy or asymmetry, an inadequate background model (requiring adjustment of spline points), or the presence of an unaccounted-for impurity phase.

Williamson-Hall Size-Strain Analysis

For refinements utilizing the TCH (Size/Strain/Aniso) profile function, the application can automatically perform a Williamson-Hall analysis to extract approximate microstructural information. This method separates the contributions of crystallite size and microstrain to the total peak broadening by analyzing their different dependencies on the diffraction angle, $\theta$.

The analysis is based on the linear Williamson-Hall equation, where $\beta$ is the total physical peak breadth (FWHM) in radians derived from the refined sample-only broadening parameters (U, X, Y), excluding instrumental contributions (V, W): $$\beta \cos(\theta) = \frac{K\lambda}{L} + 4\epsilon \sin(\theta)$$

This equation describes a straight line when plotting $\beta \cos(\theta)$ vs. $4\sin(\theta)$. The software performs a linear least-squares fit on the relevant data points (within the fitted 2θ range) to determine the y-intercept (related to crystallite size, $L$) and the slope (related to microstrain, $\epsilon$).

Two independent estimates are reported. The TCH profile already separates the two effects analytically - $Y$ broadens as $1/\cos\theta$ (size) while $X$ and $U$ broaden as $\tan\theta$ (strain) - so size and strain can be read straight off the refined parameters with no linear fit at all. The report prints these as from Y directly (nm) and from X,U directly (%) alongside the fitted Williamson-Hall values. Close agreement between the two is a good sign; a large discrepancy usually means the broadening is anisotropic and neither isotropic number should be trusted.

Reported Values

• Apparent Crystallite Size (nm): An estimate of the average size of the coherently scattering domains, calculated from the y-intercept of the Williamson-Hall plot (using $K=0.9$).
• Apparent Microstrain (%): An estimate of the root-mean-square strain within the crystallites, calculated from the slope of the plot.
• Linear Fit R²: The coefficient of determination for the linear regression. A value close to 1.0 indicates that the isotropic size/strain model is a good fit for the observed peak broadening. Values significantly less than 1.0 may suggest that broadening is anisotropic or that the model is otherwise inadequate.

Note: This Williamson-Hall analysis provides an approximation based only on the isotropic TCH parameters (U, X, Y). It does not account for anisotropic broadening effects (Stephens parameters) or instrumental contributions (V, W). For rigorous quantitative analysis, dedicated size/strain analysis software should be employed.

Data Export

One report button sits in the top-right corner of the right panel, between the light/dark toggle and the help icon. It reports whichever tab is in front, as a PDF:

A refinement is reported only from its own tab, and only for the run selected in that tab's history selector. The Plot tab reports the pattern and nothing else - otherwise the same button on the same tab would produce a Le Bail report after a Le Bail run and a Pawley one after a Pawley run, with nothing on screen saying which you were about to get.

It is never greyed out. Press it where there is nothing to report and it says what is missing, which is more than a disabled button ever conveyed. This replaces the four separate Generate PDF buttons that used to sit one per results tab, each with its own enabled/disabled rule to keep in step with the run history.

Save report on the Le Bail and Pawley tabs still writes the plain-text version: all statistical indicators, refined parameters with their ESDs, the Williamson-Hall results where applicable, the declared Lorentz–polarisation model, the background spline points, the reflection table below, and a point-by-point list of observed, calculated, background and difference intensities across the fitted range. The PDF carries the same reflection table; rows too wide for the page are set smaller rather than clipped.

The reflection table

Both the text and PDF reports list every reflection inside the fitted range with the complete chain from peak to structure factor. Le Bail and Pawley share the same columns; Pawley adds one.

A legend above the table restates the polarisation model, the formulae and the value of $K$ actually used, so every $|F_o|$ can be reproduced from the printed numbers without knowing which convention the program follows.

References & Further Reading

About This Tool

The powder5 toolkit was developed by Nita Dragoe from Université Paris-Saclay as a simple browser-based implementation of powder pattern decomposition methods. It is a long-time successor of PowderV2 (Dragoe, N. (2001). J. Appl. Cryst., 34, 535) and has been updated to include the Pawley method and modern global optimization algorithms.

All numerical work - the skyline Cholesky behind the covariance matrix and the ESDs, the banded solver for the separable Pawley intensities, and the least-squares machinery itself - is implemented locally. An earlier version depended on math.js for these; it was removed so that a refinement can never be blocked by a content-delivery network, and so the application works fully offline. Charting uses Chart.js; all pan, zoom and pinch behaviour is implemented directly against the chart's scales, with no plugin. PDF generation uses jsPDF and html2canvas. Space-group symmetry data is derived from the Computational Crystallography Toolbox (cctbx).

This document was updated with the assistance of an AI.