BR's crystallographic computing tutorials


Superposition of two molecules

The problem of superpositioning two (or more) molecules onto each other is a common task and one would assume it is also very simple. In fact, there exist several techniques to solve this problem, but there are a few details to consider which can make coding of the problem subtle. In the following I will discuss a FORTRAN program which superimposes two pdb files using the elegant quaternion method as described by S.K.Kearsley [1].

Procedure

The superposition problem can be split up into two parts: A rotation around the geometric center of the molecule to bring it into the proper orientation, and a translation that superimposes the centers of the probe and target molecule(s).

The geometric center of the molecule can be found readily by averaging the x,y, and z coordinates of the n atoms k, respectively.

This center is not exactly the center of mass (we did not care what kind the atoms are, i.e., we did not calculate the moment). But a center defined this way will be the origin of the principal axes setting up the distance ellipsoid (similar to e.g., an anisotropic thermal parameter ellipsoid).

The (orthogonal) principal axes [2] can be obtained quite easily by orthogonalization of the matrix A (i.e., the matrix filled by the sum of the matrices of the metric tensor x ) :

where

by a Jacobi transformation. The Jacobi transformation (or orthogonalization) is an iterative application of rotations to a matrix until all the off-diagonal values are zero at machine precision [3]. If the matrix A is symmetric (which A of course is) then a diagonal matrix D exists so that

where D has the eigenvalues in its diagonal :

and P the normalized (orthogonal) eigenvectors of the principal axes in its columns. This orthogonalization is in fact equivalent to the solution of the eigenvalue problem [4]

.

The feature of A is that it is real and symmetric means that a solution with real eigenvectors and eigenvalues must exist. This is good, but it is also the source of the mmm symmetry of the ellipsoid, a nasty property as we'll see shortly.

Once we have determined the center of each molecule, we can shift them to a common origin, usually (0,0,0). We are now faced with the problem to find the rotation that gives the best superposition of the molecules. What defines 'best' superposition? One might be tempted to just rotate the principal axes to an overlap. This could in principle work, and the rotation matrix to overlap the axes can be easily calculated by a Gauss-Jordan elimination [3] of three 3x3 linear systems. This does give the correct answer in cases where the molecules are pretty much pre-aligned, and the molecules can also have a different number of atoms. But as a result of the symmetry of the principal axes ellipsoid we have lost information about the orientation of the axes: we need to try 6 rotations (3 pairs of 180 deg rotated vectors) to find the actual solution. To determine the best one, we are back to a calculation of the pairwise distance r.m.s. (root mean square) between atoms which we hoped to avoid in the first place.

So do we have to solve this as a least squares problem of r.m.s. distances with all the problems associated (see, e.g. [1] for a discussion and further reading)? No, not really. The above idea using the principal axes rotation was quite elegant, but it needs to be improved. We need to find a way to resolve the degeneracy of the principal axes transformation. This can actually be done with an extension of the 3x3 problem of the metric tensor matrix to a 4x4 matrix which is constructed from a non-cyclic permutation of the pairwise differences between atom positions. The combination of elements follows the 4 constructors of the Quaternion group H, a small, non-cyclic subgroup of GL2(C) [5].

We solve our linear quaternion eigenvalue system again numerically by Jacobi transformation of the matrix Q

, where again

How the elements q are formed can be found in [1]. The smallest eigenvalue is the s.r.s.(sum of residuals squared) for the best rotation and the associated eigenvector contains element from which the 3x3 rotation matrix t for the best superposition is constructed. The largest EV contains the worst rotation, the associated eigenvalue the worst s.r.s. The Quaternion method is a very elegant way to solve the rotation problem and the r.m.s.d. least squares minimization in one shot!

The application of the rotation and the back-translation of the molecule are trivial final steps.


[1] S.K.Kearsley, On the orthogonal transformation used for structural comparisons, Acta Cryst. A45, 208 (1989)
[2] D.A.Danielson, Vectors and Tensors in Engineering and Physics, Addison-Wesley (1992), D.A.Sands, Vectors and Tensors in Crystallography, Dover Publishing (1982)
[3] W.H.Press, S.A.Teukolsky, W.T.Vettering and B.P. Flannery, Numerical Recipes, 2nd edition, Cambridge University Press (1992)
[4] S. Lipschutz, Linear Algebra, 2nd edition, Schaum Outline Series, McGraw-Hill (1991)
[5] M.Artin, Algebra, Prentice-Hall (1991)
[6] J.Hart, Quaternion demonstrator (SGI only), Stanford University


Back to Introduction
This World Wide Web site conceived and maintained by Bernhard Rupp
Last revised Dezember 27, 2009 01:40