pipe_control.pcs

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for the manipulation of pseudo-contact shift data.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi, sqrt 28 from numpy import array, float64, int32, ones, std, zeros 29 from numpy.linalg import norm 30 from random import gauss 31 import sys 32 from warnings import warn 33 34 # relax module imports. 35 from lib.alignment.pcs import ave_pcs_tensor, pcs_tensor 36 from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoPdbError, RelaxNoPCSError, RelaxNoSequenceError 37 from lib.geometry.vectors import random_unit_vector 38 from lib.io import open_write_file, read_spin_data, write_spin_data 39 from lib.physical_constants import g1H, pcs_constant 40 from lib.warnings import RelaxWarning, RelaxNoSpinWarning 41 from pipe_control import grace, pipes 42 from pipe_control.align_tensor import get_tensor_index, get_tensor_object, opt_uses_align_data, opt_uses_tensor 43 from pipe_control.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, is_pseudoatom, return_spin, spin_index_loop, spin_loop 44 45

46 -def back_calc(align_id=None):

47 """Back calculate the PCS from the alignment tensor. 48 49 @keyword align_id: The alignment tensor ID string. 50 @type align_id: str 51 """ 52 53 # Check the pipe setup. 54 check_pipe_setup(pcs_id=align_id, sequence=True, N=True, tensors=True, paramag_centre=True) 55 56 # Convert the align IDs to an array, or take all IDs. 57 if align_id: 58 align_ids = [align_id] 59 else: 60 align_ids = cdp.align_ids 61 62 # Add the ID to the PCS IDs, if needed. 63 for align_id in align_ids: 64 # Init. 65 if not hasattr(cdp, 'pcs_ids'): 66 cdp.pcs_ids = [] 67 68 # Add the ID. 69 if align_id not in cdp.pcs_ids: 70 cdp.pcs_ids.append(align_id) 71 72 # The weights. 73 weights = ones(cdp.N, float64) / cdp.N 74 75 # Unit vector data structure init. 76 unit_vect = zeros((cdp.N, 3), float64) 77 78 # Loop over the spins. 79 count = 0 80 for spin in spin_loop(): 81 # Skip spins with no position. 82 if not hasattr(spin, 'pos'): 83 continue 84 85 # Atom positions. 86 pos = spin.pos 87 if type(pos[0]) in [float, float64]: 88 pos = [pos] * cdp.N 89 90 # Loop over the alignments. 91 for id in align_ids: 92 # Vectors. 93 vect = zeros((cdp.N, 3), float64) 94 r = zeros(cdp.N, float64) 95 dj = zeros(cdp.N, float64) 96 for c in range(cdp.N): 97 # The vector. 98 vect[c] = pos[c] - cdp.paramagnetic_centre 99 100 # The length. 101 r[c] = norm(vect[c]) 102 103 # Normalise (only if the vector has length). 104 if r[c]: 105 vect[c] = vect[c] / r[c] 106 107 # Calculate the PCS constant. 108 dj[c] = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / g1H, r[c]/1e10) 109 110 # Initialise if necessary. 111 if not hasattr(spin, 'pcs_bc'): 112 spin.pcs_bc = {} 113 114 # Calculate the PCSs (in ppm). 115 spin.pcs_bc[id] = ave_pcs_tensor(dj, vect, cdp.N, cdp.align_tensors[get_tensor_index(align_id=id)].A, weights=weights) * 1e6 116 117 # Increment the counter. 118 count += 1 119 120 # No PCSs calculated. 121 if not count: 122 warn(RelaxWarning("No PCSs have been back calculated, probably due to missing spin position information."))

123 124

125 -def centre(pos=None, atom_id=None, pipe=None, verbosity=1, ave_pos=False, force=False):

