generic_fns.dipole

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2012 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 # Module docstring. 23 """Module for the manipulation of relaxation data.""" 24 25 # Python module imports. 26 from numpy import float64, zeros 27 from numpy.linalg import norm 28 import sys 29 from warnings import warn 30 31 # relax module imports. 32 from arg_check import is_float 33 from generic_fns.interatomic import create_interatom, exists_data, interatomic_loop, return_interatom 34 from generic_fns.mol_res_spin import Selection, exists_mol_res_spin_data, return_spin, spin_loop 35 from generic_fns import pipes 36 from relax_errors import RelaxError, RelaxNoInteratomError 37 from relax_io import extract_data, write_data 38 from relax_warnings import RelaxWarning, RelaxZeroVectorWarning 39 40

41 -def define(spin_id1=None, spin_id2=None, pipe=None, direct_bond=False, verbose=True):

42 """Set up the magnetic dipole-dipole interaction. 43 44 @keyword spin_id1: The spin identifier string of the first spin of the pair. 45 @type spin_id1: str 46 @keyword spin_id2: The spin identifier string of the second spin of the pair. 47 @type spin_id2: str 48 @param pipe: The data pipe to operate on. Defaults to the current data pipe. 49 @type pipe: str 50 @keyword direct_bond: A flag specifying if the two spins are directly bonded. 51 @type direct_bond: bool 52 @keyword verbose: A flag which if True will result in printouts of the created interatomoic data containers. 53 @type verbose: bool 54 """ 55 56 # The data pipe. 57 if pipe == None: 58 pipe = pipes.cdp_name() 59 60 # Get the data pipe. 61 dp = pipes.get_pipe(pipe) 62 63 # Loop over both spin selections. 64 ids = [] 65 for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True): 66 for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True): 67 # Directly bonded atoms. 68 if direct_bond: 69 # Different molecules. 70 if mol_name1 != mol_name2: 71 continue 72 73 # From structural info. 74 if hasattr(dp, 'structure') and dp.structure.get_molecule(mol_name1, model=1): 75 if not dp.structure.are_bonded(atom_id1=id1, atom_id2=id2): 76 continue 77 78 # From the residue info. 79 else: 80 # No element info. 81 if not hasattr(spin1, 'element'): 82 raise RelaxError("The spin '%s' does not have the element type set." % id1) 83 if not hasattr(spin2, 'element'): 84 raise RelaxError("The spin '%s' does not have the element type set." % id2) 85 86 # Backbone NH and CH pairs. 87 pair = False 88 if (spin1.element == 'N' and spin2.element == 'H') or (spin2.element == 'N' and spin1.element == 'H'): 89 pair = True 90 elif (spin1.element == 'C' and spin2.element == 'H') or (spin2.element == 'C' and spin1.element == 'H'): 91 pair = True 92 93 # Same residue, so skip. 94 if pair and res_num1 != None and res_num1 != res_num2: 95 continue 96 elif pair and res_num1 == None and res_name1 != res_name2: 97 continue 98 99 # Get the interatomic data object, if it exists. 100 interatom = return_interatom(id1, id2, pipe=pipe) 101 102 # Create the container if needed. 103 if interatom == None: 104 interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe) 105 106 # Check that this has not already been set up. 107 if interatom.dipole_pair: 108 raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2)) 109 110 # Set a flag indicating that a dipole-dipole interaction is present. 111 interatom.dipole_pair = True 112 113 # Store the IDs for the printout. 114 ids.append([repr(id1), repr(id2)]) 115 116 # No matches, so fail! 117 if not len(ids): 118 # Find the problem. 119 count1 = 0 120 count2 = 0 121 for spin in spin_loop(spin_id1): 122 count1 += 1 123 for spin in spin_loop(spin_id2): 124 count2 += 1 125 126 # Report the problem. 127 if count1 == 0 and count2 == 0: 128 raise RelaxError("Neither spin IDs '%s' and '%s' match any spins." % (spin_id1, spin_id2)) 129 elif count1 == 0: 130 raise RelaxError("The spin ID '%s' matches no spins." % spin_id1) 131 elif count2 == 0: 132 raise RelaxError("The spin ID '%s' matches no spins." % spin_id2) 133 else: 134 raise RelaxError("Unknown error.") 135 136 # Print out. 137 if verbose: 138 print("Magnetic dipole-dipole interactions are now defined for the following spins:\n") 139 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)

140 141

142 -def read_dist(file=None, dir=None, unit='meter', spin_id1_col=None, spin_id2_col=None, data_col=None, sep=None):

