pipe_control.structure.main

60 """Add a new atom to the structural data object. 61 62 @keyword mol_name: The name of the molecule. 63 @type mol_name: str 64 @keyword atom_name: The atom name, e.g. 'H1'. 65 @type atom_name: str or None 66 @keyword res_name: The residue name. 67 @type res_name: str or None 68 @keyword res_num: The residue number. 69 @type res_num: int or None 70 @keyword pos: The position vector of coordinates. If a rank-2 array is supplied, the length of the first dimension must match the number of models. 71 @type pos: rank-1 or rank-2 array or list of float 72 @keyword element: The element symbol. 73 @type element: str or None 74 @keyword atom_num: The atom number. 75 @type atom_num: int or None 76 @keyword chain_id: The chain identifier. 77 @type chain_id: str or None 78 @keyword segment_id: The segment identifier. 79 @type segment_id: str or None 80 @keyword pdb_record: The optional PDB record name, e.g. 'ATOM' or 'HETATM'. 81 @type pdb_record: str or None 82 """ 83 84 # Test if the current data pipe exists. 85 check_pipe() 86 87 # Place the structural object into the relax data store if needed. 88 if not hasattr(cdp, 'structure'): 89 cdp.structure = Internal() 90 91 # Add the atoms. 92 cdp.structure.add_atom(mol_name=mol_name, atom_name=atom_name, res_name=res_name, res_num=res_num, pos=pos, element=element, atom_num=atom_num, chain_id=chain_id, segment_id=segment_id, pdb_record=pdb_record, sort=True)

115 """Assemble the common atomic coordinates taking sequence alignments into account. 116 117 @keyword pipes: The data pipes to assemble the coordinates from. 118 @type pipes: None or list of str 119 @keyword models: The list of models for each data pipe. The number of elements must match the pipes argument. If set to None, then all models will be used. 120 @type models: None or list of lists of int 121 @keyword molecules: The list of molecules for each data pipe. The number of elements must match the pipes argument. 122 @type molecules: None or list of lists of str 123 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 124 @type atom_id: str or None 125 @return: The array of atomic coordinates (first dimension is the model and/or molecule, the second are the atoms, and the third are the coordinates); a list of unique IDs for each structural object, model, and molecule; the common list of molecule names; the common list of residue names; the common list of residue numbers; the common list of atom names; the common list of element names. 126 @rtype: numpy rank-3 float64 array, list of str, list of str, list of str, list of int, list of str, list of str 127 """ 128 129 # Assemble the structural objects. 130 objects, object_names, pipes = assemble_structural_objects(pipes=pipes, models=models, molecules=molecules) 131 132 # Assemble the atomic coordinates of all molecules. 133 ids, object_id_list, model_list, molecule_list, atom_pos, mol_names, res_names, res_nums, atom_names, elements, one_letter_codes, num_mols = assemble_atomic_coordinates(objects=objects, object_names=object_names, molecules=molecules, models=models, atom_id=atom_id) 134 135 # No data. 136 if mol_names == []: 137 if atom_id != None: 138 raise RelaxError("No structural data matching the atom ID string '%s' can be found." % atom_id) 139 else: 140 raise RelaxError("No structural data can be found.") 141 142 # Are all molecules the same? 143 same_mol = True 144 for mol in molecule_list: 145 if mol != molecule_list[0]: 146 same_mol = False 147 148 # Handle sequence alignments - retrieve the alignment. 149 align = None 150 if hasattr(ds, 'sequence_alignments'): 151 align = ds.sequence_alignments.find_alignment(object_ids=object_id_list, models=model_list, molecules=molecule_list, sequences=one_letter_codes) 152 if align != None: 153 # Printout. 154 print("\nSequence alignment found - common atoms will be determined based on this MSA:") 155 for i in range(len(align.object_ids)): 156 print(align.strings[i]) 157 158 # Alias the required data structures. 159 strings = align.strings 160 gaps = align.gaps 161 162 # Handle sequence alignments - no alignment required. 163 elif len(objects) == 1 and same_mol: 164 # Printout. 165 print("\nSequence alignment disabled as only models with identical molecule, residue and atomic sequences are being superimposed.") 166 167 # Set the one letter codes to be the alignment strings. 168 strings = one_letter_codes 169 170 # Create an empty gap data structure. 171 gaps = [] 172 for mol_index in range(num_mols): 173 gaps.append([]) 174 for i in range(len(one_letter_codes[mol_index])): 175 gaps[mol_index].append(0) 176 177 # Handle sequence alignments - fall back alignment based on residue numbering. 178 else: 179 # Printout. 180 print("\nSequence alignment cannot be found - falling back to a residue number based alignment.") 181 182 # Convert the residue number data structure. 183 res_num_list = [] 184 for mol_index in range(num_mols): 185 res_num_list.append([]) 186 for i in range(len(one_letter_codes[mol_index])): 187 key = list(res_nums[mol_index][i].keys())[0] 188 res_num_list[mol_index].append(res_nums[mol_index][i][key]) 189 190 # Sequence alignment. 191 strings, gaps = msa_residue_numbers(one_letter_codes, residue_numbers=res_num_list) 192 193 # Create the residue skipping data structure. 194 skip = msa_residue_skipping(strings=strings, gaps=gaps) 195 196 # Assemble and return the atomic coordinates and common atom information. 197 coord, mol_name_common, res_name_common, res_num_common, atom_name_common, element_common = assemble_coord_array(atom_pos=atom_pos, mol_names=mol_names, res_names=res_names, res_nums=res_nums, atom_names=atom_names, elements=elements, sequences=one_letter_codes, skip=skip) 198 return coord, ids, mol_name_common, res_name_common, res_num_common, atom_name_common, element_common

202 """Assemble the atomic coordinates. 203 204 @keyword pipes: The data pipes to assemble the coordinates from. 205 @type pipes: None or list of str 206 @keyword models: The list of models for each data pipe. The number of elements must match the pipes argument. If set to None, then all models will be used. 207 @type models: None or list of lists of int 208 @keyword molecules: The list of molecules for each data pipe. The number of elements must match the pipes argument. 209 @type molecules: None or list of lists of str 210 @return: The structural objects, structural object names, and data pipes list. 211 @rtype: list of lib.structure.internal.object.Internal instances, list of str, list of str 212 """ 213 214 # The data pipes to use. 215 if pipes == None: 216 pipes = [cdp_name()] 217 num_pipes = len(pipes) 218 219 # Checks. 220 for pipe in pipes: 221 check_pipe(pipe) 222 223 # Check the models and molecules arguments. 224 if models != None: 225 if len(models) != num_pipes: 226 raise RelaxError("The %i elements of the models argument does not match the %i data pipes." % (len(models), num_pipes)) 227 if molecules != None: 228 if len(molecules) != num_pipes: 229 raise RelaxError("The %i elements of the molecules argument does not match the %i data pipes." % (len(molecules), num_pipes)) 230 231 # Assemble the structural objects. 232 objects = [] 233 object_names = [] 234 for pipe_index in range(len(pipes)): 235 dp = get_pipe(pipes[pipe_index]) 236 objects.append(dp.structure) 237 object_names.append(pipes[pipe_index]) 238 239 # Return the structural objects, object names, and the new pipes list. 240 return objects, object_names, pipes

