Fix DiffractionSim coordination polyhedra for Si, SrTiO3 and elemental metals

CHANGELOG.md

@@ -6,6 +6,7 @@ new `rcN` heading when that rc is published to TestPyPI.
 
 ## Unreleased
 
+- DiffractionSim coordination polyhedra: in a compound, centres bond only to the most numerous species and get a polyhedron for 4 to 8 ligands, so SrTiO3 shows TiO6 octahedra instead of a cube of Sr around each Ti. An elemental crystal gets polyhedra on one set of a two-set neighbour network, which gives Si four SiSi4 tetrahedra inside the cell instead of interpenetrating tetrahedra on all eight atoms, and fcc, bcc and hcp metals get none.
 - Add `DiffractionSim`: rotate a unit cell and see its nanobeam, CBED or Kikuchi pattern update live, with kinematical or Bloch wave intensities computed in the browser, state save/load, and HTML export. It needs quantem 0.1.10 or later (`quantem.diffraction.Crystal` and `bloch`).
 - Add `quantem.diffsim`, a separate PyPI package with no dependencies, built from `packages/quantem.diffsim`: the browser build of DiffractionSim (`npm run build:web`) for web pages and Dash apps, with `settings()` and `save_html()` for pymatgen, ASE or CIF structures and `quantem.diffsim.dash` for Dash. A page can pass its structure as CIF text, which replaces the structure menu, and shows a "Powered by quantEM" link to the quantEM documentation. CIFs loaded in the browser now get the absorptive Weickenmeier-Kohl couplings with a Debye-Waller factor, ported from quantem, instead of an absorption fixed at 8 % of the elastic coupling, which was up to 5 times too strong for the low-order reflections of light elements and up to 4 times too weak for the high-order reflections of gold.
 - Drop abTEM. The ShowCIF potential preview takes its Lobato scattering factors from quantem (the tables agree with abTEM's to 1e-4), so the `crystal` extra is now empty and kept only so existing install lines work.

js/diffractionsim/crystal3d.test.ts

@@ -0,0 +1,43 @@
+import { describe, expect, it } from "vitest";
+
+import { cellGeometry } from "./crystal3d";
+import { parseCrystal } from "./physics";
+import { PRESETS } from "../diffractionsim-web/presets";
+
+function geometry(name: string) {
+  return cellGeometry(parseCrystal(PRESETS[name])!, [1, 1, 1], true);
+}
+
+function edges(verts: number[][]): number[] {
+  return verts.map((v, i) => {
+    const w = verts[(i + 1) % verts.length];
+    return Math.hypot(v[0] - w[0], v[1] - w[1], v[2] - w[2]);
+  });
+}
+
+describe("diffsim coordination polyhedra", () => {
+  it("draws Si tetrahedra on one fcc sublattice, inside the cell", () => {
+    const g = geometry("Si (diamond cubic)");
+    // 4 tetrahedra around the interior sublattice, with vertices on the drawn corner and face atoms
+    expect(g.faces.length).toBe(16);
+    for (const face of g.faces) {
+      for (const e of edges(face.verts)) expect(e).toBeCloseTo(5.431 / Math.SQRT2, 2);
+      for (const v of face.verts) {
+        expect(g.atoms.some((a) => Math.hypot(a.pos[0] - v[0], a.pos[1] - v[1], a.pos[2] - v[2]) < 1e-6)).toBe(true);
+      }
+    }
+  });
+
+  it("draws TiO6 octahedra in SrTiO3, without Sr", () => {
+    const g = geometry("SrTiO3 (perovskite)");
+    expect(g.faces.length).toBe(8);
+    for (const face of g.faces) {
+      expect(face.verts.length).toBe(3);
+      for (const e of edges(face.verts)) expect(e).toBeCloseTo(3.905 / Math.SQRT2, 2);
+    }
+  });
+
+  it("draws no polyhedra for an fcc metal", () => {
+    expect(geometry("Al (fcc)").faces.length).toBe(0);
+  });
+});

js/diffractionsim/crystal3d.ts

@@ -47,9 +47,9 @@ export interface CellGeometry {
 /**
  * Geometry of an na x nb x nc block of cells. Atoms on the block boundary
  * are repeated (the corner atoms of fcc all appear). Polyhedra are drawn
- * around every species except the most numerous one (cations in an oxide;
- * every atom of an elemental crystal), using neighbours within
- * 1.2 (r_i + r_j) of the covalent radii, including atoms outside the block.
+ * around the centres chosen by polyhedronSites that have 4 to 8 bonded
+ * neighbours within 1.2 (r_i + r_j) of the covalent radii, including atoms
+ * outside the block.
  */
 export function cellGeometry(c: CrystalData, nCells: [number, number, number] = [1, 1, 1], polyhedra = false): CellGeometry {
   const cell = c.cell;
@@ -101,31 +101,24 @@ export function cellGeometry(c: CrystalData, nCells: [number, number, number] =
 
   const faces: Face[] = [];
   if (polyhedra) {
-    const counts = new Map<string, number>();
-    for (const s of c.symbols) counts.set(s, (counts.get(s) || 0) + 1);
-    const species = [...counts.keys()];
-    let centres = new Set(species);
-    if (species.length > 1) {
-      const most = species.reduce((a, b) => ((counts.get(a) || 0) >= (counts.get(b) || 0) ? a : b));
-      centres = new Set(species.filter((s) => s !== most));
-    }
+    const { centre, ligand } = polyhedronSites(c, cart);
     const cutoffMax = 1.2 * 2 * Math.max(...c.radii);
     for (let i = 0; i < halo.length; i++) {
       if (!inBlock[i]) continue;
       const si = halo[i].species;
-      if (!centres.has(c.symbols[si])) continue;
+      if (!centre[si]) continue;
       const pi = halo[i].pos;
       const nbr: Vec3[] = [];
       for (let j = 0; j < halo.length; j++) {
         if (j === i) continue;
         const sj = halo[j].species;
-        if (species.length > 1 && c.symbols[sj] === c.symbols[si]) continue;
+        if (!ligand[sj]) continue;
         const pj = halo[j].pos;
         const d = Math.hypot(pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]);
         if (d > cutoffMax) continue;
         if (d <= 1.2 * (c.radii[si] + c.radii[sj])) nbr.push(pj);
       }
-      if (nbr.length >= 4 && nbr.length <= 14) {
+      if (nbr.length >= 4 && nbr.length <= 8) {
         for (const poly of convexHullFaces(nbr, pi)) faces.push({ verts: poly, color: rgb(c.colors[si]) });
       }
     }
@@ -139,6 +132,60 @@ export function cellGeometry(c: CrystalData, nCells: [number, number, number] =
   };
 }
 
+/**
+ * Sites that get a polyhedron (centre) and sites they bond to (ligand). In a
+ * compound the ligands are the most numerous species (the anions of an
+ * oxide) and every other site is a centre, so cations do not bond to each
+ * other (Ti in SrTiO3 gets TiO6, not the Sr cube). In an elemental crystal
+ * polyhedra on every atom interpenetrate, so they are drawn on one set of a
+ * two-set neighbour network (one fcc sublattice of diamond), choosing the set
+ * with fewer sites on the cell boundary so its polyhedra stay inside the
+ * cell. Other elemental crystals (fcc, bcc and hcp metals) get none.
+ */
+function polyhedronSites(c: CrystalData, cart: (f: number[]) => Vec3): { centre: boolean[]; ligand: boolean[] } {
+  const n = c.symbols.length;
+  const counts = new Map<string, number>();
+  for (const s of c.symbols) counts.set(s, (counts.get(s) || 0) + 1);
+  const species = [...counts.keys()];
+  if (species.length > 1) {
+    const most = species.reduce((a, b) => ((counts.get(a) || 0) >= (counts.get(b) || 0) ? a : b));
+    return { centre: c.symbols.map((s) => s !== most), ligand: c.symbols.map((s) => s === most) };
+  }
+  const eps = 1e-4;
+  const frac = c.positions_frac.map((f) => f.map((x) => x - Math.floor(x + eps)));
+  const pos = frac.map(cart);
+  const adj: number[][] = Array.from({ length: n }, () => []);
+  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) {
+    for (let sx = -2; sx <= 2; sx++) for (let sy = -2; sy <= 2; sy++) for (let sz = -2; sz <= 2; sz++) {
+      const p = cart([frac[j][0] + sx, frac[j][1] + sy, frac[j][2] + sz]);
+      const d = Math.hypot(p[0] - pos[i][0], p[1] - pos[i][1], p[2] - pos[i][2]);
+      if (d > 1e-6 && d <= 1.2 * (c.radii[i] + c.radii[j])) adj[i].push(j);
+    }
+  }
+  const none = { centre: new Array<boolean>(n).fill(false), ligand: new Array<boolean>(n).fill(true) };
+  const colour = new Array<number>(n).fill(-1);
+  for (let s = 0; s < n; s++) {
+    if (colour[s] >= 0) continue;
+    colour[s] = 0;
+    const queue = [s];
+    while (queue.length) {
+      const i = queue.shift()!;
+      for (const j of adj[i]) {
+        if (colour[j] < 0) {
+          colour[j] = 1 - colour[i];
+          queue.push(j);
+        } else if (colour[j] === colour[i]) {
+          return none;
+        }
+      }
+    }
+  }
+  const onBoundary = [0, 0];
+  for (let i = 0; i < n; i++) if (frac[i].some((x) => x < eps)) onBoundary[colour[i]]++;
+  const keep = onBoundary[1] < onBoundary[0] ? 1 : 0;
+  return { centre: colour.map((k) => k === keep), ligand: none.ligand };
+}
+
 /** Faces of the convex hull of a few points (brute force, n <= 14), as outward-ordered polygons. */
 function convexHullFaces(pts: Vec3[], inside: Vec3): Vec3[][] {
   const n = pts.length;