pipe_control.align

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 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 containing functions for the handling of alignment tensors.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import pi, sqrt 28 from numpy import arccos, dot, float64, linalg, zeros 29 from numpy.linalg import norm 30 import sys 31 from warnings import warn 32 33 # relax module imports. 34 from data_store.align_tensor import AlignTensorList 35 from lib.alignment.alignment_tensor import calc_chi_tensor, kappa 36 from lib.errors import RelaxError, RelaxNoTensorError, RelaxTensorError, RelaxUnknownParamCombError, RelaxUnknownParamError 37 from lib.geometry.angles import wrap_angles 38 from lib.io import write_data 39 from lib.warnings import RelaxWarning 40 import pipe_control 41 from pipe_control import pipes 42 43

44 -def align_data_exists(tensor, pipe=None):

45 """Function for determining if alignment data exists in the current data pipe. 46 47 @param tensor: The alignment tensor identification string. 48 @type tensor: str 49 @param pipe: The data pipe to search for data in. 50 @type pipe: str 51 @return: The answer to the question. 52 @rtype: bool 53 """ 54 55 # The data pipe to check. 56 if pipe == None: 57 pipe = pipes.cdp_name() 58 59 # Get the data pipe. 60 pipe = pipes.get_pipe(pipe) 61 62 # Test if an alignment tensor corresponding to the arg 'tensor' exists. 63 if hasattr(pipe, 'align_tensors'): 64 for data in pipe.align_tensors: 65 if data.name == tensor: 66 return True 67 else: 68 return False

69 70

71 -def all_tensors_fixed():

72 """Determine if all alignment tensors are fixed. 73 74 @return: True if all tensors are fixed, False otherwise. 75 @rtype: bool 76 """ 77 78 # Loop over the tensors. 79 for i in range(len(cdp.align_tensors)): 80 # Not fixed, so return False. 81 if not cdp.align_tensors[i].fixed: 82 return False 83 84 # All tensors are fixed. 85 return True

86 87

88 -def copy(tensor_from=None, pipe_from=None, tensor_to=None, pipe_to=None):

89 """Function for copying alignment tensor data from one data pipe to another. 90 91 @param tensor_from: The identification string of the alignment tensor to copy the data from. 92 @type tensor_from: str 93 @param pipe_from: The data pipe to copy the alignment tensor data from. This defaults to the current data pipe. 94 @type pipe_from: str 95 @param tensor_to: The identification string of the alignment tensor to copy the data to. If set to None, then the ID string will be set to the value of tensor_from. 96 @type tensor_to: str or None 97 @param pipe_to: The data pipe to copy the alignment tensor data to. This defaults to the current data pipe. 98 @type pipe_to: str 99 """ 100 101 # Defaults. 102 if tensor_from == tensor_to and pipe_from == None and pipe_to == None: 103 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None when the tensor names are the same.") 104 elif pipe_from == None: 105 pipe_from = pipes.cdp_name() 106 elif pipe_to == None: 107 pipe_to = pipes.cdp_name() 108 109 # The target tensor ID string. 110 if tensor_to == None: 111 tensor_to = tensor_from 112 113 # Test if the pipe_from and pipe_to data pipes exist. 114 pipes.test(pipe_from) 115 pipes.test(pipe_to) 116 117 # Get the data pipes. 118 dp_from = pipes.get_pipe(pipe_from) 119 dp_to = pipes.get_pipe(pipe_to) 120 121 # Test if pipe_from contains alignment tensor data. 122 if not align_data_exists(tensor_from, pipe_from): 123 raise RelaxNoTensorError('alignment') 124 125 # Test if pipe_to contains alignment tensor data. 126 if align_data_exists(tensor_to, pipe_to): 127 raise RelaxTensorError('alignment') 128 129 # Create the align_tensors dictionary if it doesn't yet exist. 130 if not hasattr(dp_to, 'align_tensors'): 131 dp_to.align_tensors = AlignTensorList() 132 133 # Add the tensor if it doesn't already exist. 134 if tensor_to not in dp_to.align_tensors.names(): 135 dp_to.align_tensors.add_item(tensor_to) 136 137 # Find the tensor index. 138 index_from = get_tensor_index(tensor=tensor_from, pipe=pipe_from) 139 index_to = get_tensor_index(tensor=tensor_to, pipe=pipe_to) 140 141 # Copy the data. 142 if index_to == None: 143 dp_to.align_tensors.append(deepcopy(dp_from.align_tensors[index_from])) 144 index_to = len(dp_to.align_tensors) - 1 145 else: 146 dp_to.align_tensors[index_to] = deepcopy(dp_from.align_tensors[index_from]) 147 148 # Update the tensor's name. 149 dp_to.align_tensors[index_to].set('name', tensor_to) 150 151 # Add the align ID to the target data pipe if needed. 152 align_id = dp_from.align_tensors[index_from].align_id 153 if not hasattr(dp_to, 'align_ids'): 154 dp_to.align_ids = [] 155 if align_id not in dp_to.align_ids: 156 dp_to.align_ids.append(align_id)

157 158

159 -def delete(tensor=None):

160 """Function for deleting alignment tensor data. 161 162 @param tensor: The alignment tensor identification string. 163 @type tensor: str or None 164 """ 165 166 # Test if the current data pipe exists. 167 pipes.test() 168 169 # Test if alignment tensor data exists. 170 if tensor and not align_data_exists(tensor): 171 raise RelaxNoTensorError('alignment') 172 173 # The tensor list. 174 if tensor: 175 tensors = [tensor] 176 else: 177 tensors = [] 178 for i in range(len(cdp.align_tensors)): 179 tensors.append(cdp.align_tensors[i].name) 180 181 # Loop over the tensors. 182 for tensor in tensors: 183 # Print out. 184 print("Removing the '%s' tensor." % tensor) 185 186 # Find the tensor index. 187 index = get_tensor_index(tensor=tensor) 188 189 # Delete the alignment data. 190 cdp.align_tensors.pop(index) 191 192 # Delete the alignment tensor list if empty. 193 if not len(cdp.align_tensors): 194 del(cdp.align_tensors)

195 196

197 -def display(tensor=None):