244 """Write out a correlation matrix of pairwise interatomic distance fluctuations between different structures. 245 246 @keyword pipes: The data pipes to generate the interatomic distance fluctuation correlation matrix for. 247 @type pipes: None or list of str 248 @keyword models: The list of models to generate the interatomic distance fluctuation correlation matrix for. The number of elements must match the pipes argument. If set to None, then all models will be used. 249 @type models: None or list of lists of int 250 @keyword molecules: The list of molecules to generate the interatomic distance fluctuation correlation matrix for. The number of elements must match the pipes argument. 251 @type molecules: None or list of lists of str 252 @keyword atom_id: The atom identification string of the coordinates of interest. This matches the spin ID string format. 253 @type atom_id: str or None 254 @keyword measure: The type of fluctuation to measure. This can be either 'distance' or 'angle'. 255 @type measure: str 256 @keyword file: The name of the file to write. 257 @type file: str 258 @keyword format: The output format. This can be set to "text" for text file output, or "gnuplot" for creating a gnuplot script. 259 @type format: str 260 @keyword dir: The directory where the file will be placed. If set to None, then the file will be placed in the current directory. 261 @type dir: str or None 262 @keyword force: The force flag which if True will cause the file to be overwritten. 263 @type force: bool 264 """ 265 266 # Checks. 267 check_pipe() 268 check_structure() 269 allowed_measures = ['distance', 'angle', 'parallax shift'] 270 if measure not in allowed_measures: 271 raise RelaxError("The measure '%s' must be one of %s." % (measure, allowed_measures)) 272 273 # Assemble the structural coordinates. 274 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 275 276 # The number of dimensions. 277 n = len(atom_names) 278 m = len(coord) 279 280 # Check that more than one structure is present. 281 if not m > 1: 282 raise RelaxError("Two or more structures are required.") 283 284 # The labels as spin ID strings. 285 labels = [] 286 for i in range(n): 287 labels.append(generate_spin_id_unique(mol_name=mol_names[i], res_num=res_nums[i], res_name=res_names[i], spin_name=atom_names[i])) 288 289 # Initialise the SD matrix and other structures. 290 matrix = zeros((n, n), float64) 291 dist = zeros(m, float64) 292 vectors = zeros((m, 3), float64) 293 parallax_vectors = zeros((m, 3), float64) 294 angles = zeros(m, float64) 295 296 # Generate the pairwise distance SD matrix. 297 if measure == 'distance': 298 for i in range(n): 299 for j in range(n): 300 # Only calculate the upper triangle to avoid duplicate calculations. 301 if j > i: 302 continue 303 304 # The interatomic distances between each structure. 305 for k in range(m): 306 dist[k] = norm(coord[k, i] - coord[k, j]) 307 308 # Calculate and store the corrected sample standard deviation. 309 matrix[i, j] = matrix[j, i] = std(dist, ddof=1) 310 311 # Generate the pairwise angle SD matrix. 312 elif measure == 'angle': 313 # Loop over the atom pairs. 314 for i in range(n): 315 for j in range(n): 316 # Only calculate the upper triangle to avoid duplicate calculations. 317 if j > i: 318 continue 319 320 # The interatomic vectors between each structure. 321 for k in range(m): 322 vectors[k] = coord[k, i] - coord[k, j] 323 324 # The average vector. 325 ave_vect = average(vectors, axis=0) 326 327 # The intervector angles. 328 for k in range(m): 329 angles[k] = vector_angle_atan2(ave_vect, vectors[k]) 330 331 # Calculate and store the corrected sample standard deviation. 332 matrix[i, j] = matrix[j, i] = std(angles, ddof=1) 333 334 # Generate the pairwise parallax shift SD matrix. 335 elif measure == 'parallax shift': 336 # Loop over the atom pairs. 337 for i in range(n): 338 for j in range(n): 339 # Only calculate the upper triangle to avoid duplicate calculations. 340 if j > i: 341 continue 342 343 # The interatomic vectors between each structure. 344 for k in range(m): 345 vectors[k] = coord[k, i] - coord[k, j] 346 347 # The average vector. 348 ave_vect = average(vectors, axis=0) 349 350 # Catch the zero vector. 351 length = norm(ave_vect) 352 if length == 0.0: 353 matrix[i, j] = matrix[j, i] = 0.0 354 continue 355 356 # The unit average vector. 357 unit = ave_vect / length 358 359 # The parallax shift. 360 for k in range(m): 361 # The projection onto the average vector. 362 proj = dot(vectors[k], unit) * unit 363 364 # The distance shift. 365 dist[k] = norm(vectors[k] - proj) 366 367 # Calculate and store the corrected sample standard deviation. 368 matrix[i, j] = matrix[j, i] = std(dist, ddof=1) 369 370 # Call the plotting API. 371 correlation_matrix(format=format, matrix=matrix, labels=labels, file=file, dir=dir, force=force)

375 """Calculate a mean structure. 376 377 @keyword pipes: The data pipes containing structures to average. 378 @type pipes: None or list of str 379 @keyword models: The list of models for each data pipe. The number of elements must match the pipes argument. If set to None, then all models will be used. 380 @type models: None or list of lists of int 381 @keyword molecules: The list of molecules for each data pipe. The number of elements must match the pipes argument. 382 @type molecules: None or list of lists of str 383 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 384 @type atom_id: str or None 385 @keyword set_mol_name: The molecule name for the averaged molecule. 386 @type set_mol_name: None or str 387 @keyword set_model_num: The model number for the averaged molecule. 388 @type set_model_num: None or int 389 """ 390 391 # Checks. 392 check_pipe() 393 394 # Assemble the structural coordinates. 395 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 396 397 # Calculate the mean structure. 398 struct = mean(coord, axis=0) 399 400 # Place the structural object into the relax data store if needed. 401 if not hasattr(cdp, 'structure'): 402 cdp.structure = Internal() 403 404 # Store the data. 405 cdp.structure.add_coordinates(coord=struct, mol_names=mol_names, res_names=res_names, res_nums=res_nums, atom_names=atom_names, elements=elements, set_mol_name=set_mol_name, set_model_num=set_model_num)

448 """Create the PDB representation of the diffusion tensor. 449 450 @keyword scale: The scaling factor for the diffusion tensor. 451 @type scale: float 452 @keyword file: The name of the PDB file to create. 453 @type file: str 454 @keyword dir: The name of the directory to place the PDB file into. 455 @type dir: str 456 @keyword force: Flag which if set to True will overwrite any pre-existing file. 457 @type force: bool 458 """ 459 460 # Test if the current data pipe exists. 461 check_pipe() 462 463 # Calculate the centre of mass. 464 if hasattr(cdp, 'structure') and not cdp.structure.empty(): 465 com = pipe_centre_of_mass() 466 else: 467 com = zeros(3, float64) 468 469 # Create the structural object. 470 structure = Internal() 471 472 # Create an array of data pipes to loop over (hybrid support). 473 if cdp.pipe_type == 'hybrid': 474 pipe_list = cdp.hybrid_pipes 475 else: 476 pipe_list = [pipes.cdp_name()] 477 478 # The molecule names. 479 if cdp.pipe_type == 'hybrid': 480 mol_names = [] 481 for pipe in pipe_list: 482 mol_names.append('diff_tensor_' % pipe) 483 else: 484 mol_names = ['diff_tensor'] 485 486 # Loop over the pipes. 487 for pipe_index in range(len(pipe_list)): 488 # Get the pipe container. 489 pipe = pipes.get_pipe(pipe_list[pipe_index]) 490 491 # Test if the diffusion tensor data is loaded. 492 if not hasattr(pipe, 'diff_tensor'): 493 raise RelaxNoTensorError('diffusion') 494 495 # Add a new structure. 496 structure.add_molecule(name=mol_names[pipe_index]) 497 498 # Alias the single molecule from the single model. 499 mol = structure.get_molecule(mol_names[pipe_index]) 500 501 # The diffusion tensor type. 502 diff_type = pipe.diff_tensor.type 503 if diff_type == 'spheroid': 504 diff_type = pipe.diff_tensor.spheroid_type 505 506 # Simulation info. 507 sim_num = None 508 if hasattr(pipe.diff_tensor, 'tm_sim'): 509 # The number. 510 sim_num = len(pipe.diff_tensor.tm_sim) 511 512 # Tensor axes. 513 axes = [] 514 sim_axes = [] 515 if diff_type in ['oblate', 'prolate']: 516 axes.append(pipe.diff_tensor.Dpar * pipe.diff_tensor.Dpar_unit) 517 sim_axes.append([]) 518 if sim_num != None: 519 for i in range(sim_num): 520 sim_axes[0].append(pipe.diff_tensor.Dpar_sim[i] * pipe.diff_tensor.Dpar_unit_sim[i]) 521 522 if diff_type == 'ellipsoid': 523 axes.append(pipe.diff_tensor.Dx * pipe.diff_tensor.Dx_unit) 524 axes.append(pipe.diff_tensor.Dy * pipe.diff_tensor.Dy_unit) 525 axes.append(pipe.diff_tensor.Dz * pipe.diff_tensor.Dz_unit) 526 sim_axes.append([]) 527 sim_axes.append([]) 528 sim_axes.append([]) 529 if sim_num != None: 530 for i in range(sim_num): 531 sim_axes[0].append(pipe.diff_tensor.Dx_sim[i] * pipe.diff_tensor.Dx_unit_sim[i]) 532 sim_axes[1].append(pipe.diff_tensor.Dy_sim[i] * pipe.diff_tensor.Dy_unit_sim[i]) 533 sim_axes[2].append(pipe.diff_tensor.Dz_sim[i] * pipe.diff_tensor.Dz_unit_sim[i]) 534 535 # Create the object. 536 diffusion_tensor(mol=mol, tensor=pipe.diff_tensor.tensor, tensor_diag=pipe.diff_tensor.tensor_diag, diff_type=diff_type, rotation=pipe.diff_tensor.rotation, axes=axes, sim_axes=sim_axes, com=com, scale=scale) 537 538 539 # Create the PDB file. 540 ###################### 541 542 # Print out. 543 print("\nGenerating the PDB file.") 544 545 # Create the PDB file. 546 tensor_pdb_file = open_write_file(file, dir, force=force) 547 structure.write_pdb(tensor_pdb_file) 548 tensor_pdb_file.close() 549 550 # Add the file to the results file list. 551 if not hasattr(cdp, 'result_files'): 552 cdp.result_files = [] 553 if dir == None: 554 dir = getcwd() 555 cdp.result_files.append(['diff_tensor_pdb', 'Diffusion tensor PDB', get_file_path(file, dir)]) 556 status.observers.result_file.notify()

