Subject: Summary: rotate molecules along princliple moments of inertia

From: "lala podibaer" <[email-masked]>
Date: Fri, 25 Oct 2002 08:22:15 +0000

Hi to all,

first of all I wanted to thank all the people that answered to my question 
(Alexandre Suman de Araujo, Ryszard Czerminsk, Gareth Jones, James J.P. 
Stewart, Richard Gillilan, Huub van Dam, Jordi Villà i Freixa, and Xiang-Jun 
Lu). It seems like the problem raised is quite widespread and very useful 
ways of handling this problem are around. I programmed the stuff last night 
and it works now.

Here is the summary of the most comprehensive emails I received. Please note 
that if close answers were given I mention only one !

Thanks again,
Mike

---------------
It seems that the method described in the Kabsch paper is what you are
looking for:

Wolfgang Kabsch "A solution for the best rotation
to relate two sets of vectors" Acta Cryst. (1976) A32, 922

Ryszard Czerminski
--------------
Least-square fitting can produce inversion, which is a problem if your
structure has chiral centers.  Fortunately there is an easy fix:

Given the two molecules Y and X centered on the origin, the SVD is defined
as Y'X = USV' and the best fit rotation is defined as VU'. U is the
eigenvector matrix. If the determinant of VU' is -1 then the fitting
includes inversion.  If this is the case find the eigenvector in U that
corresponds to the smallest eigenvalue and change it's sign (multiply by
-1).  With U modified, VU' will now be a rotation matrix.

Unfortunately, I don't have a reference for this handy.

Cheers, Gareth
Gareth Jones PhD
----------------------
The problem you describe arises from the fact that the eigenvectors of the 
rotation matrix are ill-defined by a phase.  In MOPAC, the diagonal terms 
are set positive and the phase is  locked such that the determinant has the 
value +1.  Please consider using the following code, iin which the 3x3 
rotation matrix is "evec":

!
!     MAKE DIAGONAL TERMS OBLIGATE POSITIVE
!
  do i = 1, 3
    if (evec(i, i) < 0.d0) then
      do j = 1, 3
        evec(j, i) = -evec(j, i)
      end do
    end if
  end do
!
!   NOW TO ORIENT THE MOLECULE SO THE CHIRALITY IS PRESERVED
!
  sum = evec(1, 1) * (evec(2, 2)*evec(3, 3)-evec(3, 2)*evec(2, 3)) + evec &
  & (1, 2) * (evec(2, 3)*evec(3, 1)-evec(2, 1)*evec(3, 3)) + evec(1, 3) * &
  & (evec(2, 1)*evec(3, 2)-evec(2, 2)*evec(3, 1))
  if (sum < 0) then
    sum = 1.d0
    do j = 1, 3
      if (evec(j, j) < sum) then
        sum = evec(j, j)
        i = j
      end if
    end do
    do j = 1, 3
      evec(j, i) = -evec(j, i)
    end do
  end if

James J.P. Stewart
---------------------------------------------
This same problem pops up when I have an ensemble of
identical, but randomly-oriented bodies (proteins)
in solution and I wish to transform them all back to
- the same - body coordinate system. You get to the body
system by diagonalizing the inertial matix, but I suspect
the problem is that eigenvalues can be returned by the
numerical routine in any order and with any sign. The
resulting body coordinate sytems then differ by
permutations and sign changes of axes subject to the
condition that the handedness must be preserved. I seem
to remember this has something to do with permutation groups
so that the required sign changes are related to the
parity of the permutation (swap axes and you must change a sign
to preserve handedness). At any rate, is there a way of putting
eigenvectors in some kind of canonical order with appropriate
sign so that the diagonalization routines always land you in the
same body system?

Richard Gillilan
MacCHESS
--------------------------------
The problem is not in the eigenvalues but with the phase factor of the
eigenvectors. If x is an eigenvector of a matrix -x is also an eigenvector 
with
the same eigenvalue. So what you need to do is tweak the sign of your
eigenvectors to remove the mirroring.

If your angles of rotation are small then simply calculate the inproduct of 
the
resulting coordinates with the starting coordinates. These overlaps should 
be
close to +1 if they are close to -1 you have to change a sign somewhere.
Alternatively you apply the transformation to the three unit vectors 
(1,0,0),
(0,1,0) and (0,0,1). Calculate the inproducts of the resulting vectors with 
the
unit vectors. I think you should find the same value for all of the 
inproducts.
Differences in signs mean you have to change a sign.
Of course you can always outer products to check whether your coordinate 
system
after transformation is still a right-hand coordinate system.

I hope this is of some help, although I realize this is not a complete
solution.

Huub van Dam
------------------------------
In a utility called "rotate" included in the program POLYRATE from Don 
Truhlar's group (http://t1.chem.umn.edu/polyrate/), we included a dirty 
trick to fix this problem. The
molecules had the same number of atoms, so you must figure out how to do 
this in a more general case. If the number of atoms is the same, we 
calculated the distance s between
the two sets of coordinates. Once the alignment was ok, we mirrored the 
coordinates according to several planes of simmetry and then recalculated 
the s, getting the alignment according to the minimum s.
It is dirty but it worked.

Jordi Villà i Freixa
-------------------------------
I believe your problem arises because of the ambiguity in selecting the
"direction" of the eigenvectors. Let the first and second eigenvectors be
whatever you have, you need to make sure the third one is the cross product
of the first two, thus ensuring a right-handed coordinate system.

"rotate_mol", a utility in 3DNA 
(http://rutchem.rutgers.edu/~xiangjun/3DNA/),
can do this and a bit more, and it is used to get the most extended views
for every structure in the updated Nucleic Acid Database (NDB)
(http://beta-ndb.rutgers.edu/).

Xiang-Jun
-------------------------------------
You raised a good point and there are indeed some technical details to go
into to get things "right". While I am not so sure about the mathematics and
physics involved, the procedures I used seem to work fine for my purpose:

[1] Pre-process xyz (n-by-3 matrix) to reset its orientation according to
     the reference frame of the first nucleic acid base; or as defined by 
the
     first 3 non-linear atoms. This will eliminate the ambiguity (e.g.,
     upside-down) associated with arbitrary initial orientation.
[2] get the co-variance matrix (real symmetric, 3-by-3).
[3] solve the eigensystem (using e.g. jacobi transformation as in Numerical
     Recipes), sort the eigenvalues, and normalize eigenvectors.
[4] check to make sure the third eigenvector is the cross product of the
     first two, otherwise reverse its sign. There are two more permutations
     which can be used to get three orthogonal views.

Does anyone out there have a more rigorous way?

Xiang-Jun Lu
-------------------------------

_________________________________________________________________
MSN - More useful every day http://www.msn.de