Adding rotation to internal coordinates I'm optimizing the geometry of a system composed of several interacting masses (a molecule). The energy of the system depends on the relative position of the masses, and all velocities are zero.
In a simple case, the energy is invariant with translation and rotation of the whole system, so I can optimized using internal coordinates, defining these internal coordinates as the different distances, angles and dihedral angles between the masses. For converting coordinates to and from Cartesian I have the matrix of coordinate derivatives, etc. J. Chem. Phys. 117, 9160 (there are small mistakes in eqs. 34 and 35).
Now the question is how to proceed when the energy is not invariant with translation and rotation. I have to include translation and rotation of the whole system to the degrees of freedom of the optimization. Translation is easy, because I can add the Cartesian coordinates of the center of mass, and the derivatives with respect to the coordinates of the masses are easy. But how can I add rotation?
I guess it will be related to the orientation of the principal axes of inertia, but I can't find a clear expression that will give me the three rotational coordinates I need, and their derivatives with respect to all the masses' coordinates.
Can anyone give me some help or pointers?
 A: EDIT: Later addition
This project is ill conceived, you should put the system in rectangular absolute coordinates, there is absolutely no gain from what you are doing, it is objectively wrong, and needlessly complicated.
However, you are specifying the location of point particles one after the other, in relative polar coordiantes, so what you can do to fix the orientation of the whole thing is add 1 fictitious new particle, which is attached to the next atom with an angle, and gives you an orientation angle for the whole thing. This is probably the simplest modification.
I have to say that you're wasting your time in writing code for point particles using relative coordinates. This is deranged.
For rigid bodies
You should use the orientation of one of the components as defining the orientation of the whole molecule. There is no analog of the center of mass for rotations.
The total energy function can be written in terms of the rotation matrices $R_i$ which rotate the parts from some initial choice of orientation to the one you are looking at, plus the position of the center of mass.
$$ H(x_1,R_1,x_2,R_2....,x_n,R_n)$$
It is invariant with respect to translations and rotations
$$ H(RR_i,x_i+c) = H(R_i,x_i)$$
The center of mass coordinates give you a center for x, but there is no analogous orientation average, because the space of orientations is not infinite in extent. So instead of using an orientation average, you just use the orientation of one object, say object 1, to fix the orientation. Then the R matrix for object 1 defines the global orientation, and the rest of the objects orientations are relative to the first.
The energy is then a function of R, the orientation of object 1, and of the relative orientations of the other objects to object 1, $R_1^{-1} R_i$ and the positions relative to the center of mass. The relative orientation between object i and object j is $R_i^{-1}R_j$, and this relative orientation is independent of whether you use the R matrix relative to the frame defined by object 1, or relative to the original frame. The whole object can be reconstructed by placing the center of mass somewhere, placing object 1 at the right position with orientation R, then the rest of the objects relative to the axes defined by object 1.
You shouldn't use the moment of inertia principal axes as your definition of the orientation of the whole thing, because these become discontinuous when two of the moments of inertia become equal, so that the object will rotate discontinuously at those times when two moments happen to become equal.
A: 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.
A: This is slightly too long, and requires a bit too much ${\LaTeX}$, for a comment.
"a given displacement of the center of mass does not transform to the same displacement to each and every mass of the system" is not quite right. The maths you are displaying say the opposite: a small displacement $\delta x_i$ on molecule $i$ effects a displacement on $x_c$ which depends on $m_i$, which is right and perfectly physically understandable. Your statement in words would describe the matrix $\frac{\partial x_i}{\partial x_c}$, for which there is not enough information. Specifically, this matrix is the inverse of the matrix you're describing, and to get that inverse, you need the whole matrix, i.e. you need the rest of the internal coordinates in order to get the effect of $\delta x_c$ on an individual position $x_i$.
To be more specific, let me draw an analogy with two masses on a line, with coordinates $x_1$ and $x_2$. You want one of your final coordinates to be the centre of mass,
$$\xi_1=x_c=\frac{m_1}{M}x_1+\frac{m_2}{M}x_2,$$
and this fixes one of the rows of the matrix $\frac{\partial \xi_j}{\partial x_i}$. The other row is still free: there is of course the canonical choice
$$\xi_2=x_r=x_2-x_1$$
but other choices are of course possible such as $(x_2-x_1)^3$, for one, or even $x_1$ or $x_2$; note that the latter are perfectly valid as a change of coordinates, but they radically change the physical content of the transformation.
To make this clearer, consider the matrix you quote, with $x_r$ chosen as $\xi_2$:
$$\frac{\partial\xi_j}{\partial x_i}=\begin{pmatrix}m_1/M &m_2/M\\ -1&1 \end{pmatrix}.$$
Then the determinant is 1 and the inverse is
$$\frac{\partial x_i}{\partial \xi_j}= \begin{pmatrix}1&-m_2/M\\1&m_1/M \end{pmatrix}.$$
If you now focus on the first column, it tells you that $\frac{\partial x_i}{\partial x_c}=1$ for $i=1,2$ and it is this that translates into your statement

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)

You can see that it does: a given displacement of the centre of mass does transform to the same displacement for each mass of the system. However, this property does not depend on us having chosen $x_c$ as one coordinate, but also on our canonical choice of $x_r$ as the other. If you choose $\xi_2=x_2$, say, (or indeed anything that is not of the form $f(x_r)$) and repeat the exercise, you'll see that the matrix $\frac{\partial x_1}{\partial\xi_j}$ has the same problem from the quote.
Sorry for the long post that doesn't deal with rotation and is certainly not an answer to your question, but I do hope it helps with the problems you're having. Particularly, I would personally advise strongly against pursuing the geometrical centre as a quantity of interest - it simply ignores the physical contribution of the masses to the geometry and cannot therefore be right. You only need one global coordinate (per spatial dimension, of course) and this must be the centre of mass. The rest of the game, as you correctly infer, is in the choice of internal coordinates. As to where to go from there, I think it would be useful if you posted some examples of what internal coordinates you've been using in your test examples, so that we can get a better idea of how your maths look like.