560 """Delete structural data. 561 562 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. If not given, then all structural data will be deleted. 563 @type atom_id: str or None 564 @keyword model: Individual structural models from a loaded ensemble can be deleted by specifying the model number. 565 @type model: None or int 566 @keyword verbosity: The amount of information to print to screen. Zero corresponds to minimal output while higher values increase the amount of output. The default value is 1. 567 @type verbosity: int 568 @keyword spin_info: A flag which if True will cause all structural information in the spin containers and interatomic data containers to be deleted as well. If False, then only the 3D structural data will be deleted. 569 @type spin_info: bool 570 """ 571 572 # Test if the current data pipe exists. 573 check_pipe() 574 575 # Run the object method. 576 if hasattr(cdp, 'structure'): 577 if verbosity: 578 print("Deleting structural data from the current pipe.") 579 selection = None 580 if atom_id != None: 581 selection = cdp.structure.selection(atom_id=atom_id) 582 cdp.structure.delete(model=model, selection=selection, verbosity=verbosity) 583 elif verbosity: 584 print("No structures are present.") 585 586 # Skip the rest. 587 if not spin_info: 588 return 589 590 # Then remove any spin specific structural info. 591 if verbosity: 592 print("Deleting all spin specific structural info.") 593 for spin in spin_loop(selection=atom_id): 594 # Delete positional information. 595 if hasattr(spin, 'pos'): 596 del spin.pos 597 598 # Then remove any interatomic vector structural info. 599 if verbosity: 600 print("Deleting all interatomic vectors.") 601 for interatom in interatomic_loop(selection1=atom_id): 602 # Delete bond vectors. 603 if hasattr(interatom, 'vector'): 604 del interatom.vector

608 """Calculate the rotational and translational displacement between structures or models. 609 610 All results will be placed into the current data pipe cdp.structure.displacements data structure. 611 612 613 @keyword pipes: The data pipes to determine the displacements for. 614 @type pipes: None or list of str 615 @keyword models: The list of models to determine the displacements for. The number of elements must match the pipes argument. If set to None, then all models will be used. 616 @type models: None or list of lists of int 617 @keyword molecules: The list of molecules to determine the displacements for. The number of elements must match the pipes argument. 618 @type molecules: None or list of lists of str 619 @keyword atom_id: The atom identification string of the coordinates of interest. This matches the spin ID string format. 620 @type atom_id: str or None 621 @keyword centroid: An alternative position of the centroid, used for studying pivoted systems. 622 @type centroid: list of float or numpy rank-1, 3D array 623 """ 624 625 # Test if the current data pipe exists. 626 check_pipe() 627 628 # Assemble the structural coordinates. 629 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 630 631 # Initialise the data structure. 632 if not hasattr(cdp.structure, 'displacments'): 633 cdp.structure.displacements = Displacements() 634 635 # Double loop over all structures, sending the data to the base container for the calculations. 636 for i in range(len(ids)): 637 for j in range(len(ids)): 638 cdp.structure.displacements._calculate(id_from=ids[i], id_to=ids[j], coord_from=coord[i], coord_to=coord[j], centroid=centroid)

642 """Find the pivoted motion of a set of structural models or structures. 643 644 The pivot will be placed into the current data pipe cdp.structure.pivot data structure. 645 646 647 @keyword pipes: The data pipes to use in the motional pivot algorithm. 648 @type pipes: None or list of str 649 @keyword models: The list of models to use. The number of elements must match the pipes argument. If set to None, then all models will be used. 650 @type models: None or list of lists of int 651 @keyword molecules: The list of molecules to find the pivoted motion for. The number of elements must match the pipes argument. 652 @type molecules: None or list of lists of str 653 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 654 @type atom_id: str or None 655 @keyword init_pos: The starting pivot position for the pivot point optimisation. 656 @type init_pos: list of float or numpy rank-1, 3D array 657 @keyword func_tol: The function tolerance which, when reached, terminates optimisation. Setting this to None turns of the check. 658 @type func_tol: None or float 659 @keyword box_limit: The simplex optimisation used in this function is constrained withing a box of +/- x Angstrom containing the pivot point using the logarithmic barrier function. This argument is the value of x. 660 @type box_limit: int 661 """ 662 663 # Test if the current data pipe exists. 664 check_pipe() 665 666 # Initialised the starting position if needed. 667 if init_pos == None: 668 init_pos = zeros(3, float64) 669 init_pos = array(init_pos) 670 671 # Assemble the structural coordinates. 672 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 673 674 # Linear constraints for the pivot position (between -1000 and 1000 Angstrom). 675 A = zeros((6, 3), float64) 676 b = zeros(6, float64) 677 for i in range(3): 678 A[2*i, i] = 1 679 A[2*i+1, i] = -1 680 b[2*i] = -box_limit 681 b[2*i+1] = -box_limit 682 683 # The target function. 684 finder = Pivot_finder(list(range(len(coord))), coord) 685 results = generic_minimise(func=finder.func, x0=init_pos, min_algor='Log barrier', min_options=('simplex',), A=A, b=b, func_tol=func_tol, print_flag=1) 686 687 # No result. 688 if results == None: 689 return 690 691 # Store the data. 692 cdp.structure.pivot = results 693 694 # Print out. 695 print("Motional pivot found at: %s" % results)

