pipe_control.interatomic

1 ############################################################################### 2 # # 3 # Copyright (C) 2012-2015 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 from pipe_control.pipes import check_pipe 41 42

43 -def copy(pipe_from=None, pipe_to=None, spin_id1=None, spin_id2=None, verbose=True):

44 """Copy the interatomic data from one data pipe to another. 45 46 @keyword pipe_from: The data pipe to copy the interatomic data from. This defaults to the current data pipe. 47 @type pipe_from: str 48 @keyword pipe_to: The data pipe to copy the interatomic data to. This defaults to the current data pipe. 49 @type pipe_to: str 50 @keyword spin_id1: The spin ID string of the first atom. 51 @type spin_id1: str 52 @keyword spin_id2: The spin ID string of the second atom. 53 @type spin_id2: str 54 @keyword verbose: A flag which if True will cause info about each spin pair to be printed out. 55 @type verbose: bool 56 """ 57 58 # Defaults. 59 if pipe_from == None and pipe_to == None: 60 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 61 elif pipe_from == None: 62 pipe_from = pipes.cdp_name() 63 elif pipe_to == None: 64 pipe_to = pipes.cdp_name() 65 66 # Test if the pipe_from and pipe_to data pipes exist. 67 check_pipe(pipe_from) 68 check_pipe(pipe_to) 69 70 # Check that the spin IDs exist. 71 if spin_id1: 72 if count_spins(selection=spin_id1, pipe=pipe_from, skip_desel=False) == 0: 73 raise RelaxNoSpinError(spin_id1, pipe_from) 74 if count_spins(selection=spin_id1, pipe=pipe_to, skip_desel=False) == 0: 75 raise RelaxNoSpinError(spin_id1, pipe_to) 76 if spin_id2: 77 if count_spins(selection=spin_id2, pipe=pipe_from, skip_desel=False) == 0: 78 raise RelaxNoSpinError(spin_id2, pipe_from) 79 if count_spins(selection=spin_id2, pipe=pipe_to, skip_desel=False) == 0: 80 raise RelaxNoSpinError(spin_id2, pipe_to) 81 82 # Check for the sequence data in the target pipe if no spin IDs are given. 83 if not spin_id1 and not spin_id2: 84 for spin, spin_id in spin_loop(pipe=pipe_from, return_id=True): 85 if not return_spin(spin_id, pipe=pipe_to): 86 raise RelaxNoSpinError(spin_id, pipe_to) 87 88 # Test if pipe_from contains interatomic data (skipping the rest of the function if it is missing). 89 if not exists_data(pipe_from): 90 return 91 92 # Loop over the interatomic data of the pipe_from data pipe. 93 ids = [] 94 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2, pipe=pipe_from): 95 # Create a new container. 96 new_interatom = create_interatom(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2, pipe=pipe_to) 97 98 # Duplicate all the objects of the container. 99 for name in dir(interatom): 100 # Skip special objects. 101 if search('^_', name): 102 continue 103 104 # Skip the spin IDs. 105 if name in ['spin_id1', 'spin_id2']: 106 continue 107 108 # Skip class methods. 109 if name in interatom.__class__.__dict__: 110 continue 111 112 # Duplicate all other objects. 113 obj = deepcopy(getattr(interatom, name)) 114 setattr(new_interatom, name, obj) 115 116 # Store the IDs for the printout. 117 ids.append([repr(interatom.spin_id1), repr(interatom.spin_id2)]) 118 119 # Print out. 120 if verbose: 121 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)

122 123

124 -def consistent_interatomic_data(pipe1=None, pipe2=None):

