pipe_control.rdc

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 RDC data.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi, sqrt 28 from numpy import array, int32, float64, ones, transpose, zeros 29 from numpy.linalg import norm 30 import sys 31 from warnings import warn 32 33 # relax module imports. 34 from lib.check_types import is_float 35 from lib.float import nan 36 from lib.alignment.rdc import ave_rdc_tensor 37 from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoJError, RelaxNoRDCError, RelaxNoSequenceError, RelaxSpinTypeError 38 from lib.io import extract_data, open_write_file, strip, write_data 39 from lib.physical_constants import dipolar_constant, return_gyromagnetic_ratio 40 from lib.warnings import RelaxWarning, RelaxSpinTypeWarning 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.interatomic import consistent_interatomic_data, create_interatom, interatomic_loop, return_interatom 44 from pipe_control.mol_res_spin import exists_mol_res_spin_data, is_pseudoatom, pseudoatom_loop, return_spin 45 46

47 -def back_calc(align_id=None):

48 """Back calculate the RDC from the alignment tensor and unit bond vectors. 49 50 @keyword align_id: The alignment tensor ID string. 51 @type align_id: str 52 """ 53 54 # Check the pipe setup. 55 check_pipe_setup(rdc_id=align_id, sequence=True, N=True, tensors=True) 56 57 # Convert the align IDs to an array, or take all IDs. 58 if align_id: 59 align_ids = [align_id] 60 else: 61 align_ids = cdp.align_ids 62 63 # Add the ID to the RDC IDs, if needed. 64 for align_id in align_ids: 65 # Init. 66 if not hasattr(cdp, 'rdc_ids'): 67 cdp.rdc_ids = [] 68 69 # Add the ID. 70 if align_id not in cdp.rdc_ids: 71 cdp.rdc_ids.append(align_id) 72 73 # The weights. 74 weights = ones(cdp.N, float64) / cdp.N 75 76 # Unit vector data structure init. 77 unit_vect = zeros((cdp.N, 3), float64) 78 79 # Loop over the interatomic data. 80 count = 0 81 for interatom in interatomic_loop(): 82 # Skip containers with no interatomic vectors. 83 if not hasattr(interatom, 'vector'): 84 continue 85 86 # Get the spins. 87 spin1 = return_spin(interatom.spin_id1) 88 spin2 = return_spin(interatom.spin_id2) 89 90 # Checks. 91 if not hasattr(spin1, 'isotope'): 92 raise RelaxSpinTypeError(interatom.spin_id1) 93 if not hasattr(spin2, 'isotope'): 94 raise RelaxSpinTypeError(interatom.spin_id2) 95 96 # Single vector. 97 if is_float(interatom.vector[0]): 98 vectors = [interatom.vector] 99 else: 100 vectors = interatom.vector 101 102 # Gyromagnetic ratios. 103 g1 = return_gyromagnetic_ratio(spin1.isotope) 104 g2 = return_gyromagnetic_ratio(spin2.isotope) 105 106 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 107 dj = 3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r) 108 109 # Unit vectors. 110 for c in range(cdp.N): 111 unit_vect[c] = vectors[c] / norm(vectors[c]) 112 113 # Initialise if necessary. 114 if not hasattr(interatom, 'rdc_bc'): 115 interatom.rdc_bc = {} 116 117 # Calculate the RDCs. 118 for id in align_ids: 119 # The signed value. 120 interatom.rdc_bc[id] = ave_rdc_tensor(dj, unit_vect, cdp.N, cdp.align_tensors[get_tensor_index(align_id=id)].A, weights=weights) 121 122 # T values. 123 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 124 if not hasattr(interatom, 'j_coupling'): 125 raise RelaxNoJError 126 127 interatom.rdc_bc[id] += interatom.j_coupling 128 129 # The absolute value. 130 if hasattr(interatom, 'absolute_rdc') and id in interatom.absolute_rdc.keys() and interatom.absolute_rdc[id]: 131 interatom.rdc_bc[id] = abs(interatom.rdc_bc[id]) 132 133 # Increment the counter. 134 count += 1 135 136 # No RDCs calculated. 137 if not count: 138 warn(RelaxWarning("No RDCs have been back calculated, probably due to missing bond vector information."))

139 140

141 -def check_pipe_setup(pipe=None, rdc_id=None, sequence=False, N=False, tensors=False, rdc=False):

142 """Check that the current data pipe has been setup sufficiently. 143 144 @keyword pipe: The data pipe to check, defaulting to the current pipe. 145 @type pipe: None or str 146 @keyword rdc_id: The RDC ID string to check for in cdp.rdc_ids. 147 @type rdc_id: None or str 148 @keyword sequence: A flag which when True will invoke the sequence data check. 149 @type sequence: bool 150 @keyword N: A flag which if True will check that cdp.N is set. 151 @type N: bool 152 @keyword tensors: A flag which if True will check that alignment tensors exist. 153 @type tensors: bool 154 @keyword rdc: A flag which if True will check that RDCs exist. 155 @type rdc: bool 156 """ 157 158 # The data pipe. 159 if pipe == None: 160 pipe = pipes.cdp_name() 161 162 # Get the data pipe. 163 dp = pipes.get_pipe(pipe) 164 165 # Test if the current data pipe exists. 166 pipes.test(pipe) 167 168 # Test if sequence data exists. 169 if sequence and not exists_mol_res_spin_data(pipe): 170 raise RelaxNoSequenceError(pipe) 171 172 # Check for dp.N. 173 if N and not hasattr(dp, 'N'): 174 raise RelaxError("The number of states N has not been set.") 175 176 # Check that alignment tensors are present. 177 if tensors and (not hasattr(dp, 'align_tensors') or len(dp.align_tensors) == 0): 178 raise RelaxNoTensorError('alignment') 179 180 # Test for the alignment ID. 181 if rdc_id and (not hasattr(dp, 'align_ids') or rdc_id not in dp.align_ids): 182 raise RelaxNoAlignError(rdc_id, pipe) 183 184 # Test if RDC data exists. 185 if rdc and not hasattr(dp, 'align_ids'): 186 raise RelaxNoAlignError() 187 if rdc and not hasattr(dp, 'rdc_ids'): 188 raise RelaxNoRDCError() 189 elif rdc and rdc_id and rdc_id not in dp.rdc_ids: 190 raise RelaxNoRDCError(rdc_id)

