comp-Harker: Technical Reference

Program Name: comp-Harker — composition-driven Harker/Patterson structure solution. Earlier versions were called Harko and sHarko; the name changed with the composition-driven route, which is what the program is now built around.
Core Architecture: Asynchronous FFT Patterson Engine + Symmetry Minimum Function (SMF) / Buerger Superposition + WebGPU replica-exchange Monte Carlo
Symmetry Database: 230 space groups in 527 settings, with Wyckoff projectors and coset operators, generated by cctbx_Harko_v1.py into sg/; scattering factors in scatters/
Structure solution: composition-driven — Wyckoff assignment search ranked by Patterson map correlation (section 7)
Version: Initial Harko version 12 Oct 2025; composition-driven solution added August 2026, when the program was renamed comp-Harker.

1. Theoretical Background & The Phase Problem

In X-ray crystallographic analyses, experimentally measured diffraction intensities $I(h,k,l)$ are directly proportional to the squared structure factor magnitudes $|F(\mathbf{h})|^2$. However, phase information $\phi(\mathbf{h})$ is inherently lost during measurement. comp-Harker bypasses the direct phase requirement by synthesizing the 3D Patterson function, which represents the real-space autocorrelation of crystal electron density:

$$ P(\mathbf{u}) = \frac{1}{V} \sum_{\mathbf{h}} |F_{\mathbf{h}}|^2 e^{2\pi i \mathbf{h}\cdot\mathbf{u}} = \frac{1}{V} \sum_{\mathbf{h}} |F_{\mathbf{h}}|^2 \cos(2\pi \mathbf{h}\cdot\mathbf{u}) $$

where $\mathbf{u} = (u,v,w)$ are fractional Patterson coordinates and $V$ represents the real-space unit cell volume.

Physical Properties of the Patterson Map

  • Interatomic Vector Representation: Patterson map peaks do not directly correspond to atomic locations $\mathbf{x}_i$. Instead, each peak located at $\mathbf{u}$ represents an interatomic displacement vector $\mathbf{u}_{ij} = \mathbf{x}_i - \mathbf{x}_j$ between atoms $i$ and $j$.
  • Weighting Factor: Peak heights scale proportionally with the product of the atomic numbers of the interacting pair ($Z_i Z_j$). As a consequence, heavy atoms generate dominant peaks that stand out clearly against lighter organic matrices.
  • Origin Peak: The point at $(0,0,0)$ corresponds to all self-vectors ($\mathbf{x}_i - \mathbf{x}_i = \mathbf{0}$) with total integrated power $\sum Z_k^2$. comp-Harker automatically masks a spherical volume around the origin to prevent scalar distortion in peak finding and colour map rendering.
  • Centrosymmetry: Regardless of whether the space group is non-centrosymmetric, the Patterson function is strictly centrosymmetric ($P(\mathbf{u}) = P(-\mathbf{u})$).

Harker Sections

Symmetry elements within space groups generate interatomic vectors constrained to specialized planes or lines termed Harker sections. For a general space group operation $\mathbf{S}\mathbf{x} = \mathbf{R}\mathbf{x} + \mathbf{t}$, symmetry-equivalent atoms generate Harker vectors defined by:

$$ \mathbf{u}_H = \mathbf{x} - (\mathbf{R}\mathbf{x} + \mathbf{t}) = (\mathbf{I} - \mathbf{R})\mathbf{x} - \mathbf{t} $$

By evaluating peak maxima on these constrained planes, comp-Harker extracts discrete fractional atomic coordinates directly from one-dimensional and two-dimensional searches.

2. Program Architecture

comp-Harker isolates computational workloads across thread boundaries to guarantee smooth UI interaction at high animation frame rates.

Multi-Threaded System Map

Execution Layer Primary Responsibility Source Components
Main Thread (UI) File parsing, UI state management, Three.js 3D structure rendering, 2D slice visualizer, PDF report generation. Harko.html, patterson3d.js, style.css, sg_engine.js
Web Worker Engine Asynchronous 3D Radix-2 FFT synthesis, Lorch filtering, peak finding, SMF & Buerger minimum superposition map generation, site consolidation. sharko_worker.js, symmetry_utils.js
WebGPU Compute Shader Hardware-accelerated structure search: the fitness evaluator for every trial structure, GPU symmetry expansion, minimum contact evaluations. The sampler itself is host-side. swarm_compute.wgsl

2.2 Web Worker Pipeline

When data files are uploaded or parameters change, comp-Harker triggers a background worker thread via a single CALCULATE execution message:

STAGE 1: ASYNCHRONOUS FFT PATTERSON SYNTHESIS
  ├─ Expand unique reflections across full Friedel sphere using space group symmetry operators
  ├─ Apply Lorch modification envelope: I *= (1-s) + s * sinc(π * d* * d_min)
  ├─ Grid determination: N = next_power_of_two( max(requested_grid, 2*h_max + 1) )
  ├─ Populate complex grid array with intensity values I(h,k,l)
  ├─ Execute 3D Radix-2 Cooley-Tukey Inverse FFT
  ├─ Map normalization: P(u,v,w) = Real(FFT) / Volume
  └─ Extract metadata: d_min resolution limit, calculated peak width σ = max(0.26 * d_min, 0.7 * dx)

STAGE 2: PEAK DETECTION, SMF & BUERGER SUPERPOSITION
  ├─ Generate origin mask: exclude all grid points within 1.1 Å radius of (0,0,0)
  ├─ Calculate map statistics (mean, sigma) over non-origin voxels
  ├─ Compute Symmetry Minimum Function (SMF) across all non-identity space group operators
  ├─ Fallback (if P1): Compute Buerger Minimum Superposition Map M(u) = min(P(u), P(u - u_top))
  └─ Extract peaks from SMF/Supermap (unmasked) as consolidated absolute candidate sites
                

2.3 WebGPU Structure Search

The third piece of the pipeline, and the only one that runs on the GPU rather than the main thread or the worker, is the search for atomic coordinates once the map and its peaks exist. Rather than duplicate that material here, this subsection is a pointer: the full treatment – the update equations, the fitness function, and importantly the two separate hardware limits that govern how large a search you can run – lives in section 6, with the hardware-limits half specifically in 6.3. The short version, for the architectural picture: the main thread builds the WGSL compute shader from swarm_compute.wgsl, sizing a couple of its constants to your actual GPU before compiling it, then dispatches one workgroup per particle every generation and reads back only the best-found fitness and position, keeping the heavy per-particle arithmetic entirely off the JavaScript thread.

Since the composition-driven route was added, each particle also carries a Wyckoff assignment, and one dispatch searches every assignment at once as independent sub-populations, each with its own ladder. See 7.7.

2.4 Symmetry & Operator Resolution

Every stage of the pipeline – reflection expansion, Harker section geometry, the SMF, and the search's own symmetry expansion of trial atoms – needs the same thing first: the actual list of symmetry operators for whichever space group was selected. That sounds like a simple lookup, but a single space-group number can correspond to more than one standard setting in the database (different origin choices or axis conventions produce different operator lists for what is nominally "the same" space group), so comp-Harker has to decide which setting you actually mean. It does this by trying to match the symbol your data file supplies against the settings on record for that number; if there is only one setting on file it uses that regardless, and if there are several and none matches the file's symbol, it falls back to the first on record and says so explicitly – because picking the wrong setting here does not cause an obvious failure downstream, it just quietly makes every Harker section, every SMF-derived site, and every search-generated symmetry copy wrong in a way that looks like a normal failed structure solution rather than a mismatched space group. If a structure refuses to solve and you're not certain the setting is right, this resolution step, not the search, is the first thing worth checking.

3. Data Input & Intensity Correction

3.1 File Formats

comp-Harker reads ASCII data exports generated by powder diffraction profile refinement tools (such as Powder 5). It also supports generic peak lists, such as the standard SHELX HKLF 4 format (e.g., lines containing h, k, l, intensity, sigma), but this functionality was not fully tested. If metadata like the Space Group or Cell parameters are found in the header, comp-Harker will auto-populate the UI. Otherwise, they can be set manually.

3.2 Intensity Classification Hierarchy

Synthesis of a physically accurate Patterson map requires pure structure factor magnitudes squared, $|F(\mathbf{h})|^2$. However, raw measured powder diffraction peak intensities ($I_{\text{obs}}$ or $I_{hkl}$) carry geometric and instrumental distortions:

  • Lorentz–Polarization Factor ($Lp$): Instrument geometry and beam polarization artificially inflate measured intensities at low and high scattering angles.
  • Multiplicity ($m$): Overlapping symmetry-equivalent reflections in powder diffraction scale the integrated peak area proportionally to the orbit multiplicity.

