lib.structure.represent.cone

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2013 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 # Python module imports. 23 from math import cos, pi, sin 24 from numpy import array, dot, eye, float64, zeros 25 26 # relax module imports. 27 from lib.geometry.rotations import two_vect_to_R 28 from pipe_control.structure.geometric import angles_uniform, get_proton_name 29 30

31 -def cone_edge(mol=None, cone=None, res_name='CON', res_num=None, chain_id='', apex=None, axis=None, R=None, scale=None, inc=None, distribution='uniform', debug=False):

32 """Add a residue to the atomic data representing a cone of the given angle. 33 34 A series of vectors totalling the number of increments and starting at the origin are equally spaced around the cone axis. The atoms representing neighbouring vectors will be directly bonded together. This will generate an object representing the outer edge of a cone. 35 36 37 @keyword mol: The molecule container. 38 @type mol: MolContainer instance 39 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the phi_max() method for returning the phi value for the given theta. 40 @type cone: class instance 41 @keyword res_name: The residue name. 42 @type res_name: str 43 @keyword res_num: The residue number. 44 @type res_num: int 45 @keyword chain_id: The chain identifier. 46 @type chain_id: str 47 @keyword apex: The apex of the cone. 48 @type apex: numpy array, len 3 49 @keyword axis: The central axis of the cone. If supplied, then this arg will be used to construct the rotation matrix. 50 @type axis: numpy array, len 3 51 @keyword R: A 3x3 rotation matrix. If the axis arg supplied, then this matrix will be ignored. 52 @type R: 3x3 numpy array 53 @keyword scale: The scaling factor to stretch all points by. 54 @type scale: float 55 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 56 @type inc: int 57 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 58 @type distribution: str 59 """ 60 61 # The atom numbers (and indices). 62 atom_num = mol.atom_num[-1]+1 63 64 # Add an atom for the cone apex. 65 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name='APX', res_name=res_name, res_num=res_num, pos=apex, segment_id=None, element='H') 66 origin_atom = atom_num 67 68 # Get the polar and azimuthal angles for the distribution. 69 if distribution == 'uniform': 70 phi, theta = angles_uniform(inc) 71 else: 72 phi, theta = angles_regular(inc) 73 74 # Initialise the rotation matrix. 75 if R == None: 76 R = eye(3) 77 78 # Get the rotation matrix. 79 if axis != None: 80 two_vect_to_R(array([0, 0, 1], float64), axis, R) 81 82 # Determine the maximum phi values and the indices of the point just above the edge. 83 phi_max = zeros(len(theta), float64) 84 edge_index = zeros(len(theta), int) 85 for i in range(len(theta)): 86 # Get the polar angle for the longitude edge atom. 87 phi_max[i] = cone.phi_max(theta[i]) 88 89 # The index. 90 for j in range(len(phi)): 91 if phi[j] <= phi_max[i]: 92 edge_index[i] = j 93 break 94 95 # Reverse edge index. 96 edge_index_rev = len(phi) - 1 - edge_index 97 98 # Move around the azimuth. 99 atom_num = atom_num + 1 100 for i in range(len(theta)): 101 # The vector in the unrotated frame. 102 x = cos(theta[i]) * sin(phi_max[i]) 103 y = sin(theta[i])* sin(phi_max[i]) 104 z = cos(phi_max[i]) 105 vector = array([x, y, z], float64) 106 107 # Rotate the vector. 108 vector = dot(R, vector) 109 110 # The atom id. 111 atom_id = 'T' + repr(i) 112 113 # The atom position. 114 pos = apex+vector*scale 115 116 # Debugging. 117 if debug: 118 print("i: %s; theta: %s" % (i, theta[i])) 119 print("%sAdding atom %s." % (" "*4, get_proton_name(atom_num))) 120 121 # Add the vector as a H atom of the cone residue. 122 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=get_proton_name(atom_num), res_name=res_name, res_num=res_num, pos=pos, segment_id=None, element='H') 123 124 # Join the longitude atom to the cone apex. 125 mol.atom_connect(index1=origin_atom-1, index2=atom_num-1) 126 127 # Increment the atom number. 128 atom_num = atom_num + 1 129 130 # Add latitude atoms for a smoother cone edge and better stitching to the cap. 131 for j in range(len(phi)): 132 # The index for the direction top to bottom! 133 k = len(phi) - 1 - j 134 135 # Debugging. 136 if debug: 137 print("%sj: %s; phi: %-20s; k: %s; phi: %-20s; phi_max: %-20s" % (" "*4, j, phi[j], k, phi[k], phi_max[i])) 138 139 # No edge. 140 skip = True 141 142 # Forward edge (skip when the latitude is phi max). 143 fwd_index = i+1 144 if i == len(theta)-1: 145 fwd_index = 0 146 if j >= edge_index[i] and j < edge_index[fwd_index] and not abs(phi_max[fwd_index] - phi[j]) < 1e-6: 147 # Debugging. 148 if debug: 149 print("%sForward edge." % (" "*8)) 150 151 # Edge found. 152 skip = False 153 154 # Find the theta value for this polar angle phi. 155 phi_val = phi[j] 156 if fwd_index == 0: 157 theta_max = theta[fwd_index] + 2*pi 158 else: 159 theta_max = theta[fwd_index] 160 theta_max = cone.theta_max(phi_val, theta_min=theta[i], theta_max=theta_max-1e-7) 161 162 # Back edge (skip when the latitude is phi max). 163 rev_index = i-1 164 if i == 0: 165 rev_index = len(theta)-1 166 if i < len(theta)-1 and j > edge_index_rev[i] and j <= edge_index_rev[i+1] and not abs(phi_max[fwd_index] - phi[k]) < 1e-6: 167 # Debugging. 168 if debug: 169 print("%sBack edge." % (" "*8)) 170 171 # Edge found. 172 skip = False 173 174 # Find the theta value for this polar angle phi. 175 phi_val = phi[k] 176 theta_max = cone.theta_max(phi_val, theta_min=theta[i-1], theta_max=theta[i]) 177 178 # Skipping. 179 if skip: 180 continue 181 182 # Debugging. 183 if debug: 184 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 185 186 # The coordinates. 187 x = cos(theta_max) * sin(phi_val) 188 y = sin(theta_max) * sin(phi_val) 189 z = cos(phi_val) 190 pos = array([x, y, z], float64) * scale 191 192 # Rotate and shift. 193 pos = apex + dot(R, pos) 194 195 # Add the atom. 196 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=get_proton_name(atom_num), res_name=res_name, chain_id=chain_id, res_num=res_num, pos=pos, segment_id=None, element='H') 197 198 # Increment the atom number. 199 atom_num = atom_num + 1 200 201 # Debugging. 202 if debug: 203 print("\nBuilding the edge.") 204 205 # Build the cone edge. 206 for i in range(origin_atom, atom_num-2): 207 # Debugging. 208 if debug: 209 print("%sCone edge, connecting %s to %s" % (" "*4, get_proton_name(i), get_proton_name(i+1))) 210 211 # Connect. 212 mol.atom_connect(index1=i, index2=i+1) 213 214 # Connect the last radial array to the first (to zip up the circle). 215 mol.atom_connect(index1=atom_num-2, index2=origin_atom)

