Tag: list

I am writing a code to generate the Wigner-Seitz cell of the reciprocal network for a given set of network translation vectors. For example, consider the cubic centered body network (BCC), whose base 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 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 network is then defined by the set of reciprocal network vectors, the set of all linear combinations of integer multiples of reciprocal base vectors, that is,

$$ 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 that contains the origin that is bounded by the perpendicular bisection planes of the reciprocal network vectors. In general, we can achieve this by considering only the first, second and perhaps the third point of the network closest 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}], two],
Rule[#] <= 2 d &];

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

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

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

(note that there is a redundancy for the origin, which simply gives True, this can be eliminated). Now, the problem will be to rewrite each of these equations as an inequality in such a way that the half space defined by the inequality contains the origin. I do not think it's too difficult, but not all equations can be solved for any of the coordinates, for example. we can not solve all the equations for $ z $, I like it

Solve[#, z] & / @ blueprints

Some of the equations must be solved for $ x $ or $ and $ Before becoming inequalities. I think I could find a brute force solution but I hope there is something more elegant.

Ultimately, I would like to obtain the inequalities that define the region in order to visualize it with RegionPlot3D and use it to Select points of a mesh.

To update:

So far I made a single Yes build to solve the equations for anyone $ z $, $ x $or $ and $,

sols = flat

The
Yes[sz=Solve[sz=Solve[sz=Resolver[sz=Solve[#, z]; sz! = {}, sz,
Yes[sx = Solve[#, x]; sx! = {}, sx, Resolve[#, y]]]& / @ blueprints]

which results in the following list of rules

{z -> -1, z -> -1 - y, z -> -1 + x, y -> -1, x -> 1 + y, x -> 1, z -> -1 - x, z -> -1 + y, x -> -1 - y, x -> 1 - y, z -> 1 + y, z -> 1 - x, x -> -1, x -> -1 + y, and -> 1, z -> 1 + x, z -> 1 - y, z -> 1}

Now I just need to convert each of these into $ LHS <RHS $ Y $ LHS> RHS $ and test which contains the origin, but I still do not know how to achieve it.