If $Lp$ and $m$ are not removed, low-angle reflections artificially dominate the Patterson synthesis, distorting vector peak heights and masking true interatomic vectors.

Because profile refinement software exports vary, comp-Harker inspects the uploaded file columns and automatically applies the highest available priority tier:

Priority Tier Detected Columns Intensity Calculation ($|F|^2$) Processing Route & Description
Tier 1 (Highest) |Fo| $|F|^2 = |F_o|^2$ Direct $|F_o|$ values supplied by the file. The refinement software has already removed $m$ and $Lp$ at high floating-point precision using refined $2\theta$ positions. No further corrections are applied.
Tier 2 I_hkl, m, Lp $|F|^2 = \frac{I_{hkl}}{m \cdot Lp}$ Explicit $m$ and $Lp$ columns present in the table header are used directly to isolate $|F|^2$ per reflection.
Tier 3 I_hkl, 2th $|F|^2 = \frac{I_{hkl}}{Lp(2\theta)}$ $Lp(2\theta)$ is calculated by comp-Harker using the beam polarisation constant $K$ read from the file header. Multiplicity $m$ is handled during symmetry orbit expansion.
Tier 4 (Fallback) I_hkl only $|F|^2 = I_{hkl}$ Neither $|F_o|$ nor $2\theta$/$Lp$ data are available. Uncorrected raw peak areas are used. Low-angle reflections will dominate the map, and a warning is displayed.

4. Patterson Map Synthesis

4.1 FFT Synthesis Engine

comp-Harker uses a 3D Radix-2 Fast Fourier Transform to calculate maps ($O(N^3 \log N)$), accelerating calculations by over 300 times compared to direct summations.

4.2 Lorch Series Termination Filter

Truncation of the reciprocal space data at resolution limit $d_{\min}$ introduces Fourier series termination ripples (ghost oscillations around heavy atom vectors). comp-Harker incorporates a tunable Lorch modification strength slider ($s \in [0, 1]$) applied per-reflection during accumulation:

$$ I(\mathbf{h})_{\text{modified}} = I(\mathbf{h}) \cdot \left[ (1 - s) + s \cdot \text{sinc}(\pi \cdot d^* \cdot d_{\min}) \right] $$

where $d^* = \frac{1}{d(\mathbf{h})}$ is the reciprocal interplanar distance in $\text{\AA}^{-1}$, $d_{\min}$ is the dataset's minimum $d$-spacing, and $\text{sinc}(x) = \frac{\sin(x)}{x}$ (with $\text{sinc}(0) = 1$).

  • $s = 0.00$ (Off): Preserves original measured intensities without filtering. Peak heights remain maximal, though termination ripples may persist.
  • $0 < s < 1.00$ (Tunable Blend): Blends raw intensities with the Lorch envelope, allowing users to dial in just enough ripple suppression without causing excessive peak flattening.
  • $s = 1.00$ (Full Lorch): Applies complete Lorch modification, effectively taming truncation ripples at the cost of slight peak broadening.

4.3 Grid Resolution & Anti-Aliasing

The FFT grid automatically expands to a power of two satisfying $N \ge 2 h_{\max} + 1$ to prevent high-index reflections from aliasing onto low-order data.

4.4 Model Gaussian Smearing

Model Patterson maps generated from trial atomic coordinates are convoluted with a Gaussian broadening kernel matching the resolution limit $d_{\min}$ and Debye-Waller temperature factor $B$.

4.5 Least-Squares Map Scaling

Trial calculated maps are scaled to observed data using least-squares scale fitting over non-origin voxels.

5. Peak Finding & Harker Analysis

5.1 3D Local Maxima Search

comp-Harker identifies peak vectors by scanning non-masked voxels using a 26-neighbor 3D local maximum search with periodic boundary wrapping.

5.2 Harker Section Deconvolution

Peaks located on Harker sections are mapped to candidate fractional atomic coordinates $(x,y,z)$ using space-group-specific analytical expressions.

5.3 Multi-Section Site Combination

Each Harker section is solved independently, so the same real atom typically produces a slightly different-looking candidate position from every section it appears in – one partial solution per section, each only as good as the peak that section happened to find. The worker still solves every section this way, and every one of those partial, per-section positions is kept and shown to you (the "Harker Solutions" column in the Peaks tab, described below), because a partial solution that agrees with the final structure is a useful sanity check even when it wasn't itself used to build that structure.

What those per-section solutions do not currently do is get cross-referenced against each other to build the "Consolidated Sites" list. There is a tolerance-based combination routine in the worker that does exactly that – averaging together whichever per-section positions land within a chosen distance of one another – and the code path for it (a COMBINE_ONLY worker message, and the underlying combineSites() function) is still there, but the slider that used to drive it is hidden in the current interface and nothing in the UI calls it any more. Instead, every entry in Consolidated Sites today comes directly from peak-finding on the SMF or Buerger superposition map described next (5.4) – a single, more reliable source, rather than an agreement vote across several partial, section-by-section solutions. If you are comparing the Harker Solutions column against the Consolidated Sites column and wondering why a position in the first doesn't obviously map onto one in the second, this is why: the second column isn't built from the first.

5.4 Symmetry Minimum Function (SMF) & Buerger Superposition

Where the SMF map shows up now. It no longer has a pane on the Maps tab — reading it as a 2D slice added little to the consolidated site list it produces, and removing it let the three remaining panes grow. Everything below still runs on every Calculate Map: the map feeds the Consolidated column on the Peaks tab whenever Consolidate sites via Harker is off, seeds the search's heavy atom, and is written out by Density (.grd) for inspection in VESTA.

Extracting atomic positions directly from a raw Patterson map is complicated by overlapping vectors and floating relative origins. comp-Harker implements a dual-mode map reduction architecture: Symmetry Minimum Function (SMF) for absolute crystallographic positioning, and Buerger Minimum Superposition as a single-shift fallback.

A. Symmetry Minimum Function (SMF)

The SMF queries the Patterson map at vector displacements predicted by all non-identity symmetry operators $\mathbf{S}_i = (\mathbf{R}_i, \mathbf{t}_i)$ of the space group simultaneously. For any trial atomic position $\mathbf{x} = (x,y,z)$, the expected interatomic Harker vector generated by operator $i$ is $\mathbf{u}_i = \mathbf{x} - (\mathbf{R}_i \mathbf{x} + \mathbf{t}_i)$. The SMF value at $\mathbf{x}$ is the minimum Patterson intensity across all active symmetry operations:

$$ \text{SMF}(\mathbf{x}) = \min_{i = 1 \dots N_{\text{ops}}} P\Big( \mathbf{x} - (\mathbf{R}_i \mathbf{x} + \mathbf{t}_i) \Big) $$

Because all symmetry operators are bound directly to the unit cell's origin, the SMF produces peaks directly in the absolute crystallographic frame (e.g., placing heavy atoms accurately on special positions such as $y=0.75$ in space group Pnma). Taking the multi-way minimum acts as a powerful logical "AND" filter, obliterating background noise and random vector overlaps while preserving genuine atomic sites.

B. Buerger Minimum Superposition Function

In space group $P1$ or when symmetry operators are unavailable, comp-Harker falls back to the classical Buerger Minimum Function. The worker selects the strongest non-origin vector peak $\mathbf{u}_{\text{top}}$ clearing the origin exclusion radius ($d_{\min}$ or $1.1\text{ \AA}$) and shifts a copy of the Patterson map:

$$ M(\mathbf{u}) = \min\left( P(\mathbf{u}), \; P(\mathbf{u} - \mathbf{u}_{\text{top}}) \right) $$

Note: Buerger superposition fixes the origin relative to the shifted atom (placing that atom at $(0,0,0)$), whereas SMF resolves positions within the true space group frame.

C. Unmasked Consolidated Peak Detection

Unlike the raw Patterson map (where the $(0,0,0)$ origin peak is masked to prevent scalar distortion), peak finding on the resulting SMF or Superposition map is performed without an origin mask. Because SMF/Superposition shifts valid atomic sites directly onto grid positions (which may include $(0,0,0)$ or special positions on mirror planes), unmasked peak detection ensures valid structural candidates are never accidentally deleted.

6. WebGPU Structure Search

6.1 Replica-Exchange Monte Carlo

comp-Harker searches for atomic coordinates with a population of independent Markov chains on the GPU. Each chain holds one trial structure. Every step it proposes a Gaussian move on the free coordinates, projects the result back onto its Wyckoff subspace, and accepts by the Metropolis rule — always if the fit improved, and with probability $\exp(\Delta f / T)$ if it did not:

$$ P(\text{accept}) \;=\; \min\!\left[\,1,\; \exp\!\left(\frac{f_{\text{new}} - f_{\text{cur}}}{T}\right)\right] $$