216 217

218 -def cone(mol=None, cone=None, start_res=1, apex=None, axis=None, R=None, inc=None, scale=30.0, distribution='regular', axis_flag=True):

219 """Create a structural representation of the given cone object. 220 221 @keyword mol: The molecule container. 222 @type mol: MolContainer instance 223 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the theta_max() method for returning the theta value for the given phi, the phi_max() method for returning the phi value for the given theta. 224 @type cone: class instance 225 @keyword start_res: The starting residue number. 226 @type start_res: str 227 @keyword apex: The apex of the cone. 228 @type apex: rank-1, 3D numpy array 229 @keyword axis: The central axis of the cone. If not supplied, the z-axis will be used. 230 @type axis: rank-1, 3D numpy array 231 @keyword R: The rotation matrix. 232 @type R: rank-2, 3D numpy array 233 @keyword inc: The increment number used to determine the number of latitude and longitude lines. 234 @type inc: int 235 @keyword scale: The scaling factor to stretch the unit cone by. 236 @type scale: float 237 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 238 @type distribution: str 239 @keyword axis_flag: A flag which if True will create the cone's axis. 240 @type axis_flag: bool 241 """ 242 243 # The cone axis default of the z-axis. 244 if not axis: 245 axis = array([0, 0, 1], float64) 246 247 # No rotation. 248 if R == None: 249 R = eye(3) 250 251 # The first atom number. 252 start_atom = 1 253 if hasattr(mol, 'atom_num'): 254 start_atom = mol.atom_num[-1]+1 255 256 # The axis. 257 if axis_flag: 258 # Add the apex. 259 mol.atom_add(pdb_record='HETATM', atom_num=start_atom, atom_name='R', res_name='APX', res_num=start_res, pos=apex, element='C') 260 261 # Generate the axis vectors. 262 print("\nGenerating the axis vectors.") 263 res_num = generate_vector_residues(mol=mol, vector=dot(R, axis), atom_name='Axis', res_name_vect='AXE', res_num=start_res+1, origin=apex, scale=scale) 264 265 # Generate the cone outer edge. 266 print("\nGenerating the cone outer edge.") 267 edge_start_atom = mol.atom_num[-1]+1 268 cone_edge(mol=mol, cone=cone, res_name='EDG', res_num=start_res+2, apex=apex, R=R, scale=scale, inc=inc, distribution=distribution) 269 270 # Generate the cone cap, and stitch it to the cone edge. 271 print("\nGenerating the cone cap.") 272 cone_start_atom = mol.atom_num[-1]+1 273 generate_vector_dist(mol=mol, res_name='CON', res_num=start_res+3, centre=apex, R=R, limit_check=cone.limit_check, scale=scale, inc=inc, distribution=distribution) 274 stitch_cone_to_edge(mol=mol, cone=cone, dome_start=cone_start_atom, edge_start=edge_start_atom+1, scale=scale, inc=inc, distribution=distribution)

275 276