699 """Load the spins from the structural object into the relax data store. 700 701 @keyword spin_id: The molecule, residue, and spin identifier string. 702 @type spin_id: str 703 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 704 @type str_id: int or str 705 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 706 @type ave_pos: bool 707 """ 708 709 # Checks. 710 check_pipe() 711 check_structure() 712 713 # The selection object. 714 selection = cdp.structure.selection(atom_id=spin_id) 715 716 # Loop over all atoms of the spin_id selection. 717 data = [] 718 for mol_name, res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(selection=selection, str_id=str_id, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True, ave=ave_pos): 719 # Remove the '+' regular expression character from the mol, res, and spin names! 720 if mol_name and search('\+', mol_name): 721 mol_name = mol_name.replace('+', '') 722 if res_name and search('\+', res_name): 723 res_name = res_name.replace('+', '') 724 if atom_name and search('\+', atom_name): 725 atom_name = atom_name.replace('+', '') 726 727 # The spin identification string. 728 id = generate_spin_id_unique(res_num=res_num, res_name=None, spin_num=atom_num, spin_name=atom_name) 729 730 # Get the spin container. 731 spin_cont = return_spin(id) 732 733 # Skip the spin if it doesn't exist. 734 if spin_cont == None: 735 continue 736 737 # Add the position vector to the spin container. 738 spin_cont.pos = pos 739 740 # Store the data for a printout at the end. 741 data.append([id, repr(pos)]) 742 743 # No positions found. 744 if not len(data): 745 raise RelaxError("No positional information matching the spin ID '%s' could be found." % spin_id) 746 747 # Update pseudo-atoms. 748 for spin in spin_loop(): 749 if hasattr(spin, 'members'): 750 # Get the spin positions. 751 positions = [] 752 for atom in spin.members: 753 # Get the spin container. 754 subspin = return_spin(atom) 755 756 # Test that the spin exists. 757 if subspin == None: 758 raise RelaxNoSpinError(atom) 759 760 # Test the position. 761 if not hasattr(subspin, 'pos') or subspin.pos == None or not len(subspin.pos): 762 raise RelaxError("Positional information is not available for the atom '%s'." % atom) 763 764 # Alias the position. 765 pos = subspin.pos 766 767 # Convert to a list of lists if not already. 768 multi_model = True 769 if type(pos[0]) in [float, float64]: 770 multi_model = False 771 pos = [pos] 772 773 # Store the position. 774 positions.append([]) 775 for i in range(len(pos)): 776 positions[-1].append(pos[i].tolist()) 777 778 # The averaging. 779 if spin.averaging == 'linear': 780 # Average pos. 781 ave = linear_ave(positions) 782 783 # Convert to the correct structure. 784 if multi_model: 785 spin.pos = ave 786 else: 787 spin.pos = ave[0] 788 789 # Print out. 790 write_data(out=sys.stdout, headings=["Spin_ID", "Position"], data=data)

794 """Load the spins from the structural object into the relax data store. 795 796 @keyword spin_id: The molecule, residue, and spin identifier string. 797 @type spin_id: str 798 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 799 @type str_id: int or str 800 @keyword from_mols: The list of similar, but not necessarily identical molecules to load spin information from. 801 @type from_mols: list of str or None 802 @keyword mol_name_target: The name of target molecule container, overriding the name of the loaded structures 803 @type mol_name_target: str or None 804 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 805 @type ave_pos: bool 806 """ 807 808 # The multi-molecule case. 809 if from_mols != None: 810 load_spins_multi_mol(spin_id=spin_id, str_id=str_id, from_mols=from_mols, mol_name_target=mol_name_target, ave_pos=ave_pos) 811 return 812 813 # Checks. 814 check_pipe() 815 check_structure() 816 817 # Print out. 818 print("Adding the following spins to the relax data store.\n") 819 820 # Init the data for printing out. 821 mol_names = [] 822 res_nums = [] 823 res_names = [] 824 spin_nums = [] 825 spin_names = [] 826 827 # Loop over all atoms of the spin_id selection. 828 selection = cdp.structure.selection(atom_id=spin_id) 829 for mol_name, res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(selection=selection, str_id=str_id, mol_name_flag=True, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True, ave=ave_pos): 830 # Override the molecule name. 831 if mol_name_target: 832 mol_name = mol_name_target 833 834 # Remove the '+' regular expression character from the mol, res, and spin names! 835 if mol_name and search('\+', mol_name): 836 mol_name = mol_name.replace('+', '') 837 if res_name and search('\+', res_name): 838 res_name = res_name.replace('+', '') 839 if atom_name and search('\+', atom_name): 840 atom_name = atom_name.replace('+', '') 841 842 # Generate a spin ID for the current atom. 843 id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 844 845 # Create the spin. 846 try: 847 spin_cont = create_spin(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=atom_num, spin_name=atom_name) 848 849 # Otherwise, get the spin container. 850 except RelaxError: 851 spin_cont = return_spin(id) 852 853 # Append all the spin ID info for printing later. 854 if mol_name_target: 855 mol_names.append(mol_name_target) 856 else: 857 mol_names.append(mol_name) 858 res_nums.append(res_num) 859 res_names.append(res_name) 860 spin_nums.append(atom_num) 861 spin_names.append(atom_name) 862 863 # Position vector. 864 spin_cont.pos = pos 865 866 # Add the element. 867 spin_cont.element = element 868 869 # Catch no data. 870 if len(mol_names) == 0: 871 warn(RelaxWarning("No spins matching the '%s' ID string could be found." % spin_id)) 872 return 873 874 # Print out. 875 write_spin_data(file=sys.stdout, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names) 876 877 # Set the number of states for use in the specific analyses. 878 cdp.N = cdp.structure.num_models()

882 """Load the spins from the structural object into the relax data store. 883 884 @keyword spin_id: The molecule, residue, and spin identifier string. 885 @type spin_id: str 886 @keyword str_id: The structure identifier. This can be the file name, model number, or structure number. 887 @type str_id: int or str 888 @keyword from_mols: The list of similar, but not necessarily identical molecules to load spin information from. 889 @type from_mols: list of str or None 890 @keyword mol_name_target: The name of target molecule container, overriding the name of the loaded structures 891 @type mol_name_target: str or None 892 @keyword ave_pos: A flag specifying if the average atom position or the atom position from all loaded structures is loaded into the SpinContainer. 893 @type ave_pos: bool 894 """ 895 896 # Checks. 897 check_pipe() 898 check_structure() 899 900 # The target molecule name must be supplied. 901 if mol_name_target == None: 902 raise RelaxError("The target molecule name must be supplied when specifying multiple molecules to load spins from.") 903 904 # Disallow multiple structural models. 905 if cdp.structure.num_models() != 1: 906 raise RelaxError("Only a single structural model is allowed when specifying multiple molecules to load spins from.") 907 908 # Split up the selection string. 909 if spin_id: 910 mol_token, res_token, spin_token = tokenise(spin_id) 911 if mol_token != None: 912 raise RelaxError("The spin ID string cannot contain molecular information when specifying multiple molecules to load spins from.") 913 914 # Print out. 915 print("Adding the following spins to the relax data store.\n") 916 917 # Initialise the data structures. 918 ids = [] 919 res_nums = {} 920 res_names = {} 921 spin_names = {} 922 positions = {} 923 elements = {} 924 925 # Loop over all target molecules. 926 for mol_name in from_mols: 927 # Add the molecule name as a key for the positions structure, and initialise as a dictionary for the spin IDs. 928 positions[mol_name] = {} 929 930 # Create a new spin ID with the molecule name. 931 new_id = '#' + mol_name 932 if spin_id != None: 933 new_id += spin_id 934 935 # Loop over all atoms of the new spin ID selection. 936 selection = cdp.structure.selection(atom_id=new_id) 937 for res_num, res_name, atom_num, atom_name, element, pos in cdp.structure.atom_loop(selection=selection, str_id=str_id, res_num_flag=True, res_name_flag=True, atom_num_flag=True, atom_name_flag=True, element_flag=True, pos_flag=True, ave=ave_pos): 938 # Remove the '+' regular expression character from the res and atom names. 939 if res_name and search('\+', res_name): 940 res_name = res_name.replace('+', '') 941 if atom_name and search('\+', atom_name): 942 atom_name = atom_name.replace('+', '') 943 944 # Generate a spin ID for the current atom. 945 id = generate_spin_id_unique(mol_name=mol_name_target, res_num=res_num, res_name=res_name, spin_name=atom_name) 946 947 # Store the position info in all cases, collapsing list of lists into single lists when needed. 948 if is_float(pos[0]): 949 positions[mol_name][id] = pos 950 else: 951 positions[mol_name][id] = pos[0] 952 953 # Not a new ID. 954 if id in ids: 955 continue 956 957 # Store the ID, residue, spin, element and position info. 958 ids.append(id) 959 res_nums[id] = res_num 960 res_names[id] = res_name 961 spin_names[id] = atom_name 962 elements[id] = element 963 964 # Catch no data. 965 if len(ids) == 0: 966 warn(RelaxWarning("No spins matching the '%s' ID string could be found." % spin_id)) 967 return 968 969 # Create the spin containers. 970 mol_names2 = [] 971 res_nums2 = [] 972 res_names2 = [] 973 spin_names2 = [] 974 for id in ids: 975 # Fetch the spin. 976 spin_cont = return_spin(id) 977 978 # Create the spin if it does not exist. 979 if spin_cont == None: 980 spin_cont = create_spin(mol_name=mol_name_target, res_num=res_nums[id], res_name=res_names[id], spin_name=spin_names[id]) 981 982 # Position vector. 983 spin_cont.pos = [] 984 for mol_name in from_mols: 985 if id in positions[mol_name]: 986 spin_cont.pos.append(positions[mol_name][id]) 987 else: 988 spin_cont.pos.append(None) 989 990 # Add the element. 991 spin_cont.element = elements[id] 992 993 # Update the structures for the printout. 994 mol_names2.append(mol_name_target) 995 res_nums2.append(res_nums[id]) 996 res_names2.append(res_names[id]) 997 spin_names2.append(spin_names[id]) 998 999 # Print out. 1000 write_spin_data(file=sys.stdout, mol_names=mol_names2, res_nums=res_nums2, res_names=res_names2, spin_names=spin_names2) 1001 1002 # Set the number of states for use in the specific analyses. 1003 cdp.N = len(from_mols)

