feat: compute the first Brillouin zone from the reciprocal lattice

src/js/lattice/reciprocal/brillouin_zone.ts

@@ -0,0 +1,199 @@
+import { Vector3DSchema } from "@mat3ra/esse/dist/js/types";
+
+/**
+ * A face of the first Brillouin zone: the polygon cut by the perpendicular bisector plane
+ * ("Bragg plane") of one reciprocal lattice vector.
+ */
+export interface BrillouinZoneFace {
+    /** Polygon vertices in reciprocal space, ordered counter-clockwise about `normal`. */
+    vertices: Vector3DSchema[];
+    /** Outward unit normal, along the reciprocal lattice vector bounding this face. */
+    normal: Vector3DSchema;
+}
+
+/**
+ * Shells of reciprocal lattice points considered when bounding the cell. Only the nearest
+ * points can contribute a face, so the search is truncated well before the triple loop below
+ * becomes expensive.
+ */
+const MAX_SHELL_INDEX = 2;
+const MAX_CANDIDATE_PLANES = 40;
+
+const SINGULAR_MATRIX_TOLERANCE = 1e-9;
+const HALF_SPACE_TOLERANCE = 1e-7;
+const COINCIDENT_POINT_TOLERANCE = 1e-6;
+const ON_PLANE_TOLERANCE = 1e-6;
+
+function crossProduct(first: Vector3DSchema, second: Vector3DSchema): Vector3DSchema {
+    return [
+        first[1] * second[2] - first[2] * second[1],
+        first[2] * second[0] - first[0] * second[2],
+        first[0] * second[1] - first[1] * second[0],
+    ];
+}
+
+function dotProduct(first: Vector3DSchema, second: Vector3DSchema): number {
+    return first[0] * second[0] + first[1] * second[1] + first[2] * second[2];
+}
+
+function subtract(first: Vector3DSchema, second: Vector3DSchema): Vector3DSchema {
+    return [first[0] - second[0], first[1] - second[1], first[2] - second[2]];
+}
+
+function scale(vector: Vector3DSchema, factor: number): Vector3DSchema {
+    return [vector[0] * factor, vector[1] * factor, vector[2] * factor];
+}
+
+function vectorLength(vector: Vector3DSchema): number {
+    return Math.sqrt(dotProduct(vector, vector));
+}
+
+function normalize(vector: Vector3DSchema): Vector3DSchema {
+    const magnitude = vectorLength(vector);
+    return magnitude === 0 ? [0, 0, 0] : scale(vector, 1 / magnitude);
+}
+
+/** Solves `matrix * x = rightHandSide` by Cramer's rule; null when the matrix is singular. */
+function solveLinearSystem(
+    matrix: [Vector3DSchema, Vector3DSchema, Vector3DSchema],
+    rightHandSide: Vector3DSchema,
+): Vector3DSchema | null {
+    const determinant = dotProduct(matrix[0], crossProduct(matrix[1], matrix[2]));
+    if (Math.abs(determinant) < SINGULAR_MATRIX_TOLERANCE) {
+        return null;
+    }
+    const determinantWithColumnReplaced = (columnIndex: 0 | 1 | 2): number => {
+        const replaced = matrix.map((row, rowIndex) => {
+            const nextRow: Vector3DSchema = [...row];
+            nextRow[columnIndex] = rightHandSide[rowIndex];
+            return nextRow;
+        }) as [Vector3DSchema, Vector3DSchema, Vector3DSchema];
+        return dotProduct(replaced[0], crossProduct(replaced[1], replaced[2]));
+    };
+    return [
+        determinantWithColumnReplaced(0) / determinant,
+        determinantWithColumnReplaced(1) / determinant,
+        determinantWithColumnReplaced(2) / determinant,
+    ];
+}
+
+/**
+ * Computes the first Brillouin zone — the Wigner-Seitz cell of the reciprocal lattice.
+ *
+ * The zone is the set of points closer to the origin than to any other reciprocal lattice
+ * point `G`, i.e. the intersection of the half-spaces `k · G <= |G|^2 / 2`. Its vertices are
+ * the points where three bounding planes meet while satisfying every other half-space, and its
+ * faces group the vertices lying on each plane.
+ *
+ * The shape follows from the lattice itself, not from its Bravais type: two materials of the
+ * same type but different axial ratios (a bulk crystal and a slab with vacuum padding, say)
+ * have differently proportioned zones.
+ *
+ * @param reciprocalVectors - the three reciprocal lattice vectors, e.g.
+ *                            `new ReciprocalLattice(material.lattice).reciprocalVectors`.
+ * @returns the zone's faces, or null when the vectors are degenerate (coplanar or zero).
+ */
+export function computeBrillouinZone(
+    reciprocalVectors: Vector3DSchema[],
+): BrillouinZoneFace[] | null {
+    if (reciprocalVectors.length !== 3) {
+        return null;
+    }
+    const [firstVector, secondVector, thirdVector] = reciprocalVectors;
+    const isDegenerate = reciprocalVectors.some(
+        (vector) => vector.length !== 3 || vector.some((component) => !Number.isFinite(component)),
+    );
+    if (isDegenerate) {
+        return null;
+    }
+
+    const latticePoints: Vector3DSchema[] = [];
+    for (let first = -MAX_SHELL_INDEX; first <= MAX_SHELL_INDEX; first += 1) {
+        for (let second = -MAX_SHELL_INDEX; second <= MAX_SHELL_INDEX; second += 1) {
+            for (let third = -MAX_SHELL_INDEX; third <= MAX_SHELL_INDEX; third += 1) {
+                if (first !== 0 || second !== 0 || third !== 0) {
+                    latticePoints.push([
+                        first * firstVector[0] + second * secondVector[0] + third * thirdVector[0],
+                        first * firstVector[1] + second * secondVector[1] + third * thirdVector[1],
+                        first * firstVector[2] + second * secondVector[2] + third * thirdVector[2],
+                    ]);
+                }
+            }
+        }
+    }
+
+    const planes = latticePoints
+        .sort((left, right) => vectorLength(left) - vectorLength(right))
+        .slice(0, MAX_CANDIDATE_PLANES)
+        .map((latticePoint) => ({
+            normal: latticePoint,
+            offset: dotProduct(latticePoint, latticePoint) / 2,
+        }));
+
+    const isInsideZone = (point: Vector3DSchema) =>
+        planes.every(
+            (plane) => dotProduct(point, plane.normal) <= plane.offset + HALF_SPACE_TOLERANCE,
+        );
+
+    const vertices: Vector3DSchema[] = [];
+    for (let first = 0; first < planes.length; first += 1) {
+        for (let second = first + 1; second < planes.length; second += 1) {
+            for (let third = second + 1; third < planes.length; third += 1) {
+                const point = solveLinearSystem(
+                    [planes[first].normal, planes[second].normal, planes[third].normal],
+                    [planes[first].offset, planes[second].offset, planes[third].offset],
+                );
+                const isZoneVertex = point !== null && isInsideZone(point);
+                const isDuplicate =
+                    isZoneVertex &&
+                    vertices.some(
+                        (existing) =>
+                            vectorLength(subtract(existing, point as Vector3DSchema)) <
+                            COINCIDENT_POINT_TOLERANCE,
+                    );
+                if (isZoneVertex && !isDuplicate) {
+                    vertices.push(point as Vector3DSchema);
+                }
+            }
+        }
+    }
+    if (vertices.length < 4) {
+        return null;
+    }
+
+    const faces: BrillouinZoneFace[] = [];
+    planes.forEach((plane) => {
+        const verticesOnPlane = vertices.filter(
+            (vertex) =>
+                Math.abs(dotProduct(vertex, plane.normal) - plane.offset) < ON_PLANE_TOLERANCE,
+        );
+        if (verticesOnPlane.length < 3) {
+            return;
+        }
+
+        // Order the polygon by angle about the face normal, in a basis lying in the face.
+        const normal = normalize(plane.normal);
+        const centroid = scale(
+            verticesOnPlane.reduce<Vector3DSchema>(
+                (sum, vertex) => [sum[0] + vertex[0], sum[1] + vertex[1], sum[2] + vertex[2]],
+                [0, 0, 0],
+            ),
+            1 / verticesOnPlane.length,
+        );
+        const inPlaneAxis = normalize(subtract(verticesOnPlane[0], centroid));
+        const inPlaneBitangent = crossProduct(normal, inPlaneAxis);
+        const angleAboutNormal = (vertex: Vector3DSchema) => {
+            const offsetFromCentroid = subtract(vertex, centroid);
+            return Math.atan2(
+                dotProduct(offsetFromCentroid, inPlaneBitangent),
+                dotProduct(offsetFromCentroid, inPlaneAxis),
+            );
+        };
+        const ordered = [...verticesOnPlane].sort(
+            (left, right) => angleAboutNormal(left) - angleAboutNormal(right),
+        );
+        faces.push({ vertices: ordered, normal });
+    });
+
+    return faces.length >= 4 ? faces : null;
+}

