1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Representations of a mechanical rotor."""
24
25
26 from math import pi
27 from numpy import array, cross, dot, float64, transpose, zeros
28 from numpy.linalg import norm
29
30
31 from lib.geometry.lines import closest_point_ax
32 from maths_fns.rotation_matrix import axis_angle_to_R
33
34
35 -def rotor_pdb(structure=None, rotor_angle=None, axis=None, axis_pt=True, centre=None, span=2e-9, blade_length=5e-10, model=None, staggered=False):
36 """Create a PDB representation of a rotor motional model.
37
38 @keyword structure: The internal structural object instance to add the rotor to as a molecule.
39 @type structure: generic_fns.structure.internal.Internal instance
40 @keyword rotor_angle: The angle of the rotor motion in radian.
41 @type rotor_angle: float
42 @keyword axis: The vector defining the rotor axis.
43 @type axis: numpy rank-1, 3D array
44 @keyword axis_pt: A point lying anywhere on the rotor axis. This is used to define the position of the axis in 3D space.
45 @type axis_pt: numpy rank-1, 3D array
46 @keyword centre: The central point of the representation. If this point is not on the rotor axis, then the closest point on the axis will be used for the centre.
47 @type centre: numpy rank-1, 3D array
48 @keyword span: The distance from the central point to the rotor blades (meters).
49 @type span: float
50 @keyword blade_length: The length of the representative rotor blades.
51 @type blade_length: float
52 @keyword model: The structural model number to add the rotor to. If not supplied, the same rotor structure will be added to all models.
53 @type model: int or None
54 @keyword staggered: A flag which if True will cause the rotor blades to be staggered. This is used to avoid blade overlap.
55 @type staggered: bool
56 """
57
58
59 axis = array(axis, float64)
60 axis_pt = array(axis_pt, float64)
61 centre = array(centre, float64)
62 span = span * 1e10
63 blade_length = blade_length * 1e10
64
65
66 axis_norm = axis / norm(axis)
67
68
69 structure.add_molecule(name='rotor')
70
71
72 for model in structure.model_loop(model):
73 print model
74
75
76 mol = structure.get_molecule('rotor', model=model.num)
77
78
79 mid_point = closest_point_ax(line_pt=axis_pt, axis=axis, point=centre)
80 mol.atom_add(pdb_record='HETATM', atom_num=1, atom_name='CTR', res_name='AX', res_num=1, pos=mid_point, element='PT')
81
82
83 prop1 = mid_point + axis_norm * span
84 prop1_index = 1
85 mol.atom_add(pdb_record='HETATM', atom_num=2, atom_name='PRP', res_name='PRC', res_num=2, pos=prop1, element='O')
86 mol.atom_connect(index1=0, index2=prop1_index)
87
88
89 prop2 = mid_point - axis_norm * span
90 prop2_index = 2
91 mol.atom_add(pdb_record='HETATM', atom_num=3, atom_name='PRP', res_name='PRC', res_num=3, pos=prop2, element='O')
92 mol.atom_connect(index1=0, index2=prop2_index)
93
94
95 rotor_propellers(mol=mol, rotor_angle=rotor_angle, centre=prop1, axis=axis, blade_length=blade_length, staggered=staggered)
96 rotor_propellers(mol=mol, rotor_angle=rotor_angle, centre=prop2, axis=-axis, blade_length=blade_length, staggered=staggered)
97
98
99 -def rotor_propellers(mol=None, rotor_angle=None, centre=None, axis=None, blade_length=5.0, staggered=False):
100 """Create a PDB representation of a rotor motional model.
101
102 @keyword mol: The internal structural object molecule container to add the atoms to.
103 @type mol: MolContainer instance
104 @keyword rotor_angle: The angle of the rotor motion in radians.
105 @type rotor_angle: float
106 @keyword centre: The central point of the propeller.
107 @type centre: numpy rank-1, 3D array
108 @keyword axis: The vector defining the rotor axis.
109 @type axis: numpy rank-1, 3D array
110 @keyword blade_length: The length of the representative rotor blades in Angstrom.
111 @type blade_length: float
112 @keyword staggered: A flag which if True will cause the rotor blades to be staggered. This is used to avoid blade overlap.
113 @type staggered: bool
114 """
115
116
117 step_angle = 2.0 / 360.0 * 2.0 * pi
118 R = zeros((3, 3), float64)
119 res_num = mol.last_residue() + 1
120
121
122 blades = zeros((4, 3), float64)
123 if abs(dot(axis, array([0, 0, 1], float64))) == 1.0:
124 blades[0] = cross(axis, array([1, 0, 0], float64))
125 else:
126 blades[0] = cross(axis, array([0, 0, 1], float64))
127 blades[0] = blades[0] / norm(blades[0])
128 blades[1] = cross(axis, blades[0])
129 blades[1] = blades[1] / norm(blades[1])
130 blades[2] = -blades[0]
131 blades[3] = -blades[1]
132
133
134 for i in range(len(blades)):
135
136 if staggered and i % 2:
137 blade_origin = centre - axis * 2
138
139
140 else:
141 blade_origin = centre
142
143
144 blade_origin_index = mol.atom_add(pdb_record='HETATM', atom_name='BLO', res_name='PRB', res_num=res_num, pos=blade_origin, element='O')
145
146
147 mid_point = blade_origin + blades[i] * blade_length
148 mid_pt_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=mid_point, element='N')
149 mol.atom_connect(index1=mid_pt_index, index2=blade_origin_index)
150
151
152 angle = 0.0
153 pos_last_index = mid_pt_index
154 neg_last_index = mid_pt_index
155 while True:
156
157 angle += step_angle
158
159
160 if angle > rotor_angle:
161 axis_angle_to_R(axis, rotor_angle, R)
162
163
164 else:
165 axis_angle_to_R(axis, angle, R)
166
167
168 pos_point = dot(R, mid_point - blade_origin) + blade_origin
169 pos_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=pos_point, element='N')
170 mol.atom_connect(index1=pos_index, index2=pos_last_index)
171 mol.atom_connect(index1=pos_index, index2=blade_origin_index)
172
173
174 neg_point = dot(transpose(R), mid_point - blade_origin) + blade_origin
175 neg_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=neg_point, element='N')
176 mol.atom_connect(index1=neg_index, index2=neg_last_index)
177 mol.atom_connect(index1=neg_index, index2=blade_origin_index)
178
179
180 pos_last_index = pos_index
181 neg_last_index = neg_index
182
183
184 if angle > rotor_angle:
185 break
186
187
188 res_num += 1
189