143 """Set up the magnetic dipole-dipole interaction. 144 145 @keyword file: The name of the file to open. 146 @type file: str 147 @keyword dir: The directory containing the file (defaults to the current directory if None). 148 @type dir: str or None 149 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 150 @type unit: str 151 @keyword spin_id1_col: The column containing the spin ID strings of the first spin. 152 @type spin_id1_col: int 153 @keyword spin_id2_col: The column containing the spin ID strings of the second spin. 154 @type spin_id2_col: int 155 @keyword data_col: The column containing the averaged distances in meters. 156 @type data_col: int or None 157 @keyword sep: The column separator which, if None, defaults to whitespace. 158 @type sep: str or None 159 """ 160 161 # Check the units. 162 if unit not in ['meter', 'Angstrom']: 163 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 164 165 # Test if the current data pipe exists. 166 pipes.test() 167 168 # Test if sequence data exists. 169 if not exists_mol_res_spin_data(): 170 raise RelaxNoSequenceError 171 172 # Extract the data from the file. 173 file_data = extract_data(file, dir, sep=sep) 174 175 # Loop over the RDC data. 176 data = [] 177 for line in file_data: 178 # Invalid columns. 179 if spin_id1_col > len(line): 180 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line)) 181 continue 182 if spin_id2_col > len(line): 183 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line)) 184 continue 185 if data_col and data_col > len(line): 186 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line)) 187 continue 188 189 # Unpack. 190 spin_id1 = line[spin_id1_col-1] 191 spin_id2 = line[spin_id2_col-1] 192 ave_dist = None 193 if data_col: 194 ave_dist = line[data_col-1] 195 196 # Convert and check the value. 197 if ave_dist != None: 198 try: 199 ave_dist = float(ave_dist) 200 except ValueError: 201 warn(RelaxWarning("The averaged distance of '%s' from the line %s is invalid." % (ave_dist, line))) 202 continue 203 204 # Unit conversion. 205 if unit == 'Angstrom': 206 ave_dist = ave_dist * 1e-10 207 208 # Get the interatomic data container. 209 interatom = return_interatom(spin_id1, spin_id2) 210 211 # No container found. 212 if interatom == None: 213 raise RelaxNoInteratomError(spin_id1=spin_id1, spin_id2=spin_id2) 214 215 # Store the averaged distance. 216 interatom.r = ave_dist 217 218 # Store the data for the printout. 219 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 220 221 # No data, so fail! 222 if not len(data): 223 raise RelaxError("No data could be extracted from the file.") 224 225 # Print out. 226 print("The following averaged distances have been read:\n") 227 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)

228 229

230 -def set_dist(spin_id1=None, spin_id2=None, ave_dist=None, unit='meter'):

231 """Set up the magnetic dipole-dipole interaction. 232 233 @keyword spin_id1: The spin identifier string of the first spin of the pair. 234 @type spin_id1: str 235 @keyword spin_id2: The spin identifier string of the second spin of the pair. 236 @type spin_id2: str 237 @keyword ave_dist: The r^-3 averaged interatomic distance. 238 @type ave_dist: float 239 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 240 @type unit: str 241 """ 242 243 # Check the units. 244 if unit not in ['meter', 'Angstrom']: 245 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 246 247 # Unit conversion. 248 if unit == 'Angstrom': 249 ave_dist = ave_dist * 1e-10 250 251 # Generate the selection objects. 252 sel_obj1 = Selection(spin_id1) 253 sel_obj2 = Selection(spin_id2) 254 255 # Loop over the interatomic containers. 256 data = [] 257 for interatom in interatomic_loop(): 258 # Get the spin info. 259 mol_name1, res_num1, res_name1, spin1 = return_spin(interatom.spin_id1, full_info=True) 260 mol_name2, res_num2, res_name2, spin2 = return_spin(interatom.spin_id2, full_info=True) 261 262 # No match, either way. 263 if not (sel_obj1.contains_spin(spin_num=spin1.num, spin_name=spin1.name, res_num=res_num1, res_name=res_name1, mol=mol_name1) and sel_obj2.contains_spin(spin_num=spin2.num, spin_name=spin2.name, res_num=res_num2, res_name=res_name2, mol=mol_name2)) and not (sel_obj2.contains_spin(spin_num=spin1.num, spin_name=spin1.name, res_num=res_num1, res_name=res_name1, mol=mol_name1) and sel_obj1.contains_spin(spin_num=spin2.num, spin_name=spin2.name, res_num=res_num2, res_name=res_name2, mol=mol_name2)): 264 continue 265 266 # Store the averaged distance. 267 interatom.r = ave_dist 268 269 # Store the data for the printout. 270 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 271 272 # No data, so fail! 273 if not len(data): 274 raise RelaxError("No data could be set.") 275 276 # Print out. 277 print("The following averaged distances have been set:\n") 278 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)