198 """Function for displaying the alignment tensor. 199 200 @keyword tensor: The alignment tensor identification string. 201 @type tensor: str or None 202 """ 203 204 # Test if the current data pipe exists. 205 pipes.test() 206 207 # Construct the tensor list. 208 tensor_list = [] 209 if not tensor: 210 for tensor_cont in cdp.align_tensors: 211 tensor_list.append(tensor_cont.name) 212 else: 213 tensor_list.append(tensor) 214 215 # Loop over the tensors. 216 for tensor in tensor_list: 217 # Test if alignment tensor data exists. 218 if not align_data_exists(tensor): 219 raise RelaxNoTensorError('alignment') 220 221 # Pull out the tensor. 222 data = get_tensor_object(tensor) 223 224 # Header. 225 head = "# Tensor: %s #" % tensor 226 print("\n\n\n" + '#' * len(head) + "\n" + head + "\n" + '#' * len(head)) 227 228 229 # The Saupe matrix. 230 ################### 231 232 title = "# Saupe order matrix." 233 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 234 235 # The parameter set {Sxx, Syy, Sxy, Sxz, Syz}. 236 print("# 5D, rank-1 notation {Sxx, Syy, Sxy, Sxz, Syz}:") 237 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Sxx, data.Syy, data.Sxy, data.Sxz, data.Syz)) 238 239 # The parameter set {Szz, Sxx-yy, Sxy, Sxz, Syz}. 240 print("# 5D, rank-1 notation {Szz, Sxx-yy, Sxy, Sxz, Syz} (the Pales default format).") 241 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Szz, data.Sxxyy, data.Sxy, data.Sxz, data.Syz)) 242 243 # 3D form. 244 print("# 3D, rank-2 notation.") 245 print("%s" % (data.S)) 246 247 248 # The alignment tensor. 249 ####################### 250 251 title = "# Alignment tensor." 252 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 253 254 # The parameter set {Axx, Ayy, Axy, Axz, Ayz}. 255 print("# 5D, rank-1 notation {Axx, Ayy, Axy, Axz, Ayz}:") 256 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Axx, data.Ayy, data.Axy, data.Axz, data.Ayz)) 257 258 # The parameter set {Azz, Axx-yy, Axy, Axz, Ayz}. 259 print("# 5D, rank-1 notation {Azz, Axx-yy, Axy, Axz, Ayz} (the Pales default format).") 260 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Azz, data.Axxyy, data.Axy, data.Axz, data.Ayz)) 261 262 # 3D form. 263 print("# 3D, rank-2 notation.") 264 print("%s" % data.A) 265 266 267 # The probability tensor. 268 ######################### 269 270 title = "# Probability tensor." 271 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 272 273 # The parameter set {Pxx, Pyy, Pxy, Pxz, Pyz}. 274 print("# 5D, rank-1 notation {Pxx, Pyy, Pxy, Pxz, Pyz}:") 275 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Pxx, data.Pyy, data.Pxy, data.Pxz, data.Pyz)) 276 277 # The parameter set {Pzz, Pxx-yy, Pxy, Pxz, Pyz}. 278 print("# 5D, rank-1 notation {Pzz, Pxx-yy, Pxy, Pxz, Pyz}.") 279 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (data.Pzz, data.Pxxyy, data.Pxy, data.Pxz, data.Pyz)) 280 281 # 3D form. 282 print("# 3D, rank-2 notation.") 283 print("%s" % data.P) 284 285 286 # The magnetic susceptibility tensor. 287 ##################################### 288 289 title = "# Magnetic susceptibility tensor." 290 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 291 chi_tensor = True 292 293 # The field strength. 294 print("# The magnetic field strength (MHz):") 295 if hasattr(cdp, 'spectrometer_frq') and tensor in cdp.spectrometer_frq: 296 print("%s\n" % (cdp.spectrometer_frq[tensor] / 1e6)) 297 else: 298 print("Not set.\n") 299 chi_tensor = False 300 301 # The temperature. 302 print("# The temperature (K):") 303 if hasattr(cdp, 'temperature') and tensor in cdp.temperature: 304 print("%s\n" % cdp.temperature[tensor]) 305 else: 306 print("Not set.\n") 307 chi_tensor = False 308 309 # No chi tensor. 310 if not chi_tensor: 311 print("# The chi tensor:\nN/A.\n") 312 313 # Calculate the chi tensor. 314 else: 315 # Conversions. 316 chi_xx = calc_chi_tensor(data.Axx, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 317 chi_xy = calc_chi_tensor(data.Axy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 318 chi_xz = calc_chi_tensor(data.Axz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 319 chi_yy = calc_chi_tensor(data.Ayy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 320 chi_yz = calc_chi_tensor(data.Ayz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 321 chi_zz = calc_chi_tensor(data.Azz, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 322 chi_xxyy = calc_chi_tensor(data.Axxyy, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 323 chi = calc_chi_tensor(data.A, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 324 325 # The parameter set {chi_xx, chi_yy, chi_xy, chi_xz, chi_yz}. 326 print("# 5D, rank-1 notation {chi_xx, chi_yy, chi_xy, chi_xz, chi_yz}:") 327 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (chi_xx, chi_yy, chi_xy, chi_xz, chi_yz)) 328 329 # The parameter set {chi_zz, chi_xx-yy, chi_xy, chi_xz, chi_yz}. 330 print("# 5D, rank-1 notation {chi_zz, chi_xx-yy, chi_xy, chi_xz, chi_yz}.") 331 print("[%25.12e, %25.12e, %25.12e, %25.12e, %25.12e]\n" % (chi_zz, chi_xxyy, chi_xy, chi_xz, chi_yz)) 332 333 # 3D form. 334 print("# 3D, rank-2 notation.") 335 print("%s" % chi) 336 337 338 # The Eigensystem. 339 ################## 340 341 title = "# Eigensystem." 342 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 343 344 # Eigenvalues. 345 print("# Saupe order matrix eigenvalues {Sxx, Syy, Szz}.") 346 print("[%25.12e, %25.12e, %25.12e]\n" % (data.S_diag[0, 0], data.S_diag[1, 1], data.S_diag[2, 2])) 347 print("# Alignment tensor eigenvalues {Axx, Ayy, Azz}.") 348 print("[%25.12e, %25.12e, %25.12e]\n" % (data.A_diag[0, 0], data.A_diag[1, 1], data.A_diag[2, 2])) 349 print("# Probability tensor eigenvalues {Pxx, Pyy, Pzz}.") 350 print("[%25.12e, %25.12e, %25.12e]\n" % (data.P_diag[0, 0], data.P_diag[1, 1], data.P_diag[2, 2])) 351 if chi_tensor: 352 chi_diag = calc_chi_tensor(data.A_diag, cdp.spectrometer_frq[tensor], cdp.temperature[tensor]) 353 print("# Magnetic susceptibility eigenvalues {chi_xx, chi_yy, chi_zz}.") 354 print("[%25.12e, %25.12e, %25.12e]\n" % (chi_diag[0, 0], chi_diag[1, 1], chi_diag[2, 2])) 355 356 # Eigenvectors. 357 print("# Eigenvector x.") 358 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_x[0], data.unit_x[1], data.unit_x[2])) 359 print("# Eigenvector y.") 360 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_y[0], data.unit_y[1], data.unit_y[2])) 361 print("# Eigenvector z.") 362 print("[%25.12f, %25.12f, %25.12f]\n" % (data.unit_z[0], data.unit_z[1], data.unit_z[2])) 363 364 # Rotation matrix. 365 print("# Rotation matrix.") 366 print("%s\n" % data.rotation) 367 368 # zyz. 369 print("# Euler angles in zyz notation {alpha, beta, gamma}.") 370 print("[%25.12f, %25.12f, %25.12f]\n" % (data.euler[0], data.euler[1], data.euler[2])) 371 372 373 # Geometric description. 374 ######################## 375 376 title = "# Geometric description." 377 print("\n\n" + title + '\n' + '#'*len(title) + '\n') 378 379 # The GDO. 380 print("# Generalized degree of order (GDO).") 381 print("GDO = %-25.12e\n" % gdo(data.A)) 382 383 # Anisotropy. 384 print("# Alignment tensor axial component (Aa = 3/2 * Azz, where Aii are the eigenvalues).") 385 print("Aa = %-25.12e\n" % data.Aa) 386 387 # Rhombicity. 388 print("# Rhombic component (Ar = Axx - Ayy, where Aii are the eigenvalues).") 389 print("Ar = %-25.12e\n" % data.Ar) 390 print("# Rhombicity (R = Ar / Aa).") 391 print("R = %-25.12f\n" % data.R) 392 print("# Asymmetry parameter (eta = (Axx - Ayy) / Azz, where Aii are the eigenvalues).") 393 print("eta = %-25.12f\n" % data.eta) 394 395 # Magnetic susceptibility tensor. 396 if chi_tensor: 397 # Chi tensor anisotropy. 398 print("# Magnetic susceptibility axial parameter (chi_ax = chi_zz - (chi_xx + chi_yy)/2, where chi_ii are the eigenvalues).") 399 print("chi_ax = %-25.12e\n" % (chi_diag[2, 2] - (chi_diag[0, 0] + chi_diag[1, 1])/2.0)) 400 401 # Chi tensor rhombicity. 402 print("# Magnetic susceptibility rhombicity parameter (chi_rh = chi_xx - chi_yy, where chi_ii are the eigenvalues).") 403 print("chi_rh = %-25.12e\n" % (chi_diag[0, 0] - chi_diag[1, 1])) 404 405 # Some white space. 406 print("\n\n\n")

407 408

409 -def fix(id=None, fixed=True):

410 """Fix the alignment tensor during optimisation. 411 412 @keyword id: The alignment tensor ID string. If set to None, then all alignment tensors will be fixed. 413 @type id: str or None 414 @keyword fixed: If True, the alignment tensor will be fixed during optimisation. If False, the alignment tensors will be optimised. 415 @type fixed: bool 416 """ 417 418 # Test if the current data pipe exists. 419 pipes.test() 420 421 # Loop over the tensors. 422 for i in range(len(cdp.align_tensors)): 423 # ID match. 424 if id and cdp.align_tensors[i].name == id: 425 cdp.align_tensors[i].set_fixed(fixed) 426 427 # Set all tensor flags. 428 if id == None: 429 cdp.align_tensors[i].set_fixed(fixed)

