In general, one cannot expect rk4 to maintain quadratic invariants of the system, it simply doesn't do that. Methods that do maintain specific invariants have to be specially devised — this usually goes by the name of geometric integration. Symplectic integrators are the most common kind of these methods, although not for your problem.
See, for example, the book Geometric Numerical Integration by Hairer, Wanner, Lubich, or the notes Solving Differential Equations on Manifolds by Hairer (free online).
Probably the simplest way to incorporate the constraint is to project the approximate solution onto the constraint manifold ("Projection methods", Example 1.3, p. 21 in the Hairer notes above). Given a system of ODEs
$$ \dot X = f(t, X), \qquad X^*X=I, $$
define $\pi(X)$ to be the closest unitary matrix to $X$. If $X = UH$ is the polar decomposition of $X$, then $\pi(X) = U$.
The simplest is to project onto the constraint manifold after each RK step, $\tilde X_{n+1} = \pi(X_{n+1})$. The second simplest is to rewrite the ODE as
$$ \dot X = f(t, \pi(X)), $$
which will incorporate the projection into the RK substeps. (This doesn't actually modify the ODE because $\pi(X)=X$ for any exact solution.)
Beyond this, the constraint manifold is called a Stiefel manifold, and there already seem to be some number of articles on google scholar on numerical integration on Stiefel manifolds.