1007 """Read structures from a Gaussian log file. 1008 1009 1010 @keyword file: The name of the Gaussian log file to read. 1011 @type file: str 1012 @keyword dir: The directory where the Gaussian log file is located. If set to None, then the file will be searched for in the current directory. 1013 @type dir: str or None 1014 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then the molecules will be automatically labelled based on the file name or other information. 1015 @type set_mol_name: None, str, or list of str 1016 @keyword set_model_num: Set the model number of the loaded molecule. 1017 @type set_model_num: None, int, or list of int 1018 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the Gaussian log file does not exist. If False, then a RelaxWarning will be trown instead. 1019 @type fail: bool 1020 @keyword verbosity: The amount of information to print to screen. Zero corresponds to minimal output while higher values increase the amount of output. The default value is 1. 1021 @type verbosity: int 1022 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the Gaussian log file does not exist. 1023 """ 1024 1025 # Test if the current data pipe exists. 1026 check_pipe() 1027 1028 # The file path. 1029 file_path = get_file_path(file, dir) 1030 1031 # Try adding '.log' to the end of the file path, if the file can't be found. 1032 if not access(file_path, F_OK): 1033 file_path_orig = file_path 1034 file_path = file_path + '.log' 1035 1036 # Test if the file exists. 1037 if not access(file_path, F_OK): 1038 if fail: 1039 raise RelaxFileError('Gaussian log', file_path_orig) 1040 else: 1041 warn(RelaxNoPDBFileWarning(file_path)) 1042 return 1043 1044 # Place the structural object into the relax data store. 1045 if not hasattr(cdp, 'structure'): 1046 cdp.structure = Internal() 1047 1048 # Load the structures. 1049 cdp.structure.load_gaussian(file_path, set_mol_name=set_mol_name, set_model_num=set_model_num, verbosity=verbosity)

1053 """The PDB loading function. 1054 1055 @keyword file: The name of the PDB file to read. 1056 @type file: str 1057 @keyword dir: The directory where the PDB file is located. If set to None, then the file will be searched for in the current directory. 1058 @type dir: str or None 1059 @keyword read_mol: The molecule(s) to read from the file, independent of model. The molecules are numbered consecutively from 1. If set to None, then all molecules will be loaded. 1060 @type read_mol: None, int, or list of int 1061 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then the molecules will be automatically labelled based on the file name or other information. 1062 @type set_mol_name: None, str, or list of str 1063 @keyword read_model: The PDB model to extract from the file. If set to None, then all models will be loaded. 1064 @type read_model: None, int, or list of int 1065 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the PDB model numbers will be preserved, if they exist. 1066 @type set_model_num: None, int, or list of int 1067 @keyword alt_loc: The PDB ATOM record 'Alternate location indicator' field value to select which coordinates to use. 1068 @type alt_loc: str or None 1069 @keyword verbosity: The amount of information to print to screen. Zero corresponds to minimal output while higher values increase the amount of output. The default value is 1. 1070 @type verbosity: int 1071 @keyword merge: A flag which if set to True will try to merge the PDB structure into the currently loaded structures. 1072 @type merge: bool 1073 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the PDB file does not exist. If False, then a RelaxWarning will be trown instead. 1074 @type fail: bool 1075 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the PDB file does not exist. 1076 """ 1077 1078 # Test if the current data pipe exists. 1079 check_pipe() 1080 1081 # The file path. 1082 file_path = get_file_path(file, dir) 1083 1084 # Try adding file extensions to the end of the file path, if the file can't be found. 1085 file_path_orig = file_path 1086 if not access(file_path, F_OK): 1087 # List of possible extensions. 1088 for ext in ['.pdb', '.gz', '.pdb.gz', '.bz2', '.pdb.bz2']: 1089 # Add the extension if the file can be found. 1090 if access(file_path+ext, F_OK): 1091 file_path = file_path + ext 1092 1093 # Test if the file exists. 1094 if not access(file_path, F_OK): 1095 if fail: 1096 raise RelaxFileError('PDB', file_path_orig) 1097 else: 1098 warn(RelaxNoPDBFileWarning(file_path)) 1099 return 1100 1101 # Place the internal structural object into the relax data store. 1102 if not hasattr(cdp, 'structure'): 1103 cdp.structure = Internal() 1104 1105 # Load the structures. 1106 cdp.structure.load_pdb(file_path, read_mol=read_mol, set_mol_name=set_mol_name, read_model=read_model, set_model_num=set_model_num, alt_loc=alt_loc, verbosity=verbosity, merge=merge) 1107 1108 # Load into Molmol (if running). 1109 molmol.molmol_obj.open_pdb()

1113 """The XYZ loading function. 1114 1115 1116 @keyword file: The name of the XYZ file to read. 1117 @type file: str 1118 @keyword dir: The directory where the XYZ file is located. If set to None, then the 1119 file will be searched for in the current directory. 1120 @type dir: str or None 1121 @keyword read_mol: The molecule(s) to read from the file, independent of model. 1122 If set to None, then all molecules will be loaded. 1123 @type read_mol: None, int, or list of int 1124 @keyword set_mol_name: Set the names of the molecules which are loaded. If set to None, then 1125 the molecules will be automatically labelled based on the file name or 1126 other information. 1127 @type set_mol_name: None, str, or list of str 1128 @keyword read_model: The XYZ model to extract from the file. If set to None, then all models 1129 will be loaded. 1130 @type read_model: None, int, or list of int 1131 @keyword set_model_num: Set the model number of the loaded molecule. If set to None, then the 1132 XYZ model numbers will be preserved, if they exist. 1133 @type set_model_num: None, int, or list of int 1134 @keyword fail: A flag which, if True, will cause a RelaxError to be raised if the XYZ 1135 file does not exist. If False, then a RelaxWarning will be trown 1136 instead. 1137 @type fail: bool 1138 @keyword verbosity: The amount of information to print to screen. Zero corresponds to 1139 minimal output while higher values increase the amount of output. The 1140 default value is 1. 1141 @type verbosity: int 1142 @raise RelaxFileError: If the fail flag is set, then a RelaxError is raised if the XYZ file 1143 does not exist. 1144 """ 1145 1146 # Test if the current data pipe exists. 1147 check_pipe() 1148 1149 # The file path. 1150 file_path = get_file_path(file, dir) 1151 1152 # Try adding '.xyz' to the end of the file path, if the file can't be found. 1153 if not access(file_path, F_OK): 1154 file_path_orig = file_path 1155 file_path = file_path + '.xyz' 1156 1157 # Test if the file exists. 1158 if not access(file_path, F_OK): 1159 if fail: 1160 raise RelaxFileError('XYZ', file_path_orig) 1161 else: 1162 warn(RelaxNoPDBFileWarning(file_path)) 1163 return 1164 1165 # Place the structural object into the relax data store. 1166 if not hasattr(cdp, 'structure'): 1167 cdp.structure = Internal() 1168 1169 # Load the structures. 1170 cdp.structure.load_xyz(file_path, read_mol=read_mol, set_mol_name=set_mol_name, read_model=read_model, set_model_num=set_model_num, verbosity=verbosity)

