  We introduce two new thermostats, one of Langevin type and one of gradient
(Brownian) type, for rigid body dynamics. We formulate rotation using the
quaternion representation of angular coordinates; both thermostats preserve the
unit length of quaternions. The Langevin thermostat also ensures that the
conjugate angular momenta stay within the tangent space of the quaternion
coordinates, as required by the Hamiltonian dynamics of rigid bodies. We have
constructed three geometric numerical integrators for the Langevin thermostat
and one for the gradient thermostat. The numerical integrators reflect key
properties of the thermostats themselves. Namely, they all preserve the unit
length of quaternions, automatically, without the need of a projection onto the
unit sphere. The Langevin integrators also ensure that the angular momenta
remain within the tangent space of the quaternion coordinates. The Langevin
integrators are quasi-symplectic and of weak order two. The numerical method
for the gradient thermostat is of weak order one. Its construction exploits
ideas of Lie-group type integrators for differential equations on manifolds. We
numerically compare the discretization errors of the Langevin integrators, as
well as the efficiency of the gradient integrator compared to the Langevin ones
when used in the simulation of rigid TIP4P water model with smoothly truncated
electrostatic interactions. We observe that the gradient integrator is
computationally less efficient than the Langevin integrators. We also compare
the relative accuracy of the Langevin integrators in evaluating various static
quantities and give recommendations as to the choice of an appropriate
integrator.
