1
2
3
4
5
6 from Numeric import matrixmultiply, transpose, sum, sqrt
7 from LinearAlgebra import singular_value_decomposition, determinant
8
9
11 """
12 SVDSuperimposer finds the best rotation and translation to put
13 two point sets on top of each other (minimizing the RMSD). This is
14 eg. useful to superimpose crystal structures.
15
16 SVD stands for Singular Value Decomposition, which is used to calculate
17 the superposition.
18
19 Reference:
20
21 Matrix computations, 2nd ed. Golub, G. & Van Loan, CF., The Johns
22 Hopkins University Press, Baltimore, 1989
23 """
26
27
28
30 self.reference_coords=None
31 self.coords=None
32 self.transformed_coords=None
33 self.rot=None
34 self.tran=None
35 self.rms=None
36 self.init_rms=None
37
38 - def _rms(self, coords1, coords2):
39 "Return rms deviations between coords1 and coords2."
40 diff=coords1-coords2
41 l=coords1.shape[0]
42 return sqrt(sum(sum(diff*diff))/l)
43
44
45
46 - def set(self, reference_coords, coords):
47 """
48 Set the coordinates to be superimposed.
49 coords will be put on top of reference_coords.
50
51 o reference_coords: an NxDIM array
52 o coords: an NxDIM array
53
54 DIM is the dimension of the points, N is the number
55 of points to be superimposed.
56 """
57
58 self._clear()
59
60 self.reference_coords=reference_coords
61 self.coords=coords
62 n=reference_coords.shape
63 m=coords.shape
64 if n!=m or not(n[1]==m[1]==3):
65 raise Exception, "Coordinate number/dimension mismatch."
66 self.n=n[0]
67
69 "Superimpose the coordinate sets."
70 if self.coords is None or self.reference_coords is None:
71 raise Exception, "No coordinates set."
72 coords=self.coords
73 reference_coords=self.reference_coords
74
75 av1=sum(coords)/self.n
76 av2=sum(reference_coords)/self.n
77 coords=coords-av1
78 reference_coords=reference_coords-av2
79
80 a=matrixmultiply(transpose(coords), reference_coords)
81 u, d, vt=singular_value_decomposition(a)
82 self.rot=transpose(matrixmultiply(transpose(vt), transpose(u)))
83
84 if determinant(self.rot)<0:
85 vt[2]=-vt[2]
86 self.rot=transpose(matrixmultiply(transpose(vt), transpose(u)))
87 self.tran=av2-matrixmultiply(av1, self.rot)
88
98
100 "Right multiplying rotation matrix and translation."
101 if self.rot is None:
102 raise Exception, "Nothing superimposed yet."
103 return self.rot, self.tran
104
106 "Root mean square deviation of untransformed coordinates."
107 if self.coords is None:
108 raise Exception, "No coordinates set yet."
109 if self.init_rms is None:
110 self.init_rms=self._rms(self.coords, self.reference_coords)
111 return self.init_rms
112
114 "Root mean square deviation of superimposed coordinates."
115 if self.rms is None:
116 transformed_coords=self.get_transformed()
117 self.rms=self._rms(transformed_coords, self.reference_coords)
118 return self.rms
119
120
121 if __name__=="__main__":
122
123 from Numeric import *
124
125
126
127 x=array([[51.65, -1.90, 50.07],
128 [50.40, -1.23, 50.65],
129 [50.68, -0.04, 51.54],
130 [50.22, -0.02, 52.85]], Float0)
131
132 y=array([[51.30, -2.99, 46.54],
133 [51.09, -1.88, 47.58],
134 [52.36, -1.20, 48.03],
135 [52.71, -1.18, 49.38]], Float0)
136
137
138 sup=SVDSuperimposer()
139
140
141
142 sup.set(x, y)
143
144
145 sup.run()
146
147
148 rms=sup.get_rms()
149
150
151 rot, tran=sup.get_rotran()
152
153
154 y_on_x1=matrixmultiply(y, rot)+tran
155
156
157 y_on_x2=sup.get_transformed()
158
159 print y_on_x1
160 print
161 print y_on_x2
162 print
163 print "%.2f" % rms
164