The velocity Verlet algorithm is broken up into two parts where one alternates between updating the positions and velocities of particles. As pointed out by Koelman and Hoogerbrugge (1993), there is an additional contribution to the stress tensor due to the constraint forces that maintain the relative positions of particles the rigid body are composed of. Clearly, we are using an algorithm that does not explicitly determine the constraint forces on each particle of the rigid body. However, since we know the position and velocity of each particle, we can effectively backout the constraint forces by mapping the individual particle's (contained on the rigid body trajectory to the velocity Verlet algorithm and then solve for the constraint forces that would be required for such motion of the particles to take place. Note that the constraint forces used to update the position and then the velocities are not the same. The difference is actually small and ignoring either contribution alone results in an error of order δt3 in the stress tensor (which should be symmetric up to that order). Although this agreement seems reasonable, a symmetric stress tensor is required to demonstrate that angular momentum is conserved. By incorporating the contributions of the constraint forces from the two steps of the velocity Verlet algorithm, it was found that the stress tensor was symmetric up to the order of precision of the computer (i.e., 16 figures for double precision).
We also compared our approach to determining the stress tensor to an entirely different but commonly used method described in Allen and Tildesley (1987). Here the constraint forces are not used to determine the stress tensor. Instead, the center of mass force each rigid body has on the other is utilized [cf. Allen and Tildesley (1987)]. This approach has the undesirable feature that it does not produce a symmetric tensor for short times but when averaged over long times approaches the correct symmetry. It was found that the time average stress tensor determined by both approaches were reasonably close in value and that they asymptotically approached each other over time.