Package pipe_control :: Module interatomic
[hide private]
[frames] | no frames]

Source Code for Module pipe_control.interatomic

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2012-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  # Module docstring. 
 23  """Module for the manipulation of the interatomic data structures in the relax data store.""" 
 24   
 25  # Python module imports. 
 26  from copy import deepcopy 
 27  from numpy import float64, zeros 
 28  from numpy.linalg import norm 
 29  from re import search 
 30  import sys 
 31  from warnings import warn 
 32   
 33  # relax module imports. 
 34  from lib.arg_check import is_float 
 35  from lib.errors import RelaxError, RelaxInteratomInconsistentError, RelaxNoInteratomError, RelaxNoSpinError 
 36  from lib.io import extract_data, strip, write_data 
 37  from lib.warnings import RelaxWarning, RelaxZeroVectorWarning 
 38  from pipe_control import pipes 
 39  from pipe_control.mol_res_spin import Selection, count_spins, exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_loop 
 40   
 41   
42 -def copy(pipe_from=None, pipe_to=None, spin_id1=None, spin_id2=None, verbose=True):
43 """Copy the interatomic data from one data pipe to another. 44 45 @keyword pipe_from: The data pipe to copy the interatomic data from. This defaults to the current data pipe. 46 @type pipe_from: str 47 @keyword pipe_to: The data pipe to copy the interatomic data to. This defaults to the current data pipe. 48 @type pipe_to: str 49 @keyword spin_id1: The spin ID string of the first atom. 50 @type spin_id1: str 51 @keyword spin_id2: The spin ID string of the second atom. 52 @type spin_id2: str 53 @keyword verbose: A flag which if True will cause info about each spin pair to be printed out. 54 @type verbose: bool 55 """ 56 57 # Defaults. 58 if pipe_from == None and pipe_to == None: 59 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 60 elif pipe_from == None: 61 pipe_from = pipes.cdp_name() 62 elif pipe_to == None: 63 pipe_to = pipes.cdp_name() 64 65 # Test if the pipe_from and pipe_to data pipes exist. 66 pipes.test(pipe_from) 67 pipes.test(pipe_to) 68 69 # Check that the spin IDs exist. 70 if spin_id1: 71 if count_spins(selection=spin_id1, pipe=pipe_from, skip_desel=False) == 0: 72 raise RelaxNoSpinError(spin_id1, pipe_from) 73 if count_spins(selection=spin_id1, pipe=pipe_to, skip_desel=False) == 0: 74 raise RelaxNoSpinError(spin_id1, pipe_to) 75 if spin_id2: 76 if count_spins(selection=spin_id2, pipe=pipe_from, skip_desel=False) == 0: 77 raise RelaxNoSpinError(spin_id2, pipe_from) 78 if count_spins(selection=spin_id2, pipe=pipe_to, skip_desel=False) == 0: 79 raise RelaxNoSpinError(spin_id2, pipe_to) 80 81 # Check for the sequence data in the target pipe if no spin IDs are given. 82 if not spin_id1 and not spin_id2: 83 for spin, spin_id in spin_loop(pipe=pipe_from, return_id=True): 84 if not return_spin(spin_id, pipe=pipe_to): 85 raise RelaxNoSpinError(spin_id, pipe_to) 86 87 # Test if pipe_from contains interatomic data (skipping the rest of the function if it is missing). 88 if not exists_data(pipe_from): 89 return 90 91 # Loop over the interatomic data of the pipe_from data pipe. 92 ids = [] 93 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2, pipe=pipe_from): 94 # Create a new container. 95 new_interatom = create_interatom(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2, pipe=pipe_to) 96 97 # Duplicate all the objects of the container. 98 for name in dir(interatom): 99 # Skip special objects. 100 if search('^_', name): 101 continue 102 103 # Skip the spin IDs. 104 if name in ['spin_id1', 'spin_id2']: 105 continue 106 107 # Skip class methods. 108 if name in list(interatom.__class__.__dict__.keys()): 109 continue 110 111 # Duplicate all other objects. 112 obj = deepcopy(getattr(interatom, name)) 113 setattr(new_interatom, name, obj) 114 115 # Store the IDs for the printout. 116 ids.append([repr(interatom.spin_id1), repr(interatom.spin_id2)]) 117 118 # Print out. 119 if verbose: 120 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)
121 122
123 -def consistent_interatomic_data(pipe1=None, pipe2=None):
124 """Check that the interatomic data is consistent between two data pipes. 125 126 @keyword pipe1: The name of the first data pipe to compare. 127 @type pipe1: str 128 @keyword pipe2: The name of the second data pipe to compare. 129 @type pipe2: str 130 @raises RelaxError: If the data is inconsistent. 131 """ 132 133 # Get the data pipes. 134 dp1 = pipes.get_pipe(pipe1) 135 dp2 = pipes.get_pipe(pipe2) 136 137 # Check the data lengths. 138 if len(dp1.interatomic) != len(dp2.interatomic): 139 raise RelaxInteratomInconsistentError(pipe1, pipe2) 140 141 # Loop over the interatomic data. 142 for i in range(len(dp1.interatomic)): 143 # Alias the containers. 144 interatom1 = dp1.interatomic[i] 145 interatom2 = dp2.interatomic[i] 146 147 # Check the spin IDs. 148 if interatom1.spin_id1 != interatom2.spin_id1: 149 raise RelaxInteratomInconsistentError(pipe1, pipe2) 150 if interatom1.spin_id2 != interatom2.spin_id2: 151 raise RelaxInteratomInconsistentError(pipe1, pipe2)
152 153
154 -def create_interatom(spin_id1=None, spin_id2=None, pipe=None, verbose=False):
155 """Create and return the interatomic data container for the two spins. 156 157 @keyword spin_id1: The spin ID string of the first atom. 158 @type spin_id1: str 159 @keyword spin_id2: The spin ID string of the second atom. 160 @type spin_id2: str 161 @keyword pipe: The data pipe to create the interatomic data container for. This defaults to the current data pipe if not supplied. 162 @type pipe: str or None 163 @keyword verbose: A flag which if True will result printouts. 164 @type verbose: bool 165 @return: The newly created interatomic data container. 166 @rtype: data.interatomic.InteratomContainer instance 167 """ 168 169 # Printout. 170 if verbose: 171 print("Creating an interatomic data container between the spins '%s' and '%s'." % (spin_id1, spin_id2)) 172 173 # The data pipe. 174 if pipe == None: 175 pipe = pipes.cdp_name() 176 177 # Get the data pipe. 178 dp = pipes.get_pipe(pipe) 179 180 # Check that the spin IDs exist. 181 spin = return_spin(spin_id1, pipe) 182 if spin == None: 183 raise RelaxNoSpinError(spin_id1) 184 spin = return_spin(spin_id2, pipe) 185 if spin == None: 186 raise RelaxNoSpinError(spin_id2) 187 188 # Check if the two spin IDs have already been added. 189 for i in range(len(dp.interatomic)): 190 if id_match(spin_id=spin_id1, interatom=dp.interatomic[i], pipe=pipe) and id_match(spin_id=spin_id2, interatom=dp.interatomic[i], pipe=pipe): 191 raise RelaxError("The spin pair %s and %s have already been added." % (spin_id1, spin_id2)) 192 193 # Add and return the data. 194 return dp.interatomic.add_item(spin_id1=spin_id1, spin_id2=spin_id2)
195 196
197 -def define(spin_id1=None, spin_id2=None, pipe=None, direct_bond=False, verbose=True):
198 """Set up the magnetic dipole-dipole interaction. 199 200 @keyword spin_id1: The spin identifier string of the first spin of the pair. 201 @type spin_id1: str 202 @keyword spin_id2: The spin identifier string of the second spin of the pair. 203 @type spin_id2: str 204 @param pipe: The data pipe to operate on. Defaults to the current data pipe. 205 @type pipe: str 206 @keyword direct_bond: A flag specifying if the two spins are directly bonded. 207 @type direct_bond: bool 208 @keyword verbose: A flag which if True will result in printouts of the created interatomoic data containers. 209 @type verbose: bool 210 """ 211 212 # The data pipe. 213 if pipe == None: 214 pipe = pipes.cdp_name() 215 216 # Get the data pipe. 217 dp = pipes.get_pipe(pipe) 218 219 # Initialise the spin ID pairs list. 220 ids = [] 221 222 # Use the structural data to find connected atoms. 223 if hasattr(dp, 'structure'): 224 # Loop over the atoms of the first spin selection. 225 for mol_name1, res_num1, res_name1, atom_num1, atom_name1, mol_index1, atom_index1 in dp.structure.atom_loop(atom_id=spin_id1, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True): 226 # Generate the first spin ID. 227 id1 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name1, res_num=res_num1, res_name=res_name1, spin_num=atom_num1, spin_name=atom_name1) 228 229 # Do the spin exist? 230 if not return_spin(id1): 231 continue 232 233 # Loop over the atoms of the second spin selection. 234 for mol_name2, res_num2, res_name2, atom_num2, atom_name2, mol_index2, atom_index2 in dp.structure.atom_loop(atom_id=spin_id2, model_num=1, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, mol_index_flag=True, index_flag=True): 235 # Directly bonded atoms. 236 if direct_bond: 237 # Different molecules. 238 if mol_name1 != mol_name2: 239 continue 240 241 # Skip non-bonded atom pairs. 242 if not dp.structure.are_bonded_index(mol_index1=mol_index1, atom_index1=atom_index1, mol_index2=mol_index2, atom_index2=atom_index2): 243 continue 244 245 # Generate the second spin ID. 246 id2 = generate_spin_id_unique(pipe_cont=dp, mol_name=mol_name2, res_num=res_num2, res_name=res_name2, spin_num=atom_num2, spin_name=atom_name2) 247 248 # Do the spin exist? 249 if not return_spin(id2): 250 continue 251 252 # Store the IDs for the printout. 253 ids.append([id1, id2]) 254 255 # No structural data present or the spin IDs are not in the structural data, so use spin loops and some basic rules. 256 if ids == []: 257 for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True): 258 for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True): 259 # Directly bonded atoms. 260 if direct_bond: 261 # Different molecules. 262 if mol_name1 != mol_name2: 263 continue 264 265 # No element info. 266 if not hasattr(spin1, 'element'): 267 raise RelaxError("The spin '%s' does not have the element type set." % id1) 268 if not hasattr(spin2, 'element'): 269 raise RelaxError("The spin '%s' does not have the element type set." % id2) 270 271 # Backbone NH and CH pairs. 272 pair = False 273 if (spin1.element == 'N' and spin2.element == 'H') or (spin2.element == 'N' and spin1.element == 'H'): 274 pair = True 275 elif (spin1.element == 'C' and spin2.element == 'H') or (spin2.element == 'C' and spin1.element == 'H'): 276 pair = True 277 278 # Same residue, so skip. 279 if pair and res_num1 != None and res_num1 != res_num2: 280 continue 281 elif pair and res_num1 == None and res_name1 != res_name2: 282 continue 283 284 # Store the IDs for the printout. 285 ids.append([id1, id2]) 286 287 # No matches, so fail! 288 if not len(ids): 289 # Find the problem. 290 count1 = 0 291 count2 = 0 292 for spin in spin_loop(spin_id1): 293 count1 += 1 294 for spin in spin_loop(spin_id2): 295 count2 += 1 296 297 # Report the problem. 298 if count1 == 0 and count2 == 0: 299 raise RelaxError("Neither spin IDs '%s' and '%s' match any spins." % (spin_id1, spin_id2)) 300 elif count1 == 0: 301 raise RelaxError("The spin ID '%s' matches no spins." % spin_id1) 302 elif count2 == 0: 303 raise RelaxError("The spin ID '%s' matches no spins." % spin_id2) 304 else: 305 raise RelaxError("Unknown error.") 306 307 # Define the interaction. 308 for id1, id2 in ids: 309 # Get the interatomic data object, if it exists. 310 interatom = return_interatom(id1, id2, pipe=pipe) 311 312 # Create the container if needed. 313 if interatom == None: 314 interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe) 315 316 # Check that this has not already been set up. 317 if interatom.dipole_pair: 318 raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2)) 319 320 # Set a flag indicating that a dipole-dipole interaction is present. 321 interatom.dipole_pair = True 322 323 # Printout. 324 if verbose: 325 # Conversion. 326 for i in range(len(ids)): 327 ids[i][0] = repr(ids[i][0]) 328 ids[i][1] = repr(ids[i][1]) 329 330 # The printout. 331 print("Interatomic interactions are now defined for the following spins:\n") 332 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)
333 334
335 -def exists_data(pipe=None):
336 """Determine if any interatomic data exists. 337 338 @keyword pipe: The data pipe in which the interatomic data will be checked for. 339 @type pipe: str 340 @return: The answer to the question about the existence of data. 341 @rtype: bool 342 """ 343 344 # The current data pipe. 345 if pipe == None: 346 pipe = pipes.cdp_name() 347 348 # Test the data pipe. 349 pipes.test(pipe) 350 351 # Get the data pipe. 352 dp = pipes.get_pipe(pipe) 353 354 # The interatomic data structure is empty. 355 if dp.interatomic.is_empty(): 356 return False 357 358 # Otherwise. 359 return True
360 361
362 -def id_match(spin_id=None, interatom=None, pipe=None):
363 """Test if the spin ID matches one of the two spins of the given container. 364 365 @keyword spin_id: The spin ID string of the first atom. 366 @type spin_id: str 367 @keyword interatom: The interatomic data container. 368 @type interatom: InteratomContainer instance 369 @keyword pipe: The data pipe containing the interatomic data container. Defaults to the current data pipe. 370 @type pipe: str or None 371 @return: True if the spin ID matches one of the two spins, False otherwise. 372 @rtype: bool 373 """ 374 375 # Get the spin containers. 376 spin1 = return_spin(interatom.spin_id1, pipe=pipe) 377 spin2 = return_spin(interatom.spin_id2, pipe=pipe) 378 379 # No spins. 380 if spin1 == None or spin2 == None: 381 return False 382 383 # Check if the ID is in the private metadata list. 384 if spin_id in spin1._spin_ids or spin_id in spin2._spin_ids: 385 return True 386 387 # Nothing found. 388 return False
389 390
391 -def interatomic_loop(selection1=None, selection2=None, pipe=None, skip_desel=True):
392 """Generator function for looping over all the interatomic data containers. 393 394 @keyword selection1: The optional spin ID selection of the first atom. 395 @type selection1: str 396 @keyword selection2: The optional spin ID selection of the second atom. 397 @type selection2: str 398 @keyword pipe: The data pipe containing the spin. Defaults to the current data pipe. 399 @type pipe: str 400 @keyword skip_desel: A flag which if True will cause only selected interatomic data containers to be returned. 401 @type skip_desel: bool 402 """ 403 404 # The data pipe. 405 if pipe == None: 406 pipe = pipes.cdp_name() 407 408 # Get the data pipe. 409 dp = pipes.get_pipe(pipe) 410 411 # Parse the spin ID selection strings. 412 select_obj = None 413 select_obj1 = None 414 select_obj2 = None 415 if selection1 and selection2: 416 select_obj1 = Selection(selection1) 417 select_obj2 = Selection(selection2) 418 elif selection1: 419 select_obj = Selection(selection1) 420 elif selection2: 421 select_obj = Selection(selection2) 422 423 # Loop over the containers, yielding them. 424 for i in range(len(dp.interatomic)): 425 # Skip deselected containers. 426 if skip_desel and not dp.interatomic[i].select: 427 continue 428 429 # Aliases. 430 interatom = dp.interatomic[i] 431 mol_index1, res_index1, spin_index1 = cdp.mol._spin_id_lookup[interatom.spin_id1] 432 mol_index2, res_index2, spin_index2 = cdp.mol._spin_id_lookup[interatom.spin_id2] 433 mol1 = cdp.mol[mol_index1] 434 res1 = cdp.mol[mol_index1].res[res_index1] 435 spin1 = cdp.mol[mol_index1].res[res_index1].spin[spin_index1] 436 mol2 = cdp.mol[mol_index2] 437 res2 = cdp.mol[mol_index2].res[res_index2] 438 spin2 = cdp.mol[mol_index2].res[res_index2].spin[spin_index2] 439 440 # The different selection combinations. 441 if select_obj: 442 sel1 = select_obj.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 443 sel2 = select_obj.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 444 if select_obj1: 445 sel11 = select_obj1.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 446 sel12 = select_obj1.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 447 if select_obj2: 448 sel21 = select_obj2.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 449 sel22 = select_obj2.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 450 451 # Check that the selections are met. 452 if select_obj: 453 if not sel1 and not sel2: 454 continue 455 if select_obj1: 456 if not (sel11 or sel12): 457 continue 458 if select_obj2: 459 if not (sel21 or sel22): 460 continue 461 462 # Return the container. 463 yield interatom
464 465
466 -def read_dist(file=None, dir=None, unit='meter', spin_id1_col=None, spin_id2_col=None, data_col=None, sep=None):
467 """Set up the magnetic dipole-dipole interaction. 468 469 @keyword file: The name of the file to open. 470 @type file: str 471 @keyword dir: The directory containing the file (defaults to the current directory if None). 472 @type dir: str or None 473 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 474 @type unit: str 475 @keyword spin_id1_col: The column containing the spin ID strings of the first spin. 476 @type spin_id1_col: int 477 @keyword spin_id2_col: The column containing the spin ID strings of the second spin. 478 @type spin_id2_col: int 479 @keyword data_col: The column containing the averaged distances in meters. 480 @type data_col: int or None 481 @keyword sep: The column separator which, if None, defaults to whitespace. 482 @type sep: str or None 483 """ 484 485 # Check the units. 486 if unit not in ['meter', 'Angstrom']: 487 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 488 489 # Test if the current data pipe exists. 490 pipes.test() 491 492 # Test if sequence data exists. 493 if not exists_mol_res_spin_data(): 494 raise RelaxNoSequenceError 495 496 # Extract the data from the file, and clean it up. 497 file_data = extract_data(file, dir, sep=sep) 498 file_data = strip(file_data, comments=True) 499 500 # Loop over the RDC data. 501 data = [] 502 for line in file_data: 503 # Invalid columns. 504 if spin_id1_col > len(line): 505 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line)) 506 continue 507 if spin_id2_col > len(line): 508 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line)) 509 continue 510 if data_col and data_col > len(line): 511 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line)) 512 continue 513 514 # Unpack. 515 spin_id1 = line[spin_id1_col-1] 516 spin_id2 = line[spin_id2_col-1] 517 ave_dist = None 518 if data_col: 519 ave_dist = line[data_col-1] 520 521 # Convert and check the value. 522 if ave_dist != None: 523 try: 524 ave_dist = float(ave_dist) 525 except ValueError: 526 warn(RelaxWarning("The averaged distance of '%s' from the line %s is invalid." % (ave_dist, line))) 527 continue 528 529 # Unit conversion. 530 if unit == 'Angstrom': 531 ave_dist = ave_dist * 1e-10 532 533 # Get the interatomic data container. 534 interatom = return_interatom(spin_id1, spin_id2) 535 536 # No container found, so create it. 537 if interatom == None: 538 interatom = create_interatom(spin_id1=spin_id1, spin_id2=spin_id2, verbose=True) 539 540 # Store the averaged distance. 541 interatom.r = ave_dist 542 543 # Store the data for the printout. 544 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 545 546 # No data, so fail! 547 if not len(data): 548 raise RelaxError("No data could be extracted from the file.") 549 550 # Print out. 551 print("The following averaged distances have been read:\n") 552 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)
553 554
555 -def return_interatom(spin_id1=None, spin_id2=None, pipe=None):
556 """Return the list of interatomic data containers for the two spins. 557 558 @keyword spin_id1: The spin ID string of the first atom. 559 @type spin_id1: str 560 @keyword spin_id2: The spin ID string of the second atom. 561 @type spin_id2: str 562 @keyword pipe: The data pipe holding the container. Defaults to the current data pipe. 563 @type pipe: str or None 564 @return: The matching interatomic data container, if it exists. 565 @rtype: data.interatomic.InteratomContainer instance or None 566 """ 567 568 # Checks. 569 if spin_id1 == None: 570 raise RelaxError("The first spin ID must be supplied.") 571 if spin_id2 == None: 572 raise RelaxError("The second spin ID must be supplied.") 573 574 # The data pipe. 575 if pipe == None: 576 pipe = pipes.cdp_name() 577 578 # Get the data pipe. 579 dp = pipes.get_pipe(pipe) 580 581 # Return the matching container. 582 for i in range(len(dp.interatomic)): 583 if id_match(spin_id=spin_id1, interatom=dp.interatomic[i], pipe=pipe) and id_match(spin_id=spin_id2, interatom=dp.interatomic[i], pipe=pipe): 584 return dp.interatomic[i] 585 586 # No matchs. 587 return None
588 589
590 -def return_interatom_list(spin_id=None, pipe=None):
591 """Return the list of interatomic data containers for the given spin. 592 593 @keyword spin_id: The spin ID string. 594 @type spin_id: str 595 @keyword pipe: The data pipe holding the container. This defaults to the current data pipe. 596 @type pipe: str or None 597 @return: The list of matching interatomic data containers, if any exist. 598 @rtype: list of data.interatomic.InteratomContainer instances 599 """ 600 601 # Check. 602 if spin_id == None: 603 raise RelaxError("The spin ID must be supplied.") 604 605 # The data pipe. 606 if pipe == None: 607 pipe = pipes.cdp_name() 608 609 # Get the data pipe. 610 dp = pipes.get_pipe(pipe) 611 612 # Initialise. 613 interatoms = [] 614 615 # Find and append all containers. 616 for i in range(len(dp.interatomic)): 617 if id_match(spin_id=spin_id, interatom=dp.interatomic[i], pipe=pipe) or id_match(spin_id=spin_id, interatom=dp.interatomic[i], pipe=pipe): 618 interatoms.append(dp.interatomic[i]) 619 620 # Return the list of containers. 621 return interatoms
622 623
624 -def set_dist(spin_id1=None, spin_id2=None, ave_dist=None, unit='meter'):
625 """Set up the magnetic dipole-dipole interaction. 626 627 @keyword spin_id1: The spin identifier string of the first spin of the pair. 628 @type spin_id1: str 629 @keyword spin_id2: The spin identifier string of the second spin of the pair. 630 @type spin_id2: str 631 @keyword ave_dist: The r^-3 averaged interatomic distance. 632 @type ave_dist: float 633 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 634 @type unit: str 635 """ 636 637 # Check the units. 638 if unit not in ['meter', 'Angstrom']: 639 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 640 641 # Unit conversion. 642 if unit == 'Angstrom': 643 ave_dist = ave_dist * 1e-10 644 645 # Generate the selection objects. 646 sel_obj1 = Selection(spin_id1) 647 sel_obj2 = Selection(spin_id2) 648 649 # Loop over the interatomic containers. 650 data = [] 651 for interatom in interatomic_loop(): 652 # Get the spin info. 653 mol_name1, res_num1, res_name1, spin1 = return_spin(interatom.spin_id1, full_info=True) 654 mol_name2, res_num2, res_name2, spin2 = return_spin(interatom.spin_id2, full_info=True) 655 656 # No match, either way. 657 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)): 658 continue 659 660 # Store the averaged distance. 661 interatom.r = ave_dist 662 663 # Store the data for the printout. 664 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 665 666 # No data, so fail! 667 if not len(data): 668 raise RelaxError("No data could be set.") 669 670 # Print out. 671 print("The following averaged distances have been set:\n") 672 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)
673 674
675 -def unit_vectors(ave=True):
676 """Extract the bond vectors from the loaded structures and store them in the spin container. 677 678 @keyword ave: A flag which if True will cause the average of all vectors to be calculated. 679 @type ave: bool 680 """ 681 682 # Test if the current data pipe exists. 683 pipes.test() 684 685 # Test if interatomic data exists. 686 if not exists_data(): 687 raise RelaxNoInteratomError 688 689 # Print out. 690 if ave: 691 print("Averaging all vectors.") 692 else: 693 print("No averaging of the vectors.") 694 695 # Loop over the interatomic data containers. 696 no_vectors = True 697 pos_info = False 698 for interatom in interatomic_loop(skip_desel=False): 699 # Get the spin info. 700 spin1 = return_spin(interatom.spin_id1) 701 spin2 = return_spin(interatom.spin_id2) 702 703 # No positional information. 704 if not hasattr(spin1, 'pos'): 705 continue 706 if not hasattr(spin2, 'pos'): 707 continue 708 709 # Positional information flag. 710 pos_info = True 711 712 # Both single positions. 713 if is_float(spin1.pos[0], raise_error=False) and is_float(spin2.pos[0], raise_error=False): 714 # The vector. 715 vector_list = [spin2.pos - spin1.pos] 716 717 # A single and multiple position pair. 718 elif is_float(spin1.pos[0], raise_error=False) or is_float(spin2.pos[0], raise_error=False): 719 # The first spin has multiple positions. 720 if is_float(spin2.pos[0], raise_error=False): 721 vector_list = [] 722 for i in range(len(spin1.pos)): 723 vector_list.append(spin2.pos - spin1.pos[i]) 724 725 # The second spin has multiple positions. 726 else: 727 vector_list = [] 728 for i in range(len(spin2.pos)): 729 vector_list.append(spin2.pos[i] - spin1.pos) 730 731 # Both spins have multiple positions. 732 else: 733 # Non-matching number of positions. 734 if len(spin1.pos) != len(spin2.pos): 735 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))) 736 737 # Calculate all vectors. 738 vector_list = [] 739 for i in range(len(spin1.pos)): 740 vector_list.append(spin2.pos[i] - spin1.pos[i]) 741 742 # Unit vectors. 743 for i in range(len(vector_list)): 744 # Normalisation factor. 745 norm_factor = norm(vector_list[i]) 746 747 # Test for zero length. 748 if norm_factor == 0.0: 749 warn(RelaxZeroVectorWarning(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2)) 750 751 # Calculate the normalised vector. 752 else: 753 vector_list[i] = vector_list[i] / norm_factor 754 755 # Average. 756 if ave: 757 ave_vector = zeros(3, float64) 758 for i in range(len(vector_list)): 759 ave_vector = ave_vector + vector_list[i] 760 vector_list = [ave_vector / len(vector_list)] 761 762 # Convert to a single vector if needed. 763 if len(vector_list) == 1: 764 vector_list = vector_list[0] 765 766 # Store the unit vector(s). 767 setattr(interatom, 'vector', vector_list) 768 769 # We have a vector! 770 no_vectors = False 771 772 # Print out. 773 num = 1 774 if not is_float(vector_list[0], raise_error=False): 775 num = len(vector_list) 776 plural = 's' 777 if num == 1: 778 plural = '' 779 if spin1.name: 780 spin1_str = spin1.name 781 else: 782 spin1_str = spin1.num 783 if spin2.name: 784 spin2_str = spin2.name 785 else: 786 spin2_str = spin2.num 787 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)) 788 789 # Catch the problem of no positional information being present. 790 if not pos_info: 791 raise RelaxError("Positional information could not be found for any spins.") 792 793 # Right, catch the problem of missing vectors to prevent massive user confusion! 794 if no_vectors: 795 raise RelaxError("No vectors could be extracted.")
796