The current isotropic initial condition uses scipy.spatial.transform.Rotation.random which seems to be constructed by sampling normal distributions in quaternion space (scipy Rotations are implemented in Cython in terms of quaternions). This does not take into account crystallite symmetry and may not produce an appropriately dispersed sampling in SO(3). The Matlab DRex code uses pseudorandom Halton sequences instead. There are other more sophisticated options like using a Bingham distribution which is axial (-x = x) or UARS distributions which are specifically designed for crystal orientations.