430 431

432 -def fold_angles(sim_index=None):

433 """Wrap the Euler angles and remove the glide reflection and translational symmetries. 434 435 Wrap the angles such that:: 436 437 0 <= alpha <= 2pi, 438 0 <= beta <= pi, 439 0 <= gamma <= 2pi. 440 441 For the simulated values, the angles are wrapped as:: 442 443 alpha - pi <= alpha_sim <= alpha + pi 444 beta - pi/2 <= beta_sim <= beta + pi/2 445 gamma - pi <= gamma_sim <= gamma + pi 446 447 448 @param sim_index: The simulation index. If set to None then the actual values will be folded 449 rather than the simulation values. 450 @type sim_index: int or None 451 """ 452 453 454 # Wrap the angles. 455 ################## 456 457 # Get the current angles. 458 alpha = cdp.align_tensors.alpha 459 beta = cdp.align_tensors.beta 460 gamma = cdp.align_tensors.gamma 461 462 # Simulated values. 463 if sim_index != None: 464 alpha_sim = cdp.align_tensors.alpha_sim[sim_index] 465 beta_sim = cdp.align_tensors.beta_sim[sim_index] 466 gamma_sim = cdp.align_tensors.gamma_sim[sim_index] 467 468 # Normal value. 469 if sim_index == None: 470 cdp.align_tensors.set(param='alpha', value=wrap_angles(alpha, 0.0, 2.0*pi)) 471 cdp.align_tensors.set(param='beta', value= wrap_angles(beta, 0.0, 2.0*pi)) 472 cdp.align_tensors.set(param='gamma', value=wrap_angles(gamma, 0.0, 2.0*pi)) 473 474 # Simulation values. 475 else: 476 cdp.align_tensors.set(param='alpha', value=wrap_angles(alpha_sim, alpha - pi, alpha + pi), category='sim', sim_index=sim_index) 477 cdp.align_tensors.set(param='beta', value= wrap_angles(beta_sim, beta - pi, beta + pi), category='sim', sim_index=sim_index) 478 cdp.align_tensors.set(param='gamma', value=wrap_angles(gamma_sim, gamma - pi, gamma + pi), category='sim', sim_index=sim_index) 479 480 481 # Remove the glide reflection and translational symmetries. 482 ########################################################### 483 484 # Normal value. 485 if sim_index == None: 486 # Fold beta inside 0 and pi. 487 if cdp.align_tensors.beta >= pi: 488 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha) 489 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta - pi) 490 491 # Simulation values. 492 else: 493 # Fold beta_sim inside beta-pi/2 and beta+pi/2. 494 if cdp.align_tensors.beta_sim[sim_index] >= cdp.align_tensors.beta + pi/2.0: 495 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha_sim[sim_index], category='sim', sim_index=sim_index) 496 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta_sim[sim_index] - pi, category='sim', sim_index=sim_index) 497 elif cdp.align_tensors.beta_sim[sim_index] <= cdp.align_tensors.beta - pi/2.0: 498 cdp.align_tensors.set(param='alpha', value=pi - cdp.align_tensors.alpha_sim[sim_index], category='sim', sim_index=sim_index) 499 cdp.align_tensors.set(param='beta', value=cdp.align_tensors.beta_sim[sim_index] + pi, category='sim', sim_index=sim_index)

500 501

502 -def gdo(A):

503 """Calculate the generalized degree of order (GDO) for the given alignment tensor. 504 505 @param A: The alignment tensor. 506 @type A: rank-2, 3D numpy array 507 @return: The GDO value. 508 @rtype: float 509 """ 510 511 # The matrix norm. 512 gdo = sqrt(3.0/2.0) * norm(A) 513 514 # Return the GDO. 515 return gdo

516 517

518 -def get_align_ids():

519 """Return the list of all alignment IDs. 520 521 @return: The list of all alignment IDs. 522 @rtype: list of str 523 """ 524 525 # No pipe. 526 if cdp == None: 527 return [] 528 529 # No tensor data. 530 if not hasattr(cdp, 'align_ids'): 531 return [] 532 533 # The tensor IDs. 534 return cdp.align_ids

535 536

537 -def get_tensor_ids():

538 """Return the list of all tensor IDs. 539 540 @return: The list of all tensor IDs. 541 @rtype: list of str 542 """ 543 544 # Init. 545 ids = [] 546 547 # No pipe. 548 if cdp == None: 549 return ids 550 551 # No tensor data. 552 if not hasattr(cdp, 'align_tensors'): 553 return ids 554 555 # Loop over the tensors. 556 for i in range(len(cdp.align_tensors)): 557 if cdp.align_tensors[i].name != None: 558 ids.append(cdp.align_tensors[i].name) 559 560 # Return the object. 561 return ids

562 563

564 -def get_tensor_index(tensor=None, align_id=None, pipe=None):

565 """Function for returning the index corresponding to the 'tensor' argument. 566 567 @keyword tensor: The alignment tensor identification string. 568 @type tensor: str or None 569 @keyword align_id: Alternative to the tensor argument, used to return the tensor index for the tensors corresponding to the alignment ID string. If more than one tensor exists, then this will fail. 570 @type align_id: str or None 571 @keyword pipe: The data pipe to search for data in. 572 @type pipe: str 573 @return: The index corresponding to the 'tensor' arg. 574 @rtype: int 575 """ 576 577 # The data pipe to check. 578 if pipe == None: 579 pipe = pipes.cdp_name() 580 581 # Get the data pipe. 582 dp = pipes.get_pipe(pipe) 583 584 # Init. 585 index = None 586 count = 0 587 588 # Loop over the tensors. 589 for i in range(len(dp.align_tensors)): 590 # Tensor name match. 591 if tensor and dp.align_tensors[i].name == tensor: 592 index = i 593 count += 1 594 595 # Alignment ID match. 596 if align_id and hasattr(dp.align_tensors[i], 'align_id') and dp.align_tensors[i].align_id == align_id: 597 index = i 598 count += 1 599 600 # No match. 601 if count == 0: 602 warn(RelaxWarning("No alignment tensors match the tensor name '%s' or alignment ID '%s' in the data pipe '%s'." % (tensor, align_id, pipe))) 603 return None 604 605 # More than one match. 606 if count > 1: 607 warn(RelaxWarning("More than one alignment tensors matches the tensor name '%s' or alignment ID '%s' in the data pipe '%s'." % (tensor, align_id, pipe))) 608 return None 609 610 # Return the index. 611 return index

612 613

614 -def get_tensor_object(tensor, pipe=None):

615 """Return the AlignTensorData instance corresponding to the tensor ID. 616 617 @param tensor: The alignment tensor identification string. 618 @type tensor: str 619 @param pipe: The data pipe to search for data in. 620 @type pipe: str 621 @return: The alignment tensor object corresponding to the 'tensor' arg. 622 @rtype: AlignTensorData instance 623 """ 624 625 # The data pipe to check. 626 if pipe == None: 627 pipe = pipes.cdp_name() 628 629 # Init. 630 data = None 631 632 # Loop over the tensors. 633 for i in range(len(cdp.align_tensors)): 634 if cdp.align_tensors[i].name == tensor: 635 data = cdp.align_tensors[i] 636 637 # Return the object. 638 return data

639 640

641 -def get_tensor_object_from_align(align_id, pipe=None):

642 """Return the unique AlignTensorData instance corresponding to the alignment ID. 643 644 @param align_id: The alignment ID for the unique tensor. 645 @type align_id: str 646 @return: The alignment tensor object corresponding to the 'tensor' arg. 647 @rtype: AlignTensorData instance 648 """ 649 650 # The data pipe to check. 651 if pipe == None: 652 pipe = pipes.cdp_name() 653 654 # Init. 655 data = None 656 657 # Loop over the tensors. 658 count = 0 659 for i in range(len(cdp.align_tensors)): 660 if hasattr(cdp.align_tensors[i], 'align_id') and cdp.align_tensors[i].align_id == align_id: 661 data = cdp.align_tensors[i] 662 count += 1 663 664 # Multiple matches. 665 if count > 1: 666 raise RelaxError("Multiple alignment tensors match the alignment ID '%s'." % align_id) 667 # Return the object. 668 return data

