lib.structure.geometric

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 Edward d'Auvergne # 4 # Copyright (C) 2008 Sebastien Morin # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Python module imports. 24 from math import cos, sin 25 from numpy import array, dot, eye, float64, zeros 26 27 # relax module imports. 28 from lib.structure.angles import angles_regular, angles_uniform 29 from lib.structure.conversion import get_proton_name 30 31

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

33 """Generate a uniformly distributed distribution of atoms on a warped sphere. 34 35 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. 36 37 38 @keyword mol: The molecule container. 39 @type mol: MolContainer instance 40 @keyword res_name: The residue name. 41 @type res_name: str 42 @keyword res_num: The residue number. 43 @type res_num: int 44 @keyword chain_id: The chain identifier. 45 @type chain_id: str 46 @keyword centre: The centre of the distribution. 47 @type centre: numpy array, len 3 48 @keyword R: The optional 3x3 rotation matrix. 49 @type R: 3x3 numpy array 50 @keyword warp: The optional 3x3 warping matrix. 51 @type warp: numpy array 52 @keyword phi_max_fn: A function with determines the limits of the distribution. It should accept the azimuthal angle theta as an argument and return the corresponding maximum allowed polar angle phi. 53 @type phi_max_fn: function 54 @keyword scale: The scaling factor to stretch all rotated and warped vectors by. 55 @type scale: float 56 @keyword inc: The number of increments or number of vectors. 57 @type inc: int 58 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 59 @type distribution: str 60 """ 61 62 # Initial atom number. 63 if len(mol.atom_num) == 0: 64 origin_num = 1 65 else: 66 origin_num = mol.atom_num[-1]+1 67 atom_num = origin_num 68 69 # Get the uniform vector distribution. 70 print(" Creating the uniform vector distribution.") 71 vectors = vect_dist_spherical_angles(inc=inc, distribution=distribution) 72 73 # Get the polar and azimuthal angles for the distribution. 74 if distribution == 'uniform': 75 phi, theta = angles_uniform(inc) 76 else: 77 phi, theta = angles_regular(inc) 78 79 # Init the arrays for stitching together. 80 edge = zeros(len(theta)) 81 edge_index = zeros(len(theta), int) 82 edge_atom = zeros(len(theta)) 83 84 # Determine the maximum phi values of the point just above the edge, and create a set of vectors to use for all points outside of the limits. 85 if phi_max_fn: 86 phi_max = zeros(len(theta), float64) 87 edge_vectors = zeros((len(theta), 3), float64) 88 for i in range(len(theta)): 89 # The maximum angle. 90 phi_max[i] = phi_max_fn(theta[i]) 91 92 # The vector in the unrotated frame. 93 edge_vectors[i, 0] = cos(theta[i]) * sin(phi_max[i]) 94 edge_vectors[i, 1] = sin(theta[i])* sin(phi_max[i]) 95 edge_vectors[i, 2] = cos(phi_max[i]) 96 97 # Loop over the radial array of vectors (change in longitude). 98 for i in range(len(theta)): 99 # Loop over the vectors of the radial array (change in latitude). 100 for j in range(len(phi)): 101 # The vector to use. 102 if phi_max_fn and phi[j] > phi_max_fn(theta[i]): 103 vector = edge_vectors[i] 104 else: 105 vector = vectors[i + j*len(theta)] 106 107 # Update the edge for this longitude. 108 if not edge[i]: 109 edge[i] = 1 110 edge_index[i] = j 111 edge_atom[i] = atom_num 112 113 # Rotate the vector into the frame. 114 vector = dot(R, vector) 115 116 # Warp the vector. 117 vector = dot(warp, vector) 118 119 # Scale the vector. 120 vector = vector * scale 121 122 # Position relative to the centre of mass. 123 pos = centre + vector 124 125 # Add the vector as a H atom of the TNS residue. 126 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') 127 128 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 129 if j > edge_index[i]: 130 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 131 132 # Connect across the radial arrays (to generate the latitudinal lines). 133 if i != 0 and j >= edge_index[i-1]: 134 # The number of atoms back to the previous longitude. 135 num = len(phi) - edge_index[i] 136 137 # Latitude line. 138 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 139 140 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 141 if i == len(theta)-1 and j >= edge_index[0]: 142 # The number of atom in the first longitude line. 143 num = origin_num + j - edge_index[0] 144 145 # Zipping up. 146 mol.atom_connect(index1=atom_num-1, index2=num-1) 147 148 # Increment the atom number. 149 atom_num = atom_num + 1

