pipe_control.rdc

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2015 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 ceil, floor, 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.periodic_table import periodic_table 40 from lib.physical_constants import dipolar_constant 41 from lib.plotting.api import write_xy_data, write_xy_header 42 from lib.warnings import RelaxWarning, RelaxSpinTypeWarning 43 from pipe_control import pipes 44 from pipe_control.align_tensor import get_tensor_index, get_tensor_object, opt_uses_align_data, opt_uses_tensor 45 from pipe_control.interatomic import consistent_interatomic_data, create_interatom, interatomic_loop, return_interatom 46 from pipe_control.mol_res_spin import exists_mol_res_spin_data, is_pseudoatom, pseudoatom_loop, return_spin 47 from pipe_control.pipes import check_pipe 48 49

50 -def back_calc(align_id=None):

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

142 143

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

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

194 195

196 -def check_rdcs(interatom):

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

272 273

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

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

306 307

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

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

442 443

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

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

568 569

570 -def delete(align_id=None):

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

608 609

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

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

624 625

626 -def opt_uses_j_couplings():

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

641 642

643 -def opt_uses_rdc(align_id):

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

679 680

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

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

909 910

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

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

1110 1111

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

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

1356 1357

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

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

1389 1390

1391 -def setup_pseudoatom_rdc():

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

1434 1435

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

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

1458 1459

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

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

1518

Source Code for Module pipe_control.rdc