Hi @Gianluca the Sivak paper (as well as many others, including e.g. here and here) deals with differences in how the stochastic forces are discretized. This is a different issue from using or not the velocity-Verlet algorithm, which is assumed the best second-order scheme for the special case without stochastic forces.
Thankfully, GJF and all schemes equivalent to it yield the same configurational sampling, which is the same fact already pointed out by @Germain . There may be some differences in the effective diffusivity (time scales), but not in the position ensemble.
Separately from positions, there is the issue of how to define instantaneous velocities in a way that is numerically stable in the high-friction limit. This is discussed in this paper, and implemented with the vhalf option of fix langevin. I and colleagues have experience this issue ourselves when looking at stress tensor components (which include a velocity term).
I have no experience with applying these schemes on granular bodies or iFLD, but hopefully with the help of others you can figure out if the same implementation applies to your setup. See also this previous discussion.
Giacomo