lib.structure.geometric

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 arccos, array, dot, eye, float64, zeros 25 26 # relax module imports. 27 from lib.structure.angles import angles_regular, angles_uniform 28 from lib.structure.conversion import get_proton_name 29 30

31 -def generate_vector_dist(mol=None, res_name=None, res_num=None, chain_id='', centre=zeros(3, float64), R=eye(3), warp=eye(3), limit_check=None, scale=1.0, inc=20, distribution='uniform', debug=False):

32 """Generate a uniformly distributed distribution of atoms on a warped sphere. 33 34 The vectors from the function vect_dist_spherical_angles() are used to generate the distribution. These vectors are rotated to the desired frame using the rotation matrix 'R', then each compressed or stretched by the dot product with the 'warp' matrix. Each vector is centred and at the head of the vector, a proton is placed. 35 36 37 @keyword mol: The molecule container. 38 @type mol: MolContainer instance 39 @keyword res_name: The residue name. 40 @type res_name: str 41 @keyword res_num: The residue number. 42 @type res_num: int 43 @keyword chain_id: The chain identifier. 44 @type chain_id: str 45 @keyword centre: The centre of the distribution. 46 @type centre: numpy array, len 3 47 @keyword R: The optional 3x3 rotation matrix. 48 @type R: 3x3 numpy array 49 @keyword warp: The optional 3x3 warping matrix. 50 @type warp: 3x3 numpy array 51 @keyword limit_check: A function with determines the limits of the distribution. It should accept 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. 52 @type limit_check: function 53 @keyword scale: The scaling factor to stretch all rotated and warped vectors by. 54 @type scale: float 55 @keyword inc: The number of increments or number of vectors. 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 # Initial atom number. 62 if len(mol.atom_num) == 0: 63 origin_num = 1 64 else: 65 origin_num = mol.atom_num[-1]+1 66 atom_num = origin_num 67 68 # Get the uniform vector distribution. 69 print(" Creating the uniform vector distribution.") 70 vectors = vect_dist_spherical_angles(inc=inc, distribution=distribution) 71 72 # Get the polar and azimuthal angles for the distribution. 73 if distribution == 'uniform': 74 phi, theta = angles_uniform(inc) 75 else: 76 phi, theta = angles_regular(inc) 77 78 # Init the arrays for stitching together. 79 edge = zeros(len(theta)) 80 edge_index = zeros(len(theta), int) 81 edge_phi = zeros(len(theta), float64) 82 edge_atom = zeros(len(theta)) 83 84 # Loop over the radial array of vectors (change in longitude). 85 for i in range(len(theta)): 86 # Debugging. 87 if debug: 88 print("i: %s; theta: %s" % (i, theta[i])) 89 90 # Loop over the vectors of the radial array (change in latitude). 91 for j in range(len(phi)): 92 # Debugging. 93 if debug: 94 print("%sj: %s; phi: %s" % (" "*4, j, phi[j])) 95 96 # Skip the vector if the point is outside of the limits. 97 if limit_check and not limit_check(phi[j], theta[i]): 98 if debug: 99 print("%sOut of limits." % (" "*8)) 100 continue 101 102 # Update the edge for this longitude. 103 if not edge[i]: 104 edge[i] = 1 105 edge_index[i] = j 106 edge_phi[i] = phi[j] 107 edge_atom[i] = atom_num 108 109 # Debugging. 110 if debug: 111 print("%sEdge detected." % (" "*8)) 112 print("%sEdge index: %s" % (" "*8, edge_index[i])) 113 print("%sEdge phi pos: %s" % (" "*8, edge_phi[i])) 114 print("%sEdge atom: %s" % (" "*8, edge_atom[i])) 115 116 # Rotate the vector into the frame. 117 vector = dot(R, vectors[i + j*len(theta)]) 118 119 # Warp the vector. 120 vector = dot(warp, vector) 121 122 # Scale the vector. 123 vector = vector * scale 124 125 # Position relative to the centre of mass. 126 pos = centre + vector 127 128 # Debugging. 129 if debug: 130 print("%sAdding atom %s." % (" "*8, get_proton_name(atom_num))) 131 132 # Add the vector as a H atom of the TNS residue. 133 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') 134 135 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 136 if j > edge_index[i]: 137 # Debugging. 138 if debug: 139 print("%sLongitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-1))) 140 141 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 142 143 # Connect across the radial arrays (to generate the latitudinal lines). 144 if i != 0 and j >= edge_index[i-1]: 145 # The number of atoms back to the previous longitude. 146 num = len(phi) - edge_index[i] 147 148 # Debugging. 149 if debug: 150 print("%sLatitude line, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(atom_num-num))) 151 152 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 153 154 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 155 if i == len(theta)-1 and j >= edge_index[0]: 156 # The number of atom in the first longitude line. 157 num = origin_num + j - edge_index[0] 158 159 # Debugging. 160 if debug: 161 print("%sZipping up, connecting %s to %s" % (" "*8, get_proton_name(atom_num), get_proton_name(num))) 162 163 mol.atom_connect(index1=atom_num-1, index2=num-1) 164 165 # Increment the atom number. 166 atom_num = atom_num + 1

167 168

169 -def generate_vector_residues(mol=None, vector=None, atom_name=None, res_name_vect='AXS', sim_vectors=None, res_name_sim='SIM', chain_id='', res_num=None, origin=None, scale=1.0, label_placement=1.1, neg=False):

170 """Generate residue representations for the vector and the MC simulationed vectors. 171 172 This is used to create a PDB representation of any vector, including its Monte Carlo simulations. 173 174 @keyword mol: The molecule container. 175 @type mol: MolContainer instance 176 @keyword vector: The vector to be represented in the PDB. 177 @type vector: numpy array, len 3 178 @keyword atom_name: The atom name used to label the atom representing the head of the vector. 179 @type atom_name: str 180 @keyword res_name_vect: The 3 letter PDB residue code used to represent the vector. 181 @type res_name_vect: str 182 @keyword sim_vectors: The optional Monte Carlo simulation vectors to be represented in the PDB. 183 @type sim_vectors: list of numpy array, each len 3 184 @keyword res_name_sim: The 3 letter PDB residue code used to represent the Monte Carlo simulation vectors. 185 @type res_name_sim: str 186 @keyword chain_id: The chain identification code. 187 @type chain_id: str 188 @keyword res_num: The residue number. 189 @type res_num: int 190 @keyword origin: The origin for the axis. 191 @type origin: numpy array, len 3 192 @keyword scale: The scaling factor to stretch the vectors by. 193 @type scale: float 194 @keyword label_placement: A scaling factor to multiply the pre-scaled vector by. This is used to place the vector labels a little further out from the vector itself. 195 @type label_placement: float 196 @keyword neg: If True, then the negative vector positioned at the origin will also be included. 197 @type neg: bool 198 @return: The new residue number. 199 @rtype: int 200 """ 201 202 # The atom numbers (and indices). 203 origin_num = len(mol.atom_num)+1 204 atom_num = len(mol.atom_num)+2 205 atom_neg_num = len(mol.atom_num)+3 206 207 # The origin atom. 208 mol.atom_add(pdb_record='HETATM', atom_num=origin_num, atom_name='R', res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin, segment_id=None, element='C') 209 210 # Create the PDB residue representing the vector. 211 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin+vector*scale, segment_id=None, element='C') 212 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 213 num = 1 214 if neg: 215 mol.atom_add(pdb_record='HETATM', atom_num=atom_neg_num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin-vector*scale, segment_id=None, element='C') 216 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 217 num = 2 218 219 # Add another atom to allow the axis labels to be shifted just outside of the vector itself. 220 mol.atom_add(pdb_record='HETATM', atom_num=atom_num+num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin+label_placement*vector*scale, segment_id=None, element='N') 221 if neg: 222 mol.atom_add(pdb_record='HETATM', atom_num=atom_neg_num+num, atom_name=atom_name, res_name=res_name_vect, chain_id=chain_id, res_num=res_num, pos=origin-label_placement*vector*scale, segment_id=None, element='N') 223 224 # Print out. 225 print(" " + atom_name + " vector (scaled + shifted to origin): " + repr(origin+vector*scale)) 226 print(" Creating the MC simulation vectors.") 227 228 # Monte Carlo simulations. 229 if sim_vectors != None: 230 for i in range(len(sim_vectors)): 231 # Increment the residue number, so each simulation is a new residue. 232 res_num = res_num + 1 233 234 # The atom numbers (and indices). 235 atom_num = mol.atom_num[-1]+1 236 atom_neg_num = mol.atom_num[-1]+2 237 238 # Create the PDB residue representing the vector. 239 mol.atom_add(pdb_record='HETATM', atom_num=atom_num, atom_name=atom_name, res_name=res_name_sim, chain_id=chain_id, res_num=res_num, pos=origin+sim_vectors[i]*scale, segment_id=None, element='C') 240 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 241 if neg: 242 mol.atom_add(pdb_record='HETATM', atom_num=atom_num+1, atom_name=atom_name, res_name=res_name_sim, chain_id=chain_id, res_num=res_num, pos=origin-sim_vectors[i]*scale, segment_id=None, element='C') 243 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 244 245 # Return the new residue number. 246 return res_num

247 248

249 -def vect_dist_spherical_angles(inc=20, distribution='uniform'):

250 """Create a distribution of vectors on a sphere using a distribution of spherical angles. 251 252 This function returns an array of unit vectors distributed within 3D space. The unit vectors are generated using the equation:: 253 254 | cos(theta) * sin(phi) | 255 vector = | sin(theta) * sin(phi) |. 256 | cos(phi) | 257 258 The vectors of this distribution generate both longitudinal and latitudinal lines. 259 260 261 @keyword inc: The number of increments in the distribution. 262 @type inc: int 263 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 264 @type distribution: str 265 @return: The distribution of vectors on a sphere. 266 @rtype: list of rank-1, 3D numpy arrays 267 """ 268 269 # Get the polar and azimuthal angles for the distribution. 270 if distribution == 'uniform': 271 phi, theta = angles_uniform(inc) 272 else: 273 phi, theta = angles_regular(inc) 274 275 # Initialise array of the distribution of vectors. 276 vectors = [] 277 278 # Loop over the longitudinal lines. 279 for j in range(len(phi)): 280 # Loop over the latitudinal lines. 281 for i in range(len(theta)): 282 # X coordinate. 283 x = cos(theta[i]) * sin(phi[j]) 284 285 # Y coordinate. 286 y = sin(theta[i])* sin(phi[j]) 287 288 # Z coordinate. 289 z = cos(phi[j]) 290 291 # Append the vector. 292 vectors.append(array([x, y, z], float64)) 293 294 # Return the array of vectors and angles. 295 return vectors

296

Source Code for Module lib.structure.geometric