126 """Specify the atom in the loaded structure corresponding to the paramagnetic centre. 127 128 @keyword pos: The atomic position. If set, the atom_id string will be ignored. 129 @type pos: list of float 130 @keyword atom_id: The atom identification string. 131 @type atom_id: str 132 @keyword pipe: An alternative data pipe to extract the paramagnetic centre from. 133 @type pipe: None or str 134 @keyword verbosity: The amount of information to print out. The bigger the number, the more information. 135 @type verbosity: int 136 @keyword ave_pos: A flag which if True causes the atomic positions from multiple models to be averaged. 137 @type ave_pos: bool 138 @keyword force: A flag which if True will cause the current PCS centre to be overwritten. 139 """ 140 141 # The data pipe. 142 if pipe == None: 143 pipe = pipes.cdp_name() 144 145 # Check the pipe setup. 146 check_pipe_setup(pipe=pipe) 147 148 # Get the data pipes. 149 source_dp = pipes.get_pipe(pipe) 150 151 # Test if the structure has been loaded. 152 if not hasattr(source_dp, 'structure'): 153 raise RelaxNoPdbError 154 155 # Test the centre has already been set. 156 if not force and hasattr(cdp, 'paramagnetic_centre'): 157 raise RelaxError("The paramagnetic centre has already been set to the coordinates " + repr(cdp.paramagnetic_centre) + ".") 158 159 # Position is supplied. 160 if pos != None: 161 centre = array(pos) 162 num_pos = 1 163 full_pos_list = [] 164 165 # Position from a loaded structure. 166 else: 167 # Get the positions. 168 centre = zeros(3, float64) 169 full_pos_list = [] 170 num_pos = 0 171 for spin, spin_id in spin_loop(atom_id, pipe=pipe, return_id=True): 172 # No atomic positions. 173 if not hasattr(spin, 'pos'): 174 continue 175 176 # Spin position list. 177 if isinstance(spin.pos[0], float) or isinstance(spin.pos[0], float64): 178 pos_list = [spin.pos] 179 else: 180 pos_list = spin.pos 181 182 # Loop over the model positions. 183 for pos in pos_list: 184 full_pos_list.append(pos) 185 centre = centre + array(pos) 186 num_pos = num_pos + 1 187 188 # No positional information! 189 if not num_pos: 190 raise RelaxError("No positional information could be found for the spin '%s'." % atom_id) 191 192 # Averaging. 193 centre = centre / float(num_pos) 194 195 # Print out. 196 if verbosity: 197 print("Paramagnetic centres located at:") 198 for pos in full_pos_list: 199 print(" [%8.3f, %8.3f, %8.3f]" % (pos[0], pos[1], pos[2])) 200 print("\nAverage paramagnetic centre located at:") 201 print(" [%8.3f, %8.3f, %8.3f]" % (centre[0], centre[1], centre[2])) 202 203 # Set the centre (place it into the current data pipe). 204 if ave_pos: 205 if verbosity: 206 print("\nUsing the average paramagnetic position.") 207 cdp.paramagnetic_centre = centre 208 else: 209 if verbosity: 210 print("\nUsing all paramagnetic positions.") 211 cdp.paramagnetic_centre = full_pos_list

212 213

214 -def check_pipe_setup(pipe=None, pcs_id=None, sequence=False, N=False, tensors=False, pcs=False, paramag_centre=False):

215 """Check that the current data pipe has been setup sufficiently. 216 217 @keyword pipe: The data pipe to check, defaulting to the current pipe. 218 @type pipe: None or str 219 @keyword pcs_id: The PCS ID string to check for in cdp.pcs_ids. 220 @type pcs_id: None or str 221 @keyword sequence: A flag which when True will invoke the sequence data check. 222 @type sequence: bool 223 @keyword N: A flag which if True will check that cdp.N is set. 224 @type N: bool 225 @keyword tensors: A flag which if True will check that alignment tensors exist. 226 @type tensors: bool 227 @keyword pcs: A flag which if True will check that PCSs exist. 228 @type pcs: bool 229 @keyword paramag_centre: A flag which if True will check that the paramagnetic centre has been set. 230 @type paramag_centre: bool 231 """ 232 233 # The data pipe. 234 if pipe == None: 235 pipe = pipes.cdp_name() 236 237 # Get the data pipe. 238 dp = pipes.get_pipe(pipe) 239 240 # Test if the current data pipe exists. 241 pipes.test(pipe) 242 243 # Test if sequence data exists. 244 if sequence and not exists_mol_res_spin_data(pipe): 245 raise RelaxNoSequenceError(pipe) 246 247 # Check for dp.N. 248 if N and not hasattr(dp, 'N'): 249 raise RelaxError("The number of states N has not been set.") 250 251 # Check that alignment tensors are present. 252 if tensors and (not hasattr(dp, 'align_tensors') or len(dp.align_tensors) == 0): 253 raise RelaxNoTensorError('alignment') 254 255 # Test for the alignment ID. 256 if pcs_id and (not hasattr(dp, 'align_ids') or pcs_id not in dp.align_ids): 257 raise RelaxNoAlignError(pcs_id, pipe) 258 259 # Test if PCS data exists. 260 if pcs and not hasattr(dp, 'align_ids'): 261 raise RelaxNoAlignError() 262 if pcs and not hasattr(dp, 'pcs_ids'): 263 raise RelaxNoPCSError() 264 elif pcs and pcs_id and pcs_id not in dp.pcs_ids: 265 raise RelaxNoPCSError(pcs_id) 266 267 # Test if the paramagnetic centre is set. 268 if paramag_centre and not hasattr(cdp, 'paramagnetic_centre'): 269 raise RelaxError("The paramagnetic centre has not been defined.")

270 271

272 -def copy(pipe_from=None, pipe_to=None, align_id=None):