669 670

671 -def init(tensor=None, align_id=None, params=None, scale=1.0, angle_units='deg', param_types=0, domain=None, errors=False):

672 """Function for initialising the alignment tensor. 673 674 @keyword tensor: The alignment tensor identification string. 675 @type tensor: str 676 @keyword align_id: The alignment ID string that the tensor corresponds to. 677 @type align_id: str or None 678 @keyword params: The alignment tensor parameters. 679 @type params: float 680 @keyword scale: The alignment tensor eigenvalue scaling value. 681 @type scale: float 682 @keyword angle_units: The units for the angle parameters (either 'deg' or 'rad'). 683 @type angle_units: str 684 @keyword param_types: The type of parameters supplied. The flag values correspond to, 0: {Axx, Ayy, Axy, Axz, Ayz}, and 1: {Azz, Axx-yy, Axy, Axz, Ayz}. 685 @type param_types: int 686 @keyword domain: The domain label. 687 @type domain: str or None 688 @keyword errors: A flag which determines if the alignment tensor data or its errors are being input. 689 @type errors: bool 690 """ 691 692 # Test if the current data pipe exists. 693 pipes.test() 694 695 # Parameter checks. 696 if align_id == None: 697 raise RelaxError("The alignment ID must be given.") 698 699 # Check the validity of the angle_units argument. 700 valid_types = ['deg', 'rad'] 701 if not angle_units in valid_types: 702 raise RelaxError("The alignment tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.") 703 704 # Test if alignment tensor data already exists. 705 if errors and (not hasattr(cdp, 'align_ids') or not align_id in cdp.align_ids): 706 raise RelaxNoTensorError('alignment') 707 708 # Check that the domain is defined. 709 if domain and (not hasattr(cdp, 'domain') or domain not in cdp.domain.keys()): 710 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % domain) 711 712 # Add the align ID to the current data pipe if needed. 713 if not hasattr(cdp, 'align_ids'): 714 cdp.align_ids = [] 715 if align_id not in cdp.align_ids: 716 cdp.align_ids.append(align_id) 717 718 # Add the align_tensors object to the data pipe. 719 tensor_obj = None 720 if not errors: 721 # Initialise the super structure. 722 if not hasattr(cdp, 'align_tensors'): 723 cdp.align_tensors = AlignTensorList() 724 725 # Add the tensor, if it doesn't already exist. 726 if tensor == None or tensor not in cdp.align_tensors.names(): 727 tensor_obj = cdp.align_tensors.add_item(tensor) 728 729 # Get the tensor. 730 if not tensor_obj: 731 if tensor: 732 tensor_obj = get_tensor_object(tensor) 733 else: 734 tensor_obj = get_tensor_object_from_align(align_id) 735 736 # {Sxx, Syy, Sxy, Sxz, Syz}. 737 if param_types == 0: 738 # Unpack the tuple. 739 Sxx, Syy, Sxy, Sxz, Syz = params 740 741 # Scaling. 742 Sxx = Sxx * scale 743 Syy = Syy * scale 744 Sxy = Sxy * scale 745 Sxz = Sxz * scale 746 Syz = Syz * scale 747 748 # Set the parameters. 749 set(tensor=tensor_obj, value=[Sxx, Syy, Sxy, Sxz, Syz], param=['Sxx', 'Syy', 'Sxy', 'Sxz', 'Syz'], errors=errors) 750 751 # {Szz, Sxx-yy, Sxy, Sxz, Syz}. 752 elif param_types == 1: 753 # Unpack the tuple. 754 Szz, Sxxyy, Sxy, Sxz, Syz = params 755 756 # Scaling. 757 Szz = Szz * scale 758 Sxxyy = Sxxyy * scale 759 Sxy = Sxy * scale 760 Sxz = Sxz * scale 761 Syz = Syz * scale 762 763 # Set the parameters. 764 set(tensor=tensor_obj, value=[Szz, Sxxyy, Sxy, Sxz, Syz], param=['Szz', 'Sxxyy', 'Sxy', 'Sxz', 'Syz'], errors=errors) 765 766 # {Axx, Ayy, Axy, Axz, Ayz}. 767 elif param_types == 2: 768 # Unpack the tuple. 769 Axx, Ayy, Axy, Axz, Ayz = params 770 771 # Scaling. 772 Axx = Axx * scale 773 Ayy = Ayy * scale 774 Axy = Axy * scale 775 Axz = Axz * scale 776 Ayz = Ayz * scale 777 778 # Set the parameters. 779 set(tensor=tensor_obj, value=[Axx, Ayy, Axy, Axz, Ayz], param=['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'], errors=errors) 780 781 # {Azz, Axx-yy, Axy, Axz, Ayz}. 782 elif param_types == 3: 783 # Unpack the tuple. 784 Azz, Axxyy, Axy, Axz, Ayz = params 785 786 # Scaling. 787 Azz = Azz * scale 788 Axxyy = Axxyy * scale 789 Axy = Axy * scale 790 Axz = Axz * scale 791 Ayz = Ayz * scale 792 793 # Set the parameters. 794 set(tensor=tensor_obj, value=[Azz, Axxyy, Axy, Axz, Ayz], param=['Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz'], errors=errors) 795 796 # {Axx, Ayy, Axy, Axz, Ayz}. 797 elif param_types == 4: 798 # Unpack the tuple. 799 Axx, Ayy, Axy, Axz, Ayz = params 800 801 # Get the bond length. 802 r = None 803 for spin in spin_loop(): 804 # First spin. 805 if r == None: 806 r = spin.r 807 808 # Different value. 809 if r != spin.r: 810 raise RelaxError("Not all spins have the same bond length.") 811 812 # Scaling. 813 scale = scale / kappa() * r**3 814 Axx = Axx * scale 815 Ayy = Ayy * scale 816 Axy = Axy * scale 817 Axz = Axz * scale 818 Ayz = Ayz * scale 819 820 # Set the parameters. 821 set(tensor=tensor_obj, value=[Axx, Ayy, Axy, Axz, Ayz], param=['Axx', 'Ayy', 'Axy', 'Axz', 'Ayz'], errors=errors) 822 823 # {Azz, Axx-yy, Axy, Axz, Ayz}. 824 elif param_types == 5: 825 # Unpack the tuple. 826 Azz, Axxyy, Axy, Axz, Ayz = params 827 828 # Get the bond length. 829 r = None 830 for spin in spin_loop(): 831 # First spin. 832 if r == None: 833 r = spin.r 834 835 # Different value. 836 if r != spin.r: 837 raise RelaxError("Not all spins have the same bond length.") 838 839 # Scaling. 840 scale = scale / kappa() * r**3 841 Azz = Azz * scale 842 Axxyy = Axxyy * scale 843 Axy = Axy * scale 844 Axz = Axz * scale 845 Ayz = Ayz * scale 846 847 # Set the parameters. 848 set(tensor=tensor_obj, value=[Azz, Axxyy, Axy, Axz, Ayz], param=['Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz'], errors=errors) 849 850 # {Pxx, Pyy, Pxy, Pxz, Pyz}. 851 elif param_types == 6: 852 # Unpack the tuple. 853 Pxx, Pyy, Pxy, Pxz, Pyz = params 854 855 # Scaling. 856 Pxx = Pxx * scale 857 Pyy = Pyy * scale 858 Pxy = Pxy * scale 859 Pxz = Pxz * scale 860 Pyz = Pyz * scale 861 862 # Set the parameters. 863 set(tensor=tensor_obj, value=[Pxx, Pyy, Pxy, Pxz, Pyz], param=['Pxx', 'Pyy', 'Pxy', 'Pxz', 'Pyz'], errors=errors) 864 865 # {Pzz, Pxx-yy, Pxy, Pxz, Pyz}. 866 elif param_types == 7: 867 # Unpack the tuple. 868 Pzz, Pxxyy, Pxy, Pxz, Pyz = params 869 870 # Scaling. 871 Pzz = Pzz * scale 872 Pxxyy = Pxxyy * scale 873 Pxy = Pxy * scale 874 Pxz = Pxz * scale 875 Pyz = Pyz * scale 876 877 # Set the parameters. 878 set(tensor=tensor_obj, value=[Pzz, Pxxyy, Pxy, Pxz, Pyz], param=['Pzz', 'Pxxyy', 'Pxy', 'Pxz', 'Pyz'], errors=errors) 879 880 # Unknown parameter combination. 881 else: 882 raise RelaxUnknownParamCombError('param_types', param_types) 883 884 # Set the domain and alignment ID. 885 if domain: 886 set_domain(tensor=tensor, domain=domain) 887 if align_id: 888 tensor_obj.set(param='align_id', value=align_id)

