generic_fns.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, 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 generic_fns import grace, pipes 36 from generic_fns.align_tensor import get_tensor_index 37 from generic_fns.mol_res_spin import exists_mol_res_spin_data, generate_spin_id_unique, return_spin, spin_index_loop, spin_loop 38 from maths_fns.pcs import ave_pcs_tensor, pcs_tensor 39 from maths_fns.vectors import random_unit_vector 40 from physical_constants import g1H, pcs_constant 41 from relax_errors import RelaxError, RelaxAlignError, RelaxNoAlignError, RelaxNoPdbError, RelaxNoPCSError, RelaxNoSequenceError, RelaxPCSError 42 from relax_io import open_write_file, read_spin_data, write_spin_data 43 from relax_warnings import RelaxWarning, RelaxNoSpinWarning 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.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(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", sets=size, set_names=set_names, linestyle=[2]+[0]*size, data_type=['pcs_bc', 'pcs'], axis_min=[-0.5, -0.5], axis_max=[0.5, 0.5], 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 q_factors(spin_id=None):

508 """Calculate the Q-factors for the PCS data. 509 510 @keyword spin_id: The spin ID string used to restrict the Q-factor calculation to a subset of all spins. 511 @type spin_id: None or str 512 """ 513 514 # Check the pipe setup. 515 check_pipe_setup(sequence=True) 516 517 # No PCSs, so no Q factors can be calculated. 518 if not hasattr(cdp, 'pcs_ids') or not len(cdp.pcs_ids): 519 warn(RelaxWarning("No PCS data exists, Q factors cannot be calculated.")) 520 return 521 522 # Q-factor dictionary. 523 cdp.q_factors_pcs = {} 524 525 # Loop over the alignments. 526 for align_id in cdp.pcs_ids: 527 # Init. 528 pcs2_sum = 0.0 529 sse = 0.0 530 531 # Spin loop. 532 spin_count = 0 533 pcs_data = False 534 pcs_bc_data = False 535 for spin in spin_loop(spin_id): 536 # Skip deselected spins. 537 if not spin.select: 538 continue 539 540 # Increment the spin counter. 541 spin_count += 1 542 543 # Data checks. 544 if hasattr(spin, 'pcs') and align_id in spin.pcs: 545 pcs_data = True 546 if hasattr(spin, 'pcs_bc') and align_id in spin.pcs_bc: 547 pcs_bc_data = True 548 549 # Skip spins without PCS data. 550 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: 551 continue 552 553 # Sum of squares. 554 sse = sse + (spin.pcs[align_id] - spin.pcs_bc[align_id])**2 555 556 # Sum the PCSs squared (for normalisation). 557 pcs2_sum = pcs2_sum + spin.pcs[align_id]**2 558 559 # The Q-factor for the alignment. 560 if pcs2_sum: 561 Q = sqrt(sse / pcs2_sum) 562 cdp.q_factors_pcs[align_id] = Q 563 564 # Warnings (and then exit). 565 if not spin_count: 566 warn(RelaxWarning("No spins have been used in the calculation, skipping the PCS Q factor calculation.")) 567 return 568 if not pcs_data: 569 warn(RelaxWarning("No PCS data can be found for the alignment ID '%s', skipping the PCS Q factor calculation for this alignment." % align_id)) 570 continue 571 if not pcs_bc_data: 572 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)) 573 continue 574 575 # The total Q-factor. 576 cdp.q_pcs = 0.0 577 for id in cdp.q_factors_pcs: 578 cdp.q_pcs = cdp.q_pcs + cdp.q_factors_pcs[id]**2 579 cdp.q_pcs = cdp.q_pcs / len(cdp.q_factors_pcs) 580 cdp.q_pcs = sqrt(cdp.q_pcs)

581 582

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

584 """Read the PCS data from file. 585 586 @param align_id: The alignment tensor ID string. 587 @type align_id: str 588 @param file: The name of the file to open. 589 @type file: str 590 @param dir: The directory containing the file (defaults to the current directory if None). 591 @type dir: str or None 592 @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. 593 @type file_data: list of lists 594 @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. 595 @type spin_id_col: int or None 596 @keyword mol_name_col: The column containing the molecule name information. If supplied, spin_id_col must be None. 597 @type mol_name_col: int or None 598 @keyword res_name_col: The column containing the residue name information. If supplied, spin_id_col must be None. 599 @type res_name_col: int or None 600 @keyword res_num_col: The column containing the residue number information. If supplied, spin_id_col must be None. 601 @type res_num_col: int or None 602 @keyword spin_name_col: The column containing the spin name information. If supplied, spin_id_col must be None. 603 @type spin_name_col: int or None 604 @keyword spin_num_col: The column containing the spin number information. If supplied, spin_id_col must be None. 605 @type spin_num_col: int or None 606 @keyword data_col: The column containing the PCS data in Hz. 607 @type data_col: int or None 608 @keyword error_col: The column containing the PCS errors. 609 @type error_col: int or None 610 @keyword sep: The column separator which, if None, defaults to whitespace. 611 @type sep: str or None 612 @keyword spin_id: The spin ID string used to restrict data loading to a subset of all spins. 613 @type spin_id: None or str 614 """ 615 616 # Check the pipe setup. 617 check_pipe_setup(sequence=True) 618 619 # Test if sequence data exists. 620 if not exists_mol_res_spin_data(): 621 raise RelaxNoSequenceError 622 623 # Either the data or error column must be supplied. 624 if data_col == None and error_col == None: 625 raise RelaxError("One of either the data or error column must be supplied.") 626 627 628 # Spin specific data. 629 ##################### 630 631 # Loop over the PCS data. 632 mol_names = [] 633 res_nums = [] 634 res_names = [] 635 spin_nums = [] 636 spin_names = [] 637 values = [] 638 errors = [] 639 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): 640 # Unpack. 641 if data_col and error_col: 642 mol_name, res_num, res_name, spin_num, spin_name, value, error = data 643 elif data_col: 644 mol_name, res_num, res_name, spin_num, spin_name, value = data 645 error = None 646 else: 647 mol_name, res_num, res_name, spin_num, spin_name, error = data 648 value = None 649 650 # Test the error value (cannot be 0.0). 651 if error == 0.0: 652 raise RelaxError("An invalid error value of zero has been encountered.") 653 654 # Get the corresponding spin container. 655 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) 656 spin = return_spin(id) 657 if spin == None and spin_id[0] == '@': # Allow spin IDs of atom names to be used to specify multi column data. 658 spin = return_spin(id+spin_id) 659 if spin == None: 660 warn(RelaxNoSpinWarning(id)) 661 continue 662 663 # Add the data. 664 if data_col: 665 # Initialise. 666 if not hasattr(spin, 'pcs'): 667 spin.pcs = {} 668 669 # Append the value. 670 spin.pcs[align_id] = value 671 672 # Add the error. 673 if error_col: 674 # Initialise. 675 if not hasattr(spin, 'pcs_err'): 676 spin.pcs_err = {} 677 678 # Append the error. 679 spin.pcs_err[align_id] = error 680 681 # Append the data for printout. 682 mol_names.append(mol_name) 683 res_nums.append(res_num) 684 res_names.append(res_name) 685 spin_nums.append(spin_num) 686 spin_names.append(spin_name) 687 values.append(value) 688 errors.append(error) 689 690 # Print out. 691 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') 692 693 694 # Global (non-spin specific) data. 695 ################################## 696 697 # No data, so return. 698 if not len(values): 699 return 700 701 # Initialise. 702 if not hasattr(cdp, 'align_ids'): 703 cdp.align_ids = [] 704 if not hasattr(cdp, 'pcs_ids'): 705 cdp.pcs_ids = [] 706 707 # Add the PCS id string. 708 if align_id not in cdp.align_ids: 709 cdp.align_ids.append(align_id) 710 if align_id not in cdp.pcs_ids: 711 cdp.pcs_ids.append(align_id)

712 713

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

715 """Set the PCS errors if not already present. 716 717 @keyword align_id: The optional alignment tensor ID string. 718 @type align_id: str 719 @keyword spin_id: The optional spin ID string. 720 @type spin_id: None or str 721 @keyword sd: The PCS standard deviation in ppm. 722 @type sd: float or int. 723 """ 724 725 # Check the pipe setup. 726 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 727 728 # Convert the align IDs to an array, or take all IDs. 729 if align_id: 730 align_ids = [align_id] 731 else: 732 align_ids = cdp.pcs_ids 733 734 # Loop over the spins. 735 for spin in spin_loop(spin_id): 736 # Skip deselected spins. 737 if not spin.select: 738 continue 739 740 # Skip spins with no PCSs. 741 if not hasattr(spin, 'pcs') or (align_id and not align_id in spin.pcs): 742 continue 743 744 # No data structure. 745 if not hasattr(spin, 'pcs_err'): 746 spin.pcs_err = {} 747 748 # Set the error. 749 for id in align_ids: 750 spin.pcs_err[id] = sd

751 752

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

754 """Determine the PCS error due to structural noise via simulation. 755 756 For the simulation the following must already be set up in the current data pipe: 757 758 - The position of the paramagnetic centre. 759 - The alignment and magnetic susceptibility tensor. 760 761 The protocol for the simulation is as follows: 762 763 - 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. 764 - 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. 765 - The PCS for the randomised position will be back calculated. 766 - The PCS standard deviation will be calculated from the N randomised PCS values. 767 768 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. 769 770 If the alignment ID string is not supplied, the procedure will be applied to the PCS data from all alignments. 771 772 773 @keyword align_id: The alignment tensor ID string. 774 @type align_id: str 775 @keyword rmsd: The atomic position RMSD, in Angstrom, to randomise the spin positions with for the simulations. 776 @type rmsd: float 777 @keyword sim_num: The number of simulations, N, to perform to determine the structural noise component of the PCS errors. 778 @type sim_num: int 779 @keyword file: The optional name of the Grace file to plot the structural errors verses the paramagnetic centre to spin distances. 780 @type file: None or str 781 @keyword dir: The directory name to place the Grace file into. 782 @type dir: None or str 783 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 784 @type force: bool 785 """ 786 787 # Check the pipe setup. 788 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True, paramag_centre=True) 789 790 # Convert the align IDs to an array, or take all IDs. 791 if align_id: 792 align_ids = [align_id] 793 else: 794 align_ids = cdp.align_ids 795 796 # Initialise some numpy data structures for use in the simulations. 797 grace_data = [] 798 unit_vect = zeros(3, float64) 799 pcs = {} 800 for id in align_ids: 801 grace_data.append([]) 802 pcs[id] = zeros(sim_num, float64) 803 804 # Print out. 805 print("Executing %i simulations for each spin system." % sim_num) 806 807 # Loop over the spins. 808 for spin, spin_id in spin_loop(return_id=True): 809 # Deselected spins. 810 if not spin.select: 811 continue 812 813 # Skip spins with no PCS or position. 814 if not hasattr(spin, 'pcs'): 815 continue 816 if not hasattr(spin, 'pos'): 817 continue 818 819 # Print out. 820 print(spin_id) 821 822 # Average the atom position. 823 if type(spin.pos[0]) in [float, float64]: 824 pos = spin.pos 825 else: 826 pos = zeros(3, float64) 827 for i in range(len(spin.pos)): 828 pos += spin.pos[i] 829 pos = pos / len(spin.pos) 830 831 # The original vector length (for the Grace plot). 832 orig_vect = pos - cdp.paramagnetic_centre 833 orig_r = norm(orig_vect) 834 835 # Loop over the N randomisations. 836 for i in range(sim_num): 837 # The random unit vector. 838 random_unit_vector(unit_vect) 839 840 # The random displacement (in Angstrom). 841 disp = gauss(0, rmsd) 842 843 # Move the atom. 844 new_pos = pos + disp*unit_vect 845 846 # The vector and length. 847 vect = new_pos - cdp.paramagnetic_centre 848 r = norm(vect) 849 vect = vect / r 850 851 # Loop over the alignments. 852 for id in align_ids: 853 # No PCS value, so skip. 854 if id not in spin.pcs: 855 continue 856 857 # Calculate the PCS constant. 858 dj = pcs_constant(cdp.temperature[id], cdp.frq[id] * 2.0 * pi / g1H, r/1e10) 859 860 # Calculate the PCS value (in ppm). 861 pcs[id][i] = pcs_tensor(dj, vect, cdp.align_tensors[get_tensor_index(id)].A) * 1e6 862 863 # Initialise if necessary. 864 if not hasattr(spin, 'pcs_struct_err'): 865 spin.pcs_struct_err = {} 866 867 # Loop over the alignments. 868 align_index = 0 869 for id in align_ids: 870 # No PCS value, so skip. 871 if id not in spin.pcs: 872 align_index += 1 873 continue 874 875 # The PCS standard deviation. 876 sd = std(pcs[id]) 877 878 # Remove the previous error. 879 if id in spin.pcs_struct_err: 880 warn(RelaxWarning("Removing the previous structural error value from the PCS error of the spin '%s' for the alignment ID '%s'." % (spin_id, id))) 881 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 - spin.pcs_struct_err[id]**2) 882 883 # Store the structural error. 884 spin.pcs_struct_err[id] = sd 885 886 # Add it to the PCS error (with variance addition). 887 spin.pcs_err[id] = sqrt(spin.pcs_err[id]**2 + sd**2) 888 889 # Store the data for the Grace plot. 890 grace_data[align_index].append([orig_r, sd]) 891 892 # Increment the alignment index. 893 align_index += 1 894 895 # The Grace output. 896 if file: 897 # Open the Grace file for writing. 898 file = open_write_file(file, dir, force) 899 900 # The header. 901 grace.write_xy_header(file=file, title="PCS structural noise", subtitle="%s Angstrom structural noise"%rmsd, sets=len(align_ids), set_names=align_ids, symbol_sizes=[0.5]*len(align_ids), linetype=[0]*len(align_ids), data_type=['pcs_bc', 'pcs'], axis_labels=["Ln\\S3+\\N to spin distance (Angstrom)", "PCS standard deviation (ppm)"]) 902 903 # The main data. 904 grace.write_xy_data(data=[grace_data], file=file, graph_type='xy')

905 906

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

908 """Set optimisation weights on the PCS data. 909 910 @keyword align_id: The alignment tensor ID string. 911 @type align_id: str 912 @keyword spin_id: The spin ID string. 913 @type spin_id: None or str 914 @keyword weight: The optimisation weight. The higher the value, the more importance the PCS will have. 915 @type weight: float or int. 916 """ 917 918 # Check the pipe setup. 919 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 920 921 # Loop over the spins. 922 for spin in spin_loop(spin_id): 923 # No data structure. 924 if not hasattr(spin, 'pcs_weight'): 925 spin.pcs_weight = {} 926 927 # Set the weight. 928 spin.pcs_weight[align_id] = weight

929 930

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

932 """Display the PCS data corresponding to the alignment ID. 933 934 @keyword align_id: The alignment tensor ID string. 935 @type align_id: str 936 @keyword file: The file name or object to write to. 937 @type file: str or file object 938 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 939 @type dir: str 940 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 941 @type bc: bool 942 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 943 @type force: bool 944 """ 945 946 # Check the pipe setup. 947 check_pipe_setup(sequence=True, pcs_id=align_id, pcs=True) 948 949 # Open the file for writing. 950 file = open_write_file(file, dir, force) 951 952 # Loop over the spins and collect the data. 953 mol_names = [] 954 res_nums = [] 955 res_names = [] 956 spin_nums = [] 957 spin_names = [] 958 values = [] 959 errors = [] 960 for spin, mol_name, res_num, res_name in spin_loop(full_info=True): 961 # Skip deselected spins. 962 if not spin.select: 963 continue 964 965 # Skip spins with no PCSs. 966 if not bc and (not hasattr(spin, 'pcs') or not align_id in spin.pcs.keys()): 967 continue 968 elif bc and (not hasattr(spin, 'pcs_bc') or align_id not in spin.pcs_bc.keys()): 969 continue 970 971 # Append the spin data. 972 mol_names.append(mol_name) 973 res_nums.append(res_num) 974 res_names.append(res_name) 975 spin_nums.append(spin.num) 976 spin_names.append(spin.name) 977 978 # The value. 979 if bc: 980 values.append(spin.pcs_bc[align_id]) 981 else: 982 values.append(spin.pcs[align_id]) 983 984 # The error. 985 if hasattr(spin, 'pcs_err') and align_id in spin.pcs_err.keys(): 986 errors.append(spin.pcs_err[align_id]) 987 else: 988 errors.append(None) 989 990 # Write out. 991 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')

992

Source Code for Module generic_fns.pcs