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.
$$\vec{G} = n_1 \vec{b}_1 + n_2 \vec{b}_2 + n_3 \vec{b}_3, \qquad n_i \in \mathbb{Z}$$
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.