889 890

891 -def matrix_angles(basis_set=0, tensors=None):

892 """Function for calculating the 5D angles between the alignment tensors. 893 894 The basis set used for the 5D vector construction changes the angles calculated. 895 896 @param basis_set: The basis set to use for constructing the 5D vectors. If set to 0, the 897 basis set is {Sxx, Syy, Sxy, Sxz, Syz}. If 1, then the basis set is {Szz, 898 Sxxyy, Sxy, Sxz, Syz}. 899 @type basis_set: int 900 @param tensors: An array of tensors to apply SVD to. If None, all tensors will be used. 901 @type tensors: None or array of str 902 """ 903 904 # Test that alignment tensor data exists. 905 if not hasattr(cdp, 'align_tensors') or len(cdp.align_tensors) == 0: 906 raise RelaxNoTensorError('alignment') 907 908 # Count the number of tensors. 909 tensor_num = 0 910 for tensor in cdp.align_tensors: 911 if tensors and tensor.name not in tensors: 912 continue 913 tensor_num = tensor_num + 1 914 915 # Create the matrix which contains the 5D vectors. 916 matrix = zeros((tensor_num, 5), float64) 917 918 # Loop over the tensors. 919 i = 0 920 for tensor in cdp.align_tensors: 921 # Skip tensors. 922 if tensors and tensor.name not in tensors: 923 continue 924 925 # Unitary basis set. 926 if basis_set == 0: 927 # Pack the elements. 928 matrix[i, 0] = tensor.Sxx 929 matrix[i, 1] = tensor.Syy 930 matrix[i, 2] = tensor.Sxy 931 matrix[i, 3] = tensor.Sxz 932 matrix[i, 4] = tensor.Syz 933 934 # Geometric basis set. 935 elif basis_set == 1: 936 # Pack the elements. 937 matrix[i, 0] = tensor.Szz 938 matrix[i, 1] = tensor.Sxxyy 939 matrix[i, 2] = tensor.Sxy 940 matrix[i, 3] = tensor.Sxz 941 matrix[i, 4] = tensor.Syz 942 943 # Normalisation. 944 norm = linalg.norm(matrix[i]) 945 matrix[i] = matrix[i] / norm 946 947 # Increment the index. 948 i = i + 1 949 950 # Initialise the matrix of angles. 951 cdp.align_tensors.angles = zeros((tensor_num, tensor_num), float64) 952 953 # Header printout. 954 sys.stdout.write("\nData pipe: " + repr(pipes.cdp_name()) + "\n") 955 sys.stdout.write("\n5D angles in deg between the vectors ") 956 if basis_set == 0: 957 sys.stdout.write("{Sxx, Syy, Sxy, Sxz, Syz}") 958 elif basis_set == 1: 959 sys.stdout.write("{Szz, Sxx-yy, Sxy, Sxz, Syz}") 960 sys.stdout.write(":\n") 961 962 # Initialise the table of data. 963 table = [] 964 965 # The table header. 966 table.append(['']) 967 for i in range(tensor_num): 968 if cdp.align_tensors[i].name == None: 969 table[0].append(repr(i)) 970 else: 971 table[0].append(cdp.align_tensors[i].name) 972 973 # First loop over the rows. 974 for i in range(tensor_num): 975 # Add the tensor name. 976 if cdp.align_tensors[i].name == None: 977 table.append([repr(i)]) 978 else: 979 table.append([cdp.align_tensors[i].name]) 980 981 # Second loop over the columns. 982 for j in range(tensor_num): 983 # Dot product. 984 delta = dot(matrix[i], matrix[j]) 985 986 # Check. 987 if delta > 1: 988 delta = 1 989 990 # The angle (in rad). 991 cdp.align_tensors.angles[i, j] = arccos(delta) 992 993 # Add to the table as degrees. 994 table[i+1].append("%8.1f" % (cdp.align_tensors.angles[i, j]*180.0/pi)) 995 996 # Write out the table. 997 write_data(out=sys.stdout, data=table)

998 999

1000 -def num_tensors(skip_fixed=True):

1001 """Count the number of tensors. 1002 1003 @keyword skip_fixed: If set to True, then only the tensors without the fixed flag will be counted. If set to False, then all tensors will be counted. 1004 @type skip_fixed: bool 1005 @return: The number of tensors (excluding fixed tensors by default). 1006 @rtype: int 1007 """ 1008 1009 # Init. 1010 count = 0 1011 1012 # Loop over the tensors. 1013 for tensor_cont in cdp.align_tensors: 1014 # Skip fixed tensors. 1015 if skip_fixed and tensor_cont.fixed: 1016 continue 1017 1018 # Increment. 1019 count += 1 1020 1021 # Return the count. 1022 return count

1023 1024

1025 -def opt_uses_align_data(align_id=None):

1026 """Determine if the PCS or RDC data for the given alignment ID is needed for optimisation. 1027 1028 @keyword align_id: The optional alignment ID string. 1029 @type align_id: str 1030 @return: True if alignment data is to be used for optimisation, False otherwise. 1031 @rtype: bool 1032 """ 1033 1034 # No alignment IDs. 1035 if not hasattr(cdp, 'align_ids'): 1036 return False 1037 1038 # Convert the align IDs to an array, or take all IDs. 1039 if align_id: 1040 align_ids = [align_id] 1041 else: 1042 align_ids = cdp.align_ids 1043 1044 # Check the PCS and RDC. 1045 for align_id in align_ids: 1046 if pipe_control.pcs.opt_uses_pcs(align_id) or pipe_control.rdc.opt_uses_rdc(align_id): 1047 return True 1048 1049 # No alignment data is used for optimisation. 1050 return False

1051 1052

1053 -def opt_uses_tensor(tensor):

1054 """Determine if the given tensor is to be optimised. 1055 1056 @param tensor: The alignment tensor to check. 1057 @type tensor: AlignmentTensor object. 1058 @return: True if the tensor is to be optimised, False otherwise. 1059 @rtype: bool 1060 """ 1061 1062 # Combine all RDC and PCS IDs. 1063 ids = [] 1064 if hasattr(cdp, 'rdc_ids'): 1065 ids += cdp.rdc_ids 1066 if hasattr(cdp, 'pcs_ids'): 1067 ids += cdp.pcs_ids 1068 1069 # No RDC or PCS data for the alignment, so skip the tensor as it will not be optimised. 1070 if tensor.align_id not in ids: 1071 return False 1072 1073 # Fixed tensor. 1074 if tensor.fixed: 1075 return False 1076 1077 # The tensor is to be optimised. 1078 return True

1079 1080

1081 -def reduction(full_tensor=None, red_tensor=None):

1082 """Specify which tensor is a reduction of which other tensor. 1083 1084 @param full_tensor: The full alignment tensor. 1085 @type full_tensor: str 1086 @param red_tensor: The reduced alignment tensor. 1087 @type red_tensor: str 1088 """ 1089 1090 # Tensor information. 1091 match_full = False 1092 match_red = False 1093 i = 0 1094 for tensor_cont in cdp.align_tensors: 1095 # Test the tensor names. 1096 if tensor_cont.name == full_tensor: 1097 match_full = True 1098 index_full = i 1099 if tensor_cont.name == red_tensor: 1100 match_red = True 1101 index_red = i 1102 1103 # Increment. 1104 i = i + 1 1105 1106 # No match. 1107 if not match_full: 1108 raise RelaxNoTensorError('alignment', full_tensor) 1109 if not match_red: 1110 raise RelaxNoTensorError('alignment', red_tensor) 1111 1112 # Store. 1113 if not hasattr(cdp.align_tensors, 'reduction'): 1114 cdp.align_tensors.reduction = [] 1115 cdp.align_tensors.reduction.append([index_full, index_red])

