This work presents a novel methodology for the dynamic analysis of general non-linear multibody systems composed of rigid and deformable bodies, the latter under the small strain assumption. In Part I we developed the 6-D compact representation and parameterization of motion for constrained bodies. Part II is devoted to the design of a class of modified Runge–Kutta (RK) methods dedicated to non-linear dynamics. These are capable of integrating on the configuration manifold and of preserving linear and angular momenta. Within this class of methods, two second-order algorithms are designed under the requirement of attaining non-linear unconditional stability: the energy preserving (EP) and energy decaying (ED) methods. These schemes are associated with an algorithmic law of conservation and dissipation, respectively, of the total mechanical energy of the system, together with the vanishing of the algorithmic work done by ideal, time-independent constraints. Their performances are assessed with the aid of some representative numerical applications which confirm the non-conventional properties predicted in the analysis.