pipe_control.rdc

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2015,2017 Edward d'Auvergne # 4 # Copyright (C) 2006 Chris MacRaild # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 20 # # 21 ############################################################################### 22 23 # Module docstring. 24 """Module for the manipulation of RDC data.""" 25 26 # Python module imports. 27 from copy import deepcopy 28 from math import ceil, floor, pi, sqrt 29 from numpy import array, int32, float64, ones, transpose, zeros 30 from numpy.linalg import norm 31 import sys 32 from warnings import warn 33 34 # relax module imports. 35 from lib.check_types import is_float 36 from lib.float import nan 37 from lib.alignment.rdc import ave_rdc_tensor 38 from lib.errors import RelaxError, RelaxNoAlignError, RelaxNoJError, RelaxNoRDCError, RelaxNoSequenceError, RelaxSpinTypeError 39 from lib.io import extract_data, open_write_file, strip, write_data 40 from lib.periodic_table import periodic_table 41 from lib.physical_constants import dipolar_constant 42 from lib.plotting.api import write_xy_data, write_xy_header 43 from lib.warnings import RelaxWarning, RelaxSpinTypeWarning 44 from pipe_control import pipes 45 from pipe_control.align_tensor import get_tensor_index, get_tensor_object, opt_uses_align_data, opt_uses_tensor 46 from pipe_control.interatomic import create_interatom, interatomic_loop, return_interatom 47 from pipe_control.mol_res_spin import exists_mol_res_spin_data, is_pseudoatom, pseudoatom_loop, return_spin 48 from pipe_control.pipes import check_pipe 49 50

51 -def back_calc(align_id=None):

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

143 144

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

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

195 196

197 -def check_rdcs(interatom):

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

275 276

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

278 """Convert the RDC based on the 'D' or '2D' data type. 279 280 @param value: The value or error to convert. 281 @type value: float or None 282 @param data_type: The RDC data type. Either 'D', '2D' or 'T'. 283 @type data_type: str 284 @param align_id: The alignment tensor ID string. 285 @type align_id: str 286 @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. 287 @type to_intern: bool 288 @return: The converted value. 289 @rtype: float or None 290 """ 291 292 # Handle values of None. 293 if value == None: 294 return None 295 296 # The conversion factor. 297 factor = 1.0 298 if data_type == '2D': 299 # Convert from 2D to D. 300 if to_intern: 301 factor = 2.0 302 303 # Convert from D to 2D. 304 else: 305 factor = 0.5 306 307 # Return the converted value. 308 return value * factor

309 310

311 -def copy(pipe_from=None, pipe_to=None, align_id=None, back_calc=True):

312 """Copy the RDC data from one data pipe to another. 313 314 @keyword pipe_from: The data pipe to copy the RDC data from. This defaults to the current data pipe. 315 @type pipe_from: str 316 @keyword pipe_to: The data pipe to copy the RDC data to. This defaults to the current data pipe. 317 @type pipe_to: str 318 @keyword align_id: The alignment ID string. 319 @type align_id: str 320 @keyword back_calc: A flag which if True will cause any back-calculated RDCs present to also be copied with the real values and errors. 321 @type back_calc: bool 322 """ 323 324 # Defaults. 325 if pipe_from == None and pipe_to == None: 326 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 327 elif pipe_from == None: 328 pipe_from = pipes.cdp_name() 329 elif pipe_to == None: 330 pipe_to = pipes.cdp_name() 331 332 # Check the pipe setup. 333 check_pipe_setup(pipe=pipe_from, rdc_id=align_id, sequence=True, rdc=True) 334 check_pipe_setup(pipe=pipe_to, sequence=True) 335 336 # Get the data pipes. 337 dp_from = pipes.get_pipe(pipe_from) 338 dp_to = pipes.get_pipe(pipe_to) 339 340 # The IDs. 341 if align_id == None: 342 align_ids = dp_from.align_ids 343 else: 344 align_ids = [align_id] 345 346 # Init target pipe global structures. 347 if not hasattr(dp_to, 'align_ids'): 348 dp_to.align_ids = [] 349 if not hasattr(dp_to, 'rdc_ids'): 350 dp_to.rdc_ids = [] 351 352 # Loop over the align IDs. 353 for align_id in align_ids: 354 # Printout. 355 print("\nCoping RDCs for the alignment ID '%s'." % align_id) 356 357 # Copy the global data. 358 if align_id not in dp_to.align_ids and align_id not in dp_to.align_ids: 359 dp_to.align_ids.append(align_id) 360 if align_id in dp_from.rdc_ids and align_id not in dp_to.rdc_ids: 361 dp_to.rdc_ids.append(align_id) 362 363 # Loop over the interatomic data of the source data pipe. 364 data = [] 365 for interatom_from in interatomic_loop(pipe=pipe_from): 366 # Find the matching interatomic data container in the target data pipe. 367 interatom_to = [] 368 for interatom in interatomic_loop(selection1=interatom_from.spin_id1, selection2=interatom_from.spin_id2, pipe=pipe_to, skip_desel=False): 369 interatom_to.append(interatom) 370 371 # No matching interatomic data container. 372 if interatom_to == []: 373 warn(RelaxWarning("The interatomic data container between the spins '%s' and '%s' cannot be found in the target data pipe." % (interatom_from.spin_id1, interatom_from.spin_id2))) 374 continue 375 376 # Too many containers. 377 elif len(interatom_to) != 1: 378 raise RelaxError("Too many interatomic data containers between the spins '%s' and '%s' exist in the target data pipe." % (interatom_from.spin_id1, interatom_from.spin_id2)) 379 380 # Collapse the container. 381 interatom_to = interatom_to[0] 382 383 # No data or errors. 384 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): 385 continue 386 387 # Initialise the data structures if necessary. 388 if hasattr(interatom_from, 'rdc') and not hasattr(interatom_to, 'rdc'): 389 interatom_to.rdc = {} 390 if back_calc and hasattr(interatom_from, 'rdc_bc') and not hasattr(interatom_to, 'rdc_bc'): 391 interatom_to.rdc_bc = {} 392 if hasattr(interatom_from, 'rdc_err') and not hasattr(interatom_to, 'rdc_err'): 393 interatom_to.rdc_err = {} 394 if hasattr(interatom_from, 'rdc_data_types') and not hasattr(interatom_to, 'rdc_data_types'): 395 interatom_to.rdc_data_types = {} 396 if hasattr(interatom_from, 'absolute_rdc') and not hasattr(interatom_to, 'absolute_rdc'): 397 interatom_to.absolute_rdc = {} 398 399 # Copy the value and error from pipe_from. 400 value = None 401 error = None 402 value_bc = None 403 data_type = None 404 absolute_rdc = None 405 if hasattr(interatom_from, 'rdc'): 406 value = interatom_from.rdc[align_id] 407 interatom_to.rdc[align_id] = value 408 if back_calc and hasattr(interatom_from, 'rdc_bc'): 409 value_bc = interatom_from.rdc_bc[align_id] 410 interatom_to.rdc_bc[align_id] = value_bc 411 if hasattr(interatom_from, 'rdc_err'): 412 error = interatom_from.rdc_err[align_id] 413 interatom_to.rdc_err[align_id] = error 414 if hasattr(interatom_from, 'rdc_data_types'): 415 data_type = interatom_from.rdc_data_types[align_id] 416 interatom_to.rdc_data_types[align_id] = data_type 417 if hasattr(interatom_from, 'absolute_rdc'): 418 absolute_rdc = interatom_from.absolute_rdc[align_id] 419 interatom_to.absolute_rdc[align_id] = absolute_rdc 420 421 # Append the data for printout. 422 data.append([interatom_from.spin_id1, interatom_from.spin_id2]) 423 if is_float(value): 424 data[-1].append("%20.15f" % value) 425 else: 426 data[-1].append("%20s" % value) 427 if back_calc: 428 if is_float(value_bc): 429 data[-1].append("%20.15f" % value_bc) 430 else: 431 data[-1].append("%20s" % value_bc) 432 if is_float(error): 433 data[-1].append("%20.15f" % error) 434 else: 435 data[-1].append("%20s" % error) 436 data[-1].append("%20s" % data_type) 437 data[-1].append("%20s" % absolute_rdc) 438 439 # Printout. 440 print("The following RDCs have been copied:\n") 441 if back_calc: 442 write_data(out=sys.stdout, headings=["Spin_ID1", "Spin_ID2", "Value", "Back-calculated", "Error", "Data_type", "Absolute_RDC"], data=data) 443 else: 444 write_data(out=sys.stdout, headings=["Spin_ID1", "Spin_ID2", "Value", "Error", "Data_type", "Absolute_RDC"], data=data)