273 """Copy the PCS data from one data pipe to another. 274 275 @keyword pipe_from: The data pipe to copy the PCS data from. This defaults to the current data pipe. 276 @type pipe_from: str 277 @keyword pipe_to: The data pipe to copy the PCS data to. This defaults to the current data pipe. 278 @type pipe_to: str 279 @param align_id: The alignment ID string. 280 @type align_id: str 281 """ 282 283 # Defaults. 284 if pipe_from == None and pipe_to == None: 285 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 286 elif pipe_from == None: 287 pipe_from = pipes.cdp_name() 288 elif pipe_to == None: 289 pipe_to = pipes.cdp_name() 290 291 # Check the pipe setup. 292 check_pipe_setup(pipe=pipe_from, pcs_id=align_id, sequence=True, pcs=True) 293 check_pipe_setup(pipe=pipe_to, sequence=True) 294 295 # Get the data pipes. 296 dp_from = pipes.get_pipe(pipe_from) 297 dp_to = pipes.get_pipe(pipe_to) 298 299 # The IDs. 300 if align_id == None: 301 align_ids = dp_from.align_ids 302 else: 303 align_ids = [align_id] 304 305 # Init target pipe global structures. 306 if not hasattr(dp_to, 'align_ids'): 307 dp_to.align_ids = [] 308 if not hasattr(dp_to, 'pcs_ids'): 309 dp_to.pcs_ids = [] 310 311 # Loop over the align IDs. 312 for align_id in align_ids: 313 # Copy the global data. 314 if align_id not in dp_to.align_ids and align_id not in dp_to.align_ids: 315 dp_to.align_ids.append(align_id) 316 if align_id in dp_from.pcs_ids and align_id not in dp_to.pcs_ids: 317 dp_to.pcs_ids.append(align_id) 318 319 # Spin loop. 320 for mol_index, res_index, spin_index in spin_index_loop(): 321 # Alias the spin containers. 322 spin_from = dp_from.mol[mol_index].res[res_index].spin[spin_index] 323 spin_to = dp_to.mol[mol_index].res[res_index].spin[spin_index] 324 325 # No data or errors. 326 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): 327 continue 328 329 # Initialise the spin data if necessary. 330 if hasattr(spin_from, 'pcs') and not hasattr(spin_to, 'pcs'): 331 spin_to.pcs = {} 332 if hasattr(spin_from, 'pcs_err') and not hasattr(spin_to, 'pcs_err'): 333 spin_to.pcs_err = {} 334 335 # Copy the value and error from pipe_from. 336 if hasattr(spin_from, 'pcs'): 337 spin_to.pcs[align_id] = spin_from.pcs[align_id] 338 if hasattr(spin_from, 'pcs_err'): 339 spin_to.pcs_err[align_id] = spin_from.pcs_err[align_id]

340 341

342 -def corr_plot(format=None, file=None, dir=None, force=False):

343 """Generate a correlation plot of the measured vs. back-calculated PCSs. 344 345 @keyword format: The format for the plot file. The following values are accepted: 'grace', a Grace plot; None, a plain text file. 346 @type format: str or None 347 @keyword file: The file name or object to write to. 348 @type file: str or file object 349 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 350 @type dir: str 351 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 352 @type force: bool 353 """ 354 355 # Check the pipe setup. 356 check_pipe_setup(sequence=True) 357 358 # Does PCS data exist? 359 if not hasattr(cdp, 'pcs_ids') or not cdp.pcs_ids: 360 warn(RelaxWarning("No PCS data exists, skipping file creation.")) 361 return 362 363 # Open the file for writing. 364 file = open_write_file(file, dir, force) 365 366 # Init. 367 data = [] 368 369 # The diagonal. 370 data.append([[-100, -100, 0], [100, 100, 0]]) 371 372 # The spin types. 373 types = [] 374 for spin in spin_loop(): 375 if spin.element not in types: 376 types.append(spin.element) 377 378 # Loop over the PCS data. 379 for align_id in cdp.pcs_ids: 380 # Loop over the spin types. 381 for i in range(len(types)): 382 # Append a new list for this alignment. 383 data.append([]) 384 385 # Errors present? 386 err_flag = False 387 for spin in spin_loop(): 388 # Skip deselected spins. 389 if not spin.select: 390 continue 391 392 # Error present. 393 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 394 err_flag = True 395 break 396 397 # Loop over the spins. 398 for spin, spin_id in spin_loop(return_id=True): 399 # Skip deselected spins. 400 if not spin.select: 401 continue 402 403 # Incorrect spin type. 404 if spin.element != types[i]: 405 continue 406 407 # Skip if data is missing. 408 if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc') or not align_id in spin.pcs.keys() or not align_id in spin.pcs_bc.keys(): 409 continue 410 411 # Append the data. 412 data[-1].append([spin.pcs_bc[align_id], spin.pcs[align_id]]) 413 414 # Errors. 415 if err_flag: 416 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 417 data[-1][-1].append(spin.pcs_err[align_id]) 418 else: 419 data[-1][-1].append(None) 420 421 # Label. 422 data[-1][-1].append(spin_id) 423 424 # The data size. 425 size = len(data) 426 427 # Only one data set. 428 data = [data] 429 430 # Graph type. 431 if err_flag: 432 graph_type = 'xydy' 433 else: 434 graph_type = 'xy' 435 436 # Grace file. 437 if format == 'grace': 438 # The set names. 439 set_names = [None] 440 for i in range(len(cdp.pcs_ids)): 441 for j in range(len(types)): 442 set_names.append("%s (%s)" % (cdp.pcs_ids[i], types[j])) 443 444 # The header. 445 grace.write_xy_header(file=file, title="PCS correlation plot", data_type=['pcs_bc', 'pcs'], sets=[size], set_names=[set_names], linestyle=[[2]+[0]*size], legend_pos=[[1, 0.5]]) 446 447 # The main data. 448 grace.write_xy_data(data=data, file=file, graph_type=graph_type)