191 192

193 -def check_rdcs(interatom):

194 """Check if all data required for calculating the RDC is present. 195 196 @param interatom: The interatomic data container. 197 @type interatom: InteratomContainer instance 198 @return: True if all data required for calculating the RDC is present, False otherwise. 199 @rtype: bool 200 """ 201 202 # Skip deselected interatomic data containers. 203 if not interatom.select: 204 return False 205 206 # Only use interatomic data containers with RDC data. 207 if not hasattr(interatom, 'rdc'): 208 return False 209 210 # Only use interatomic data containers with RDC and J coupling data. 211 if opt_uses_j_couplings() and not hasattr(interatom, 'j_coupling'): 212 return False 213 214 # Get the spins. 215 spin1 = return_spin(interatom.spin_id1) 216 spin2 = return_spin(interatom.spin_id2) 217 218 # Spin information checks. 219 if not hasattr(spin1, 'isotope'): 220 warn(RelaxSpinTypeWarning(interatom.spin_id1)) 221 return False 222 if not hasattr(spin2, 'isotope'): 223 warn(RelaxSpinTypeWarning(interatom.spin_id2)) 224 return False 225 if is_pseudoatom(spin1) and is_pseudoatom(spin2): 226 warn(RelaxWarning("Support for both spins being in a dipole pair being pseudo-atoms is not implemented yet.")) 227 return False 228 229 # Pseudo-atom checks. 230 if is_pseudoatom(spin1) or is_pseudoatom(spin2): 231 # Alias the pseudo and normal atoms. 232 pseudospin = spin1 233 base_spin_id = interatom.spin_id2 234 pseudospin_id = interatom.spin_id1 235 if is_pseudoatom(spin2): 236 pseudospin = spin2 237 base_spin_id = interatom.spin_id1 238 pseudospin_id = interatom.spin_id2 239 240 # Loop over the atoms of the pseudo-atom. 241 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 242 # Get the corresponding interatomic data container. 243 pseudo_interatom = return_interatom(spin_id1=spin_id, spin_id2=base_spin_id) 244 245 # Unit vector check. 246 if not hasattr(pseudo_interatom, 'vector'): 247 warn(RelaxWarning("The interatomic vector is missing for the spin pair '%s' and '%s' of the pseudo-atom '%s', skipping the RDC for the spin pair '%s' and '%s'." % (pseudo_interatom.spin_id1, pseudo_interatom.spin_id2, pseudospin_id, base_spin_id, pseudospin_id))) 248 return False 249 250 # Check. 251 if not hasattr(pseudo_interatom, 'r'): 252 warn(RelaxWarning("The averaged interatomic distance between spins '%s' and '%s' for the pseudo-atom '%s' has not been set yet." % (spin_id, base_spin_id, pseudospin_id))) 253 return False 254 255 # Normal atoms checks. 256 else: 257 # Unit vector check. 258 if not hasattr(interatom, 'vector'): 259 warn(RelaxWarning("The interatomic vector is missing, skipping the spin pair '%s' and '%s'." % (interatom.spin_id1, interatom.spin_id2))) 260 return False 261 262 # Distance information check. 263 if not hasattr(interatom, 'r'): 264 warn(RelaxWarning("The averaged interatomic distance between spins '%s' and '%s' has not been set yet." % (interatom.spin_id1, interatom.spin_id2))) 265 return False 266 267 # Everything is ok. 268 return True

269 270

271 -def convert(value, data_type, align_id, to_intern=False):

272 """Convert the RDC based on the 'D' or '2D' data type. 273 274 @param value: The value or error to convert. 275 @type value: float or None 276 @param data_type: The RDC data type. Either 'D', '2D' or 'T'. 277 @type data_type: str 278 @param align_id: The alignment tensor ID string. 279 @type align_id: str 280 @keyword to_intern: A flag which if True will convert to the internal D notation if needed, or if False will convert from the internal D notation to the external D or 2D format. 281 @type to_intern: bool 282 @return: The converted value. 283 @rtype: float or None 284 """ 285 286 # Handle values of None. 287 if value == None: 288 return None 289 290 # The conversion factor. 291 factor = 1.0 292 if data_type == '2D': 293 # Convert from 2D to D. 294 if to_intern: 295 factor = 2.0 296 297 # Convert from D to 2D. 298 else: 299 factor = 0.5 300 301 # Return the converted value. 302 return value * factor

303 304

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

