Xjyuk

I'm writing a code to generate the Wigner-Seitz cell of the reciprocal lattice for a given set of lattice translation vectors. For example, consider the Body Centered Cubic (BCC) lattice whose basis translation vectors are given by

a1 = -1, 1, 1/2;
a2 = 1, -1, 1/2;
a3 = 1, 1, -1/2;

The reciprocal basis vectors are then defined according to

d = 2 Pi;
v = a1.(a2[Cross]a3);
b1 = d/v (a2[Cross]a3);
b2 = d/v (a3[Cross]a1);
b3 = d/v (a1[Cross]a2);

The reciprocal lattice is then defined by the set of reciprocal lattice vectors, the set of all linear combinations of integer multiples of reciprocal basis vectors, i.e.

$$vecG = n_1 vecb_1 + n_2 vecb_2 + n_3 vecb_3, qquad n_i in mathbbZ$$

The Wigner-Seitz cell (in this case the First Brillouin Zone) is defined as the region containing the origin which is bounded by the perpendicular bisecting planes of the reciprocal lattice vectors. We generally can accomplish this by only considering the first, second, and maybe third closest reciprocal lattice points to the origin. In the case of BCC, for example, the following vectors will suffice:

recipvecs = 
 Select[Flatten[
 Table[n1 b1 + n2 b2 + n3 b3, n1, -1, 1, n2, -1, 1, n3, -1, 1], 2],
 Norm[#] <= 2 d &];

Question: Given these vectors, how can I construct the Wigner-Seitz cell?

For example, one possibility is to construct the equations for all the planes

planes = (x, y, z - (#/2)).# == 0 & /@ reciplattice

(note there is a redundancy for the origin, which just gives True, this can be removed). Now the issue is going to be to rewrite each of these equations as an inequality such that the half-space defined by the inequality contains the origin. I don't think that would be too difficult, but not every one of the equations can be solved for any one of the coordinates, e.g. we cannot solve every equation for $z$, like

Solve[#, z] & /@ planes

Some of the equations will have to be solved for $x$ or $y$ before being turned into inequalities. I think I could find a brute force solution but I'm hoping there's something more elegant.

Ultimately I'd like to obtain the inequalities that define the region so that I can visualize it with RegionPlot3D and use it to Select points from a mesh.

Unfortunately, VoronoiMesh does not work in 3D. So we do it manually.

the crystal lattice vectors:

a1 = -1, 1, 1/2;
a2 = 1, -1, 1/2;
a3 = 1, 1, -1/2;

the reciprocal lattice vectors: (Inverse is easier than using cross products, but ultimately the same thing)

B = b1, b2, b3 = 2π*Inverse[Transpose[a1, a2, a3]];

an inequality defining the perpendicular bisecting plane of a reciprocal lattice point v:

pbp[0, 0, 0, r_] = True;
pbp[v_, r_] := v.r/v.v <= 1/2

make a list of such inequalities, And them, and simplify: (here you may have to go to larger s to get all the constraints, as you said)

With[s = 1,
 WS[x_, y_, z_] = FullSimplify[
 And @@ Flatten[Table[pbp[n1,n2,n3.B, x,y,z], n1,-s,s, n2,-s,s, n3,-s,s]]]]

-2 π <= y + z <= 2 π && z <= 2 π + x && x <= 2 π + z && y <= 2 π + x && x <= 2 π + y && -2 π <= x + z <= 2 π && z <= 2 π + y && y <= 2 π + z && -2 π <= x + y <= 2 π

make a 3D plot of the Wigner-Seitz cell: (use more PlotPoints to make it prettier)

With[t = 2π,
 RegionPlot3D[WS[x, y, z], x, -t, t, y, -t, t, z, -t, t]]

You can also check if a point is in the Wigner-Seitz cell or not:

WS[0.1, 0.2, 0.3]
(* True *)
WS[3.1, 3.2, 0.3]
(* False *)

It is really unfortunate that we don't have a 3D implementation of VoronoiMesh.

Borrowing quite a lot from Roman, the following tries to compute the extremal points of the Wigner-Seitz cells and applies ConvexHullMesh to the result in order to obtain the precise polyhedron.

a1 = -1, 1, 1/2;
a2 = 1, -1, 1/2;
a3 = 1, 1, -1/2;
B = b1, b2, b3 = 2 π*Inverse[Transpose[a1, a2, a3]];
pts = Flatten[Table[b1, b2, b3.n1, n2, n3, n1, -1, 1, n2, -1, 1, n3, -1, 1], 2];
G = NearestNeighborGraph[pts, VertexCoordinates -> pts];
neighbors = Rest[VertexOutComponent[G, 0, 0, 0, 1]];
rhs = MapThread[Dot, neighbors, neighbors]/2;
subsets = Subsets[Range[Length[neighbors]], 3];

q = Module[A, x,
 Table[
 A = neighbors[[s]];
 If[Det[A] != 0,
 x = LinearSolve[A, rhs[[s]]];
 If[And @@ Thread[neighbors.x <= rhs], x, Nothing],
 Nothing
 ],
 s, subsets]
 ];
R = ConvexHullMesh[q]

The other answers are great and very enlightening, I had already found a brute force solution but I took elements of both @Henrik Schumacher and @Roman's answers to produce this nice minimal one for what I wanted. I think both of their answers are better in that they provide more functionality.

d = 2 Pi;
a1 = -1, 1, 1/2;
a2 = 1, -1, 1/2;
a3 = 1, 1, -1/2;
b1, b2, b3 = d*Inverse[Transpose[a1, a2, a3]];
reciplattice = 
 Select[Flatten[
 Table[n1 b1 + n2 b2 + n3 b3, n1, -1, 1, n2, -1, 1, n3, -1, 1], 
 2], 0 < Norm[#] <= 2 d &];
region = And@@FullSimplify[(x, y, z - (#/2)).# <= 0 & /@ reciplattice] 

And plotting it with

e = d + 0.1;
fbz = RegionPlot3D[region, x, -e, e, y, -e, e, z, -e, e, 
 PlotPoints -> 60]