From f4815581f581ac984e56573f3b663ff46a869a2e Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 16 Aug 2026 20:32:16 +0000 Subject: [PATCH 1/2] feat: compute the first Brillouin zone from the reciprocal lattice MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds ReciprocalLattice.brillouinZone, returning the Wigner-Seitz cell of the reciprocal lattice as polygonal faces: the intersection of the half-spaces k.G <= |G|^2/2, with vertices at the plane triple intersections that satisfy every other half-space. The zone follows from the lattice's own vectors rather than its Bravais type, so materials sharing a type but differing in axial ratios — a bulk crystal and a slab padded with vacuum, say — yield correctly differing zones. Consumers currently illustrate the zone with one static image per lattice type, which cannot express that difference (and ships as an asset outside the package). Validated against the material fixtures: silicon gives a truncated octahedron, Na4Cl4 a cube, graphene a hexagonal prism, each closed with V - E + F = 2 and enclosing no reciprocal lattice point but the origin. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01FFFpkKD3zHysJXgAbT2uYV --- .../js/lattice/reciprocal/brillouin_zone.d.ts | 28 +++ dist/js/lattice/reciprocal/brillouin_zone.js | 142 +++++++++++++ .../reciprocal/lattice_reciprocal.d.ts | 12 ++ .../lattice/reciprocal/lattice_reciprocal.js | 14 ++ src/js/lattice/reciprocal/brillouin_zone.ts | 199 ++++++++++++++++++ .../lattice/reciprocal/lattice_reciprocal.ts | 15 ++ tests/js/lattice/brillouin_zone.ts | 103 +++++++++ 7 files changed, 513 insertions(+) create mode 100644 dist/js/lattice/reciprocal/brillouin_zone.d.ts create mode 100644 dist/js/lattice/reciprocal/brillouin_zone.js create mode 100644 src/js/lattice/reciprocal/brillouin_zone.ts create mode 100644 tests/js/lattice/brillouin_zone.ts diff --git a/dist/js/lattice/reciprocal/brillouin_zone.d.ts b/dist/js/lattice/reciprocal/brillouin_zone.d.ts new file mode 100644 index 000000000..795849c4e --- /dev/null +++ b/dist/js/lattice/reciprocal/brillouin_zone.d.ts @@ -0,0 +1,28 @@ +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; +} +/** + * 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 declare function computeBrillouinZone(reciprocalVectors: Vector3DSchema[]): BrillouinZoneFace[] | null; diff --git a/dist/js/lattice/reciprocal/brillouin_zone.js b/dist/js/lattice/reciprocal/brillouin_zone.js new file mode 100644 index 000000000..94baf98f7 --- /dev/null +++ b/dist/js/lattice/reciprocal/brillouin_zone.js @@ -0,0 +1,142 @@ +"use strict"; +Object.defineProperty(exports, "__esModule", { value: true }); +exports.computeBrillouinZone = computeBrillouinZone; +/** + * 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, second) { + 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, second) { + return first[0] * second[0] + first[1] * second[1] + first[2] * second[2]; +} +function subtract(first, second) { + return [first[0] - second[0], first[1] - second[1], first[2] - second[2]]; +} +function scale(vector, factor) { + return [vector[0] * factor, vector[1] * factor, vector[2] * factor]; +} +function vectorLength(vector) { + return Math.sqrt(dotProduct(vector, vector)); +} +function normalize(vector) { + 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, rightHandSide) { + const determinant = dotProduct(matrix[0], crossProduct(matrix[1], matrix[2])); + if (Math.abs(determinant) < SINGULAR_MATRIX_TOLERANCE) { + return null; + } + const determinantWithColumnReplaced = (columnIndex) => { + const replaced = matrix.map((row, rowIndex) => { + const nextRow = [...row]; + nextRow[columnIndex] = rightHandSide[rowIndex]; + return nextRow; + }); + 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). + */ +function computeBrillouinZone(reciprocalVectors) { + 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 = []; + 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) => planes.every((plane) => dotProduct(point, plane.normal) <= plane.offset + HALF_SPACE_TOLERANCE); + const vertices = []; + 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)) < + COINCIDENT_POINT_TOLERANCE); + if (isZoneVertex && !isDuplicate) { + vertices.push(point); + } + } + } + } + if (vertices.length < 4) { + return null; + } + const faces = []; + 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((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) => { + 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; +} diff --git a/dist/js/lattice/reciprocal/lattice_reciprocal.d.ts b/dist/js/lattice/reciprocal/lattice_reciprocal.d.ts index 20e8ffc16..92753c640 100644 --- a/dist/js/lattice/reciprocal/lattice_reciprocal.d.ts +++ b/dist/js/lattice/reciprocal/lattice_reciprocal.d.ts @@ -1,5 +1,6 @@ import { Vector3DSchema } from "@mat3ra/esse/dist/js/types"; import { Lattice } from "../lattice"; +import { BrillouinZoneFace } from "./brillouin_zone"; export type KPointCoordinates = number[]; export type KPointPath = Array<{ point: string; @@ -40,6 +41,17 @@ export declare class ReciprocalLattice extends Lattice { * @return {SymmetryPoint[]} */ get symmetryPoints(): SymmetryPoint[]; + /** + * 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; /** * Get the default path in reciprocal space for the current lattice. * @return {Array<{point: string; steps: number}>} diff --git a/dist/js/lattice/reciprocal/lattice_reciprocal.js b/dist/js/lattice/reciprocal/lattice_reciprocal.js index 38886ba87..b1103f5df 100644 --- a/dist/js/lattice/reciprocal/lattice_reciprocal.js +++ b/dist/js/lattice/reciprocal/lattice_reciprocal.js @@ -8,6 +8,7 @@ const constants_1 = require("@mat3ra/code/dist/js/constants"); const utils_1 = require("@mat3ra/utils"); const lodash_1 = __importDefault(require("lodash")); const lattice_1 = require("../lattice"); +const brillouin_zone_1 = require("./brillouin_zone"); const paths_1 = require("./paths"); const symmetry_points_1 = require("./symmetry_points"); class ReciprocalLattice extends lattice_1.Lattice { @@ -54,6 +55,19 @@ class ReciprocalLattice extends lattice_1.Lattice { get symmetryPoints() { return (0, symmetry_points_1.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() { + return (0, brillouin_zone_1.computeBrillouinZone)(this.reciprocalVectors); + } /** * Get the default path in reciprocal space for the current lattice. * @return {Array<{point: string; steps: number}>} diff --git a/src/js/lattice/reciprocal/brillouin_zone.ts b/src/js/lattice/reciprocal/brillouin_zone.ts new file mode 100644 index 000000000..65929caef --- /dev/null +++ b/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( + (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; +} diff --git a/src/js/lattice/reciprocal/lattice_reciprocal.ts b/src/js/lattice/reciprocal/lattice_reciprocal.ts index 298145204..b52c7a3f1 100644 --- a/src/js/lattice/reciprocal/lattice_reciprocal.ts +++ b/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}>} diff --git a/tests/js/lattice/brillouin_zone.ts b/tests/js/lattice/brillouin_zone.ts new file mode 100644 index 000000000..a325d32e5 --- /dev/null +++ b/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>): number { + const keys = new Set(); + faces.forEach((face) => + face.vertices.forEach((vertex) => + keys.add(vertex.map((component) => component.toFixed(5)).join(",")), + ), + ); + return keys.size; +} + +function countEdges(faces: NonNullable>): number { + return faces.reduce((sum, face) => sum + face.vertices.length, 0) / 2; +} + +function extentAlongThirdAxis(faces: NonNullable>): 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); + }); +}); From b544032dcdd0d400ad9d90ea68902106c9d5b881 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 27 Aug 2026 03:36:17 +0000 Subject: [PATCH 2/2] chore: stop committing dist/ (build output; see mat3ra/agents#5) --- .../js/lattice/reciprocal/brillouin_zone.d.ts | 28 ---- dist/js/lattice/reciprocal/brillouin_zone.js | 142 ------------------ 2 files changed, 170 deletions(-) delete mode 100644 dist/js/lattice/reciprocal/brillouin_zone.d.ts delete mode 100644 dist/js/lattice/reciprocal/brillouin_zone.js diff --git a/dist/js/lattice/reciprocal/brillouin_zone.d.ts b/dist/js/lattice/reciprocal/brillouin_zone.d.ts deleted file mode 100644 index 795849c4e..000000000 --- a/dist/js/lattice/reciprocal/brillouin_zone.d.ts +++ /dev/null @@ -1,28 +0,0 @@ -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; -} -/** - * 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 declare function computeBrillouinZone(reciprocalVectors: Vector3DSchema[]): BrillouinZoneFace[] | null; diff --git a/dist/js/lattice/reciprocal/brillouin_zone.js b/dist/js/lattice/reciprocal/brillouin_zone.js deleted file mode 100644 index 94baf98f7..000000000 --- a/dist/js/lattice/reciprocal/brillouin_zone.js +++ /dev/null @@ -1,142 +0,0 @@ -"use strict"; -Object.defineProperty(exports, "__esModule", { value: true }); -exports.computeBrillouinZone = computeBrillouinZone; -/** - * 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, second) { - 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, second) { - return first[0] * second[0] + first[1] * second[1] + first[2] * second[2]; -} -function subtract(first, second) { - return [first[0] - second[0], first[1] - second[1], first[2] - second[2]]; -} -function scale(vector, factor) { - return [vector[0] * factor, vector[1] * factor, vector[2] * factor]; -} -function vectorLength(vector) { - return Math.sqrt(dotProduct(vector, vector)); -} -function normalize(vector) { - 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, rightHandSide) { - const determinant = dotProduct(matrix[0], crossProduct(matrix[1], matrix[2])); - if (Math.abs(determinant) < SINGULAR_MATRIX_TOLERANCE) { - return null; - } - const determinantWithColumnReplaced = (columnIndex) => { - const replaced = matrix.map((row, rowIndex) => { - const nextRow = [...row]; - nextRow[columnIndex] = rightHandSide[rowIndex]; - return nextRow; - }); - 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). - */ -function computeBrillouinZone(reciprocalVectors) { - 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 = []; - 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) => planes.every((plane) => dotProduct(point, plane.normal) <= plane.offset + HALF_SPACE_TOLERANCE); - const vertices = []; - 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)) < - COINCIDENT_POINT_TOLERANCE); - if (isZoneVertex && !isDuplicate) { - vertices.push(point); - } - } - } - } - if (vertices.length < 4) { - return null; - } - const faces = []; - 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((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) => { - 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; -}