306 """Copy the RDC data from one data pipe to another. 307 308 @keyword pipe_from: The data pipe to copy the RDC data from. This defaults to the current data pipe. 309 @type pipe_from: str 310 @keyword pipe_to: The data pipe to copy the RDC data to. This defaults to the current data pipe. 311 @type pipe_to: str 312 @param align_id: The alignment ID string. 313 @type align_id: str 314 """ 315 316 # Defaults. 317 if pipe_from == None and pipe_to == None: 318 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 319 elif pipe_from == None: 320 pipe_from = pipes.cdp_name() 321 elif pipe_to == None: 322 pipe_to = pipes.cdp_name() 323 324 # Check the pipe setup. 325 check_pipe_setup(pipe=pipe_from, rdc_id=align_id, sequence=True, rdc=True) 326 check_pipe_setup(pipe=pipe_to, sequence=True) 327 328 # Get the data pipes. 329 dp_from = pipes.get_pipe(pipe_from) 330 dp_to = pipes.get_pipe(pipe_to) 331 332 # Test that the interatomic data is consistent between the two data pipe. 333 consistent_interatomic_data(pipe1=pipe_to, pipe2=pipe_from) 334 335 # The IDs. 336 if align_id == None: 337 align_ids = dp_from.align_ids 338 else: 339 align_ids = [align_id] 340 341 # Init target pipe global structures. 342 if not hasattr(dp_to, 'align_ids'): 343 dp_to.align_ids = [] 344 if not hasattr(dp_to, 'rdc_ids'): 345 dp_to.rdc_ids = [] 346 347 # Loop over the align IDs. 348 for align_id in align_ids: 349 # Copy the global data. 350 if align_id not in dp_to.align_ids and align_id not in dp_to.align_ids: 351 dp_to.align_ids.append(align_id) 352 if align_id in dp_from.rdc_ids and align_id not in dp_to.rdc_ids: 353 dp_to.rdc_ids.append(align_id) 354 355 # Loop over the interatomic data. 356 for i in range(len(dp_from.interatomic)): 357 # Alias the containers. 358 interatom_from = dp_from.interatomic[i] 359 interatom_to = dp_to.interatomic[i] 360 361 # No data or errors. 362 if (not hasattr(interatom_from, 'rdc') or not align_id in interatom_from.rdc) and (not hasattr(interatom_from, 'rdc_err') or not align_id in interatom_from.rdc_err): 363 continue 364 365 # Initialise the data structures if necessary. 366 if hasattr(interatom_from, 'rdc') and not hasattr(interatom_to, 'rdc'): 367 interatom_to.rdc = {} 368 if hasattr(interatom_from, 'rdc_err') and not hasattr(interatom_to, 'rdc_err'): 369 interatom_to.rdc_err = {} 370 371 # Copy the value and error from pipe_from. 372 if hasattr(interatom_from, 'rdc'): 373 interatom_to.rdc[align_id] = interatom_from.rdc[align_id] 374 if hasattr(interatom_from, 'rdc_err'): 375 interatom_to.rdc_err[align_id] = interatom_from.rdc_err[align_id]

376 377

378 -def corr_plot(format=None, title=None, subtitle=None, file=None, dir=None, force=False):

379 """Generate a correlation plot of the measured vs. back-calculated RDCs. 380 381 @keyword format: The format for the plot file. The following values are accepted: 'grace', a Grace plot; None, a plain text file. 382 @type format: str or None 383 @keyword title: The title for the plot, overriding the default. 384 @type title: None or str 385 @keyword subtitle: The subtitle for the plot, overriding the default. 386 @type subtitle: None or str 387 @keyword file: The file name or object to write to. 388 @type file: str or file object 389 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 390 @type dir: str 391 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 392 @type force: bool 393 """ 394 395 # Check the pipe setup. 396 check_pipe_setup(sequence=True) 397 398 # Does RDC data exist? 399 if not hasattr(cdp, 'rdc_ids') or not cdp.rdc_ids: 400 warn(RelaxWarning("No RDC data exists, skipping file creation.")) 401 return 402 403 # Open the file for writing. 404 file = open_write_file(file, dir, force) 405 406 # Init. 407 data = [] 408 orig_title = title 409 if orig_title == None: 410 title = "RDC correlation plot" 411 axis_labels = ["Measured RDC (Hz)", "Back-calculated RDC (Hz)"] 412 413 # The diagonal. 414 data.append([[-100, -100, 0], [100, 100, 0]]) 415 416 # Loop over the RDC data. 417 for align_id in cdp.rdc_ids: 418 # Append a new list for this alignment. 419 data.append([]) 420 421 # Errors present? 422 err_flag = False 423 for interatom in interatomic_loop(): 424 # Error present. 425 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 426 err_flag = True 427 break 428 429 # Loop over the interatomic data. 430 for interatom in interatomic_loop(): 431 # Skip if data is missing. 432 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'rdc_bc') or not align_id in interatom.rdc.keys() or not align_id in interatom.rdc_bc.keys(): 433 continue 434 435 # Convert between the 2D and D notation. 436 rdc_bc = convert(interatom.rdc_bc[align_id], interatom.rdc_data_types[align_id], align_id) 437 rdc = convert(interatom.rdc[align_id], interatom.rdc_data_types[align_id], align_id) 438 439 # T=J+D type data. 440 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 441 # Convert T=J+D type data to D, if not absolute. 442 if not interatom.absolute_rdc[align_id]: 443 rdc_bc -= interatom.j_coupling 444 rdc -= interatom.j_coupling 445 446 # A different title and axis labels. 447 else: 448 if orig_title == None: 449 title = "T = J+D correlation plot" 450 axis_labels = ["Measured T = J+D (Hz)", "Back-calculated T = J+D (Hz)"] 451 452 # Append the data. 453 data[-1].append([rdc, rdc_bc]) 454 455 # Errors. 456 if err_flag: 457 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 458 data[-1][-1].append(convert(interatom.rdc_err[align_id], interatom.rdc_data_types[align_id], align_id)) 459 else: 460 data[-1][-1].append(None) 461 462 # Label. 463 data[-1][-1].append("%s-%s" % (interatom.spin_id1, interatom.spin_id2)) 464 465 # The data size. 466 size = len(data) 467 468 # Only one data set. 469 data = [data] 470 471 # Graph type. 472 if err_flag: 473 graph_type = 'xydy' 474 else: 475 graph_type = 'xy' 476 477 # Grace file. 478 if format == 'grace': 479 # The header. 480 grace.write_xy_header(file=file, title=title, subtitle=subtitle, world=[[-50, -50, 50, 50]], sets=[size], set_names=[[None]+cdp.rdc_ids], linestyle=[[2]+[0]*size], data_type=['rdc', 'rdc_bc'], axis_labels=[axis_labels], tick_major_spacing=[[10, 10]], tick_minor_count=[[9, 9]], legend_pos=[[1, 0.5]]) 481 482 # The main data. 483 grace.write_xy_data(data=data, file=file, graph_type=graph_type, autoscale=False)

484 485

486 -def delete(align_id=None):

