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

Source Code for Module pipe_control.pcs

   1  ############################################################################### 
   2  #                                                                             # 
   3  # Copyright (C) 2003-2016 Edward d'Auvergne                                   # 
   4  # Copyright (C) 2006 Chris MacRaild                                           # 
   5  #                                                                             # 
   6  # This file is part of the program relax (http://www.nmr-relax.com).          # 
   7  #                                                                             # 
   8  # This program is free software: you can redistribute it and/or modify        # 
   9  # it under the terms of the GNU General Public License as published by        # 
  10  # the Free Software Foundation, either version 3 of the License, or           # 
  11  # (at your option) any later version.                                         # 
  12  #                                                                             # 
  13  # This program is distributed in the hope that it will be useful,             # 
  14  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
  15  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
  16  # GNU General Public License for more details.                                # 
  17  #                                                                             # 
  18  # You should have received a copy of the GNU General Public License           # 
  19  # along with this program.  If not, see <http://www.gnu.org/licenses/>.       # 
  20  #                                                                             # 
  21  ############################################################################### 
  22   
  23  # Module docstring. 
  24  """Module for the manipulation of pseudo-contact shift data.""" 
  25   
  26  # Python module imports. 
  27  from copy import deepcopy 
  28  from math import ceil, floor, pi, sqrt 
  29  from numpy import array, eye, float64, int32, ones, std, zeros 
  30  from numpy.linalg import norm 
  31  from numpy.random import multivariate_normal 
  32  import sys 
  33  from warnings import warn 
  34   
  35  # relax module imports. 
  36  from lib.alignment.pcs import ave_pcs_tensor, pcs_tensor 
  37  from lib.check_types import is_float 
  38  from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoPdbError, RelaxNoPCSError, RelaxNoSequenceError 
  39  from lib.io import open_write_file, write_data 
  40  from lib.periodic_table import periodic_table 
  41  from lib.physical_constants import pcs_constant 
  42  from lib.plotting.api import write_xy_data, write_xy_header 
  43  from lib.sequence import read_spin_data, write_spin_data 
  44  from lib.warnings import RelaxWarning, RelaxNoSpinWarning 
  45  from pipe_control import pipes 
  46  from pipe_control.align_tensor import get_tensor_index, get_tensor_object, opt_uses_align_data, opt_uses_tensor 
  47  from pipe_control.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, is_pseudoatom, return_spin, spin_loop 
  48  from pipe_control.pipes import check_pipe 
  49   
  50   