125 """Check that the interatomic data is consistent between two data pipes. 126 127 @keyword pipe1: The name of the first data pipe to compare. 128 @type pipe1: str 129 @keyword pipe2: The name of the second data pipe to compare. 130 @type pipe2: str 131 @raises RelaxError: If the data is inconsistent. 132 """ 133 134 # Get the data pipes. 135 dp1 = pipes.get_pipe(pipe1) 136 dp2 = pipes.get_pipe(pipe2) 137 138 # Check the data lengths. 139 if len(dp1.interatomic) != len(dp2.interatomic): 140 raise RelaxInteratomInconsistentError(pipe1, pipe2) 141 142 # Loop over the interatomic data. 143 for i in range(len(dp1.interatomic)): 144 # Alias the containers. 145 interatom1 = dp1.interatomic[i] 146 interatom2 = dp2.interatomic[i] 147 148 # Check the spin IDs. 149 if interatom1.spin_id1 != interatom2.spin_id1: 150 raise RelaxInteratomInconsistentError(pipe1, pipe2) 151 if interatom1.spin_id2 != interatom2.spin_id2: 152 raise RelaxInteratomInconsistentError(pipe1, pipe2)

153 154

155 -def create_interatom(spin_id1=None, spin_id2=None, pipe=None, verbose=False):

156 """Create and return the interatomic data container for the two spins. 157 158 @keyword spin_id1: The spin ID string of the first atom. 159 @type spin_id1: str 160 @keyword spin_id2: The spin ID string of the second atom. 161 @type spin_id2: str 162 @keyword pipe: The data pipe to create the interatomic data container for. This defaults to the current data pipe if not supplied. 163 @type pipe: str or None 164 @keyword verbose: A flag which if True will result printouts. 165 @type verbose: bool 166 @return: The newly created interatomic data container. 167 @rtype: data.interatomic.InteratomContainer instance 168 """ 169 170 # Printout. 171 if verbose: 172 print("Creating an interatomic data container between the spins '%s' and '%s'." % (spin_id1, spin_id2)) 173 174 # The data pipe. 175 if pipe == None: 176 pipe = pipes.cdp_name() 177 178 # Get the data pipe. 179 dp = pipes.get_pipe(pipe) 180 181 # Check that the spin IDs exist. 182 spin = return_spin(spin_id1, pipe) 183 if spin == None: 184 raise RelaxNoSpinError(spin_id1) 185 spin = return_spin(spin_id2, pipe) 186 if spin == None: 187 raise RelaxNoSpinError(spin_id2) 188 189 # Check if the two spin IDs have already been added. 190 for i in range(len(dp.interatomic)): 191 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): 192 raise RelaxError("The spin pair %s and %s have already been added." % (spin_id1, spin_id2)) 193 194 # Add and return the data. 195 return dp.interatomic.add_item(spin_id1=spin_id1, spin_id2=spin_id2)

196 197

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