487 """Delete the RDC data corresponding to the alignment ID. 488 489 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted. 490 @type align_id: str or None 491 """ 492 493 # Check the pipe setup. 494 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 495 496 # The IDs. 497 if not align_id: 498 ids = deepcopy(cdp.rdc_ids) 499 else: 500 ids = [align_id] 501 502 # Loop over the alignments, removing all the corresponding data. 503 for id in ids: 504 # The RDC ID. 505 cdp.rdc_ids.pop(cdp.rdc_ids.index(id)) 506 507 # The interatomic data. 508 for interatom in interatomic_loop(): 509 # The data. 510 if hasattr(interatom, 'rdc') and id in interatom.rdc: 511 interatom.rdc.pop(id) 512 513 # The error. 514 if hasattr(interatom, 'rdc_err') and id in interatom.rdc_err: 515 interatom.rdc_err.pop(id) 516 517 # The data type. 518 if hasattr(interatom, 'rdc_data_types') and id in interatom.rdc_data_types: 519 interatom.rdc_data_types.pop(id) 520 521 # Clean the global data. 522 if not hasattr(cdp, 'pcs_ids') or id not in cdp.pcs_ids: 523 cdp.align_ids.pop(cdp.align_ids.index(id))

524 525

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

527 """Display the RDC data corresponding to the alignment ID. 528 529 @keyword align_id: The alignment tensor ID string. 530 @type align_id: str 531 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed. 532 @type bc: bool 533 """ 534 535 # Check the pipe setup. 536 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 537 538 # Call the write method with sys.stdout as the file. 539 write(align_id=align_id, file=sys.stdout, bc=bc)

540 541

542 -def opt_uses_j_couplings():

543 """Determine of J couplings are needed for optimisation. 544 545 @return: True if J couplings are required, False otherwise. 546 @rtype: bool 547 """ 548 549 # Loop over the alignments. 550 for align_id in cdp.align_ids: 551 for interatom in interatomic_loop(): 552 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 553 return True 554 555 # No J values required. 556 return False

557 558

559 -def opt_uses_rdc(align_id):

560 """Determine if the RDC data for the given alignment ID is needed for optimisation. 561 562 @param align_id: The alignment ID string. 563 @type align_id: str 564 @return: True if the RDC data is to be used for optimisation, False otherwise. 565 @rtype: bool 566 """ 567 568 # No alignment IDs. 569 if not hasattr(cdp, 'rdc_ids'): 570 return False 571 572 # No RDC data for the alignment. 573 if align_id not in cdp.rdc_ids: 574 return False 575 576 # Is the tensor optimised? 577 tensor_flag = opt_uses_tensor(get_tensor_object(align_id)) 578 579 # Is the paramagnetic position optimised? 580 pos_flag = False 581 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 582 pos_flag = True 583 584 # Are the populations optimised? 585 prob_flag = False 586 if cdp.model == 'population': 587 prob_flag = True 588 589 # Not used. 590 if not tensor_flag and not pos_flag and not prob_flag: 591 return False 592 593 # The RDC data is to be used for optimisation. 594 return True

595 596

597 -def q_factors(spin_id=None):

598 """Calculate the Q factors for the RDC data. 599 600 @keyword spin_id: The spin ID string used to restrict the Q factor calculation to a subset of all spins. 601 @type spin_id: None or str 602 """ 603 604 # Check the pipe setup. 605 check_pipe_setup(sequence=True) 606 607 # No RDCs, so no Q factors can be calculated. 608 if not hasattr(cdp, 'rdc_ids') or not len(cdp.rdc_ids): 609 warn(RelaxWarning("No RDC data exists, Q factors cannot be calculated.")) 610 return 611 612 # Q factor dictonaries. 613 cdp.q_factors_rdc = {} 614 cdp.q_factors_rdc_norm2 = {} 615 616 # Loop over the alignments. 617 for align_id in cdp.rdc_ids: 618 # Init. 619 D2_sum = 0.0 620 sse = 0.0 621 622 # Interatomic data loop. 623 dj = None 624 N = 0 625 interatom_count = 0 626 rdc_data = False 627 rdc_bc_data = False 628 norm2_flag = True 629 for interatom in interatomic_loop(): 630 # Increment the counter. 631 interatom_count += 1 632 633 # Data checks. 634 if hasattr(interatom, 'rdc') and align_id in interatom.rdc: 635 rdc_data = True 636 if hasattr(interatom, 'rdc_bc') and align_id in interatom.rdc_bc: 637 rdc_bc_data = True 638 j_flag = False 639 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 640 j_flag = True 641 if not hasattr(interatom, 'j_coupling'): 642 raise RelaxNoJError 643 644 # Skip containers without RDC data. 645 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'rdc_bc') or not align_id in interatom.rdc or interatom.rdc[align_id] == None or not align_id in interatom.rdc_bc or interatom.rdc_bc[align_id] == None: 646 continue 647 648 # Get the spins. 649 spin1 = return_spin(interatom.spin_id1) 650 spin2 = return_spin(interatom.spin_id2) 651 652 # Sum of squares. 653 sse = sse + (interatom.rdc[align_id] - interatom.rdc_bc[align_id])**2 654 655 # Sum the RDCs squared (for one type of normalisation). 656 if j_flag: 657 D2_sum = D2_sum + (interatom.rdc[align_id] - interatom.j_coupling)**2 658 else: 659 D2_sum = D2_sum + interatom.rdc[align_id]**2 660 661 # Gyromagnetic ratios. 662 g1 = return_gyromagnetic_ratio(spin1.isotope) 663 g2 = return_gyromagnetic_ratio(spin2.isotope) 664 665 # Skip the 2Da^2(4 + 3R)/5 normalised Q factor if pseudo-atoms are present. 666 if norm2_flag and (is_pseudoatom(spin1) or is_pseudoatom(spin2)): 667 warn(RelaxWarning("Pseudo-atoms are present, skipping the Q factor normalised with 2Da^2(4 + 3R)/5.")) 668 norm2_flag = False 669 670 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 671 if norm2_flag: 672 dj_new = 3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r) 673 if dj != None and dj_new != dj: 674 warn(RelaxWarning("The dipolar constant is not the same for all RDCs, skipping the Q factor normalised with 2Da^2(4 + 3R)/5.")) 675 norm2_flag = False 676 else: 677 dj = dj_new 678 679 # Increment the number of data sets. 680 N = N + 1 681 682 # Warnings (and then exit). 683 if not interatom_count: 684 warn(RelaxWarning("No interatomic data containers have been used in the calculation, skipping the RDC Q factor calculation.")) 685 return 686 if not rdc_data: 687 warn(RelaxWarning("No RDC data can be found for the alignment ID '%s', skipping the RDC Q factor calculation for this alignment." % align_id)) 688 continue 689 if not rdc_bc_data: 690 warn(RelaxWarning("No back-calculated RDC data can be found for the alignment ID '%s', skipping the RDC Q factor calculation for this alignment." % align_id)) 691 continue 692 693 # Normalisation factor of 2Da^2(4 + 3R)/5. 694 if norm2_flag: 695 D = dj * cdp.align_tensors[cdp.align_ids.index(align_id)].A_diag 696 Da = 1.0/3.0 * (D[2, 2] - (D[0, 0]+D[1, 1])/2.0) 697 Dr = 1.0/3.0 * (D[0, 0] - D[1, 1]) 698 if Da == 0: 699 R = nan 700 else: 701 R = Dr / Da 702 norm = 2.0 * (Da)**2 * (4.0 + 3.0*R**2)/5.0 703 if Da == 0.0: 704 norm = 1e-15 705 706 # The Q factor for the alignment. 707 cdp.q_factors_rdc[align_id] = sqrt(sse / N / norm) 708 709 else: 710 cdp.q_factors_rdc[align_id] = 0.0 711 712 # The second Q factor definition. 713 cdp.q_factors_rdc_norm2[align_id] = sqrt(sse / D2_sum) 714 715 # The total Q factor. 716 cdp.q_rdc = 0.0 717 cdp.q_rdc_norm2 = 0.0 718 for id in cdp.q_factors_rdc: 719 cdp.q_rdc = cdp.q_rdc + cdp.q_factors_rdc[id]**2 720 for id in cdp.q_factors_rdc_norm2: 721 cdp.q_rdc_norm2 = cdp.q_rdc_norm2 + cdp.q_factors_rdc_norm2[id]**2 722 cdp.q_rdc = sqrt(cdp.q_rdc / len(cdp.q_factors_rdc)) 723 cdp.q_rdc_norm2 = sqrt(cdp.q_rdc_norm2 / len(cdp.q_factors_rdc_norm2))

