lib.structure.represent.rotor

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Representations of a mechanical rotor.""" 24 25 # Python module imports. 26 from math import pi 27 from numpy import array, cross, dot, float64, transpose, zeros 28 from numpy.linalg import norm 29 30 # relax module imports. 31 from lib.geometry.lines import closest_point_ax 32 from lib.geometry.rotations 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: lib.structure.internal.object.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 # Convert the arguments to numpy arrays, radians and Angstrom. 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 # Normalise. 66 axis_norm = axis / norm(axis) 67 68 # Add a structure (handling up to 3 rotors). 69 if structure.has_molecule(name='rotor') and structure.has_molecule(name='rotor2'): 70 structure.add_molecule(name='rotor3') 71 mol_name = 'rotor3' 72 elif structure.has_molecule(name='rotor'): 73 structure.add_molecule(name='rotor2') 74 mol_name = 'rotor2' 75 else: 76 structure.add_molecule(name='rotor') 77 mol_name = 'rotor' 78 79 # Loop over the models. 80 for model in structure.model_loop(model): 81 # Alias the single molecule from the single model. 82 mol = structure.get_molecule(mol_name, model=model.num) 83 84 # The central point. 85 mid_point = closest_point_ax(line_pt=axis_pt, axis=axis, point=centre) 86 pos_index = mol.atom_add(pdb_record='HETATM', atom_name='CTR', res_name='AX', res_num=1, pos=mid_point, element='PT') 87 88 # Centre of the propellers. 89 prop1 = mid_point + axis_norm * span 90 prop1_index = pos_index + 1 91 mol.atom_add(pdb_record='HETATM', atom_name='PRP', res_name='PRC', res_num=2, pos=prop1, element='O') 92 mol.atom_connect(index1=pos_index, index2=prop1_index) 93 94 # Centre of the propellers. 95 prop2 = mid_point - axis_norm * span 96 prop2_index = pos_index + 2 97 mol.atom_add(pdb_record='HETATM', atom_name='PRP', res_name='PRC', res_num=3, pos=prop2, element='O') 98 mol.atom_connect(index1=pos_index, index2=prop2_index) 99 100 # Create the rotor propellers. 101 rotor_propellers(mol=mol, rotor_angle=rotor_angle, centre=prop1, axis=axis, blade_length=blade_length, staggered=staggered) 102 rotor_propellers(mol=mol, rotor_angle=rotor_angle, centre=prop2, axis=-axis, blade_length=blade_length, staggered=staggered)

103 104

105 -def rotor_propellers(mol=None, rotor_angle=None, centre=None, axis=None, blade_length=5.0, staggered=False):

106 """Create a PDB representation of a rotor motional model. 107 108 @keyword mol: The internal structural object molecule container to add the atoms to. 109 @type mol: MolContainer instance 110 @keyword rotor_angle: The angle of the rotor motion in radians. 111 @type rotor_angle: float 112 @keyword centre: The central point of the propeller. 113 @type centre: numpy rank-1, 3D array 114 @keyword axis: The vector defining the rotor axis. 115 @type axis: numpy rank-1, 3D array 116 @keyword blade_length: The length of the representative rotor blades in Angstrom. 117 @type blade_length: float 118 @keyword staggered: A flag which if True will cause the rotor blades to be staggered. This is used to avoid blade overlap. 119 @type staggered: bool 120 """ 121 122 # Init. 123 step_angle = 2.0 / 360.0 * 2.0 * pi 124 R = zeros((3, 3), float64) 125 res_num = mol.last_residue() + 1 126 127 # Blade vectors. 128 blades = zeros((4, 3), float64) 129 if abs(dot(axis, array([0, 0, 1], float64))) == 1.0: # Avoid failures in artificial situations. 130 blades[0] = cross(axis, array([1, 0, 0], float64)) 131 else: 132 blades[0] = cross(axis, array([0, 0, 1], float64)) 133 blades[0] = blades[0] / norm(blades[0]) 134 blades[1] = cross(axis, blades[0]) 135 blades[1] = blades[1] / norm(blades[1]) 136 blades[2] = -blades[0] 137 blades[3] = -blades[1] 138 139 # Create the 4 blades. 140 for i in range(len(blades)): 141 # Staggering. 142 if staggered and i % 2: 143 blade_origin = centre - axis * 2 144 145 # Non-staggered. 146 else: 147 blade_origin = centre 148 149 # Add an atom for the blage origin. 150 blade_origin_index = mol.atom_add(pdb_record='HETATM', atom_name='BLO', res_name='PRB', res_num=res_num, pos=blade_origin, element='O') 151 152 # The centre edge point of the blade. 153 mid_point = blade_origin + blades[i] * blade_length 154 mid_pt_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=mid_point, element='N') 155 mol.atom_connect(index1=mid_pt_index, index2=blade_origin_index) 156 157 # Build the blade. 158 angle = 0.0 159 pos_last_index = mid_pt_index 160 neg_last_index = mid_pt_index 161 while True: 162 # Increase the angle. 163 angle += step_angle 164 165 # The edge rotation. 166 if angle > rotor_angle: 167 axis_angle_to_R(axis, rotor_angle, R) 168 169 # The normal rotation matrix. 170 else: 171 axis_angle_to_R(axis, angle, R) 172 173 # The positive edge. 174 pos_point = dot(R, mid_point - blade_origin) + blade_origin 175 pos_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=pos_point, element='N') 176 mol.atom_connect(index1=pos_index, index2=pos_last_index) 177 mol.atom_connect(index1=pos_index, index2=blade_origin_index) 178 179 # The negative edge. 180 neg_point = dot(transpose(R), mid_point - blade_origin) + blade_origin 181 neg_index = mol.atom_add(pdb_record='HETATM', atom_name='BLD', res_name='PRB', res_num=res_num, pos=neg_point, element='N') 182 mol.atom_connect(index1=neg_index, index2=neg_last_index) 183 mol.atom_connect(index1=neg_index, index2=blade_origin_index) 184 185 # Update the indices. 186 pos_last_index = pos_index 187 neg_last_index = neg_index 188 189 # Finish. 190 if angle > rotor_angle: 191 break 192 193 # Increment the residue number. 194 res_num += 1

195

Source Code for Module lib.structure.represent.rotor