449 450

451 -def delete(align_id=None):

452 """Delete the PCS data corresponding to the alignment ID. 453 454 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted. 455 @type align_id: str or None 456 """ 457 458 # Check the pipe setup. 459 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 460 461 # The IDs. 462 if not align_id: 463 ids = deepcopy(cdp.pcs_ids) 464 else: 465 ids = [align_id] 466 467 # Loop over the alignments, removing all the corresponding data. 468 for id in ids: 469 # The PCS ID. 470 cdp.pcs_ids.pop(cdp.pcs_ids.index(id)) 471 472 # The data type. 473 if hasattr(cdp, 'pcs_data_types') and id in cdp.pcs_data_types: 474 cdp.pcs_data_types.pop(id) 475 476 # The spin data. 477 for spin in spin_loop(): 478 # The data. 479 if hasattr(spin, 'pcs') and id in spin.pcs: 480 spin.pcs.pop(id) 481 482 # The error. 483 if hasattr(spin, 'pcs_err') and id in spin.pcs_err: 484 spin.pcs_err.pop(id) 485 486 # Clean the global data. 487 if not hasattr(cdp, 'rdc_ids') or id not in cdp.rdc_ids: 488 cdp.align_ids.pop(cdp.align_ids.index(id))

489 490

491 -def display(align_id=None, bc=False):

492 """Display the PCS data corresponding to the alignment ID. 493 494 @keyword align_id: The alignment tensor ID string. 495 @type align_id: str 496 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed. 497 @type bc: bool 498 """ 499 500 # Check the pipe setup. 501 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 502 503 # Call the write method with sys.stdout as the file. 504 write(align_id=align_id, file=sys.stdout, bc=bc)

505 506

507 -def opt_uses_pcs(align_id):

508 """Determine if the PCS data for the given alignment ID is needed for optimisation. 509 510 @param align_id: The alignment ID string. 511 @type align_id: str 512 @return: True if the PCS data is to be used for optimisation, False otherwise. 513 @rtype: bool 514 """ 515 516 # No alignment IDs. 517 if not hasattr(cdp, 'pcs_ids'): 518 return False 519 520 # No PCS data for the alignment. 521 if align_id not in cdp.pcs_ids: 522 return False 523 524 # Is the tensor optimised? 525 tensor_flag = opt_uses_tensor(get_tensor_object(align_id)) 526 527 # Is the paramagnetic position optimised? 528 pos_flag = False 529 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 530 pos_flag = True 531 532 # Are the populations optimised? 533 prob_flag = False 534 if cdp.model == 'population': 535 prob_flag = True 536 537 # Not used. 538 if not tensor_flag and not pos_flag and not prob_flag: 539 return False 540 541 # The PCS data is to be used for optimisation. 542 return True

543 544

545 -def q_factors(spin_id=None):

546 """Calculate the Q-factors for the PCS data. 547 548 @keyword spin_id: The spin ID string used to restrict the Q-factor calculation to a subset of all spins. 549 @type spin_id: None or str 550 """ 551 552 # Check the pipe setup. 553 check_pipe_setup(sequence=True) 554 555 # No PCSs, so no Q factors can be calculated. 556 if not hasattr(cdp, 'pcs_ids') or not len(cdp.pcs_ids): 557 warn(RelaxWarning("No PCS data exists, Q factors cannot be calculated.")) 558 return 559 560 # Q-factor dictionary. 561 cdp.q_factors_pcs = {} 562 563 # Loop over the alignments. 564 for align_id in cdp.pcs_ids: 565 # Init. 566 pcs2_sum = 0.0 567 sse = 0.0 568 569 # Spin loop. 570 spin_count = 0 571 pcs_data = False 572 pcs_bc_data = False 573 for spin in spin_loop(spin_id): 574 # Skip deselected spins. 575 if not spin.select: 576 continue 577 578 # Increment the spin counter. 579 spin_count += 1 580 581 # Data checks. 582 if hasattr(spin, 'pcs') and align_id in spin.pcs: 583 pcs_data = True 584 if hasattr(spin, 'pcs_bc') and align_id in spin.pcs_bc: 585 pcs_bc_data = True 586 587 # Skip spins without PCS data. 588 if not hasattr(spin, 'pcs') or not hasattr(spin, 'pcs_bc') or not align_id in spin.pcs or spin.pcs[align_id] == None or not align_id in spin.pcs_bc or spin.pcs_bc[align_id] == None: 589 continue 590 591 # Sum of squares. 592 sse = sse + (spin.pcs[align_id] - spin.pcs_bc[align_id])**2 593 594 # Sum the PCSs squared (for normalisation). 595 pcs2_sum = pcs2_sum + spin.pcs[align_id]**2 596 597 # The Q-factor for the alignment. 598 if pcs2_sum: 599 Q = sqrt(sse / pcs2_sum) 600 cdp.q_factors_pcs[align_id] = Q 601 602 # Warnings (and then exit). 603 if not spin_count: 604 warn(RelaxWarning("No spins have been used in the calculation, skipping the PCS Q factor calculation.")) 605 return 606 if not pcs_data: 607 warn(RelaxWarning("No PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id)) 608 continue 609 if not pcs_bc_data: 610 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)) 611 continue 612 613 # The total Q-factor. 614 cdp.q_pcs = 0.0 615 for id in cdp.q_factors_pcs: 616 cdp.q_pcs = cdp.q_pcs + cdp.q_factors_pcs[id]**2 617 cdp.q_pcs = cdp.q_pcs / len(cdp.q_factors_pcs) 618 cdp.q_pcs = sqrt(cdp.q_pcs)