445 446

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

448 """Generate a correlation plot of the measured vs. back-calculated RDCs. 449 450 @keyword format: The format for the plot file. The following values are accepted: 'grace', a Grace plot; None, a plain text file. 451 @type format: str or None 452 @keyword title: The title for the plot, overriding the default. 453 @type title: None or str 454 @keyword subtitle: The subtitle for the plot, overriding the default. 455 @type subtitle: None or str 456 @keyword file: The file name or object to write to. 457 @type file: str or file object 458 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 459 @type dir: str 460 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 461 @type force: bool 462 """ 463 464 # Check the pipe setup. 465 check_pipe_setup(sequence=True) 466 467 # Does RDC data exist? 468 if not hasattr(cdp, 'rdc_ids') or not cdp.rdc_ids: 469 warn(RelaxWarning("No RDC data exists, skipping file creation.")) 470 return 471 472 # Open the file for writing. 473 file = open_write_file(file, dir, force) 474 475 # Init. 476 data = [] 477 orig_title = title 478 if orig_title == None: 479 title = "RDC correlation plot" 480 axis_labels = ["Back-calculated RDC (Hz)", "Measured RDC (Hz)"] 481 482 # The diagonal. 483 data.append([[-100, -100, 0], [100, 100, 0]]) 484 485 # Loop over the RDC data. 486 min_rdc = 1e100 487 max_rdc = -1e100 488 for align_id in cdp.rdc_ids: 489 # Append a new list for this alignment. 490 data.append([]) 491 492 # Errors present? 493 err_flag = False 494 for interatom in interatomic_loop(): 495 # Error present. 496 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 497 err_flag = True 498 break 499 500 # Loop over the interatomic data. 501 for interatom in interatomic_loop(skip_desel=True): 502 # Skip if data is missing. 503 if not hasattr(interatom, 'rdc') or not hasattr(interatom, 'rdc_bc') or not align_id in interatom.rdc or not align_id in interatom.rdc_bc: 504 continue 505 if interatom.rdc[align_id] == None or interatom.rdc_bc[align_id] == None: 506 continue 507 508 # Convert between the 2D and D notation. 509 rdc_bc = convert(interatom.rdc_bc[align_id], interatom.rdc_data_types[align_id], align_id) 510 rdc = convert(interatom.rdc[align_id], interatom.rdc_data_types[align_id], align_id) 511 512 # T=J+D type data. 513 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 514 # Convert T=J+D type data to D, if not absolute. 515 if not interatom.absolute_rdc[align_id]: 516 rdc_bc -= interatom.j_coupling 517 rdc -= interatom.j_coupling 518 519 # A different title and axis labels. 520 else: 521 if orig_title == None: 522 title = "T = J+D correlation plot" 523 axis_labels = ["Measured T = J+D (Hz)", "Back-calculated T = J+D (Hz)"] 524 525 # Append the data. 526 data[-1].append([rdc_bc, rdc]) 527 528 # The minimum of all data sets. 529 if rdc < min_rdc: 530 min_rdc = rdc 531 if rdc_bc < min_rdc: 532 min_rdc = rdc_bc 533 534 # The maximum of all data sets. 535 if rdc > max_rdc: 536 max_rdc = rdc 537 if rdc_bc > max_rdc: 538 max_rdc = rdc_bc 539 540 # Errors. 541 if err_flag: 542 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 543 data[-1][-1].append(convert(interatom.rdc_err[align_id], interatom.rdc_data_types[align_id], align_id)) 544 else: 545 data[-1][-1].append(None) 546 547 # Label. 548 data[-1][-1].append("%s-%s" % (interatom.spin_id1, interatom.spin_id2)) 549 550 # The data size. 551 size = len(data) 552 553 # Round the data limits. 554 max_rdc = float(ceil(max_rdc)) 555 min_rdc = float(floor(min_rdc)) 556 557 # Only one data set. 558 data = [data] 559 560 # Graph type. 561 if err_flag: 562 graph_type = 'xydy' 563 else: 564 graph_type = 'xy' 565 566 # The header. 567 write_xy_header(format=format, file=file, title=title, subtitle=subtitle, world=[[min_rdc, min_rdc, max_rdc, max_rdc]], 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]]) 568 569 # The main data. 570 write_xy_data(format=format, data=data, file=file, graph_type=graph_type, autoscale=False)