Accepting a worse structure some of the time is the whole point: it is how a chain climbs out of a hole it has fallen into. $T$ decides how readily. The proposal width of each chain is adapted individually toward a $30\%$ acceptance rate, so a twelve-parameter assignment and a six-parameter one each settle at whatever step size their own subspace can take — a single global step size handicaps the flexible assignments, which are the ones that matter.

The temperature ladder

Temperatures are fixed, not cooled. The chains of each assignment are grouped into ladders of eight, holding a geometric sequence from $T = 5\times10^{-4}$ to $0.05$ in correlation units. Every ten steps, adjacent rungs attempt to exchange their structures:

$$ P(\text{swap}) \;=\; \min\!\left[\,1,\; \exp\!\big( (\beta_{\text{cold}} - \beta_{\text{hot}})(f_{\text{hot}} - f_{\text{cold}}) \big)\right], \qquad \beta = 1/T $$

A hot replica that finds something better hands the structure down the ladder; the cold replica's worse structure goes up to be knocked about further. The configurations trade places and the temperatures stay with the rungs, so the population always holds the same spread of behaviour: hot chains roaming for the whole run, cold chains refining. Exchanges are only ever attempted between chains of the same assignment — two assignments are different problems in different subspaces and their coordinates are not interchangeable.

This replaced particle swarm optimisation, and the reason is worth recording. PSO's social term pulls every particle of an assignment toward that assignment's best, so the population collapses onto whichever basin was found first. On real PbSO$_4$ data that showed as a run-to-run coin flip between structures containing a Pb–O contact of $1.01\,\text{\AA}$ — not a worse answer, not an answer at all — while the correct structure, when constraints were supplied, scored best on both the correlation and the $R$ factor. The objective had never been the problem; the sampling had. Independent chains cannot collapse, because nothing shares a direction.

6.1b The Final Quench

What an assignment stores at the end of a run is the best structure ever proposed for it – and a proposal is judged by one chain, at the temperature of the moment, under a partial reflection set. Nothing has gone back and refined that particular structure. So after the restarts, every assignment's best becomes the starting point for a greedy descent: $T = 0$, improvements only, the step shrinking geometrically to well below the precision the coordinates are printed at, at full resolution and full penalty weight.

Path length, not breadth. A greedy descent walks down a valley one move at a time, so many short chains and one long chain are not the same purchase. Measured on real PbSO$_4$ data, starting from a structure at $R = 11.9\%$: $23\times150$ steps reached $10.5\%$, $23\times600$ reached $10.0\%$, and a single chain of $10{,}000$ steps reached $9.9\%$ on a third of the evaluations. Chains only try different directions from the same point; they cannot walk further for you.

Every assignment is quenched, not only those that survived pruning. The candidate table ranks on $R$, and refining some structures and not others would put "did it survive the prune" into that ranking — the same bias against high-dimensional assignments that the particle weighting and the late prune point exist to remove.

6.2 Vector Fitness & Sigma Normalization

Superseded. The vector fitness below is no longer what the search maximises. It is retained here because it explains the quantity the older sigma-scaled readouts reported, and because it remains a good description of why a Patterson vector sum makes a smooth funnel. The current fitness is the correlation between the observed and calculated Patterson maps — dimensionless, bounded by 1, so the $\sigma$ conversion does not apply and the chart is labelled in CC. See 7.4 for what replaced it and why.

comp-Harker evaluates trial structural candidate configurations by measuring how strongly the observed Patterson map supports the predicted interatomic vectors.

A. Raw Vector Fitness Calculation

For $N$ independent atoms expanded across all space group symmetry operations, the raw fitness $\text{Fitness}_{\text{raw}}$ computes the pair-weight-normalized mean Patterson intensity at all predicted vector positions $\mathbf{u}_{ij} = \mathbf{x}_i - \mathbf{x}_j$, minus steric collision penalties:

$$ \text{Fitness}_{\text{raw}} = \frac{1}{\sum_{i < j} 2 Z_i Z_j} \left( \sum_{i < j} 2 Z_i Z_j \cdot P(\mathbf{u}_{ij}) \;-\; \text{Penalty}_{\text{bump}} \right) $$

where $P(\mathbf{u}_{ij})$ is trilinearly interpolated from the observed map. Normalizing by the constant sum of all pair weights ensures the search is not artificially biased toward structures sitting on special positions.

B. Sigma ($\sigma$) Scale Conversion

Raw map values depend entirely on the scale of the source file intensities. To provide a standard metric, comp-Harker converts the raw fitness into standard deviation units ($\sigma$):

$$ \text{Fitness}_{\sigma} = \frac{\text{Fitness}_{\text{raw}}}{\sigma_{\text{map}}} $$

where $\sigma_{\text{map}}$ is the standard deviation of non-origin voxels in the observed map. Quoting the fitness in $\sigma$ places the convergence score on the exact same physical scale used by the 2D color bars, peak table height values, and PDF reports.

6.3 Restarts & GPU Hardware Limits

A chain that has cooled is finished — at the bottom rung it only accepts improvements, so it stays in whatever basin it is in. A restart throws every chain back to a fresh random position and re-runs the ladder, keeping every structure already found. Restarts and Iterations spend the same budget differently: restarts buy more independent attempts, iterations buy a longer look at each one. Underneath that, comp-Harker has to fit two very different things inside the fixed budget every GPU gives a compute shader, and it now measures your actual hardware for both rather than assuming a number that might be too small for a good GPU or too large for a modest one.

A. How many particles can run at once

Every chain (one trial structure) needs its own slice of three large storage buffers holding the fractional coordinates of every atom in the asymmetric unit. comp-Harker reads the browser's maxStorageBufferBindingSize and maxComputeWorkgroupsPerDimension limits the moment WebGPU becomes available, works out how many particles' worth of buffer space that leaves after reserving room for the map textures, and sets the Chains slider's maximum accordingly – rounded down to a multiple of 64, never below 64, never above 8192. Because the buffers scale with the number of atoms in your asymmetric unit, adding another atom to the unit shrinks the particle ceiling; you can see the current number, and the raw binding-size and workgroup limits it was computed from, by clicking the ? beside the Chains slider in section 4 of the left panel.

B. How many symmetry-expanded atoms one particle can hold

Separately from the particle count, evaluating a single particle's fitness means generating every symmetry-equivalent copy of every asymmetric-unit atom – all of them, for every operator in the space group – and holding all of those positions in the compute shader's per-workgroup memory at once so every pair of them can be compared. That per-workgroup memory is a small, fixed pool (WebGPU guarantees at least 16 KB per workgroup, and most desktop GPUs offer more), so there is a hard ceiling on (atoms in the asymmetric unit) × (symmetry operators in the space group) that any one particle can carry.