1174 """Calculate the RMSD between the loaded models. 1175 1176 The RMSD value will be placed into the current data pipe cdp.structure.rmsd data structure. 1177 1178 1179 @keyword pipes: The data pipes to determine the RMSD for. 1180 @type pipes: None or list of str 1181 @keyword models: The list of models to determine the RMSD for. The number of elements must match the pipes argument. If set to None, then all models will be used. 1182 @type models: None or list of lists of int 1183 @keyword molecules: The list of molecules to determine the RMSD for. The number of elements must match the pipes argument. 1184 @type molecules: None or list of lists of str 1185 @keyword atom_id: The atom identification string of the coordinates of interest. This matches the spin ID string format. 1186 @type atom_id: str or None 1187 @return: The RMSD value. 1188 @rtype: float 1189 """ 1190 1191 # Test if the current data pipe exists. 1192 check_pipe() 1193 1194 # Assemble the structural coordinates. 1195 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 1196 1197 # Calculate the RMSD. 1198 cdp.structure.rmsd = atomic_rmsd(coord, verbosity=1) 1199 1200 # Return the RMSD. 1201 return cdp.structure.rmsd

1205 """Rotate the structural data about the origin by the specified forwards rotation. 1206 1207 @keyword R: The forwards rotation matrix. 1208 @type R: numpy 3D, rank-2 array or a 3x3 list of floats 1209 @keyword origin: The origin of the rotation. If not supplied, the origin will be set to [0, 0, 0]. 1210 @type origin: numpy 3D, rank-1 array or list of len 3 or None 1211 @keyword model: The model to rotate. If None, all models will be rotated. 1212 @type model: int 1213 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 1214 @type atom_id: str or None 1215 @keyword pipe_name: The name of the data pipe containing the structures to translate. This defaults to the current data pipe. 1216 @type pipe_name: None or str 1217 """ 1218 1219 # Defaults. 1220 if pipe_name == None: 1221 pipe_name = cdp_name() 1222 1223 # Checks. 1224 check_pipe(pipe_name) 1225 check_structure(pipe_name) 1226 1227 # Get the data pipe. 1228 dp = get_pipe(pipe_name) 1229 1230 # Set the origin if not supplied. 1231 if origin == None: 1232 origin = [0, 0, 0] 1233 1234 # Convert the args to numpy float data structures. 1235 R = array(R, float64) 1236 origin = array(origin, float64) 1237 1238 # Call the specific code. 1239 selection = dp.structure.selection(atom_id=atom_id) 1240 dp.structure.rotate(R=R, origin=origin, model=model, selection=selection) 1241 1242 # Final printout. 1243 if model != None: 1244 print("Rotated %i atoms of model %i." % (selection.count_atoms(), model)) 1245 else: 1246 print("Rotated %i atoms." % selection.count_atoms())

1250 """Superimpose a set of related, but not identical structures. 1251 1252 @keyword pipes: The data pipes to include in the alignment and superimposition. 1253 @type pipes: None or list of str 1254 @keyword models: The list of models to for each data pipe superimpose. The number of elements must match the pipes argument. If set to None, then all models will be used. 1255 @type models: list of lists of int or None 1256 @keyword molecules: The molecule names to include in the alignment and superimposition. The number of elements must match the pipes argument. 1257 @type molecules: None or list of str 1258 @keyword msa_algorithm: The multiple sequence alignment (MSA) algorithm to use. 1259 @type msa_algorithm: str 1260 @keyword pairwise_algorithm: The pairwise sequence alignment algorithm to use. 1261 @type pairwise_algorithm: str 1262 @keyword matrix: The substitution matrix to use. 1263 @type matrix: str 1264 @keyword gap_open_penalty: The penalty for introducing gaps, as a positive number. 1265 @type gap_open_penalty: float 1266 @keyword gap_extend_penalty: The penalty for extending a gap, as a positive number. 1267 @type gap_extend_penalty: float 1268 @keyword end_gap_open_penalty: The optional penalty for opening a gap at the end of a sequence. 1269 @type end_gap_open_penalty: float 1270 @keyword end_gap_extend_penalty: The optional penalty for extending a gap at the end of a sequence. 1271 @type end_gap_extend_penalty: float 1272 """ 1273 1274 # Assemble the structural objects. 1275 objects, object_names, pipes = assemble_structural_objects(pipes=pipes, models=models, molecules=molecules) 1276 1277 # Assemble the atomic coordinates of all molecules. 1278 ids, object_id_list, model_list, molecule_list, atom_pos, mol_names, res_names, res_nums, atom_names, elements, one_letter_codes, num_mols = assemble_atomic_coordinates(objects=objects, object_names=object_names, molecules=molecules, models=models) 1279 1280 # Convert the residue number data structure. 1281 res_num_list = [] 1282 for mol_index in range(num_mols): 1283 res_num_list.append([]) 1284 for i in range(len(one_letter_codes[mol_index])): 1285 key = list(res_nums[mol_index][i].keys())[0] 1286 res_num_list[mol_index].append(res_nums[mol_index][i][key]) 1287 1288 # MSA. 1289 strings, gaps = msa_general(one_letter_codes, residue_numbers=res_num_list, msa_algorithm=msa_algorithm, pairwise_algorithm=pairwise_algorithm, matrix=matrix, gap_open_penalty=gap_open_penalty, gap_extend_penalty=gap_extend_penalty, end_gap_open_penalty=end_gap_open_penalty, end_gap_extend_penalty=end_gap_extend_penalty) 1290 1291 # Set up the data store object. 1292 if not hasattr(ds, 'sequence_alignments'): 1293 ds.sequence_alignments = Sequence_alignments() 1294 1295 # Set some unused arguments to None for storage. 1296 if msa_algorithm == 'residue number': 1297 pairwise_algorithm = None 1298 matrix = None 1299 gap_open_penalty = None 1300 gap_extend_penalty = None 1301 end_gap_open_penalty = None 1302 end_gap_extend_penalty = None 1303 1304 # Store the alignment. 1305 ds.sequence_alignments.add(object_ids=object_id_list, models=model_list, molecules=molecule_list, sequences=one_letter_codes, strings=strings, gaps=gaps, msa_algorithm=msa_algorithm, pairwise_algorithm=pairwise_algorithm, matrix=matrix, gap_open_penalty=gap_open_penalty, gap_extend_penalty=gap_extend_penalty, end_gap_open_penalty=end_gap_open_penalty, end_gap_extend_penalty=end_gap_extend_penalty)