1116 1117

1118 -def return_tensor(index, skip_fixed=True):

1119 """Return the tensor container for the given index, skipping fixed tensors if required. 1120 1121 @param index: The index of the tensor (if skip_fixed is True, then fixed tensors are not included in the index count). 1122 @type index: int 1123 @keyword skip_fixed: A flag which if True will exclude fixed tensors from the indexation. 1124 @type skip_fixed: bool 1125 @return: The tensor corresponding to the index. 1126 @rtype: data.align_tensor.AlignTensorData instance 1127 """ 1128 1129 # Init. 1130 count = 0 1131 1132 # Loop over the tensors. 1133 for tensor_cont in cdp.align_tensors: 1134 # Skip fixed tensors. 1135 if skip_fixed and tensor_cont.fixed: 1136 continue 1137 1138 # Index match, so return the container. 1139 if index == count: 1140 return tensor_cont 1141 1142 # Increment. 1143 count += 1 1144 1145 # Return False if the container was not found. 1146 return False

1147 1148

1149 -def set(tensor=None, value=None, param=None, errors=False):

1150 """Set the tensor. 1151 1152 @keyword tensor: The alignment tensor object. 1153 @type tensor: AlignTensorData instance 1154 @keyword value: The list of values to set the parameters to. 1155 @type value: list of float 1156 @keyword param: The list of parameter names. 1157 @type param: list of str 1158 @keyword errors: A flag which determines if the alignment tensor data or its errors are being 1159 input. 1160 @type errors: bool 1161 """ 1162 1163 # Initialise. 1164 geo_params = [] 1165 geo_values = [] 1166 orient_params = [] 1167 orient_values = [] 1168 1169 # Loop over the parameters. 1170 for i in range(len(param)): 1171 # Unknown parameter. 1172 if param[i] == None: 1173 raise RelaxUnknownParamError("alignment tensor", 'None') 1174 1175 # Default value. 1176 if value[i] == None: 1177 value[i] = 0.0 1178 1179 # Geometric parameter. 1180 if param[i] in ['Sxx', 'Syy', 'Szz', 'Sxxyy', 'Sxy', 'Sxz', 'Syz', 'Axx', 'Ayy', 'Azz', 'Axxyy', 'Axy', 'Axz', 'Ayz', 'Pxx', 'Pyy', 'Pzz', 'Pxxyy', 'Pxy', 'Pxz', 'Pyz']: 1181 geo_params.append(param[i]) 1182 geo_values.append(value[i]) 1183 1184 # Orientational parameter. 1185 if param[i] in ['alpha', 'beta', 'gamma']: 1186 orient_params.append(param[i]) 1187 orient_values.append(value[i]) 1188 1189 1190 # Geometric parameters. 1191 ####################### 1192 1193 # A single geometric parameter. 1194 if len(geo_params) == 1: 1195 # Saupe order matrix. 1196 ##################### 1197 1198 # The single parameter Sxx. 1199 if geo_params[0] == 'Sxx': 1200 if errors: 1201 tensor.set(param='Sxx', value=geo_values[0], category='err') 1202 else: 1203 tensor.set(param='Sxx', value=geo_values[0]) 1204 1205 # The single parameter Syy. 1206 elif geo_params[0] == 'Syy': 1207 if errors: 1208 tensor.set(param='Syy', value=geo_values[0], category='err') 1209 else: 1210 tensor.set(param='Syy', value=geo_values[0]) 1211 1212 # The single parameter Sxy. 1213 elif geo_params[0] == 'Sxy': 1214 if errors: 1215 tensor.set(param='Sxy', value=geo_values[0], category='err') 1216 else: 1217 tensor.set(param='Sxy', value=geo_values[0]) 1218 1219 # The single parameter Sxz. 1220 elif geo_params[0] == 'Sxz': 1221 if errors: 1222 tensor.set(param='Sxz', value=geo_values[0], category='err') 1223 else: 1224 tensor.set(param='Sxz', value=geo_values[0]) 1225 1226 # The single parameter Syz. 1227 elif geo_params[0] == 'Syz': 1228 if errors: 1229 tensor.set(param='Syz', value=geo_values[0], category='err') 1230 else: 1231 tensor.set(param='Syz', value=geo_values[0]) 1232 1233 1234 # Alignment tensor. 1235 ################### 1236 1237 # The single parameter Axx. 1238 elif geo_params[0] == 'Axx': 1239 if errors: 1240 tensor.set(param='Sxx', value=3.0/2.0 * geo_values[0], category='err') 1241 else: 1242 tensor.set(param='Sxx', value=3.0/2.0 * geo_values[0]) 1243 1244 # The single parameter Ayy. 1245 elif geo_params[0] == 'Ayy': 1246 if errors: 1247 tensor.set(param='Syy', value=3.0/2.0 * geo_values[0], category='err') 1248 else: 1249 tensor.set(param='Syy', value=3.0/2.0 * geo_values[0]) 1250 1251 # The single parameter Axy. 1252 elif geo_params[0] == 'Axy': 1253 if errors: 1254 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0], category='err') 1255 else: 1256 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0]) 1257 1258 # The single parameter Axz. 1259 elif geo_params[0] == 'Axz': 1260 if errors: 1261 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0], category='err') 1262 else: 1263 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0]) 1264 1265 # The single parameter Ayz. 1266 elif geo_params[0] == 'Ayz': 1267 if errors: 1268 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0], category='err') 1269 else: 1270 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0]) 1271 1272 1273 # Probability tensor. 1274 ##################### 1275 1276 # The single parameter Pxx. 1277 elif geo_params[0] == 'Pxx': 1278 if errors: 1279 tensor.set(param='Sxx', value=3.0/2.0 * (geo_values[0] - 1.0/3.0), category='err') 1280 else: 1281 tensor.set(param='Sxx', value=3.0/2.0 * (geo_values[0] - 1.0/3.0)) 1282 1283 # The single parameter Pyy. 1284 elif geo_params[0] == 'Pyy': 1285 if errors: 1286 tensor.set(param='Syy', value=3.0/2.0 * (geo_values[0] - 1.0/3.0), category='err') 1287 else: 1288 tensor.set(param='Syy', value=3.0/2.0 * (geo_values[0] - 1.0/3.0)) 1289 1290 # The single parameter Pxy. 1291 elif geo_params[0] == 'Pxy': 1292 if errors: 1293 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0], category='err') 1294 else: 1295 tensor.set(param='Sxy', value=3.0/2.0 * geo_values[0]) 1296 1297 # The single parameter Pxz. 1298 elif geo_params[0] == 'Pxz': 1299 if errors: 1300 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0], category='err') 1301 else: 1302 tensor.set(param='Sxz', value=3.0/2.0 * geo_values[0]) 1303 1304 # The single parameter Pyz. 1305 elif geo_params[0] == 'Pyz': 1306 if errors: 1307 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0], category='err') 1308 else: 1309 tensor.set(param='Syz', value=3.0/2.0 * geo_values[0]) 1310 1311 # Cannot set the single parameter. 1312 else: 1313 raise RelaxError("The geometric alignment parameter " + repr(geo_params[0]) + " cannot be set.") 1314 1315 # 5 geometric parameters. 1316 elif len(geo_params) == 5: 1317 # The geometric parameter set {Sxx, Syy, Sxy, Sxz, Syz}. 1318 if geo_params.count('Sxx') == 1 and geo_params.count('Syy') == 1 and geo_params.count('Sxy') == 1 and geo_params.count('Sxz') == 1 and geo_params.count('Syz') == 1: 1319 # The parameters. 1320 Sxx = geo_values[geo_params.index('Sxx')] 1321 Syy = geo_values[geo_params.index('Syy')] 1322 Sxy = geo_values[geo_params.index('Sxy')] 1323 Sxz = geo_values[geo_params.index('Sxz')] 1324 Syz = geo_values[geo_params.index('Syz')] 1325 1326 # Set the internal parameter values. 1327 if errors: 1328 tensor.set(param='Axx', value=2.0/3.0 * Sxx, category='err') 1329 tensor.set(param='Ayy', value=2.0/3.0 * Syy, category='err') 1330 tensor.set(param='Axy', value=2.0/3.0 * Sxy, category='err') 1331 tensor.set(param='Axz', value=2.0/3.0 * Sxz, category='err') 1332 tensor.set(param='Ayz', value=2.0/3.0 * Syz, category='err') 1333 else: 1334 tensor.set(param='Axx', value=2.0/3.0 * Sxx) 1335 tensor.set(param='Ayy', value=2.0/3.0 * Syy) 1336 tensor.set(param='Axy', value=2.0/3.0 * Sxy) 1337 tensor.set(param='Axz', value=2.0/3.0 * Sxz) 1338 tensor.set(param='Ayz', value=2.0/3.0 * Syz) 1339 1340 # The geometric parameter set {Szz, Sxxyy, Sxy, Sxz, Syz}. 1341 elif geo_params.count('Szz') == 1 and geo_params.count('Sxxyy') == 1 and geo_params.count('Sxy') == 1 and geo_params.count('Sxz') == 1 and geo_params.count('Syz') == 1: 1342 # The parameters. 1343 Szz = geo_values[geo_params.index('Szz')] 1344 Sxxyy = geo_values[geo_params.index('Sxxyy')] 1345 Sxy = geo_values[geo_params.index('Sxy')] 1346 Sxz = geo_values[geo_params.index('Sxz')] 1347 Syz = geo_values[geo_params.index('Syz')] 1348 1349 # Set the internal parameter values. 1350 if errors: 1351 tensor.set(param='Axx', value=2.0/3.0 * -0.5*(Szz-Sxxyy), category='err') 1352 tensor.set(param='Ayy', value=2.0/3.0 * -0.5*(Szz+Sxxyy), category='err') 1353 tensor.set(param='Axy', value=2.0/3.0 * Sxy, category='err') 1354 tensor.set(param='Axz', value=2.0/3.0 * Sxz, category='err') 1355 tensor.set(param='Ayz', value=2.0/3.0 * Syz, category='err') 1356 else: 1357 tensor.set(param='Axx', value=2.0/3.0 * -0.5*(Szz-Sxxyy)) 1358 tensor.set(param='Ayy', value=2.0/3.0 * -0.5*(Szz+Sxxyy)) 1359 tensor.set(param='Axy', value=2.0/3.0 * Sxy) 1360 tensor.set(param='Axz', value=2.0/3.0 * Sxz) 1361 tensor.set(param='Ayz', value=2.0/3.0 * Syz) 1362 1363 # The geometric parameter set {Axx, Ayy, Axy, Axz, Ayz}. 1364 elif geo_params.count('Axx') == 1 and geo_params.count('Ayy') == 1 and geo_params.count('Axy') == 1 and geo_params.count('Axz') == 1 and geo_params.count('Ayz') == 1: 1365 # The parameters. 1366 Axx = geo_values[geo_params.index('Axx')] 1367 Ayy = geo_values[geo_params.index('Ayy')] 1368 Axy = geo_values[geo_params.index('Axy')] 1369 Axz = geo_values[geo_params.index('Axz')] 1370 Ayz = geo_values[geo_params.index('Ayz')] 1371 1372 # Set the internal parameter values. 1373 if errors: 1374 tensor.set(param='Axx', value=Axx, category='err') 1375 tensor.set(param='Ayy', value=Ayy, category='err') 1376 tensor.set(param='Axy', value=Axy, category='err') 1377 tensor.set(param='Axz', value=Axz, category='err') 1378 tensor.set(param='Ayz', value=Ayz, category='err') 1379 else: 1380 tensor.set(param='Axx', value=Axx) 1381 tensor.set(param='Ayy', value=Ayy) 1382 tensor.set(param='Axy', value=Axy) 1383 tensor.set(param='Axz', value=Axz) 1384 tensor.set(param='Ayz', value=Ayz) 1385 1386 # The geometric parameter set {Azz, Axxyy, Axy, Axz, Ayz}. 1387 elif geo_params.count('Azz') == 1 and geo_params.count('Axxyy') == 1 and geo_params.count('Axy') == 1 and geo_params.count('Axz') == 1 and geo_params.count('Ayz') == 1: 1388 # The parameters. 1389 Azz = geo_values[geo_params.index('Azz')] 1390 Axxyy = geo_values[geo_params.index('Axxyy')] 1391 Axy = geo_values[geo_params.index('Axy')] 1392 Axz = geo_values[geo_params.index('Axz')] 1393 Ayz = geo_values[geo_params.index('Ayz')] 1394 1395 # Set the internal parameter values. 1396 if errors: 1397 tensor.set(param='Axx', value=-0.5*(Azz-Axxyy), category='err') 1398 tensor.set(param='Ayy', value=-0.5*(Azz+Axxyy), category='err') 1399 tensor.set(param='Axy', value=Axy, category='err') 1400 tensor.set(param='Axz', value=Axz, category='err') 1401 tensor.set(param='Ayz', value=Ayz, category='err') 1402 else: 1403 tensor.set(param='Axx', value=-0.5*(Azz-Axxyy)) 1404 tensor.set(param='Ayy', value=-0.5*(Azz+Axxyy)) 1405 tensor.set(param='Axy', value=Axy) 1406 tensor.set(param='Axz', value=Axz) 1407 tensor.set(param='Ayz', value=Ayz) 1408 1409 # The geometric parameter set {Pxx, Pyy, Pxy, Pxz, Pyz}. 1410 elif geo_params.count('Pxx') == 1 and geo_params.count('Pyy') == 1 and geo_params.count('Pxy') == 1 and geo_params.count('Pxz') == 1 and geo_params.count('Pyz') == 1: 1411 # The parameters. 1412 Pxx = geo_values[geo_params.index('Pxx')] 1413 Pyy = geo_values[geo_params.index('Pyy')] 1414 Pxy = geo_values[geo_params.index('Pxy')] 1415 Pxz = geo_values[geo_params.index('Pxz')] 1416 Pyz = geo_values[geo_params.index('Pyz')] 1417 1418 # Set the internal parameter values. 1419 if errors: 1420 tensor.set(param='Axx', value=Pxx - 1.0/3.0, category='err') 1421 tensor.set(param='Ayy', value=Pyy - 1.0/3.0, category='err') 1422 tensor.set(param='Axy', value=Pxy, category='err') 1423 tensor.set(param='Axz', value=Pxz, category='err') 1424 tensor.set(param='Ayz', value=Pyz, category='err') 1425 else: 1426 tensor.set(param='Axx', value=Pxx - 1.0/3.0) 1427 tensor.set(param='Ayy', value=Pyy - 1.0/3.0) 1428 tensor.set(param='Axy', value=Pxy) 1429 tensor.set(param='Axz', value=Pxz) 1430 tensor.set(param='Ayz', value=Pyz) 1431 1432 # The geometric parameter set {Pzz, Pxxyy, Pxy, Pxz, Pyz}. 1433 elif geo_params.count('Pzz') == 1 and geo_params.count('Pxxyy') == 1 and geo_params.count('Pxy') == 1 and geo_params.count('Pxz') == 1 and geo_params.count('Pyz') == 1: 1434 # The parameters. 1435 Pzz = geo_values[geo_params.index('Pzz')] 1436 Pxxyy = geo_values[geo_params.index('Pxxyy')] 1437 Pxy = geo_values[geo_params.index('Pxy')] 1438 Pxz = geo_values[geo_params.index('Pxz')] 1439 Pyz = geo_values[geo_params.index('Pyz')] 1440 1441 # Set the internal parameter values. 1442 if errors: 1443 tensor.set(param='Axx', value=-0.5*(Pzz-Pxxyy) - 1.0/3.0, category='err') 1444 tensor.set(param='Ayy', value=-0.5*(Pzz+Pxxyy) - 1.0/3.0, category='err') 1445 tensor.set(param='Axy', value=Pxy, category='err') 1446 tensor.set(param='Axz', value=Pxz, category='err') 1447 tensor.set(param='Ayz', value=Pyz, category='err') 1448 else: 1449 tensor.set(param='Axx', value=-0.5*(Pzz-Pxxyy) - 1.0/3.0) 1450 tensor.set(param='Ayy', value=-0.5*(Pzz+Pxxyy) - 1.0/3.0) 1451 tensor.set(param='Axy', value=Pxy) 1452 tensor.set(param='Axz', value=Pxz) 1453 tensor.set(param='Ayz', value=Pyz) 1454 1455 # Unknown parameter combination. 1456 else: 1457 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1458 1459 1460 # Unknown geometric parameters. 1461 else: 1462 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1463 1464 1465 # Orientational parameters. 1466 ########################### 1467 1468 # A single orientational parameter. 1469 if len(orient_params) == 1: 1470 # The single parameter alpha. 1471 if orient_params[0] == 'alpha': 1472 if errors: 1473 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err') 1474 else: 1475 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1476 1477 # The single parameter beta. 1478 elif orient_params[0] == 'beta': 1479 if errors: 1480 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err') 1481 else: 1482 tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1483 1484 # The single parameter gamma. 1485 elif orient_params[0] == 'gamma': 1486 if errors: 1487 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err') 1488 else: 1489 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1490 1491 # Two orientational parameters. 1492 elif len(orient_params) == 2: 1493 # The orientational parameter set {alpha, beta}. 1494 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1495 if errors: 1496 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err') 1497 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err') 1498 else: 1499 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1500 tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1501 1502 # The orientational parameter set {alpha, gamma}. 1503 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1: 1504 if errors: 1505 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err') 1506 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err') 1507 else: 1508 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1509 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1510 1511 # The orientational parameter set {beta, gamma}. 1512 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1: 1513 if errors: 1514 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err') 1515 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err') 1516 else: 1517 tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1518 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1519 1520 # Unknown parameter combination. 1521 else: 1522 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1523 1524 # Three orientational parameters. 1525 elif len(orient_params) == 3: 1526 # The orientational parameter set {alpha, beta, gamma}. 1527 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1528 if errors: 1529 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')], category='err') 1530 tensor.set(param='beta', value=orient_values[orient_params.index('beta')], category='err') 1531 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')], category='err') 1532 else: 1533 tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1534 tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1535 tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1536 1537 # Unknown parameter combination. 1538 else: 1539 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1540 1541 # More than three orientational parameters. 1542 elif len(orient_params) > 3: 1543 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1544 1545 1546 # Fold the angles in. 1547 ##################### 1548 1549 if orient_params: 1550 fold_angles()