This used to be a single number wired into the shader by hand. It now measures the connected GPU's maxComputeWorkgroupStorageSize at startup, works out the largest atom×operator total that fits with headroom to spare, rounds that down to a whole multiple of the space group's operator count (so an atom's symmetry orbit is never split – carrying five of an atom's six symmetry images and dropping the sixth would silently bias that atom's contribution to the fitness), and compiles that number directly into the shader before the run starts. The floor is 384 combinations, which is generous for most structures – eight independent atoms in a 48-operator space group, or two atoms in one of the 192-operator cubic groups – and capable GPUs comfortably clear a thousand or more.

If your asymmetric unit and space group between them ask for more symmetry-expanded positions than even this device-sized limit allows, comp-Harker does not silently drop the excess atoms and quietly hand back a fitness computed on an incomplete structure. It refuses to start the run at all and reports the exact numbers involved – how many atoms, how many operators, the total that implies, and the GPU's actual limit – so the fix (remove an atom, or accept that this particular structure needs a beefier GPU) is obvious rather than something you'd have to reverse-engineer from a suspiciously bad-looking result.

6.4 Minimum Contact Distance

The Minimum contact distance control is a hard floor under every interatomic contact in the cell, symmetry mates and cell-boundary neighbours included. It runs from $0$ to $3\,\text{\AA}$ and starts at $1.00\,\text{\AA}$. The same number is enforced by the GPU kernel during the search and by the filter that drops candidates before ranking, so the table cannot show a structure the search was not charging for.

Why a floor is needed at all. The correlation does not know that atoms are solid. A model with lead and oxygen $0.32\,\text{\AA}$ apart correlated at $0.9428$ and ranked fourth in a real run before any floor existed — not a worse answer than one at $2.6\,\text{\AA}$, but not an answer at all. The original zero floor was there so that independent atoms of the same element could merge onto a special position without penalty; the composition-driven route removed that need, because coset_ops generates exactly the right number of atoms for each Wyckoff position (7.2), so there is nothing to merge.

⚠ This changed, and the old behaviour was worth reporting as a bug. Between those two designs the control was a ceiling on automatic per-pair floors: the search enforced $\min(0.65\times\Sigma r_{\text{cov}},\ \text{slider})$, so raising the slider above a pair's radii value did nothing. Set to $1.10\,\text{\AA}$ on PbSO$_4$ it still returned candidates with O–O at $0.94\,\text{\AA}$, because the radii floor for that pair is $0.86$ and $\min()$ kept it. A control named "minimum contact distance" that cannot raise a minimum contact distance is misnamed, so the rule is now simply $\text{floor}(A,B) = \text{slider}$, for every pair.

Choosing a value. $1.00\,\text{\AA}$ sits below any real bond — the shortest in ordinary inorganics is a little over $1.1\,\text{\AA}$ — and well above what a wrong Wyckoff assignment produces, so it discards nonsense without putting a true structure out of reach. $0$ turns it off. Raise it with care: too high does not merely reject the right answer, it charges it more than a wrong one, because the clash penalty grows with overlap depth, and the search will spend itself avoiding the region the answer lives in.

Pair-specific distances belong in Distance constraints (7.8). A stated dmin is a hard floor for that pair alone — S O - 1.30/ holds S–O at $1.30\,\text{\AA}$ while everything else keeps the global floor. Constraints are evaluated by the rule machinery, not folded into the global floor, so the same approach is never charged twice.

6.5 Sampler Constants

⚠ Not exposed in the interface. The panel offers only the three controls that change the answer in a way you can reason about — Chains, Iterations and Restarts. Everything below is fixed in SW_DEFAULTS in swarm_wyckoff.js. It is documented because the values govern how the search moves, not because there is a slider for them.
  • Ladder: $8$ rungs, geometric from $T_{\text{cold}} = 5\times10^{-4}$ to $T_{\text{hot}} = 0.05$, exchange attempted every $10$ steps. The units are correlation units, so $0.05$ means a proposal a twentieth worse is taken about a third of the time, and $5\times10^{-4}$ means only a thousandth of a correlation unit is negotiable.
  • Step adaptation: Robbins–Monro toward a target acceptance of $0.3$, with gain $0.08$, clamped to $[2\times10^{-4},\, 0.5]$ in fractional coordinates. Each rung self-tunes: hot rungs settle on large moves, cold rungs on small ones.
  • Resolution ramp: quantised into $24$ blocks. This one is not cosmetic. A Metropolis test compares a proposal against the current state's energy, and both have to come from the same objective; adding a reflection group every generation silently lowers everything's correlation, so a state measured one generation ago looks too good and almost every proposal is refused. Held over a block, with the current state re-measured when the block changes, the comparison is honest at one extra dispatch per block.
  • Quench: $600$ greedy steps, move size $0.02 \to 5\times10^{-5}$.

Reading the acceptance rate

The Acceptance row in the search panel, and the restart lines in the log, report the fraction of proposals the chains took. Healthy is $20$–$35\%$ across the ladder, with exchange between rungs at $20$–$45\%$. Single-digit acceptance early means the chains are frozen and the budget is being wasted; above $90\%$ means they are random-walking and nothing will settle. An exchange rate near zero means the rungs are too far apart in temperature and there are too few of them.

7. Composition-Driven Solution (Wyckoff Search)

This is how comp-Harker now solves structures. You state a formula, and the program decides how many atoms there are, which Wyckoff positions hold them, and where. Run Search requires a formula and will not start without one.

Section 6 still describes the machinery — the particle dynamics, the hardware limits, the contact rules and the tuning parameters all apply unchanged. What has been replaced is what is optimised and what the parameters are: the fitness is now a map correlation rather than a vector sum in $\sigma$ (7.4), and the free parameters are the coordinates a Wyckoff assignment leaves free rather than three per independent atom. The list of independent atoms with quantities that older versions asked for is gone: it stated the one thing this search exists to work out.

7.1 Searching Wyckoff Space

A general-position search optimises $3N$ free coordinates for $N$ independent atoms. Symmetry then generates whatever it generates: in $Pnma$ five independent atoms become forty in the cell, so a model built from one Pb, one S and three O is $\text{Pb}_8\text{S}_8\text{O}_{24}$, not $\text{PbSO}_4$. Multiplicity is measured afterwards, by testing whether symmetry images happen to coincide, rather than imposed.

The composition-driven route inverts this. A trial structure is described by

$$ (\text{Wyckoff assignment}) \;\times\; (\text{free parameters of that assignment}) $$

where an assignment chooses, for each element, a multiset of Wyckoff positions whose multiplicities sum to the required atom count. Two consequences follow:

  • Stoichiometry is exact by construction. Every structure the search can express has the requested formula. Nothing has to be rescued afterwards.
  • Dimensionality collapses. A site on $4c$ in $Pnma$ is $(x, \tfrac{1}{4}, z)$ and carries two parameters; a site on $4a$ is $(0,0,0)$ and carries none. $\text{PbSO}_4$ goes from 15 loose coordinates to 11 constrained ones; $\text{Y}_2\text{O}_3$ in $Ia\bar{3}$ goes from 15 to 4.
The projection trick. Each Wyckoff position carries an idempotent affine projector $\mathbf{x} \mapsto \mathbf{P}\mathbf{x} + \mathbf{t}$ taken from cctbx. Applying it is what keeps a trial structure on its subspace, and under the Monte Carlo sampler it is the only thing that does. A particle swarm stayed on the subspace largely for free, because its update is a linear combination of points already on it; a Gaussian proposal is not, and lands off the subspace almost surely. So every proposal is projected before it is evaluated — not as a correction for drift, but as part of how a move is made. A move is therefore a step in the ambient cell followed by a projection, which is a symmetric proposal on the subspace and leaves the Metropolis acceptance rule valid.

7.2 The Space-Group Database

The generator cctbx_Harko_v1.py writes a per-space-group database into sg/. Beyond the symbol, operators and Harker geometry, each Wyckoff position now carries machine-usable fields:

FieldMeaning
P_num, P_den, T_num, T_den The special operator as exact rationals — the projector above.
n_free $\text{rank}(\mathbf{P})$: the true number of continuous parameters. Not recoverable from the special_op string, which cctbx may print as 1/3x+1/3y+1/3z, ... — three symbols for a rank-1 subspace.
coset_ops Operator indices giving exactly multiplicity distinct images. Since a Wyckoff position has one site-symmetry group throughout, these depend on the position and not on where along it the atom sits.

coset_ops removes the shader's distance-based special-position collapse entirely: the correct number of atoms is generated by construction, with no tolerance to tune and no class of bug where multiplicity depends on a threshold. All 527 settings verify exact.

Three fetch routes are tried in order: sg/all.json (one request), sg/index.json plus the per-group files in parallel, then the legacy monolithic file. The console reports which won and how many settings carry projections — worth checking, because a stale legacy file silently outranking a fresh sg/ shows up only as every Wyckoff position appearing to lack n_free.

7.3 Assignment Enumeration & Ranking

Enumeration is an integer-partition problem per element, with one cross-element rule: a position with no free parameter is a single point in the cell and cannot be shared. In $Pnma$, $\text{PbSO}_4$ at $Z=4$ yields 35 assignments; $\text{Y}_2\text{O}_3$ in $Ia\bar{3}$ at $Z=16$ yields 17. Both include the known answer.

$Z$ may be left blank. The estimate spans 11–22 Å$^3$ per non-hydrogen atom and flags values dividing the group order — a range, not a determination. The familiar 18 Å$^3$ rule is organic-chemistry lore and returns $Z=3$ for $\text{PbSO}_4$, whose $Z$ is 4.

Assignments are ordered before the search by Harker residual: the distance from the nearest consolidated Harker site to the subspace an assignment gives the heavy atom. This is free evidence from work already done. It only separates assignments that place the heavy atom differently, so when every candidate agrees on it the column reads the same value throughout — expected, and stated explicitly in the results header rather than left as a row of identical zeros. Only the correlation ranks the final list.

7.4 The Map-Correlation Fitness

Candidates are ranked by the agreement between the observed and calculated Patterson maps. That comparison is not made in real space, for a reason worth stating carefully.

By Parseval's theorem, the correlation of two Patterson maps equals the cosine similarity of their intensity lists — verified numerically to six decimal places. So the calculated map never has to be built. More importantly, a flat intensity $I(\mathbf{h}) = c$ transforms to a delta at the origin, so subtracting the mean intensity removes the Patterson origin peak exactly — which is precisely what a Pearson correlation does. No real-space mask can achieve that. Measured on a scrambled structure:

ComparisonCorrect structureScrambled structure
CC of the raw maps1.0000.89 — useless
CC of the maps, origin masked1.0000.47 — mediocre
CC of the intensities (Pearson)1.000−0.01 — correct

The observation is the multiplicity-weighted $|F_o|^2$, since Parseval sums over the full sphere where each unique reflection appears $m$ times. Weighting by $m \cdot Lp$ instead — comparing raw peak areas — gives a higher number for the true structure (0.99 against 0.96) and worse discrimination (3.2$\sigma$ against 5.9$\sigma$ over 200 random trials), because $Lp$ spans a factor of 47 across a pattern and a handful of strong low-angle lines come to dominate both.

Reflections a powder pattern cannot resolve are summed into groups and compared as one observation, since Pawley's division of intensity between overlapped lines is a fitting artefact rather than a measurement. Groups are ordered by $d^*$, which makes the resolution ramp a single bound on the group index: the search starts on the low-order quarter of the data, where the $|F|^2$ landscape is smoothest, and reaches the full set at 60% of the run.

On the real PbSO4 data in this repository: the true anglesite structure scores CC = 0.966 against a mean of 0.300 for random structures and 0.475 for the best of thirty — a separation of 5.9$\sigma$.

7.5 Observations, Wavelength & Wilson B

Everything reduces to $|F_o|^2$ and the powder multiplicity $m$ per reflection. How they are obtained depends on the file, and the routes are not interchangeable — applying a correction twice is as wrong as not applying it, and both look plausible.

File providesRoute
$|F_o|$ columnUsed directly; $m$ from the file or the Laue group.
$m$, $Lp$, $I$$|F_o|^2 = I / (m \cdot Lp)$.
$2\theta$, $I$$Lp$ from $2\theta$, $m$ from the Laue group.
$d$, $I$ (+ wavelength)$2\theta$ from $d$ and $\lambda$, then as above.
$d$, $I$, no wavelengthRefused, naming the UI field.

The reconstruction is exact where it can be checked: against a file stating both columns, all 282 multiplicities matched and $Lp$ agreed to 0.067% (the residual being the refined zero shift, which no reconstruction can know). Constant factors — the K$\alpha_2$ doublet ratio, the Pawley scale, the sample volume — cancel entirely, because a correlation is scale-free. Wavelength and Polarisation K are now read from the report header and pushed into the UI.

An overall temperature factor is not optional: on the real data it is worth 0.31 in correlation (0.65 at $B=0$ against 0.97 at the optimum). It is fitted once, before the search, from a Wilson plot:

$$ \ln \frac{\langle |F_o|^2 \rangle}{\sum_j f_j(s)^2} \;=\; \ln K \;-\; 2B\left(\frac{s}{2}\right)^{2} $$

This is structure-independent, so $B$ is fixed thereafter and written into the $B$ slider so the displayed difference map describes the model the search actually scored. A negative slope is reported as unphysical rather than passed through: intensity rising with angle almost always means $Lp$ or a sharpening was applied twice.

7.6 Scattering Factors

cctbx_scatterers_v1.py writes Waasmaier–Kirfel and International Tables coefficients, plus neutron scattering lengths, into scatters/ — 206 entries covering neutral atoms and common ions. The resolution fall-off belongs on the model:

$$ f(s_{tol}) = c + \sum_i a_i \exp\left(-b_i\, s_{tol}^2\right), \qquad s_{tol} = \frac{\sin\theta}{\lambda} = \frac{1}{2d} $$
Two traps, both guarded. First, $s_{tol}$ is half of $s = 1/d$; passing $s$ leaves $f(0)$ untouched, so the error survives a casual check, while making every atom fall off four times too fast — $f(\text{O})$ at $d = 1$ Å becomes 1.38 instead of 2.34. Second, cctbx's table lookups fall back to the neutral atom for an unknown ion rather than failing, so a naive probe of 1128 labels "succeeds" 1128 times and 922 results are the neutral element wearing an ion's name. Both sides now assert $f(0) = Z - \text{charge}$ before the table is used.

An earlier version normalised the observed intensities shell by shell to match point atoms. That is backwards and measurably so: rescaling the data while leaving the model unscaled put the true structure's correlation at 0.88 instead of 1.000.

7.7 Multi-Assignment GPU Search

Every assignment is searched simultaneously in one dispatch, as independent sub-swarms. Particles on different assignments must never attract each other — each keeps its own global best, and chains of different assignments never exchange.

  • Particles are shared by dimensionality, not equally. An assignment with no free parameters is solved by its first particle; an eleven-parameter one is still searching when the run ends and reports a correlation below what it could reach. An equal split therefore ranks assignments by how hard they were to search rather than how well they fit. Weighting by parameter count gives the 11-parameter case roughly 40% more particles than the 7-parameter one, and more than double after pruning.
  • Waves. A wave holds only as many assignments as leaves room to weight. Packing in $\text{particles}/\text{minimum}$ of them spends the whole budget on the floor and the weighting does nothing.
  • Restarts re-seed the particles while keeping every assignment's best, which is how a high-dimensional assignment escapes a bad first draw.
  • Successive halving retires the worst half once, at 45% of the run. A fixed interval fired twelve times in a long run and cut the field to the floor on evidence gathered before the large assignments had converged.
  • Heavy-atom seeding starts half the particles with the heaviest site on a consolidated Harker peak, projected onto that assignment's subspace first. Half rather than all: if the Harker site is wrong, seeding everything would make the assignment unsolvable and hide the error.

A sanity check runs at generation 0. If every particle scores identically the run stops and says so — a kernel that dispatches but writes nothing leaves the fitness buffer at zero, and the resulting list of valid Wyckoff assignments looks entirely reasonable while being nothing but the enumeration order.

7.8 Contacts & Distance Windows

Minimum contact distance is one control used twice: the search is pushed away from anything shorter, and candidates containing a shorter contact are dropped before ranking. It is a hard floor on every element pair, $0$–$3\,\text{\AA}$, default $1.00$ (6.4). The clash penalty scales with how deep the overlap runs — a flat charge meant $0.3\,\text{\AA}$ cost exactly what $2.1\,\text{\AA}$ did, leaving the chains no gradient telling them which way out.

Distance windows are entered one per line, S O 1.35 1.65. The lower bound is a true per-pair constraint, checked on every pair. The upper bound cannot be — most Pb–O pairs in a cell are legitimately far apart — so it is read as a nearest-neighbour condition: every S must have some O within $1.65\,\text{\AA}$. That is the honest reading of the statement, and it is what builds a tetrahedron.

Windows are enforced, not weighed. A model with a $1.8\,\text{\AA}$ O–O contact is not a slightly worse answer; it is not an answer. Candidates breaking a window are removed from the list rather than marked down, and the count is reported. If every candidate is rejected the list is shown unfiltered with a warning, because an empty table says less than a bad one. The penalty weights inside the search are frankly arbitrary and therefore take no part in the ranking.

7.9 Judging a Candidate

Top candidates typically separate by about $0.01$ in correlation, which decides nothing. The Structure quality panel below the coordinate table separates them by a great deal more, and two of its three measures assume nothing at all.

  • Coordination number counts the first shell only — everything within $1.35\times$ the shortest contact. A flat radius counts the four oxygens of a sulfate and whatever lies beyond, reporting a coordination of nine at a mean of $2.7\,\text{\AA}$, which describes nothing.
  • Bond-length spread is the scatter of those first-shell distances. A real sulfate's four S–O bonds agree to $0.014\,\text{\AA}$; the same structure with sulfur displaced gives $0.618\,\text{\AA}$. A factor of forty, where the correlation moved by one part in a hundred. No parameters, no tables — only the fact that coordination polyhedra are regular.
  • Bond-valence sum $\sum_j \exp\!\left((R_0 - d_j)/0.37\right)$, reported against whichever oxidation state fits best, since a site whose sum lands near its own formal charge is evidence in itself. This is the one measure needing tabulated $R_0$; values are quoted from Brese & O'Keeffe (1991) for common cation–oxygen pairs, and unlisted pairs report nothing rather than guessing.
  • R factor on $|F|$ after least-squares scaling, shown in the panel header. It weights differently from the correlation — a weak reflection predicted strongly counts as heavily as the reverse — so models a hundredth apart in CC are often several percent apart in $R$.
Why you need more than one. A candidate correlating at 0.9647 in a real run had three oxygens around its sulfur all at $2.16\,\text{\AA}$: a spread of $0.004\,\text{\AA}$, perfectly regular, and not a sulfate. Its bond-valence sum was 0.89 where 6 is expected. Spread catches misplaced atoms; bond valence catches atoms that are not bonded at all. The verdict column states the reading in words.

Clicking any row in the Candidate structures table loads that structure into the 3D view, the coordinate table, the density readout and the calculated and difference maps. Clicking a row in the coordinate table lists that site's neighbours. Both tables are read-only: the coordinates are a search result and everything else on the page is derived from them.

7.11 Ranking Candidates: $R$ and Hamilton's Test

The search maximises the Patterson map correlation, which is the right thing to maximise and the wrong thing to compare assignments with. That map is dominated by the heavy-atom vectors, so two models agreeing on the lead and differing on the oxygens score within a thousandth of each other. On real PbSO$_4$ data the correct assignment once came fourth at $0.9804$, behind two wrong models at $0.9872$.

The candidate table therefore ranks on $R$ against $|F|$, where every reflection counts rather than every grid point — which is where the light atoms are. The residual is formed on the powder lines, not the reflections: an overlapped group contributed one measured intensity constraining several $|F|$, so the group sums are compared. On this data $282$ reflections are $182$ lines.

The significance column

A model with more free coordinates fits at least as well as one with fewer, so a lower $R$ across assignments of different dimension is not by itself evidence of anything. Hamilton's $R$-ratio test asks the question that matters: is this fit enough better to justify the freedom that bought it. Under the hypothesis that the smaller model is correct,

$$ \frac{wR_{\text{small}}}{wR_{\text{large}}} \;\sim\; \left[\, 1 + \frac{b}{\nu}\,F(b, \nu) \right]^{1/2} $$

with $b$ the number of parameters given up and $\nu = n - m$ the degrees of freedom of the larger model. $n$ is the number of lines, and $m$ counts the free coordinates plus two, for the scale factor and the fitted overall $B$. The column reports $p$, the probability of a difference this large if the smaller model were the true one:

  • $p < 0.05$ — the top row is genuinely better.
  • n.s. — these two are not distinguishable by this data. Two rows with identical correlations are usually this, and saying so is the point: the table stops implying a precision it does not have.
  • equal dim — same parameter count, so the models are not nested and the test does not apply.
  • dominated — fits worse while spending more parameters. No test needed.
⚠ Read $p$ as indicative. Hamilton's distribution theory is exact only when the weights are the inverse covariance of the observations. The weights here are unity, which is the honest default when the file carries no sigmas, but it is an assumption.

7.10 Files & Versions

FileRole
wyckoff_assign.jsEnumeration, projection, particle allocation, formula parsing.
observations.jsIntensity routes, multiplicity reconstruction, Wilson $B$.
scatterers.jsLoads scatters/; asserts $f(0) = Z$.
contacts.jsContacts, coordination, bond valence, $R$ factor.
swarm_wyckoff.jsOrchestration: packing, buffers, the Metropolis / replica-exchange loop, the quench, waves, pruning.
swarm_cc.wgslThe correlation kernel.

Each module declares a version string that Harko.html checks at load and before every run. A mismatch is reported by name, because a cached module throws at line numbers that no longer exist and the obvious reading is that the fix never landed. Script tags carry the version as a query string. Save the modules and Harko.html together.

8. UI & Controls Guided Tour

Welcome to comp-Harker! The UI is designed to keep controls organized on the left while leaving the bulk of your screen open for interactive 2D maps, 3D atomic structures, and peak analysis tables on the right.

💡 Pro-Tip: Resizing Your Workspace
You can drag the vertical divider handle between the left control panel and the right results area to adjust their widths. Every map, chart, and 3D structure automatically resizes on the fly to fit your screen.

Top-Right Quick Controls

Floating cleanly in the upper right corner, you'll find three subtle, lightweight icons:

  • Info Button (ⓘ): Hover your mouse over this icon to see a small tooltip naming the program and the date it was built, without taking up any permanent space on the page.
  • Theme Toggle (☀️ / 🌙): Instantly switches between Light Mode and Dark Mode. Your preference is saved locally, and all 3D cell frames and 2D map canvases dynamically adapt their color palettes.
  • Help Button (?): Opens this documentation page in a fresh browser tab.

8.2 Left Controls Panel

The left panel handles data loading, structural setup, map calculation parameters, the search solver's settings and the run log. It used to be split into two tabs, but a tab hides everything except the one you're on, and several of these settings are ones you want to glance at together – the peak width $\sigma$ slider next to the temperature factor $B$ slider, say, since both broaden the calculated map and you're usually adjusting them in response to the same difference-map residual. So instead the panel is five collapsible sections, stacked top to bottom: 1. Data Input, 2. Structure, 3. Patterson Map Synthesis, 4. Structure Search and 5. Log (the first four open by default, the log closed). Click a section's title bar to fold it away once it is set the way you want; folding is purely to reclaim vertical space and each section remembers its own state.

One way to ask a control anything: click its question mark. The ? beside a caption opens that control's help box; click it again, or click anywhere else, to dismiss it. Nothing in the panel opens on hover. Where a control has state worth reporting – the grid actually used, which intensity route the file took, how many sites combined, the shortest contact found – that line sits at the foot of the same box, under a hairline rule, so the explanation and the current reading are read together. Those readings used to be permanent lines under their controls, which left some controls two lines taller than others and broke the panel's rhythm; the panel is now one uniform column of caption-and-control pairs.

1. Data Input

  • Select Data File: loads a powder export or generic peak list (.txt, .hkl, .dat).
  • Wavelength (Å): used to convert between $2\theta$ and $d$ when the file needs that conversion. Defaults to Cu-K$\alpha$, $1.54056$ Å.
  • Intensity Correction → Divide out Lorentz–polarisation: checked by default, and it only does anything when the file leaves comp-Harker no better option. A file supplying $|F_o|$ directly (Tier 1) or an explicit $Lp$ column (Tier 2) takes priority and the toggle has nothing left to correct; it matters for Tier 3 files that give $2\theta$ but not $Lp$. The route actually taken is reported at the foot of the toggle's help box. See 3.2 for the full priority order.

2. Structure

  • Select Space Group »: opens a searchable modal listing all 230 groups; search by number or symbol (14, P 21/c). The chosen group is named in the line underneath.
  • Unit cell parameters ($a, b, c, \alpha, \beta, \gamma$): filled from the file header where it carries them, editable otherwise.
  • Target Composition → Formula and $Z$: a formula such as PbSO4, with $Z$ optional. This is what makes the composition-driven route available: the formula gives both the elements and how many of each, so the search solves for which Wyckoff positions carry that composition rather than for a fixed list of independent atoms. Leave $Z$ blank and the cell volume suggests a range, reported at the foot of the Target Composition help box.
  • Max sites / element and Max reuse / position: caps on the enumeration — how many distinct Wyckoff positions one element may occupy, and how many independent atoms may share one position at different coordinates. Blank means automatic, which is whatever the composition arithmetically requires; raise them by hand only to allow assignments with more sites than the minimum.
  • Distance constraints: one restraint per line, A B N dmin/dmax — for example S O 4 1.4/1.9 for four oxygens around every sulfur between $1.4$ and $1.9\,\text{\AA}$. dmin is a hard floor, not merely the lower edge of a counting window, and a coordination number is directed. See 7.8.

3. Patterson Map Synthesis

  • High resolution grid: off is a $64^3$ grid, on is $128^3$ – eight times the voxels, memory and time. The FFT needs a power of two, so these are the only two settings exposed; if the data reach an index high enough that $128^3$ would alias one reflection onto another, comp-Harker steps the grid up on its own; the grid actually used is reported in the toggle's help box.
  • Consolidate sites via Harker: chooses which of the two routes from map to atom list is shown and which feeds the search's heavy-atom seeding. Off is the superposition / SMF route, which works from the map itself and needs no symmetry bookkeeping but inherits the map's artefacts; on is the Harker combination, which is sharper when the sections are clean and returns nothing when either tolerance is too tight. Both are always computed, so switching is instant.
  • Harker section tolerance (Å): how close a peak must lie to a Harker plane to be read as coming from it. Too tight and real peaks are discarded; too loose and peaks that merely pass near a plane are solved as if they lay on it. $0$ is automatic: the larger of the peak width $\sigma$ and one and a half grid steps; its help box reports how many peaks were assigned and how many near-missed.
  • Harker combination tolerance (Å): how closely two sections must agree on a shared coordinate before their partial positions are merged into one atom. This slider starts disabled and stays disabled until at least two sections have yielded partial sites — there is nothing to combine before that, and in a group whose sections give fewer than two it stays disabled for the whole session. It re-combines immediately when moved; no map recalculation is needed.
  • Lorch filter strength ($s \in [0,1]$): $s=0$ is off; higher values damp series termination ripples at the cost of broader peaks. See 4.2.
  • Peak width $\sigma$ (Å): Gaussian broadening applied only to the calculated map so it can be compared fairly with the resolution-limited observed one. Seeded to $0.26\,d_{\min}$ whenever a map is calculated.
  • Overall temperature factor $B$ (Å$^2$): a second, independent broadening of the calculated map, combined in quadrature with the peak width. Set it here or turn on Fit B on the Search Results tab. Like $\sigma$, it affects only the calculated and difference maps, never the observed map and never the search fitness.
  • Peak Search Filters — Max Peaks and Min $\sigma$ Threshold: keep at most this many peaks, ignore anything weaker than this many standard deviations above the map mean. Defaults $50$ and $3.0\sigma$.
  • Calculate Map: runs the FFT synthesis, the Lorch filter, the peak search, the Harker section analysis and the SMF / Buerger superposition in one background pass.
⚠ None of the settings above recalculate anything by themselves. Changing the grid toggle, the Lorch slider or the peak filters records the new value and marks the map on screen as out of date – the status line names what changed and the Calculate Map button takes an amber ring. This is deliberate: recalculating on every slider tick meant an accidental drag could kick off a full $128^3$ transform. The two Harker tolerances are the exception in the other direction: they act on peaks already in memory, so they re-derive the site lists as you move them.

4. Structure Search

Four sliders, one group, no sub-headings — they are all parameters of the same search and were previously split under captions that implied otherwise.

  • Minimum contact distance (Å): a hard floor under every interatomic contact, symmetry included — no two atoms may come closer than this, whatever it costs the correlation. Range $0$–$3$, default $1.00\,\text{\AA}$: below any real bond, and well above what a wrong Wyckoff assignment produces. $0$ turns it off. It applies to every element pair alike; to hold one pair further apart than the rest, state a dmin in Distance constraints. Raise it with care — too high does not merely reject the right answer, it charges it more than a wrong one, because the clash penalty grows with overlap depth. See 6.4.
  • Chains: independent trial structures, each its own Markov chain (6.1). The ceiling is measured from your hardware; click the ? to see the limits your GPU reported.
    More chains is not the same as more search. A particle swarm shared a best position, so adding particles added eyes on one basin. These do not talk, so doubling the chains doubles the number of places looked at while each one still looks just as shallowly. Measured on a rugged eleven-parameter test surface at a fixed budget, $2048$ chains of $100$ steps found the answer $6$ times in $10$ while $256$ chains of $800$ steps found it every time.
    But do spend enough chains to collapse the waves. See the callout below — it is worth a factor of four in wall time on a typical run.
  • Iterations: steps taken by each chain. This is the control that matters most: a chain of $100$ steps has barely moved in a structure with ten free coordinates, and no number of chains fixes that. If the fitness chart is still climbing when the run ends, or the acceptance rate is still high, this is the number to raise. Default $400$, up to $4000$.
  • Restarts: independent searches from fresh random positions, keeping every structure already found (6.3).
  • Run Search / Stop: a formula is required. While the search runs, the progress bar, the step counter, the best fitness, the acceptance rate and an elapsed / rate line update in place below the buttons, and a line under those names the GPU doing the work.
⚠ Waves are not a control — they are a symptom of too few chains. An assignment needs a floor of chains to behave as a search rather than a lottery, and the allocator packs only $\lfloor \text{chains} / 48 \rfloor$ assignments into one wave, leaving headroom to weight the rest by free-parameter count. With $512$ chains that is $10$ assignments per wave, so a structure offering $35$ assignments runs in 4 waves — and each wave re-runs the entire restart and iteration budget. Since a run is dominated by the per-generation round trip to the GPU rather than by GPU work, that is close to four times the wall clock for the same number of evaluations. On real PbSO$_4$ data, $512$ chains took $4$ m $30$ s and $2816$ chains took $34$ s for the same search. Raise Chains until the log says "over 1 wave(s)", then spend the rest of the budget on Iterations.
Reading the acceptance rate. Healthy is $20$–$35\%$, with exchange between ladder rungs at $20$–$45\%$. Single digits in the first half of a run means the chains are frozen and the budget is being wasted; above $90\%$ means they are still wandering and nothing will settle. The row turns amber in either case. Both numbers also appear on every restart line in the log.

5. Log

The run log: the intensity route taken, the fitted $B$, the number of assignments enumerated, the generation-0 fitness spread and any warnings. Clear empties it and Save writes it to a text file.

While anything is calculating. The map transform and the structure search both lock the rest of the interface for their duration: loading a new file, recalculating the map, starting a second run and every export button are disabled until the current job finishes. Loading new data out from under a running GPU search would corrupt the buffers it is reading from, so the lock makes that impossible rather than merely unlikely. The exception is Stop, which stays live throughout for the obvious reason; everything else re-enables the instant the run ends or is stopped.

8.3 Right Results Panel

The right side is three tabs: Maps, Peaks and Swarm Results.

1. Maps

A single row of three equal square panes, all sliced at the same axis and section:

  • Observed: the experimental Patterson synthesis computed from the intensities.
  • Calculated: the model Patterson of the structure currently selected.
  • Difference (observed − calculated): missing density and misplaced atoms.
The SMF pane is gone; the SMF map is not. A fourth pane used to show the symmetry minimum function. It was removed because reading it as a slice added little to what the consolidated site list already says, and dropping it lets the remaining three panes be half again as large. The map itself is still computed on every Calculate Map, still drives the Consolidated Sites column on the Peaks tab, still feeds heavy-atom seeding when Consolidate sites via Harker is off, and is still written out by Density (.grd). See 5.4.

Each pane carries its own colour bar with its own minimum and maximum, so the three are individually readable rather than forced onto one scale. Below them: Auto range, the Axis dropdown (⊥ c, ⊥ b or ⊥ a), the Section slider that moves through the cell, and the Density (.grd) export. Auto range left unchecked (the default) scales each pane over its whole volume, so any two sections are directly comparable – a faint slice really does look faint. Checking it rescales to the current slice only, useful when hunting a weak feature, at the cost of making a faint section look as vivid as a strong one.

2. Peaks

Three columns, left to right, so what is derived from what stays visible as one pipeline:

  1. Patterson peaks: raw 3D local maxima $(u,v,w)$, ordered by height in $\sigma$.
  2. Harker solutions: positions $(x,y,z)$ deconvoluted from individual Harker sections – one candidate per section, each only as reliable as the peak that section contributed.
  3. Consolidated: the final candidate list. Which route produced it depends on the Consolidate sites via Harker toggle in section 3 of the left panel: the superposition / SMF map, or positions agreed across sections within the combination tolerance.

Peaks (.csv) at the head of the tab writes all three columns to one file.

3. Search Results

Three columns across the top, then everything the selected structure implies underneath.

  • Fitness evolution: the run's best correlation against cumulative generation, across every wave and restart.
  • Proposed structure: a Three.js view of all symmetry-generated atoms in the cell. Drag to rotate, right-drag to pan, scroll to zoom.
  • Atoms in cell: total atom count, shortest contact and cell mass. The theoretical density is quoted beside the Asymmetric unit heading further down.
  • Fit B: scans $B$ for the value minimising the residual between calculated and observed maps and writes it back to the $B$ slider. Runs immediately if a structure is on screen, and again at the end of every subsequent run.
  • Candidate structures: the ranked list. The ranking is on $R$ against $|F|$, not on the correlation the search maximises — see 7.11 for why. Columns: $R$, the significance of the difference from the top row, the correlation, the search penalty, the shortest contact, whether the distance constraints are met, the Wyckoff assignment, the free-parameter count and, where relevant, the Harker residual. Clicking a row loads that structure into everything below it: the site list, the 3D view, the density, and the calculated and difference maps.
  • Asymmetric unit: the independent sites of the selected structure — element, Wyckoff label, fractional $x, y, z$, and an eye icon that hides a site from the 3D view. The multiplicity had its own column and no longer does: it is the number in the Wyckoff label, so the column said the same thing twice; it survives as the label cell's tooltip. The coordinates are read-only – they are a search result, and the contact list, the calculated map and the density are all derived from them. Clicking a row lists that site's neighbours; a Structure quality panel below reports coordination, bond-length spread, bond valence and $R$ (7.9).
Candidates appear while the search is still running. Every time the run improves on its own best, that structure is kept as a row in the Candidate structures table — correlation, the restart and generation it came from, and its assignment – with the best eight retained. They are clickable like finished candidates, so an intermediate solution that looked right can be gone back to and compared against the final answer instead of being overwritten and lost. Two things to know: a snapshot you select mid-run will be replaced by the next improvement, since the live display owns the site table while the search is going; and the whole live list is replaced by the properly filtered, chemistry-checked ranking the moment the run finishes.

Three actions sit below: Revert edits, Save CIF and PDF Report.

8.4 Exporting CIF, GRD, CSV & Summary Reports

  • Save CIF (Search Results, bottom action row): unit cell, space group and the solved coordinates of the selected candidate, with occupancy $1.0$ and the Wyckoff multiplicity.
  • PDF Report (same row): crystal parameters, the map slices, the fitness chart, the peak and site tables and a structure snapshot. Like Save CIF it needs an actual swarm result before it will click.
  • Density (.grd) (Maps, in the controls row under the panes): writes every populated map as a VESTA-readable volume — observed, calculated, difference and the SMF map, which is still computed even though it no longer has a pane of its own. Panes with no data yet are skipped and named in the status line.
  • Peaks (.csv) (Peaks, at the head of the tab): the Patterson peaks, the Harker solutions and the consolidated sites in one file.

9. Step-by-Step Workflows

The three walkthroughs below cover the same ground at different levels of granularity: 9.1 and 9.2 together are the detailed, one-decision-at-a-time version of what 9.3 compresses into five lines. If this is your first run, follow 9.1 and 9.2; once the workflow is familiar, 9.3 is enough of a reminder.

9.1 Deterministic Harker & SMF Analysis

Everything in this half of the workflow is deterministic – no swarm, no GPU, just the FFT and the peak search.

  1. Load Data: Click Select Data File in the "1. Data Input" section of the left panel to load your .txt, .hkl, or .dat file. comp-Harker will attempt to parse reflections and any available metadata headers.
  2. Define Structure: In "2. Structure", click Select Space Group » to choose your symmetry, and verify the populated cell parameters ($a, b, c, \alpha, \beta, \gamma$).
  3. Calculate Map: In "3. Patterson Map Synthesis", set the Lorch strength, grid resolution, and peak-search filters the way you want them, then click Calculate Map to generate the Patterson synthesis. If you change any of those settings afterward, remember to press Calculate Map again – they no longer trigger a recalculation on their own (see 7.2).
  4. Inspect the Maps: Switch to the Maps tab on the right, choose the slicing Axis and drag the Section slider through the cell. Only the observed pane means anything at this stage – the calculated and difference panes stay blank until a structure exists to calculate from.
  5. Set the Harker tolerances. Leave the section tolerance at $0$ (automatic) for the first pass and read its note: if it reports near misses, widen it. The combination tolerance stays disabled until at least two sections have produced partial sites, which is normal and not a fault.
  6. Check the candidate sites: the Peaks tab, Consolidated column. Flip Consolidate sites via Harker in the left panel to see the other route's answer — both are always computed, so the switch is instant, and agreement between them is worth more than either list alone.

9.2 Composition-Driven Global Search

Picking up from the Consolidated Sites you just inspected, this half hands the search to the GPU. A formula is required: Run Search will not start without one.

  1. State the composition. In 2. Structure, type the formula — PbSO4 — and leave $Z$ blank unless you know it. The line below suggests a range from the cell volume and flags the values that divide the group order. Check it: a wrong $Z$ fails completely rather than partially.
  2. Add what you know about the chemistry. A window such as S O 1.35 1.65 in the Distance windows box tells the search that sulfates exist. Without one the search has no reason to build a polyanion, and a candidate whose sulfur has no oxygen nearer than $2.16\,\text{\AA}$ can still correlate at $0.96$. Set Minimum contact distance a little under the shortest bond you expect.
  3. Set the search size. Chains, Iterations and Restarts in 4. Structure Search. Raise Chains until the log reports a single wave, then put the rest of the budget into Iterations. Assignments with many free parameters need a longer run to be judged fairly against small ones — if a plausible assignment ranks low, raise Iterations before concluding anything. Heavy-atom seeding from the consolidated sites of 9.1 is automatic and needs no toggle: half the particles are seeded, the peak first projected onto the assignment's Wyckoff subspace.
  4. Run. Click Run Search. The fitness curve on the Search Results tab plots the run's best correlation; the generation counter, restart and rate update live on the left. The Candidate structures table fills as the run goes: each improvement is kept as a clickable snapshot, so a structure that looked right at generation 300 is still there at the end. Open 5. Log afterwards for the intensity route, the fitted $B$, the number of assignments and the generation-0 fitness spread.
  5. Judge, do not just read the ranking. The top candidates will sit within about $0.01$ of each other, which decides nothing. Click each row in turn and read the Structure quality panel: bond-length spread, bond valence and $R$ separate them far more sharply than the correlation. A cation whose valence sums far below its formal charge has no real bonds, whatever the map says.
  6. Export. Save CIF writes the selected candidate with occupancy $1.0$ and the Wyckoff multiplicity. PDF Report includes every candidate with full coordinates, not only the one on screen.

9.3 Combined Hybrid Resolution Workflow

  1. Load your diffraction data file and click Calculate Map. comp-Harker computes the FFT Patterson map, applies the optional Lorch filter, runs the peak search and the Harker section analysis, and constructs the Symmetry Minimum Function (SMF) map.
  2. Inspect the observed pane in the Maps tab, section by section.
  3. Review the Consolidated column in the Peaks tab, on both settings of Consolidate sites via Harker.
  4. Enter the formula in Target Composition (section "2. Structure") and any distance constraints you can justify. Then open "4. Structure Search", set the contact floor and the three search-size sliders, and click Run Search.
  5. Work through the candidates in the Search Results tab — the ranking alone will not separate them — then export your CIF file (and, if you'd like, the PDF report right beside it).

10. Troubleshooting & Technical FAQ

Symptom / Question Root Cause / Explanation Resolution
Lorch filter does not change peak list coordinates True interatomic vectors are robust physical features. Lorch modification tames ripples and broadens peaks, but does not displace true vector centers. This is normal and expected. Use Lorch strength ($0.1 - 0.3$) if termination ripples create noise peaks in the peak table.
SMF map looks sparse/different from the raw Patterson map SMF evaluates an $N$-way minimum across all space group operators. It acts as a logical "AND" filter, suppressing background noise and leaving only true atomic positions in the absolute crystallographic frame. This is the expected behavior of SMF. Use the extracted SMF peaks as absolute candidate sites for seeding.
"WebGPU is not supported..." WebGPU disabled or unsupported browser. Use a Chromium v113+ browser with updated graphics drivers.
Atoms refusing to merge on special positions Minimum contact distance filter set too high. Lower the minimum contact distance slider down toward $0.00\text{ \AA}$ to permit overlapping symmetry-equivalent atoms.
"… generated positions, above this GPU's limit of …. Search for fewer atoms at a time." Evaluating one particle means holding every symmetry-generated copy of every asymmetric-unit atom in the compute shader's per-workgroup memory at once, and that memory is a small, fixed pool on any GPU. The number of copies is (atoms you've entered) × (space group operators), and this message means that product is larger than what your specific GPU's workgroup memory can hold – comp-Harker sizes this limit to your actual hardware rather than assuming a fixed number, but there is still a ceiling somewhere. See 6.3.B for the full explanation of why this is a hard refusal rather than a silent, partially-wrong result. Remove an atom from the asymmetric unit if you can (the message tells you exactly how many positions were requested and what your GPU's actual limit is, so you can see how far over you are), or accept that this particular combination of structure size and space group needs a GPU with more workgroup storage than the one you're currently running on.

