This is what I've managed so far. In general, to convert to and from Cartesian coordinates, we need the the (Wilson's) $B$ matrix, a rectangular matrix whose elements are:
$$B_{ij}=\frac{\partial q_i}{\partial x_j}$$
where $q_i$ are the internal or collective coordinates, and $x_j$ are the Cartesian coordinates ($3N$ in total, $x$, $y$ and $z$ for each of the $N$ masses). This matrix can be used to convert derivatives and displacements. So from the vector of first derivatives (gradient) in Cartesian coordinates $g_x$, we can obtain the corresponding gradient in internal coordinates $g_q=(B^T)^+g_x$ (where $(B^T)^+$ is the Moore-Penrose pseudo inverse of the transpose of the $B$ matrix). A displacement in internal coordinates $\delta q$ can be obtained from a displacement of Cartesian coordinates: $\delta q=B\delta x$. For higher-order derivatives, the matrix of second derivatives $B'$ is needed:
$$B'_{ijk}=\frac{\partial^2 q_i}{\partial x_jx_k}$$
OK, now for the case of purely internal coordinates (distances, angles, dihedral angles), the expressions for the above $B$ and $B'$ elements are given in J. Chem. Phys. 117, 9160. But this assumes that the energy (function to minimize) is invariant with respect to translations and rotations, so the internal coordinates can describe all the degrees of freedom needed. If there is an external potential, field, charges, etc. this is not the case, and I must add translation and rotation of the whole system to the set of internal coordinates.
As I said, translation is easy, but not so much. At first I thought I would include the center of mass coordinates.
$$x_c=\sum\frac{m_ix_i}{M} \qquad M=\sum m_i$$
(now I use $x$ only for the $x$ coordinates, similar expressions would apply to $y_c$ and $z_c$) which gives:
$$\frac{\partial x_c}{\partial x_i} = \frac{m_i}{M} \qquad \frac{\partial x_c}{\partial y_i} = 0 \qquad \frac{\partial x_c}{\partial z_i} = 0$$
and similarly for $y_c$ and $z_c$. This is all fine, but there is a problem: it does not give the same displacement for all masses, i.e. a given displacement of the center of mass does not transform to the same displacement to each and every mass of the system (as I'd expect)[*]. In order to get this behaviour, I have to add not the center of mass but the geometrical center:
$$x_c=\sum\frac{x_i}{N} \qquad \frac{\partial x_c}{\partial x_i}=\frac{1}{N} \qquad \frac{\partial^2 x_c}{\partial x_i\partial x_i}=0$$
This ensures that when a displacement is converted to cartesian coordinates ($\delta x=B^+\delta q$) all masses will be displaced in the same way.
Now for rotation. Since I've used the geometrical center above, I consider rotation around this same center (so all Cartesian coordinates should be assumed to have $x_c,y_c,z_c$ subtracted). For a given mass, rotation around the $z$ axis can be intuitively given as $\arctan\frac{y}{x}$, or:
$$\frac{\partial R_z}{\partial x_i} = \frac{-y_i}{x_i^2+y_i^2} \qquad \frac{\partial R_z}{\partial y_i} = \frac{x_i}{x_i^2+y_i^2} \qquad \frac{\partial R_z}{\partial z_i} = 0$$
and symmetrically (cyclic) for $R_y$ and $R_x$. The second derivatives:
$$\frac{\partial^2 R_z}{\partial x_i\partial x_i} = \frac{2x_iy_i}{(x_i^2+y_i^2)^2} \qquad \frac{\partial^2 R_z}{\partial x_i\partial y_i} = \frac{y_i^2-x_i^2}{(x_i^2+y_i^2)^2} \qquad \frac{\partial^2 R_z}{\partial x_i\partial z_i} = 0$$
or, by analogy with translation, I could divide them all by $N$.
Applying these expressions in some test cases seemed to work, but I'm not sure if it was just luck. Besides, I have no idea what the rotational coordinates themselves, $R_x$, $R_y$, $R_z$, would be, I just used their derivatives. For any single mass I can get the angle with the axes, but for the whole system of $N$ masses? Computing the principal axes of inertia and the Euler angles of those is a possibility, but is this related to the derivatives above?
Does any of the above make sense?
[*] After reading Emilio Pisanty's answer, I think this is not a problem. To follow his example with two masses (assuming $m_1=1$ and $m_2=2$), if I remove $\xi_2=x_2-x_1$ as a coordinate, the $\frac{\partial x_i}{\partial \xi_j}$ matrix becomes the column vector $(0.6, 1.2)$. In this case, a given displacement of the center of mass $\Delta\xi_1$ will transform in a Cartesian displacement $\Delta x_1=0.6\Delta\xi_1$ and $\Delta x_2=1.2\Delta\xi_1$. Now suppose initially $x_1=0$, $x_2=1$, initially the center of mass would be at $\xi_1=\frac{2}{3}$; if the desired displacement is $\Delta\xi_1=0.5$, this gives $\Delta x_1=0.3$ and $\Delta x_2=0.6$, so the final coordinates would be $x_1=0.3$ and $x_2=1.6$, and the final center of mass is indeed $\xi_1=\frac{2}{3}+0.5=\frac{7}{6}$. But the distance between $x_1$ and $x_2$ has changed, which should be considered OK, as I didn't include the distance in the coordinates, which is like saying I don't care what happens with it. So, I guess my problem was that the system is underdetermined, and this shouldn't happen if I have at least the needed $3N$ coordinates.