1551 1552

1553 -def set_align_id(tensor=None, align_id=None):

1554 """Set the align ID string for the given tensor. 1555 1556 @keyword tensor: The alignment tensor label. 1557 @type tensor: str 1558 @keyword align_id: The alignment ID string. 1559 @type align_id: str 1560 """ 1561 1562 # Loop over the tensors. 1563 match = False 1564 for tensor_cont in cdp.align_tensors: 1565 # Find the matching tensor and then store the ID. 1566 if tensor_cont.name == tensor: 1567 tensor_cont.align_id = align_id 1568 match = True 1569 1570 # The tensor label doesn't exist. 1571 if not match: 1572 raise RelaxNoTensorError('alignment', tensor)

1573 1574

1575 -def set_domain(tensor=None, domain=None):

1576 """Set the domain label for the given tensor. 1577 1578 @param tensor: The alignment tensor label. 1579 @type tensor: str 1580 @param domain: The domain label. 1581 @type domain: str 1582 """ 1583 1584 # Check that the domain is defined. 1585 if not hasattr(cdp, 'domain') or domain not in cdp.domain.keys(): 1586 raise RelaxError("The domain '%s' has not been defined. Please use the domain user function." % domain) 1587 1588 # Loop over the tensors. 1589 match = False 1590 for tensor_cont in cdp.align_tensors: 1591 # Find the matching tensor and then store the domain label. 1592 if tensor_cont.name == tensor: 1593 tensor_cont.set(param='domain', value=domain) 1594 match = True 1595 1596 # The tensor label doesn't exist. 1597 if not match: 1598 raise RelaxNoTensorError('alignment', tensor)