199 """Set up the magnetic dipole-dipole interaction. 200 201 @keyword spin_id1: The spin identifier string of the first spin of the pair. 202 @type spin_id1: str 203 @keyword spin_id2: The spin identifier string of the second spin of the pair. 204 @type spin_id2: str 205 @param pipe: The data pipe to operate on. Defaults to the current data pipe. 206 @type pipe: str 207 @keyword direct_bond: A flag specifying if the two spins are directly bonded. 208 @type direct_bond: bool 209 @keyword spin_selection: Define the interatomic data container selection based on the spin selection. If either spin is deselected, the interatomic container will also be deselected. Otherwise the container will be selected. 210 @type spin_selection: bool 211 @keyword verbose: A flag which if True will result in printouts of the created interatomoic data containers. 212 @type verbose: bool 213 """ 214 215 # The data pipe. 216 if pipe == None: 217 pipe = pipes.cdp_name() 218 219 # Get the data pipe. 220 dp = pipes.get_pipe(pipe) 221 222 # Initialise the spin ID pairs list. 223 ids = [] 224 spin_selections = [] 225 226 # Use the structural data to find connected atoms. 227 if hasattr(dp, 'structure'): 228 # The selection objects. 229 selection1 = cdp.structure.selection(atom_id=spin_id1) 230 selection2 = cdp.structure.selection(atom_id=spin_id2) 231 232 # Loop over the atoms of the first spin selection. 233 for mol_name1, res_num1, res_name1, atom_num1, atom_name1, mol_index1, atom_index1 in dp.structure.atom_loop(selection=selection1, 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): 234 # Generate the first spin ID. 235 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) 236 237 # Do the spin exist? 238 spin1 = return_spin(id1) 239 if not spin1: 240 continue 241 242 # Loop over the atoms of the second spin selection. 243 for mol_name2, res_num2, res_name2, atom_num2, atom_name2, mol_index2, atom_index2 in dp.structure.atom_loop(selection=selection2, 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): 244 # Directly bonded atoms. 245 if direct_bond: 246 # Different molecules. 247 if mol_name1 != mol_name2: 248 continue 249 250 # Skip non-bonded atom pairs. 251 if not dp.structure.are_bonded_index(mol_index1=mol_index1, atom_index1=atom_index1, mol_index2=mol_index2, atom_index2=atom_index2): 252 continue 253 254 # Generate the second spin ID. 255 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) 256 257 # Do the spin exist? 258 spin2 = return_spin(id2) 259 if not spin2: 260 continue 261 262 # Store the IDs for the printout. 263 ids.append([id1, id2]) 264 spin_selections.append([spin1.select, spin2.select]) 265 266 # No structural data present or the spin IDs are not in the structural data, so use spin loops and some basic rules. 267 if ids == []: 268 for spin1, mol_name1, res_num1, res_name1, id1 in spin_loop(spin_id1, pipe=pipe, full_info=True, return_id=True): 269 for spin2, mol_name2, res_num2, res_name2, id2 in spin_loop(spin_id2, pipe=pipe, full_info=True, return_id=True): 270 # Directly bonded atoms. 271 if direct_bond: 272 # Different molecules. 273 if mol_name1 != mol_name2: 274 continue 275 276 # No element info. 277 if not hasattr(spin1, 'element'): 278 raise RelaxError("The spin '%s' does not have the element type set." % id1) 279 if not hasattr(spin2, 'element'): 280 raise RelaxError("The spin '%s' does not have the element type set." % id2) 281 282 # Backbone NH and CH pairs. 283 pair = False 284 if (spin1.element == 'N' and spin2.element == 'H') or (spin2.element == 'N' and spin1.element == 'H'): 285 pair = True 286 elif (spin1.element == 'C' and spin2.element == 'H') or (spin2.element == 'C' and spin1.element == 'H'): 287 pair = True 288 289 # Same residue, so skip. 290 if pair and res_num1 != None and res_num1 != res_num2: 291 continue 292 elif pair and res_num1 == None and res_name1 != res_name2: 293 continue 294 295 # Store the IDs for the printout. 296 ids.append([id1, id2]) 297 spin_selections.append([spin1.select, spin2.select]) 298 299 # No matches, so fail! 300 if not len(ids): 301 # Find the problem. 302 count1 = 0 303 count2 = 0 304 for spin in spin_loop(spin_id1): 305 count1 += 1 306 for spin in spin_loop(spin_id2): 307 count2 += 1 308 309 # Report the problem. 310 if count1 == 0 and count2 == 0: 311 raise RelaxError("Neither spin IDs '%s' and '%s' match any spins." % (spin_id1, spin_id2)) 312 elif count1 == 0: 313 raise RelaxError("The spin ID '%s' matches no spins." % spin_id1) 314 elif count2 == 0: 315 raise RelaxError("The spin ID '%s' matches no spins." % spin_id2) 316 else: 317 raise RelaxError("Unknown error.") 318 319 # Define the interaction. 320 for i in range(len(ids)): 321 # Unpack. 322 id1, id2 = ids[i] 323 # Get the interatomic data object, if it exists. 324 interatom = return_interatom(id1, id2, pipe=pipe) 325 326 # Create the container if needed. 327 if interatom == None: 328 interatom = create_interatom(spin_id1=id1, spin_id2=id2, pipe=pipe) 329 330 # Check that this has not already been set up. 331 if interatom.dipole_pair: 332 raise RelaxError("The magnetic dipole-dipole interaction already exists between the spins '%s' and '%s'." % (id1, id2)) 333 334 # Set a flag indicating that a dipole-dipole interaction is present. 335 interatom.dipole_pair = True 336 337 # Set the selection. 338 if spin_selection: 339 interatom.select = False 340 if spin_selections[i][0] and spin_selections[i][1]: 341 interatom.select = True 342 343 # Printout. 344 if verbose: 345 # Conversion. 346 for i in range(len(ids)): 347 ids[i][0] = repr(ids[i][0]) 348 ids[i][1] = repr(ids[i][1]) 349 350 # The printout. 351 print("Interatomic interactions are now defined for the following spins:\n") 352 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2"], data=ids)