571 572

573 -def delete(align_id=None):

574 """Delete the RDC data corresponding to the alignment ID. 575 576 @keyword align_id: The alignment tensor ID string. If not specified, all data will be deleted. 577 @type align_id: str or None 578 """ 579 580 # Check the pipe setup. 581 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 582 583 # The IDs. 584 if not align_id: 585 ids = deepcopy(cdp.rdc_ids) 586 else: 587 ids = [align_id] 588 589 # Loop over the alignments, removing all the corresponding data. 590 for id in ids: 591 # The RDC ID. 592 cdp.rdc_ids.pop(cdp.rdc_ids.index(id)) 593 594 # The interatomic data. 595 for interatom in interatomic_loop(): 596 # The data. 597 if hasattr(interatom, 'rdc') and id in interatom.rdc: 598 interatom.rdc.pop(id) 599 600 # The error. 601 if hasattr(interatom, 'rdc_err') and id in interatom.rdc_err: 602 interatom.rdc_err.pop(id) 603 604 # The data type. 605 if hasattr(interatom, 'rdc_data_types') and id in interatom.rdc_data_types: 606 interatom.rdc_data_types.pop(id) 607 608 # Clean the global data. 609 if not hasattr(cdp, 'pcs_ids') or id not in cdp.pcs_ids: 610 cdp.align_ids.pop(cdp.align_ids.index(id))

611 612

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

614 """Display the RDC data corresponding to the alignment ID. 615 616 @keyword align_id: The alignment tensor ID string. 617 @type align_id: str 618 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be displayed. 619 @type bc: bool 620 """ 621 622 # Check the pipe setup. 623 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 624 625 # Call the write method with sys.stdout as the file. 626 write(align_id=align_id, file=sys.stdout, bc=bc)

627 628

629 -def opt_uses_j_couplings():

630 """Determine of J couplings are needed for optimisation. 631 632 @return: True if J couplings are required, False otherwise. 633 @rtype: bool 634 """ 635 636 # Loop over the alignments. 637 for align_id in cdp.align_ids: 638 for interatom in interatomic_loop(): 639 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 640 return True 641 642 # No J values required. 643 return False

644 645

646 -def opt_uses_rdc(align_id):

647 """Determine if the RDC data for the given alignment ID is needed for optimisation. 648 649 @param align_id: The alignment ID string. 650 @type align_id: str 651 @return: True if the RDC data is to be used for optimisation, False otherwise. 652 @rtype: bool 653 """ 654 655 # No alignment IDs. 656 if not hasattr(cdp, 'rdc_ids'): 657 return False 658 659 # No RDC data for the alignment. 660 if align_id not in cdp.rdc_ids: 661 return False 662 663 # Is the tensor optimised? 664 tensor_flag = opt_uses_tensor(get_tensor_object(align_id)) 665 666 # Is the paramagnetic position optimised? 667 pos_flag = False 668 if hasattr(cdp, 'paramag_centre_fixed') and not cdp.paramag_centre_fixed: 669 pos_flag = True 670 671 # Are the populations optimised? 672 prob_flag = False 673 if cdp.model == 'population': 674 prob_flag = True 675 676 # Not used. 677 if not tensor_flag and not pos_flag and not prob_flag: 678 return False 679 680 # The RDC data is to be used for optimisation. 681 return True

682 683

684 -def q_factors(spin_id=None, sim_index=None, verbosity=1):