724 725

726 -def read(align_id=None, file=None, dir=None, file_data=None, data_type='D', spin_id1_col=None, spin_id2_col=None, data_col=None, error_col=None, sep=None, neg_g_corr=False, absolute=False):

727 """Read the RDC data from file. 728 729 @keyword align_id: The alignment tensor ID string. 730 @type align_id: str 731 @keyword file: The name of the file to open. 732 @type file: str 733 @keyword dir: The directory containing the file (defaults to the current directory if None). 734 @type dir: str or None 735 @keyword file_data: An alternative to 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. 736 @type file_data: list of lists 737 @keyword data_type: A string which is set to 'D' means that the splitting in the aligned sample was assumed to be J + D, or if set to '2D' then the splitting was taken as J + 2D. If set to 'T', then the data will be marked as being J+D values. 738 @keyword spin_id1_col: The column containing the spin ID strings of the first spin. 739 @type spin_id1_col: int 740 @keyword spin_id2_col: The column containing the spin ID strings of the second spin. 741 @type spin_id2_col: int 742 @keyword data_col: The column containing the RDC data in Hz. 743 @type data_col: int or None 744 @keyword error_col: The column containing the RDC errors. 745 @type error_col: int or None 746 @keyword sep: The column separator which, if None, defaults to whitespace. 747 @type sep: str or None 748 @keyword neg_g_corr: A flag which is used to correct for the negative gyromagnetic ratio of 15N. If True, a sign inversion will be applied to all RDC values to be loaded. 749 @type neg_g_corr: bool 750 @keyword absolute: A flag which if True indicates that the RDCs to load are signless. All RDCs will then be converted to positive values. 751 @type absolute: bool 752 """ 753 754 # Check the pipe setup. 755 check_pipe_setup(sequence=True) 756 757 # Either the data or error column must be supplied. 758 if data_col == None and error_col == None: 759 raise RelaxError("One of either the data or error column must be supplied.") 760 761 # Check the data types. 762 rdc_types = ['D', '2D', 'T'] 763 if data_type not in rdc_types: 764 raise RelaxError("The RDC data type '%s' must be one of %s." % (data_type, rdc_types)) 765 766 # Spin specific data. 767 ##################### 768 769 # Extract the data from the file, and remove comments and blank lines. 770 file_data = extract_data(file, dir, sep=sep) 771 file_data = strip(file_data, comments=True) 772 773 # Loop over the RDC data. 774 data = [] 775 for line in file_data: 776 # Invalid columns. 777 if spin_id1_col > len(line): 778 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line)) 779 continue 780 if spin_id2_col > len(line): 781 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line)) 782 continue 783 if data_col and data_col > len(line): 784 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line)) 785 continue 786 if error_col and error_col > len(line): 787 warn(RelaxWarning("The data %s is invalid, no error column can be found." % line)) 788 continue 789 790 # Unpack. 791 spin_id1 = line[spin_id1_col-1] 792 spin_id2 = line[spin_id2_col-1] 793 value = None 794 if data_col: 795 value = line[data_col-1] 796 error = None 797 if error_col: 798 error = line[error_col-1] 799 800 # Convert the spin IDs. 801 if spin_id1[0] in ["\"", "\'"]: 802 spin_id1 = eval(spin_id1) 803 if spin_id2[0] in ["\"", "\'"]: 804 spin_id2 = eval(spin_id2) 805 806 # Convert and check the value. 807 if value == 'None': 808 value = None 809 if value != None: 810 try: 811 value = float(value) 812 except ValueError: 813 warn(RelaxWarning("The RDC value of the line %s is invalid." % line)) 814 continue 815 816 # Convert and check the error. 817 if error == 'None': 818 error = None 819 if error != None: 820 try: 821 error = float(error) 822 except ValueError: 823 warn(RelaxWarning("The error value of the line %s is invalid." % line)) 824 continue 825 826 # Get the spins. 827 spin1 = return_spin(spin_id1) 828 spin2 = return_spin(spin_id2) 829 830 # Check the spin IDs. 831 if not spin1: 832 warn(RelaxWarning("The spin ID '%s' cannot be found in the current data pipe, skipping the data %s." % (spin_id1, line))) 833 continue 834 if not spin2: 835 warn(RelaxWarning("The spin ID '%s' cannot be found in the current data pipe, skipping the data %s." % (spin_id2, line))) 836 continue 837 838 # Get the interatomic data container. 839 interatom = return_interatom(spin_id1, spin_id2) 840 841 # Create the container if needed. 842 if interatom == None: 843 interatom = create_interatom(spin_id1=spin_id1, spin_id2=spin_id2) 844 845 # Test the error value (a value of 0.0 will cause the interatomic container to be deselected). 846 if error == 0.0: 847 interatom.select = False 848 warn(RelaxWarning("An error value of zero has been encountered, deselecting the interatomic container between spin '%s' and '%s'." % (spin_id1, spin_id2))) 849 continue 850 851 # Store the data type as global data (need for the conversion of RDC data). 852 if not hasattr(interatom, 'rdc_data_types'): 853 interatom.rdc_data_types = {} 854 if not align_id in interatom.rdc_data_types: 855 interatom.rdc_data_types[align_id] = data_type 856 857 # Convert and add the data. 858 if data_col: 859 # Data conversion. 860 value = convert(value, data_type, align_id, to_intern=True) 861 862 # Correction for the negative gyromagnetic ratio of 15N. 863 if neg_g_corr and value != None: 864 value = -value 865 866 # Absolute values. 867 if absolute: 868 # Force the value to be positive. 869 value = abs(value) 870 871 # Initialise. 872 if not hasattr(interatom, 'rdc'): 873 interatom.rdc = {} 874 875 # Add the value. 876 interatom.rdc[align_id] = value 877 878 # Store the absolute value flag. 879 if not hasattr(interatom, 'absolute_rdc'): 880 interatom.absolute_rdc = {} 881 interatom.absolute_rdc[align_id] = absolute 882 883 # Convert and add the error. 884 if error_col: 885 # Data conversion. 886 error = convert(error, data_type, align_id, to_intern=True) 887 888 # Initialise. 889 if not hasattr(interatom, 'rdc_err'): 890 interatom.rdc_err = {} 891 892 # Append the error. 893 interatom.rdc_err[align_id] = error 894 895 # Append the data for printout. 896 data.append([spin_id1, spin_id2]) 897 if is_float(value): 898 data[-1].append("%20.15f" % value) 899 else: 900 data[-1].append("%20s" % value) 901 if is_float(error): 902 data[-1].append("%20.15f" % error) 903 else: 904 data[-1].append("%20s" % error) 905 906 # No data, so fail hard! 907 if not len(data): 908 raise RelaxError("No RDC data could be extracted.") 909 910 # Print out. 911 print("The following RDCs have been loaded into the relax data store:\n") 912 write_data(out=sys.stdout, headings=["Spin_ID1", "Spin_ID2", "Value", "Error"], data=data) 913 914 # Initialise some global structures. 915 if not hasattr(cdp, 'align_ids'): 916 cdp.align_ids = [] 917 if not hasattr(cdp, 'rdc_ids'): 918 cdp.rdc_ids = [] 919 920 # Add the RDC id string. 921 if align_id not in cdp.align_ids: 922 cdp.align_ids.append(align_id) 923 if align_id not in cdp.rdc_ids: 924 cdp.rdc_ids.append(align_id)