353 354

355 -def exists_data(pipe=None):

356 """Determine if any interatomic data exists. 357 358 @keyword pipe: The data pipe in which the interatomic data will be checked for. 359 @type pipe: str 360 @return: The answer to the question about the existence of data. 361 @rtype: bool 362 """ 363 364 # The current data pipe. 365 if pipe == None: 366 pipe = pipes.cdp_name() 367 368 # Test the data pipe. 369 check_pipe(pipe) 370 371 # Get the data pipe. 372 dp = pipes.get_pipe(pipe) 373 374 # The interatomic data structure is empty. 375 if dp.interatomic.is_empty(): 376 return False 377 378 # Otherwise. 379 return True

380 381

382 -def id_match(spin_id=None, interatom=None, pipe=None):

383 """Test if the spin ID matches one of the two spins of the given container. 384 385 @keyword spin_id: The spin ID string of the first atom. 386 @type spin_id: str 387 @keyword interatom: The interatomic data container. 388 @type interatom: InteratomContainer instance 389 @keyword pipe: The data pipe containing the interatomic data container. Defaults to the current data pipe. 390 @type pipe: str or None 391 @return: True if the spin ID matches one of the two spins, False otherwise. 392 @rtype: bool 393 """ 394 395 # Get the spin containers. 396 spin1 = return_spin(interatom.spin_id1, pipe=pipe) 397 spin2 = return_spin(interatom.spin_id2, pipe=pipe) 398 399 # No spins. 400 if spin1 == None or spin2 == None: 401 return False 402 403 # Check if the ID is in the private metadata list. 404 if spin_id in spin1._spin_ids or spin_id in spin2._spin_ids: 405 return True 406 407 # Nothing found. 408 return False

409 410

411 -def interatomic_loop(selection1=None, selection2=None, pipe=None, skip_desel=True):

412 """Generator function for looping over all the interatomic data containers. 413 414 @keyword selection1: The optional spin ID selection of the first atom. 415 @type selection1: str 416 @keyword selection2: The optional spin ID selection of the second atom. 417 @type selection2: str 418 @keyword pipe: The data pipe containing the spin. Defaults to the current data pipe. 419 @type pipe: str 420 @keyword skip_desel: A flag which if True will cause only selected interatomic data containers to be returned. 421 @type skip_desel: bool 422 """ 423 424 # The data pipe. 425 if pipe == None: 426 pipe = pipes.cdp_name() 427 428 # Get the data pipe. 429 dp = pipes.get_pipe(pipe) 430 431 # Parse the spin ID selection strings. 432 select_obj = None 433 select_obj1 = None 434 select_obj2 = None 435 if selection1 and selection2: 436 select_obj1 = Selection(selection1) 437 select_obj2 = Selection(selection2) 438 elif selection1: 439 select_obj = Selection(selection1) 440 elif selection2: 441 select_obj = Selection(selection2) 442 443 # Loop over the containers, yielding them. 444 for i in range(len(dp.interatomic)): 445 # Skip deselected containers. 446 if skip_desel and not dp.interatomic[i].select: 447 continue 448 449 # Aliases. 450 interatom = dp.interatomic[i] 451 mol_index1, res_index1, spin_index1 = dp.mol._spin_id_lookup[interatom.spin_id1] 452 mol_index2, res_index2, spin_index2 = dp.mol._spin_id_lookup[interatom.spin_id2] 453 mol1 = dp.mol[mol_index1] 454 res1 = dp.mol[mol_index1].res[res_index1] 455 spin1 = dp.mol[mol_index1].res[res_index1].spin[spin_index1] 456 mol2 = dp.mol[mol_index2] 457 res2 = dp.mol[mol_index2].res[res_index2] 458 spin2 = dp.mol[mol_index2].res[res_index2].spin[spin_index2] 459 460 # The different selection combinations. 461 if select_obj: 462 sel1 = select_obj.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 463 sel2 = select_obj.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 464 if select_obj1: 465 sel11 = select_obj1.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 466 sel12 = select_obj1.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 467 if select_obj2: 468 sel21 = select_obj2.contains_spin(spin_name=spin1.name, spin_num=spin1.num, res_name=res1.name, res_num=res1.num, mol=mol1.name) 469 sel22 = select_obj2.contains_spin(spin_name=spin2.name, spin_num=spin2.num, res_name=res2.name, res_num=res2.num, mol=mol2.name) 470 471 # Check that the selections are met. 472 if select_obj: 473 if not sel1 and not sel2: 474 continue 475 if select_obj1: 476 if not (sel11 or sel12): 477 continue 478 if select_obj2: 479 if not (sel21 or sel22): 480 continue 481 482 # Return the container. 483 yield interatom

