r14916 - /1.3/generic_fns/structure/api_base.py -- October 25, 2011 - 16:41

Author: bugman
Date: Tue Oct 25 16:41:43 2011
New Revision: 14916

URL: http://svn.gna.org/viewcvs/relax?rev=14916&view=rev
Log:
Fixes for the structure.displacement user function.

The rotations are now correctly calculated - the problem was that the numpy 
SVD is not transposing
the last matrix v.  The technique is now labelled as the Kabsch algorithm.


Modified:
    1.3/generic_fns/structure/api_base.py

Modified: 1.3/generic_fns/structure/api_base.py
URL: 
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/api_base.py?rev=14916&r1=14915&r2=14916&view=diff
==============================================================================
--- 1.3/generic_fns/structure/api_base.py (original)
+++ 1.3/generic_fns/structure/api_base.py Tue Oct 25 16:41:43 2011
@@ -30,8 +30,8 @@
 
 # Python module imports.
 from math import pi
-from numpy import array, dot, float64, outer, transpose, zeros
-from numpy.linalg import norm, svd
+from numpy import array, diag, dot, float64, outer, sign, transpose, zeros
+from numpy.linalg import det, norm, svd
 from os import sep
 from re import match
 from types import MethodType
@@ -883,8 +883,8 @@
         @rtype:                 numpy rank-2, 3D array
         """
 
-        # Initialise the H matrix.
-        H = zeros((3, 3), float64)
+        # Initialise the covariance matrix A.
+        A = zeros((3, 3), float64)
 
         # Loop over the atoms.
         for i in range(coord_from.shape[0]):
@@ -893,13 +893,17 @@
             orig_to = coord_to[i] - centroid_to
 
             # The outer product.
-            H += outer(orig_from, orig_to)
+            A += outer(orig_from, orig_to)
 
         # SVD.
-        U, S, V = svd(H)
+        U, S, V = svd(A)
+
+        # The handedness of the covariance matrix.
+        d = sign(det(A))
+        D = diag([1, 1, d])
 
         # The rotation.
-        R = dot(V, transpose(U))
+        R = dot(transpose(V), dot(D, transpose(U)))
 
         # Return the rotation.
         return R
@@ -907,6 +911,9 @@
 
     def _calculate(self, model_from=None, model_to=None, coord_from=None, 
coord_to=None):
         """Calculate the rotational and translational displacements using 
the given coordinate sets.
+
+        This uses the Kabsch algorithm 
(http://en.wikipedia.org/wiki/Kabsch_algorithm).
+
 
         @keyword model_from:    The model number of the starting structure.
         @type model_from:       int