1322 """Generator function for looping over all internal structural objects, models and molecules. 1323 1324 @keyword pipes: The data pipes to loop over. 1325 @type pipes: None or list of str 1326 @keyword models: The list of models for each data pipe. The number of elements must match the pipes argument. If set to None, then all models will be used. 1327 @type models: None or list of lists of int 1328 @keyword molecules: The list of molecules for each data pipe. The number of elements must match the pipes argument. 1329 @type molecules: None or list of lists of str 1330 @keyword atom_id: The molecule, residue, and atom identifier string of the coordinates of interest. This matches the spin ID string format. 1331 @type atom_id: None or str 1332 @return: The data pipe index, model number, and molecule name. 1333 @rtype: int, int or None, str 1334 """ 1335 1336 # The data pipes to use. 1337 if pipes == None: 1338 pipes = [cdp_name()] 1339 num_pipes = len(pipes) 1340 1341 # Checks. 1342 for pipe in pipes: 1343 check_pipe(pipe) 1344 1345 # Check the models and molecules arguments. 1346 if models != None: 1347 if len(models) != num_pipes: 1348 raise RelaxError("The %i elements of the models argument does not match the %i data pipes." % (len(models), num_pipes)) 1349 if molecules != None: 1350 if len(molecules) != num_pipes: 1351 raise RelaxError("The %i elements of the molecules argument does not match the %i data pipes." % (len(molecules), num_pipes)) 1352 1353 # Assemble the structural objects. 1354 objects = [] 1355 for pipe_index in range(len(pipes)): 1356 dp = get_pipe(pipes[pipe_index]) 1357 objects.append(dp.structure) 1358 1359 # Call the library method to do all of the work. 1360 for pipe_index, model_num, mol_name in loop_coord_structures(objects=objects, models=models, molecules=molecules, atom_id=atom_id): 1361 yield pipe_index, model_num, mol_name

1365 """Superimpose a set of structures. 1366 1367 @keyword pipes: The data pipes to include in the alignment and superimposition. 1368 @type pipes: None or list of str 1369 @keyword models: The list of models to for each data pipe superimpose. The number of elements must match the pipes argument. If set to None, then all models will be used. 1370 @type models: list of lists of int or None 1371 @keyword molecules: The molecule names to include in the alignment and superimposition. The number of elements must match the pipes argument. 1372 @type molecules: None or list of str 1373 @keyword atom_id: The molecule, residue, and atom identifier string. This matches the spin ID string format. 1374 @type atom_id: str or None 1375 @keyword displace_id: The atom ID string for restricting the displacement to a subset of all atoms. If not set, then all atoms will be translated and rotated. This can be a list of atom IDs with each element corresponding to one of the structures. 1376 @type displace_id: None, str, or list of str 1377 @keyword method: The superimposition method. It must be one of 'fit to mean' or 'fit to first'. 1378 @type method: str 1379 @keyword centre_type: The type of centre to superimpose over. This can either be the standard centroid superimposition or the CoM could be used instead. 1380 @type centre_type: str 1381 @keyword centroid: An alternative position of the centroid to allow for different superpositions, for example of pivot point motions. 1382 @type centroid: list of float or numpy rank-1, 3D array 1383 """ 1384 1385 # Check the method. 1386 allowed = ['fit to mean', 'fit to first'] 1387 if method not in allowed: 1388 raise RelaxError("The superimposition method '%s' is unknown. It must be one of %s." % (method, allowed)) 1389 1390 # Check the type. 1391 allowed = ['centroid', 'CoM'] 1392 if centre_type not in allowed: 1393 raise RelaxError("The superimposition centre type '%s' is unknown. It must be one of %s." % (centre_type, allowed)) 1394 1395 # The data pipes to use. 1396 if pipes == None: 1397 pipes = [cdp_name()] 1398 1399 # Assemble the structural coordinates. 1400 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 1401 1402 # Catch missing data. 1403 if len(coord[0]) == 0: 1404 raise RelaxError("No common atoms could be found between the structures.") 1405 1406 # The different algorithms. 1407 if method == 'fit to mean': 1408 T, R, pivot = fit_to_mean(models=list(range(len(ids))), coord=coord, centre_type=centre_type, elements=elements, centroid=centroid) 1409 elif method == 'fit to first': 1410 T, R, pivot = fit_to_first(models=list(range(len(ids))), coord=coord, centre_type=centre_type, elements=elements, centroid=centroid) 1411 1412 # Loop over all pipes, models, and molecules. 1413 i = 0 1414 for pipe_index, model_num, mol_name in structure_loop(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id): 1415 # The current displacement ID. 1416 curr_displace_id = None 1417 if isinstance(displace_id, str): 1418 curr_displace_id = displace_id 1419 elif isinstance(displace_id, list): 1420 if len(displace_id) <= i: 1421 raise RelaxError("Not enough displacement ID strings have been provided.") 1422 curr_displace_id = displace_id[i] 1423 1424 # Add the molecule name to the displacement ID if required. 1425 id = curr_displace_id 1426 if id == None or (mol_name and not search('#', id)): 1427 if curr_displace_id == None: 1428 id = '#%s' % mol_name 1429 elif search('#', curr_displace_id): 1430 id = curr_displace_id 1431 else: 1432 id = '#%s%s' % (mol_name, curr_displace_id) 1433 1434 # Translate the molecule first (the rotational pivot is defined in the first model). 1435 translate(T=T[i], model=model_num, pipe_name=pipes[pipe_index], atom_id=id) 1436 1437 # Rotate the molecule. 1438 rotate(R=R[i], origin=pivot[i], model=model_num, pipe_name=pipes[pipe_index], atom_id=id) 1439 1440 # Increment the index. 1441 i += 1

1445 """Shift the structural data by the specified translation vector. 1446 1447 @keyword T: The translation vector. 1448 @type T: numpy rank-1, 3D array or list of float 1449 @keyword model: The model to translate. If None, all models will be rotated. 1450 @type model: int or None 1451 @keyword atom_id: The molecule, residue, and atom identifier string. Only atoms matching this selection will be used. 1452 @type atom_id: str or None 1453 @keyword pipe_name: The name of the data pipe containing the structures to translate. This defaults to the current data pipe. 1454 @type pipe_name: None or str 1455 """ 1456 1457 # Defaults. 1458 if pipe_name == None: 1459 pipe_name = cdp_name() 1460 1461 # Checks. 1462 check_pipe(pipe_name) 1463 check_structure(pipe_name) 1464 1465 # Get the data pipe. 1466 dp = get_pipe(pipe_name) 1467 1468 # Convert the args to numpy float data structures. 1469 T = array(T, float64) 1470 1471 # Call the specific code. 1472 selection = dp.structure.selection(atom_id=atom_id) 1473 dp.structure.translate(T=T, model=model, selection=selection) 1474 1475 # Final printout. 1476 if model != None: 1477 print("Translated %i atoms of model %i." % (selection.count_atoms(), model)) 1478 else: 1479 print("Translated %i atoms." % selection.count_atoms())