484 485

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

487 """Set up the magnetic dipole-dipole interaction. 488 489 @keyword file: The name of the file to open. 490 @type file: str 491 @keyword dir: The directory containing the file (defaults to the current directory if None). 492 @type dir: str or None 493 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 494 @type unit: str 495 @keyword spin_id1_col: The column containing the spin ID strings of the first spin. 496 @type spin_id1_col: int 497 @keyword spin_id2_col: The column containing the spin ID strings of the second spin. 498 @type spin_id2_col: int 499 @keyword data_col: The column containing the averaged distances in meters. 500 @type data_col: int or None 501 @keyword sep: The column separator which, if None, defaults to whitespace. 502 @type sep: str or None 503 """ 504 505 # Check the units. 506 if unit not in ['meter', 'Angstrom']: 507 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 508 509 # Test if the current data pipe exists. 510 check_pipe() 511 512 # Test if sequence data exists. 513 if not exists_mol_res_spin_data(): 514 raise RelaxNoSequenceError 515 516 # Extract the data from the file, and clean it up. 517 file_data = extract_data(file, dir, sep=sep) 518 file_data = strip(file_data, comments=True) 519 520 # Loop over the RDC data. 521 data = [] 522 for line in file_data: 523 # Invalid columns. 524 if spin_id1_col > len(line): 525 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line)) 526 continue 527 if spin_id2_col > len(line): 528 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line)) 529 continue 530 if data_col and data_col > len(line): 531 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line)) 532 continue 533 534 # Unpack. 535 spin_id1 = line[spin_id1_col-1] 536 spin_id2 = line[spin_id2_col-1] 537 ave_dist = None 538 if data_col: 539 ave_dist = line[data_col-1] 540 541 # Convert and check the value. 542 if ave_dist != None: 543 try: 544 ave_dist = float(ave_dist) 545 except ValueError: 546 warn(RelaxWarning("The averaged distance of '%s' from the line %s is invalid." % (ave_dist, line))) 547 continue 548 549 # Unit conversion. 550 if unit == 'Angstrom': 551 ave_dist = ave_dist * 1e-10 552 553 # Get the interatomic data container. 554 interatom = return_interatom(spin_id1, spin_id2) 555 556 # No container found, so create it. 557 if interatom == None: 558 interatom = create_interatom(spin_id1=spin_id1, spin_id2=spin_id2, verbose=True) 559 560 # Store the averaged distance. 561 interatom.r = ave_dist 562 563 # Store the data for the printout. 564 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 565 566 # No data, so fail! 567 if not len(data): 568 raise RelaxError("No data could be extracted from the file.") 569 570 # Print out. 571 print("The following averaged distances have been read:\n") 572 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)

573 574

575 -def return_interatom(spin_id1=None, spin_id2=None, pipe=None):

