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, sin 24 from numpy import 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), phi_max_fn=None, scale=1.0, inc=20, distribution='uniform'):

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: numpy array 51 @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. 52 @type phi_max_fn: 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_atom = zeros(len(theta)) 82 83 # 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. 84 if phi_max_fn: 85 phi_max = zeros(len(theta), float64) 86 edge_vectors = zeros((len(theta), 3), float64) 87 for i in range(len(theta)): 88 # The maximum angle. 89 phi_max[i] = phi_max_fn(theta[i]) 90 91 # The vector in the unrotated frame. 92 edge_vectors[i, 0] = cos(theta[i]) * sin(phi_max[i]) 93 edge_vectors[i, 1] = sin(theta[i])* sin(phi_max[i]) 94 edge_vectors[i, 2] = cos(phi_max[i]) 95 96 # Loop over the radial array of vectors (change in longitude). 97 for i in range(len(theta)): 98 # Loop over the vectors of the radial array (change in latitude). 99 for j in range(len(phi)): 100 # The vector to use. 101 if phi_max_fn and phi[j] > phi_max_fn(theta[i]): 102 vector = edge_vectors[i] 103 else: 104 vector = vectors[i + j*len(theta)] 105 106 # Update the edge for this longitude. 107 if not edge[i]: 108 edge[i] = 1 109 edge_index[i] = j 110 edge_atom[i] = atom_num 111 112 # Rotate the vector into the frame. 113 vector = dot(R, vector) 114 115 # Warp the vector. 116 vector = dot(warp, vector) 117 118 # Scale the vector. 119 vector = vector * scale 120 121 # Position relative to the centre of mass. 122 pos = centre + vector 123 124 # Add the vector as a H atom of the TNS residue. 125 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') 126 127 # Connect to the previous atom to generate the longitudinal lines (except for the first point). 128 if j > edge_index[i]: 129 mol.atom_connect(index1=atom_num-1, index2=atom_num-2) 130 131 # Connect across the radial arrays (to generate the latitudinal lines). 132 if i != 0 and j >= edge_index[i-1]: 133 # The number of atoms back to the previous longitude. 134 num = len(phi) - edge_index[i] 135 136 # Latitude line. 137 mol.atom_connect(index1=atom_num-1, index2=atom_num-1-num) 138 139 # Connect the last radial array to the first (to zip up the geometric object and close the latitudinal lines). 140 if i == len(theta)-1 and j >= edge_index[0]: 141 # The number of atom in the first longitude line. 142 num = origin_num + j - edge_index[0] 143 144 # Zipping up. 145 mol.atom_connect(index1=atom_num-1, index2=num-1) 146 147 # Increment the atom number. 148 atom_num = atom_num + 1

149 150

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

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

229 230

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

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

278

Source Code for Module lib.structure.geometric