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