685 """Calculate the Q factors for the RDC data. 686 687 @keyword spin_id: The spin ID string used to restrict the Q factor calculation to a subset of all spins. 688 @type spin_id: None or str 689 @keyword sim_index: The optional Monte Carlo simulation index. 690 @type sim_index: None or int 691 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 692 @type verbosity: int 693 """ 694 695 # Check the pipe setup. 696 check_pipe_setup(sequence=True) 697 698 # Simulation flag. 699 sim_flag = (sim_index != None) 700 701 # No RDCs, so no Q factors can be calculated. 702 if not hasattr(cdp, 'rdc_ids') or not len(cdp.rdc_ids): 703 warn(RelaxWarning("No RDC data exists, Q factors cannot be calculated.")) 704 return 705 706 # Q factor dictonaries. 707 if not sim_flag: 708 cdp.q_factors_rdc_norm_tensor_size = {} 709 cdp.q_factors_rdc_norm_squared_sum = {} 710 else: 711 if not hasattr(cdp, 'q_factors_rdc_norm_tensor_size_sim'): 712 cdp.q_factors_rdc_norm_tensor_size_sim = {} 713 if not hasattr(cdp, 'q_factors_rdc_norm_squared_sum_sim'): 714 cdp.q_factors_rdc_norm_squared_sum_sim = {} 715 716 # Loop over the alignments. 717 for align_id in cdp.rdc_ids: 718 # Init. 719 D2_sum = 0.0 720 sse = 0.0 721 if sim_flag and align_id not in cdp.q_factors_rdc_norm_tensor_size_sim: 722 cdp.q_factors_rdc_norm_tensor_size_sim[align_id] = [None] * cdp.sim_number 723 if sim_flag and align_id not in cdp.q_factors_rdc_norm_squared_sum_sim: 724 cdp.q_factors_rdc_norm_squared_sum_sim[align_id] = [None] * cdp.sim_number 725 726 # Interatomic data loop. 727 dj = None 728 N = 0 729 interatom_count = 0 730 rdc_data = False 731 rdc_bc_data = False 732 norm2_flag = True 733 for interatom in interatomic_loop(): 734 rdc_data_align_id = False 735 rdc_bc_data_align_id = False 736 737 # Increment the counter. 738 interatom_count += 1 739 740 # Data checks. 741 if not sim_flag: 742 if hasattr(interatom, 'rdc') and align_id in interatom.rdc and interatom.rdc[align_id] != None: 743 rdc_data_align_id = True 744 if hasattr(interatom, 'rdc_bc') and align_id in interatom.rdc_bc and interatom.rdc_bc[align_id] != None: 745 rdc_bc_data_align_id = True 746 else: 747 if hasattr(interatom, 'rdc_sim') and align_id in interatom.rdc_sim and interatom.rdc_sim[align_id][sim_index] != None: 748 rdc_data_align_id = True 749 if hasattr(interatom, 'rdc_sim_bc') and align_id in interatom.rdc_sim_bc and interatom.rdc_sim_bc[align_id][sim_index] != None: 750 rdc_bc_data_align_id = True 751 j_flag = False 752 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 753 j_flag = True 754 if not hasattr(interatom, 'j_coupling'): 755 raise RelaxNoJError 756 757 # The global flags. 758 rdc_data = rdc_data or rdc_data_align_id 759 rdc_bc_data = rdc_bc_data or rdc_bc_data_align_id 760 761 # Skip containers without RDC data. 762 if not rdc_data_align_id or not rdc_bc_data_align_id: 763 continue 764 765 # Get the spins. 766 spin1 = return_spin(spin_hash=interatom._spin_hash1) 767 spin2 = return_spin(spin_hash=interatom._spin_hash2) 768 769 # Alias the RDC data. 770 if not sim_flag: 771 rdc = interatom.rdc[align_id] 772 rdc_bc = interatom.rdc_bc[align_id] 773 else: 774 rdc = interatom.rdc_sim[align_id][sim_index] 775 rdc_bc = interatom.rdc_sim_bc[align_id][sim_index] 776 777 # Sum of squares. 778 sse = sse + (rdc - rdc_bc)**2 779 780 # Sum the RDCs squared (for one type of normalisation). 781 if j_flag: 782 D2_sum = D2_sum + (rdc - interatom.j_coupling)**2 783 else: 784 D2_sum = D2_sum + rdc**2 785 786 # Skip the 2Da^2(4 + 3R)/5 normalised Q factor if no tensor is present. 787 if norm2_flag and not hasattr(cdp, 'align_tensors'): 788 if not sim_flag: 789 warn(RelaxWarning("No alignment tensors are present for the alignment '%s', skipping the Q factor normalised with 2Da^2(4 + 3R)/5." % align_id)) 790 norm2_flag = False 791 792 # Skip the 2Da^2(4 + 3R)/5 normalised Q factor if pseudo-atoms are present. 793 if norm2_flag and (is_pseudoatom(spin1) or is_pseudoatom(spin2)): 794 if not sim_flag: 795 warn(RelaxWarning("Pseudo-atoms are present for the alignment '%s', skipping the Q factor normalised with 2Da^2(4 + 3R)/5." % align_id)) 796 norm2_flag = False 797 798 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 799 if norm2_flag: 800 # Data checks. 801 if not hasattr(spin1, 'isotope'): 802 raise RelaxSpinTypeError(spin_id=interatom.spin_id1) 803 if not hasattr(spin2, 'isotope'): 804 raise RelaxSpinTypeError(spin_id=interatom.spin_id2) 805 806 # Gyromagnetic ratios. 807 g1 = periodic_table.gyromagnetic_ratio(spin1.isotope) 808 g2 = periodic_table.gyromagnetic_ratio(spin2.isotope) 809 810 # Calculate the dipolar constant. 811 dj_new = 3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r) 812 if dj != None and dj_new != dj: 813 if not sim_flag: 814 warn(RelaxWarning("The dipolar constant is not the same for all RDCs for the alignment '%s', skipping the Q factor normalised with 2Da^2(4 + 3R)/5." % align_id)) 815 norm2_flag = False 816 else: 817 dj = dj_new 818 819 # Increment the number of data sets. 820 N = N + 1 821 822 # Warnings (and then exit). 823 if not interatom_count: 824 warn(RelaxWarning("No interatomic data containers have been used in the calculation, skipping the RDC Q factor calculation.")) 825 return 826 if not rdc_data: 827 warn(RelaxWarning("No RDC data can be found for the alignment ID '%s', skipping the RDC Q factor calculation for this alignment." % align_id)) 828 continue 829 if not rdc_bc_data: 830 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)) 831 continue 832 833 # Set up the sim Q factor structure if needed. 834 if sim_flag and not align_id in cdp.q_factors_rdc_norm_tensor_size_sim: 835 cdp.q_factors_rdc_norm_tensor_size_sim[align_id] = [None] * cdp.sim_number 836 837 # Normalisation factor of 2Da^2(4 + 3R)/5. 838 if norm2_flag: 839 A = cdp.align_tensors[cdp.align_ids.index(align_id)] 840 if not sim_flag: 841 D = dj * A.A_diag 842 else: 843 D = dj * A.A_diag_sim[sim_index] 844 Da = 1.0/3.0 * (D[2, 2] - (D[0, 0]+D[1, 1])/2.0) 845 Dr = 1.0/3.0 * (D[0, 0] - D[1, 1]) 846 if Da == 0: 847 R = nan 848 else: 849 R = Dr / Da 850 norm = 2.0 * (Da)**2 * (4.0 + 3.0*R**2)/5.0 851 if Da == 0.0: 852 norm = 1e-15 853 854 # The Q factor for the alignment. 855 if sim_flag: 856 cdp.q_factors_rdc_norm_tensor_size_sim[align_id][sim_index] = sqrt(sse / N / norm) 857 else: 858 cdp.q_factors_rdc_norm_tensor_size[align_id] = sqrt(sse / N / norm) 859 860 else: 861 if sim_flag: 862 cdp.q_factors_rdc_norm_tensor_size_sim[align_id][sim_index] = 0.0 863 else: 864 cdp.q_factors_rdc_norm_tensor_size[align_id] = 0.0 865 866 # The second Q factor definition. 867 if sim_flag: 868 cdp.q_factors_rdc_norm_squared_sum_sim[align_id][sim_index] = sqrt(sse / D2_sum) 869 else: 870 cdp.q_factors_rdc_norm_squared_sum[align_id] = sqrt(sse / D2_sum) 871 872 # ID and RDC Q factor printout. 873 if verbosity: 874 print("\nRDC Q factors normalised by the tensor size (2Da^2(4 + 3R)/5):") 875 for align_id in cdp.rdc_ids: 876 if align_id in cdp.q_factors_rdc_norm_tensor_size: 877 if sim_flag: 878 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_rdc_norm_tensor_size_sim[align_id][sim_index])) 879 else: 880 print(" Alignment ID '%s': %.3f" % (align_id, cdp.q_factors_rdc_norm_tensor_size[align_id])) 881 print("\nRDC Q factors normalised by the sum of RDCs squared:") 882 for align_id in cdp.rdc_ids: 883 if align_id in cdp.q_factors_rdc_norm_squared_sum: 884 if sim_flag: 885 print(" Alignment ID '%s': %.13f" % (align_id, cdp.q_factors_rdc_norm_squared_sum_sim[align_id][sim_index])) 886 else: 887 print(" Alignment ID '%s': %.13f" % (align_id, cdp.q_factors_rdc_norm_squared_sum[align_id])) 888 889 # The total Q factor. 890 if sim_flag: 891 if not hasattr(cdp, 'q_rdc_norm_tensor_size_sim'): 892 cdp.q_rdc_norm_tensor_size_sim = [None] * cdp.sim_number 893 if not hasattr(cdp, 'q_rdc_norm_squared_sum_sim'): 894 cdp.q_rdc_norm_squared_sum_sim = [None] * cdp.sim_number 895 cdp.q_rdc_norm_tensor_size_sim[sim_index] = 0.0 896 cdp.q_rdc_norm_squared_sum_sim[sim_index] = 0.0 897 for id in cdp.q_factors_rdc_norm_tensor_size_sim: 898 cdp.q_rdc_norm_tensor_size_sim[sim_index] = cdp.q_rdc_norm_tensor_size_sim[sim_index] + cdp.q_factors_rdc_norm_tensor_size_sim[id][sim_index]**2 899 for id in cdp.q_factors_rdc_norm_squared_sum_sim: 900 cdp.q_rdc_norm_squared_sum_sim[sim_index] = cdp.q_rdc_norm_squared_sum_sim[sim_index] + cdp.q_factors_rdc_norm_squared_sum_sim[id][sim_index]**2 901 cdp.q_rdc_norm_tensor_size_sim[sim_index] = sqrt(cdp.q_rdc_norm_tensor_size_sim[sim_index] / len(cdp.q_factors_rdc_norm_tensor_size_sim)) 902 cdp.q_rdc_norm_squared_sum_sim[sim_index] = sqrt(cdp.q_rdc_norm_squared_sum_sim[sim_index] / len(cdp.q_factors_rdc_norm_squared_sum_sim)) 903 else: 904 cdp.q_rdc_norm_tensor_size = 0.0 905 cdp.q_rdc_norm_squared_sum = 0.0 906 for id in cdp.q_factors_rdc_norm_tensor_size: 907 cdp.q_rdc_norm_tensor_size = cdp.q_rdc_norm_tensor_size + cdp.q_factors_rdc_norm_tensor_size[id]**2 908 for id in cdp.q_factors_rdc_norm_squared_sum: 909 cdp.q_rdc_norm_squared_sum = cdp.q_rdc_norm_squared_sum + cdp.q_factors_rdc_norm_squared_sum[id]**2 910 cdp.q_rdc_norm_tensor_size = sqrt(cdp.q_rdc_norm_tensor_size / len(cdp.q_factors_rdc_norm_tensor_size)) 911 cdp.q_rdc_norm_squared_sum = sqrt(cdp.q_rdc_norm_squared_sum / len(cdp.q_factors_rdc_norm_squared_sum))

