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 all of the frame order specific user functions."""
24
25
26 from copy import deepcopy
27 from math import pi
28 from numpy import array, dot, eye, float64, transpose, zeros
29 import sys
30
31
32 from lib.arg_check import is_float_array
33 from lib.errors import RelaxError
34 from lib.frame_order.rotor_axis import create_rotor_axis_alpha, create_rotor_axis_euler, create_rotor_axis_spherical
35 from lib.geometry.rotations import euler_to_R_zyz, two_vect_to_R
36 from lib.io import open_write_file
37 from lib.order import order_parameters
38 from lib.structure.cones import Iso_cone, Pseudo_elliptic
39 from lib.structure.geometric import generate_vector_residues
40 from lib.structure.internal.object import Internal
41 from lib.structure.represent.cone import cone
42 from lib.structure.represent.rotor import rotor_pdb
43 from lib.text.sectioning import subsection
44 from pipe_control import pipes
45 from pipe_control.structure.mass import pipe_centre_of_mass
46 from specific_analyses.frame_order.data import domain_moving, translation_fixed
47 from specific_analyses.frame_order.parameters import update_model
48
49
51 """Set up the mechanics of the average domain position.
52
53 @keyword pivot: What to use as the motional pivot. This can be 'com' for the centre of mass of the moving domain, or 'motional' to link the pivot of the motion to the rotation of the average domain position.
54 @type pivot: str
55 @keyword translation: If True, translation to the average domain position will be allowed. If False, then translation will not occur.
56 @type translation: bool
57 """
58
59
60 pipes.test()
61
62
63 if pivot not in ['com', 'motional']:
64 raise RelaxError("The pivot for the rotation to the average domain position must be either 'com' or 'motional'.")
65
66
67 cdp.ave_pos_pivot = pivot
68 cdp.ave_pos_translation = translation
69
70
72 """Set the number of integration points to use in the quasi-random Sobol' sequence.
73
74 @keyword num: The number of integration points.
75 @type num: int
76 """
77
78
79 pipes.test()
80
81
82 cdp.num_int_pts = num
83
84
86 """Create a PDB file of the molecule with the moving domains shifted to the average position.
87
88 @keyword file: The name of the file for the average molecule structure.
89 @type file: str
90 @keyword dir: The name of the directory to place the PDB file into.
91 @type dir: str
92 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten.
93 @type force: bool
94 """
95
96
97 subsection(file=sys.stdout, text="Creating a PDB file with the moving domains shifted to the average position.")
98
99
100 structure = deepcopy(cdp.structure)
101
102
103 R = zeros((3, 3), float64)
104 if hasattr(cdp, 'ave_pos_alpha'):
105 euler_to_R_zyz(cdp.ave_pos_alpha, cdp.ave_pos_beta, cdp.ave_pos_gamma, R)
106 else:
107 euler_to_R_zyz(0.0, cdp.ave_pos_beta, cdp.ave_pos_gamma, R)
108 if cdp.ave_pos_pivot == 'com':
109 origin = pipe_centre_of_mass(atom_id=domain_moving(), verbosity=0)
110 else:
111 origin = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z])
112 structure.rotate(R=R, origin=origin, atom_id=domain_moving())
113
114
115 if not translation_fixed():
116 structure.translate(T=[cdp.ave_pos_x, cdp.ave_pos_y, cdp.ave_pos_z], atom_id=domain_moving())
117
118
119 file = open_write_file(file_name=file, dir=dir, force=force)
120 structure.write_pdb(file=file)
121 file.close()
122
123
125 """Create a PDB file of a distribution of positions coving the full dynamics of the moving domain.
126
127 @keyword file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model.
128 @type file: str
129 @keyword dir: The name of the directory to place the PDB file into.
130 @type dir: str
131 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten.
132 @type force: bool
133 """
134
135
136 subsection(file=sys.stdout, text="Creating a PDB file of a distribution of positions coving the full dynamics of the moving domain.")
137
138
139 -def pdb_geometric_rep(file=None, dir=None, size=30.0, inc=36, force=False, neg_cone=True):
140 """Create a PDB file containing a geometric object representing the frame order dynamics.
141
142 @keyword file: The name of the file of the PDB representation of the frame order dynamics to create.
143 @type file: str
144 @keyword dir: The name of the directory to place the PDB file into.
145 @type dir: str
146 @keyword size: The size of the geometric object in Angstroms.
147 @type size: float
148 @keyword inc: The number of increments for the filling of the cone objects.
149 @type inc: int
150 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten.
151 @type force: bool
152 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation. This is ignored for the rotor models.
153 @type neg_cone: bool
154 """
155
156
157 subsection(file=sys.stdout, text="Creating a PDB file containing a geometric object representing the frame order dynamics.")
158
159
160 sim = False
161 num_sim = 0
162 if hasattr(cdp, 'sim_number'):
163 sim = True
164 num_sim = cdp.sim_number
165
166
167 inv_mat = -eye(3)
168
169
170 structure = Internal()
171
172
173 if cdp.model in ['rotor', 'free rotor']:
174 neg_cone = False
175 model = structure.add_model(model=1)
176 model_num = 1
177 if neg_cone:
178 model_neg = structure.add_model(model=2)
179 model_num = 2
180
181
182 pivot = array([cdp.pivot_x, cdp.pivot_y, cdp.pivot_z])
183
184
185 if cdp.model in ['rotor', 'free rotor', 'iso cone', 'iso cone, free rotor', 'pseudo-ellipse']:
186
187 if cdp.model in ['free rotor', 'iso cone, free rotor']:
188 rotor_angle = pi
189 else:
190 rotor_angle = cdp.cone_sigma_max
191
192
193 com = pivot
194 if cdp.model in ['rotor', 'free rotor']:
195 com = pipe_centre_of_mass(verbosity=0)
196
197
198 if cdp.model in ['rotor']:
199 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com)
200 elif cdp.model in ['free rotor', 'iso cone', 'iso cone, free rotor']:
201 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi)
202 else:
203 axis = create_rotor_axis_euler(alpha=cdp.eigen_alpha, beta=cdp.eigen_beta, gamma=cdp.eigen_gamma)
204
205
206 if cdp.model in ['rotor', 'free rotor']:
207 span = 20e-10
208 else:
209 span = 35e-10
210
211
212 rotor_pdb(structure=structure, rotor_angle=rotor_angle, axis=axis, axis_pt=pivot, centre=com, span=span, blade_length=5e-10, staggered=True)
213
214
215 structure.add_molecule(name='rest')
216
217
218 mol = model.mol[0]
219 if neg_cone:
220 mol_neg = model_neg.mol[0]
221
222
223
224
225
226
227 structure.add_atom(mol_name=cdp.model, pdb_record='HETATM', atom_num=1, atom_name='R', res_name='PIV', res_num=1, pos=pivot, element='C')
228
229
230
231
232
233
234 if cdp.model in ['iso cone', 'free rotor', 'iso cone, torsionless', 'iso cone, free rotor', 'rotor']:
235
236 print("\nGenerating the z-axis system.")
237
238
239 if cdp.model in ['rotor']:
240 axis = create_rotor_axis_alpha(alpha=cdp.axis_alpha, pivot=pivot, point=com)
241 else:
242 axis = create_rotor_axis_spherical(theta=cdp.axis_theta, phi=cdp.axis_phi)
243 print(("Central axis: %s." % axis))
244
245
246 axis_pos = axis
247 axis_neg = dot(inv_mat, axis)
248
249
250 axis_sim_pos = None
251 axis_sim_neg = None
252 if sim:
253
254 axis_sim = zeros((cdp.sim_number, 3), float64)
255
256
257 for i in range(cdp.sim_number):
258 if cdp.model in ['rotor']:
259 axis_sim[i] = create_rotor_axis_alpha(alpha=cdp.axis_alpha_sim[i], pivot=pivot, point=com)
260 else:
261 axis_sim[i] = create_rotor_axis_spherical(theta=cdp.axis_theta_sim[i], phi=cdp.axis_phi_sim[i])
262
263
264 axis_sim_pos = axis_sim
265 axis_sim_neg = transpose(dot(inv_mat, transpose(axis_sim_pos)))
266
267
268 print("\nGenerating the axis vectors.")
269 res_num = generate_vector_residues(mol=mol, vector=axis_pos, atom_name='z-ax', res_name_vect='AXE', sim_vectors=axis_sim_pos, res_num=2, origin=pivot, scale=size)
270
271
272 if neg_cone:
273 res_num = generate_vector_residues(mol=mol_neg, vector=axis_neg, atom_name='z-ax', res_name_vect='AXE', sim_vectors=axis_sim_neg, res_num=2, origin=pivot, scale=size)
274
275
276 else:
277
278 print("\nGenerating the full axis system.")
279
280
281 axes = zeros((3, 3), float64)
282 euler_to_R_zyz(cdp.eigen_alpha, cdp.eigen_beta, cdp.eigen_gamma, axes)
283 print(("Axis system:\n%s" % axes))
284
285
286 axes_pos = axes
287 axes_neg = dot(inv_mat, axes)
288
289
290 axes_sim_pos = None
291 axes_sim_neg = None
292 if sim:
293
294 axes_sim_pos = zeros((cdp.sim_number, 3, 3), float64)
295 axes_sim_neg = zeros((cdp.sim_number, 3, 3), float64)
296
297
298 for i in range(cdp.sim_number):
299
300 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_pos[i])
301
302
303 euler_to_R_zyz(cdp.eigen_alpha_sim[i], cdp.eigen_beta_sim[i], cdp.eigen_gamma_sim[i], axes_sim_neg[i])
304 axes_sim_neg[i] = dot(inv_mat, axes_sim_neg[i])
305
306
307 print("\nGenerating the axis vectors.")
308 label = ['x', 'y', 'z']
309 for j in range(len(label)):
310
311 axis_sim_pos = None
312 axis_sim_neg = None
313 if sim:
314 axis_sim_pos = axes_sim_pos[:,:, j]
315 axis_sim_neg = axes_sim_neg[:,:, j]
316
317
318 res_num = generate_vector_residues(mol=mol, vector=axes_pos[:, j], atom_name='%s-ax'%label[j], res_name_vect='AXE', sim_vectors=axis_sim_pos, res_num=2, origin=pivot, scale=size)
319 if neg_cone:
320 res_num = generate_vector_residues(mol=mol_neg, vector=axes_neg[:, j], atom_name='%s-ax'%label[j], res_name_vect='AXE', sim_vectors=axis_sim_neg, res_num=2, origin=pivot, scale=size)
321
322
323
324
325
326
327 if cdp.model not in ['rotor', 'free rotor']:
328
329 if cdp.model not in ['iso cone', 'iso cone, torsionless', 'iso cone, free rotor']:
330 R = axes
331 else:
332 R = zeros((3, 3), float64)
333 two_vect_to_R(array([0, 0, 1], float64), axis, R)
334
335
336 R_pos = R
337 R_neg = dot(inv_mat, R)
338
339
340 if cdp.model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor']:
341 cone_obj = Pseudo_elliptic(cdp.cone_theta_x, cdp.cone_theta_y)
342
343
344 else:
345
346 if hasattr(cdp, 'cone_theta'):
347 cone_theta = cdp.cone_theta
348 elif hasattr(cdp, 'cone_s1'):
349 cone_theta = order_parameters.iso_cone_S_to_theta(cdp.cone_s1)
350
351
352 cone_obj = Iso_cone(cone_theta)
353
354
355 cone(mol=mol, cone_obj=cone_obj, start_res=mol.res_num[-1]+1, apex=pivot, R=R_pos, inc=inc, distribution='regular', axis_flag=False)
356
357
358 if neg_cone:
359 cone(mol=mol_neg, cone_obj=cone_obj, start_res=mol_neg.res_num[-1]+1, apex=pivot, R=R_neg, inc=inc, distribution='regular', axis_flag=False)
360
361
362
363
364
365
366 print("\nGenerating the PDB file.")
367
368
369 pdb_file = open_write_file(file, dir, force=force)
370 structure.write_pdb(pdb_file)
371 pdb_file.close()
372
373
374 -def pdb_model(ave_pos_file="ave_pos.pdb", rep_file="frame_order.pdb", dist_file="domain_distribution.pdb", dir=None, size=30.0, inc=36, force=False, neg_cone=True):
375 """Create 3 different PDB files for representing the frame order dynamics of the system.
376
377 @keyword ave_pos_file: The name of the file for the average molecule structure.
378 @type ave_pos_file: str or None
379 @keyword rep_file: The name of the file of the PDB representation of the frame order dynamics to create.
380 @type rep_file: str or None
381 @keyword dist_file: The name of the file which will contain multiple models spanning the full dynamics distribution of the frame order model.
382 @type dist_file: str or None
383 @keyword dir: The name of the directory to place the PDB file into.
384 @type dir: str
385 @keyword size: The size of the geometric object in Angstroms.
386 @type size: float
387 @keyword inc: The number of increments for the filling of the cone objects.
388 @type inc: int
389 @keyword force: Flag which if set to True will cause any pre-existing file to be overwritten.
390 @type force: bool
391 @keyword neg_cone: A flag which if True will cause the negative cone to be added to the representation.
392 @type neg_cone: bool
393 """
394
395
396 if not ave_pos_file and not rep_file and not dist_file:
397 raise RelaxError("Minimally one PDB file name must be supplied.")
398
399
400 pipes.test()
401
402
403 if ave_pos_file:
404 pdb_ave_pos(file=ave_pos_file, dir=dir, force=force)
405
406
407 if cdp.model == 'rigid':
408 return
409
410
411 if rep_file:
412 pdb_geometric_rep(file=rep_file, dir=dir, size=size, inc=inc, force=force, neg_cone=neg_cone)
413
414
415 if dist_file:
416 pdb_distribution(file=dist_file, dir=dir, force=force)
417
418
419 -def pivot(pivot=None, order=1, fix=False):
420 """Set the pivot point for the 2 body motion.
421
422 @keyword pivot: The pivot point of the two bodies (domains, etc.) in the structural coordinate system.
423 @type pivot: list of num
424 @keyword order: The ordinal number of the pivot point. The value of 1 is for the first pivot point, the value of 2 for the second pivot point, and so on.
425 @type order: int
426 @keyword fix: A flag specifying if the pivot point should be fixed during optimisation.
427 @type fix: bool
428 """
429
430
431 pipes.test()
432
433
434 cdp.pivot_fixed = fix
435
436
437 if pivot == None:
438 if hasattr(cdp, 'model'):
439 update_model()
440 return
441
442
443 pivot = array(pivot, float64)
444
445
446 is_float_array(pivot, name='pivot point', size=3)
447
448
449 if order == 1:
450 cdp.pivot_x = pivot[0]
451 cdp.pivot_y = pivot[1]
452 cdp.pivot_z = pivot[2]
453 else:
454
455 name_x = 'pivot_x_%i' % order
456 name_y = 'pivot_y_%i' % order
457 name_z = 'pivot_z_%i' % order
458
459
460 setattr(cdp, name_x, pivot[0])
461 setattr(cdp, name_y, pivot[1])
462 setattr(cdp, name_z, pivot[2])
463
464
465 if hasattr(cdp, 'model'):
466 update_model()
467
468
470 """Turn the high precision Scipy quadratic numerical integration on or off.
471
472 @keyword flag: The flag which if True will perform high precision numerical integration via the scipy.integrate quad(), dblquad() and tplquad() integration methods rather than the rough quasi-random numerical integration.
473 @type flag: bool
474 """
475
476
477 pipes.test()
478
479
480 cdp.quad_int = flag
481
482
483 -def ref_domain(ref=None):
484 """Set the reference domain for the frame order, multi-domain models.
485
486 @param ref: The reference domain.
487 @type ref: str
488 """
489
490
491 pipes.test()
492
493
494 if not hasattr(cdp, 'domain') or ref not in list(cdp.domain.keys()):
495 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % ref)
496
497
498 exists = False
499 for tensor_cont in cdp.align_tensors:
500 if hasattr(tensor_cont, 'domain') and tensor_cont.domain == ref:
501 exists = True
502 if not exists:
503 raise RelaxError("The reference domain cannot be found within any of the loaded tensors.")
504
505
506 cdp.ref_domain = ref
507
508
509 if hasattr(cdp, 'model'):
510 update_model()
511
512
514 """Select the Frame Order model.
515
516 @param model: The Frame Order model. This can be one of 'pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor', 'double rotor'.
517 @type model: str
518 """
519
520
521 pipes.test()
522
523
524 if not model in ['pseudo-ellipse', 'pseudo-ellipse, torsionless', 'pseudo-ellipse, free rotor', 'iso cone', 'iso cone, torsionless', 'iso cone, free rotor', 'line', 'line, torsionless', 'line, free rotor', 'rotor', 'rigid', 'free rotor', 'double rotor']:
525 raise RelaxError("The model name " + repr(model) + " is invalid.")
526
527
528 cdp.model = model
529
530
531 cdp.params = []
532
533
534 if not hasattr(cdp, 'quad_int'):
535
536 if cdp.model in []:
537 cdp.quad_int = True
538
539
540 else:
541 cdp.quad_int = False
542
543
544 update_model()
545