1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module for handling all types of structural superimpositions."""
24
25
26 from copy import deepcopy
27 from math import pi
28 from numpy import diag, dot, eye, float64, outer, sign, transpose, zeros
29 from numpy.linalg import det, norm, svd
30
31
32 from lib.structure.mass import centre_of_mass
33 from lib.structure.statistics import calc_mean_structure
34 from lib.geometry.rotations import R_to_axis_angle, R_to_euler_zyz
35
36
38 """Calculate the centroid of the structural coordinates.
39
40 @keyword coord: The atomic coordinates.
41 @type coord: numpy rank-2, Nx3 array
42 @return: The centroid.
43 @rtype: numpy rank-1, 3D array
44 """
45
46
47 centroid = coords.sum(0) / coords.shape[0]
48
49
50 return centroid
51
52
53 -def fit_to_first(models=None, coord=None, centre_type="centroid", elements=None, centroid=None):
54 """Superimpose a set of structural models using the fit to first algorithm.
55
56 @keyword models: The list of models to superimpose.
57 @type models: list of int
58 @keyword coord: The list of coordinates of all models to superimpose. The first index is the models, the second is the atomic positions, and the third is the xyz coordinates.
59 @type coord: list of numpy rank-2, Nx3 arrays
60 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions.
61 @type centroid: list of float or numpy rank-1, 3D array
62 @keyword centre_type: The type of centre to superimpose over. This can either be the standard centroid superimposition or the CoM could be used instead.
63 @type centre_type: str
64 @keyword elements: The list of elements corresponding to the atoms.
65 @type elements: list of str
66 @return: The lists of translation vectors, rotation matrices, and rotation pivots.
67 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays
68 """
69
70
71 print("\nSuperimposition of structural models %s using the 'fit to first' algorithm." % models)
72
73
74 T_list = [zeros(3, float64)]
75 R_list = [eye(3, dtype=float64)]
76 pivot_list = [zeros(3, float64)]
77
78
79 for i in range(1, len(models)):
80
81 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='model %s'%models[i], coord_from=coord[i], coord_to=coord[0], centre_type=centre_type, elements=elements, centroid=centroid)
82
83
84 T_list.append(trans_vect)
85 R_list.append(R)
86 pivot_list.append(pivot)
87
88
89 return T_list, R_list, pivot_list
90
91
92 -def fit_to_mean(models=None, coord=None, centre_type="centroid", elements=None, centroid=None, verbosity=1):
93 """Superimpose a set of structural models using the fit to first algorithm.
94
95 @keyword models: The list of models to superimpose.
96 @type models: list of int
97 @keyword coord: The list of coordinates of all models to superimpose. The first index is the models, the second is the atomic positions, and the third is the xyz coordinates.
98 @type coord: list of numpy rank-2, Nx3 arrays
99 @keyword centre_type: The type of centre to superimpose over. This can either be the standard centroid superimposition or the CoM could be used instead.
100 @type centre_type: str
101 @keyword elements: The list of elements corresponding to the atoms.
102 @type elements: list of str
103 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions.
104 @type centroid: list of float or numpy rank-1, 3D array
105 @keyword verbosity: The amount of information to print out. If 0, nothing will be printed.
106 @type verbosity: int
107 @return: The lists of translation vectors, rotation matrices, and rotation pivots.
108 @rtype: list of numpy rank-1 3D arrays, list of numpy rank-2 3D arrays, list of numpy rank-1 3D arrays
109 """
110
111
112 if verbosity:
113 print("\nSuperimposition of structural models %s using the 'fit to mean' algorithm." % models)
114
115
116 orig_coord = deepcopy(coord)
117
118
119 T_list = []
120 R_list = []
121 pivot_list = []
122
123
124 N = len(coord[0])
125 mean = zeros((N, 3), float64)
126
127
128 converged = False
129 iter = 0
130 while not converged:
131
132 if verbosity:
133 print("\nIteration %i of the algorithm." % iter)
134 print("%-10s%-25s%-25s" % ("Model", "Translation (Angstrom)", "Rotation (deg)"))
135
136
137 calc_mean_structure(coord, mean)
138
139
140 converged = True
141 for i in range(len(models)):
142
143 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[0], name_to='mean', coord_from=coord[i], coord_to=mean, centre_type=centre_type, elements=elements, centroid=centroid, verbosity=0)
144
145
146 if verbosity:
147 print("%-10i%25.3g%25.3g" % (i, trans_dist, (angle / 2.0 / pi * 360.0)))
148
149
150 for j in range(N):
151
152 coord[i][j] += trans_vect
153
154
155 coord[i][j] -= pivot
156
157
158 coord[i][j] = dot(R, coord[i][j])
159
160
161 coord[i][j] += pivot
162
163
164 if trans_dist > 1e-10 or angle > 1e-10:
165 converged = False
166
167
168 iter += 1
169
170
171 for i in range(len(models)):
172
173 trans_vect, trans_dist, R, axis, angle, pivot = kabsch(name_from='model %s'%models[i], name_to='the mean structure', coord_from=orig_coord[i], coord_to=mean, centre_type=centre_type, elements=elements, centroid=centroid, verbosity=0)
174
175
176 T_list.append(trans_vect)
177 R_list.append(R)
178 pivot_list.append(pivot)
179
180
181 return T_list, R_list, pivot_list
182
183
184 -def kabsch(name_from=None, name_to=None, coord_from=None, coord_to=None, centre_type="centroid", elements=None, centroid=None, verbosity=1):
185 """Calculate the rotational and translational displacements between the two given coordinate sets.
186
187 This uses the U{Kabsch algorithm<http://en.wikipedia.org/wiki/Kabsch_algorithm>}.
188
189
190 @keyword name_from: The name of the starting structure, used for the printouts.
191 @type name_from: str
192 @keyword name_to: The name of the ending structure, used for the printouts.
193 @type name_to: str
194 @keyword coord_from: The list of atomic coordinates for the starting structure.
195 @type coord_from: numpy rank-2, Nx3 array
196 @keyword coord_to: The list of atomic coordinates for the ending structure.
197 @type coord_to: numpy rank-2, Nx3 array
198 @keyword centre_type: The type of centre to superimpose over. This can either be the standard centroid superimposition or the CoM could be used instead.
199 @type centre_type: str
200 @keyword elements: The list of elements corresponding to the atoms.
201 @type elements: list of str
202 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems.
203 @type centroid: list of float or numpy rank-1, 3D array
204 @return: The translation vector T, translation distance d, rotation matrix R, rotation axis r, rotation angle theta, and the rotational pivot defined as the centroid of the ending structure.
205 @rtype: numpy rank-1 3D array, float, numpy rank-2 3D array, numpy rank-1 3D array, float, numpy rank-1 3D array
206 """
207
208
209 if centroid != None:
210 centroid_from = centroid
211 centroid_to = centroid
212 elif centre_type == 'centroid':
213 centroid_from = find_centroid(coord_from)
214 centroid_to = find_centroid(coord_to)
215 else:
216 centroid_from, mass_from = centre_of_mass(pos=coord_from, elements=elements)
217 centroid_to, mass_to = centre_of_mass(pos=coord_to, elements=elements)
218
219
220 trans_vect = centroid_to - centroid_from
221 trans_dist = norm(trans_vect)
222
223
224 R = kabsch_rotation(coord_from=coord_from, coord_to=coord_to, centroid_from=centroid_from, centroid_to=centroid_to)
225 axis, angle = R_to_axis_angle(R)
226 a, b, g = R_to_euler_zyz(R)
227
228
229 if verbosity >= 1:
230 print("\n\nCalculating the rotational and translational displacements from %s to %s using the Kabsch algorithm.\n" % (name_from, name_to))
231 if centre_type == 'centroid':
232 print("Start centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_from[0], centroid_from[1], centroid_from[2]))
233 print("End centroid: [%20.15f, %20.15f, %20.15f]" % (centroid_to[0], centroid_to[1], centroid_to[2]))
234 else:
235 print("Start CoM: [%20.15f, %20.15f, %20.15f]" % (centroid_from[0], centroid_from[1], centroid_from[2]))
236 print("End CoM: [%20.15f, %20.15f, %20.15f]" % (centroid_to[0], centroid_to[1], centroid_to[2]))
237 print("Translation vector: [%20.15f, %20.15f, %20.15f]" % (trans_vect[0], trans_vect[1], trans_vect[2]))
238 print("Translation distance: %.15f" % trans_dist)
239 print("Rotation matrix:")
240 print(" [[%20.15f, %20.15f, %20.15f]" % (R[0, 0], R[0, 1], R[0, 2]))
241 print(" [%20.15f, %20.15f, %20.15f]" % (R[1, 0], R[1, 1], R[1, 2]))
242 print(" [%20.15f, %20.15f, %20.15f]]" % (R[2, 0], R[2, 1], R[2, 2]))
243 print("Rotation axis: [%20.15f, %20.15f, %20.15f]" % (axis[0], axis[1], axis[2]))
244 print("Rotation euler angles: [%20.15f, %20.15f, %20.15f]" % (a, b, g))
245 print("Rotation angle (deg): %.15f" % (angle / 2.0 / pi * 360.0))
246
247
248 return trans_vect, trans_dist, R, axis, angle, centroid_to
249
250
251 -def kabsch_rotation(coord_from=None, coord_to=None, centroid_from=None, centroid_to=None):
252 """Calculate the rotation via SVD.
253
254 @keyword coord_from: The list of atomic coordinates for the starting structure.
255 @type coord_from: numpy rank-2, Nx3 array
256 @keyword coord_to: The list of atomic coordinates for the ending structure.
257 @type coord_to: numpy rank-2, Nx3 array
258 @keyword centroid_from: The starting centroid.
259 @type centroid_from: numpy rank-1, 3D array
260 @keyword centroid_to: The ending centroid.
261 @type centroid_to: numpy rank-1, 3D array
262 @return: The rotation matrix.
263 @rtype: numpy rank-2, 3D array
264 """
265
266
267 A = zeros((3, 3), float64)
268
269
270 for i in range(coord_from.shape[0]):
271
272 orig_from = coord_from[i] - centroid_from
273 orig_to = coord_to[i] - centroid_to
274
275
276 A += outer(orig_from, orig_to)
277
278
279 U, S, V = svd(A)
280
281
282 d = sign(det(A))
283 D = diag([1, 1, d])
284
285
286 R = dot(transpose(V), dot(D, transpose(U)))
287
288
289 return R
290