912 913

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

915 """Read the RDC data from file. 916 917 @keyword align_id: The alignment tensor ID string. 918 @type align_id: str 919 @keyword file: The name of the file to open. 920 @type file: str 921 @keyword dir: The directory containing the file (defaults to the current directory if None). 922 @type dir: str or None 923 @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. 924 @type file_data: list of lists 925 @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. 926 @keyword spin_id1_col: The column containing the spin ID strings of the first spin. 927 @type spin_id1_col: int 928 @keyword spin_id2_col: The column containing the spin ID strings of the second spin. 929 @type spin_id2_col: int 930 @keyword data_col: The column containing the RDC data in Hz. 931 @type data_col: int or None 932 @keyword error_col: The column containing the RDC errors. 933 @type error_col: int or None 934 @keyword sep: The column separator which, if None, defaults to whitespace. 935 @type sep: str or None 936 @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. 937 @type neg_g_corr: bool 938 @keyword absolute: A flag which if True indicates that the RDCs to load are signless. All RDCs will then be converted to positive values. 939 @type absolute: bool 940 """ 941 942 # Check the pipe setup. 943 check_pipe_setup(sequence=True) 944 945 # Either the data or error column must be supplied. 946 if data_col == None and error_col == None: 947 raise RelaxError("One of either the data or error column must be supplied.") 948 949 # Check the data types. 950 rdc_types = ['D', '2D', 'T'] 951 if data_type not in rdc_types: 952 raise RelaxError("The RDC data type '%s' must be one of %s." % (data_type, rdc_types)) 953 954 # Spin specific data. 955 ##################### 956 957 # Extract the data from the file, and remove comments and blank lines. 958 file_data = extract_data(file, dir, sep=sep) 959 file_data = strip(file_data, comments=True) 960 961 # Loop over the RDC data. 962 data = [] 963 for line in file_data: 964 # Invalid columns. 965 if spin_id1_col > len(line): 966 warn(RelaxWarning("The data %s is invalid, no first spin ID column can be found." % line)) 967 continue 968 if spin_id2_col > len(line): 969 warn(RelaxWarning("The data %s is invalid, no second spin ID column can be found." % line)) 970 continue 971 if data_col and data_col > len(line): 972 warn(RelaxWarning("The data %s is invalid, no data column can be found." % line)) 973 continue 974 if error_col and error_col > len(line): 975 warn(RelaxWarning("The data %s is invalid, no error column can be found." % line)) 976 continue 977 978 # Unpack. 979 spin_id1 = line[spin_id1_col-1] 980 spin_id2 = line[spin_id2_col-1] 981 value = None 982 if data_col: 983 value = line[data_col-1] 984 error = None 985 if error_col: 986 error = line[error_col-1] 987 988 # Convert the spin IDs. 989 if spin_id1[0] in ["\"", "\'"]: 990 spin_id1 = eval(spin_id1) 991 if spin_id2[0] in ["\"", "\'"]: 992 spin_id2 = eval(spin_id2) 993 994 # Convert and check the value. 995 if value == 'None': 996 value = None 997 if value != None: 998 try: 999 value = float(value) 1000 except ValueError: 1001 warn(RelaxWarning("The RDC value of the line %s is invalid." % line)) 1002 continue 1003 1004 # Convert and check the error. 1005 if error == 'None': 1006 error = None 1007 if error != None: 1008 try: 1009 error = float(error) 1010 except ValueError: 1011 warn(RelaxWarning("The error value of the line %s is invalid." % line)) 1012 continue 1013 1014 # Get the spins. 1015 spin1 = return_spin(spin_id=spin_id1) 1016 spin2 = return_spin(spin_id=spin_id2) 1017 1018 # Check the spin IDs. 1019 if not spin1: 1020 warn(RelaxWarning("The spin ID '%s' cannot be found in the current data pipe, skipping the data %s." % (spin_id1, line))) 1021 continue 1022 if not spin2: 1023 warn(RelaxWarning("The spin ID '%s' cannot be found in the current data pipe, skipping the data %s." % (spin_id2, line))) 1024 continue 1025 1026 # Get the interatomic data container. 1027 interatom = return_interatom(spin_hash1=spin1._hash, spin_hash2=spin2._hash) 1028 1029 # Create the container if needed. 1030 if interatom == None: 1031 interatom = create_interatom(spin_id1=spin_id1, spin_id2=spin_id2) 1032 1033 # Test the error value (a value of 0.0 will cause the interatomic container to be deselected). 1034 if error == 0.0: 1035 interatom.select = False 1036 warn(RelaxWarning("An error value of zero has been encountered, deselecting the interatomic container between spin '%s' and '%s'." % (spin_id1, spin_id2))) 1037 continue 1038 1039 # Store the data type as global data (need for the conversion of RDC data). 1040 if not hasattr(interatom, 'rdc_data_types'): 1041 interatom.rdc_data_types = {} 1042 if not align_id in interatom.rdc_data_types: 1043 interatom.rdc_data_types[align_id] = data_type 1044 1045 # Convert and add the data. 1046 if data_col: 1047 # Data conversion. 1048 value = convert(value, data_type, align_id, to_intern=True) 1049 1050 # Correction for the negative gyromagnetic ratio of 15N. 1051 if neg_g_corr and value != None: 1052 value = -value 1053 1054 # Absolute values. 1055 if absolute: 1056 # Force the value to be positive. 1057 value = abs(value) 1058 1059 # Initialise. 1060 if not hasattr(interatom, 'rdc'): 1061 interatom.rdc = {} 1062 1063 # Add the value. 1064 interatom.rdc[align_id] = value 1065 1066 # Store the absolute value flag. 1067 if not hasattr(interatom, 'absolute_rdc'): 1068 interatom.absolute_rdc = {} 1069 interatom.absolute_rdc[align_id] = absolute 1070 1071 # Convert and add the error. 1072 if error_col: 1073 # Data conversion. 1074 error = convert(error, data_type, align_id, to_intern=True) 1075 1076 # Initialise. 1077 if not hasattr(interatom, 'rdc_err'): 1078 interatom.rdc_err = {} 1079 1080 # Append the error. 1081 interatom.rdc_err[align_id] = error 1082 1083 # Append the data for printout. 1084 data.append([interatom.spin_id1, interatom.spin_id2]) 1085 if is_float(value): 1086 data[-1].append("%20.15f" % value) 1087 else: 1088 data[-1].append("%20s" % value) 1089 if is_float(error): 1090 data[-1].append("%20.15f" % error) 1091 else: 1092 data[-1].append("%20s" % error) 1093 1094 # No data, so fail hard! 1095 if not len(data): 1096 raise RelaxError("No RDC data could be extracted.") 1097 1098 # Print out. 1099 print("The following RDCs have been loaded into the relax data store:\n") 1100 write_data(out=sys.stdout, headings=["Spin_ID1", "Spin_ID2", "Value", "Error"], data=data) 1101 1102 # Initialise some global structures. 1103 if not hasattr(cdp, 'align_ids'): 1104 cdp.align_ids = [] 1105 if not hasattr(cdp, 'rdc_ids'): 1106 cdp.rdc_ids = [] 1107 1108 # Add the RDC id string. 1109 if align_id not in cdp.align_ids: 1110 cdp.align_ids.append(align_id) 1111 if align_id not in cdp.rdc_ids: 1112 cdp.rdc_ids.append(align_id)