1599 1600

1601 -def svd(basis_set=0, tensors=None):

1602 """Function for calculating the singular values of all the loaded tensors. 1603 1604 The matrix on which SVD will be performed is:: 1605 1606 | Sxx1 Syy1 Sxy1 Sxz1 Syz1 | 1607 | Sxx2 Syy2 Sxy2 Sxz2 Syz2 | 1608 | Sxx3 Syy3 Sxy3 Sxz3 Syz3 | 1609 | . . . . . | 1610 | . . . . . | 1611 | . . . . . | 1612 | SxxN SyyN SxyN SxzN SyzN | 1613 1614 This is the default unitary basis set (selected when basis_set is 0). Alternatively a geometric 1615 basis set consisting of the stretching and skewing parameters Szz and Sxx-yy respectively 1616 replacing Sxx and Syy can be chosen by setting basis_set to 1. The matrix in this case is:: 1617 1618 | Szz1 Sxxyy1 Sxy1 Sxz1 Syz1 | 1619 | Szz2 Sxxyy2 Sxy2 Sxz2 Syz2 | 1620 | Szz3 Sxxyy3 Sxy3 Sxz3 Syz3 | 1621 | . . . . . | 1622 | . . . . . | 1623 | . . . . . | 1624 | SzzN SxxyyN SxyN SxzN SyzN | 1625 1626 The relationships between the geometric and unitary basis sets are:: 1627 1628 Szz = - Sxx - Syy, 1629 Sxxyy = Sxx - Syy, 1630 1631 The SVD values and condition number are dependendent upon the basis set chosen. 1632 1633 1634 @param basis_set: The basis set to create the 5 by n matrix on which to perform SVD. 1635 @type basis_set: int 1636 @param tensors: An array of tensors to apply SVD to. If None, all tensors will be used. 1637 @type tensors: None or array of str 1638 """ 1639 1640 # Test that alignment tensor data exists. 1641 if not hasattr(cdp, 'align_tensors') or len(cdp.align_tensors) == 0: 1642 raise RelaxNoTensorError('alignment') 1643 1644 # Count the number of tensors used in the SVD. 1645 tensor_num = 0 1646 for tensor in cdp.align_tensors: 1647 if tensors and tensor.name not in tensors: 1648 continue 1649 tensor_num = tensor_num + 1 1650 1651 # Create the matrix to apply SVD on. 1652 matrix = zeros((tensor_num, 5), float64) 1653 1654 # Pack the elements. 1655 i = 0 1656 for tensor in cdp.align_tensors: 1657 # Skip tensors. 1658 if tensors and tensor.name not in tensors: 1659 continue 1660 1661 # Unitary basis set. 1662 if basis_set == 0: 1663 matrix[i, 0] = tensor.Sxx 1664 matrix[i, 1] = tensor.Syy 1665 matrix[i, 2] = tensor.Sxy 1666 matrix[i, 3] = tensor.Sxz 1667 matrix[i, 4] = tensor.Syz 1668 1669 # Geometric basis set. 1670 elif basis_set == 1: 1671 matrix[i, 0] = tensor.Szz 1672 matrix[i, 1] = tensor.Sxxyy 1673 matrix[i, 2] = tensor.Sxy 1674 matrix[i, 3] = tensor.Sxz 1675 matrix[i, 4] = tensor.Syz 1676 1677 # Increment the index. 1678 i = i + 1 1679 1680 # SVD. 1681 u, s, vh = linalg.svd(matrix) 1682 1683 # Store the singular values. 1684 cdp.align_tensors.singular_vals = s 1685 1686 # Calculate and store the condition number. 1687 cdp.align_tensors.cond_num = s[0] / s[-1] 1688 1689 # Print out. 1690 print("\nData pipe: " + repr(pipes.cdp_name())) 1691 print("\nSingular values:") 1692 for val in s: 1693 print(" %.4e" % val) 1694 print("\nCondition number: %.2f" % cdp.align_tensors.cond_num)

1695

Source Code for Module pipe_control.align_tensor