lib.structure.represent.cone

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 # 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 lib.structure.angles import angles_regular, angles_uniform 29 from lib.structure.conversion import get_proton_name 30 from lib.structure.geometric import generate_vector_dist 31 32

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

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

218 219

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

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

277 278

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

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

405

Source Code for Module lib.structure.represent.cone