1113 1114

1115 -def return_rdc_data(sim_index=None, verbosity=0):

1116 """Set up the data structures for optimisation using RDCs as base data sets. 1117 1118 @keyword sim_index: The index of the simulation to optimise. This should be None if normal optimisation is desired. 1119 @type sim_index: None or int 1120 @keyword verbosity: A flag specifying the amount of information to print. The higher the value, the greater the verbosity. 1121 @type verbosity: int 1122 @return: The assembled data structures for using RDCs as the base data for optimisation. These include: 1123 - rdc, the RDC values. 1124 - rdc_err, the RDC errors. 1125 - rdc_weight, the RDC weights. 1126 - vectors, the interatomic vectors (pseudo-atom dependent). 1127 - rdc_const, the dipolar constants (pseudo-atom dependent). 1128 - absolute, the absolute value flags (as 1's and 0's). 1129 - T_flags, the flags for T = J+D type data (as 1's and 0's). 1130 - j_couplings, the J coupling values if the RDC data type is set to T = J+D. 1131 - pseudo_flags, the list of flags indicating if the interatomic data contains a pseudo-atom (as 1's and 0's). 1132 @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) 1133 """ 1134 1135 # Initial printout. 1136 if verbosity: 1137 print("\nRDC data counts:") 1138 1139 # Sort out pseudo-atoms first. This only needs to be called once. 1140 setup_pseudoatom_rdc() 1141 1142 # Initialise. 1143 rdc = [] 1144 rdc_err = [] 1145 rdc_weight = [] 1146 unit_vect = [] 1147 rdc_const = [] 1148 absolute = [] 1149 T_flags = [] 1150 j_couplings = [] 1151 pseudo_flags = [] 1152 1153 # The unit vectors, RDC constants, and J couplings. 1154 for interatom in interatomic_loop(): 1155 # Get the spins. 1156 spin1 = return_spin(spin_hash=interatom._spin_hash1) 1157 spin2 = return_spin(spin_hash=interatom._spin_hash2) 1158 1159 # RDC checks. 1160 if not check_rdcs(interatom): 1161 continue 1162 1163 # Gyromagnetic ratios. 1164 g1 = periodic_table.gyromagnetic_ratio(spin1.isotope) 1165 g2 = periodic_table.gyromagnetic_ratio(spin2.isotope) 1166 1167 # Pseudo atoms. 1168 if is_pseudoatom(spin1) and is_pseudoatom(spin2): 1169 raise RelaxError("Support for both spins being in a dipole pair being pseudo-atoms is not implemented yet.") 1170 if is_pseudoatom(spin1) or is_pseudoatom(spin2): 1171 # Set the flag. 1172 pseudo_flags.append(1) 1173 1174 # Alias the pseudo and normal atoms. 1175 if is_pseudoatom(spin1): 1176 pseudospin = spin1 1177 base_spin = spin2 1178 pseudospin_id = interatom.spin_id1 1179 base_spin_id = interatom.spin_id2 1180 base_spin_hash = interatom._spin_hash2 1181 else: 1182 pseudospin = spin2 1183 base_spin = spin1 1184 pseudospin_id = interatom.spin_id2 1185 base_spin_id = interatom.spin_id1 1186 base_spin_hash = interatom._spin_hash1 1187 1188 # Loop over the atoms of the pseudo-atom, storing the data. 1189 pseudo_unit_vect = [] 1190 pseudo_rdc_const = [] 1191 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 1192 # Get the corresponding interatomic data container. 1193 pseudo_interatom = return_interatom(spin_hash1=spin._hash, spin_hash2=base_spin_hash) 1194 1195 # Check. 1196 if pseudo_interatom == None: 1197 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)) 1198 1199 # Add the vectors. 1200 if is_float(interatom.vector[0]): 1201 pseudo_unit_vect.append([pseudo_interatom.vector]) 1202 else: 1203 pseudo_unit_vect.append(pseudo_interatom.vector) 1204 1205 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 1206 pseudo_rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, pseudo_interatom.r)) 1207 1208 # Reorder the unit vectors so that the structure and pseudo-atom dimensions are swapped. 1209 pseudo_unit_vect = transpose(array(pseudo_unit_vect, float64), (1, 0, 2)) 1210 1211 # Block append the pseudo-data. 1212 unit_vect.append(pseudo_unit_vect) 1213 rdc_const.append(pseudo_rdc_const) 1214 1215 # Normal atom. 1216 else: 1217 # Set the flag. 1218 pseudo_flags.append(0) 1219 1220 # Add the vectors. 1221 if is_float(interatom.vector[0]): 1222 unit_vect.append([interatom.vector]) 1223 else: 1224 unit_vect.append(interatom.vector) 1225 1226 # Calculate the RDC dipolar constant (in Hertz, and the 3 comes from the alignment tensor), and append it to the list. 1227 rdc_const.append(3.0/(2.0*pi) * dipolar_constant(g1, g2, interatom.r)) 1228 1229 # Sanity check, to prevent cryptic Python errors. 1230 indices = [] 1231 for i in range(len(unit_vect[-1])): 1232 if unit_vect[-1][i] is None: 1233 raise RelaxError("Unit vectors of None have been detected between the spins '%s' and '%s' %s." % (interatom.spin_id1, interatom.spin_id2, unit_vect[-1])) 1234 1235 # Store the measured J coupling. 1236 if opt_uses_j_couplings(): 1237 j_couplings.append(interatom.j_coupling) 1238 1239 # Fix the unit vector data structure. 1240 num = None 1241 for rdc_index in range(len(unit_vect)): 1242 # Convert to numpy structures. 1243 unit_vect[rdc_index] = array(unit_vect[rdc_index], float64) 1244 1245 # Number of vectors. 1246 if num == None: 1247 if unit_vect[rdc_index] is not None: 1248 num = len(unit_vect[rdc_index]) 1249 continue 1250 1251 # Check. 1252 if unit_vect[rdc_index] is not None and len(unit_vect[rdc_index]) != num: 1253 raise RelaxError("The number of interatomic vectors for all no match:\n%s" % unit_vect[rdc_index]) 1254 1255 # Missing unit vectors. 1256 if num == None: 1257 raise RelaxError("No interatomic vectors could be found.") 1258 1259 # Update None entries. 1260 for i in range(len(unit_vect)): 1261 if unit_vect[i] is None: 1262 unit_vect[i] = [[None, None, None]]*num 1263 1264 # The RDC data. 1265 for i in range(len(cdp.align_ids)): 1266 # Alias the ID. 1267 align_id = cdp.align_ids[i] 1268 1269 # Skip non-optimised data. 1270 if not opt_uses_align_data(align_id): 1271 continue 1272 1273 # Append empty arrays to the RDC structures. 1274 rdc.append([]) 1275 rdc_err.append([]) 1276 rdc_weight.append([]) 1277 absolute.append([]) 1278 T_flags.append([]) 1279 1280 # Interatom loop. 1281 j = 0 1282 for interatom in interatomic_loop(): 1283 # Get the spins. 1284 spin1 = return_spin(spin_hash=interatom._spin_hash1) 1285 spin2 = return_spin(spin_hash=interatom._spin_hash2) 1286 1287 # RDC checks. 1288 if not check_rdcs(interatom): 1289 continue 1290 1291 # T-type data. 1292 if align_id in interatom.rdc_data_types and interatom.rdc_data_types[align_id] == 'T': 1293 T_flags[-1].append(True) 1294 else: 1295 T_flags[-1].append(False) 1296 1297 # Check for J couplings if the RDC data type is T = J+D. 1298 if T_flags[-1][-1] and not hasattr(interatom, 'j_coupling'): 1299 continue 1300 1301 # Defaults of None. 1302 value = None 1303 error = None 1304 1305 # Normal set up. 1306 if align_id in interatom.rdc: 1307 # The RDC. 1308 if sim_index != None: 1309 value = interatom.rdc_sim[align_id][sim_index] 1310 else: 1311 value = interatom.rdc[align_id] 1312 1313 # The error. 1314 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 1315 # T values. 1316 if T_flags[-1][-1]: 1317 error = sqrt(interatom.rdc_err[align_id]**2 + interatom.j_coupling_err**2) 1318 1319 # D values. 1320 else: 1321 error = interatom.rdc_err[align_id] 1322 1323 # Increment the RDC count. 1324 if value != None: 1325 j += 1 1326 1327 # Append the RDCs to the list. 1328 rdc[-1].append(value) 1329 1330 # Append the RDC errors. 1331 rdc_err[-1].append(error) 1332 1333 # Append the weight. 1334 if hasattr(interatom, 'rdc_weight') and align_id in interatom.rdc_weight: 1335 rdc_weight[-1].append(interatom.rdc_weight[align_id]) 1336 else: 1337 rdc_weight[-1].append(1.0) 1338 1339 # Append the absolute value flag. 1340 if hasattr(interatom, 'absolute_rdc') and align_id in interatom.absolute_rdc: 1341 absolute[-1].append(interatom.absolute_rdc[align_id]) 1342 else: 1343 absolute[-1].append(False) 1344 1345 # ID and RDC count printout. 1346 if verbosity: 1347 print(" Alignment ID '%s': %i" % (align_id, j)) 1348 1349 # Convert to numpy objects. 1350 rdc = array(rdc, float64) 1351 rdc_err = array(rdc_err, float64) 1352 rdc_weight = array(rdc_weight, float64) 1353 absolute = array(absolute, int32) 1354 T_flags = array(T_flags, int32) 1355 if not opt_uses_j_couplings(): 1356 j_couplings = None 1357 pseudo_flags = array(pseudo_flags, int32) 1358 1359 # Return the data structures. 1360 return rdc, rdc_err, rdc_weight, unit_vect, rdc_const, absolute, T_flags, j_couplings, pseudo_flags

