Tal — Crystallographic Application

1. Introduction

Tal is a sophisticated desktop application for crystallographic visualization and analysis. Built with .NET 8 and MonoGame, it offers a high-performance 3D rendering environment for exploring molecular structures, analyzing intermolecular interactions, and performing advanced computational tasks like cavity detection and Hirshfeld surface analysis.

The application uses a flexible docking interface, allowing you to customize the layout of various tool windows around the central 3D viewport.

Key Features:

• Loading and visualization of common crystallographic file formats (CIF, RES, XYZ, PDB, MOL, MOL2, GJF).
• Multiple display styles: CPK, Sticks, Balls and Sticks, Ellipsoids (ORTEP), and Polyhedrals.
• Hydrogen normalization to neutron-standard X–H distances (manual tool + automatic Hirshfeld path).
• Advanced coloring modes, including by element, symmetry, molecule, charge, and hydrophobicity.
• Comprehensive analysis tools for cavities, Hirshfeld surfaces, short contacts, and hydrogen bonds.
• Powerful overlay feature to superimpose multiple structures.
• Tools for creating publication-quality images (screenshots, POV-Ray export).
• Powder X-ray diffraction simulation with 2D fingerprint-style plotting.
• Crystal plane construction (hkl and mean plane via least-squares).

2. File Management

Operations related to opening, saving, and exporting data are found in the File menu (MainForm.Menu.cs / MainForm.Events.Menu.cs).

2.1. Open... (Ctrl+O)

Opens a structure file. Tal supports:

• CIF (*.cif) — multi-block via CifDivider / CifParser.ParseMultiBlock
• SHELX (*.res, *.ins)
• XYZ (*.xyz)
• PDB (*.pdb, *.ent)
• MOL/SDF (*.mol, *.sdf, *.sd)
• MOL2 (*.mol2, *.ml2)
• Gaussian (*.gjf, *.com, *.gau)

If a file contains multiple structures (e.g. multi-record CIF), Tal prompts to open all, only the first, or none. Loaded structures are added to StructureStorage / Structure List.

2.2. Open Recent

Recent files are managed by the singleton RecentFilesManager and persisted to %LocalAppData%\Tal\RecentFiles.json.

• Maximum 10 entries (MaxEntries).
• Deduplicated case-insensitively; re-opening moves a path to the top.
• Missing files are filtered out on load; unresolvable paths (e.g. broken UNC) are skipped.

2.3. Save (Ctrl+S) / Save As...

Both menu items call the same handler OnSaveAs: a save dialog is always shown. The currently selected structure from the Structure List is written via StructureWriter.Save(), dispatched by extension.

Supported write formats

CRYSTAL17 export and standard setting

Before writing .d12, Tal calls SpaceGroupTransformer.GetTransformToStandard(). If the structure is non-standard, a dialog is shown:

• YES — transform a copy to the standard setting (original in StructureStorage unchanged) and save.
• NO — save without transformation; a warning is written into the file comment.
• CANCEL — abort export.

2.4. Save Screenshot... (Ctrl+Shift+S)

Opens ScreenshotDialog: captures the active MonoGame viewport at a chosen supersampled resolution.

Final size: ResultW = round(baseWidth × Scale), clamped to MAX_DIMENSION = 8192 px per axis. Base size = current backbuffer of the active document.

2.5. Export to POV-Ray...

Exports the scene (atoms, bonds, unit cell, visible cavities, polyhedra, planes) to a *.pov file via PovRayGenerator. Uses current camera, lighting, background, and display flags. See section 13.

2.6. Close (Ctrl+W) / Exit (Alt+F4)

• Close — closes the active 3D document window.
• Exit — shuts down the application; layout and Flags settings are saved on exit.

2b. Edit, View, Window & Help Menus

Edit menu

For atom/molecule selection prefer the Selection menu and picking modes.

View menu

• Reset Layout — restores the default docking layout of tool windows (OnResetLayout).

Window menu

• 📋 Structure List — toggles the structure list dock window.
• 🎨 Display Options — toggles DisplayOptionsWindow (moved here from the Display menu).
• 🐞 Debug — Symmetry — opens the debug / symmetry log window (moved here from Help).

Help menu

• Contents (F1) — opens this help.html documentation.
• About — shows Flags.VERSION (currently 3.2 beta), author, and third-party credits (.NET 8, MonoGame, DockPanel Suite, SharpDX).

3. Selection & Manipulation

The Selection menu, picking toolbar, mouse modes, and the right-click context menu control interaction with atoms and molecules.

3.1. Picking Mode (toolbar)

The picking strip combo box (MainForm.Toolbars.cs) offers:

Also on the picking strip: Clear Measurements, and Show Labels for: All atoms / Not hydrogen atoms / Selected atoms / None.

3.2. Selection menu

• Select Atoms — all visible atoms in the current structure.
• Select Molecules — atoms of a molecule in a packing view.
• Select Unit Cell — unit-cell frame.
• Invert Selection / Clear Selection.

3.3. Context menu (right-click)

Built in MainForm.ContextMenu.cs, updated from the atom under the cursor:

• Hide selected / Hide reverse selected
• Show all
• Hide this atom / Hide this molecule / Hide this symmetry equivalent
• Hide short contacts / Remove hanging contacts

These call into AtomVisibilityManager (see section 16). After bulk hides, use Remove hanging contacts so contact lines are not drawn to invisible atoms.

3.4. Programmatic selection helpers

MeanPlaneCalculator (namespace Tal.Core.PowderDiffraction) exposes:

• CollectSelectedAtomIndices(Molecule) — indices with selected == true.
• CountSelectedAtoms(Molecule) — count of selected atoms.

Used by the Crystal Planes tool when building mean planes.

4. Visual Display

The Display menu, style toolbar, view toolbar, and docked Display Options window control appearance.

4.1. Style (toolbar + Display → Style)

Global style combo (exact labels from code):

• CPK — spheres scaled to van der Waals radii.
• Sticks — bonds only (no atom spheres).
• Balls and Sticks — spheres + cylindrical bonds.
• Ellipsoids — ORTEP-style thermal ellipsoids; isotropic atoms as spheres.
• Polyhedrals — coordination polyhedra for configured centres.

Per-atom overrides remain available via Atom Display Settings / colour tools.

4.2. Colour atom

Opens the atom colouring window. Modes include:

• Atomic Number (default element colours)
• Symmetry Equivalent
• Molecule Number
• Symmetry + Molecule
• Individual Color (Selected Atoms)
• Charge
• Hydrophobicity

4.3. Labels

Opens LabelsWindow (🏷️ Labels Settings): show/hide, font, shadow, default offset, optional property beside the name. Individual repositioning: picking mode Move Labels.

4.4. Show Hydrogens / Show Unit Cell / Reset View

• Show Hydrogens — OnShowHydrogens.
• Show Unit Cell — OnShowUnitCell.
• Reset View — OnResetView (camera to default).

4.5. View toolbar

• ▶ Animate — auto-rotate (see also AutoRotateWindow).
• View along a, b, c, a*, b*, c*.
• Rotate ±1° about X/Y/Z; rotate 90° (←90 →90 ↑90 ↓90).
• Pan arrows; Zoom+ / Zoom−.
• Persp. / Ortho. — toggle perspective vs orthographic.

4.6. Background... (Ctrl+B)

BackgroundColorWindow — solid colour or two-colour gradient; presets or custom colours.

4.7. Lighting...

LightingWindow — primary and secondary directional lights (direction, diffuse, ambient, specular colour/power/intensity); optional per-pixel lighting.

4.8. Atom Display Settings

Window title 🎨 Atom Display Settings (AtomDisplayWindow):

• BALL SIZES: Non-Metal, Metal, Hydrogen (%) of ball size.
• STICK SIZES: Radius, Length Factor.
• UNIT CELL LINES: Thickness, Length Factor.
• ELLIPSOID PARAMETERS: Probability (%), Disk Thickness, Base Radius; Show Hydrogen as Sphere; Show Disks.

Buttons: Default, Reset, Apply. Probability is the source of truth for ellipsoid scale (χ² table).

4.9. Atom Colors (global)...

Edits the persistent per-element palette (atomic number index) via AtomPropertiesWindow / related helpers. Overrides are stored in Flags.AtomColorOverrides and applied on startup. Reset Selected / Reset All restore defaults.

4.10. Short Contacts...

Opens ShortContactsWindow (“Define Short Contacts”):

• Fixed distance (Å) or Sum of vdW radii (with ± Å offset controls).
• Intermolecular / Intramolecular (separated by > 3 bonds).
• Select whole molecule (Expand Contacts).

Detected contacts are drawn as dashed lines. Calculate → Contacts... reuses this path (OnCalcContacts → OnShortContacts).

4.11. Scene Quality...

SceneQualityWindow (labels from Designer; defaults from Flags):

• SSAA: Off (1×), 2×, 4× — default Flags.SSAASamples = 2.
• MSAA on scene RT: Off, 2×, 4×, 8× — default MSAAOnSceneRT = 4.
• Anisotropy of textures — default AnisotropyLevel = 4.
• Progressive AA when the camera is static: Flags.ProgressiveAASamples (default 4).

4.12. Display Options (Window menu)

DisplayOptionsWindow quick toggles:

• Display group: Packing, Short Contact, Asymmetric Unit, H-Bond, Auto centre.
• Options: Show hydrogens, Show cell axes, Label atoms, Depth cue, Z-Clipping, Stereo.
• Buttons: Reset, Contacts..., More Info, Powder pattern.

Depth cue details: DepthCueWindow. Stereo and Z-clipping are experimental depending on build.

4b. Symmetry Elements

The Symmetry Elements window visualizes the symmetry operations of the crystal structure in the style of Mercury and the ITA (International Tables for Crystallography) conventions. It shows rotation axes, screw axes, mirror planes, glide planes, inversion centers, and rotoinversion axes directly in the 3D viewport.

4b.1. Opening the window

The window is opened from Display → 🔷 Symmetry Elements... or with the keyboard shortcut Ctrl+Shift+E (OnSymmetryElements in MainForm.Events.Menu.cs). It is a dockable DockContent (SymmetryElementsWindow) and can be placed anywhere in the layout, but is shown by default on the left dock with the standard tool-window width.

4b.2. Master switch

The Show Symmetry Elements checkbox at the top of the window controls Flags.ShowSymmetryElements. When it is off, the renderer skips the entire pass and nothing is drawn. This is the only global toggle; every other control in the window affects one specific kind of element.

4b.3. Element types & colors

The window groups every symmetry element into a two-column grid: the left column is a checkbox that controls visibility, the right column is a small color swatch that opens a standard ColorDialog. Each row is independent: you can hide, for example, only the 3₂ screw axes while keeping 3₁ visible, and each one has its own color.

4b.4. Rendering style (ITA markers)

Elements are drawn in the style of the International Tables for Crystallography:

• Rotation axis 2 — a rod with a conical arrow head at the end (representing the 180° flip direction).
• Rotation axes 3, 4, 6 — a rod with the classic ITA marker at the end: a filled triangle (3), a filled rhombus (4), or a filled hexagon with an inner ring (6).
• Screw axes — the same rod, but dashed. The dash length encodes the subscript: 2₁, 3₁, 4₁, 6₁ have dense dashes; 3₂, 4₂, 6₂ have wider dashes; 6₃ has the widest. This matches the "gęsty vs rzadszy" convention from the on-line reference style.
• Rotoinversion axes (−n) — a rod with the ITA polygon marker for the order, plus a small ring drawn exactly at the inversion center (which is not always in the middle of the visible segment — a −4 axis can have its inversion point anywhere along its line, e.g. at z = ¼).
• Mirror plane (m) — the polygon of intersection between the plane and the unit cell, drawn as a translucent filled surface with rod edges.
• Glide plane — same as mirror, plus a small arrow showing the glide vector (the component of the symmetry operation that lies inside the plane). The letter (a/b/c/n/d/g) is derived from that vector, not from the translation alone.
• Inversion center (−1) — a small sphere plus a three-axis "crosshair" (three short rods along X, Y, Z), so it stays visible even when the sphere would be hidden inside a molecule.

4b.5. Sizes and transparency

The bottom group of the window controls the sizes and opacities. These are the same parameters for every element of that kind — the per-element controls are only for visibility and color.

All five sliders are clamped to their declared range, so an old tal_settings.json with values outside the range (for example a leftover 0.185 for axis thickness) cannot crash the window — it is silently clamped to the nearest legal value.

4b.6. Live updates and colors

Colors and sizes are read live from Flags.Instance every frame by the renderer CellSymmetryRenderer. The pipeline is:

1. A control changes a Flags field, calls Flags.MarkDirty(), and calls the shared RefreshView() helper.
2. RefreshView() bumps Flags.VisualRevision, calls MainForm.InvalidateSymmetryRenderers() (which sets a rebuild flag in every renderer), and calls MainForm.RefreshAllDocuments().
3. On the next frame, CellSymmetryRenderer.Draw recomputes the color for each element via GetColorForElement(kind, order, subscript), which reads the current Flags values. The rod thickness, marker size, inversion radius and plane alphas are also read from Flags inside the drawing code.
4. The rebuild flag triggers a full rebuild of the cached _elements list from the current Molecule if the geometry needs to change (for example after a different symmetry operation became visible). For simple color changes, the rebuild is skipped and only the colors change.

This makes the behavior identical to what you see in Atom Color and Atom Display Settings: no Apply button, no OK/Close required, no need to reopen the window . The change appears as soon as the next frame is rendered.

4b.7. Refresh and Reset buttons

The window has two buttons at the bottom:

• 🔄 Refresh — forces an immediate redraw without changing any flags. Useful after loading a new structure, after switching documents, or after a Transform space group operation that changed the underlying Molecule.SymmetryOperations.
• ↺ Reset to Defaults — restores every color and every size to the built-in defaults, re-enables every element type, and refreshes the view. It does not touch tal_settings.json; changes are persisted on exit along with every other flag.

4b.8. What is shown and where it comes from

Every element copy that intersects the unit cell [0,1]³ is generated. The pipeline is:

1. The full symmetry group is closed from the generators stored in Molecule.SymmetryOperations. If the CIF (or RES) file only lists generators, the closure expands them up to the full group (max 192 operations).
2. Every operation is classified by trace and determinant of its linear part: det = +1 with trace −1, 0, 1, 2 is a rotation of order 2/3/4/6; det = −1 with trace 1 is a mirror plane; det = −1 with trace −3 is an inversion center; and the remaining det = −1 cases are rotoinversion axes −3, −4, −6.
3. For each operation and for every lattice translation L ∈ {−2 … 2}³, the intersection of the element with the unit cell is computed in fractional coordinates, then converted to Cartesian using the crystal's fractional-to-Cartesian matrix.
4. Duplicate copies (same axis line, same plane, same point) are deduplicated with a geometric key, so a 4-fold axis in a cell with four equivalent copies is not drawn four times on top of itself.
5. Screw subscripts are derived from the intrinsic translation along the axis: w = (1/n)·Σ Rᵏ·t, projected onto the axis direction and quantized to 1/n. The rotation sense (right-hand rule in Cartesian space) decides whether a 3-fold becomes 3₁ or 3₂.
6. Glide letters are derived from the glide vector w = (t + R·t) / 2. A single half-translation along a/b/c gives a/b/c; two half-translations give n; quarter-translations give d; anything else is g.

4b.9. Interaction with other tools

• The window respects Flags.ShowSymmetryElements as its master switch: unchecking it disables everything in the window but does not close the window.
• Symmetry elements are not exported to POV-Ray: the POV-Ray scene contains atoms, bonds, the unit cell, planes, polyhedra, and cavities, but not the auxiliary ITA markers.
• Symmetry elements are affected by the unit-cell repetition settings. If you display a supercell (for example 0 → 2 along a), the symmetry elements are still computed only for the base cell — this is intentional, because drawing every copy in a supercell would clutter the view without adding information.
• The colors and sizes are part of tal_settings.json and are saved automatically when the application exits. There is no per-project persistence.

4b.10. Keyboard shortcuts and menu map

5. Calculations

The Calculate menu (MainForm.Menu.cs). Several entries are fully wired; a few are still stubs.

5.1. Centroids...

Not implemented. OnCalcCentroids only sets the status bar. Reserved for the geometric centroid of selected atoms.

5.2. Planes...

Opens PlanesWindow (dock left). Define planes from selected atoms (mean plane via MeanPlaneCalculator) or Miller indices (hkl). Visibility, colour, transparency; distance-to-plane measurements. See section 9.

5.3. Packing...

Opens PackingWindow (“📦 Packing”). Controls fractional extent of the displayed block:

• Per-axis Min / Max for X, Y, Z (numeric up-downs).
• Shortcuts: +0.5, 0→2, 0→3 per axis; All +0.5, All 0→2, All 0→3; Reset.
• Packing checkbox sets Flags.DrawPacking and syncs Display Options when present.

Typical presets in code: asymmetric-style 0–0.5, 0–1 unit cell, 0–2 / 0–3 supercells.

5.4. Contacts...

Alias of Display → Short Contacts (OnCalcContacts → OnShortContacts). See §4.10 and Expand Contacts picking mode.

5.5. Powder Pattern...

Opens PowderDiffractionWindow when an active document with atoms exists. See section 8.

5.6. Molecular Shell...

Not implemented. OnCalcMolecularShell only updates the status bar.

5.7. Graph Sets...

Not implemented. OnCalcGraphSets only updates the status bar (reserved for H-bond graph-set analysis).

5.8. Cavity Settings... (Ctrl+Shift+C)

Opens CavitySettingsWindow (title 🕳️ Cavity Settings). Full algorithm and parameters: section 6. UI highlights: Enable Cavity Display; Resolution (Å); VdW Multiplier; Smoothing; Transparency; Use thermal motion (Uij/Uiso); Detect / Clear / Export Report; display modes Hidden / Normal Cavity / Property Map / Topology Map / Channel Paths / Sphere Fits; property types Charge, Hydrophobicity, Atomic Number, VdW Radius, Occupancy.

