Author: bugman
Date: Wed Oct 26 10:38:16 2011
New Revision: 14919
URL: http://svn.gna.org/viewcvs/relax?rev=14919&view=rev
Log:
Shifted out the Kabsch superimposition code into its own module.
Added:
1.3/generic_fns/structure/superimpose.py
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=14919&r1=14918&r2=14919&view=diff
==============================================================================
--- 1.3/generic_fns/structure/api_base.py (original)
+++ 1.3/generic_fns/structure/api_base.py Wed Oct 26 10:38:16 2011
@@ -29,9 +29,7 @@
"""
# Python module imports.
-from math import pi
-from numpy import array, diag, dot, float64, outer, sign, transpose, zeros
-from numpy.linalg import det, norm, svd
+from numpy import float64
from os import sep
from re import match
from types import MethodType
@@ -40,7 +38,7 @@
# relax module import.
from data.relax_xml import fill_object_contents, node_value_to_python,
xml_to_object
from float import floatAsByteArray, packBytesAsPyFloat
-from maths_fns.rotation_matrix import R_to_axis_angle
+from generic_fns.structure.superimpose import kabsch
from relax_errors import RelaxError, RelaxFileError,
RelaxFromXMLNotEmptyError, RelaxImplementError
from relax_io import file_root
from relax_warnings import RelaxWarning
@@ -852,63 +850,6 @@
self._rotation_angle = {}
- def _calc_centriod(self, coords):
- """Calculate the centroid of the structure.
-
- @keyword coord: The atomic coordinates.
- @type coord: numpy rank-2, Nx3 array
- @return: The centroid.
- @rtype: numpy rank-1, 3D array
- """
-
- # The sum.
- centroid = coords.sum(0) / coords.shape[0]
-
- # Return.
- return centroid
-
-
- def _calc_rotation(self, coord_from=None, coord_to=None,
centroid_from=None, centroid_to=None):
- """Calculate the rotation via SVD.
-
- @keyword coord_from: The list of atomic coordinates for the
starting structure.
- @type coord_from: numpy rank-2, Nx3 array
- @keyword coord_to: The list of atomic coordinates for the
ending structure.
- @type coord_to: numpy rank-2, Nx3 array
- @keyword centroid_from: The starting centroid.
- @type centroid_from: numpy rank-1, 3D array
- @keyword centroid_to: The ending centroid.
- @type centroid_to: numpy rank-1, 3D array
- @return: The rotation matrix.
- @rtype: numpy rank-2, 3D array
- """
-
- # Initialise the covariance matrix A.
- A = zeros((3, 3), float64)
-
- # Loop over the atoms.
- for i in range(coord_from.shape[0]):
- # The positions shifted to the origin.
- orig_from = coord_from[i] - centroid_from
- orig_to = coord_to[i] - centroid_to
-
- # The outer product.
- A += outer(orig_from, orig_to)
-
- # SVD.
- U, S, V = svd(A)
-
- # The handedness of the covariance matrix.
- d = sign(det(A))
- D = diag([1, 1, d])
-
- # The rotation.
- R = dot(transpose(V), dot(D, transpose(U)))
-
- # Return the rotation.
- return R
-
-
def _calculate(self, model_from=None, model_to=None, coord_from=None,
coord_to=None, centroid=None):
"""Calculate the rotational and translational displacements using
the given coordinate sets.
@@ -926,35 +867,6 @@
@keyword centroid: An alternative position of the centroid,
used for studying pivoted systems.
@type centroid: list of float or numpy rank-1, 3D array
"""
-
- # Calculate the centroids.
- if centroid != None:
- centroid_from = centroid
- centroid_to = centroid
- else:
- centroid_from = self._calc_centriod(coord_from)
- centroid_to = self._calc_centriod(coord_to)
-
- # The translation.
- trans_vect = centroid_to - centroid_from
- trans_dist = norm(trans_vect)
-
- # Calculate the rotation.
- R = self._calc_rotation(coord_from=coord_from, coord_to=coord_to,
centroid_from=centroid_from, centroid_to=centroid_to)
- axis, angle = R_to_axis_angle(R)
-
- # Print out.
- print("\n\nCalculating the rotational and translational
displacements from model %s to model %s.\n" % (model_from, model_to))
- print("Start centroid: [%20.15f, %20.15f, %20.15f]" %
(centroid_from[0], centroid_from[1], centroid_from[2]))
- print("End centroid: [%20.15f, %20.15f, %20.15f]" %
(centroid_to[0], centroid_to[1], centroid_to[2]))
- print("Translation vector: [%20.15f, %20.15f, %20.15f]" %
(trans_vect[0], trans_vect[1], trans_vect[2]))
- print("Translation distance: %.15f" % trans_dist)
- print("Rotation matrix:")
- print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0,
2]))
- print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1,
2]))
- print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2,
2]))
- print("Rotation axis: [%20.15f, %20.15f, %20.15f]" %
(axis[0], axis[1], axis[2]))
- print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0))
# Initialise structures if necessary.
if not self._translation_vector.has_key(model_from):
@@ -968,6 +880,9 @@
if not self._rotation_angle.has_key(model_from):
self._rotation_angle[model_from] = {}
+ # The Kabsch algorithm.
+ trans_vect, trans_dist, R, axis, angle = kabsch(name_from='model
%s'%model_from, name_to='model %s'%model_to, coord_from=coord_from,
coord_to=coord_to, centroid=centroid)
+
# Store the data.
self._translation_vector[model_from][model_to] = trans_vect
self._translation_distance[model_from][model_to] = trans_dist
Added: 1.3/generic_fns/structure/superimpose.py
URL:
http://svn.gna.org/viewcvs/relax/1.3/generic_fns/structure/superimpose.py?rev=14919&view=auto
==============================================================================
--- 1.3/generic_fns/structure/superimpose.py (added)
+++ 1.3/generic_fns/structure/superimpose.py Wed Oct 26 10:38:16 2011
@@ -1,0 +1,143 @@
+###############################################################################
+#
#
+# Copyright (C) 2011 Edward d'Auvergne
#
+#
#
+# This file is part of the program relax.
#
+#
#
+# relax is free software; you can redistribute it and/or modify
#
+# it under the terms of the GNU General Public License as published by
#
+# the Free Software Foundation; either version 2 of the License, or
#
+# (at your option) any later version.
#
+#
#
+# relax is distributed in the hope that it will be useful,
#
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
#
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#
+# GNU General Public License for more details.
#
+#
#
+# You should have received a copy of the GNU General Public License
#
+# along with relax; if not, write to the Free Software
#
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
+#
#
+###############################################################################
+
+# Module docstring.
+"""Module for handling all types of structural superimpositions."""
+
+# Python module imports.
+from math import pi
+from numpy import diag, dot, float64, outer, sign, transpose, zeros
+from numpy.linalg import det, norm, svd
+
+# relax module import.
+from maths_fns.rotation_matrix import R_to_axis_angle
+
+
+def find_centroid(coords):
+ """Calculate the centroid of the structural coordinates.
+
+ @keyword coord: The atomic coordinates.
+ @type coord: numpy rank-2, Nx3 array
+ @return: The centroid.
+ @rtype: numpy rank-1, 3D array
+ """
+
+ # The sum.
+ centroid = coords.sum(0) / coords.shape[0]
+
+ # Return.
+ return centroid
+
+
+def kabsch(name_from=None, name_to=None, coord_from=None, coord_to=None,
centroid=None, verbosity=1):
+ """Calculate the rotational and translational displacements between the
two given coordinate sets.
+
+ This uses the Kabsch algorithm
(http://en.wikipedia.org/wiki/Kabsch_algorithm).
+
+
+ @keyword name_from: The name of the starting structure, used for the
print outs.
+ @type name_from: str
+ @keyword name_to: The name of the ending structure, used for the
print outs.
+ @type name_to: str
+ @keyword coord_from: The list of atomic coordinates for the starting
structure.
+ @type coord_from: numpy rank-2, Nx3 array
+ @keyword coord_to: The list of atomic coordinates for the ending
structure.
+ @type coord_to: numpy rank-2, Nx3 array
+ @keyword centroid: An alternative position of the centroid, used
for studying pivoted systems.
+ @type centroid: list of float or numpy rank-1, 3D array
+ @return: The translation vector T, translation distance
d, rotation matrix R, rotation axis r, and rotation angle theta.
+ @rtype: numpy rank-1 3D array, float, numpy rank-2 3D
array, numpy rank-1 3D array, float
+ """
+
+ # Calculate the centroids.
+ if centroid != None:
+ centroid_from = centroid
+ centroid_to = centroid
+ else:
+ centroid_from = find_centroid(coord_from)
+ centroid_to = find_centroid(coord_to)
+
+ # The translation.
+ trans_vect = centroid_to - centroid_from
+ trans_dist = norm(trans_vect)
+
+ # Calculate the rotation.
+ R = kabsch_rotation(coord_from=coord_from, coord_to=coord_to,
centroid_from=centroid_from, centroid_to=centroid_to)
+ axis, angle = R_to_axis_angle(R)
+
+ # Print out.
+ if verbosity >= 1:
+ print("\n\nCalculating the rotational and translational
displacements from model %s to model %s.\n" % (name_from, name_to))
+ print("Start centroid: [%20.15f, %20.15f, %20.15f]" %
(centroid_from[0], centroid_from[1], centroid_from[2]))
+ print("End centroid: [%20.15f, %20.15f, %20.15f]" %
(centroid_to[0], centroid_to[1], centroid_to[2]))
+ print("Translation vector: [%20.15f, %20.15f, %20.15f]" %
(trans_vect[0], trans_vect[1], trans_vect[2]))
+ print("Translation distance: %.15f" % trans_dist)
+ print("Rotation matrix:")
+ print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0,
2]))
+ print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1,
2]))
+ print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2,
2]))
+ print("Rotation axis: [%20.15f, %20.15f, %20.15f]" %
(axis[0], axis[1], axis[2]))
+ print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0))
+
+ # Return the data.
+ return trans_vect, trans_dist, R, axis, angle
+
+
+def kabsch_rotation(coord_from=None, coord_to=None, centroid_from=None,
centroid_to=None):
+ """Calculate the rotation via SVD.
+
+ @keyword coord_from: The list of atomic coordinates for the starting
structure.
+ @type coord_from: numpy rank-2, Nx3 array
+ @keyword coord_to: The list of atomic coordinates for the ending
structure.
+ @type coord_to: numpy rank-2, Nx3 array
+ @keyword centroid_from: The starting centroid.
+ @type centroid_from: numpy rank-1, 3D array
+ @keyword centroid_to: The ending centroid.
+ @type centroid_to: numpy rank-1, 3D array
+ @return: The rotation matrix.
+ @rtype: numpy rank-2, 3D array
+ """
+
+ # Initialise the covariance matrix A.
+ A = zeros((3, 3), float64)
+
+ # Loop over the atoms.
+ for i in range(coord_from.shape[0]):
+ # The positions shifted to the origin.
+ orig_from = coord_from[i] - centroid_from
+ orig_to = coord_to[i] - centroid_to
+
+ # The outer product.
+ A += outer(orig_from, orig_to)
+
+ # SVD.
+ U, S, V = svd(A)
+
+ # The handedness of the covariance matrix.
+ d = sign(det(A))
+ D = diag([1, 1, d])
+
+ # The rotation.
+ R = dot(transpose(V), dot(D, transpose(U)))
+
+ # Return the rotation.
+ return R
r14919 - in /1.3/generic_fns/structure: api_base.py superimpose.py