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, 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 file: The file name or object to write to. 384 @type file: str or file object 385 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 386 @type dir: str 387 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 388 @type force: bool 389 """ 390 391 # Check the pipe setup. 392 check_pipe_setup(sequence=True) 393 394 # Does RDC data exist? 395 if not hasattr(cdp, 'rdc_ids') or not cdp.rdc_ids: 396 warn(RelaxWarning("No RDC data exists, skipping file creation.")) 397 return 398 399 # Open the file for writing. 400 file = open_write_file(file, dir, force) 401 402 # Init. 403 data = [] 404 405 # The diagonal. 406 data.append([[-100, -100, 0], [100, 100, 0]]) 407 408 # Loop over the RDC data. 409 for align_id in cdp.rdc_ids: 410 # Append a new list for this alignment. 411 data.append([]) 412 413 # Errors present? 414 err_flag = False 415 for interatom in interatomic_loop(): 416 # Error present. 417 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 418 err_flag = True 419 break 420 421 # Loop over the interatomic data. 422 for interatom in interatomic_loop(): 423 # Skip if data is missing. 424 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(): 425 continue 426 427 # Append the data. 428 rdc_bc = convert(interatom.rdc_bc[align_id], interatom.rdc_data_types[align_id], align_id) 429 rdc = convert(interatom.rdc[align_id], interatom.rdc_data_types[align_id], align_id) 430 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 431 rdc_bc -= interatom.j_coupling 432 rdc -= interatom.j_coupling 433 data[-1].append([rdc_bc, rdc]) 434 435 # Errors. 436 if err_flag: 437 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 438 data[-1][-1].append(convert(interatom.rdc_err[align_id], interatom.rdc_data_types[align_id], align_id)) 439 else: 440 data[-1][-1].append(None) 441 442 # Label. 443 data[-1][-1].append("%s-%s" % (interatom.spin_id1, interatom.spin_id2)) 444 445 # The data size. 446 size = len(data) 447 448 # Only one data set. 449 data = [data] 450 451 # Graph type. 452 if err_flag: 453 graph_type = 'xydy' 454 else: 455 graph_type = 'xy' 456 457 # Grace file. 458 if format == 'grace': 459 # The header. 460 grace.write_xy_header(file=file, title="RDC correlation plot", sets=[size], set_names=[[None]+cdp.rdc_ids], linestyle=[[2]+[0]*size], data_type=['rdc_bc', 'rdc'], legend_pos=[[1, 0.5]]) 461 462 # The main data. 463 grace.write_xy_data(data=data, file=file, graph_type=graph_type)

464 465

466 -def delete(align_id=None):

467 """Delete the RDC data corresponding to the alignment ID. 468 469 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted. 470 @type align_id: str or None 471 """ 472 473 # Check the pipe setup. 474 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 475 476 # The IDs. 477 if not align_id: 478 ids = deepcopy(cdp.rdc_ids) 479 else: 480 ids = [align_id] 481 482 # Loop over the alignments, removing all the corresponding data. 483 for id in ids: 484 # The RDC ID. 485 cdp.rdc_ids.pop(cdp.rdc_ids.index(id)) 486 487 # The interatomic data. 488 for interatom in interatomic_loop(): 489 # The data. 490 if hasattr(interatom, 'rdc') and id in interatom.rdc: 491 interatom.rdc.pop(id) 492 493 # The error. 494 if hasattr(interatom, 'rdc_err') and id in interatom.rdc_err: 495 interatom.rdc_err.pop(id) 496 497 # The data type. 498 if hasattr(interatom, 'rdc_data_types') and id in interatom.rdc_data_types: 499 interatom.rdc_data_types.pop(id) 500 501 # Clean the global data. 502 if not hasattr(cdp, 'pcs_ids') or id not in cdp.pcs_ids: 503 cdp.align_ids.pop(cdp.align_ids.index(id))

504 505

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

507 """Display the RDC data corresponding to the alignment ID. 508 509 @keyword align_id: The alignment tensor ID string. 510 @type align_id: str 511 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed. 512 @type bc: bool 513 """ 514 515 # Check the pipe setup. 516 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 517 518 # Call the write method with sys.stdout as the file. 519 write(align_id=align_id, file=sys.stdout, bc=bc)

520 521

522 -def opt_uses_j_couplings():

523 """Determine of J couplings are needed for optimisation. 524 525 @return: True if J couplings are required, False otherwise. 526 @rtype: bool 527 """ 528 529 # Loop over the alignments. 530 for align_id in cdp.align_ids: 531 for interatom in interatomic_loop(): 532 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 533 return True 534 535 # No J values required. 536 return False

537 538

539 -def opt_uses_rdc(align_id):

540 """Determine if the RDC data for the given alignment ID is needed for optimisation. 541 542 @param align_id: The alignment ID string. 543 @type align_id: str 544 @return: True if the RDC data is to be used for optimisation, False otherwise. 545 @rtype: bool 546 """ 547 548 # No alignment IDs. 549 if not hasattr(cdp, 'rdc_ids'): 550 return False 551 552 # No RDC data for the alignment. 553 if align_id not in cdp.rdc_ids: 554 return False 555 556 # Is the tensor optimised? 557 tensor_flag = opt_uses_tensor(get_tensor_object(align_id)) 558 559 # Is the paramagnetic position optimised? 560 pos_flag = False 561 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 562 pos_flag = True 563 564 # Are the populations optimised? 565 prob_flag = False 566 if cdp.model == 'population': 567 prob_flag = True 568 569 # Not used. 570 if not tensor_flag and not pos_flag and not prob_flag: 571 return False 572 573 # The RDC data is to be used for optimisation. 574 return True

575 576

577 -def q_factors(spin_id=None):

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

704 705

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

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

905 906

907 -def return_rdc_data(sim_index=None):

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

1130 1131

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

1133 """Set the RDC errors if not already present. 1134 1135 @keyword align_id: The optional alignment tensor ID string. 1136 @type align_id: str 1137 @keyword spin_id1: The optional spin ID string of the first spin. 1138 @type spin_id1: None or str 1139 @keyword spin_id2: The optional spin ID string of the second spin. 1140 @type spin_id2: None or str 1141 @keyword sd: The RDC standard deviation in Hz. 1142 @type sd: float or int. 1143 """ 1144 1145 # Check the pipe setup. 1146 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1147 1148 # Convert the align IDs to an array, or take all IDs. 1149 if align_id: 1150 align_ids = [align_id] 1151 else: 1152 align_ids = cdp.rdc_ids 1153 1154 # Loop over the interatomic data. 1155 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2): 1156 # No data structure. 1157 if not hasattr(interatom, 'rdc_err'): 1158 interatom.rdc_err = {} 1159 1160 # Set the error. 1161 for id in align_ids: 1162 interatom.rdc_err[id] = sd

1163 1164

1165 -def setup_pseudoatom_rdc():

1166 """Make sure that the interatom system is properly set up for pseudo-atoms and RDCs. 1167 1168 Interatomic data containers between the non-pseudo-atom and the pseudo-atom members will be deselected. 1169 """ 1170 1171 # Loop over all interatomic data containers. 1172 for interatom in interatomic_loop(): 1173 # Get the spins. 1174 spin1 = return_spin(interatom.spin_id1) 1175 spin2 = return_spin(interatom.spin_id2) 1176 1177 # Checks. 1178 flag1 = is_pseudoatom(spin1) 1179 flag2 = is_pseudoatom(spin2) 1180 1181 # No pseudo-atoms, so do nothing. 1182 if not (flag1 or flag2): 1183 continue 1184 1185 # Both are pseudo-atoms. 1186 if flag1 and flag2: 1187 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))) 1188 interatom.select = False 1189 1190 # Alias the pseudo and normal atoms. 1191 pseudospin = spin1 1192 base_spin_id = interatom.spin_id2 1193 pseudospin_id = interatom.spin_id1 1194 if flag2: 1195 pseudospin = spin2 1196 base_spin_id = interatom.spin_id1 1197 pseudospin_id = interatom.spin_id2 1198 1199 # Loop over the atoms of the pseudo-atom. 1200 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 1201 # Get the corresponding interatomic data container. 1202 pseudo_interatom = return_interatom(spin_id1=spin_id, spin_id2=base_spin_id) 1203 1204 # Deselect if needed. 1205 if pseudo_interatom.select: 1206 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))) 1207 pseudo_interatom.select = False

1208 1209

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

1211 """Set optimisation weights on the RDC data. 1212 1213 @keyword align_id: The alignment tensor ID string. 1214 @type align_id: str 1215 @keyword spin_id: The spin ID string. 1216 @type spin_id: None or str 1217 @keyword weight: The optimisation weight. The higher the value, the more importance the RDC will have. 1218 @type weight: float or int. 1219 """ 1220 1221 # Check the pipe setup. 1222 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1223 1224 # Loop over the interatomic data. 1225 for interatom in interatomic_loop(): 1226 # No data structure. 1227 if not hasattr(interatom, 'rdc_weight'): 1228 interatom.rdc_weight = {} 1229 1230 # Set the weight. 1231 interatom.rdc_weight[align_id] = weight

1232 1233

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

1235 """Display the RDC data corresponding to the alignment ID. 1236 1237 @keyword align_id: The alignment tensor ID string. 1238 @type align_id: str 1239 @keyword file: The file name or object to write to. 1240 @type file: str or file object 1241 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 1242 @type dir: str 1243 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 1244 @type bc: bool 1245 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1246 @type force: bool 1247 """ 1248 1249 # Check the pipe setup. 1250 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1251 1252 # Open the file for writing. 1253 file = open_write_file(file, dir, force) 1254 1255 # Loop over the interatomic data containers and collect the data. 1256 data = [] 1257 for interatom in interatomic_loop(): 1258 # Skip deselected containers. 1259 if not interatom.select: 1260 continue 1261 1262 # Skip containers with no RDCs. 1263 if not bc and (not hasattr(interatom, 'rdc') or align_id not in interatom.rdc.keys()): 1264 continue 1265 elif bc and (not hasattr(interatom, 'rdc_bc') or align_id not in interatom.rdc_bc.keys()): 1266 continue 1267 1268 # Append the spin data. 1269 data.append([]) 1270 data[-1].append(repr(interatom.spin_id1)) 1271 data[-1].append(repr(interatom.spin_id2)) 1272 1273 # Handle the missing rdc_data_types variable. 1274 data_type = None 1275 if hasattr(interatom, 'rdc_data_types'): 1276 data_type = interatom.rdc_data_types[align_id] 1277 1278 # The value. 1279 if bc: 1280 data[-1].append(repr(convert(interatom.rdc_bc[align_id], data_type, align_id))) 1281 else: 1282 data[-1].append(repr(convert(interatom.rdc[align_id], data_type, align_id))) 1283 1284 # The error. 1285 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err.keys(): 1286 data[-1].append(repr(convert(interatom.rdc_err[align_id], data_type, align_id))) 1287 else: 1288 data[-1].append(repr(None)) 1289 1290 # Write out. 1291 write_data(out=file, headings=["Spin_ID1", "Spin_ID2", "RDCs", "RDC_error"], data=data)

1292

Source Code for Module pipe_control.rdc