5.9. Transform space group...

Transforms the active structure between equivalent settings (e.g. P2₁/c ↔ P2₁/n) via SpaceGroupTransformer. Updates cell, coordinates, ADP tensors, and symmetry in place. See section 13.

5.10. Normalize Hydrogens...

Opens HydrogenNormalizationWindow (🧪 Hydrogen Normalization). Loads the active molecule, then docks the tool window. Uses HydrogenNormalizer (neutron-standard X–H lengths, optional skip of bridging H, Apply / Reset / cache of originals). Details and Hirshfeld auto-path: section 7.2 and the dedicated section below if present.

6. Cavity Analysis

Tal includes a powerful and robust engine for the detection and analysis of voids and cavities in crystal structures. This is a key feature for studying porosity, host–guest chemistry, and gas storage potential.

6.1. Cavity Detection & Display Window

This is the main control panel for all cavity-related operations, opened from Calculate → Cavity Settings... (Ctrl+Shift+C).

Detection Parameters

Controls in CavitySettingsWindow (Designer labels) and corresponding Flags / detector fields:

• Enable Cavity Display (chkEnableCavity) — master switch for cavity visualization.
• Resolution (Å) (numResolution, Flags.Resolution, default often 0.10–0.15 Å) — voxel size of the SDF grid. Smaller is more accurate and much more expensive (voxel count ~ 1/res³).
• VdW Multiplier (numVdwMultiplier, Flags.VdWValue, default 1.00) — scales every atom’s van der Waals radius. See section 15.
• Use thermal motion (Uij/Uiso) (chkUseThermalMotion, Flags.UseThermalMotion) — when enabled, the effective atomic shape is an ellipsoid from ADP data instead of a pure sphere. See section 6.2. The thermal term is the 1σ (RMS) amplitude of the ADP; it is independent of the ellipsoid probability used for the ORTEP display (see ADP conventions).
• Transparency (trackTransparency, Flags.transparency) — cavity mesh opacity.
• Smoothing (numSmoothing, Flags.smuth, default 3) — Gaussian / Laplacian smoothing iterations before Marching Cubes.

Actions

• Detect Cavities: Runs the main detection algorithm on the active structure. The process uses a Signed Distance Function (SDF) and a flood-fill algorithm to identify regions of empty space that are not connected to the outside of the crystal.
• Clear Cavities: Removes all detected cavities from the current document, freeing memory and hiding the visualization.
• Export Report: Saves a text file containing a complete single-run report — crystal data, grid info, per-cavity table, channels, sphere fits, gas/solvent compatibility, single-run QSPR descriptors, and refinement metadata. See section 6.7.

Display Modes

Once cavities have been detected, you can choose how to visualize them:

• Hidden: Cavities are not shown.
• Normal Cavity: Renders the cavities as semi-transparent isosurfaces (meshes).
• Property Map: Colors the cavity surface based on a molecular property of the nearby atoms. The UI offers: Charge, Hydrophobicity, Atomic Number, VdW Radius, Occupancy.
• Topology Map: Colors the cavity surface based on topological features (number of wall contacts).
• Channel Paths: Analyzes the cavity network to find and visualize interconnected 1D channels that percolate through the crystal.
• Sphere Fits: Finds the largest spheres that can fit inside each cavity and renders them, giving a quick visual estimate of the cavity's internal dimensions.

6.2. The Cavity Detection Algorithm

Tal's cavity detector is a sophisticated implementation designed for speed and accuracy, even for large and complex structures.

1. Atom generation: The scanner first generates all atoms of the unit cell together with their neighbours (the atoms in the 26 periodic images that are within ScanMargin = 2 Å of the cell boundary). Every atom is stored with its effective radius:

   r_eff(i) = r_vdw(i) · VdWValue + ProbeRadius

   This radius is baked into the atom and into the SDF that is computed from it. There is no separate "base SDF" cached across VdW values — see section 15.
2. Signed Distance Function (SDF) calculation: A high-resolution 3D grid (the "voxel grid") is created over the unit cell, extended by a margin of ScanMargin on each side. For every voxel the algorithm computes the signed distance to the surface of any atom:

   d_i       = distance from the voxel to the surface of atom i
               = |p − c_i| − r_eff(i)      (spherical case)
   SDF(p)    = min_i d_i
   

   SDF is negative inside an atom and positive outside. It is truncated at MaxDistFine = 5 Å — voxels further than 5 Å from any atom are set to that value, which is more than enough for the flood-fill to distinguish "empty" from "occupied".
3. Handling of anisotropic displacement parameters (ADPs): When Use thermal motion (Uij/Uiso) is enabled, the distance is not computed against a sphere but against an ellipsoid whose principal half-axes are:

   semi-axis_i = r_vdw(i) · VdWValue + ProbeRadius + √|λ_i|

   where √|λ_i| is the i-th row length of the Un matrix produced by the Eigen decomposition of the ADP tensor (see section 6.2). The half-axes are stored per atom and used inside the SDF hot loop.
   Symmetry-equivalent atoms. The ADP of every symmetry-generated atom is rotated together with its position:

   Un' = Un · R_cartᵀ        Ue' = Ue · R_cartᵀ
   R_cart = A · R_frac · A⁻¹

   where A is the fractional → Cartesian matrix (x_cart = A · x_frac, column convention) and R_frac is the 3 × 3 linear part of the symmetry operation (translations do not rotate the ADP). This is valid for every crystal system, including hexagonal, trigonal and rhombohedral settings; no diagonal-only or sign-only shortcut is used on this path.
   The split between the two consumers is deliberate:
   • Cavity detector / VoxelScanner reads the rows of Un' as the ellipsoid axes. Because Un' = Un · R_cartᵀ, each row i has length √|λ_i| (the 1σ half-axis is preserved) and direction R_cart · e_i (the eigenvector is rotated). Transformation is ignored by this path (ReExpandSymmetry = false).
   • ORTEP / Draw multiplies on the left by Transformation = R_cart, so the rendered quadric is R_cart · Un · R_cartᵀ · prob — the same form as the old S · Un · S for diagonal sign-only operations, generalized to full R_cart.
   The fallback path (TransformUnSignOnly + Transformation = diag(sign1)) is used only when the cell matrix A is not invertible; it is bit-compatible with the historical CalculateSymmetryEquivalents. See the header of Molecule.Symmetry.cs for the full derivation.
4. Periodicity and flood-fill: The grid is extended into a margin around the unit cell. A flood-fill algorithm starts from the outermost boundary of this extended grid (the "exterior") and fills all connected empty voxels. Any empty voxels not reached by this fill are considered to be part of a "buried" cavity.
5. Wall detection: For every detected cavity the algorithm records which unit-cell walls it touches (WallXMin/Max, WallYMin/Max, WallZMin/Max). This is what tells you later whether a cavity is closed or part of a percolating channel.
6. Marching Cubes mesh generation: The Marching Cubes algorithm is applied to the SDF grid to generate a smooth triangle mesh for the surface of each cavity. The quality of the mesh is controlled via the Smoothing parameter (see section 6.8).
7. Post-processing per cavity: For each cavity the engine also computes:
   • volume (voxel count × voxel volume) and its error estimate,
   • surface area and its error estimate,
   • sphericity,
   • center of mass,
   • fractional coordinates of every voxel (used by property mapping),
   • per-voxel charge, hydrophobicity and Z-effective (used by the property map display mode).

q(p) = ex²·invRx² + ey²·invRy² + ez²·invRz²

6.3. Resolution, memory and performance

The number of voxels in the supercell grid grows cubically with the resolution:

N_voxels ≈ (cell_size / resolution)³ · (1 + 2 · margin_fraction)³

Tal enforces a hard limit Flags.MAX_VOXELS on the total voxel count. If the grid exceeds this limit, detection is aborted with an error message and a suggested coarser resolution. Two additional thresholds trigger warnings:

• 15 million voxels: a "large grid" warning — detection will take a while but is allowed.
• Grid dimension > 32767 per axis: the underlying VoxelGrid stores resolutions as short, so anything larger cannot be represented. Detection is aborted.

The memory footprint is roughly 4 bytes × N_voxels for the density grid, plus about the same again for the periodic copy, gradient arrays, and the voxel mask. A 512³ grid (≈ 134 M voxels) already needs ~1.5 GB of RAM, which is why the default resolution is capped around 0.10–0.15 Å for typical organic crystals.

6.4. Why the void volume shrinks: an erosion model

Both the VdW multiplier and the anisotropic (ADP) mode make every atom larger and therefore reduce the free space. The reduction is often much larger than the relative growth of the atoms , and the number of cavities usually changes as well. This section gives a simple model that predicts the size of the effect from two numbers that Tal already reports (void volume V and void surface S), and explains when a void splits or disappears. Throughout, κADP = V(ADP = YES) / V(ADP = NO); this ratio is not the same quantity as the κ = V(1.5)/V(1.0) of the VdW sweep (section 15.6).

The core idea: inflating the atoms erodes the void

Let d(p) be the distance from a point p to the nearest atom surface (this is the SDF of section 6.2). A point belongs to the void when d(p) > 0. If every atom is inflated by the same amount δ, the point stays in the void only when d(p) > δ. In other words, the new void is the old void eroded by a ball of radius δ: every wall of the void moves inwards by δ. In ADP mode δ is not the same for every wall: it is the displacement amplitude of the nearest atom in the direction of the wall normal (see "What is δ for an ellipsoid?" below).

Effect 1 — loss of volume

Because every wall moves inwards by δ, the volume lost per unit of δ is exactly the area of the wall. For the eroded void this gives the exact relation

dV/dδ = − S(δ)          V(δ) = V₀ − ∫₀^δ S(t) dt

where S(δ) is the void surface after the erosion. To first order in δ:

κ_ADP ≈ 1 − δ·S₀/V₀ = 1 − δ/ℓ          ℓ = V₀/S₀

The characteristic length ℓ = V/S of the void measured without ADP decides everything: a void is robust when δ ≪ ℓ and collapses when δ is comparable to ℓ. For simple shapes with inscribed radius ρ the erosion can be written exactly:

All three are of the form (1 − δ/ρ)m with m = 3, 2, 1 for pockets, channels and layers. Real voids are mixtures of these. A void of large volume but small ℓ (thin, branched, low sphericity) loses most of its volume at a very small δ, even if its equivalent-sphere radius is large. Note that a void of inscribed radius ρ (the "max sphere radius" of the sphere fit) disappears completely at δ = ρ, regardless of its volume.

Effect 2 — change of connectivity

A neck of free width w (the diameter of the sphere that can pass through it) closes when its two walls have both moved by δ, i.e. when δ ≥ w/2. A void network therefore splits at the bottleneck radius δc = wmin/2 of the narrowest neck on its connecting path (MinWidth of the channel analysis, section 15.10, is a local diameter, so δc = MinWidth/2). Splitting changes the number of cavities but, taken alone, hardly changes the total volume (a neck contributes very little volume). The total volume falls additionally through fragmentation only in two ways:

• fragments whose volume is below Min cavity volume are discarded, and
• each fragment is eroded separately, and small fragments have a small ℓ and vanish quickly (δ ≥ ρ).

Consequently, the number of cavities, the total volume and κADP depend on Min cavity volume. Always quote that parameter together with κADP. The effects of volume loss and of connectivity change are two observables of the same erosion; they are not two independent factors that multiply each other.

What is δ for an ellipsoid?

An atom with the displacement tensor U (Cartesian, Ų) is represented by the 1σ ellipsoid with half-axes √λi. A wall with outward normal n is moved by the extent of this ellipsoid in the direction n, that is, by

δ(n) = √(nᵀ·U·n)         √λ_min ≤ δ(n) ≤ √λ_max

So it is the vibration towards the void that matters, not the mean amplitude √Ueq. Typical values of √Ueq for C, N, O atoms are ≈ 0.12–0.16 Å near 100 K (Ueq ≈ 0.015–0.025 Ų) and ≈ 0.20–0.25 Å at room temperature (Ueq ≈ 0.04–0.06 Ų). Here U is the 1σ tensor of section 6.2; a thermal term scaled by a factor s (e.g. 1.538 for ellipsoids drawn at 50%) would scale δ by the same factor.

Illustration: a lattice without necks (FCC)

The following didactic model contains only the volume effect and no necks. It is not a molecular crystal. Take an FCC lattice with a = 3.615 Å (4 atoms per cell, Vcell = a³ = 47.24 ų, nearest-neighbour distance dnn = a/√2 = 2.556 Å), a base radius rbase = 1.00 Å (deliberately smaller than a real VdW radius, so that free space exists) and an isotropic increase δ. As long as the spheres do not overlap (δ < δtouch = dnn/2 − rbase = 0.278 Å):

V_free(δ) = a³ − (16π/3)·(r_base + δ)³

At contact the spheres fill 74.0% of the cell (the close-packed value), so the free volume is 26.0%, which is a useful check of the table. The model illustrates how fast the free volume falls when the radius grows by only δ/rbase = 10–25%. The numbers are specific to this lattice and are not transferable to a real crystal (real VdW radii are larger, so the same δ is a smaller relative change; real voids have a different ℓ). Also note that this lattice has no necks to close: for a point probe the free space of FCC stays connected through the triangular windows up to r = dnn/√3 = 1.476 Å, i.e. well beyond the contact. Fragmentation is therefore described by Effect 2 above, not by this model.

6.5. Meaning of κ_ADP and relation to packing coefficients

The packing coefficient of Kitaigorodskii, k = Vmolecules / Vcell (for hard spheres at most 0.74), is a static measure: it counts the cell volume occupied by the VdW envelopes of the molecules. Since 1 − k is the whole space not covered by the VdW envelopes, including the narrow interstices that Tal does not report as cavities, the void fill fraction of Tal (probe 0, VdW multiplier 1.0) is at most 1 − k. The two numbers are not expected to be equal.

κADP adds a different piece of information: how much of the reported void lies within a distance of about δ of the walls. By section 6.4, this depends on ℓ = V/S of the void, hence on its shape (low sphericity means a small ℓ). A high κADP (close to 1) means that the voids are compact and survive the thermal inflation; a low value means that the reported volume consists largely of thin regions near the walls.

6.6. Worked example: ADP versus spheres

Consider a molecular crystal (P2₁/n, cell 11.80 × 10.39 × 19.46 Å, β = 103.73°, V_cell = 2317.08 ų, 224 atoms in the cell, data collected at 293 K), detected twice with identical parameters (resolution 0.15 Å, VdW multiplier 1.0, probe radius 0.0, Min cavity volume = 0.1 ų, the default). The only difference is the Use Anisotropic VdW flag.

Reading the result with the erosion model

• The void at ADP = NO is far from compact. Its volume corresponds to a sphere of radius 4.3 Å, but its sphericity is only 0.164 and ℓ = 0.23 Å (a compact sphere of that size would have ℓ = R/3 ≈ 1.42 Å). The estimate for a single sphere, (1 − δ/R)³ with δ = 0.15 Å, would predict κ_ADP ≈ 0.90. The observed value is 0.21. The ℓ-based estimate is the appropriate one here.
• First-order estimate. With ℓ = 0.232 Å we obtain κADP ≈ 1 − δ/ℓ = 0.57, 0.35 and 0.14 for δ = 0.10, 0.15 and 0.20 Å, respectively. Conversely, the observed κ_ADP = 0.21 corresponds to an effective δeff = (1 − κADP)·ℓ ≈ 0.18 Å, which is in the range of the RMS amplitudes of C, N, O atoms expected for the data collected near room temperature. Because δeff is comparable with ℓ, the first-order formula is only an order-of-magnitude check here (the higher-order, curvature terms and the vanishing of small pieces are not included).
• Sphericity increases (0.164 → 0.399) mostly because of selection: the thin, high-S/V parts of the void disappear first, and the survivors are compact pockets. It does not mean that the original void became rounder.
• One void became six (equivalent radii 1.0–1.8 Å, volumes 4.8–24.5 ų): the necks of the original void are narrower than about 2δ.

Quick estimate before running ADP mode

1. Run the detection with ADP = NO and read V and S from the report.
2. Compute ℓ = V/S.
3. Take δ ≈ √Ueq of the atoms lining the void (see above).
4. If δ ≪ ℓ, the void is robust. If δ is comparable to or larger than ℓ, expect a strong decrease of volume and a split into several cavities. Cavities whose max sphere radius is smaller than δ will disappear.

Which run is the right one?

The two runs answer different questions.

• ADP = NO is the geometric void of the model with the atoms at their mean positions. It is the quantity to compare with published geometric void analyses. It is not a 0 K model: the coordinates and the cell are those of the experiment, i.e. of the temperature of the measurement.
• ADP = YES is the space that remains outside the 1σ displacement envelope of all atoms, that is, the space that atoms practically never visit. It is a conservative (lower) estimate of the space available to a guest, not the space seen at one instant: at any moment atoms occupy only one point of their envelope, and a guest may pass through necks that are closed in the envelope (gating). For adsorption and diffusion at room temperature the ADP run is therefore the more realistic of the two, but still a model, and not a prediction of accessibility.

6.7. Single-run report

After a successful detection, the Export Report button writes a complete text report. It is deliberately self-contained: everything needed to reproduce or interpret the numbers is in the file. Sections:

• Header: timestamp, structure name, space group, and the full set of detection parameters (resolution, VdW multiplier, probe radius, thermal-motion flag, min cavity volume, scan margin, base SDF radius, smoothing level).
• Crystal data: unit cell (a, b, c, α, β, γ, volume), crystal system, symmetry count, atom count in asymmetric unit and in the cell, molecule count, density, formula, formula weight, F(000), and the charge/hydrophobicity ranges used for property mapping.
• Grid / SDF data: periodic cell resolution (ResX × ResY × ResZ), total voxels, margins, supercell resolution (with margin), total supercell voxels, estimated memory, voxel size, voxel volume, max SDF distance, and non-orthogonal flag.
• Cavity summary: total count, primary/duplicate/periodic/closed splits, count with channels, with mesh, with sphere fits; total void volume and error; total void surface and error; volume-to-surface ratio; cell volume; void fill fraction; min/avg/max void volume; average sphericity; largest sphere radius and diameter.
• Per-cavity table: ID, volume, surface, sphericity, V/S ratio, wall count, voxel count, max sphere radius, effective diameter, channel directions.
• Crystal channels: per-cavity per-direction table (length, avg/min/max width, effective diameter, helicity, cell crossings), direction summary, network type (1D / 2D / 3D / isolated).
• Sphere fits: per-cavity count, max radius, effective diameter, packing efficiency, usable volume, number of spheres on the cell wall.
• Single-run QSPR descriptors: I_frag, porosity, V/S ratio, mean charge variance, mean hydrophobicity variance, mean Z-effective.
• Refinement / experiment / publication metadata: R1 (observed / all), wR2 (all / observed), GooF, reflections, parameters, restraints, shift/su, Δρ max/min/RMS, weighting scheme; reflection statistics; temperature, wavelength, radiation type, crystal colour/shape/size, absorption coefficient.

6.8. Smoothing and the Marching Cubes pipeline

Smoothing is applied in two stages, both controlled by the Smoothing slider:

1. Grid-level smoothing (applied before Marching Cubes):
   • Gaussian blur — a separable 3-tap convolution with a Gaussian kernel whose σ is derived from the smoothing level. This softens the SDF field itself.
   • Laplacian smoothing — a few iterations of the 6-neighbour average, weighted by a smoothing factor (0.25–0.40 depending on the level). This removes the residual staircase from the Gaussian step.
   • Gradient scaling — a small multiplier (0.80–1.00) applied to the deviation from the mean density. Lower values produce flatter, smoother surfaces at the cost of some volume shrinkage.
2. Mesh-level snapping (applied after Marching Cubes):
   • Vertex welding — vertices that are within SnapTolerance of each other are merged. The tolerance is 0.0004–0.0008 Å depending on the level.
   • Normal recomputation — vertex normals are recomputed from the SDF gradient (not from the triangle geometry), which gives smoother shading on the final mesh.

Smoothing level 0 disables both stages. Levels 1–5 follow the presets shown in section 4. Higher levels produce smoother surfaces but also systematically shrink the cavity volume by 5–15%, so for quantitative work it is recommended to use a low smoothing level (0 or 1) and only enable smoothing for publication figures.

6.9. Clearing and resetting

Clear Cavities removes all detected cavities and their meshes, but it does not clear the cached atoms or the voxel grid — a subsequent detection at the same parameters reuses the cached data and is nearly instant. To fully reset the detector (free all memory, drop the atom cache, drop the spatial index), use Calculate → Cavity Detection → Force Reset. This is also the recommended action when switching to a different crystal structure.

6.10. Space-group transformation — consistency test

The same crystal structure can be expressed in several equivalent settings of the same space group — for example, P2₁/n and P2₁/c are the two standard monoclinic settings of space group #14. They describe the same crystal, with the same unit-cell volume and the same number of atoms per cell, but with different cell parameters and different symmetry-operation strings. Because Tal can transform between settings in place (Calculate → Transform space group...), a natural question is: does the cavity detector give the same answer for both settings? The short answer is yes — but only if the ADP tensors of the symmetry equivalents are rotated correctly. This section describes the consistency test and what to expect from it.

Why this test matters

When the atoms carry anisotropic displacement parameters (ADPs), the cavity detector replaces each atom with its 1σ thermal ellipsoid. The shape and orientation of that ellipsoid are stored as the 3×3 tensor Un, and the detector reads the rows of Un as the three half-axes of the ellipsoid (see section 6.2, "Anisotropic VdW — exact math"). If the transformation between two space-group settings rotates the tensor in the wrong convention — for example, full rotation R·U·RT instead of the row-wise U·RT that the detector expects — the row lengths no longer match √λi, the ellipsoids come out physically wrong, and the cavity volumes differ between settings. The test below detects exactly this class of error.

Test procedure

1. Load a structure and record the cell parameters (a, b, c, α, β, γ, V, Z, atoms in asymmetric unit).
2. Run a normal cavity detection (Calculate → Cavity Settings → Detect) with Thermal motion = YES, at a chosen resolution (for example 0.15 Å). Export the single-run report.
3. Transform the structure to the other equivalent setting (Calculate → Transform space group...) and confirm that:
   • the unit-cell volume V is unchanged,
   • the number of atoms in the asymmetric unit is unchanged,
   • the number of symmetry operations is unchanged.
4. Run cavity detection again, at the same resolution and the same settings as in step 2. Export the second report.
5. Compare the two reports on the quantities listed in the table below. For a fully consistent implementation, ΔV should be well below 1 % at resolution 0.15 Å, and it should decrease when the resolution is reduced to 0.06 Å.
6. Optionally, repeat the whole cycle with Thermal motion = NO (spheres). This isolates the contribution of the ADP rotation: any remaining difference in V between settings must be due to the rasterization of the two unit cells, not to the ADP handling.

What to compare

Worked example — structure 488_auto

