We present a hard constraint, linear complementarity based, method for the simulation of stiff multibody dynamics with contact, joints, and friction. The approach uses a linearization of the modified trapezoidal method, incorporates a Poisson restitution model at collision, and solves only one linear complementarity problem per time step when no collisions are encountered. We prove that the method has order two, a fact that is also demonstrated by our numerical simulations. When we use a special approximation of the Jacobian matrix for the case where the stiff forces originate in springs and dampers attached to two points in the system, we prove that the linear complementarity problem can be solved for any value of the time step and that the method is stiffly stable. In particular, the numerical solution converges to the numerical solution of the constraint problem, where the stiff springs and dampers have been replaced by rigid links, a fact that we also demonstrate by numerical simulations. The method was implemented in UMBRA, an industrial-grade virtual prototyping software.