150 151

152 -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):

153 """Generate residue representations for the vector and the MC simulationed vectors. 154 155 This is used to create a PDB representation of any vector, including its Monte Carlo simulations. 156 157 @keyword mol: The molecule container. 158 @type mol: MolContainer instance 159 @keyword vector: The vector to be represented in the PDB. 160 @type vector: numpy array, len 3 161 @keyword atom_name: The atom name used to label the atom representing the head of the vector. 162 @type atom_name: str 163 @keyword res_name_vect: The 3 letter PDB residue code used to represent the vector. 164 @type res_name_vect: str 165 @keyword sim_vectors: The optional Monte Carlo simulation vectors to be represented in the PDB. 166 @type sim_vectors: list of numpy array, each len 3 167 @keyword res_name_sim: The 3 letter PDB residue code used to represent the Monte Carlo simulation vectors. 168 @type res_name_sim: str 169 @keyword chain_id: The chain identification code. 170 @type chain_id: str 171 @keyword res_num: The residue number. 172 @type res_num: int 173 @keyword origin: The origin for the axis. 174 @type origin: numpy array, len 3 175 @keyword scale: The scaling factor to stretch the vectors by. 176 @type scale: float 177 @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. 178 @type label_placement: float 179 @keyword neg: If True, then the negative vector positioned at the origin will also be included. 180 @type neg: bool 181 @return: The new residue number. 182 @rtype: int 183 """ 184 185 # The atom numbers (and indices). 186 origin_num = len(mol.atom_num)+1 187 atom_num = len(mol.atom_num)+2 188 atom_neg_num = len(mol.atom_num)+3 189 190 # The origin atom. 191 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') 192 193 # Create the PDB residue representing the vector. 194 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') 195 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 196 num = 1 197 if neg: 198 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') 199 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 200 num = 2 201 202 # Add another atom to allow the axis labels to be shifted just outside of the vector itself. 203 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') 204 if neg: 205 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') 206 207 # Print out. 208 print(" " + atom_name + " vector (scaled + shifted to origin): " + repr(origin+vector*scale)) 209 print(" Creating the MC simulation vectors.") 210 211 # Monte Carlo simulations. 212 if sim_vectors != None: 213 for i in range(len(sim_vectors)): 214 # Increment the residue number, so each simulation is a new residue. 215 res_num = res_num + 1 216 217 # The atom numbers (and indices). 218 atom_num = mol.atom_num[-1]+1 219 atom_neg_num = mol.atom_num[-1]+2 220 221 # Create the PDB residue representing the vector. 222 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') 223 mol.atom_connect(index1=atom_num-1, index2=origin_num-1) 224 if neg: 225 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') 226 mol.atom_connect(index1=atom_neg_num-1, index2=origin_num-1) 227 228 # Return the new residue number. 229 return res_num

230 231

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

233 """Create a distribution of vectors on a sphere using a distribution of spherical angles. 234 235 This function returns an array of unit vectors distributed within 3D space. The unit vectors are generated using the equation:: 236 237 | cos(theta) * sin(phi) | 238 vector = | sin(theta) * sin(phi) |. 239 | cos(phi) | 240 241 The vectors of this distribution generate both longitudinal and latitudinal lines. 242 243 244 @keyword inc: The number of increments in the distribution. 245 @type inc: int 246 @keyword distribution: The type of point distribution to use. This can be 'uniform' or 'regular'. 247 @type distribution: str 248 @return: The distribution of vectors on a sphere. 249 @rtype: list of rank-1, 3D numpy arrays 250 """ 251 252 # Get the polar and azimuthal angles for the distribution. 253 if distribution == 'uniform': 254 phi, theta = angles_uniform(inc) 255 else: 256 phi, theta = angles_regular(inc) 257 258 # Initialise array of the distribution of vectors. 259 vectors = [] 260 261 # Loop over the longitudinal lines. 262 for j in range(len(phi)): 263 # Loop over the latitudinal lines. 264 for i in range(len(theta)): 265 # X coordinate. 266 x = cos(theta[i]) * sin(phi[j]) 267 268 # Y coordinate. 269 y = sin(theta[i])* sin(phi[j]) 270 271 # Z coordinate. 272 z = cos(phi[j]) 273 274 # Append the vector. 275 vectors.append(array([x, y, z], float64)) 276 277 # Return the array of vectors and angles. 278 return vectors

279

Source Code for Module lib.structure.geometric