619 620

621 -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):

622 """Read the PCS data from file. 623 624 @param align_id: The alignment tensor ID string. 625 @type align_id: str 626 @param file: The name of the file to open. 627 @type file: str 628 @param dir: The directory containing the file (defaults to the current directory if None). 629 @type dir: str or None 630 @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. 631 @type file_data: list of lists 632 @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. 633 @type spin_id_col: int or None 634 @keyword mol_name_col: The column containing the molecule name information. If supplied, spin_id_col must be None. 635 @type mol_name_col: int or None 636 @keyword res_name_col: The column containing the residue name information. If supplied, spin_id_col must be None. 637 @type res_name_col: int or None 638 @keyword res_num_col: The column containing the residue number information. If supplied, spin_id_col must be None. 639 @type res_num_col: int or None 640 @keyword spin_name_col: The column containing the spin name information. If supplied, spin_id_col must be None. 641 @type spin_name_col: int or None 642 @keyword spin_num_col: The column containing the spin number information. If supplied, spin_id_col must be None. 643 @type spin_num_col: int or None 644 @keyword data_col: The column containing the PCS data in Hz. 645 @type data_col: int or None 646 @keyword error_col: The column containing the PCS errors. 647 @type error_col: int or None 648 @keyword sep: The column separator which, if None, defaults to whitespace. 649 @type sep: str or None 650 @keyword spin_id: The spin ID string used to restrict data loading to a subset of all spins. 651 @type spin_id: None or str 652 """ 653 654 # Check the pipe setup. 655 check_pipe_setup(sequence=True) 656 657 # Test if sequence data exists. 658 if not exists_mol_res_spin_data(): 659 raise RelaxNoSequenceError 660 661 # Either the data or error column must be supplied. 662 if data_col == None and error_col == None: 663 raise RelaxError("One of either the data or error column must be supplied.") 664 665 666 # Spin specific data. 667 ##################### 668 669 # Loop over the PCS data. 670 mol_names = [] 671 res_nums = [] 672 res_names = [] 673 spin_nums = [] 674 spin_names = [] 675 values = [] 676 errors = [] 677 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): 678 # Unpack. 679 if data_col and error_col: 680 mol_name, res_num, res_name, spin_num, spin_name, value, error = data 681 elif data_col: 682 mol_name, res_num, res_name, spin_num, spin_name, value = data 683 error = None 684 else: 685 mol_name, res_num, res_name, spin_num, spin_name, error = data 686 value = None 687 688 # Test the error value (cannot be 0.0). 689 if error == 0.0: 690 raise RelaxError("An invalid error value of zero has been encountered.") 691 692 # Get the corresponding spin container. 693 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) 694 spin = return_spin(id) 695 if spin == None and spin_id[0] == '@': # Allow spin IDs of atom names to be used to specify multi column data. 696 spin = return_spin(id+spin_id) 697 if spin == None: 698 warn(RelaxNoSpinWarning(id)) 699 continue 700 701 # Add the data. 702 if data_col: 703 # Initialise. 704 if not hasattr(spin, 'pcs'): 705 spin.pcs = {} 706 707 # Append the value. 708 spin.pcs[align_id] = value 709 710 # Add the error. 711 if error_col: 712 # Initialise. 713 if not hasattr(spin, 'pcs_err'): 714 spin.pcs_err = {} 715 716 # Append the error. 717 spin.pcs_err[align_id] = error 718 719 # Append the data for printout. 720 mol_names.append(mol_name) 721 res_nums.append(res_num) 722 res_names.append(res_name) 723 spin_nums.append(spin_num) 724 spin_names.append(spin_name) 725 values.append(value) 726 errors.append(error) 727 728 # Print out. 729 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') 730 731 732 # Global (non-spin specific) data. 733 ################################## 734 735 # No data, so return. 736 if not len(values): 737 return 738 739 # Initialise. 740 if not hasattr(cdp, 'align_ids'): 741 cdp.align_ids = [] 742 if not hasattr(cdp, 'pcs_ids'): 743 cdp.pcs_ids = [] 744 745 # Add the PCS id string. 746 if align_id not in cdp.align_ids: 747 cdp.align_ids.append(align_id) 748 if align_id not in cdp.pcs_ids: 749 cdp.pcs_ids.append(align_id)