The table below shows the full set of measurements that were taken on the same crystal of 488_auto (space group #14, Z = 4, 38 atoms in the asymmetric unit, 152 atoms per cell) in the two settings P2₁/n and P2₁/c. The test was performed twice at each of two resolutions (0.15 Å and 0.06 Å), once with Thermal motion = YES and once with Thermal motion = NO.

Two things are worth pointing out in this example:

• The failing case was not a numerical artefact. Going from 0.15 Å to 0.06 Å changed the difference from +9.1 % to +9.6 % — it grew slightly instead of shrinking. A pure rasterization difference would have decreased by a factor of about 4 (in V) when the voxel volume shrinks by 16×. The stability of the offset was the smoking gun for a convention bug.
• The passing case (0.06 Å, ADP = NO) confirmed the diagnosis. With spheres only, both settings gave essentially the same void volume (0.06 % apart). That proved the atomic positions were correct in both settings, and that the difference had to come from the ADP path.

Known residual differences after the fix

Even after the ADP rotation is correct, two differences between the two settings remain and are expected:

• Surface area. The total void surface differs by roughly 5–10 % between two non-orthogonal cells. This is intrinsic to Marching Cubes: the same iso-surface sampled on two grids whose axes are not aligned with the cell walls produces slightly different triangle counts. The difference does not disappear with resolution; it only stabilises.
• Number of voxels inside the cavity. Because the two grids have different voxel volumes (the periodic cell grid has different ResX × ResY × ResZ for the two cells, and the volumes are not equal), the same cavity is counted with different numbers of voxels. This is a bookkeeping difference; it does not affect the total volume.

7. Hirshfeld Surface Analysis

This advanced tool, opened from the Advanced menu, is essential for understanding and quantifying intermolecular interactions in molecular crystals. It is based on the Hirshfeld partitioning scheme, which defines a molecule's space in a crystal based on the electron density.

7.1. The Workflow

1. Open the Hirshfeld Surface Analysis window from the Advanced menu.
2. The window automatically builds the fresh unit cell and enumerates all molecules found in it. The Target molecule dropdown is filled with the detected molecules (each entry shows the number of atoms).
3. Select a target molecule from the dropdown, or click Use current selection to pick the molecule that contains the currently selected atoms in the main viewport.
4. Choose the density model (isotropic or anisotropic) and the surface coloring mode.
5. Click the Generate button. The application computes the isosurface, the per-vertex properties (d_e, d_i, d_norm, charge, hydrophobicity) and the contact-type decomposition.
6. The resulting surface is displayed in the 3D view, colored according to the selected property. The corresponding 2D fingerprint plot can be opened with Show Fingerprint....

7.2. Hydrogen Normalization (automatic)

Hirshfeld surface analysis, by convention (CrystalExplorer, Tonto and all crystallographic tools), requires neutron-normalized hydrogen positions. X-ray diffraction locates the maximum of the electron density of an H atom, which is systematically shifted toward the heavy atom relative to the true nuclear position. Using X-ray H positions in a Hirshfeld analysis introduces a bias of typically 0.10–0.20 Å on every X–H bond, which distorts both the surface and the fingerprint decomposition. Tal therefore always normalizes H before generating a Hirshfeld surface — silently, without a dialog.

What "normalized" means

Each hydrogen atom is placed at its neutron bond length along the direction defined by the heavy-atom coordinate and the current H coordinate. The bond lengths used are tabulated in HydrogenNormalizer.NeutronLengths:

Detection and application

When Generate is pressed, the pipeline calls Molecule.CheckHydrogenState, which uses the pure function HydrogenNormalizer.Check to analyze every H in the asymmetric part without modifying anything. Each H is classified into one of five states:

• NoStructure — the asymmetric part is empty; the pipeline proceeds normally (the case never occurs for a valid molecule).
• NoHydrogens — the structure contains no hydrogen atoms. No work needed.
• AlreadyNormalized — every H is within 0.005 Å of its neutron target. The pipeline logs a success message and skips the normalization step.
• DisabledByUser — the user clicked Reset in the Hydrogen Normalization window at some point in the session. The pipeline overrides this flag for the Hirshfeld run: it logs "H normalization was disabled for 3D view — but Hirshfeld REQUIRES normalized H. Overriding." and proceeds with normalization. X-ray H is physically incompatible with a correct Hirshfeld surface, so the user's global 3D-view preference is temporarily overridden for this specific analysis. The user's Reset is not lost — the flag is restored afterwards.
• NeedsNormalization — the normal path. The pipeline calls Molecule.ApplyHydrogenNormalization, which

1. Finds the nearest heavy neighbor for every H. The search uses a 1.3 Å cutoff and a minimum-image fractional difference , so H and heavy atoms in different unit cells are matched correctly.
2. Skips bridging hydrogens (H that is within the cutoff of more than one heavy atom at the same time, e.g. a hydrogen bond donor that also sits near a second acceptor). Skipping avoids forcing a wrong bond direction on a shared H. Skipped H are reported with reason "bridging".
3. Applies a Cartesian delta: cartShift = unit · (target − len), where unit is the normalized direction from the heavy atom to the current H. This shift is then converted to a fractional delta via Molecule.CartDeltaToFracDelta and added to the H's fractional coordinates. The H stays in the same copy of the unit cell as its heavy neighbor (the fractional coordinates are never wrapped), so cross-boundary molecules remain intact.
4. Rebuilds the whole 3D scene from the modified asymmetric part via Molecule.RebuildFromAsymetricPart. This regenerates the symmetry equivalents, the 3D atom array, the bond graph, the hydrogen bond graph, and the MoleculeShift.

The pipeline then closes any open fingerprint window (its per-vertex data is now stale), clears the previous Hirshfeld surface (it was built on X-ray H), and continues with the normal detection algorithm.

7.3. Hirshfeld Surface Properties

• Hirshfeld Surface (d_i, d_e, d_norm): A surface that represents the boundary of the molecule's space in the crystal. It is defined by the w(r) = 0.5 isosurface of the Hirshfeld weight function, which is the ratio of the electron density of a promolecule (the selected molecule) to the total electron density of the crystal (procrystal). The surface can be colored by various properties.
  • d_i: The distance from the surface point to the nearest nucleus inside the surface.
  • d_e: The distance from the surface point to the nearest nucleus outside the surface.
  • d_norm: A normalized contact distance, calculated as d_norm = (d_i − r_ivdW) / r_ivdW + (d_e − r_evdW) / r_evdW . Areas with d_norm < 0 (typically colored red) are regions of strong interaction (e.g., hydrogen bonds). Blue areas represent longer, weaker contacts.
• Charge: The average partial charge at each vertex, computed as a Gaussian-weighted average over all nearby atoms. Uses the per-atom charge field from the source file. The color ramp matches the CrystalExplorer convention: green (negative), yellow (neutral), red (positive). This is a powerful visualization for identifying electrophilic and nucleophilic regions on the surface.
• Hydrophobicity: The average hydrophobicity at each vertex, computed as a Gaussian-weighted average over all nearby atoms. Uses the per-atom Hydrophobicity field from the source file. The color ramp is: blue (hydrophilic), white (neutral), brown (hydrophobic). Useful for studying host–guest chemistry and protein–ligand interactions.

7.4. Electron Density Model — Slater-Type Orbitals (STO)

The promolecule density ρ_target(r) and the surrounding density ρ_outer(r) are computed from a Slater-Type Orbital basis. Each atom in the crystal is represented by the sum of the squared radial orbitals of its electronic configuration, with the coefficients and exponents taken from the tabulated Roothaan–Hartree–Fock STO basis of Koga and coworkers.

STO radial function

For an orbital with principal quantum number n and exponent ξ, the normalized radial part is

R_n(r) = N · r^(n-1) · exp(−ξ·r)
N     = (2ξ)^(n+1/2) / sqrt( (2n)! )

The spherical density of a single shell is the sum of the squared contractions over all orbitals of that shell:

ρ_shell(r) = Σ_orbital  occ(orbital) · | Σ_term  c_term · R_term(r) |²

where occ is the occupancy of the orbital in the atom's configuration (alpha and beta spin summed). The total density of the atom is

ρ_atom(r) = (ρ_s + ρ_p + ρ_d + ρ_f) / (4π)

Dataset — Koga Roothaan–Hartree–Fock

Tal ships with three JSON files that contain the complete STO basis for every element that appears in the neutral, cation and anion datasets. The files are located next to the executable (or in one of the search paths listed in StoDatabase.EnsureLoaded):

Each JSON entry contains the atom symbol, total, kinetic and potential energy, the config string (e.g. "CARBON 1S(2)2S(2)2P(2), 3P"), and the four shells s, p, d, f, each of which is a list of terms {n, ξ, coefs[]}. The shell also carries the array of per-orbital occupancies OccAlpha and OccBeta (parsed on the fly by ElectronConfigParser).

Density LUT

Because evaluating the sum of squared contractions at every voxel is expensive, the density of each symbol is precomputed once per session into a StoDensityLut. The LUT is a 1-D array of ρ(r) sampled on a uniform grid with step 0.03 Å up to a hard limit of 15 Å. Three global parameters control the size and shape of the LUT:

These three values can be changed at runtime from the Advanced panel (see 7.5). Any change requires a rebuild of the LUT cache, which is done via the "Rebuild LUT cache & Regenerate" button. Because the LUT is shared by all atoms of the same symbol, the cache must be cleared before a regeneration picks up the new parameters.

7.5. Advanced Panel — STO Tuning

The right-side control panel contains an Advanced checkbox. Ticking it reveals a small group with three numeric up-down controls and a rebuild button:

7.6. 2D Fingerprint Plots

This is the cornerstone of Hirshfeld analysis. The fingerprint plot is a 2D histogram of the d_i versus d_e values for every point on the Hirshfeld surface.

• Each point on the plot represents a point on the 3D surface.
• Areas with a high density of points (shown as "hot" colors) represent a large fraction of the total surface area involved in those specific types of contacts.
• The plot can be filtered by the chemical element of the atom outside the surface, which allows you to see the contribution of specific interactions, like O···H, C···H, or H···H contacts, to the overall crystal packing.
• The shape and position of the features on the plot provide a detailed signature of the intermolecular interaction network.

7.7. Hirshfeld Window (HirshfeldDocument)

The main 3D view for the Hirshfeld analysis is the HirshfeldDocument — a DockContent that hosts a MonoGame viewport and a right-side control panel. It is created from the Hirshfeld Surface Analysis menu item and reused on subsequent invocations.

Layout

• Left: 3D viewport (_panel, fill). Renders the molecule (or isolated target), the Hirshfeld isosurface, optional wireframe, axes and unit cell. Background color follows Flags.Instance.color1.
• Right: control panel (_controlsHost, 300 px wide, docked Right). A FlowLayoutPanel with all controls, sorted top-down.

Control panel — Target selection

Control panel — Detection parameters

Control panel — Property & display

Control panel — Actions

Progress bar & status label

• _progressBar: ProgressBarStyle.Continuous, 0 – 100. Updated through a thread-safe BeginInvoke adapter.
• _lblStatus: 4-line status text (monospace). Shows "Ready", "Generating…", final vertex/triangle counts, and the current blur state.

GenerateAsync — nine-step algorithm

1. STEP 0/9 — Hydrogen normalization. NormalizeHydrogensForHirshfeldSilently runs first (see 7.2). If any H needed normalization, the pipeline rebuilds the 3D scene, closes the fingerprint window, and clears the previous surface before proceeding.
2. STEP 1/9 — Building HirshfeldScene. Calls HirshfeldScene.Build(targetMolId, neighborCutoff). If the target molecule is not found, the operation aborts.
3. STEP 2/9 — Preparing atom entries. Each SceneAtom is wrapped either as an AtomStoParams (if the STO database is loaded and a LUT exists for the element) or, as a fallback, as an AtomGaussParams via HirshfeldParameters.BuildFromAtomSimple. The log line "Total: STO=N, Gaussian=M" shows how many atoms went through each path.
4. STEP 3/9 — Building density grid. A cubic grid is created around the target center. Resolution n = ceil(2·(r_max + margin) / voxelSize). Fails gracefully if n³ exceeds Flags.MAX_VOXELS.
5. STEP 4/9 — Filling density grid. Parallel per-thread accumulation of density contributions. Uses all available processor cores (not ProcessorCount − 1); the atoms are distributed round-robin across the threads so that load is even. The hot loop calls StoDensityLut.SampleSq for STO atoms and HirshfeldExpLut.FastExpNeg for the Gaussian fallback.
6. STEP 5/9 — Computing w(r). In-place reduction: W = RhoPro / (RhoPro + RhoProm_Outer) with a strict fallback of 0 for empty voxels. The distribution is written to the debug log as w(r) [after compute].
7. STEP 6/9 — Marching Cubes. Calls HirshfeldMarchingCubes.Extract, which delegates to HirshfeldMarchingCubesGenerator.Extract. It reads the flat W array, subtracts isoLevel inline, and produces a triangle mesh. Duplicate vertices are merged with a 27-neighbour hash grid.
8. STEP 7/9 — Metadata. Populates TargetMoleculeId, TargetAtomCount, OuterAtomCount, Center, UsedAnisotropic, and calls ComputeBoundingRadius().
9. STEP 8/9 — Per-vertex properties. HirshfeldProperties.Compute populates DE, DI, DNorm, Charge, Hydrophobicity, and stores the nearest-atom indices and names. The original values are also cloned into DEOriginal, DIOriginal, DNormOriginal.
10. STEP 9/9 — Post-processing. Wireframe buffer, axes buffer, cell buffer, and view matrix are rebuilt; the MoleculeShift sentinel is updated; the fingerprint window (if open) is refreshed with the new analysis.

Mouse & keyboard

• Left drag — rotate the camera. The scene sorting is invalidated and triangles are re-sorted next frame.
• Mouse wheel — zoom the whole scene (multiplicative 1.05 / 0.95, clamped 0.2 – 8.0). This affects _zoomFactor (scene scale), not the camera zoom.
• Ctrl+S — triggers SavePng().

MoleculeShift compensation

The surface vertices are stored in "raw" Cartesian coordinates (no −MoleculeShift). At render time the document subtracts the current MoleculeShift from each vertex, so the surface always aligns with the atoms rendered by the main viewport. When the user changes the presentation mode (and hence the shift), the buffers are rebuilt automatically — the change is detected via the _lastMoleculeShift sentinel field.

7.8. Hirshfeld Pipeline — Computational Classes

The Hirshfeld surface is produced by a self-contained pipeline that does not depend on the main document's presentation state. It is built on top of the same Molecule data used by the rest of Tal, but with its own cell reconstruction step.

HirshfeldScene

Central coordinator of the pipeline. Given a Molecule and a target molecule ID, it reconstructs a compact unit cell, runs the molecule finder, marks the host atoms and collects all outer (guest) atoms within a neighbour cutoff.

• Inputs: target Molecule, targetMoleculeId, cutoff (Å). Default cutoff = HirshfeldParameters.CUT_OFF = 10 Å.
• Public properties: TargetAtoms, OuterAtoms, Fragments, TargetCenter, TargetBBoxMin/Max, TargetRadius, TargetMoleculeId, LastCutoff, LastNeighborCount.
• SceneAtom struct: stores SymmetryAtomIndex, CellShiftX/Y/Z, FracPosition, CartPosition, AtomicNumber, IsTarget, IsAnisotropic, Un (ADP matrix), Occupancy, Name, MoleculeId, SourceAsymIndex, and (v14) Charge and Hydrophobicity.

HirshfeldUnitCell

Independent representation of the unit cell, built exclusively from Molecule.AsymetricPart and Molecule.SymmetryOperations. It deliberately does not read Molecule.SymetryEquivalents or Molecule.atom[], because those are mutated by the presentation-mode code.

1. For each atom of the asymmetric part and each symmetry operation, generate a raw fractional position.
2. Deduplicate against the growing list with a fractional tolerance DEDUP_TOL_FRAC = 0.005.
3. For each symmetry operation ops[k], k ≠ 0, find the anchor atom (first non-hydrogen of ops[0]) and its counterpart in ops[k]. Compute a single integer translation from the anchor to its counterpart, then apply that same translation to all atoms of ops[k]. This keeps symmetry-generated molecules intact.
4. Deduplicate again to remove post-translation collisions.
5. v14: the per-atom Charge and Hydrophobicity values are copied from AtomSimple so they propagate through the whole pipeline.

HirshfeldDetector

Detects the surface and returns a HirshfeldSurface. Its main entry point is GenerateAsync, which drives the nine-step algorithm described above. Important public properties:

HirshfeldGrid

Wraps a VoxelScanner.VoxelGrid and holds flat density arrays:

float[] W;                // w(r) = rho_pro / (rho_pro + rho_outer)
float[] RhoPro;           // promolecule (target)
float[] RhoProm_Outer;    // surrounding atoms

Indexing is idx = i + Nx·(j + Ny·k) — the deepest loop index k gives sequential memory access. ComputeRelativeDensity is parallelized and iterates linearly, with the historical fallback of 0 for empty voxels.

HirshfeldParameters & HirshfeldAtomicRadii

HirshfeldParameters.BuildFromAtomSimple converts an atom into an AtomGaussParams:

• Isotropic: ρ(r) = Z_eff / (√(2π)σ)³ · exp(−|r−r₀|² / (2σ²)). The per-element σ comes from HirshfeldAtomicRadii.GetSigma(Z), multiplied by the user's Sigma multiplier.
• Anisotropic: the ADP matrix Un is row-normalized to get the three eigenvectors; the row lengths before normalization are the principal half-axes. Each half-axis is combined with the electronic σ as a_k = sqrt(λ_k² + σ_e²).
• Cutoff radius: CUT_OFF × max(σ or a_k) with CUT_OFF = 10.

This class is only used by the Gaussian fallback path. When the STO database is loaded and contains an entry for the element, the density comes from the LUT instead, and none of these parameters are involved.

HirshfeldExpLut

Hot loop replacement for MathF.Exp(-x). A static LUT of 4096 entries over x ∈ [0, 25] with linear interpolation. Maximum relative error around 1e-4, invisible in the final rendering. Used only by the Gaussian fallback path.

HirshfeldMarchingCubes & HirshfeldMarchingCubesGenerator

A dedicated copy of Tal's standard Marching Cubes generator, adapted for flat density arrays and Hirshfeld iso conventions. It never modifies the input array; the iso offset is subtracted inline. Duplicate vertices are merged with a 27-neighbour spatial hash with tolerance 0.0001.

HirshfeldProperties

For each mesh vertex, computes:

d_e     = min distance to any outer atom
d_i     = min distance to any target atom
d_norm  = soft-min average of (d_e − r_vdW) / r_vdW  +  (d_i − r_vdW) / r_vdW
charge  = gaussian-weighted average of atom.charge
hydro   = gaussian-weighted average of atom.Hydrophobicity

The d_norm value uses a soft-min over the K = 8 nearest neighbours with a Gaussian weight exp(−β·d) and β = 8. This removes the Voronoi discontinuities that a hard nearest-neighbour approach produces on the surface. The result is clamped to [−3, +3].

The charge and hydrophobicity values use a plain Gaussian average over all target and outer atoms with σ = 0.6 Å; the outer atoms are weighted by an extra factor 0.7. If the source file does not provide these per-atom values, the corresponding surface field is uniform and a warning is logged.

HirshfeldSmoothing

A geometry-only Taubin (λ | μ) smoothing applied to the mesh itself. Unlike a plain Laplacian, the negative μ step compensates for shrinkage, so the surface volume is approximately preserved. The original vertex positions and the adjacency list are cached per HirshfeldSurface. In the current version this class is not wired to a UI slider — it remains available for programmatic use.

HirshfeldMeshRelaxation

An optional pre-color smoothing step that moves vertices tangentially along the mesh and slightly along their normals to remove the "staircase" artefacts of Marching Cubes without losing the physical shape. The maximum displacement is clamped to 5% of the grid unit to avoid overshooting. In the current version this class is not wired to a UI slider — it remains available for programmatic use.

HirshfeldTargetSelector

Simple helper that lists all molecules in a Molecule by their IndividualMoleculeInCell ID. Used by the _cmbMolecule dropdown and by UseCurrentSelection().

7.9. Fingerprint 2D Window (HirshfeldFingerprintWindow)

A separate DockContent that renders a 2D heat map of the d_e × d_i distribution. Each bin corresponds to a small angular sector of the surface; the number of vertices in a bin is mapped to a log-scaled color ramp.

Layout

• Left: 2D plot panel (_panel, fill). The heat map, axes, grid, diagonal and highlighted bins are drawn here with a SpriteBatch in PointClamp mode.
• Right: control panel (_controlsHost, 360 px wide, docked Right). A FlowLayoutPanel with all controls, sorted top-down.

Side panel — Plot parameters

Side panel — Contact decompose

This is the essential feature that makes the fingerprint usable for real chemistry. It consists of a CheckedListBox that lists every contact type present on the surface, a pair of buttons for mass (de)selection, a read-only TextBox with the full decompose table, and a Copy table button.

Contact   Area (Ų)      %    Verts
C···H         245.3    18.4      9823
O···H         185.7    13.9      7412
H···H         412.5    30.9     16840
...

Bin construction & re-bucketing

• Default range [0.0 Å, 3.4 Å], default bin size 0.005 Å. Both can be changed live from the side panel.
• RebuildHeatmapFromRaw() recomputes a full histogram and a per-vertex (binDe, binDi) table. Vertices outside the range are tagged (−1, −1).
• The new per-vertex bins are sent back to the HirshfeldDocument via the OnVertexBinsUpdated callback, so the 3D highlighting stays in sync even after re-bucketing.
• Range and bin are clamped to safe limits: binsInRange ≤ MAX_TOTAL_CELLS = 8192 and texture bytes ≤ 512 MiB.

Rendering

• Background: Flags.Instance.color1 (same as the 3D viewport).
• Heat map: single Texture2D of size texSize × texSize, drawn with SamplerState.PointClamp (crisp bins).
• Color ramp: 4-segment palette white → cyan → blue → magenta/red, driven by log(1 + count) / log(1 + maxCount) and an intensity^0.6 compression. Minimum visible intensity is 0.15 so isolated bins are still visible.
• Optional overlays: major grid, minor grid, diagonal d_e = d_i line, axes and labels for d_e (Å) and d_i (Å).

Mouse interaction

Highlight sync

When a bin is selected the window calls OnRegionClicked(binDe, binDi). The HirshfeldDocument stores _highlightedBin and rebuilds the vertex buffer. Every vertex whose (binDe, binDi) matches is colored bright yellow (255, 255, 0, alpha). When the user changes the bin size or the range, the new bins are pushed back through OnVertexBinsUpdated so the highlight remains correct.

Save PNG

_btnSavePng renders the same content into a 2048 × 2048 RenderTarget2D and saves it as PNG. The label offsets, axis ranges, and highlight state are preserved.

7.10. Surface Data Model and Diagnostics

HirshfeldSurface

The data container for a computed surface. All vertex-based arrays share the same length.

HirshfeldAnalysisResults

A compact result set that is passed to the fingerprint window. It contains:

• FingerprintHistogram — Dictionary<(int binDe, int binDi), int>.
• FingerprintBinSize, FingerprintRangeMin, FingerprintRangeMax.
• ContactTypes — a string[] of length VertexCount with the type of each contact (e.g. "C···H").
• ContactTypeCounts — a Dictionary<string, int> of vertex counts per contact type.
• ContactTypeAreas — a Dictionary<string, float> of total surface area per contact type.
• ContactTypePercentages — a Dictionary<string, float> of area percentages per contact type.
• TotalSurfaceArea — total surface area in Ų.
• DecomposeTableText — the pre-formatted text table ready for copy/paste.
• VertexBins — a per-vertex (binDe, binDi) table.

HirshfeldAtomDistanceDiagnostics

A diagnostic object built after each generation. It contains:

• Target and outer atom count, min/max/median/p95/p99 of distances from the target center.
• Grid metadata: TargetCenter, RTarget, RGrid, Margin, VoxelSize, GridResX/Y/Z, TotalAtoms.
• A histogram (bin size 0.5 Å) of outer atom distances.
• A cutoff analysis table showing, for several R_target + delta thresholds (2–8 Å), how many atoms would be rejected.

The diagnostic is shown in the DebugWindow via the Show Diagnostics in Debug… button. The formatting is handled by HirshfeldDiagnosticsFormatter, which builds a list of (text, System.Drawing.Color) lines.

Export formats

• PLY — ASCII, one line per vertex (x y z nx ny nz r g b), then a face list. Compatible with MeshLab, Blender, CloudCompare.
• PNG — a high-resolution screenshot rendered with SSAA = 2.

8. Powder Diffraction

The Powder Diffraction window simulates the theoretical powder X-ray diffraction pattern for the loaded crystal structure. It is opened from Calculate → Powder Pattern....

The underlying engine (PowderDiffractionEngine) enumerates allowed hkl reflections in the 2θ range, computes the structure factor F(hkl) from the full unit cell, applies Lorentz–polarization and Debye–Waller corrections, groups symmetry-equivalent reflections, and renders the resulting profile in a WinForms plot.

8.1. Window Layout

The window is a DockContent composed of a large plot panel (panelPlot, a BufferedPanel) and a bottom control strip (panelControls) divided into group boxes.

8.2. Controls and Parameters

Range 2θ (groupBoxRange)

Peak parameters (groupBoxPeak)

Debye-Waller (groupBoxDebye)

Radiation source (groupBoxSource)

Available sources (from XRaySource enum):

Display options (groupBoxDisplay)

Info and actions (groupBoxHover, buttons)

8.3. Mouse Interaction

• Wheel: Zooms the 2θ range about the cursor position (factor 0.8 in, 1.25 out).
• Left drag: Rubber-band zoom. The selected pixel range is converted into a new 2θ range.
• Double-click: Resets the range to 0–90°.
• Hover: Updates the readout labels. The "closest peak" tolerance is 2% of the current 2θ span.

8.4. Algorithm and Formulas

1. d-spacing: Computed from the metric tensor of the triclinic cell using the standard reciprocal-lattice formula.
2. Angle: 2θ = 2 · asin(λ / (2d)), rejected if sinθ > 1.
3. Centering prefilter: Determined from the first letter of the actual Hermann–Mauguin symbol (P/A/B/C/F/I/R), not from the space-group number. This is only a fast performance filter — systematic absences from screw axes and glide planes are handled naturally by F(hkl).
4. Structure factor: F(hkl) = Σ_j f_j(s) · exp(−B·s²) · occ_j · exp(2πi(h·x_j + k·y_j + l·z_j)), summed over all atoms in the unit cell (using SymetryEquivalents if available). The atomic scattering factors f_j(s) use the Cromer–Mann 4-Gaussian + constant approximation (International Tables Vol. C, Table 6.1.1.4).
5. Lorentz-polarization correction: L·P = (1 / (2 sin²θ cosθ)) · ((1 + cos²2θ) / 2).
6. Intensity: I = |F|² · L·P.
7. Grouping: Reflections with the same rounded 2θ (4 decimal places) are merged; intensities are summed, multiplicities added. This is valid for every crystal system without hard-coded symmetry tables.
8. Profile: The rendered profile is a sum of pseudo-Voigt profiles centered at each peak 2θ_k with the chosen FWHM:
   • Gaussian: exp(−4 ln 2 · x²) where x = (2θ − 2θ_k) / FWHM.
   • Lorentzian: 1 / (1 + 4x²).
   • Pseudo-Voigt: η · G + (1 − η) · L with η = PseudoVoigtEta (default 0.5).

8.5. Export

The ExportData() dialog offers two paths:

• 📊 CSV data — one row per peak: 2θ(°), d(Å), Q(Å⁻¹), h, k, l, Intensity, Multiplicity, |F|. Written as UTF-8 with a leading comment line.
• 🖼️ Image — saves the current plot bitmap as PNG, JPEG (quality 95, via EncoderParameters), or BMP.

btnCopy places the bitmap on the clipboard as both a DataFormats.Bitmap and a raw PNG stream.

8.6. Settings Persistence

The PowderDiffractionSettings class holds all numeric and boolean values as properties with defaults defined in constants (single source of truth). Clone() and ResetToDefaults() are available for programmatic use. The window synchronizes controls to settings via LoadSettingsToUI() and SaveUIToSettings().

9. Crystal Planes & Mean Planes

Tal can define, edit and visualize two kinds of crystallographic planes:

• Mean plane — the least-squares best-fit plane through a set of selected atoms.
• hkl plane — a crystallographic plane defined by its Miller indices.

9.1. Planes Window (list manager)

Opened from Calculate → Planes... or from the main Planes toolbar button. It is a DockContent named PlanesWindow, docked left/right or floatable.

Controls

The window maintains a color icon cache (_colorImageCache) keyed by ARGB value, so repeated planes of the same color reuse the same ImageList index.

9.2. Plane Properties Window

PlanePropertiesWindow is a modal-ish dialog (border fixed, center-parent) with a tab control for the two plane types and a common options group.

Tab: Mean plane

• Pick atom to select: three radio buttons — Picked atom (default), Picked molecule, Ring.
• Or select: three buttons — All, Non-hydrogen, Hydrogen.
• Selection info: a status label that shows the current count. If fewer than 3 atoms are selected the label turns firebrick and shows "Selected: N atoms (need at least 3)".

When the Mean plane tab is active, the window enters plane selection mode in the main form. A 300 ms timer monitors the selection count and refreshes the label. Atom clicks are routed through HandleAtomPick(int atomIndex) according to the current selection mode:

• Picked atom — toggles the single atom's selected flag and updates SelectionOrder and NumberSelected.
• Picked molecule — selects all atoms with the same IndividualMoleculeInCell.
• Ring — calls RingFinder.FindSmallestRing() with parameters maxRingSize = 12, requirePlanar = true, planarityTolerance = 0.15 Å. If no ring is found, falls back to selecting the single clicked atom.

Tab: hkl

• Numeric up-down controls for h, k, l. Range −20 to +20. Default all zero; if all three are zero the dialog refuses to save and shows a warning.

Common options (groupBoxCommon)

On OK, SaveUIToPlane() copies the UI state to the working plane and calls CalculatePlaneGeometry(). The result is exposed via the Result property.

9.3. Mean Plane Calculation

The least-squares mean plane is computed by MeanPlaneCalculator.CalculateMeanPlane():

1. Collect the selected atom positions in raw Cartesian space. Because Molecule.atom[i].Coordinate is stored with −MoleculeShift applied, the shift is added back so the plane is independent of re-centering.
2. Compute the centroid.
3. Build the 3×3 symmetric covariance matrix of the centered positions.
4. Eigendecompose it with the Jacobi algorithm (max 50 sweeps, ε = 1e-12). Eigenvectors are returned as columns.
5. The eigenvector with the smallest eigenvalue is the plane normal.
6. Corners are built by PlaneGeometry.BuildRectangle(center, normal, size), which constructs an orthonormal (u, v) basis inside the plane and returns the four corners center ± u·size ± v·size.

The plane equation coefficients A, B, C, D are stored (Ax + By + Cz + D = 0). The corners are stored in raw Cartesian (without shift), so the draw layer subtracts the current shift at render time.

9.4. hkl Plane Calculation

PlaneGeometry.CalculateHKLPlane() uses the molecule's fractional-to-Cartesian matrix (transposed, to match Molecule.PrepareUnitCell) to convert intercept points into Cartesian space. The intercepts depend on how many Miller indices are zero:

• All zero — invalid, empty result.
• One non-zero — the plane is a cell face; returns the 4 corners of that face at x = 1/h, y = 1/k, or z = 1/l.
• Two non-zero — the plane is parallel to the axis with the zero index; returns 4 corners that sweep the cut edge along that axis.
• Three non-zero — general case; returns the 3 intercepts with the axes, and a rectangle is built around the centroid with the chosen size (default 10 Å).

9.5. Plane Intersection

PlaneGeometry.ComputeIntersection(p1, p2) returns a PlaneIntersection struct with a point on the line, the normalized direction, and a half-length. The algorithm:

1. Direction = n₁ × n₂ (normalized). If nearly zero, the planes are parallel and the result is invalid.
2. A starting point is found by two Newton-style projections (onto plane 1, then along the projected n₂, then again onto plane 1).
3. For each plane, the line (point + t·dir) is projected into the plane's local (u, v) basis. The u/v extents are computed by scanning all 4 corners, so any corner ordering works. The permitted t-interval is the intersection of the u and v slab constraints.
4. The final point is at the midpoint of the interval; the half-length is half the interval length (minimum 0.1 Å). A bounding-box fallback is used if the precise computation fails.

9.6. Ring Finder

RingFinder finds the smallest ring containing a clicked atom:

1. For each bonded neighbor of the clicked atom, run a BFS back to the clicked atom.
2. The BFS forbids passing through the start atom as an intermediate and forbids the direct neighbor → start edge (which would give a trivial 2-ring).
3. Candidate rings shorter than maxRingSize and free of duplicate atoms are tested for planarity by PCA (smallest eigenvalue of the covariance matrix).
4. The smallest planar ring is returned. If none is planar, the smallest non-planar ring is not returned — the caller falls back to single-atom selection.

9.7. CrystalPlane Data Model

10. Morphology

Opened from the Advanced menu, this tool is used for analyzing and visualizing the macroscopic crystal habit (morphology) based on the internal crystal structure.

The morphology is predicted using the Bravais-Friedel-Donnay-Harker (BFDH) method, which relates the importance of a crystal face to its interplanar spacing (dhkl) and the lattice geometry. The tool allows you to:

• Generate a 3D polyhedron representing the predicted crystal habit.
• Color the faces based on their Miller index (hkl) or other properties.
• Overlay the predicted morphology with the 3D structure of the molecule.

10.1. Morphology Window

The main control window is the MorphologyWindow, a DockContent that can be docked left, right, or floated. It is opened from Advanced → Crystal Morphology.... The window is divided into three vertical zones:

• Top panel — calculation parameters and action buttons.
• Info panel — a summary of the predicted habitus and unit cell.
• Face table — an owner-drawn ListView listing every face that survived the BFDH filtering.

10.2. Calculation Parameters

All parameters are exposed as WinForms controls in the top panel:

10.3. Face Table

The face table is a ListView with View.Details, FullRowSelect, and custom owner-drawing (DrawColumnHeader / DrawItem / DrawSubItem). The columns are:

Clicking a row highlights the corresponding face in the 3D view (all triangles with the same FaceName get a yellow outline). This is implemented via MorphologyDocument.HighlightFace(string faceName), which rebuilds a dedicated VertexBuffer (_highlightVB) containing the face's triangle edges.

10.4. The BFDH Algorithm

The BFDH (Bravais–Friedel–Donnay–Harker) method predicts crystal morphology from the crystal structure by assuming that the growth rate of a face (hkl) is inversely proportional to its interplanar spacing d_hkl:

RelativeGrowthRate(hkl) ∝ 1 / d(hkl)

Faces with large d_hkl grow slowly and therefore dominate the final crystal habit. The algorithm implemented in CrystalMorphologyFinder proceeds as follows:

1. Plane enumeration: All integer triples (h, k, l) with |h|, |k|, |l| ≤ maxIndex are generated. The trivial (0,0,0) is rejected, as are all triples whose GCD is greater than 1 (they represent the same plane as the primitive triple). This is done by GenerateUniquePlanes.
2. d-spacing calculation: For each candidate, the d-spacing is computed from the triclinic metric tensor by CalculateDSpacing. Planes with d < minDSpacing are discarded.
3. Multiplicity calculation: The real multiplicity is computed from the crystal's space group by CountEquivalentPlanesReal, which builds the orbit of (hkl) under the linear parts of all symmetry operations and counts distinct canonical forms. If the space group is unavailable or invalid, the code falls back to a cubic approximation (CountEquivalentPlanesCubic), which is clearly marked as a fallback.
4. Relative growth rate: RelativeGrowthRate = 1 / d.
5. DistFrac computation: DistFrac is defined as RelativeGrowthRate / min(RelativeGrowthRate), so the largest face has DistFrac = 1.
6. Softening (optional): The DistFrac values are compressed according to the selected SofteningMode (see section 10.5).
7. Relative area: RelativeArea = 1 / DistFrac, normalized to the maximum. Note that this is a visual heuristic — the code explicitly does not multiply by multiplicity, because the physical area of a face depends on the depth of the cut (DistFrac), not on the number of symmetry-equivalent planes.
8. Sorting: The face list is sorted by DistFrac ascending, so the largest faces come first.
9. Habitus classification: ClassifyHabitus inspects the relative areas of the {100}, {010}, {001}, {110}, {101}, {011}, and {111} families and returns one of: cubic / octahedral, plate, needle / prismatic, rhombic / dodecahedral, equant (isometric), prismatic, or unknown.

10.5. Softening Modes

Pure BFDH (SofteningMode.None) is physically correct but often produces crystals with many very small, almost invisible faces. The softening modes are visual heuristics that compress the DistFrac distribution to produce more aesthetically pleasing or recognizable shapes. They are not physically based models and should be used with that in mind.

When the mode is changed in the combo box, the numeric up-down is automatically updated to the suggested strength for that mode (via GetSuggestedSoftening). The user can of course override this value.

10.6. Morphology 3D View (MorphologyDocument)

The 3D crystal is displayed in a MorphologyDocument, a DockContent that hosts a MonoGame viewport. It is created automatically the first time "Calculate" is clicked in the Morphology window, and reused on subsequent clicks (the title is updated to reflect the current structure name).

Construction

The polyhedron is built by MorphologyCrystal.Build:

1. Start with a large cube (ConvexPolyhedronFrac.FromBox) of half-edge 3.0 · scale Å in fractional space.
2. Sort the surviving faces by DistFrac ascending (largest faces first).
3. For each face, compute its unit normal in Cartesian space and the scaled plane offset dScaled = DistFrac · ‖g‖ · faceOffset, where ‖g‖ is the norm of the reciprocal-lattice vector of (hkl).
4. Clip the current polyhedron by the half-space defined by this plane (ConvexPolyhedronFrac.ClipByPlane). The clipping is a Sutherland–Hodgman-style operation on each triangle, followed by a cap construction from the intersection points (a fan triangulation of the cap polygon, oriented to match the plane normal).
5. If a clip would destroy the polyhedron entirely, it is rolled back. A maximum of MaxSuccessfulClips = 14 successful clips is enforced to keep the polyhedron manageable.
6. The resulting fractional-space triangles are converted to Cartesian space using the molecule's fractional-to-Cartesian matrix (transposed to match Molecule.PrepareUnitCell).
7. Face colors are assigned deterministically from the face name (NameToColor) using an FNV-1a hash mapped to HSV. The starting cube's "cell-*" faces are rendered in neutral grey.

Rendering and interaction

• Camera: a shared Camera instance; the default position is (0, 0, 100) looking at the origin.
• Rotation: left mouse drag. The rotation matrix is exposed via GetRotationMatrix() so the Wulff plot can follow it.
• Zoom: mouse wheel. The zoom factor is clamped to 0.2–5.0 and applied on top of an automatic scene scale (40 / BoundingRadius) that keeps the crystal nicely framed.
• Reset View: restores the default camera position, identity rotation, and 1× zoom.
• Projection: switches between orthographic (default) and perspective based on Flags.Instance.PerspectiveMode.

Side panel controls

• Alpha (0–255): sets the transparency of the crystal faces. The alpha value is applied to every triangle's color at render time (a TrackBar with a live numeric label).
• Show Edges: toggles the drawing of black edges along the boundaries between faces. The edge buffer is rebuilt whenever the visibility, alpha, or crystal changes.
• Show Face Labels: toggles the drawing of (hkl) labels at the centroid of each face. Uses the shared SpriteFont1 and respects Flags.LabelScale, Flags.LabelColor, Flags.LabelShadow, and Flags.LabelShadowColor. Faces with duplicate names share a single label at their combined centroid.
• Monochrome / Color swatch: enables a single user-selected color for all faces instead of the per-face deterministic colors. Clicking the swatch opens a ColorDialog.
• Reset View: resets the camera and the zoom factor.
• Zoom info label: shows the current zoom factor in the form "Zoom: X.XX× (mouse wheel)".

Edge detection algorithm

Edges are drawn only along the true boundaries between two different faces. The algorithm is deliberately careful about numerical robustness:

1. All triangle vertices are merged into unique positions using a 27-neighbour hash grid with a tolerance EPS_world = BoundingRadius · 1e-4 (in the crystal's original scale, not after any zoom or scale transform).
2. The set of unique raw edges (pairs of vertex indices) is built from all triangles.
3. T-junction removal: for each raw edge, all other unique vertices are tested for lying on the segment (parameter t ∈ (tEps, 1 − tEps), perpendicular distance < 4 · EPS_world). Vertices that pass are sorted by t and the edge is split into a chain. This removes the visual gaps that appear at T-junctions where two triangles of the same face share a boundary with a third triangle of a different face.
4. For each triangle, the chain is walked and each sub-segment is tagged with the face ID of the triangle. If a sub-segment ends up tagged with two or more distinct face IDs, it is drawn once as a boundary edge.
5. Face IDs are computed by GetFaceId: if the triangle has a non-empty FaceName, that name is used; otherwise a synthetic ID is built from the sign of the triangle normal.

Face highlighting

When the user selects a row in the face table, MorphologyWindow calls MorphologyDocument.HighlightFace(faceName) with the selected face name (or null to clear). This rebuilds the highlight vertex buffer by:

1. Collecting all triangles whose FaceName equals the requested name.
2. Merging their vertices into a small local hash grid (tolerance derived from the bounding radius and the current zoom).
3. Emitting a yellow (255, 255, 0, 255) line segment for every edge that appears in exactly one triangle of the selected face.
4. The highlight is drawn with DepthStencilState.None so it always appears on top of the crystal, regardless of depth ordering.

10.7. Wulff Plot (Stereographic Projection)

The Wulff plot is opened by clicking the Wulff Plot button in the Morphology window. It opens a WulffPlotDocument, a DockContent that renders a stereographic projection of the crystal faces. This is the standard graphical tool in crystallography for visualizing the angular relationships between faces and the symmetry of the crystal.

Stereographic projection

Each face is represented by its unit normal vector n in Cartesian space. The projection is defined with the +c axis at the center of the plot:

θ = acos(n.z)                    // polar angle from +c
φ = atan2(n.y, n.x)               // azimuthal angle in the ab plane
r = tan(θ / 2)                    // stereographic radius (gnomonic-scaled)
x_plot = r · cos(φ)
y_plot = r · sin(φ)

The factor tan(θ/2) is the standard Wulff-net stereographic projection (equal-angle). Faces close to +c project near the center; faces in the ab plane (θ = 90°) project onto the outer circle. Lower-hemisphere faces (n.z < 0) are mirrored onto the upper hemisphere for display (their z component is negated before projection), and are drawn as hollow circles; upper-hemisphere faces are drawn as filled circles.

Plot layout

• The plot area is always square-ish: the effective radius is min(screenW, screenH) / 2 · 0.85, multiplied by the current zoom factor. This keeps angular relationships correct regardless of the window aspect ratio.
• The outer circle (the θ = 90° equator) is drawn in light grey.
• An optional polar grid (every 10° of θ, plus radial lines every 30° of φ) is drawn in a dark blue-grey.
• An optional bold cyan equator is drawn on top of the grid when Show Equator is enabled.
• Each face is drawn as a small filled circle (upper hemisphere) or a ring (lower hemisphere). The size is controlled by the Point Size slider (4–24 px, default 10).
• Face labels (hkl strings) are drawn above each point using the shared SpriteFont1, respecting the global label settings from Flags.

Side panel controls

• Follow 3D Rotation: when checked, the plot rotates in sync with the 3D crystal view. The rotation matrix is obtained from MorphologyDocument.GetRotationMatrix() and multiplied by the same view matrix used by the 3D renderer (Matrix.CreateLookAt((0,0,100), Zero, Up)) so that the two views are visually consistent.
• Show Labels (hkl): toggles the (hkl) text labels.
• Show ALL faces: when unchecked (default), only faces present in the current 3D crystal polyhedron are shown. When checked, every face from the BFDH calculation is shown, including those that were clipped away during polyhedron construction. This is useful for understanding which planes the algorithm considered.
• Show Grid (10°): toggles the polar grid.
• Show Equator (θ=90°): toggles the bold equator circle.
• Color by d-spacing: when checked, replaces the deterministic per-face color with a rainbow gradient from red (large d) to blue (small d). The gradient is computed from t = d / maxD as (255·t, 100 + 100·(1 − |2t − 1|), 255·(1 − t)).
• Hemisphere group: three radio buttons — Upper only, Lower only, Both (default). Filters the projected points by the sign of n.z.
• Point Size: a slider (4–24, default 10) controlling the radius of the plotted points.
• Zoom info label + Reset Zoom: shows the current zoom factor; the reset button sets it back to 1.0.
• Refresh: forces a re-extraction of plot positions and a rebuild of the vertex buffers. Useful when the underlying crystal or its visibility has changed outside the Wulff plot.

Mouse interaction

The only mouse interaction is the wheel: scrolling up multiplies the zoom factor by 1.15, scrolling down by 0.87. The zoom factor is clamped to 0.3–5.0.

Vertex buffer strategy

Three separate VertexBuffers are maintained and rebuilt whenever the plot needs to be refreshed:

• _circleVB — the outer circle (128 segments, LineList).
• _gridVB — polar grid, radial lines, and equator (LineList).
• _pointsVB — the face markers, drawn as triangle fans (upper hemisphere) or triangle strips (lower hemisphere rings) (TriangleList).

The plot background is taken from Flags.Instance.color1, matching the 3D viewport so that screenshots and the on-screen plot have consistent backgrounds.

11. Structure Overlay

Opened from the Advanced menu, the 🔬 Structure Overlay window is a powerful tool for comparing multiple molecular structures.

Key Features:

• Add Structures: You can add any of the structures from your Structure List to the overlay.
• Manual Alignment: Each structure in the overlay can be independently translated and rotated using sliders, allowing you to manually superimpose them.
• Automatic Alignment: A powerful algorithm can automatically align the structures based on a common substructure or by minimizing the Root Mean Square Deviation (RMSD) between their atoms.
• Visualization: Each structure in the overlay can be given a unique color to easily distinguish them. You can also change their display style (e.g., show one as balls & sticks and another as a wireframe).

This tool is invaluable for comparing different conformations of a flexible molecule, or comparing a calculated structure with an experimental one.

11.1. Overlay Manager Window

The main control window is the MoleculeOverlayWindow (a DockContent named "🧬 Overlay Manager"), opened from Advanced → Overlay Manager.... It can be docked left, right, or floated, and holds an OverlayManager instance that is the shared data model for the whole overlay feature. The window is laid out as three vertical panels:

• Top panel — buttons for adding/removing molecules, the list of overlay entries, the visibility checkbox, and the atom counter.
• Middle panel — the ListBox of overlay entries (custom owner-drawn so each entry can be colored with the molecule's actual color).
• Bottom panel — the transformation radio buttons (None, Mirror, Inversion) and the status label.

11.2. The Overlay Data Model

Every entry in the overlay is an instance of OverlayMoleculeData, which holds:

The class exposes two important methods: InitializeWorkingData() (deep-copies OriginalAtoms into WorkingAtoms and stores the original ADP matrices) and ResetToOriginal() (reverses the process, restoring the working atoms from the originals while preserving the current color and style).

11.3. The Overlay Manager (OverlayManager.cs)

OverlayManager owns the list of OverlayMoleculeData entries and exposes a small API used by the window. The public surface is:

• IReadOnlyList<OverlayMoleculeData> OverlayMolecules — read-only view of the current entries.
• int TotalAtomCount — sum of WorkingAtoms.Count across all entries.
• event EventHandler MoleculesChanged — raised whenever the list or any entry's state changes, so the window and any overlay views can refresh.
• AddSelectedMolecule(doc, mode) — adds the whole molecule containing the selected atoms as a single entry.
• AddSelectedAtoms(doc, mode) — adds just the selected atoms as a single entry (marked IsSelection = true).
• AddSelectedStructure(doc, mode) — adds every atom of the document's molecule (including symmetry-generated ones) as a single entry (marked IsStructure = true).
• RemoveMolecule(id) — removes an entry by ID and raises MoleculesChanged.
• ClearOverlayMolecules() — removes all entries.
• ToggleVisibility(id) — flips the IsVisible flag.
• SetMoleculeColor(id, color) — sets the per-entry color and propagates it to every CustomColor field of the working atoms.
• SetMoleculeStyle(id, styleIndex) — sets the per-entry style (0 = CPK, 1 = Sticks, 2 = Balls and Sticks, 3 = Ellipsoids).
• SetTransformationForMolecule(id, mode) — the main transformation entry point (see below).
• FitMoleculeToPoints(source, sourcePoints, targetPoints) — the fitting entry point (see below).
• IsBaseMolecule(id) / GetBaseMoleculeId() — helpers for the "base molecule" concept (see section 11.6).

11.4. Adding Molecules — Workflow

The Overlay Manager exposes three "Add" buttons, each backed by a different method of OverlayManager. All three share the same underlying extraction logic, but differ in which atoms are collected:

1. Add Molecule (AddSelectedMolecule): collects all atoms of the molecule that contains the currently selected atoms. The molecule ID is taken from IndividualMoleculeInCell. If the source document has no molecule numbering (single-molecule file), all atoms are added as one entry. This is the most common case.
2. Add Selection (AddSelectedAtoms): collects only the currently selected atoms. This produces an entry marked IsSelection = true and is useful for comparing, e.g., just a ring system or a metal coordination sphere. Note: the current version of the Molecule Overlay window does not expose a button for this method — the button is commented out in the source, but the API remains available for programmatic use.
3. Add Structure (AddSelectedStructure): collects every atom of the source document's molecule, including all symmetry-generated atoms. This produces an entry marked IsStructure = true. Use this when comparing two whole crystal packings.

In all three cases, the first 3 selected atoms are stored as reference points. If fewer than 3 atoms are selected, the entry is still created, but HasReferencePoints is false and the automatic fitting (section 11.5) is skipped.

11.5. Automatic Fitting (3-Point Rigid-Body Superposition)

The core of the overlay feature is the automatic fitting algorithm, implemented as OverlayManager.IterativeFit(sourceMolecule, sourcePoints, targetPoints). Despite the name "iterative", the current implementation is a single-pass, closed-form rigid-body transformation that maps a source reference frame (defined by the 3 source points) onto a target reference frame (defined by the 3 target points).

The algorithm

1. Build a local orthonormal basis for the source points:

   xs = normalize(P1_source − P0_source)
   zs = normalize(cross(xs, P2_source − P0_source))
   ys = cross(zs, xs)   // already unit length

   The three vectors (xs, ys, zs) form a right-handed orthonormal frame anchored at P0_source.
2. Build the corresponding local basis for the target points using the exact same construction (same cross-product order, same sign conventions), anchored at P0_target. Because the construction is identical, the two bases have the same handedness, and the map between them is guaranteed to be a proper rotation (determinant +1, no reflection).
3. Compute the rigid-body transformation as a change of basis. For any atom at position p (in the source frame), its new position is:

   v = p − P0_source
   a = dot(v, xs), b = dot(v, ys), c = dot(v, zs)
   p_new = P0_target + a · xt + b · yt + c · zt

   This is the classic "re-express in the source basis, then rebuild in the target basis" trick. It implements a pure rotation + translation, with no scaling and no reflection. The three source points are mapped exactly onto the three target points by construction.
4. Apply the same mapping to every working atom and to the working reference points.
5. Recalculate bonds via RecalculateBondsForMolecule, which rebuilds the per-atom bond lists, distances, angles, and precomputed bond matrices from the new coordinates. This is necessary because the overlay atoms are rendered with the same bond-drawing code as regular molecules.
6. Report the fit quality: the initial and final sum of distances between the three source points and the three target points is logged, along with the improvement (in Å and as a percentage). For a successful rigid-body fit, the final error should be 0.0000 Å (to floating-point precision) — any non-zero residual indicates either that the three source points were collinear (in which case the fit is aborted with a "punkty współliniowe" error) or that the fit was called on a degenerate reference frame.

When is fitting triggered?

• Automatically on add: when a new molecule is added to an overlay that already contains at least one other entry, AddSelectedMolecule (and its siblings) calls AutoFitMolecule at the end. The new molecule is fitted onto the base molecule (the first entry in the list).
• Automatically on transformation change: when the user changes the transformation radio button for a non-base molecule, RbTransformation_CheckedChanged calls AutoFitSelectedMolecule(), which re-fits the molecule onto the base using the (possibly mirrored) reference points.
• Manually: the user can re-trigger a fit at any time by re-selecting the same transformation or by removing and re-adding the molecule. There is no dedicated "Fit" button in the current UI.

11.6. Base Molecule and Locked Operations

The first entry in the overlay list is special: it is the base molecule (marked with ⭐ BASE in the list). All subsequent molecules are fitted onto the base. To preserve a stable reference, the following operations are disabled for the base molecule when the overlay contains more than one entry:

• Remove — the base cannot be removed (it would invalidate all fits).
• Center — the base is already at the origin by construction.
• Transformation radio buttons — mirror/inversion on the base would invalidate the reference frame.

The only operation that remains enabled for the base is Color: changing the color of the base is purely cosmetic and does not affect any fit.

When the overlay contains exactly one entry, all operations are enabled for that single molecule (it is not yet "base" in any meaningful sense). When the overlay is empty, all operations are disabled.

11.7. Transformations (Mirror and Inversion)

Before fitting, the user can apply one of three transformations to a non-base molecule. These are exposed as radio buttons in the bottom panel:

Transformation is applied before the automatic fit. The workflow is:

1. The user selects a transformation radio button.
2. RbTransformation_CheckedChanged verifies the molecule is not the base (otherwise it restores the previous radio state and shows a "Cannot transform BASE molecule" status).
3. SetTransformationForMolecule(id, mode) is called. It resets the working atoms from the originals, preserves the current color and style, and applies the appropriate transform matrix.
4. AutoFitSelectedMolecule() is called, which re-fits the (now transformed) molecule onto the base using the reference points.
5. Both the list and the Overlay View are refreshed.

11.8. Colors and Styles

Each overlay entry has its own color, chosen from a random hue at creation time (via GetRandomColor, which picks a random hue and converts to RGB with a sine-based palette). The color is applied to every atom's CustomColor field, so the whole molecule is rendered in a single, distinct color.

To change the color of the selected molecule, click the colored button in the top panel (its background is the current molecule color). A standard ColorDialog opens; on OK, the color is applied via SetMoleculeColor and the list is refreshed. Color can be changed even for the base molecule — it is the only operation that remains enabled for the base.

The style of an overlay entry can be changed programmatically via SetMoleculeStyle(id, styleIndex), but this is not currently exposed as a UI control in the Overlay Manager window. The available styles are 0 = CPK, 1 = Sticks, 2 = Balls and Sticks (the default), 3 = Ellipsoids.

11.9. Overlay View Window

The Show button in the Overlay Manager opens (or brings to front) the OverlayViewWindow, a DockContent hosting a MonoGame viewport that renders all visible overlay atoms. The window is subscribed to OverlayManager.MoleculesChanged, so any change to the overlay (add, remove, transform, color change, visibility toggle) is reflected immediately.

Rendering

OverlayViewWindow.UpdateDisplay() collects the WorkingAtoms of every visible entry and passes them to the underlying GameWindow via SetOverlayAtoms(). The atoms are rendered with the same pipeline as regular molecules, so all display styles (CPK, balls & sticks, ellipsoids) and per-atom colors work as expected. Bond detection is done per-molecule (the overlay does not merge bonds across molecules), so an intermolecular contact in the overlay is not drawn as a bond — this is intentional, since the overlay is meant for comparing molecular geometries, not crystal packings.

Lifecycle

• The window is created lazily: the first time the user clicks "Show", a new OverlayViewWindow is created with the current OverlayManager and docked as a document.
• If a window already exists (e.g., after a layout save/restore), FixExistingOverlayView() detects it and injects the current manager if the existing one is null. This prevents the classic "Overlay View is empty after reloading the layout" bug.
• On close, the window unsubscribes from MoleculesChanged and clears its atom buffer.

11.10. Typical Workflow — Step by Step

1. Load the base structure. Open a CIF/RES/etc. file. The structure appears in the main 3D view.
2. Pick 3 reference atoms on the base. Switch to Pick Atoms mode (or use the default rotate mode with Ctrl+click). Click three atoms that define a unique, non-collinear frame — for example, three atoms of a phenyl ring, or a central atom plus two substituents. You should see them highlighted in yellow.
3. Open the Overlay Manager (Advanced → Overlay Manager...). The list is initially empty.
4. Add the base as the first entry. Click Add Molecule. The molecule containing the 3 selected atoms is added as the first entry, marked ⭐ BASE. Its 3 reference points are stored.
5. Load the target structure. Open a second file. The main 3D view switches to the new structure.
6. Pick 3 corresponding atoms on the target. These must be the "same" 3 atoms as on the base — same chemical environment, same relative geometry. For example, the equivalent three atoms of the same ring.
7. Add the target. Click Add Molecule again. The new molecule is added to the list and automatically fitted onto the base using the 3+3 reference points. The fit status is shown in the status label ("Fitted to base").
8. Inspect the overlay. Click Show to open the Overlay View. Both molecules should overlap at the 3 reference atoms; any deviation elsewhere is a real conformational or symmetry difference.
9. If the fit is wrong, try a mirror. If the two molecules are related by reflection (e.g., enantiomers), select the target entry in the list and click the Mirror radio button. The molecule is re-fitted with the mirror applied first, and should now overlap correctly.
10. Try a different set of reference atoms if needed. Remove the target entry (❌ Remove), pick a different set of 3 atoms, and re-add. The fitting is deterministic — the same 3+3 pair always produces the same result.
11. Color and style. Use the color swatch button to distinguish the molecules. Use Show All to toggle visibility of the whole overlay.
12. Clear when done. Click 🗑️ Clear to remove all entries, or close the Overlay Manager — the overlay state is kept in memory as long as the window is open.

11.11. Limitations and Notes

• Rigid-body only. The fit is a pure rotation + translation. It does not allow bond stretching, angle bending, or torsional flexibility. If the two molecules differ in conformation, the overlay will show a real, unfittable discrepancy at the flexible parts.
• Exactly 3 reference points. The algorithm requires exactly 3. Using more points would require a least-squares (Kabsch) fit, which is not implemented. If you need to fit on more atoms, pick the 3 that best define the rigid core.
• Handedness matters. If the two structures are related by reflection, the direct fit will fail (the fit will be "correct" in the sense of a rigid transform, but the molecule will appear inside-out). Use the Mirror transformation to pre-reflect the source.
• Symmetry-generated atoms. "Add Molecule" uses the atoms as they are in the document, which may include symmetry-generated copies (depending on the current packing mode). If you want just the asymmetric unit, set the packing mode accordingly before adding.
• No persistence. The overlay list is not saved to disk. Closing the Overlay Manager and reopening it starts from an empty list.

12. Advanced Tools

12.1. H-Bond Settings

Opened from the menu (or the corresponding toolbar / context command), the H-Bond Settings window controls how hydrogen bonds are detected and drawn in the 3D view.

Contact distance range

Two modes are available:

• Actual distance — absolute limits in Ångströms (e.g. Minimum = 1.50 Å, Maximum = 3.50 Å).
• VdW distance (default) — limits relative to the sum of van der Waals radii of the donor and acceptor:
  Min/Max = rvdW(donor) + rvdW(acceptor) ± value
  The spin boxes set the offset (e.g. “radii minus 0.50” tightens the contact).

D–H···A angle

Minimum angle (degrees) — geometry filter. A candidate is accepted only if the angle ∠D–H···A is at least this value (default 63°). Typical strong H-bonds are closer to 180°; lower thresholds also accept weaker / bent interactions.

Options

• Show H-Bonds in 3D view — draw detected bonds in the scene.
• Require hydrogen atom to be present — only D–H···A with an explicit H atom; otherwise D···A contacts without H are ignored.
• Draw H-bond from hydrogen — the graphical segment starts at H (not at the donor heavy atom).
• Intermolecular / Intramolecular — include contacts between different molecules and/or within the same molecule.

Donors / Acceptors / Color

Tick which element types may act as donors and acceptors (e.g. N, O, F, P). Only checked pairs are considered. H-Bond color sets the drawing colour in the 3D view (default green).

Use Default to restore factory settings, Apply to recalculate without closing the window, or OK to apply and close.

H-bond list (example table)

After detection, the Structure Info / H-bond panel lists each bond in a table. Column meanings:

Example (illustrative):

Row 1 is an intermolecular O–H···N bond to a symmetry-related molecule in a neighbouring cell. Row 2 is an intramolecular O–H···N contact within the same asymmetric unit (both partners at x,y,z, cell (0,0,0)).

12.2. Atom Display Settings

Refer to section 4.6 for a detailed description of this window.

12.3. Debug Window

Opened from the Window menu, this window shows a detailed log of symmetry operations found in the loaded structure. It's a diagnostic tool for developers and advanced users to troubleshoot issues with symmetry or file parsing.

12.4. BufferedPanel (helper control)

Tal.Forms.BufferedPanel is a System.Windows.Forms.Panel subclass that enables DoubleBuffered plus AllPaintingInWmPaint, OptimizedDoubleBuffer and UserPaint control styles. It is used wherever smooth custom drawing is required (e.g., the powder diffraction plot).

12.5. Structure Info Window

Opened from Tools → Structure Info. The StructureInfoWindow is a dockable panel (a DockContent) that acts as the central, always-available reference for every numeric and textual property of the currently loaded structure. It replaces the need to dig through the raw file or the debug log to find a specific value: every quantity the parsers extracted (cell parameters, symmetry operations, refinement statistics, atomic coordinates with uncertainties, contacts, H-bonds, measured distances, etc.) is laid out in a clean, scrollable table.

12.5.1. Window Layout

The window is divided into four regions:

• Top panel (panelTop): a read-only txtCurrentStructure showing the identifier of the structure currently loaded into the window (e.g., the CIF data block name or the file name), plus a btnCustomise button reserved for future per-window customization.
• Left panel (panelButtons): a vertical stack of section-switching buttons. Exactly one is highlighted (back-coloured with a light blue) at any time to indicate the active section.
• Center panel (panelTableScroll): an auto-scrolling container hosting a single TableLayoutPanel (table). The panel is a BufferedPanel, so the table redraws without flicker when the user switches sections.
• Bottom panel (panelBottom): a Close button.

12.5.2. Section Buttons and Tables

Each button populates the table with a dedicated Section object (a small inner class that holds a list of column headers and a list of string rows). The current sections are:

12.5.3. Error Estimation

For every geometric quantity, the window reports an estimated standard uncertainty (e.s.d.). The formulas are consistent with the ones used by the Measure tool, and propagate the per-atom positional errors stored in Coordinate_error (which already include the cell-parameter contribution):

• Bond length error: if the two atoms have positions p1, p2 with errors e1, e2 and d = |p1 − p2|, then

  σ_d² = (dx/d)²·(e1x² + e2x²)
       + (dy/d)²·(e1y² + e2y²)
       + (dz/d)²·(e1z² + e2z²)

  where (dx, dy, dz) = p1 − p2.
• Angle error: an upper-bound estimate

  σ_angle ≈ √2 · mean(|e1|, |e2|, |e3|) / min(l1, l2)   [rad]
         → σ_angle · 180/π                                  [deg]

  clamped to [0.001°, 10°] to avoid absurd values for near-degenerate geometries.
• Torsion error: an upper-bound estimate

  σ_torsion ≈ √2 · mean(|e1|, |e2|, |e3|, |e4|) / l_BC · 180/π

  where l_BC is the length of the central bond, clamped to [0.001°, 30°]. This reflects the fact that torsions are most sensitive to errors when the central bond is short.

When an error is too small to be meaningful (or the source file did not provide positional errors), the cell shows an em-dash (—) instead of a number.

12.5.4. Torsion Angle Formula

The dihedral angle A-B-C-D is computed by:

a = A − B
b = C − B
c = D − C
d = B − C

n1 = normalize(cross(a, b))
n2 = normalize(cross(c, d))

angle = acos(clamp(dot(n1, n2), -1, 1))

if dot(cross(n1, n2), normalize(b)) < 0:
    angle = −angle

The sign convention matches IUPAC: a positive torsion is a clockwise rotation of the far bond when looking along the central bond from B to C. The result is in radians internally, and converted to degrees for display.

12.5.5. Dynamic Behaviour

• The window caches Section objects in a Dictionary<string, Section>. When the structure changes (a new file is loaded or the active document switches), LoadStructure re-populates every section from the new Molecule, then re-fills the currently visible table.
• The Angles and Torsions sections are special: they are populated lazily, only when the corresponding button is clicked. This avoids the cost of enumerating up to 50,000 rows when the user is only interested in, say, the unit cell parameters.
• Switching between sections does not re-run the parsers — all data is extracted once in PopulateSections and simply re-displayed.

12.5.6. Number Formatting

The window uses the classic crystallographic "value(esd)" notation. For example, a cell parameter a = 5.6789 with a_e = 0.0012 is rendered as 5.6789(12). The number of decimal places is chosen from the magnitude of the e.s.d.: if error < 1, the number of decimals is the position of the first significant digit of the error. Volume is computed from the cell parameters and its e.s.d. is propagated analytically from the partial derivatives with respect to a, b, c, α, β, γ.

12.5.7. Source of the Populated Data (CIF / RES)

All quantities displayed by the Structure Info Window are extracted by the parsers from the source file. The mapping is fully deterministic — a value appears in the table only if the corresponding field is present in the CIF (or RES/INS) file. Missing fields are silently omitted; the window never invents data.

Cells, symmetry and atoms

• Unit cell a, b, c, α, β, γ and their e.s.d.'s — from the CIF _cell_length_* and _cell_angle_* tags, or from the RES CELL/ZERR cards.
• Cell volume — from _cell_volume if present; otherwise computed from a, b, c, α, β, γ.
• Space group number, H-M name — from _space_group_IT_number, _space_group_name_H-M_alt / _symmetry_space_group_name_H-M; RES reads TITL and LATT.
• Symmetry operations — from the _space_group_symop_operation_xyz loop; RES reads the SYMM cards. The full list is displayed verbatim in the Symmetry section.
• Atoms — from the _atom_site_* loop (label, type symbol, fractional x/y/z, x/y/z esd's, occupancy, U_iso, ADP type, disorder group, and the six U^ij components when _atom_site_adp_type = Uani); RES reads the atom cards with the SFAC mapping.

Refinement statistics

• _refine_ls_R_factor_gt → R1 observed.
• _refine_ls_R_factor_all → R1 all.
• _refine_ls_wR_factor_ref / _refine_ls_wR_factor_all → wR2 (all / observed).
• _refine_ls_goodness_of_fit_ref → GooF.
• _refine_ls_number_reflns, _refine_ls_number_parameters, _refine_ls_number_restraints → reflection/parameter/restraint counts.
• _refine_ls_shift/su_max, _refine_ls_shift/su_mean → shift/su max and mean.
• _refine_diff_density_max, _refine_diff_density_min, _refine_diff_density_rms → residual electron density.
• _refine_ls_weighting_scheme and _refine_ls_weighting_details → weighting scheme.
• For RES files the same values are read from the REM lines (R1, wR2, GooF, Nref, Npar, Nres, shift/su max/mean, Δρ max/min).

Reflection statistics

• _diffrn_reflns_number → total reflections.
• _diffrn_reflns_av_R_equivalents and the _reflns_number_gt / _reflns_number_total tags → observed reflections and threshold.
• _diffrn_reflns_theta_min, _diffrn_reflns_theta_max → θ_min, θ_max.
• _diffrn_reflns_resolution_high (or _reflns_d_resolution_high), _diffrn_reflns_resolution_low → resolution limits d_min, d_max.
• _reflns_threshold_expression → the threshold column.

Experiment

• _diffrn_radiation_wavelength, _diffrn_radiation_type → wavelength and radiation type.
• _diffrn_measurement_device_type → measurement device.
• _diffrn_ambient_temperature → temperature (K).
• _exptl_crystal_colour, _exptl_crystal_description → crystal colour and shape.
• _exptl_crystal_size_max, _exptl_crystal_size_mid, _exptl_crystal_size_min → crystal size.
• _exptl_absorpt_coefficient_mu → absorption coefficient μ.
• _exptl_absorpt_correction_type, _exptl_absorpt_correction_T_min, _exptl_absorpt_correction_T_max → absorption correction type and transmission range.
• _exptl_crystal_density_diffrn → calculated density.
• _chemical_formula_sum, _chemical_formula_weight, _chemical_name_common / _chemical_name_systematic → formula, formula weight and compound name.
• _cell_formula_units_Z, _cell_formula_units_Z' → Z and Z'.
• _exptl_crystal_F_000 or _cell_formula_units_Z × _chemical_formula_weight → F(000) if present.

Audit and publication

• _audit_creation_date, _audit_creation_method → audit lines.
• _publ_section_title → publication title.
• _journal_name_full / _journal_name_abbrev → journal name.
• _journal_volume, _journal_page_first, _journal_year → journal reference.
• _journal_doi or _publ_doi → DOI.
• _chemical_melting_point, _chemical_compound_source etc. — read if present; otherwise omitted.

13. Export & Format Integration

13.1. StructureWriter — dispatcher

StructureWriter.Save(structure, path) inspects the file extension and dispatches to the appropriate writer. Returns true on success. Exceptions are caught and logged.

13.2. CifWriter

Writes a StructureData to a CIF file using the IUCr data-name rules. It emits: data_ block, audit lines, chemical metadata, cell parameters (with esd where available), _space_group_IT_number, H-M name, symmetry-operation loop (never guessed — via SymmetryHelper.GetOperationStrings), atom-site loop (label, type, fract coords, occupancy, U_iso, adp_type, optional disorder group), and an anisotropic displacement loop when any atom has AdpType == "Uani".

Number formatting uses CultureInfo.InvariantCulture. The file is written as UTF-8 without BOM.

13.3. ResWriter

Writes a SHELX-format .res file: TITL, CELL, ZERR, LATT, SYMM cards (identity excluded; inversion sign of LATT inferred from SymmetryHelper.ContainsInversion), SFAC/UNIT, atom lines, HKLF 4, END. Occupancy is written as occ + 10 (SHELX free-variable convention). Anisotropic atoms emit the full 6-element U line with the second-line convention.

13.4. XyzWriter

Writes atom count, comment line (structure name), and one Element x y z line per atom. Coordinates are written as stored (Cartesian for non-crystal formats). UTF-8 encoding.

13.5. PdbWriter

Writes HEADER, TITLE, an optional CRYST1 record (when a cell is present, with the H-M name truncated to 11 chars and the space-group number), one HETATM record per atom (serial, 4-char name, "UNK" residue, residue 1, x/y/z, occupancy, B-factor derived from U_iso as B = 8π²·U, element), and END.

13.6. MolWriter

Writes an MDL molfile V2000: header (title truncated to 80 chars, " TAL 3D", blank), counts line (atoms, bonds = 0), atom lines in fixed format (%10.4f%10.4f%10.4f %-3s + zeroed flags), and M END.

13.7. Mol2Writer

Writes a Tripos MOL2 file with @<TRIPOS>MOLECULE, @<TRIPOS>ATOM, and optionally @<TRIPOS>CRYSIN when a cell is available (a, b, c, α, β, γ, space-group number, setting 1). Charges are not written (charge type NO_CHARGES).

13.8. GaussianWriter

Writes a Gaussian input file with a %chk= link, a #p B3LYP/6-31G(d) Opt route, the structure title, a default charge/multiplicity line 0 1, and the Cartesian atom block. ASCII encoding.

13.9. CrystalWriter (.d12)

Writes a CRYSTAL17 geometry + SCF input deck. The writer performs the non-standard-setting check via SpaceGroupTransformer.GetTransformToStandard() and, on user confirmation, transforms a deep copy of the structure (the original in StructureStorage is untouched). The output includes the title line, CRYSTAL, the geometry flag 0 0 0, the space-group number, the cell line, the number of atoms, and one atom line per atom (Z, x, y, z, optional label). The settings block emits OPTGEOM / FULLOPTG / ENDOPT, BASISSET POB-TZVP, DFT B3LYP-D3 XLGRID, SCFDIR, memory sizes, SHRINK 8 8, LEVSHIFT 6 1.0, and TOLINTEG 7 7 7 7 25.

CrystalSpaceGroupValidator.Validate(structure) returns a warning string when the structure is not in the standard setting; otherwise null. The comparison is always done on the normalized H-M name.

13.10. PovRayGenerator

PovRayGenerator(Molecule) produces a complete POV-Ray 3.7 scene file.

Scene structure

• Header: commented metadata (author, date, molecule, atoms, cell, style), POV-Ray version, global_settings { assumed_gamma 2.2; max_trace_level 10 }.
• Camera: position, look-at, up, fixed right, and an angle converted from XNA's vertical FOV to POV-Ray's horizontal angle via angle = degrees(2·atan(tan(vFov/2)·aspect)). The CAMERA_FOV_STRETCH constant (1.0) is applied after conversion.
• Lights: up to two parallel light sources with direction and diffuse intensity taken from Flags; a global ambient is declared as #declare GlobalAmbient.
• Background: a screen-aligned quad placed far behind the scene (quadDist = 500, halfSize = 400) with a vertical gradient between color2 (bottom) and color1 (top). A sky_sphere fallback is used if no camera is available.
• Atoms: style-dependent. CPK spheres scale with the vdW radius / RuntimeIcosphere.BASE_RADIUS. Balls & sticks draw two cylinders per bond (each half colored by its atom). Ellipsoids are drawn as a unit sphere transformed by Transformation · Un · Scale(probability) · Translation(position), with the model scaled by rModel.
• Ellipsoid disks: three thin ellipsoids aligned with the principal axes of the ADP tensor. Disk thickness is constant in Å; the radius scales with the ellipsoid.
• Hydrogen: if ShowHydrogenIsoSpheres is false, H atoms are drawn as small isotropic spheres; otherwise as full ellipsoids.
• Hydrogen bonds: solid or dashed cylinders, colored from Flags.HBondColor; dash/gap lengths are read from flags.
• Unit cell: cylinders along the cell edges, with the edge direction reconstructed from AngleZ and AngleXY and the color saturated by _saturationBoost.
• Polyhedrals: when Presentation == 4, each Polihedrals.VertexesList entry is written as a mesh2 with per-vertex normals, transparent faces, and an inside_vector.
• Crystal planes: visible planes from Molecule.Planes are written as double-sided mesh2 quads (two triangles) with thin cylinder edges. Corners are shifted by −MoleculeShift to align with the atom coordinates.
• Cavities: each rendered cavity with a valid marching-cubes mesh is written as a mesh2. If the draw mode is CavityColorMap and per-vertex colors are available, the mesh uses texture_list and per-face texture indices; otherwise a single transparent pigment is used.
• World transform: the XNA world matrix Scale(cam.OneOverZoom) · Rotate · Translate is emitted once at the end of the union.

Safety guards

• IsValidSegment rejects any cylinder shorter than MIN_SEGMENT_LENGTH = 1e-4 to avoid POV-Ray "degenerate cylinder" errors.
• IsFiniteVector rejects NaN/Inf coordinates.
• Sphere radii smaller than 1e-4 are clamped to 0.01.

Contrast / saturation constants

13.11. SpaceGroupTransformer

Transforms a StructureData from one setting of a space group to another equivalent setting (e.g. P2₁/c ↔ P2₁/n).

Public API

• TransformTo(structure, targetSettingName) — main entry point. Returns a TransformResult with success flag, applied transform, target setting name, and message.
• ApplyTransformToStructure(structure, transform, targetName, wrapAtoms = false) — applies the transform in place: fractional coordinates, ADP tensors (U' = R·U·Rᵀ), cell parameters (M' = M·P), symmetry operations (regenerated for the target setting) and the space-group name.
• ApplyTransform(structure, transform) — low-level version; atoms + ADP + cell only, no symmetry/name update.
• ResolveTransform(structure, targetSettingName) — finds the transform matrix. It first tries a hard-coded lookup and verifies it against the operation sets; if that fails it falls back to a brute-force search over 3×3 integer matrices with determinant ±1 and 27 candidate shifts.
• GetStandardName(group) — returns the standard H-M name for a space group number (first entry with orientation 0).
• GetAvailableSettings(group) — returns all distinct setting names for a group.
• GetCurrentSetting(structure) — returns the setting name currently stored in the structure (matched by normalized name against the internal table).
• GetOperationStringsForSetting(group, settingName) — returns the list of operation strings for a given setting.

Algorithm

1. Resolve the source and target operation sets (ReadSymm) from the internal SymmetryElements table.
2. If the sets are already equivalent, return null (nothing to do).
3. Try a hard-coded transformation; verify it by conjugating the source set and comparing with the target set (3 test points, tolerance 1e-3).
4. If verification fails, brute-force over all integer matrices with determinant ±1 (cached) and a small set of candidate shifts (0, ¼, ½). The first matrix that maps the source set onto the target set is used.
5. Apply the transform to the atoms (optionally wrapping to the unit cell), rotate the ADP tensors, update the cell by M' = M · P where P = Q⁻¹, regenerate the symmetry operations, and set the new name.

Refinement after transform

RefineMoleculePositionsAfterTransform(structure, bonds, maxIterations) performs an iterative relaxation of atoms around their covalent bonds (minimum-image convention), followed by an "attractor" step that pulls smaller connected fragments toward the largest one using integer fractional shifts of up to ±2. Finally, each molecule is translated so that its heaviest non-hydrogen anchor sits near (0.5, 0.5, 0.5) — a whole-molecule translation, not a wrap.

Helper methods: DetectBondsInStructure (minimum-image covalent-bond detection with element-pair thresholds), DetectMoleculesSimple (BFS over bonds), and VerifyRoundTrip (round-trip diagnostic for A→B→A).

13.12. Parsers

The application exposes a set of parsers, each of which populates a data container and then converts to a StructureData:

All parsers call a Clean() method on their container before populating it, and set a completed flag on success.

13.13. CifDivider

CifDivider.Divide(content, filePath) splits a multi-block CIF into a list of CifBlock objects (Name, RawData). It preserves empty lines to keep line numbering intact, scans for lines whose trimmed form starts with data_, and falls back to a single block named after the file when no data_ directive is found.

13.14. RecentFilesManager

Singleton that maintains the recent-files list persisted in %LocalAppData%\Tal\RecentFiles.json. See section 2.2.

13.15. Complex (Tal.Utils)

A lightweight struct Complex (double-precision) used by the powder-diffraction structure-factor calculation. It provides the standard arithmetic operators, Magnitude, MagnitudeSquared, Phase, Conjugate, and a FromPolar factory. The BCL System.Numerics.Complex is intentionally not used to avoid name clashes and to keep the operations minimal.

14. Appendix: Technical Details

14.1. Rendering Pipeline

Tal uses a modern rendering pipeline with the following key stages:

• Geometry Pass: The scene is rendered to a high-resolution off-screen render target (RT). This is where Supersampling Anti-Aliasing (SSAA) takes place. The size of this RT is determined by the Scene Quality settings (e.g., 2x means the RT is twice the width and height of the window).
• Progressive AA Pass: If enabled, the previous frame is blended with the current frame when the camera is stationary. This accumulation of frames acts as a temporal anti-aliasing filter, producing a very clean and smooth final image.
• Sprite Layer Pass: All 2D elements, such as atom labels, measurement annotations, and the color scale legend, are drawn directly onto the high-resolution RT. This ensures that they are perfectly integrated with the 3D scene and benefit from the same anti-aliasing.
• Final Composite Pass: The high-resolution RT is then downsampled and drawn to the screen. This is where the SSAA algorithm performs its final filtering, averaging the sub-pixels to create the final, smooth image.

14.2. Supersampling (SSAA)

Supersampling is the primary method of anti-aliasing in Tal. By rendering the scene at a higher resolution (e.g., 2x, 4x) and then scaling it down, the jagged edges of polygons are smoothed. The SSAASamples flag determines the multiplier. A value of 2 means the scene is rendered at 2x resolution, a value of 4 at 4x, etc. This technique provides excellent image quality, especially for thin features like bonds and ellipsoid disks.

14.3. Progressive Accumulation AA

This technique builds upon SSAA. When the user is not interacting with the scene (i.e., the camera is not moving), the application renders the scene multiple times with a sub-pixel jitter on the camera's projection matrix. The results of these jittered frames are accumulated (averaged) into a single buffer. This effectively multiplies the number of samples per pixel, producing an extremely smooth, high-quality final image. The number of accumulated frames is controlled by the ProgressiveAASamples flag (e.g., 16 frames). This method is perfect for producing publication-quality screenshots.

14.4. Van der Waals Radii and Color Tables

The VdwParameters static class is the single source of truth for all van der Waals radii and the default color palette used in the application. It contains a comprehensive, hardcoded table of radii for elements 1 through 118, along with a default color for each. The Atom Colors (global)... window allows the user to modify the color palette at runtime, and these overrides are stored in the application's settings file (Flags) and applied on startup.

14.5. Settings and Persistence

All application settings, window layouts, and user preferences are managed by the Flags singleton. Key settings are saved to a JSON file on application exit and loaded on startup. This includes window positions, lighting settings, background colors, atom color overrides, and camera parameters. The Flags.Instance.MarkDirty() method is used to flag that settings have changed and need to be saved.

14.6. Camera and POV-Ray FOV Conversion

XNA's Matrix.CreatePerspectiveFieldOfView takes a vertical field of view. POV-Ray's angle parameter is the horizontal angle. The conversion used by PovRayGenerator is:

halfV = radians(fovVertical) / 2
halfH = atan(tan(halfV) * aspect)
angle = degrees(2 * halfH)

The POV-Ray right vector must remain constant (1.0) — POV-Ray scales the image itself from the render resolution (+W/+H), and setting right = aspect would squeeze the image horizontally. The render resolution must therefore have the same aspect ratio as the TAL viewport.

14.7. Powder Diffraction: Cromer–Mann Form Factors

Atomic form factors in StructureFactorCalculator use the 4-Gaussian + constant approximation from International Tables for Crystallography, Vol. C, Table 6.1.1.4:

f0(s) = Σ aᵢ·exp(−bᵢ·s²) + c      (i = 1..4)
s = sinθ / λ

Coefficients are stored as 9-element arrays { a1, b1, a2, b2, a3, b3, a4, b4, c } keyed by atomic number. Elements without tabulated coefficients fall back to a carbon-like default.

14.8. Hydrogen U_iso Recaclulation (SHELX AFIX)

In SHELX RES/INS files, a negative U_iso value on a hydrogen atom encodes an AFIX multiplier: the hydrogen's U_iso is to be derived from its bonded neighbor's equivalent isotropic displacement as

U_iso(H) = |AFIX| · U_eq(neighbor)

ResParser performs a second pass over the parsed atoms. For each hydrogen with afix_multiplier > 0, it finds the nearest non-hydrogen atom within 1.3 Å (minimum-image convention, using the cell's fractional-to-Cartesian matrix), computes U_eq = (U11 + U22 + U33) / 3 for anisotropic neighbors or uses U_iso for isotropic ones, and writes the resulting value into iso. Hydrogens without an AFIX multiplier keep the value given directly by SHELX.

15. Advanced Cavity Analysis (VdW Scaling)

The VdW Scaling analyzer performs a systematic sweep of cavity detection at multiple van der Waals radii. For each scale it collects geometric, electrostatic and channel information, and then derives a set of QSPR descriptors that summarize the structural stability of the framework, its porosity and its potential for host-guest chemistry. This is the tool to reach for when you want to move beyond a single "static" cavity picture and characterize a porous material.

15.1. How to run the analysis

1. Open a crystal structure (CIF, RES, etc.) so that a valid unit cell with symmetry is available.
2. Run at least one normal cavity detection first (Calculate → Cavity Detection) so that the scanner has generated its atoms, its spatial index and at least one density grid. The VdW analysis re-uses the scanner and re-runs the SDF computation for every scale, but it needs a fully initialised scanner and a valid cell.
3. Invoke the analysis programmatically via VdWScaleAnalyzer.RunAnalysisAsync(molecule, detector, scales) (wired from MainForm.RunVdWAnalysisAsync). In the current UI the entry is not a separate top-level Calculate menu item; single-run cavity reports still refer to it as “Calculate → VdW Scaling Analysis…”. Prefer launching it from the cavity workflow after a successful detection, or from a debug/status path that calls RunVdWAnalysisAsync.
4. Optionally provide a custom scale list. If none is given, the default list is used:

   2.00, 1.95, 1.90, 1.85, 1.80, 1.75, 1.70, 1.65, 1.60, 1.55,
   1.50, 1.45, 1.40, 1.35, 1.30, 1.25, 1.20, 1.15, 1.10, 1.05, 1.00

   Larger scales first (aggressive shrinking of the pores), then down to 1.00. This ordering is important for the "collapse scale" descriptor — the analyzer walks from the most aggressive scale down to the physical one, and records the first scale at which total cavity volume drops below the collapse threshold.
5. Watch the progress events (ProgressChanged, StatusChanged) or the status bar. Each scale reports its cavity count and (when available) the largest fitted sphere.
6. When the loop finishes, the analyzer restores the original VdW multiplier and SDF flag state, and stores the VdWAnalysisResults object in VdWScaleAnalyzer.CurrentResults. The report can be saved as TXT or HTML.

r_eff(i, mult) = r_vdw(i) · mult + ProbeRadius

15.2. VdWAnalysisResults — data model

The result object is the single source of truth for the report. Its main fields are:

15.3. Per-cavity snapshot (VdWAnalysisResult)

For every (cavity, scale) pair the analyzer stores both geometric and electrostatic data in a single VdWAnalysisResult:

• Geometric: Volume, VolumeError, SurfaceArea, SurfaceAreaError, Sphericity, VolumeToSurfaceRatio, RatioError, TriangleCount, VoxelCount, WallConnections.
• Electrostatic: MeanCharge, ChargeVariance, MeanHydrophobicity, HydrophobicityVariance, MeanZEffective, DipoleMomentX/Y/Z, DipoleMomentMagnitude. The dipole is computed as the charge-weighted first moment of the cavity voxels.
• Channel information: HasChannelsAtThisScale, ChannelDirections. This reflects the channel state at VdW = 1.0 (channels are only analyzed at that scale to keep the runtime manageable).

The most important derived quantity is VolumeToSurfaceRatio (V/S), which behaves like a generalized radius: it grows with the cavity size and shrinks when a cavity becomes rough or fragmented.

15.4. Cavity tracking and critical events

After the sweep, the analyzer builds a CavityTrackData for every cavity ID it saw at least once. The tracking object contains one CavitySnapshot per scale, with the volume/surface/electrostatic data frozen at that point. The DisappearanceScale is the first scale at which the volume drops below EPSILON after having been positive earlier.

CriticalEvent records all interesting transitions:

• Disappearance: a cavity that existed at scale s_i is absent at scale s_i+1.
• ChannelOpen / ChannelClose: the channel state (present vs. not present) changed for a given cavity between two consecutive scales.
• Merge / Split / TopologyChange: reserved for future use.

15.5. QSPR descriptors

Six descriptors summarize the whole sweep:

15.6. Interpreting the shape of the V(mult) curve

The VdW sweep is not just a way to compute f_collapse. The shape of the total volume curve V(mult) is itself diagnostic and allows you to distinguish qualitatively different kinds of cavities:

The volume of a spherical cavity of radius R scales as V ∝ R³. The VdW multiplier shrinks R by δ = (mult − 1) · r_vdw, so the expected profile for a compact cavity is a cubic polynomial in mult. Deviations from that cubic — in particular a sharp drop followed by a long tail — indicate that the cavity is not a compact object but a network of narrow necks. This is exactly the kind of qualitative distinction that the VdW sweep is designed to make visible.

15.7. Sphere packing analysis (SphereFitAnalyzer)

For each cavity detected at VdW = 1.0 the analyzer builds a SpherePackingResult. The sphere fits themselves are found by SphereFitAnalyzer, which operates entirely on the Cavity object and does not need a graphics device.

Small cavities

If the cavity has fewer than MIN_VOXELS_SMALL = 50 voxels, the analyzer runs a FindSingleMaximum routine: a coarse grid scan (step = min(resX,resY,resZ)/20) locates the maximum of the SDF, then a local refinement over a radius of 2·step voxels finds the true apex. A single SphereFitResult is created with MaxSphereRadius = SDF_max.

Large cavities

For larger cavities the analyzer enumerates all local maxima of the SDF in parallel and then selects a non-overlapping subset:

1. FindLocalMaximaParallel: scan with step = min(resX,resY,resZ)/15, test the 6 face neighbours, and refine each candidate over a 3-voxel cube.
2. Filter: discard every local maximum whose SDF is below MIN_SDF_VALUE_LARGE = 0.5 Å.
3. SelectSpheres: sort by descending SDF; add a candidate to the selected set only if its distance to every already-selected sphere is at least max( (rcand + rsel) · SEPARATION_FACTOR / 10, MIN_DIST_FRAC = 0.05 ) in fractional coordinates. This yields a maximal packing of spheres.
4. Result: one SphereFitResult per selected sphere, with the cavity ID, sphere radius, sphere center, fractional coordinates, wall flag, and distance to the cavity center.

SpherePackingResult

15.8. Gas compatibility

For every cavity the analyzer cross-checks the fitted spheres against a fixed library of gases. Each gas is a GuestMolecule with a kinetic diameter, polarizability, dipole moment, and a flag IsGas = true. The library includes:

• H₂ (2.9 Å), He (2.6 Å), CH₄ (3.8 Å), CO₂ (4.2 Å), N₂ (3.6 Å), O₂ (3.7 Å), Ar (3.9 Å), NH₃ (4.0 Å).

For each gas the result is a GasGuestResult:

• CanFit: gasRadius ≤ 0.95 · maxSphereRadius.
• VolumeRatio: gasVolume / cavityVolume.
• Score: 0.70 · sizeMatch + 0.30 · polarityMatch, where sizeMatch = 1 − |gasRadius / maxRadius − 0.5|·2.
• Status: GOOD (score ≥ 0.7), FAIR (score ≥ 0.4), POOR.
• BestSphereLocation: textual reference to the closest-sized sphere (index + fractional coordinates).
• RecommendedCount: floor(usableVolume / gasVolume).
• IsDiffusionLimited: true if the cavity has wall connections but no channel (the gas can enter, but diffusion is restricted).

15.9. Solvent compatibility

The same library (with IsGas = false) is used to score solvents. For each solvent the analyzer builds a SolventGuestResult whose score combines four indicators:

• sizeMatch — how well the solvent radius matches the local sphere radius.
• polarMatch — closeness of the solvent polarizability to a reference value of 0.5.
• chargeMatch — penalizes strongly polar solvents when the cavity is not particularly polar.
• hydroMatch — closeness of the solvent hydrophobicity to zero.

The final score is a weighted sum (0.25 size + 0.15 polar + 0.15 charge + 0.15 hydro + 0.10 shape + 0.10 voxel) and VoxelClusterMatch is a derived 0.5–1.0 confidence that the solvent fits the shape of the cavity's voxel cloud, not just the largest sphere.

15.10. Semi-periodic channels (CavityPathFinder)

Channel analysis is run at VdW = 1.0 for every cavity that touches at least one cell wall (WallConnections > 0). The CavityPathFinder performs a semi-periodic A* search along the three cell axes and stores its results in cavity.CrystalChannelResults.

Path result

Algorithm

1. Setup: the density grid, the periodic grid and the inverse fractional-to-Cartesian matrix are flattened into contiguous arrays. The maximum SDF value inside the cavity volume is recorded for the heuristic.
2. Start points: up to MaxStartPoints = 20 candidates are enumerated in expanding shells around the center of the periodic cell. A candidate is accepted only if the local voxel is not a wall and has at least 2 units of free space in the 6 axis directions.
3. A* search: for each start point, two separate A* searches are run — one toward the minimum wall and one toward the maximum wall. Each search uses a 26-neighbour stencil and a cost function distance + distance · (1 − SDF / SDF_max) · SDFWeight so that the path prefers to stay near the center of the cavity.
4. Path assembly: the two half-paths are concatenated. Paths with fewer than 3 points are discarded.
5. Measurement: total length, helicity (angle between the first and second half), cell crossings, and the diameter at each of up to 20 sample points.
6. Smoothing (optional): a Laplacian smoothing pass (iterations controlled by SmoothingIterations, default 3). The RMS displacement is reported as SmoothingError.
7. Statistics: the analyzer aggregates the results of all start points into a ChannelPathStatistics object (min/max/avg length, diameter, helicity).
8. XNA conversion: WorldPathXna is filled from WorldPath so the renderer can consume it directly.

15.11. Reports

VdWAnalysisReport.GenerateTxtReport(results) produces a plain-text report; GenerateHtmlReport(results) produces a self-contained dark-theme HTML report; SaveReport(results, path) dispatches to the right one by extension. Both reports are composed of the same sections:

• Header (cell parameters, analysis parameters, summary).
• QSPR descriptors.
• Scale summary (cavity count, total volume/surface, V/S, Vol/Voxels, Vol/Atom, channels per scale).
• Scale dependence (dV/dscale, dS/dscale, ΔV%, ΔS%, V/S).
• Sphere packing (per-cavity table + a detailed sphere list for the largest packing).
• Gas compatibility (the top 5 gases).
• Solvent compatibility (the top 8 solvents, with all four match metrics).
• Electrostatic & hydrophobic summary (mean/variance of charge and hydrophobicity, dipole, surface area, per scale).
• Critical events (disappearance, channel open/close).
• Cavity survival tracking (Volume and Surface at 1.0, 1.1, 1.2, 1.3, 1.4, 1.5).
• Channel analysis (channel summary, statistics per direction, network type: 1D / 2D / 3D).

16. Atom Visibility Manager

The AtomVisibilityManager is a small static helper class in Tal.DrawMethods that centralizes the "show / hide" operations for atoms. It is not bound to any particular tool window; any command that needs to hide or show atoms should go through it so that all the ancillary structures (contacts, hydrogen bonds, etc.) stay consistent.

The class is stateless: it does not cache any data, it just manipulates the Molecule.Atom.show flag and the atomic contact arrays. This follows the same pattern as ShortContactsFinder and HBondFinder.

16.1. Hide operations

16.2. Hanging contacts cleanup

Contacts in Tal are stored twice — once on each of the two atoms that participate in the contact. If one end is later hidden or removed, the other end is left with a dangling reference to an atom that is no longer visible. RemoveHangingContacts(molecule) walks every atom, filters its Contact* arrays, and drops every entry whose target atom is out of range or has show == false. The filtered arrays are:

• ContactDistance
• ContactAtom
• ContactAngleZ
• ContactAngleXY
• PrecomputedContactMatrix
• ContactIsNew

ContactCount is updated to match the new list length. The method returns the total number of contacts removed and writes a single line to the debug log.

16.3. Statistics

• CountVisible(molecule) — number of atoms with show == true.
• CountHidden(molecule) — NumberOfAtoms − CountVisible(molecule).

16.4. Typical workflows

Isolate a single molecule in a packing view:

1. Pick one atom of the target molecule (Pick Atoms mode).
2. Call HideReverseSelected(molecule) first (so that only the picked atom's selection is kept), then HideMolecule(molecule, atomIndex) to hide every copy of that molecule except the selected one.
3. Call RemoveHangingContacts(molecule) to clean up dangling contacts.

Hide all copies of a single symmetry-equivalent atom:

1. Pick one atom.
2. HideSymmetryEquivalent(molecule, atomIndex).

Undo all hides:

1. ShowAll(molecule).

17. Toolbars

Toolbars are built in MainForm.Toolbars.cs. They sit under the menu strip and provide style, camera, and picking controls without opening dock windows.

17.1. Style strip

17.2. View strip

17.3. Picking strip

• Clear Measurements — removes distance/angle/torsion annotations.
• Show Labels for: All atoms · Not hydrogen atoms · Selected atoms · None.
• Extra buttons on the strip include 👆 Pick and lasso shortcuts where present.

18. Hydrogen Bonds (H-Bond Settings)

Detection and drawing are implemented by HBondFinder (Tal/Calculate/HBondFinder.cs). The settings UI is HBondSettingsWindow (title H-Bond Settings), opened via OnHBondSettings / Display Options → H-Bond related paths.

18.1. Contact distance range

18.2. D–H⋯A angle

• Minimum angle (degrees) — numAngle. Bonds with a smaller D–H⋯A angle are rejected.

18.3. Options

• Show H-Bonds in 3D view
• Intramolecular / Intermolecular
• Draw H-bond from hydrogen — cylinder starts at H rather than donor heavy atom.
• Require hydrogen atom to be present — skip putative D⋯A contacts without an explicit H.

18.4. Donors / Acceptors / Color

• Editable donor and acceptor element lists (classification may be incomplete for non-CSD-style atom typing — see on-form warning).
• H-Bond color — maps to Flags.HBondColor for the 3D cylinders (solid or dashed; dash geometry from flags).

18.5. Actions

• OK — apply and close.
• Apply — apply without closing.
• Default — restore built-in criteria.

19. Hydrogen Normalization

Menu: Calculate → Normalize Hydrogens... → OnNormalizeHydrogens. Window: HydrogenNormalizationWindow (title Hydrogen Normalization). Engine: HydrogenNormalizer.

19.1. When to use it

• X-ray structures often have shortened X–H distances; neutron-standard lengths improve Hirshfeld surfaces and H-bond geometry.
• The Hirshfeld pipeline can also normalize silently before surface generation (see main help §7.2).

19.2. Window controls

19.3. Length table

Standard distances live in HydrogenNormalizer.NeutronLengths (element-pair keys). Analysis uses HydrogenNormalizer.Check over the asymmetric unit (and related helpers for Hirshfeld).

19.4. Workflow

1. Open a structure with explicit hydrogen atoms.
2. Calculate → Normalize Hydrogens... — molecule is loaded before the window is shown.
3. Review the table; adjust Skip bridging if needed.
4. Apply — updates coordinates in the active molecule.
5. Optional Reset to restore the pre-normalization geometry from the session cache.

20. Display Options, Depth Cue & Auto-Rotate

20.1. Display Options window

Dock window DisplayOptionsWindow, toggled from Window → 🎨 Display Options (OnToggleDisplayOptions).

Display group

• Packing — same flag as Packing window / Flags.DrawPacking.
• Short Contact — show short-contact lines.
• Asymmetric Unit — restrict view to the asymmetric unit when applicable.
• H-Bond — show hydrogen bonds.
• Auto centre — keep the structure centred on interaction.

Options group

• Show hydrogens
• Show cell axes
• Label atoms
• Depth cue — enables depth cueing (see below).
• Z-Clipping
• Stereo

Buttons

• Reset — default display flags.
• Contacts... — opens Short Contacts settings.
• More Info — structure / info path.
• Powder pattern — shortcut to powder diffraction.

20.2. Depth Cue window

DepthCueWindow fine-tunes fog / attenuation when Depth cue is enabled. Opened from display options or related commands; parameters are stored in Flags and applied in the draw pipeline.

20.3. Auto-Rotate window

AutoRotateWindow complements the toolbar ▶ Animate button with explicit speed / axis options when available. Animation advances the camera or model rotation each frame while active.

21. Complete Menu Map

Built by InitializeMenu() in MainForm.Menu.cs.

File

• Open... (Ctrl+O)
• Open Recent
• Save (Ctrl+S) · Save As...
• Save Screenshot... (Ctrl+Shift+S)
• Export to POV-Ray...
• Close (Ctrl+W) · Exit (Alt+F4)

Edit

• Undo (Ctrl+Z) · Redo (Ctrl+Y)
• Cut · Copy · Paste
• Select All (Ctrl+A)

Selection

• Select Atoms · Select Molecules · Select Unit Cell
• Invert Selection · Clear Selection

Display

• Style · Colour atom
• Labels · Show Hydrogens · Show Unit Cell
• Short Contacts...
• Reset View · Background... (Ctrl+B) · Lighting...
• Atom Display Settings · Atom Colors (global)...
• Scene Quality...

Calculate

• Centroids... (stub)
• Planes... · Packing... · Contacts...
• Molecular Shell... (stub) · Graph Sets... (stub)
• Powder Pattern...
• Cavity Settings... (Ctrl+Shift+C)
• Transform space group...
• Normalize Hydrogens...

Advanced

• Hirshfeld Analysis...
• Morphology...
• Molecule Property Mapping...
• Structure Overlay...

View

• Reset Layout

Window

• Structure List · Display Options · Debug — Symmetry

Help

• Contents (F1) · About

22. License

TAL — Crystallographic Application
Copyright © 2026 Daniel M. Kamiński, Lublin, Poland
All rights reserved.

Website: https://dkami.umcs.pl/tal/
Contact: daniel.kaminski@umcs.pl

Academic / Non-Commercial License

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to use, copy, modify, and distribute the Software for non-commercial, academic, and educational purposes only , subject to the conditions below.

1. Non-commercial use only

The Software may be used solely for non-commercial purposes, including but not limited to academic research, teaching, and personal study. Commercial use of any kind is strictly prohibited without a separate commercial license from the copyright holder.

“Commercial use” includes, without limitation:

• use in a commercial product or service,
• use by a for-profit company or organization in the course of its business,
• use that generates revenue (directly or indirectly),
• use in consulting, contract research, or industrial R&D performed for a paying client.

2. Attribution

Any distribution, publication, or derivative work that uses or is based on the Software must include the following notice:

TAL - Crystallographic Application
Copyright (c) 2026 Daniel M. Kamiński
Contact: daniel.kaminski@umcs.pl
Website: https://dkami.umcs.pl/tal/

3. No warranty

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

4. Commercial licensing

For commercial licensing inquiries, please contact:
Daniel M. Kamiński
E-mail: daniel.kaminski@umcs.pl
Lublin, Poland

5. Third-party components

This Software incorporates third-party libraries distributed under their own licenses (including but not limited to .NET, MonoGame, DockPanel Suite, SharpDX). Those licenses remain in force and are unaffected by this license.