src/js/lattice/reciprocal/lattice_reciprocal.ts

@@ -4,6 +4,7 @@ import { Utils } from "@mat3ra/utils";
 import lodash from "lodash";
 
 import { Lattice } from "../lattice";
+import { BrillouinZoneFace, computeBrillouinZone } from "./brillouin_zone";
 import { paths } from "./paths";
 import { symmetryPoints } from "./symmetry_points";
 
@@ -96,6 +97,20 @@ export class ReciprocalLattice extends Lattice {
         return symmetryPoints(this);
     }
 
+    /**
+     * Get the first Brillouin zone — the Wigner-Seitz cell of the reciprocal lattice — as a
+     * list of polygonal faces, ready to be projected and drawn.
+     *
+     * The zone follows from this lattice's own vectors rather than from its Bravais type, so
+     * materials sharing a type but differing in axial ratios (a bulk crystal and a slab with
+     * vacuum padding, say) yield correctly differing zones.
+     *
+     * @return {BrillouinZoneFace[] | null} null for a degenerate lattice.
+     */
+    get brillouinZone(): BrillouinZoneFace[] | null {
+        return computeBrillouinZone(this.reciprocalVectors);
+    }
+
     /**
      * Get the default path in reciprocal space for the current lattice.
      * @return {Array<{point: string; steps: number}>}

tests/js/lattice/brillouin_zone.ts

@@ -0,0 +1,103 @@
+import "../setup";
+
+import { LatticeSchema } from "@mat3ra/esse/dist/js/types";
+import { expect } from "chai";
+
+import { computeBrillouinZone } from "../../../src/js/lattice/reciprocal/brillouin_zone";
+import { ReciprocalLattice } from "../../../src/js/lattice/reciprocal/lattice_reciprocal";
+import { Graphene, Na4Cl4, Silicon, SiSlab } from "../fixtures";
+
+/** Distinct vertices across all faces, keyed by rounded coordinates. */
+function countVertices(faces: NonNullable<ReturnType<typeof computeBrillouinZone>>): number {
+    const keys = new Set<string>();
+    faces.forEach((face) =>
+        face.vertices.forEach((vertex) =>
+            keys.add(vertex.map((component) => component.toFixed(5)).join(",")),
+        ),
+    );
+    return keys.size;
+}
+
+function countEdges(faces: NonNullable<ReturnType<typeof computeBrillouinZone>>): number {
+    return faces.reduce((sum, face) => sum + face.vertices.length, 0) / 2;
+}
+
+function extentAlongThirdAxis(faces: NonNullable<ReturnType<typeof computeBrillouinZone>>): number {
+    const coordinates = faces.flatMap((face) => face.vertices.map((vertex) => vertex[2]));
+    return Math.max(...coordinates) - Math.min(...coordinates);
+}
+
+describe("Brillouin Zone", () => {
+    it("should be a truncated octahedron for a face-centered cubic lattice", () => {
+        const faces = new ReciprocalLattice(Silicon.lattice as LatticeSchema).brillouinZone;
+        expect(faces).to.not.be.null;
+        // 8 hexagons on the <111> planes and 6 squares on the <200> planes.
+        expect(faces).to.have.lengthOf(14);
+        expect(faces!.filter((face) => face.vertices.length === 6)).to.have.lengthOf(8);
+        expect(faces!.filter((face) => face.vertices.length === 4)).to.have.lengthOf(6);
+        expect(countVertices(faces!)).to.be.equal(24);
+    });
+
+    it("should be a cube for a simple cubic lattice", () => {
+        const faces = new ReciprocalLattice(Na4Cl4.lattice as LatticeSchema).brillouinZone;
+        expect(faces).to.not.be.null;
+        expect(faces).to.have.lengthOf(6);
+        expect(countVertices(faces!)).to.be.equal(8);
+        faces!.forEach((face) => expect(face.vertices).to.have.lengthOf(4));
+    });
+
+    it("should be a hexagonal prism for a hexagonal lattice", () => {
+        const faces = new ReciprocalLattice(Graphene.lattice as LatticeSchema).brillouinZone;
+        expect(faces).to.not.be.null;
+        expect(faces).to.have.lengthOf(8);
+        expect(faces!.filter((face) => face.vertices.length === 6)).to.have.lengthOf(2);
+        expect(faces!.filter((face) => face.vertices.length === 4)).to.have.lengthOf(6);
+    });
+
+    it("should follow the lattice itself, not only its type", () => {
+        // A slab pads the cell with vacuum along the third axis, which shrinks the
+        // corresponding reciprocal vector and flattens the zone.
+        const bulk = new ReciprocalLattice(Silicon.lattice as LatticeSchema).brillouinZone;
+        const slab = new ReciprocalLattice(SiSlab.lattice as LatticeSchema).brillouinZone;
+        expect(bulk).to.not.be.null;
+        expect(slab).to.not.be.null;
+        expect(extentAlongThirdAxis(slab!)).to.be.lessThan(extentAlongThirdAxis(bulk!));
+    });
+
+    it("should be a closed convex polyhedron", () => {
+        [Silicon, Na4Cl4, Graphene, SiSlab].forEach((material) => {
+            const faces = new ReciprocalLattice(material.lattice as LatticeSchema).brillouinZone;
+            expect(faces).to.not.be.null;
+            // Euler characteristic of a convex polyhedron: V - E + F = 2.
+            expect(countVertices(faces!) - countEdges(faces!) + faces!.length).to.be.equal(2);
+        });
+    });
+
+    it("should enclose the origin and no other reciprocal lattice point", () => {
+        const lattice = new ReciprocalLattice(Silicon.lattice as LatticeSchema);
+        const faces = lattice.brillouinZone!;
+        const [firstVector] = lattice.reciprocalVectors;
+        faces.forEach((face) => {
+            face.vertices.forEach((vertex) => {
+                const distanceToOrigin = Math.hypot(...vertex);
+                const distanceToNeighbour = Math.hypot(
+                    vertex[0] - firstVector[0],
+                    vertex[1] - firstVector[1],
+                    vertex[2] - firstVector[2],
+                );
+                expect(distanceToOrigin).to.be.at.most(distanceToNeighbour + 1e-6);
+            });
+        });
+    });
+
+    it("should return null for degenerate reciprocal vectors", () => {
+        expect(
+            computeBrillouinZone([
+                [1, 0, 0],
+                [1, 0, 0],
+                [0, 0, 1],
+            ]),
+        ).to.be.equal(null);
+        expect(computeBrillouinZone([[1, 0, 0]])).to.be.equal(null);
+    });
+});