925 926

927 -def return_rdc_data(sim_index=None):

928 """Set up the data structures for optimisation using RDCs as base data sets. 929 930 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 931 @type sim_index: None or int 932 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 933 - rdc, the RDC values. 934 - rdc_err, the RDC errors. 935 - rdc_weight, the RDC weights. 936 - vectors, the interatomic vectors (pseudo-atom dependent). 937 - rdc_const, the dipolar constants (pseudo-atom dependent). 938 - absolute, the absolute value flags (as 1's and 0's). 939 - T_flags, the flags for T = J+D type data (as 1's and 0's). 940 - j_couplings, the J coupling values if the RDC data type is set to T = J+D. 941 - pseudo_flags, the list of flags indicating if the interatomic data contains a pseudo-atom (as 1's and 0's). 942 @rtype: tuple of (numpy rank-2 float64 array, numpy rank-2 float64 array, numpy rank-2 float64 array, list of numpy rank-3 float64 arrays, list of lists of floats, numpy rank-2 int32 array, numpy rank-2 int32 array, numpy rank-2 float64 array, numpy rank-1 int32 array) 943 """ 944 945 # Sort out pseudo-atoms first. This only needs to be called once. 946 setup_pseudoatom_rdc() 947 948 # Initialise. 949 rdc = [] 950 rdc_err = [] 951 rdc_weight = [] 952 unit_vect = [] 953 rdc_const = [] 954 absolute = [] 955 T_flags = [] 956 j_couplings = [] 957 pseudo_flags = [] 958 959 # The unit vectors, RDC constants, and J couplings. 960 for interatom in interatomic_loop(): 961 # Get the spins. 962 spin1 = return_spin(interatom.spin_id1) 963 spin2 = return_spin(interatom.spin_id2) 964 965 # RDC checks. 966 if not check_rdcs(interatom): 967 continue 968 969 # Gyromagnetic ratios. 970 g1 = return_gyromagnetic_ratio(spin1.isotope) 971 g2 = return_gyromagnetic_ratio(spin2.isotope) 972 973 # Pseudo atoms. 974 if is_pseudoatom(spin1) and is_pseudoatom(spin2): 975 raise RelaxError("Support for both spins being in a dipole pair being pseudo-atoms is not implemented yet.") 976 if is_pseudoatom(spin1) or is_pseudoatom(spin2): 977 # Set the flag. 978 pseudo_flags.append(1) 979 980 # Alias the pseudo and normal atoms. 981 if is_pseudoatom(spin1): 982 pseudospin = spin1 983 base_spin = spin2 984 pseudospin_id = interatom.spin_id1 985 base_spin_id = interatom.spin_id2 986 else: 987 pseudospin = spin2 988 base_spin = spin1 989 pseudospin_id = interatom.spin_id2 990 base_spin_id = interatom.spin_id1 991 992 # Loop over the atoms of the pseudo-atom, storing the data. 993 pseudo_unit_vect = [] 994 pseudo_rdc_const = [] 995 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 996 # Get the corresponding interatomic data container. 997 pseudo_interatom = return_interatom(spin_id1=spin_id, spin_id2=base_spin_id) 998 999 # Check. 1000 if pseudo_interatom == None: 1001 raise RelaxError("The interatomic data container between the spins '%s' and '%s' for the pseudo-atom '%s' is not defined." % (base_spin_id, spin_id, pseudospin_id)) 1002 1003 # Add the vectors. 1004 if is_float(interatom.vector[0]): 1005 pseudo_unit_vect.append([pseudo_interatom.vector]) 1006 else: 1007 pseudo_unit_vect.append(pseudo_interatom.vector) 1008 1009 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 1010 pseudo_rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, pseudo_interatom.r)) 1011 1012 # Reorder the unit vectors so that the structure and pseudo-atom dimensions are swapped. 1013 pseudo_unit_vect = transpose(array(pseudo_unit_vect, float64), (1, 0, 2)) 1014 1015 # Block append the pseudo-data. 1016 unit_vect.append(pseudo_unit_vect) 1017 rdc_const.append(pseudo_rdc_const) 1018 1019 # Normal atom. 1020 else: 1021 # Set the flag. 1022 pseudo_flags.append(0) 1023 1024 # Add the vectors. 1025 if is_float(interatom.vector[0]): 1026 unit_vect.append([interatom.vector]) 1027 else: 1028 unit_vect.append(interatom.vector) 1029 1030 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 1031 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r)) 1032 1033 # Store the measured J coupling. 1034 if opt_uses_j_couplings(): 1035 j_couplings.append(interatom.j_coupling) 1036 1037 # Fix the unit vector data structure. 1038 num = None 1039 for rdc_index in range(len(unit_vect)): 1040 # Convert to numpy structures. 1041 unit_vect[rdc_index] = array(unit_vect[rdc_index], float64) 1042 1043 # Number of vectors. 1044 if num == None: 1045 if unit_vect[rdc_index] != None: 1046 num = len(unit_vect[rdc_index]) 1047 continue 1048 1049 # Check. 1050 if unit_vect[rdc_index] != None and len(unit_vect[rdc_index]) != num: 1051 raise RelaxError("The number of interatomic vectors for all no match:\n%s" % unit_vect[rdc_index]) 1052 1053 # Missing unit vectors. 1054 if num == None: 1055 raise RelaxError("No interatomic vectors could be found.") 1056 1057 # Update None entries. 1058 for i in range(len(unit_vect)): 1059 if unit_vect[i] == None: 1060 unit_vect[i] = [[None, None, None]]*num 1061 1062 # The RDC data. 1063 for i in range(len(cdp.align_ids)): 1064 # Alias the ID. 1065 align_id = cdp.align_ids[i] 1066 1067 # Skip non-optimised data. 1068 if not opt_uses_align_data(align_id): 1069 continue 1070 1071 # Append empty arrays to the RDC structures. 1072 rdc.append([]) 1073 rdc_err.append([]) 1074 rdc_weight.append([]) 1075 absolute.append([]) 1076 T_flags.append([]) 1077 1078 # Interatom loop. 1079 for interatom in interatomic_loop(): 1080 # Get the spins. 1081 spin1 = return_spin(interatom.spin_id1) 1082 spin2 = return_spin(interatom.spin_id2) 1083 1084 # RDC checks. 1085 if not check_rdcs(interatom): 1086 continue 1087 1088 # T-type data. 1089 if align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 1090 T_flags[-1].append(True) 1091 else: 1092 T_flags[-1].append(False) 1093 1094 # Check for J couplings if the RDC data type is T = J+D. 1095 if T_flags[-1][-1] and not hasattr(interatom, 'j_coupling'): 1096 continue 1097 1098 # Defaults of None. 1099 value = None 1100 error = None 1101 1102 # Normal set up. 1103 if align_id in interatom.rdc.keys(): 1104 # The RDC. 1105 if sim_index != None: 1106 value = interatom.rdc_sim[align_id][sim_index] 1107 else: 1108 value = interatom.rdc[align_id] 1109 1110 # The error. 1111 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 1112 # T values. 1113 if T_flags[-1][-1]: 1114 error = sqrt(interatom.rdc_err[align_id]**2 + interatom.j_coupling_err**2) 1115 1116 # D values. 1117 else: 1118 error = interatom.rdc_err[align_id] 1119 1120 # Append the RDCs to the list. 1121 rdc[-1].append(value) 1122 1123 # Append the RDC errors. 1124 rdc_err[-1].append(error) 1125 1126 # Append the weight. 1127 if hasattr(interatom, 'rdc_weight') and align_id in interatom.rdc_weight.keys(): 1128 rdc_weight[-1].append(interatom.rdc_weight[align_id]) 1129 else: 1130 rdc_weight[-1].append(1.0) 1131 1132 # Append the absolute value flag. 1133 if hasattr(interatom, 'absolute_rdc') and align_id in interatom.absolute_rdc.keys(): 1134 absolute[-1].append(interatom.absolute_rdc[align_id]) 1135 else: 1136 absolute[-1].append(False) 1137 1138 # Convert to numpy objects. 1139 rdc = array(rdc, float64) 1140 rdc_err = array(rdc_err, float64) 1141 rdc_weight = array(rdc_weight, float64) 1142 absolute = array(absolute, int32) 1143 T_flags = array(T_flags, int32) 1144 if not opt_uses_j_couplings(): 1145 j_couplings = None 1146 pseudo_flags = array(pseudo_flags, int32) 1147 1148 # Return the data structures. 1149 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const, absolute, T_flags, j_couplings, pseudo_flags