1483 """Extract the bond vectors from the loaded structures and store them in the spin container. 1484 1485 @keyword spin_id1: The spin identifier string of the first spin of the pair. 1486 @type spin_id1: str 1487 @keyword spin_id2: The spin identifier string of the second spin of the pair. 1488 @type spin_id2: str 1489 @keyword model: The model to extract the vector from. If None, all vectors will be extracted. 1490 @type model: str 1491 @keyword verbosity: The higher the value, the more information is printed to screen. 1492 @type verbosity: int 1493 @keyword ave: A flag which if True will cause the average of all vectors to be extracted. 1494 @type ave: bool 1495 @keyword unit: A flag which if True will cause the function to calculate the unit vectors. 1496 @type unit: bool 1497 """ 1498 1499 # Checks. 1500 check_pipe() 1501 check_structure() 1502 check_mol_res_spin_data() 1503 1504 # Print out. 1505 if verbosity: 1506 # Number of models. 1507 num_models = cdp.structure.num_models() 1508 1509 # Multiple models loaded. 1510 if num_models > 1: 1511 if model: 1512 print("Extracting vectors for model '%s'." % model) 1513 else: 1514 print("Extracting vectors for all %s models." % num_models) 1515 if ave: 1516 print("Averaging all vectors.") 1517 1518 # Single model loaded. 1519 else: 1520 print("Extracting vectors from the single model.") 1521 1522 # Unit vectors. 1523 if unit: 1524 print("Calculating the unit vectors.") 1525 1526 # Loop over the spins. 1527 no_vectors = True 1528 for spin, mol_name, res_num, res_name in spin_loop(selection=spin_id, full_info=True): 1529 # Skip deselected spins. 1530 if not spin.select: 1531 continue 1532 1533 # The spin identification string. The residue name and spin num is not included to allow molecules with point mutations to be used as different models. 1534 id = generate_spin_id_unique(res_num=res_num, res_name=None, spin_name=spin.name, spin_num=spin.num) 1535 1536 # Test that the spin number or name are set (one or both are essential for the identification of the atom). 1537 if spin.num == None and spin.name == None: 1538 warn(RelaxWarning("Either the spin number or name must be set for the spin " + repr(id) + " to identify the corresponding atom in the molecule.")) 1539 continue 1540 1541 # The bond vector already exists. 1542 if hasattr(spin, 'vector'): 1543 obj = getattr(spin, 'vector') 1544 if obj != None: 1545 warn(RelaxWarning("The bond vector for the spin " + repr(id) + " already exists.")) 1546 continue 1547 1548 # Get the bond info. 1549 bond_vectors, attached_name, warnings = cdp.structure.bond_vectors(attached_atom=attached, model_num=model, res_num=res_num, spin_name=spin.name, spin_num=spin.num, return_name=True, return_warnings=True) 1550 id2 = generate_spin_id_unique(res_num=res_num, res_name=None, spin_name=spin.name) 1551 1552 # No attached atom. 1553 if not bond_vectors: 1554 # Warning messages. 1555 if warnings: 1556 warn(RelaxWarning(warnings + " (atom ID " + repr(id) + ").")) 1557 1558 # Skip the spin. 1559 continue 1560 1561 # Set the attached atom name. 1562 if not hasattr(spin, 'attached_atom'): 1563 spin.attached_atom = attached_name 1564 elif spin.attached_atom != attached_name: 1565 raise RelaxError("The " + repr(spin.attached_atom) + " atom already attached to the spin does not match the attached atom " + repr(attached_name) + ".") 1566 1567 # Initialise the average vector. 1568 if ave: 1569 ave_vector = zeros(3, float64) 1570 1571 # Loop over the individual vectors. 1572 for i in range(len(bond_vectors)): 1573 # Unit vector. 1574 if unit: 1575 # Normalisation factor. 1576 norm_factor = norm(bond_vectors[i]) 1577 1578 # Test for zero length. 1579 if norm_factor == 0.0: 1580 warn(RelaxZeroVectorWarning(spin_id1=id, spin_id2=id2)) 1581 1582 # Calculate the normalised vector. 1583 else: 1584 bond_vectors[i] = bond_vectors[i] / norm_factor 1585 1586 # Sum the vectors. 1587 if ave: 1588 ave_vector = ave_vector + bond_vectors[i] 1589 1590 # Average. 1591 if ave: 1592 vector = ave_vector / float(len(bond_vectors)) 1593 else: 1594 vector = bond_vectors 1595 1596 # Convert to a single vector if needed. 1597 if len(vector) == 1: 1598 vector = vector[0] 1599 1600 # Set the vector. 1601 setattr(spin, 'vector', vector) 1602 1603 # We have a vector! 1604 no_vectors = False 1605 1606 # Print out of modified spins. 1607 if verbosity: 1608 # The number of vectors. 1609 num = len(bond_vectors) 1610 plural = 's' 1611 if num == 1: 1612 plural = '' 1613 1614 if spin.name: 1615 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.name, attached_name, plural, id)) 1616 else: 1617 print("Extracted %s %s-%s vector%s for the spin '%s'." % (num, spin.num, attached_name, plural, id)) 1618 1619 # Right, catch the problem of missing vectors to prevent massive user confusion! 1620 if no_vectors: 1621 raise RelaxError("No vectors could be extracted.")

1625 """Create a PDB representation of the motion between a set of models. 1626 1627 This will create a PDB file containing the atoms of all models, with identical atoms links using CONECT records. This function only supports the internal structural object. 1628 1629 1630 @keyword pipes: The data pipes to generate the web between. 1631 @type pipes: None or list of str 1632 @keyword models: The list of models to generate the web between. The number of elements must match the pipes argument. If set to None, then all models will be used. 1633 @type models: None or list of lists of int 1634 @keyword molecules: The list of molecules to generate the web between. The number of elements must match the pipes argument. 1635 @type molecules: None or list of lists of str 1636 @keyword atom_id: The atom identification string of the coordinates of interest. This matches the spin ID string format. 1637 @type atom_id: str or None 1638 @keyword file: The name of the PDB file to write. 1639 @type file: str 1640 @keyword dir: The directory where the PDB file will be placed. If set to None, then the file will be placed in the current directory. 1641 @type dir: str or None 1642 @keyword force: The force flag which if True will cause the file to be overwritten. 1643 @type force: bool 1644 """ 1645 1646 # Checks. 1647 check_pipe() 1648 check_structure() 1649 1650 # Assemble the structural coordinates. 1651 coord, ids, mol_names, res_names, res_nums, atom_names, elements = assemble_structural_coordinates(pipes=pipes, models=models, molecules=molecules, atom_id=atom_id) 1652 1653 # Check that more than one structure is present. 1654 if not len(coord) > 1: 1655 raise RelaxError("Two or more structures are required.") 1656 1657 # Initialise the structural object. 1658 web = Internal() 1659 1660 # Loop over the atoms. 1661 for atom_index in range(len(atom_names)): 1662 # Loop over the structures. 1663 for struct_index in range(len(ids)): 1664 # Add the atom. 1665 web.add_atom(mol_name=mol_names[atom_index], atom_name=atom_names[atom_index], res_name=res_names[atom_index], res_num=res_nums[atom_index], pos=coord[struct_index, atom_index], element=elements[atom_index], sort=False) 1666 1667 # Loop over the structures again, this time twice. 1668 for k in range(len(ids)): 1669 for l in range(len(ids)): 1670 # Skip identical atoms. 1671 if k == l: 1672 continue 1673 1674 # The atom index. 1675 index1 = atom_index*len(ids) + k 1676 index2 = atom_index*len(ids) + l 1677 1678 # Connect to the previous atoms. 1679 web.connect_atom(mol_name=mol_names[atom_index], index1=index1, index2=index2) 1680 1681 # Append the PDB extension if needed. 1682 if isinstance(file, str): 1683 # The file path. 1684 file = get_file_path(file, dir) 1685 1686 # Add '.pdb' to the end of the file path if it isn't there yet. 1687 if not search(".pdb$", file): 1688 file += '.pdb' 1689 1690 # Open the file for writing. 1691 file = open_write_file(file, force=force) 1692 1693 # Write the structure. 1694 web.write_pdb(file)

1698 """The PDB writing function. 1699 1700 @keyword file: The name of the PDB file to write. This can also be a file instance. 1701 @type file: str or file instance 1702 @keyword dir: The directory where the PDB file will be placed. If set to None, then the file will be placed in the current directory. 1703 @type dir: str or None 1704 @keyword model_num: The model to place into the PDB file. If not supplied, then all models will be placed into the file. 1705 @type model_num: None or int 1706 @keyword compress_type: The compression type. The integer values correspond to the compression type: 0, no compression; 1, Bzip2 compression; 2, Gzip compression. 1707 @type compress_type: int 1708 @keyword force: The force flag which if True will cause the file to be overwritten. 1709 @type force: bool 1710 """ 1711 1712 # Test if the current data pipe exists. 1713 check_pipe() 1714 1715 # Check if the structural object exists. 1716 if not hasattr(cdp, 'structure'): 1717 raise RelaxError("No structural data is present in the current data pipe.") 1718 1719 # Path handling. 1720 if isinstance(file, str): 1721 # The file path. 1722 file = get_file_path(file, dir) 1723 1724 # Add '.pdb' to the end of the file path if it isn't there yet. 1725 if not search(".pdb$", file): 1726 file = file + '.pdb' 1727 1728 # Open the file for writing. 1729 file = open_write_file(file, compress_type=compress_type, force=force) 1730 1731 # Write the structures. 1732 cdp.structure.write_pdb(file, model_num=model_num)

Source Code for Module pipe_control.structure.main