51 -def back_calc(align_id=None):
52 """Back calculate the PCS from the alignment tensor. 53 54 @keyword align_id: The alignment tensor ID string. 55 @type align_id: str 56 """ 57 58 # Check the pipe setup. 59 check_pipe_setup(pcs_id=align_id, sequence=True, N=True, tensors=True, paramag_centre=True) 60 61 # Convert the align IDs to an array, or take all IDs. 62 if align_id: 63 align_ids = [align_id] 64 else: 65 align_ids = cdp.align_ids 66 67 # Add the ID to the PCS IDs, if needed. 68 for align_id in align_ids: 69 # Init. 70 if not hasattr(cdp, 'pcs_ids'): 71 cdp.pcs_ids = [] 72 73 # Add the ID. 74 if align_id not in cdp.pcs_ids: 75 cdp.pcs_ids.append(align_id) 76 77 # The weights. 78 weights = ones(cdp.N, float64) / cdp.N 79 80 # Unit vector data structure init. 81 unit_vect = zeros((cdp.N, 3), float64) 82 83 # Loop over the spins. 84 count = 0 85 for spin in spin_loop(): 86 # Skip spins with no position. 87 if not hasattr(spin, 'pos'): 88 continue 89 90 # Atom positions. 91 pos = spin.pos 92 if type(pos[0]) in [float, float64]: 93 pos = [pos] * cdp.N 94 95 # Loop over the alignments. 96 for id in align_ids: 97 # Vectors. 98 vect = zeros((cdp.N, 3), float64) 99 r = zeros(cdp.N, float64) 100 dj = zeros(cdp.N, float64) 101 for c in range(cdp.N): 102 # The vector. 103 vect[c] = pos[c] - cdp.paramagnetic_centre 104 105 # The length. 106 r[c] = norm(vect[c]) 107 108 # Normalise (only if the vector has length). 109 if r[c]: 110 vect[c] = vect[c] / r[c] 111 112 # Calculate the PCS constant. 113 dj[c] = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H'), r[c]/1e10) 114 115 # Initialise if necessary. 116 if not hasattr(spin, 'pcs_bc'): 117 spin.pcs_bc = {} 118 119 # Calculate the PCSs (in ppm). 120 spin.pcs_bc[id] = ave_pcs_tensor(dj, vect, cdp.N, cdp.align_tensors[get_tensor_index(align_id=id)].A, weights=weights) * 1e6 121 122 # Increment the counter. 123 count += 1 124 125 # No PCSs calculated. 126 if not count: 127 warn(RelaxWarning("No PCSs have been back calculated, probably due to missing spin position information."))
128 129
130 -def centre(pos=None, atom_id=None, pipe=None, verbosity=1, ave_pos=False, force=False):
131 """Specify the atom in the loaded structure corresponding to the paramagnetic centre. 132 133 @keyword pos: The atomic position. If set, the atom_id string will be ignored. 134 @type pos: list of float 135 @keyword atom_id: The atom identification string. 136 @type atom_id: str 137 @keyword pipe: An alternative data pipe to extract the paramagnetic centre from. 138 @type pipe: None or str 139 @keyword verbosity: The amount of information to print out. The bigger the number, the more information. 140 @type verbosity: int 141 @keyword ave_pos: A flag which if True causes the atomic positions from multiple models to be averaged. 142 @type ave_pos: bool 143 @keyword force: A flag which if True will cause the current PCS centre to be overwritten. 144 """ 145 146 # The data pipe. 147 if pipe == None: 148 pipe = pipes.cdp_name() 149 150 # Check the pipe setup. 151 check_pipe_setup(pipe=pipe) 152 153 # Get the data pipes. 154 source_dp = pipes.get_pipe(pipe) 155 156 # Test if the structure has been loaded. 157 if not hasattr(source_dp, 'structure'): 158 raise RelaxNoPdbError 159 160 # Test the centre has already been set. 161 if not force and hasattr(cdp, 'paramagnetic_centre'): 162 raise RelaxError("The paramagnetic centre has already been set to the coordinates " + repr(cdp.paramagnetic_centre) + ".") 163 164 # Position is supplied. 165 if pos != None: 166 centre = array(pos) 167 num_pos = 1 168 full_pos_list = [] 169 170 # Position from a loaded structure. 171 else: 172 # Get the positions. 173 centre = zeros(3, float64) 174 full_pos_list = [] 175 num_pos = 0 176 for spin, spin_id in spin_loop(atom_id, pipe=pipe, return_id=True): 177 # No atomic positions. 178 if not hasattr(spin, 'pos'): 179 continue 180 181 # Spin position list. 182 if isinstance(spin.pos[0], float) or isinstance(spin.pos[0], float64): 183 pos_list = [spin.pos] 184 else: 185 pos_list = spin.pos 186 187 # Loop over the model positions. 188 for pos in pos_list: 189 full_pos_list.append(pos) 190 centre = centre + array(pos) 191 num_pos = num_pos + 1 192 193 # No positional information! 194 if not num_pos: 195 raise RelaxError("No positional information could be found for the spin '%s'." % atom_id) 196 197 # Averaging. 198 centre = centre / float(num_pos) 199 200 # Print out. 201 if verbosity: 202 print("Paramagnetic centres located at:") 203 for pos in full_pos_list: 204 print(" [%8.3f, %8.3f, %8.3f]" % (pos[0], pos[1], pos[2])) 205 print("\nAverage paramagnetic centre located at:") 206 print(" [%8.3f, %8.3f, %8.3f]" % (centre[0], centre[1], centre[2])) 207 208 # Set the centre (place it into the current data pipe). 209 if ave_pos: 210 if verbosity: 211 print("\nUsing the average paramagnetic position.") 212 cdp.paramagnetic_centre = centre 213 else: 214 if verbosity: 215 print("\nUsing all paramagnetic positions.") 216 cdp.paramagnetic_centre = full_pos_list
217 218
219 -def check_pipe_setup(pipe=None, pcs_id=None, sequence=False, N=False, tensors=False, pcs=False, paramag_centre=False):
220 """Check that the current data pipe has been setup sufficiently. 221 222 @keyword pipe: The data pipe to check, defaulting to the current pipe. 223 @type pipe: None or str 224 @keyword pcs_id: The PCS ID string to check for in cdp.pcs_ids. 225 @type pcs_id: None or str 226 @keyword sequence: A flag which when True will invoke the sequence data check. 227 @type sequence: bool 228 @keyword N: A flag which if True will check that cdp.N is set. 229 @type N: bool 230 @keyword tensors: A flag which if True will check that alignment tensors exist. 231 @type tensors: bool 232 @keyword pcs: A flag which if True will check that PCSs exist. 233 @type pcs: bool 234 @keyword paramag_centre: A flag which if True will check that the paramagnetic centre has been set. 235 @type paramag_centre: bool 236 """ 237 238 # The data pipe. 239 if pipe == None: 240 pipe = pipes.cdp_name() 241 242 # Get the data pipe. 243 dp = pipes.get_pipe(pipe) 244 245 # Test if the current data pipe exists. 246 check_pipe(pipe) 247 248 # Test if sequence data exists. 249 if sequence and not exists_mol_res_spin_data(pipe): 250 raise RelaxNoSequenceError(pipe) 251 252 # Check for dp.N. 253 if N and not hasattr(dp, 'N'): 254 raise RelaxError("The number of states N has not been set.") 255 256 # Check that alignment tensors are present. 257 if tensors and (not hasattr(dp, 'align_tensors') or len(dp.align_tensors) == 0): 258 raise RelaxNoTensorError('alignment') 259 260 # Test for the alignment ID. 261 if pcs_id and (not hasattr(dp, 'align_ids') or pcs_id not in dp.align_ids): 262 raise RelaxNoAlignError(pcs_id, pipe) 263 264 # Test if PCS data exists. 265 if pcs and not hasattr(dp, 'align_ids'): 266 raise RelaxNoAlignError() 267 if pcs and not hasattr(dp, 'pcs_ids'): 268 raise RelaxNoPCSError() 269 elif pcs and pcs_id and pcs_id not in dp.pcs_ids: 270 raise RelaxNoPCSError(pcs_id) 271 272 # Test if the paramagnetic centre is set. 273 if paramag_centre and not hasattr(cdp, 'paramagnetic_centre'): 274 raise RelaxError("The paramagnetic centre has not been defined.")
275 276
277 -def copy(pipe_from=None, pipe_to=None, align_id=None, back_calc=True):
278 """Copy the PCS data from one data pipe to another. 279 280 @keyword pipe_from: The data pipe to copy the PCS data from. This defaults to the current data pipe. 281 @type pipe_from: str 282 @keyword pipe_to: The data pipe to copy the PCS data to. This defaults to the current data pipe. 283 @type pipe_to: str 284 @keyword align_id: The alignment ID string. 285 @type align_id: str 286 @keyword back_calc: A flag which if True will cause any back-calculated RDCs present to also be copied with the real values and errors. 287 @type back_calc: bool 288 """ 289 290 # Defaults. 291 if pipe_from == None and pipe_to == None: 292 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 293 elif pipe_from == None: 294 pipe_from = pipes.cdp_name() 295 elif pipe_to == None: 296 pipe_to = pipes.cdp_name() 297 298 # Check the pipe setup. 299 check_pipe_setup(pipe=pipe_from, pcs_id=align_id, sequence=True, pcs=True) 300 check_pipe_setup(pipe=pipe_to, sequence=True) 301 302 # Get the data pipes. 303 dp_from = pipes.get_pipe(pipe_from) 304 dp_to = pipes.get_pipe(pipe_to) 305 306 # The IDs. 307 if align_id == None: 308 align_ids = dp_from.align_ids 309 else: 310 align_ids = [align_id] 311 312 # Init target pipe global structures. 313 if not hasattr(dp_to, 'align_ids'): 314 dp_to.align_ids = [] 315 if not hasattr(dp_to, 'pcs_ids'): 316 dp_to.pcs_ids = [] 317 318 # Loop over the align IDs. 319 for align_id in align_ids: 320 # Printout. 321 print("\nCoping PCSs for the alignment ID '%s'." % align_id) 322 323 # Copy the global data. 324 if align_id not in dp_to.align_ids and align_id not in dp_to.align_ids: 325 dp_to.align_ids.append(align_id) 326 if align_id in dp_from.pcs_ids and align_id not in dp_to.pcs_ids: 327 dp_to.pcs_ids.append(align_id) 328 329 # Spin loop. 330 data = [] 331 for spin_from, spin_id in spin_loop(return_id=True, skip_desel=True, pipe=pipe_from): 332 # Find the matching spin container in the target data pipe. 333 spin_to = return_spin(spin_id=spin_id, pipe=pipe_to) 334 335 # No matching spin container. 336 if spin_to == None: 337 warn(RelaxWarning("The spin container for the spin '%s' cannot be found in the target data pipe." % spin_id)) 338 continue 339 340 # No data or errors. 341 if (not hasattr(spin_from, 'pcs') or not align_id in spin_from.pcs) and (not hasattr(spin_from, 'pcs_err') or not align_id in spin_from.pcs_err): 342 continue 343 344 # Initialise the spin data if necessary. 345 if hasattr(spin_from, 'pcs') and not hasattr(spin_to, 'pcs'): 346 spin_to.pcs = {} 347 if back_calc and hasattr(spin_from, 'pcs_bc') and not hasattr(spin_to, 'pcs_bc'): 348 spin_to.pcs_bc = {} 349 if hasattr(spin_from, 'pcs_err') and not hasattr(spin_to, 'pcs_err'): 350 spin_to.pcs_err = {} 351 352 # Copy the value and error from pipe_from. 353 value = None 354 error = None 355 value_bc = None 356 if hasattr(spin_from, 'pcs'): 357 value = spin_from.pcs[align_id] 358 spin_to.pcs[align_id] = value 359 if back_calc and hasattr(spin_from, 'pcs_bc'): 360 value_bc = spin_from.pcs_bc[align_id] 361 spin_to.pcs_bc[align_id] = value_bc 362 if hasattr(spin_from, 'pcs_err'): 363 error = spin_from.pcs_err[align_id] 364 spin_to.pcs_err[align_id] = error 365 366 # Append the data for printout. 367 data.append([spin_id]) 368 if is_float(value): 369 data[-1].append("%20.15f" % value) 370 else: 371 data[-1].append("%20s" % value) 372 if back_calc: 373 if is_float(value_bc): 374 data[-1].append("%20.15f" % value_bc) 375 else: 376 data[-1].append("%20s" % value_bc) 377 if is_float(error): 378 data[-1].append("%20.15f" % error) 379 else: 380 data[-1].append("%20s" % error) 381 382 # Printout. 383 print("The following PCSs have been copied:\n") 384 if back_calc: 385 write_data(out=sys.stdout, headings=["Spin_ID", "Value", "Back-calculated", "Error"], data=data) 386 else: 387 write_data(out=sys.stdout, headings=["Spin_ID", "Value", "Error"], data=data)
388 389
390 -def corr_plot(format=None, title=None, subtitle=None, file=None, dir=None, force=False):
391 """Generate a correlation plot of the measured vs. back-calculated PCSs. 392 393 @keyword format: The format for the plot file. The following values are accepted: 'grace', a Grace plot; None, a plain text file. 394 @type format: str or None 395 @keyword title: The title for the plot, overriding the default. 396 @type title: None or str 397 @keyword subtitle: The subtitle for the plot, overriding the default. 398 @type subtitle: None or str 399 @keyword file: The file name or object to write to. 400 @type file: str or file object 401 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 402 @type dir: str 403 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 404 @type force: bool 405 """ 406 407 # Check the pipe setup. 408 check_pipe_setup(sequence=True) 409 410 # Does PCS data exist? 411 if not hasattr(cdp, 'pcs_ids') or not cdp.pcs_ids: 412 warn(RelaxWarning("No PCS data exists, skipping file creation.")) 413 return 414 415 # Open the file for writing. 416 file = open_write_file(file, dir, force) 417 418 # Init. 419 data = [] 420 orig_title = title 421 if orig_title == None: 422 title = "PCS correlation plot" 423 axis_labels = ["Back-calculated PCS (ppm)", "Measured PCS (ppm)"] 424 425 # The diagonal. 426 data.append([[-100, -100, 0], [100, 100, 0]]) 427 428 # The spin types. 429 types = [] 430 for spin in spin_loop(): 431 if not hasattr(spin, 'element'): 432 if None not in types: 433 types.append(None) 434 elif spin.element not in types: 435 types.append(spin.element) 436 437 # Loop over the PCS data. 438 min_pcs = 1e100 439 max_pcs = -1e100 440 for align_id in cdp.pcs_ids: 441 # Loop over the spin types. 442 for i in range(len(types)): 443 # Append a new list for this alignment. 444 data.append([]) 445 446 # Errors present? 447 err_flag = False 448 for spin in spin_loop(): 449 # Skip deselected spins. 450 if not spin.select: 451 continue 452 453 # Error present. 454 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 455 err_flag = True 456 break 457 458 # Loop over the spins. 459 for spin, spin_id in spin_loop(return_id=True): 460 # Skip deselected spins. 461 if not spin.select: 462 continue 463 464 # Incorrect spin type. 465 if hasattr(spin, 'element') and spin.element != types[i]: 466 continue 467 468 # Skip if data is missing. 469 if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc') or not align_id in spin.pcs or not align_id in spin.pcs_bc: 470 continue 471 if spin.pcs[align_id] == None or spin.pcs_bc[align_id] == None: 472 continue 473 474 # Append the data. 475 data[-1].append([spin.pcs_bc[align_id], spin.pcs[align_id]]) 476 477 # The minimum of all data sets. 478 if spin.pcs[align_id] < min_pcs: 479 min_pcs = spin.pcs[align_id] 480 if spin.pcs_bc[align_id] < min_pcs: 481 min_pcs = spin.pcs_bc[align_id] 482 483 # The maximum of all data sets. 484 if spin.pcs[align_id] > max_pcs: 485 max_pcs = spin.pcs[align_id] 486 if spin.pcs_bc[align_id] > max_pcs: 487 max_pcs = spin.pcs_bc[align_id] 488 489 # Errors. 490 if err_flag: 491 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 492 data[-1][-1].append(spin.pcs_err[align_id]) 493 else: 494 data[-1][-1].append(None) 495 496 # Label. 497 data[-1][-1].append(spin_id) 498 499 # The data size. 500 size = len(data) 501 502 # Round the data limits. 503 max_pcs = float(ceil(max_pcs)) 504 min_pcs = float(floor(min_pcs)) 505 506 # Only one data set. 507 data = [data] 508 509 # Graph type. 510 if err_flag: 511 graph_type = 'xydy' 512 else: 513 graph_type = 'xy' 514 515 # Grace file. 516 if format == 'grace': 517 # The set names. 518 set_names = [None] 519 for i in range(len(cdp.pcs_ids)): 520 for j in range(len(types)): 521 set_names.append("%s (%s)" % (cdp.pcs_ids[i], types[j])) 522 523 # The header. 524 write_xy_header(format=format, file=file, title=title, subtitle=subtitle, world=[[min_pcs, min_pcs, max_pcs, max_pcs]], sets=[size], set_names=[set_names], linestyle=[[2]+[0]*size], data_type=['pcs_bc', 'pcs'], axis_labels=[axis_labels], tick_major_spacing=[[1, 1]], tick_minor_count=[[9, 9]], legend_pos=[[1, 0.5]]) 525 526 # The main data. 527 write_xy_data(format=format, data=data, file=file, graph_type=graph_type, autoscale=False)
528 529
530 -def delete(align_id=None):
531 """Delete the PCS data corresponding to the alignment ID. 532 533 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted. 534 @type align_id: str or None 535 """ 536 537 # Check the pipe setup. 538 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 539 540 # The IDs. 541 if not align_id: 542 ids = deepcopy(cdp.pcs_ids) 543 else: 544 ids = [align_id] 545 546 # Loop over the alignments, removing all the corresponding data. 547 for id in ids: 548 # The PCS ID. 549 cdp.pcs_ids.pop(cdp.pcs_ids.index(id)) 550 551 # The data type. 552 if hasattr(cdp, 'pcs_data_types') and id in cdp.pcs_data_types: 553 cdp.pcs_data_types.pop(id) 554 555 # The spin data. 556 for spin in spin_loop(): 557 # The data. 558 if hasattr(spin, 'pcs') and id in spin.pcs: 559 spin.pcs.pop(id) 560 561 # The error. 562 if hasattr(spin, 'pcs_err') and id in spin.pcs_err: 563 spin.pcs_err.pop(id) 564 565 # Clean the global data. 566 if not hasattr(cdp, 'rdc_ids') or id not in cdp.rdc_ids: 567 cdp.align_ids.pop(cdp.align_ids.index(id))
568 569
570 -def display(align_id=None, bc=False):
571 """Display the PCS data corresponding to the alignment ID. 572 573 @keyword align_id: The alignment tensor ID string. 574 @type align_id: str 575 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed. 576 @type bc: bool 577 """ 578 579 # Check the pipe setup. 580 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 581 582 # Call the write method with sys.stdout as the file. 583 write(align_id=align_id, file=sys.stdout, bc=bc)
584 585
586 -def opt_uses_pcs(align_id):
587 """Determine if the PCS data for the given alignment ID is needed for optimisation. 588 589 @param align_id: The alignment ID string. 590 @type align_id: str 591 @return: True if the PCS data is to be used for optimisation, False otherwise. 592 @rtype: bool 593 """ 594 595 # No alignment IDs. 596 if not hasattr(cdp, 'pcs_ids'): 597 return False 598 599 # No PCS data for the alignment. 600 if align_id not in cdp.pcs_ids: 601 return False 602 603 # Is the tensor optimised? 604 tensor_flag = opt_uses_tensor(get_tensor_object(align_id)) 605 606 # Is the paramagnetic position optimised? 607 pos_flag = False 608 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 609 pos_flag = True 610 611 # Are the populations optimised? 612 prob_flag = False 613 if cdp.model == 'population': 614 prob_flag = True 615 616 # Not used. 617 if not tensor_flag and not pos_flag and not prob_flag: 618 return False 619 620 # The PCS data is to be used for optimisation. 621 return True
622 623
624 -def q_factors(spin_id=None, sim_index=None, verbosity=1):
625 """Calculate the Q factors for the PCS data. 626 627 @keyword spin_id: The spin ID string used to restrict the Q factor calculation to a subset of all spins. 628 @type spin_id: None or str 629 @keyword sim_index: The optional Monte Carlo simulation index. 630 @type sim_index: None or int 631 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 632 @type verbosity: int 633 """ 634 635 # Check the pipe setup. 636 check_pipe_setup(sequence=True) 637 638 # Simulation flag. 639 sim_flag = (sim_index != None) 640 641 # No PCSs, so no Q factors can be calculated. 642 if not hasattr(cdp, 'pcs_ids') or not len(cdp.pcs_ids): 643 warn(RelaxWarning("No PCS data exists, Q factors cannot be calculated.")) 644 return 645 646 # Q factor dictionary. 647 if not sim_flag: 648 cdp.q_factors_pcs_norm_squared_sum = {} 649 else: 650 if not hasattr(cdp, 'q_factors_pcs_norm_squared_sum_sim'): 651 cdp.q_factors_pcs_norm_squared_sum_sim = {} 652 653 # Loop over the alignments. 654 for align_id in cdp.pcs_ids: 655 # Init. 656 pcs2_sum = 0.0 657 sse = 0.0 658 if sim_flag and align_id not in cdp.q_factors_pcs_norm_squared_sum_sim: 659 cdp.q_factors_pcs_norm_squared_sum_sim[align_id] = [None] * cdp.sim_number 660 661 # Spin loop. 662 spin_count = 0 663 pcs_data = False 664 pcs_bc_data = False 665 for spin in spin_loop(spin_id): 666 pcs_data_align_id = False 667 pcs_bc_data_align_id = False 668 669 # Skip deselected spins. 670 if not spin.select: 671 continue 672 673 # Increment the spin counter. 674 spin_count += 1 675 676 # Data checks. 677 if not sim_flag: 678 if hasattr(spin, 'pcs') and align_id in spin.pcs and spin.pcs[align_id] != None: 679 pcs_data_align_id = True 680 if hasattr(spin, 'pcs_bc') and align_id in spin.pcs_bc and spin.pcs_bc[align_id] != None: 681 pcs_bc_data_align_id = True 682 else: 683 if hasattr(spin, 'pcs_sim') and align_id in spin.pcs_sim and spin.pcs_sim[align_id][sim_index] != None: 684 pcs_data_align_id = True 685 if hasattr(spin, 'pcs_sim_bc') and align_id in spin.pcs_sim_bc and spin.pcs_sim_bc[align_id][sim_index] != None: 686 pcs_bc_data_align_id = True 687 688 # The global flags. 689 pcs_data = pcs_data or pcs_data_align_id 690 pcs_bc_data = pcs_bc_data or pcs_bc_data_align_id 691 692 # Skip spins without PCS data. 693 if not pcs_data_align_id or not pcs_bc_data_align_id: 694 continue 695 696 # Alias the PCS data. 697 if not sim_flag: 698 pcs = spin.pcs[align_id] 699 pcs_bc = spin.pcs_bc[align_id] 700 else: 701 pcs = spin.pcs_sim[align_id][sim_index] 702 pcs_bc = spin.pcs_sim_bc[align_id][sim_index] 703 704 # Sum of squares. 705 sse = sse + (pcs - pcs_bc)**2 706 707 # Sum the PCSs squared (for normalisation). 708 pcs2_sum = pcs2_sum + pcs**2 709 710 # The Q factor for the alignment. 711 if pcs2_sum: 712 Q = sqrt(sse / pcs2_sum) 713 if sim_flag: 714 cdp.q_factors_pcs_norm_squared_sum_sim[align_id][sim_index] = Q 715 else: 716 cdp.q_factors_pcs_norm_squared_sum[align_id] = Q 717 718 # Warnings (and then exit). 719 if not spin_count: 720 if not sim_flag: 721 warn(RelaxWarning("No spins have been used in the calculation, skipping the PCS Q factor calculation.")) 722 return 723 if not pcs_data: 724 if not sim_flag: 725 warn(RelaxWarning("No PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id)) 726 continue 727 if not pcs_bc_data: 728 if not sim_flag: 729 warn(RelaxWarning("No back-calculated PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id)) 730 continue 731 732 # ID and PCS Q factor printout. 733 if verbosity: 734 print("\nPCS Q factors normalised by the sum of PCSs squared:") 735 for align_id in cdp.pcs_ids: 736 if align_id in cdp.q_factors_pcs_norm_squared_sum: 737 if sim_flag: 738 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_pcs_norm_squared_sum_sim[align_id][sim_index])) 739 else: 740 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_pcs_norm_squared_sum[align_id])) 741 742 # The total Q factor. 743 if sim_flag: 744 if not hasattr(cdp, 'q_pcs_norm_squared_sum_sim'): 745 cdp.q_pcs_norm_squared_sum_sim = [None] * cdp.sim_number 746 cdp.q_pcs_norm_squared_sum_sim[sim_index] = 0.0 747 for id in cdp.q_factors_pcs_norm_squared_sum_sim: 748 cdp.q_pcs_norm_squared_sum_sim[sim_index] = cdp.q_pcs_norm_squared_sum_sim[sim_index] + cdp.q_factors_pcs_norm_squared_sum_sim[id][sim_index]**2 749 cdp.q_pcs_norm_squared_sum_sim[sim_index] = cdp.q_pcs_norm_squared_sum_sim[sim_index] / len(cdp.q_factors_pcs_norm_squared_sum_sim) 750 cdp.q_pcs_norm_squared_sum_sim[sim_index] = sqrt(cdp.q_pcs_norm_squared_sum_sim[sim_index]) 751 else: 752 cdp.q_pcs_norm_squared_sum = 0.0 753 for id in cdp.q_factors_pcs_norm_squared_sum: 754 cdp.q_pcs_norm_squared_sum = cdp.q_pcs_norm_squared_sum + cdp.q_factors_pcs_norm_squared_sum[id]**2 755 cdp.q_pcs_norm_squared_sum = cdp.q_pcs_norm_squared_sum / len(cdp.q_factors_pcs_norm_squared_sum) 756 cdp.q_pcs_norm_squared_sum = sqrt(cdp.q_pcs_norm_squared_sum)
757 758
759 -def read(align_id=None, file=None, dir=None, file_data=None, spin_id_col=None, mol_name_col=None, res_num_col=None, res_name_col=None, spin_num_col=None, spin_name_col=None, data_col=None, error_col=None, sep=None, spin_id=None):
760 """Read the PCS data from file. 761 762 @param align_id: The alignment tensor ID string. 763 @type align_id: str 764 @param file: The name of the file to open. 765 @type file: str 766 @param dir: The directory containing the file (defaults to the current directory if None). 767 @type dir: str or None 768 @param file_data: An alternative opening a file, if the data already exists in the correct format. The format is a list of lists where the first index corresponds to the row and the second the column. 769 @type file_data: list of lists 770 @keyword spin_id_col: The column containing the spin ID strings. If supplied, the mol_name_col, res_name_col, res_num_col, spin_name_col, and spin_num_col arguments must be none. 771 @type spin_id_col: int or None 772 @keyword mol_name_col: The column containing the molecule name information. If supplied, spin_id_col must be None. 773 @type mol_name_col: int or None 774 @keyword res_name_col: The column containing the residue name information. If supplied, spin_id_col must be None. 775 @type res_name_col: int or None 776 @keyword res_num_col: The column containing the residue number information. If supplied, spin_id_col must be None. 777 @type res_num_col: int or None 778 @keyword spin_name_col: The column containing the spin name information. If supplied, spin_id_col must be None. 779 @type spin_name_col: int or None 780 @keyword spin_num_col: The column containing the spin number information. If supplied, spin_id_col must be None. 781 @type spin_num_col: int or None 782 @keyword data_col: The column containing the PCS data in Hz. 783 @type data_col: int or None 784 @keyword error_col: The column containing the PCS errors. 785 @type error_col: int or None 786 @keyword sep: The column separator which, if None, defaults to whitespace. 787 @type sep: str or None 788 @keyword spin_id: The spin ID string used to restrict data loading to a subset of all spins. 789 @type spin_id: None or str 790 """ 791 792 # Check the pipe setup. 793 check_pipe_setup(sequence=True) 794 795 # Test if sequence data exists. 796 if not exists_mol_res_spin_data(): 797 raise RelaxNoSequenceError 798 799 # Either the data or error column must be supplied. 800 if data_col == None and error_col == None: 801 raise RelaxError("One of either the data or error column must be supplied.") 802 803 804 # Spin specific data. 805 ##################### 806 807 # Loop over the PCS data. 808 mol_names = [] 809 res_nums = [] 810 res_names = [] 811 spin_nums = [] 812 spin_names = [] 813 values = [] 814 errors = [] 815 for data in read_spin_data(file=file, dir=dir, file_data=file_data, spin_id_col=spin_id_col, mol_name_col=mol_name_col, res_num_col=res_num_col, res_name_col=res_name_col, spin_num_col=spin_num_col, spin_name_col=spin_name_col, data_col=data_col, error_col=error_col, sep=sep, spin_id=spin_id): 816 # Unpack. 817 if data_col and error_col: 818 mol_name, res_num, res_name, spin_num, spin_name, value, error = data 819 elif data_col: 820 mol_name, res_num, res_name, spin_num, spin_name, value = data 821 error = None 822 else: 823 mol_name, res_num, res_name, spin_num, spin_name, error = data 824 value = None 825 826 # Test the error value (cannot be 0.0). 827 if error == 0.0: 828 raise RelaxError("An invalid error value of zero has been encountered.") 829 830 # Get the corresponding spin container. 831 id = generate_spin_id_unique(mol_name=mol_name, res_num=res_num, res_name=res_name, spin_num=spin_num, spin_name=spin_name) 832 spin = return_spin(spin_id=id) 833 if spin == None and spin_id and spin_id[0] == '@': # Allow spin IDs of atom names to be used to specify multi column data. 834 spin = return_spin(spin_id=id+spin_id) 835 if spin == None: 836 warn(RelaxNoSpinWarning(id)) 837 continue 838 839 # Add the data. 840 if data_col: 841 # Initialise. 842 if not hasattr(spin, 'pcs'): 843 spin.pcs = {} 844 845 # Append the value. 846 spin.pcs[align_id] = value 847 848 # Add the error. 849 if error_col: 850 # Initialise. 851 if not hasattr(spin, 'pcs_err'): 852 spin.pcs_err = {} 853 854 # Append the error. 855 spin.pcs_err[align_id] = error 856 857 # Append the data for printout. 858 mol_names.append(mol_name) 859 res_nums.append(res_num) 860 res_names.append(res_name) 861 spin_nums.append(spin_num) 862 spin_names.append(spin_name) 863 values.append(value) 864 errors.append(error) 865 866 # Print out. 867 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, data=values, data_name='PCSs', error=errors, error_name='PCS_error') 868 869 870 # Global (non-spin specific) data. 871 ################################## 872 873 # No data, so return. 874 if not len(values): 875 return 876 877 # Initialise. 878 if not hasattr(cdp, 'align_ids'): 879 cdp.align_ids = [] 880 if not hasattr(cdp, 'pcs_ids'): 881 cdp.pcs_ids = [] 882 883 # Add the PCS id string. 884 if align_id not in cdp.align_ids: 885 cdp.align_ids.append(align_id) 886 if align_id not in cdp.pcs_ids: 887 cdp.pcs_ids.append(align_id)
888 889
890 -def return_pcs_data(sim_index=None, verbosity=0):
891 """Set up the data structures for optimisation using PCSs as base data sets. 892 893 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 894 @type sim_index: None or int 895 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 896 @type verbosity: int 897 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 898 - the PCS values. 899 - the unit vectors connecting the paramagnetic centre (the electron spin) to the spin. 900 - the PCS weight. 901 - the experimental temperatures. 902 - the spectrometer frequencies. 903 - pseudo_flags, the list of flags indicating if the interatomic data contains a pseudo-atom (as 1's and 0's). 904 @rtype: tuple of (numpy rank-2 float64 array, numpy rank-2 float64 array, numpy rank-2 float64 array, list of float, list of float, numpy rank-1 int32 array) 905 """ 906 907 # Initial printout. 908 if verbosity: 909 print("\nPCS data counts:") 910 911 # Data setup tests. 912 if not hasattr(cdp, 'paramagnetic_centre') and (hasattr(cdp, 'paramag_centre_fixed') and cdp.paramag_centre_fixed): 913 raise RelaxError("The paramagnetic centre has not yet been specified.") 914 if not hasattr(cdp, 'temperature'): 915 raise RelaxError("The experimental temperatures have not been set.") 916 if not hasattr(cdp, 'spectrometer_frq'): 917 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 918 919 # Sort out pseudo-atoms first. This only needs to be called once. 920 setup_pseudoatom_pcs() 921 922 # Initialise. 923 pcs = [] 924 pcs_err = [] 925 pcs_weight = [] 926 temp = [] 927 frq = [] 928 pseudo_flags = [] 929 930 # The PCS data. 931 for i in range(len(cdp.align_ids)): 932 # Alias the ID. 933 align_id = cdp.align_ids[i] 934 935 # Skip non-optimised data. 936 if not opt_uses_align_data(align_id): 937 continue 938 939 # Append empty arrays to the PCS structures. 940 pcs.append([]) 941 pcs_err.append([]) 942 pcs_weight.append([]) 943 944 # Get the temperature for the PCS constant. 945 if align_id in cdp.temperature: 946 temp.append(cdp.temperature[align_id]) 947 948 # The temperature must be given! 949 else: 950 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id) 951 952 # Get the spectrometer frequency in Tesla units for the PCS constant. 953 if align_id in cdp.spectrometer_frq: 954 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H')) 955 956 # The frequency must be given! 957 else: 958 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id) 959 960 # Spin loop. 961 j = 0 962 for spin in spin_loop(): 963 # Skip deselected spins. 964 if not spin.select: 965 continue 966 967 # Skip spins without PCS data. 968 if not hasattr(spin, 'pcs'): 969 continue 970 971 # Append the PCSs to the list. 972 if align_id in spin.pcs: 973 if sim_index != None: 974 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 975 else: 976 pcs[-1].append(spin.pcs[align_id]) 977 978 # Increment the PCS count. 979 if pcs[-1][-1] != None: 980 j += 1 981 982 # No data. 983 else: 984 pcs[-1].append(None) 985 986 # Append the PCS errors. 987 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 988 pcs_err[-1].append(spin.pcs_err[align_id]) 989 else: 990 pcs_err[-1].append(None) 991 992 # Append the weight. 993 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight: 994 pcs_weight[-1].append(spin.pcs_weight[align_id]) 995 else: 996 pcs_weight[-1].append(1.0) 997 998 # ID and PCS count printout. 999 if verbosity: 1000 print(" Alignment ID '%s': %i" % (align_id, j)) 1001 1002 # Pseudo-atom. 1003 for spin in spin_loop(): 1004 if is_pseudoatom(spin): 1005 pseudo_flags.append(1) 1006 else: 1007 pseudo_flags.append(0) 1008 1009 # Convert to numpy objects. 1010 pcs = array(pcs, float64) 1011 pcs_err = array(pcs_err, float64) 1012 pcs_weight = array(pcs_weight, float64) 1013 pseudo_flags = array(pseudo_flags, int32) 1014 1015 # Convert the PCS from ppm to no units. 1016 pcs = pcs * 1e-6 1017 pcs_err = pcs_err * 1e-6 1018 1019 # Return the data structures. 1020 return pcs, pcs_err, pcs_weight, temp, frq, pseudo_flags
1021 1022
1023 -def set_errors(align_id=None, spin_id=None, sd=None):
1024 """Set the PCS errors if not already present. 1025 1026 @keyword align_id: The optional alignment tensor ID string. 1027 @type align_id: str 1028 @keyword spin_id: The optional spin ID string. 1029 @type spin_id: None or str 1030 @keyword sd: The PCS standard deviation in ppm. 1031 @type sd: float or int. 1032 """ 1033 1034 # Check the pipe setup. 1035 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 1036 1037 # Convert the align IDs to an array, or take all IDs. 1038 if align_id: 1039 align_ids = [align_id] 1040 else: 1041 align_ids = cdp.pcs_ids 1042 1043 # Loop over the spins. 1044 for spin in spin_loop(spin_id): 1045 # Skip deselected spins. 1046 if not spin.select: 1047 continue 1048 1049 # Skip spins with no PCSs. 1050 if not hasattr(spin, 'pcs') or (align_id and not align_id in spin.pcs): 1051 continue 1052 1053 # No data structure. 1054 if not hasattr(spin, 'pcs_err'): 1055 spin.pcs_err = {} 1056 1057 # Set the error. 1058 for id in align_ids: 1059 spin.pcs_err[id] = sd
1060 1061
1062 -def setup_pseudoatom_pcs():
1063 """Make sure that the spin systems are properly set up for pseudo-atoms and PCSs. 1064 1065 All spin data containers which are a member of a pseudo-atom will be deselected. 1066 """ 1067 1068 # Loop over all spin data containers. 1069 for pseudospin, pseudospin_id in spin_loop(return_id=True): 1070 # No pseudo-atom, so do nothing. 1071 if not is_pseudoatom(pseudospin): 1072 return 1073 1074 # Loop over the atoms of the pseudo-atom. 1075 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 1076 # Deselect if needed. 1077 if spin.select: 1078 warn(RelaxWarning("Deselecting the '%s' spin as it is a member of the '%s' pseudo-atom system." % (spin_id, pseudospin_id))) 1079 spin.select = False
1080 1081
1082 -def structural_noise(align_id=None, rmsd=0.2, sim_num=1000, file=None, dir=None, force=False):
1083 """Determine the PCS error due to structural noise via simulation. 1084 1085 For the simulation the following must already be set up in the current data pipe: 1086 1087 - The position of the paramagnetic centre. 1088 - The alignment and magnetic susceptibility tensor. 1089 1090 The protocol for the simulation is as follows: 1091 1092 - The lanthanide or paramagnetic centre position will be fixed. Its motion is assumed to be on the femto- to pico- and nanosecond timescales. Hence the motion is averaged over the evolution of the PCS and can be ignored. 1093 - The positions of the nuclear spins will be randomised N times using a multivariate normal distribution. 1094 - The PCS for the randomised position will be back calculated. 1095 - The PCS standard deviation will be calculated from the N randomised PCS values. 1096 1097 The standard deviation will both be stored in the spin container data structure in the relax data store as well as being added to the already present PCS error (using variance addition). This will then be used in any optimisations involving the PCS. 1098 1099 If the alignment ID string is not supplied, the procedure will be applied to the PCS data from all alignments. 1100 1101 1102 @keyword align_id: The alignment tensor ID string. 1103 @type align_id: str 1104 @keyword rmsd: The atomic position RMSD, in Angstrom, to randomise the spin positions with for the simulations. 1105 @type rmsd: float 1106 @keyword sim_num: The number of simulations, N, to perform to determine the structural noise component of the PCS errors. 1107 @type sim_num: int 1108 @keyword file: The optional name of the Grace file to plot the structural errors verses the paramagnetic centre to spin distances. 1109 @type file: None or str 1110 @keyword dir: The directory name to place the Grace file into. 1111 @type dir: None or str 1112 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1113 @type force: bool 1114 """ 1115 1116 # Check the pipe setup. 1117 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True, paramag_centre=True) 1118 1119 # Convert the align IDs to an array, or take all IDs. 1120 if align_id: 1121 align_ids = [align_id] 1122 else: 1123 align_ids = cdp.align_ids 1124 1125 # Initialise some numpy data structures for use in the simulations. 1126 grace_data = [] 1127 unit_vect = zeros(3, float64) 1128 pcs = {} 1129 for id in align_ids: 1130 grace_data.append([]) 1131 pcs[id] = zeros(sim_num, float64) 1132 1133 # Pre-calculate the covariance matrix as the spherical covariance. 1134 cov = rmsd**2 * eye(3) 1135 1136 # Print out. 1137 print("Executing %i simulations for each spin system." % sim_num) 1138 1139 # Loop over the spins. 1140 for spin, spin_id in spin_loop(return_id=True): 1141 # Deselected spins. 1142 if not spin.select: 1143 continue 1144 1145 # Skip spins with no PCS or position. 1146 if not hasattr(spin, 'pcs'): 1147 continue 1148 if not hasattr(spin, 'pos'): 1149 continue 1150 1151 # Print out. 1152 print(spin_id) 1153 1154 # Average the atom position. 1155 if type(spin.pos[0]) in [float, float64]: 1156 pos = spin.pos 1157 else: 1158 pos = zeros(3, float64) 1159 for i in range(len(spin.pos)): 1160 pos += spin.pos[i] 1161 pos = pos / len(spin.pos) 1162 1163 # The original vector length (for the Grace plot). 1164 orig_vect = pos - cdp.paramagnetic_centre 1165 orig_r = norm(orig_vect) 1166 1167 # Sample from the multivariate normal distribution. 1168 new_pos = multivariate_normal(pos, cov, sim_num) 1169 1170 # Loop over the N randomisations. 1171 for i in range(sim_num): 1172 # The vector and length. 1173 vect = new_pos[i] - cdp.paramagnetic_centre 1174 r = norm(vect) 1175 vect = vect / r 1176 1177 # Loop over the alignments. 1178 for id in align_ids: 1179 # No PCS value, so skip. 1180 if id not in spin.pcs: 1181 continue 1182 1183 # Calculate the PCS constant. 1184 dj = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / periodic_table.gyromagnetic_ratio('1H'), r/1e10) 1185 1186 # Calculate the PCS value (in ppm). 1187 pcs[id][i] = pcs_tensor(dj, vect, cdp.align_tensors[get_tensor_index(id)].A) * 1e6 1188 1189 # Initialise if necessary. 1190 if not hasattr(spin, 'pcs_struct_err'): 1191 spin.pcs_struct_err = {} 1192 1193 # Loop over the alignments. 1194 align_index = 0 1195 for id in align_ids: 1196 # No PCS value, so skip. 1197 if id not in spin.pcs or spin.pcs[id] == None: 1198 align_index += 1 1199 continue 1200 1201 # The PCS standard deviation. 1202 sd = std(pcs[id]) 1203 1204 # Remove the previous error. 1205 if id in spin.pcs_struct_err: 1206 warn(RelaxWarning("Removing the previous structural error value from the PCS error of the spin '%s' for the alignment ID '%s'." % (spin_id, id))) 1207 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 - spin.pcs_struct_err[id]**2) 1208 1209 # Store the structural error. 1210 spin.pcs_struct_err[id] = sd 1211 1212 # Add it to the PCS error (with variance addition). 1213 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 + sd**2) 1214 1215 # Store the data for the Grace plot. 1216 grace_data[align_index].append([orig_r, sd, spin_id]) 1217 1218 # Increment the alignment index. 1219 align_index += 1 1220 1221 # The Grace output. 1222 if file: 1223 # Open the Grace file for writing. 1224 file = open_write_file(file, dir, force) 1225 1226 # The header. 1227 write_xy_header(format='grace', file=file, title="PCS structural noise", subtitle="%s Angstrom structural noise"%rmsd, data_type=['pcs_bc', 'pcs'], sets=[len(align_ids)], set_names=[align_ids], symbol_sizes=[[0.5]*len(align_ids)], linetype=[[0]*len(align_ids)], axis_labels=[["Ln\\S3+\\N to spin distance (Angstrom)", "PCS standard deviation (ppm)"]]) 1228 1229 # The main data. 1230 write_xy_data(format='grace', data=[grace_data], file=file, graph_type='xy')
1231 1232
1233 -def weight(align_id=None, spin_id=None, weight=1.0):
1234 """Set optimisation weights on the PCS data. 1235 1236 @keyword align_id: The alignment tensor ID string. 1237 @type align_id: str 1238 @keyword spin_id: The spin ID string. 1239 @type spin_id: None or str 1240 @keyword weight: The optimisation weight. The higher the value, the more importance the PCS will have. 1241 @type weight: float or int. 1242 """ 1243 1244 # Check the pipe setup. 1245 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 1246 1247 # Loop over the spins. 1248 for spin in spin_loop(spin_id): 1249 # No data structure. 1250 if not hasattr(spin, 'pcs_weight'): 1251 spin.pcs_weight = {} 1252 1253 # Set the weight. 1254 spin.pcs_weight[align_id] = weight
1255 1256
1257 -def write(align_id=None, file=None, dir=None, bc=False, force=False):
1258 """Display the PCS data corresponding to the alignment ID. 1259 1260 @keyword align_id: The alignment tensor ID string. 1261 @type align_id: str 1262 @keyword file: The file name or object to write to. 1263 @type file: str or file object 1264 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 1265 @type dir: str 1266 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 1267 @type bc: bool 1268 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1269 @type force: bool 1270 """ 1271 1272 # Check the pipe setup. 1273 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 1274 1275 # Open the file for writing. 1276 file = open_write_file(file, dir, force) 1277 1278 # Loop over the spins and collect the data. 1279 mol_names = [] 1280 res_nums = [] 1281 res_names = [] 1282 spin_nums = [] 1283 spin_names = [] 1284 values = [] 1285 errors = [] 1286 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 1287 # Skip deselected spins. 1288 if not spin.select: 1289 continue 1290 1291 # Skip spins with no PCSs. 1292 if not bc and (not hasattr(spin, 'pcs') or not align_id in spin.pcs): 1293 continue 1294 elif bc and (not hasattr(spin, 'pcs_bc') or align_id not in spin.pcs_bc): 1295 continue 1296 1297 # Append the spin data. 1298 mol_names.append(mol_name) 1299 res_nums.append(res_num) 1300 res_names.append(res_name) 1301 spin_nums.append(spin.num) 1302 spin_names.append(spin.name) 1303 1304 # The value. 1305 if bc: 1306 values.append(spin.pcs_bc[align_id]) 1307 else: 1308 values.append(spin.pcs[align_id]) 1309 1310 # The error. 1311 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err: 1312 errors.append(spin.pcs_err[align_id]) 1313 else: 1314 errors.append(None) 1315 1316 # Write out. 1317 write_spin_data(file=file, mol_names=mol_names, res_nums=res_nums, res_names=res_names, spin_nums=spin_nums, spin_names=spin_names, data=values, data_name='PCSs', error=errors, error_name='PCS_error')
1318