1150 1151

1152 -def set_errors(align_id=None, spin_id1=None, spin_id2=None, sd=None):

1153 """Set the RDC errors if not already present. 1154 1155 @keyword align_id: The optional alignment tensor ID string. 1156 @type align_id: str 1157 @keyword spin_id1: The optional spin ID string of the first spin. 1158 @type spin_id1: None or str 1159 @keyword spin_id2: The optional spin ID string of the second spin. 1160 @type spin_id2: None or str 1161 @keyword sd: The RDC standard deviation in Hz. 1162 @type sd: float or int. 1163 """ 1164 1165 # Check the pipe setup. 1166 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1167 1168 # Convert the align IDs to an array, or take all IDs. 1169 if align_id: 1170 align_ids = [align_id] 1171 else: 1172 align_ids = cdp.rdc_ids 1173 1174 # Loop over the interatomic data. 1175 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2): 1176 # No data structure. 1177 if not hasattr(interatom, 'rdc_err'): 1178 interatom.rdc_err = {} 1179 1180 # Set the error. 1181 for id in align_ids: 1182 interatom.rdc_err[id] = sd

1183 1184

1185 -def setup_pseudoatom_rdc():

1186 """Make sure that the interatom system is properly set up for pseudo-atoms and RDCs. 1187 1188 Interatomic data containers between the non-pseudo-atom and the pseudo-atom members will be deselected. 1189 """ 1190 1191 # Loop over all interatomic data containers. 1192 for interatom in interatomic_loop(): 1193 # Get the spins. 1194 spin1 = return_spin(interatom.spin_id1) 1195 spin2 = return_spin(interatom.spin_id2) 1196 1197 # Checks. 1198 flag1 = is_pseudoatom(spin1) 1199 flag2 = is_pseudoatom(spin2) 1200 1201 # No pseudo-atoms, so do nothing. 1202 if not (flag1 or flag2): 1203 continue 1204 1205 # Both are pseudo-atoms. 1206 if flag1 and flag2: 1207 warn(RelaxWarning("Support for both spins being in a dipole pair being pseudo-atoms is not implemented yet, deselecting the interatomic data container for the spin pair '%s' and '%s'." % (interatom.spin_id1, interatom.spin_id2))) 1208 interatom.select = False 1209 1210 # Alias the pseudo and normal atoms. 1211 pseudospin = spin1 1212 base_spin_id = interatom.spin_id2 1213 pseudospin_id = interatom.spin_id1 1214 if flag2: 1215 pseudospin = spin2 1216 base_spin_id = interatom.spin_id1 1217 pseudospin_id = interatom.spin_id2 1218 1219 # Loop over the atoms of the pseudo-atom. 1220 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 1221 # Get the corresponding interatomic data container. 1222 pseudo_interatom = return_interatom(spin_id1=spin_id, spin_id2=base_spin_id) 1223 1224 # Deselect if needed. 1225 if pseudo_interatom.select: 1226 warn(RelaxWarning("Deselecting the interatomic data container for the spin pair '%s' and '%s' as it is part of the pseudo-atom system of the spin pair '%s' and '%s'." % (pseudo_interatom.spin_id1, pseudo_interatom.spin_id2, base_spin_id, pseudospin_id))) 1227 pseudo_interatom.select = False