750 751

752 -def return_pcs_data(sim_index=None):

753 """Set up the data structures for optimisation using PCSs as base data sets. 754 755 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 756 @type sim_index: None or int 757 @return: The assembled data structures for using PCSs as the base data for optimisation. These include: 758 - the PCS values. 759 - the unit vectors connecting the paramagnetic centre (the electron spin) to the spin. 760 - the PCS weight. 761 - the experimental temperatures. 762 - the spectrometer frequencies. 763 - pseudo_flags, the list of flags indicating if the interatomic data contains a pseudo-atom (as 1's and 0's). 764 @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) 765 """ 766 767 # Data setup tests. 768 if not hasattr(cdp, 'paramagnetic_centre') and (hasattr(cdp, 'paramag_centre_fixed') and cdp.paramag_centre_fixed): 769 raise RelaxError("The paramagnetic centre has not yet been specified.") 770 if not hasattr(cdp, 'temperature'): 771 raise RelaxError("The experimental temperatures have not been set.") 772 if not hasattr(cdp, 'spectrometer_frq'): 773 raise RelaxError("The spectrometer frequencies of the experiments have not been set.") 774 775 # Sort out pseudo-atoms first. This only needs to be called once. 776 setup_pseudoatom_pcs() 777 778 # Initialise. 779 pcs = [] 780 pcs_err = [] 781 pcs_weight = [] 782 temp = [] 783 frq = [] 784 pseudo_flags = [] 785 786 # The PCS data. 787 for i in range(len(cdp.align_ids)): 788 # Alias the ID. 789 align_id = cdp.align_ids[i] 790 791 # Skip non-optimised data. 792 if not opt_uses_align_data(align_id): 793 continue 794 795 # Append empty arrays to the PCS structures. 796 pcs.append([]) 797 pcs_err.append([]) 798 pcs_weight.append([]) 799 800 # Get the temperature for the PCS constant. 801 if align_id in cdp.temperature: 802 temp.append(cdp.temperature[align_id]) 803 804 # The temperature must be given! 805 else: 806 raise RelaxError("The experimental temperature for the alignment ID '%s' has not been set." % align_id) 807 808 # Get the spectrometer frequency in Tesla units for the PCS constant. 809 if align_id in cdp.spectrometer_frq: 810 frq.append(cdp.spectrometer_frq[align_id] * 2.0 * pi / g1H) 811 812 # The frequency must be given! 813 else: 814 raise RelaxError("The spectrometer frequency for the alignment ID '%s' has not been set." % align_id) 815 816 # Spin loop. 817 j = 0 818 for spin in spin_loop(): 819 # Skip deselected spins. 820 if not spin.select: 821 continue 822 823 # Skip spins without PCS data. 824 if not hasattr(spin, 'pcs'): 825 continue 826 827 # Append the PCSs to the list. 828 if align_id in spin.pcs.keys(): 829 if sim_index != None: 830 pcs[-1].append(spin.pcs_sim[align_id][sim_index]) 831 else: 832 pcs[-1].append(spin.pcs[align_id]) 833 else: 834 pcs[-1].append(None) 835 836 # Append the PCS errors. 837 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 838 pcs_err[-1].append(spin.pcs_err[align_id]) 839 else: 840 pcs_err[-1].append(None) 841 842 # Append the weight. 843 if hasattr(spin, 'pcs_weight') and align_id in spin.pcs_weight.keys(): 844 pcs_weight[-1].append(spin.pcs_weight[align_id]) 845 else: 846 pcs_weight[-1].append(1.0) 847 848 # Spin index. 849 j = j + 1 850 851 # Pseudo-atom. 852 for spin in spin_loop(): 853 if is_pseudoatom(spin): 854 pseudo_flags.append(1) 855 else: 856 pseudo_flags.append(0) 857 858 # Convert to numpy objects. 859 pcs = array(pcs, float64) 860 pcs_err = array(pcs_err, float64) 861 pcs_weight = array(pcs_weight, float64) 862 pseudo_flags = array(pseudo_flags, int32) 863 864 # Convert the PCS from ppm to no units. 865 pcs = pcs * 1e-6 866 pcs_err = pcs_err * 1e-6 867 868 # Return the data structures. 869 return pcs, pcs_err, pcs_weight, temp, frq, pseudo_flags

870 871

872 -def set_errors(align_id=None, spin_id=None, sd=None):