1361 1362

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

1364 """Set the RDC errors if not already present. 1365 1366 @keyword align_id: The optional alignment tensor ID string. 1367 @type align_id: str 1368 @keyword spin_id1: The optional spin ID string of the first spin. 1369 @type spin_id1: None or str 1370 @keyword spin_id2: The optional spin ID string of the second spin. 1371 @type spin_id2: None or str 1372 @keyword sd: The RDC standard deviation in Hz. 1373 @type sd: float or int. 1374 """ 1375 1376 # Check the pipe setup. 1377 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1378 1379 # Convert the align IDs to an array, or take all IDs. 1380 if align_id: 1381 align_ids = [align_id] 1382 else: 1383 align_ids = cdp.rdc_ids 1384 1385 # Loop over the interatomic data. 1386 for interatom in interatomic_loop(selection1=spin_id1, selection2=spin_id2): 1387 # No data structure. 1388 if not hasattr(interatom, 'rdc_err'): 1389 interatom.rdc_err = {} 1390 1391 # Set the error. 1392 for id in align_ids: 1393 interatom.rdc_err[id] = sd

1394 1395

1396 -def setup_pseudoatom_rdc():

1397 """Make sure that the interatom system is properly set up for pseudo-atoms and RDCs. 1398 1399 Interatomic data containers between the non-pseudo-atom and the pseudo-atom members will be deselected. 1400 """ 1401 1402 # Loop over all interatomic data containers. 1403 for interatom in interatomic_loop(): 1404 # Get the spins. 1405 spin1 = return_spin(spin_hash=interatom._spin_hash1) 1406 spin2 = return_spin(spin_hash=interatom._spin_hash2) 1407 1408 # Checks. 1409 flag1 = is_pseudoatom(spin1) 1410 flag2 = is_pseudoatom(spin2) 1411 1412 # No pseudo-atoms, so do nothing. 1413 if not (flag1 or flag2): 1414 continue 1415 1416 # Both are pseudo-atoms. 1417 if flag1 and flag2: 1418 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))) 1419 interatom.select = False 1420 1421 # Alias the pseudo and normal atoms. 1422 pseudospin = spin1 1423 base_spin_id = interatom.spin_id2 1424 base_spin_hash = interatom._spin_hash2 1425 pseudospin_id = interatom.spin_id1 1426 if flag2: 1427 pseudospin = spin2 1428 base_spin_id = interatom.spin_id1 1429 base_spin_hash = interatom._spin_hash1 1430 pseudospin_id = interatom.spin_id2 1431 1432 # Loop over the atoms of the pseudo-atom. 1433 for spin, spin_id in pseudoatom_loop(pseudospin, return_id=True): 1434 # Get the corresponding interatomic data container. 1435 pseudo_interatom = return_interatom(spin_hash1=spin._hash, spin_hash2=base_spin_hash) 1436 1437 # Deselect if needed. 1438 if pseudo_interatom.select: 1439 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))) 1440 pseudo_interatom.select = False