277 -def stitch_cone_to_edge(mol=None, cone=None, chain_id='', dome_start=None, edge_start=None, scale=1.0, inc=None, distribution='uniform', debug=False):

278 """Function for stitching the cone dome to its edge, in the PDB representations. 279 280 @keyword mol: The molecule container. 281 @type mol: MolContainer instance 282 @keyword cone: The cone object. This should provide the limit_check() method with determines the limits of the distribution accepting two arguments, the polar angle phi and the azimuthal angle theta, and return True if the point is in the limits or False if outside. It should also provide the theta_max() method for returning the theta value for the given phi. 283 @type cone: class instance 284 @keyword chain_id: The chain identifier. 285 @type chain_id: str 286 @keyword dome_start: The starting atom number of the cone dome residue. 287 @type dome_start: int 288 @keyword edge_start: The starting atom number of the cone edge residue. 289 @type edge_start: int 290 @keyword scale: The scaling factor to stretch all points by. 291 @type scale: float 292 @keyword inc: The number of increments or number of vectors used to generate the outer edge of the cone. 293 @type inc: int 294 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 295 @type distribution: str 296 """ 297 298 # Get the polar and azimuthal angles for the distribution. 299 if distribution == 'uniform': 300 phi, theta = angles_uniform(inc) 301 else: 302 phi, theta = angles_regular(inc) 303 304 # Determine the maximum phi values and the indices of the point just above the edge. 305 phi_max = zeros(len(theta), float64) 306 edge_index = zeros(len(theta), int) 307 for i in range(len(theta)): 308 # Get the polar angle for the longitude edge atom. 309 phi_max[i] = cone.phi_max(theta[i]) 310 311 # The index. 312 for j in range(len(phi)): 313 if phi[j] <= phi_max[i]: 314 edge_index[i] = j 315 break 316 317 # Reverse edge index. 318 edge_index_rev = len(phi) - 1 - edge_index 319 320 # Debugging. 321 if debug: 322 print("\nDome start: %s" % dome_start) 323 print("Edge start: %s" % edge_start) 324 print("Edge indices: %s" % edge_index) 325 print("Edge indices rev: %s" % edge_index_rev) 326 327 # Move around the azimuth. 328 dome_atom = dome_start 329 edge_atom = edge_start 330 for i in range(len(theta)): 331 # Debugging. 332 if debug: 333 print("i: %s; theta: %s" % (i, theta[i])) 334 print("%sDome atom: %s" % (" "*4, get_proton_name(dome_atom))) 335 print("%sStitching longitudinal line to edge - %s to %s." % (" "*4, get_proton_name(edge_atom), get_proton_name(dome_atom))) 336 337 # Connect the two atoms (to stitch up the 2 objects). 338 mol.atom_connect(index1=dome_atom-1, index2=edge_atom-1) 339 340 # Update the edge atom. 341 edge_atom = edge_atom + 1 342 343 # Stitch up the latitude atoms. 344 for j in range(len(phi)): 345 # The index for the direction top to bottom! 346 k = len(phi) - 1 - j 347 348 # Debugging. 349 if debug: 350 print("%sj: %s; phi: %-20s; k: %s; phi: %-20s; phi_max: %-20s" % (" "*4, j, phi[j], k, phi[k], phi_max[i])) 351 352 # No edge. 353 skip = True 354 355 # Forward edge (skip when the latitude is phi max). 356 fwd_index = i+1 357 if i == len(theta)-1: 358 fwd_index = 0 359 if j >= edge_index[i] and j < edge_index[fwd_index] and not abs(phi_max[fwd_index] - phi[j]) < 1e-6: 360 # Debugging. 361 if debug: 362 print("%sForward edge." % (" "*8)) 363 364 # Edge found. 365 skip = False 366 forward = True 367 368 # Back edge (skip when the latitude is phi max). 369 rev_index = i-1 370 if i == 0: 371 rev_index = len(theta)-1 372 if i < len(theta)-1 and j > edge_index_rev[i] and j <= edge_index_rev[i+1] and not abs(phi_max[fwd_index] - phi[k]) < 1e-6: 373 # Debugging. 374 if debug: 375 print("%sBack edge." % (" "*8)) 376 377 # Edge found. 378 skip = False 379 forward = False 380 381 # Skipping. 382 if skip: 383 continue 384 385 # The dome atom to stitch to (current dome atom + one latitude line to shift across). 386 if forward: 387 atom = dome_atom + j - edge_index[i] 388 else: 389 atom = dome_atom + (edge_index_rev[i]+1) + k - edge_index[fwd_index] 390 391 # Debugging. 392 if debug: 393 print("%sStitching latitude line to edge - %s to %s." % (" "*8, get_proton_name(edge_atom), get_proton_name(atom))) 394 395 # Connect the two atoms (to stitch up the 2 objects). 396 mol.atom_connect(index1=atom-1, index2=edge_atom-1) 397 398 # Increment the cone edge atom number. 399 edge_atom = edge_atom + 1 400 401 # Update the cone dome atom. 402 dome_atom = dome_atom + (len(phi) - edge_index[i])

403

Source Code for Module lib.structure.represent.cone