873 """Set the PCS errors if not already present. 874 875 @keyword align_id: The optional alignment tensor ID string. 876 @type align_id: str 877 @keyword spin_id: The optional spin ID string. 878 @type spin_id: None or str 879 @keyword sd: The PCS standard deviation in ppm. 880 @type sd: float or int. 881 """ 882 883 # Check the pipe setup. 884 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 885 886 # Convert the align IDs to an array, or take all IDs. 887 if align_id: 888 align_ids = [align_id] 889 else: 890 align_ids = cdp.pcs_ids 891 892 # Loop over the spins. 893 for spin in spin_loop(spin_id): 894 # Skip deselected spins. 895 if not spin.select: 896 continue 897 898 # Skip spins with no PCSs. 899 if not hasattr(spin, 'pcs') or (align_id and not align_id in spin.pcs): 900 continue 901 902 # No data structure. 903 if not hasattr(spin, 'pcs_err'): 904 spin.pcs_err = {} 905 906 # Set the error. 907 for id in align_ids: 908 spin.pcs_err[id] = sd

909 910

911 -def setup_pseudoatom_pcs():

912 """Make sure that the spin systems are properly set up for pseudo-atoms and PCSs. 913 914 All spin data containers which are a member of a pseudo-atom will be deselected. 915 """ 916 917 # Loop over all spin data containers. 918 for pseudospin, pseudospin_id in spin_loop(return_id=True): 919 # No pseudo-atom, so do nothing. 920 if not is_pseudoatom(pseudospin): 921 return 922 923 # Loop over the atoms of the pseudo-atom. 924 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 925 # Deselect if needed. 926 if spin.select: 927 warn(RelaxWarning("Deselecting the '%s' spin as it is a member of the '%s' pseudo-atom system." % (spin_id, pseudospin_id))) 928 spin.select = False

929 930

931 -def structural_noise(align_id=None, rmsd=0.2, sim_num=1000, file=None, dir=None, force=False):

932 """Determine the PCS error due to structural noise via simulation. 933 934 For the simulation the following must already be set up in the current data pipe: 935 936 - The position of the paramagnetic centre. 937 - The alignment and magnetic susceptibility tensor. 938 939 The protocol for the simulation is as follows: 940 941 - 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. 942 - The positions of the nuclear spins will be randomised N times. For each simulation a random unit vector will be generated. Then a random distance along the unit vector will be generated by sampling from a Gaussian distribution centered at zero, the original spin position, with a standard deviation set to the given RMSD. Both positive and negative displacements will be used. 943 - The PCS for the randomised position will be back calculated. 944 - The PCS standard deviation will be calculated from the N randomised PCS values. 945 946 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. 947 948 If the alignment ID string is not supplied, the procedure will be applied to the PCS data from all alignments. 949 950 951 @keyword align_id: The alignment tensor ID string. 952 @type align_id: str 953 @keyword rmsd: The atomic position RMSD, in Angstrom, to randomise the spin positions with for the simulations. 954 @type rmsd: float 955 @keyword sim_num: The number of simulations, N, to perform to determine the structural noise component of the PCS errors. 956 @type sim_num: int 957 @keyword file: The optional name of the Grace file to plot the structural errors verses the paramagnetic centre to spin distances. 958 @type file: None or str 959 @keyword dir: The directory name to place the Grace file into. 960 @type dir: None or str 961 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 962 @type force: bool 963 """ 964 965 # Check the pipe setup. 966 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True, paramag_centre=True) 967 968 # Convert the align IDs to an array, or take all IDs. 969 if align_id: 970 align_ids = [align_id] 971 else: 972 align_ids = cdp.align_ids 973 974 # Initialise some numpy data structures for use in the simulations. 975 grace_data = [] 976 unit_vect = zeros(3, float64) 977 pcs = {} 978 for id in align_ids: 979 grace_data.append([]) 980 pcs[id] = zeros(sim_num, float64) 981 982 # Print out. 983 print("Executing %i simulations for each spin system." % sim_num) 984 985 # Loop over the spins. 986 for spin, spin_id in spin_loop(return_id=True): 987 # Deselected spins. 988 if not spin.select: 989 continue 990 991 # Skip spins with no PCS or position. 992 if not hasattr(spin, 'pcs'): 993 continue 994 if not hasattr(spin, 'pos'): 995 continue 996 997 # Print out. 998 print(spin_id) 999 1000 # Average the atom position. 1001 if type(spin.pos[0]) in [float, float64]: 1002 pos = spin.pos 1003 else: 1004 pos = zeros(3, float64) 1005 for i in range(len(spin.pos)): 1006 pos += spin.pos[i] 1007 pos = pos / len(spin.pos) 1008 1009 # The original vector length (for the Grace plot). 1010 orig_vect = pos - cdp.paramagnetic_centre 1011 orig_r = norm(orig_vect) 1012 1013 # Loop over the N randomisations. 1014 for i in range(sim_num): 1015 # The random unit vector. 1016 random_unit_vector(unit_vect) 1017 1018 # The random displacement (in Angstrom). 1019 disp = gauss(0, rmsd) 1020 1021 # Move the atom. 1022 new_pos = pos + disp*unit_vect 1023 1024 # The vector and length. 1025 vect = new_pos - cdp.paramagnetic_centre 1026 r = norm(vect) 1027 vect = vect / r 1028 1029 # Loop over the alignments. 1030 for id in align_ids: 1031 # No PCS value, so skip. 1032 if id not in spin.pcs: 1033 continue 1034 1035 # Calculate the PCS constant. 1036 dj = pcs_constant(cdp.temperature[id], cdp.spectrometer_frq[id] * 2.0 * pi / g1H, r/1e10) 1037 1038 # Calculate the PCS value (in ppm). 1039 pcs[id][i] = pcs_tensor(dj, vect, cdp.align_tensors[get_tensor_index(id)].A) * 1e6 1040 1041 # Initialise if necessary. 1042 if not hasattr(spin, 'pcs_struct_err'): 1043 spin.pcs_struct_err = {} 1044 1045 # Loop over the alignments. 1046 align_index = 0 1047 for id in align_ids: 1048 # No PCS value, so skip. 1049 if id not in spin.pcs: 1050 align_index += 1 1051 continue 1052 1053 # The PCS standard deviation. 1054 sd = std(pcs[id]) 1055 1056 # Remove the previous error. 1057 if id in spin.pcs_struct_err: 1058 warn(RelaxWarning("Removing the previous structural error value from the PCS error of the spin '%s' for the alignment ID '%s'." % (spin_id, id))) 1059 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 - spin.pcs_struct_err[id]**2) 1060 1061 # Store the structural error. 1062 spin.pcs_struct_err[id] = sd 1063 1064 # Add it to the PCS error (with variance addition). 1065 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 + sd**2) 1066 1067 # Store the data for the Grace plot. 1068 grace_data[align_index].append([orig_r, sd]) 1069 1070 # Increment the alignment index. 1071 align_index += 1 1072 1073 # The Grace output. 1074 if file: 1075 # Open the Grace file for writing. 1076 file = open_write_file(file, dir, force) 1077 1078 # The header. 1079 grace.write_xy_header(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)"]]) 1080 1081 # The main data. 1082 grace.write_xy_data(data=[grace_data], file=file, graph_type='xy')

