  We present a code to construct initial data for binary neutron star systems
in which the stars are rotating. Our code, based on a formalism developed by
Tichy, allows for arbitrary rotation axes of the neutron stars and is able to
achieve rotation rates near rotational breakup. We compute the neutron star
angular momentum through quasi-local angular momentum integrals. When
constructing irrotational binary neutron stars, we find a very small residual
dimensionless spin of $\sim 2\times 10^{-4}$. Evolutions of rotating neutron
star binaries show that the magnitude of the stars' angular momentum is
conserved, and that the spin- and orbit-precession of the stars is well
described by post-Newtonian approximation. We demonstrate that orbital
eccentricity of the binary neutron stars can be controlled to $\sim 0.1\%$. The
neutron stars show quasi-normal mode oscillations at an amplitude which
increases with the rotation rate of the stars.