576 """Return the list of interatomic data containers for the two spins. 577 578 @keyword spin_id1: The spin ID string of the first atom. 579 @type spin_id1: str 580 @keyword spin_id2: The spin ID string of the second atom. 581 @type spin_id2: str 582 @keyword pipe: The data pipe holding the container. Defaults to the current data pipe. 583 @type pipe: str or None 584 @return: The matching interatomic data container, if it exists. 585 @rtype: data.interatomic.InteratomContainer instance or None 586 """ 587 588 # Checks. 589 if spin_id1 == None: 590 raise RelaxError("The first spin ID must be supplied.") 591 if spin_id2 == None: 592 raise RelaxError("The second spin ID must be supplied.") 593 594 # The data pipe. 595 if pipe == None: 596 pipe = pipes.cdp_name() 597 598 # Get the data pipe. 599 dp = pipes.get_pipe(pipe) 600 601 # Return the matching container. 602 for i in range(len(dp.interatomic)): 603 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): 604 return dp.interatomic[i] 605 606 # No matchs. 607 return None

608 609

610 -def return_interatom_list(spin_id=None, pipe=None):

611 """Return the list of interatomic data containers for the given spin. 612 613 @keyword spin_id: The spin ID string. 614 @type spin_id: str 615 @keyword pipe: The data pipe holding the container. This defaults to the current data pipe. 616 @type pipe: str or None 617 @return: The list of matching interatomic data containers, if any exist. 618 @rtype: list of data.interatomic.InteratomContainer instances 619 """ 620 621 # Check. 622 if spin_id == None: 623 raise RelaxError("The spin ID must be supplied.") 624 625 # The data pipe. 626 if pipe == None: 627 pipe = pipes.cdp_name() 628 629 # Get the data pipe. 630 dp = pipes.get_pipe(pipe) 631 632 # Initialise. 633 interatoms = [] 634 635 # Find and append all containers. 636 for i in range(len(dp.interatomic)): 637 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): 638 interatoms.append(dp.interatomic[i]) 639 640 # Return the list of containers. 641 return interatoms

642 643

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

645 """Set up the magnetic dipole-dipole interaction. 646 647 @keyword spin_id1: The spin identifier string of the first spin of the pair. 648 @type spin_id1: str 649 @keyword spin_id2: The spin identifier string of the second spin of the pair. 650 @type spin_id2: str 651 @keyword ave_dist: The r^-3 averaged interatomic distance. 652 @type ave_dist: float 653 @keyword unit: The measurement unit. This can be either 'meter' or 'Angstrom'. 654 @type unit: str 655 """ 656 657 # Check the units. 658 if unit not in ['meter', 'Angstrom']: 659 raise RelaxError("The measurement unit of '%s' must be one of 'meter' or 'Angstrom'." % unit) 660 661 # Unit conversion. 662 if unit == 'Angstrom': 663 ave_dist = ave_dist * 1e-10 664 665 # Generate the selection objects. 666 sel_obj1 = Selection(spin_id1) 667 sel_obj2 = Selection(spin_id2) 668 669 # Loop over the interatomic containers. 670 data = [] 671 for interatom in interatomic_loop(): 672 # Get the spin info. 673 mol_name1, res_num1, res_name1, spin1 = return_spin(interatom.spin_id1, full_info=True) 674 mol_name2, res_num2, res_name2, spin2 = return_spin(interatom.spin_id2, full_info=True) 675 676 # No match, either way. 677 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)): 678 continue 679 680 # Store the averaged distance. 681 interatom.r = ave_dist 682 683 # Store the data for the printout. 684 data.append([repr(interatom.spin_id1), repr(interatom.spin_id2), repr(ave_dist)]) 685 686 # No data, so fail! 687 if not len(data): 688 raise RelaxError("No data could be set.") 689 690 # Print out. 691 print("The following averaged distances have been set:\n") 692 write_data(out=sys.stdout, headings=["Spin_ID_1", "Spin_ID_2", "Ave_distance(meters)"], data=data)

693 694

695 -def unit_vectors(ave=True):