1083 1084

1085 -def weight(align_id=None, spin_id=None, weight=1.0):

1086 """Set optimisation weights on the PCS data. 1087 1088 @keyword align_id: The alignment tensor ID string. 1089 @type align_id: str 1090 @keyword spin_id: The spin ID string. 1091 @type spin_id: None or str 1092 @keyword weight: The optimisation weight. The higher the value, the more importance the PCS will have. 1093 @type weight: float or int. 1094 """ 1095 1096 # Check the pipe setup. 1097 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 1098 1099 # Loop over the spins. 1100 for spin in spin_loop(spin_id): 1101 # No data structure. 1102 if not hasattr(spin, 'pcs_weight'): 1103 spin.pcs_weight = {} 1104 1105 # Set the weight. 1106 spin.pcs_weight[align_id] = weight

1107 1108

1109 -def write(align_id=None, file=None, dir=None, bc=False, force=False):

1110 """Display the PCS data corresponding to the alignment ID. 1111 1112 @keyword align_id: The alignment tensor ID string. 1113 @type align_id: str 1114 @keyword file: The file name or object to write to. 1115 @type file: str or file object 1116 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 1117 @type dir: str 1118 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 1119 @type bc: bool 1120 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1121 @type force: bool 1122 """ 1123 1124 # Check the pipe setup. 1125 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 1126 1127 # Open the file for writing. 1128 file = open_write_file(file, dir, force) 1129 1130 # Loop over the spins and collect the data. 1131 mol_names = [] 1132 res_nums = [] 1133 res_names = [] 1134 spin_nums = [] 1135 spin_names = [] 1136 values = [] 1137 errors = [] 1138 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 1139 # Skip deselected spins. 1140 if not spin.select: 1141 continue 1142 1143 # Skip spins with no PCSs. 1144 if not bc and (not hasattr(spin, 'pcs') or not align_id in spin.pcs.keys()): 1145 continue 1146 elif bc and (not hasattr(spin, 'pcs_bc') or align_id not in spin.pcs_bc.keys()): 1147 continue 1148 1149 # Append the spin data. 1150 mol_names.append(mol_name) 1151 res_nums.append(res_num) 1152 res_names.append(res_name) 1153 spin_nums.append(spin.num) 1154 spin_names.append(spin.name) 1155 1156 # The value. 1157 if bc: 1158 values.append(spin.pcs_bc[align_id]) 1159 else: 1160 values.append(spin.pcs[align_id]) 1161 1162 # The error. 1163 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 1164 errors.append(spin.pcs_err[align_id]) 1165 else: 1166 errors.append(None) 1167 1168 # Write out. 1169 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')

1170

Source Code for Module pipe_control.pcs