1228 1229

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

1231 """Set optimisation weights on the RDC data. 1232 1233 @keyword align_id: The alignment tensor ID string. 1234 @type align_id: str 1235 @keyword spin_id: The spin ID string. 1236 @type spin_id: None or str 1237 @keyword weight: The optimisation weight. The higher the value, the more importance the RDC will have. 1238 @type weight: float or int. 1239 """ 1240 1241 # Check the pipe setup. 1242 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1243 1244 # Loop over the interatomic data. 1245 for interatom in interatomic_loop(): 1246 # No data structure. 1247 if not hasattr(interatom, 'rdc_weight'): 1248 interatom.rdc_weight = {} 1249 1250 # Set the weight. 1251 interatom.rdc_weight[align_id] = weight

1252 1253

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

1255 """Display the RDC data corresponding to the alignment ID. 1256 1257 @keyword align_id: The alignment tensor ID string. 1258 @type align_id: str 1259 @keyword file: The file name or object to write to. 1260 @type file: str or file object 1261 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 1262 @type dir: str 1263 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 1264 @type bc: bool 1265 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1266 @type force: bool 1267 """ 1268 1269 # Check the pipe setup. 1270 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1271 1272 # Open the file for writing. 1273 file = open_write_file(file, dir, force) 1274 1275 # Loop over the interatomic data containers and collect the data. 1276 data = [] 1277 for interatom in interatomic_loop(): 1278 # Skip deselected containers. 1279 if not interatom.select: 1280 continue 1281 1282 # Skip containers with no RDCs. 1283 if not bc and (not hasattr(interatom, 'rdc') or align_id not in interatom.rdc.keys()): 1284 continue 1285 elif bc and (not hasattr(interatom, 'rdc_bc') or align_id not in interatom.rdc_bc.keys()): 1286 continue 1287 1288 # Append the spin data. 1289 data.append([]) 1290 data[-1].append(repr(interatom.spin_id1)) 1291 data[-1].append(repr(interatom.spin_id2)) 1292 1293 # Handle the missing rdc_data_types variable. 1294 data_type = None 1295 if hasattr(interatom, 'rdc_data_types'): 1296 data_type = interatom.rdc_data_types[align_id] 1297 1298 # The value. 1299 if bc: 1300 data[-1].append(repr(convert(interatom.rdc_bc[align_id], data_type, align_id))) 1301 else: 1302 data[-1].append(repr(convert(interatom.rdc[align_id], data_type, align_id))) 1303 1304 # The error. 1305 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 1306 data[-1].append(repr(convert(interatom.rdc_err[align_id], data_type, align_id))) 1307 else: 1308 data[-1].append(repr(None)) 1309 1310 # Write out. 1311 write_data(out=file, headings=["Spin_ID1", "Spin_ID2", "RDCs", "RDC_error"], data=data)

1312

Source Code for Module pipe_control.rdc