Composition-Driven Runs

SymptomCause & Remedy
"Setting X carries no Wyckoff table" An old database is being served. Regenerate with cctbx_Harko_v1.py and check the console line naming which file won and how many settings carry projections. A stale legacy file in lib/ silently outranks a fresh sg/ unless the loader is told otherwise.
"Not loaded: observations.js" or a stale-module warning The modules must sit beside Harko.html and be served by the same server. If versions disagree, reload with Ctrl+Shift+R — save the modules and Harko.html together.
"Every particle scored exactly ..." The kernel dispatched but wrote nothing. Check the console for WGSL validation errors, and hard-reload: swarm_cc.wgsl is fetched without a cache-busting query.
Every candidate has an unbonded cation No distance windows were given — the header line says so explicitly. The search has nothing telling it that a polyanion should exist. Add S O 1.35 1.65 or equivalent.
Candidates all within 0.01 in CC Expected. The map cannot separate them; the chemistry can. Use the Structure quality panel (7.9). Also raise Iterations and Restarts: high-dimensional assignments need a longer run to be judged fairly.
"No combination of Wyckoff multiplicities sums to N" $Z$ is wrong, or the composition needs partial occupancy, which this search does not model. Try the neighbouring $Z$ values the suggestion line offers.
Wilson $B$ reported as unphysical Intensity rising with angle. Usually $Lp$ or a sharpening applied twice — check which intensity route the Log reports.