(* First choose three lattice vectors that form a parallelepied with unit volume *)
u1={1,1,1};
u2 = {3,0,1};
u3= {1, 0, 0};
U = {u1, u2, u3};
If[Det[U]==1, Print["Okay."], Print["Lattice vectors are invalid."]]

(* Shear the parallelepiped into a cuboid *)
α=-Dot[u1,u2]/Dot[u1,u1];
γ=-(α Dot[u1,u3] + Dot[u2,u3])/(α Dot[u1,u2] + Dot[u2,u2]);
β=-(Dot[u1,u3] + γ Dot[u1,u2])/Dot[u1,u1];
e1=u1
e2=u2 + α u1
e3=u3 + β u1 + γ u2

(* Check that all edges are in fact orthogonal *)
If[Dot[e1,e2]==0 && Dot[e2,e3]==0 && Dot[e3,e1]==0, Print["Okay."], Print["Edges are not orthogonal."]]

(* Compute the dimensions of the cuboid *)
L1=Norm[e1];
L2=Norm[e2];
L3=Norm[e3];
Print[L1,"×",L2,"×",L3]
Print[L1//N," × ",L2//N," × ",L3//N]

(* Get the minimum and maximum coordinates of the cuboid in each direction *)
v={{0,0,0},e1,e2,e1+e2,e3,e1+e3,e2+e3,e1+e2+e3};
xvals=Table[v[[i,1]],{i,1,8}];
yvals=Table[v[[i,2]],{i,1,8}];
zvals=Table[v[[i,3]],{i,1,8}];
xmin=Min[xvals];
xmax=Max[xvals];
ymin=Min[yvals];
ymax=Max[yvals];
zmin=Min[zvals];
zmax=Max[zvals];

(* Compute the axis-angle rotation from an axis-aligned cuboid to our cuboid *)
n1=e1/L1;
n2=e2/L2;
n3=e3/L3;
R={n1,n2,n3};
θ=ArcCos[((Tr[R]-1.)/2)];
ω=(1/(2 Sin[θ])){R[[3,2]]-R[[2,3]], R[[1,3]]-R[[3,1]], R[[2,1]] -R[[1,2]]};

(* Draw the cuboid, with the origin and first edge highlighted, and the lattice grid in gray *)
cbd=Cuboid[{0,0,0},{L1,L2,L3}];
rotcbd=Rotate[cbd,-θ,ω];
grid=Table[Cuboid[{ix,iy,iz}],{ix,Floor[xmin],Ceiling[xmax-1]},{iy,Floor[ymin],Ceiling[ymax-1]},{iz,Floor[zmin],Ceiling[zmax-1]}];
origin={0,0,0};
edge1=Line[{origin,e1}];
edge2=Line[{origin,e2}];
edge3=Line[{origin,e3}];
Graphics3D[{Opacity[0.4],rotcbd,Opacity[1],Thick,edge1,PointSize[Large],Point[origin],FaceForm[],EdgeForm[Opacity[0.1]],grid}, {Boxed->False}]