Integration of angular velocity of nonspherical particles is wrong
As explained in "Particlescale computational approaches to model dry and saturated granular flows of nonBrownian, noncohesive, and nonspherical rigid bodies" by Wachs (2019), the equation of rotational rigid particle motion, formulated in the body frame, generally reads (Eq. 27):
\frac{d (J \omega)}{d t} = J \frac{d (\omega)}{d t} + \omega \times J \omega = T
The second term featuring the cross product is nonzero if the particle is nonspherical. However, this term is missing in all integrator kernels of PE and mesapd (e.g. used in combination with DEMlike collision models).
On the other hand, the HCSITS uses own integration procedures and should also be checked to ensure correctness. There, a difference in the peHCSITS (pe/cr/HCSITS.impl.h, line 308):
bodyCache.w_[j] = body>getAngularVel() + dt * ( body>getInvInertia() * ( ( body>getInertia() * body>getAngularVel() ) % body>getAngularVel() ) );
and the mesapdHCSITS (mesa_pd/kernel/InitParticlesForHCSITS.h)
ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * ac.getAngularVelocity(j) ) % ac.getAngularVelocity(j) ) );
has been noted. While the pe uses the worldframe (WF) consistently, mesa_pd mixes WF and body frame (BF) which seems wrong.
Possible solution (assuming that pe HCSITS variant is correct)

add getInertia
andgetInvInertia
utility functions to mesa_pd (adapt generation files). See pe/rigidbody/RigidBody.h for implementation. 
Use in XWithShape.h integrators (adapt generation files) according to pe formulation. For VelocityVerlet, it could be that more advanced algorithm needs to be use, see e.g. https://yadedem.org/doc/formulation.html#orientationaspherical 
Fix InitParticlesForHCSITS.h accordingly.