I prefer to use a rotation matrix
approach to do coordinate transformations. But you
still need the appropriate angles to plug in.

I am using a rotation matrix. Why did you think I was not?

Semi major axis, eccentricity, inclination, longitude of ascending node,
longitude of perihelion, mean longitude [at epoch].

So, for your circular orbit the longitude of perihelion
will be undefined (no perihelion). If it happens to be an
equitorial orbit, you lose the ascending node, too.

Just to be clear, I am going from the position and velocity at a given time,
to the orbital elements.

The mean longitude at epoch will help. For the circular
orbit it should be the same as the true longitude at epoch,
which is the angle between the Vernal Equinox and the
radius vector at epoch, measured eastwards to the
ascending node, and then in the orbital plane to the
radius vector.

....

The problem is the lack of the perihelion; there's
no point of reference for the true anomaly for a circular
orbit without explicity defining one. Often the ascending
node is chosen. For circular equatorial orbits, the
Vernal Equinox (I axis in the IJK frame) is often used.

The ijk frame you mention, is that in the 'q-plane' of the orbit (the orbit
before the rotations that would be applied by I,w and Om) or the 'real' axis
system?

let's take your case of a circular orbit where you are
given the semimajor axis (orbit radius, a), inclination (i),
longitude of ascending node (Om), and mean longitude at epoch
(M).

Sorry, going this way round is easy. It's going back that's the problem.
If I'm given a position and a velocity at a point on a elliptical orbit,
everything's fine. If the orbit turns out to be circular though, then I
can't calculate those undefined parameters!
Presumably, in this case (e=0) I could just set them to zero?