696 """Extract the bond vectors from the loaded structures and store them in the spin container. 697 698 @keyword ave: A flag which if True will cause the average of all vectors to be calculated. 699 @type ave: bool 700 """ 701 702 # Test if the current data pipe exists. 703 check_pipe() 704 705 # Test if interatomic data exists. 706 if not exists_data(): 707 raise RelaxNoInteratomError 708 709 # Print out. 710 if ave: 711 print("Averaging all vectors.") 712 else: 713 print("No averaging of the vectors.") 714 715 # Loop over the interatomic data containers. 716 no_vectors = True 717 pos_info = False 718 for interatom in interatomic_loop(skip_desel=False): 719 # Get the spin info. 720 spin1 = return_spin(interatom.spin_id1) 721 spin2 = return_spin(interatom.spin_id2) 722 723 # No positional information. 724 if not hasattr(spin1, 'pos'): 725 continue 726 if not hasattr(spin2, 'pos'): 727 continue 728 729 # Positional information flag. 730 pos_info = True 731 732 # Both single positions. 733 if is_float(spin1.pos[0], raise_error=False) and is_float(spin2.pos[0], raise_error=False): 734 # The vector. 735 vector_list = [spin2.pos - spin1.pos] 736 737 # A single and multiple position pair. 738 elif is_float(spin1.pos[0], raise_error=False) or is_float(spin2.pos[0], raise_error=False): 739 # The first spin has multiple positions. 740 if is_float(spin2.pos[0], raise_error=False): 741 vector_list = [] 742 for i in range(len(spin1.pos)): 743 vector_list.append(spin2.pos - spin1.pos[i]) 744 745 # The second spin has multiple positions. 746 else: 747 vector_list = [] 748 for i in range(len(spin2.pos)): 749 vector_list.append(spin2.pos[i] - spin1.pos) 750 751 # Both spins have multiple positions. 752 else: 753 # Non-matching number of positions. 754 if len(spin1.pos) != len(spin2.pos): 755 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))) 756 757 # Calculate all vectors. 758 vector_list = [] 759 for i in range(len(spin1.pos)): 760 # No structural information. 761 if spin1.pos[i] == None or spin2.pos[i] == None: 762 warn(RelaxWarning("No structural information for state %i can be found between spins '%s' and '%s'." % (i, interatom.spin_id1, interatom.spin_id2))) 763 vector_list.append(None) 764 765 # All data is present. 766 else: 767 vector_list.append(spin2.pos[i] - spin1.pos[i]) 768 769 # Unit vectors. 770 for i in range(len(vector_list)): 771 # No vector. 772 if vector_list[i] == None: 773 continue 774 775 # Normalisation factor. 776 norm_factor = norm(vector_list[i]) 777 778 # Test for zero length. 779 if norm_factor == 0.0: 780 warn(RelaxZeroVectorWarning(spin_id1=interatom.spin_id1, spin_id2=interatom.spin_id2)) 781 782 # Calculate the normalised vector. 783 else: 784 vector_list[i] = vector_list[i] / norm_factor 785 786 # Average. 787 if ave: 788 ave_vector = zeros(3, float64) 789 count = 0 790 for i in range(len(vector_list)): 791 if vector_list[i] != None: 792 ave_vector = ave_vector + vector_list[i] 793 count += 1 794 vector_list = [ave_vector / count] 795 796 # Convert to a single vector if needed. 797 if len(vector_list) == 1: 798 vector_list = vector_list[0] 799 800 # Store the unit vector(s). 801 setattr(interatom, 'vector', vector_list) 802 803 # We have a vector! 804 no_vectors = False 805 806 # Printout. 807 num = 1 808 if not is_float(vector_list[0], raise_error=False): 809 num = len(vector_list) 810 plural = 's' 811 if num == 1: 812 plural = '' 813 if spin1.name: 814 spin1_str = spin1.name 815 else: 816 spin1_str = spin1.num 817 if spin2.name: 818 spin2_str = spin2.name 819 else: 820 spin2_str = spin2.num 821 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)) 822 823 # Catch the problem of no positional information being present. 824 if not pos_info: 825 raise RelaxError("Positional information could not be found for any spins.") 826 827 # Right, catch the problem of missing vectors to prevent massive user confusion! 828 if no_vectors: 829 raise RelaxError("No vectors could be extracted.")

830

Source Code for Module pipe_control.interatomic