# 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?

-
Is the system invariant or not under translations and rotations? i.e.: do you have an invariant system and are hoping to include a description of translation and rotation? or do you have interactions that break the translational and rotational invariance? if so, what coordinates do these interactions depend on? –  Emilio Pisanty Jul 9 '12 at 14:09
@EmilioPisanty In general, the system is not invariant under translations and rotations. The dependence is not simple: there's a large number of fixed external point charges ("fixed" meaning they are not affected by rotation or translation). Note that I'm not asking for a method or algorithm to optimize the geometry, but for a way to add translation and rotation degrees of freedom of the whole system to a set of redundant internal coordinates. –  Jellby Jul 9 '12 at 14:46
If you are looking for a way to describe an orientation of an object, Euler angles might be what you are looking for. –  Alexey Bobrick Jul 9 '12 at 23:50
What's wrong with just using cartesian coordinates for everything? You can add SO(3) coordinates for the system, but you can also just write down the force laws in cartesian coordinates, or using SO(3) coordinates plus translations for the internal stuff. –  Ron Maimon Jul 10 '12 at 20:06
@RonMaimon In chemical systems it is often more efficient (fewer iterations) to use internal coordinates. Internal coordinates also make it easier to add constraints to distances and angles, and are more closely related to the "meaning" of the process. But of course, using Cartesian coordinates is always a possibility. –  Jellby Jul 11 '12 at 7:59

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.

-
I'd need two objects (not aligned with the center of rotation) to define the orientation, I guess. Otherwise the orientation around the axis determined by the first object is undefined, right? Now, when it comes to defining all the derivatives $B_{ij}$, does it mean all derivatives with respect to objects other that 1 are zero? Or are all as in the equations in my answer? –  Jellby Jul 16 '12 at 8:59
@Jellby: I assumed the "objects" are irregular, so that they each have an orientation as rigid bodies. Are these rigid bodies or point particles? If they are point particles, ignore this answer, but from your description, it seemed you had orientation matrices for pairs, so they were rigid bodies, like amino acids. The derivatives "B_ij" for the rotation matrices are the angular velocities, which are given in antisymmetric matrix form by $R^{-1}\dot{R}$ where dot is a time derivative. –  Ron Maimon Jul 16 '12 at 9:15
No, my "objects" are just point particles (nuclei), and there is no velocity either, as there's no time. The derivatives I need are of internal (or global) coordinates with respect to Cartesian coordinates of the nuclei, or vice versa, which are used to convert between the two kinds of coordinates. –  Jellby Jul 16 '12 at 9:37
@Jellby: Ok, then you need to fix four arbitrary point particles to define the frame. You should have said this, because if you are doing point particles, it is absolutely 100% wrong to use any sort of relative coordinates, you should use cartesian coordinates period, and stop wasting time. –  Ron Maimon Jul 16 '12 at 17:27
I'm sorry that's what you think. But as I said, in the comments to my question, Cartesian coordinates, while possible, are generally not the most efficient. It is true this is not a life-or-death question, I can have the job done by using Cartesian coordinates, but since I already have a program that works very nicely with internal coordinates when there is translational and rotational invariance, I just needed to add these degrees of freedom. –  Jellby Jul 16 '12 at 17:29

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.

-
I hadn't thought about how the presence of other internal coordinates would affect the inverse conversion matrix. I guess I thought that, if I remove all degrees of freedom except center of mass translation, the system will move as a rigid object, but now I believe that's not the case, rather internal degrees of freedom will change in undetermined ways. I've added a footnote to my answer. As for the physics of the problem, since I'm optimizing a geometry and the (potential) energy is independent of the masses, I don't think there's any problem with using the geometric center. –  Jellby Jul 16 '12 at 8:53
And about examples of what internal coordinates, I'm using in general redundant internal coordinates, that is (a subset af all possible) distances between two points, angles between three points and dihedral formed by four points. In methanol, for instance, I'd have 3 CH distances, 1 CO, 1 OH, 3 HCH angles, 3 HCO, 1 COH, 3 HCOH dihedral. That's a total of 15 internal coordinates, which is more than the 3*6-6=12 needed (that's why they are redundant), but still don't contain the translational and rotational degrees of freedom. –  Jellby Jul 16 '12 at 10:40
For methanol, I'd use the Euler angles of the CO bond. You should be careful because at equilibrium the two smaller moments of inertia are equal and, as Ron points out, this can bring trouble. –  Emilio Pisanty Jul 16 '12 at 12:50
But I'm after a solution that would work for any other system, even linear. I guess I'll have to choose some specific atoms to define the orientation (or rather think of an algorithm to choose them automatically). –  Jellby Jul 16 '12 at 13:09
Linear molecules are a problem because two of their moments of inertia are basically zero. (The electronic angular momentum is already accounted for and nuclear rotation only contributes to hyperfine structure.) Choosing the orientation of any given bond (or possibly two adjacent bonds) will most likely work. Having a universal algorithm for choosing which bond to use sounds unlikely. –  Emilio Pisanty Jul 16 '12 at 13:50

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.

-