1441 1442

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

1444 """Set optimisation weights on the RDC data. 1445 1446 @keyword align_id: The alignment tensor ID string. 1447 @type align_id: str 1448 @keyword spin_id: The spin ID string. 1449 @type spin_id: None or str 1450 @keyword weight: The optimisation weight. The higher the value, the more importance the RDC will have. 1451 @type weight: float or int. 1452 """ 1453 1454 # Check the pipe setup. 1455 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1456 1457 # Loop over the interatomic data. 1458 for interatom in interatomic_loop(): 1459 # No data structure. 1460 if not hasattr(interatom, 'rdc_weight'): 1461 interatom.rdc_weight = {} 1462 1463 # Set the weight. 1464 interatom.rdc_weight[align_id] = weight

1465 1466

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

1468 """Display the RDC data corresponding to the alignment ID. 1469 1470 @keyword align_id: The alignment tensor ID string. 1471 @type align_id: str 1472 @keyword file: The file name or object to write to. 1473 @type file: str or file object 1474 @keyword dir: The name of the directory to place the file into (defaults to the current directory). 1475 @type dir: str 1476 @keyword bc: The back-calculation flag which if True will cause the back-calculated rather than measured data to be written. 1477 @type bc: bool 1478 @keyword force: A flag which if True will cause any pre-existing file to be overwritten. 1479 @type force: bool 1480 """ 1481 1482 # Check the pipe setup. 1483 check_pipe_setup(sequence=True, rdc_id=align_id, rdc=True) 1484 1485 # Open the file for writing. 1486 file = open_write_file(file, dir, force) 1487 1488 # Loop over the interatomic data containers and collect the data. 1489 data = [] 1490 for interatom in interatomic_loop(): 1491 # Skip deselected containers. 1492 if not interatom.select: 1493 continue 1494 1495 # Skip containers with no RDCs. 1496 if not bc and (not hasattr(interatom, 'rdc') or align_id not in interatom.rdc): 1497 continue 1498 elif bc and (not hasattr(interatom, 'rdc_bc') or align_id not in interatom.rdc_bc): 1499 continue 1500 1501 # Append the spin data. 1502 data.append([]) 1503 data[-1].append(repr(interatom.spin_id1)) 1504 data[-1].append(repr(interatom.spin_id2)) 1505 1506 # Handle the missing rdc_data_types variable. 1507 data_type = None 1508 if hasattr(interatom, 'rdc_data_types') and align_id in interatom.rdc_data_types: 1509 data_type = interatom.rdc_data_types[align_id] 1510 1511 # The value. 1512 if bc: 1513 data[-1].append(repr(convert(interatom.rdc_bc[align_id], data_type, align_id))) 1514 else: 1515 data[-1].append(repr(convert(interatom.rdc[align_id], data_type, align_id))) 1516 1517 # The error. 1518 if hasattr(interatom, 'rdc_err') and align_id in interatom.rdc_err: 1519 data[-1].append(repr(convert(interatom.rdc_err[align_id], data_type, align_id))) 1520 else: 1521 data[-1].append(repr(None)) 1522 1523 # Write out. 1524 write_data(out=file, headings=["Spin_ID1", "Spin_ID2", "RDCs", "RDC_error"], data=data)

1525

Source Code for Module pipe_control.rdc