279 280

281 -def unit_vectors(ave=True):

282 """Extract the bond vectors from the loaded structures and store them in the spin container. 283 284 @keyword ave: A flag which if True will cause the average of all vectors to be calculated. 285 @type ave: bool 286 """ 287 288 # Test if the current data pipe exists. 289 pipes.test() 290 291 # Test if interatomic data exists. 292 if not exists_data(): 293 raise RelaxNoInteratomError 294 295 # Print out. 296 if ave: 297 print("Averaging all vectors.") 298 else: 299 print("No averaging of the vectors.") 300 301 # Loop over the interatomic data containers. 302 no_vectors = True 303 pos_info = False 304 for interatom in interatomic_loop(): 305 # Get the spin info. 306 spin1 = return_spin(interatom.spin_id1) 307 spin2 = return_spin(interatom.spin_id2) 308 309 # No positional information. 310 if not hasattr(spin1, 'pos'): 311 continue 312 if not hasattr(spin2, 'pos'): 313 continue 314 315 # Positional information flag. 316 pos_info = True 317 318 # Both single positions. 319 if is_float(spin1.pos[0], raise_error=False) and is_float(spin2.pos[0], raise_error=False): 320 # The vector. 321 vector_list = [spin2.pos - spin1.pos] 322 323 # A single and multiple position pair. 324 elif is_float(spin1.pos[0], raise_error=False) or is_float(spin2.pos[0], raise_error=False): 325 # The first spin has multiple positions. 326 if is_float(spin2.pos[0], raise_error=False): 327 vector_list = [] 328 for i in range(len(spin1.pos)): 329 vector_list.append(spin2.pos - spin1.pos[i]) 330 331 # The second spin has multiple positions. 332 else: 333 vector_list = [] 334 for i in range(len(spin2.pos)): 335 vector_list.append(spin2.pos[i] - spin1.pos) 336 337 # Both spins have multiple positions. 338 else: 339 # Non-matching number of positions. 340 if len(spin1.pos) != len(spin2.pos): 341 raise RelaxError("The spin '%s' consists of %s positions whereas the spin '%s' consists of %s - these numbers must match." % (interatom.spin_id1, len(spin1.pos), interatom.spin_id1, len(spin1.pos))) 342 343 # Calculate all vectors. 344 vector_list = [] 345 for i in range(len(spin1.pos)): 346 vector_list.append(spin2.pos[i] - spin1.pos[i]) 347 348 # Unit vectors. 349 for i in range(len(vector_list)): 350 # Normalisation factor. 351 norm_factor = norm(vector_list[i]) 352 353 # Test for zero length. 354 if norm_factor == 0.0: 355 warn(RelaxZeroVectorWarning(id)) 356 357 # Calculate the normalised vector. 358 else: 359 vector_list[i] = vector_list[i] / norm_factor 360 361 # Average. 362 if ave: 363 ave_vector = zeros(3, float64) 364 for i in range(len(vector_list)): 365 ave_vector = ave_vector + vector_list[i] 366 vector_list = [ave_vector / len(vector_list)] 367 368 # Convert to a single vector if needed. 369 if len(vector_list) == 1: 370 vector_list = vector_list[0] 371 372 # Store the unit vector(s). 373 setattr(interatom, 'vector', vector_list) 374 375 # We have a vector! 376 no_vectors = False 377 378 # Print out. 379 num = 1 380 if not is_float(vector_list[0], raise_error=False): 381 num = len(vector_list) 382 plural = 's' 383 if num == 1: 384 plural = '' 385 if spin1.name: 386 spin1_str = spin1.name 387 else: 388 spin1_str = spin1.num 389 if spin2.name: 390 spin2_str = spin2.name 391 else: 392 spin2_str = spin2.num 393 print("Calculated %s %s-%s unit vector%s between the spins '%s' and '%s'." % (num, spin1_str, spin2_str, plural, interatom.spin_id1, interatom.spin_id2)) 394 395 # Catch the problem of no positional information being present. 396 if not pos_info: 397 raise RelaxError("Positional information could not be found for any spins.") 398 399 # Right, catch the problem of missing vectors to prevent massive user confusion! 400 if no_vectors: 401 raise RelaxError("No vectors could be extracted.")

402

Source Code for Module generic_fns.dipole_pair