The next function we describe is the clever mathematical part of the whole coordinate
alignment operation. It takes two arrays of coordinates, and an optional array of weights,
then calculates the rotation that best superimposes equivalent pairs of coordinates
(corresponding atoms from each structure). The underlying mathematical theory is
somewhat involved, and we will not discuss this in great detail, but by reading through the
Python code it should be clear as to what the overall effect of the commands is.
The function is defined, and if we do not pass in an array of weights a new one is
created using ones(), to give a series of 1’s the same length as the coordinates (one for
each atom).
def alignCoords(coordsA, coordsB, weights=None):
n = len(coordsA)
if weights is None:
weights = ones(n)
The alignment routine is basically a minimisation of the differences in positions
between corresponding coordinates, which we can achieve by finding an optimum rotation
transformation. We define an initial matrix as the weighted dot product (matrix
multiplication) of the one coordinate array and the transpose of the other. We need to make
the transpose (switch rows for columns) to align the long, atom axes of the arrays so we
end up with a 3×3 matrix. You can imagine each coordinate array as a transformation
between atom number and spatial position. Accordingly, by applying one array to the
other we get the transformation of positions, through a common set of atoms, to new
positions, thus defining a spatial transformation.
rMat = dot(coordsB.transpose()*weights, coordsA)
The mathematical magic that performs the actual minimisation is an operation called
singular value decomposition or SVD for short. Luckily the linear algebra module of
NumPy (linalg) has a convenient svd() function that will do all the hard work for us.
Effectively this takes the 3×3 transformation matrix defined above and splits (factorises) it
into three components, the combination of which would perform the same transformation.
These three components are as follows: a rotation matrix, an array of linear scaling factors
and an opposing rotation matrix.
15
The detailed explanation of why this is done is
probably very confusing if you don’t already have a good understanding of what are
known as eigenvectors. However, it is sufficient to say that the matrices generated
represent very special directions, which allows us to get directly at the rotation that we
need to apply to align the coordinates.
rMat1, scales, rMat2 = linalg.svd(rMat)
Before using the extracted rotation matrices to align the coordinates, we must first
check whether the SVD has given transformations that would cause a mirror image; the
decomposition does not distinguish between normal and reflected solutions. To address
this problem we calculate what is known as the determinant of the matrices: a single
number that represents the overall scaling factor for the transformation. Helpfully, the
linalg.det() function easily calculates the determinant, which will be −1 or +1 for our
matrices. Given that the SVD operation has factored out a change in size, all that remains
is a change in sign. The two determinants are multiplied, so that if they have opposite sign
the result is −1.
sign = linalg.det(rMat1) * linalg.det(rMat2)
We check whether the sign is negative, and if so we know to remove the reflection. This
is done by flipping the sign of the numbers in the last column in the matrix.
if sign < 0:
rMat1[:,2] *= -1
Then the optimised rotation matrix, which will allow us to transform the coordinates to
do the superimposition, is calculated by multiplying the two matrices obtained by the
linalg.svd decomposition. Effectively this is reconstructing most of the original rMat
transformation, except that the scaling, extracted by the SVD factorisation, has been
eliminated so that only a pure rotation remains.
rotation = dot(rMat1, rMat2)
We use the final rotation matrix to transform one set of coordinates, so that it is
repositioned closer to the other set. Naturally we only move one coordinate array, and the
other is left alone. The dot() function is used to perform the matrix multiplication on the
arrays: apply the coordinate transformation. Lastly, we return the calculated rotation
matrix and the updated coordinates:
coordsB = dot(coordsB, rotation)
return rotation, coordsB
Do'stlaringiz bilan baham: 
