generic_fns.diffusion

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2012 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 support of diffusion tensors.""" 24 25 # Python module imports. 26 from copy import deepcopy 27 from math import cos, pi, sin 28 from numpy import cross, float64, int32, ones, transpose, zeros 29 from numpy.linalg import norm, svd 30 from re import search 31 import string 32 33 # relax module imports. 34 from data.diff_tensor import DiffTensorData 35 from generic_fns import pipes 36 from generic_fns.angles import fold_spherical_angles, wrap_angles 37 from generic_fns.interatomic import return_interatom_list 38 from generic_fns.mol_res_spin import get_molecule_names, return_spin, spin_loop 39 from maths_fns.coord_transform import cartesian_to_spherical 40 from maths_fns.rotation_matrix import R_to_euler_zyz 41 from physical_constants import element_from_isotope, number_from_isotope 42 from relax_errors import RelaxError, RelaxNoTensorError, RelaxStrError, RelaxTensorError, RelaxUnknownParamCombError, RelaxUnknownParamError 43 from user_functions.data import Uf_tables; uf_tables = Uf_tables() 44 from user_functions.objects import Desc_container 45 46

47 -def bmrb_read(star):

48 """Read the relaxation data from the NMR-STAR dictionary object. 49 50 @param star: The NMR-STAR dictionary object. 51 @type star: NMR_STAR instance 52 """ 53 54 # Get the diffusion tensor data. 55 found = 0 56 for data in star.tensor.loop(): 57 # No data. 58 if data == None: 59 continue 60 61 # Not a diffusion tensor. 62 if data['tensor_type'] != 'diffusion': 63 continue 64 65 # Found. 66 found = found + 1 67 68 # No diffusion tensor data. 69 if not found: 70 return 71 72 # Check. 73 if found != 1: 74 raise RelaxError("More than one diffusion tensor found.") 75 76 # Rebuild the tensor. 77 tensor = zeros((3, 3), float64) 78 tensor[0, 0] = data['tensor_11'][0] 79 tensor[0, 1] = data['tensor_12'][0] 80 tensor[0, 2] = data['tensor_13'][0] 81 tensor[1, 0] = data['tensor_21'][0] 82 tensor[1, 1] = data['tensor_22'][0] 83 tensor[1, 2] = data['tensor_23'][0] 84 tensor[2, 0] = data['tensor_31'][0] 85 tensor[2, 1] = data['tensor_32'][0] 86 tensor[2, 2] = data['tensor_33'][0] 87 88 # Eigenvalues. 89 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 90 91 # X-Y eigenvalue comparison. 92 xy_match = False 93 epsilon = 1e-1 94 if abs(Di[0] - Di[1]) < epsilon: 95 xy_match = True 96 97 # Y-Z eigenvalue comparison. 98 yz_match = False 99 if abs(Di[1] - Di[2]) < epsilon: 100 yz_match = True 101 102 # Determine the tensor type. 103 if xy_match and yz_match: 104 shape = ['sphere'] 105 elif xy_match: 106 shape = ['spheroid', 'prolate spheroid'] 107 type = 'prolate' 108 Dpar = Di[2] 109 Dper = Di[0] 110 elif yz_match: 111 shape = ['spheroid', 'oblate spheroid'] 112 type = 'oblate' 113 Dpar = Di[0] 114 Dper = Di[2] 115 else: 116 shape = ['ellipsoid'] 117 118 # Check the shape. 119 if data['geometric_shape'] not in shape: 120 raise RelaxError("The tensor with eigenvalues %s does not match the %s geometric shape." % (Di, shape[0])) 121 122 # Add the diff_tensor object to the data pipe. 123 cdp.diff_tensor = DiffTensorData() 124 125 # Set the fixed flag. 126 cdp.diff_tensor.set_fixed(True) 127 128 # Sphere. 129 if data['geometric_shape'] == 'sphere': 130 sphere(params=Di[0], d_scale=1.0, param_types=1) 131 132 # Spheroid. 133 elif data['geometric_shape'] in ['spheroid', 'oblate spheroid', 'prolate spheroid']: 134 # The spherical angles. 135 theta = data['axial_sym_axis_polar_angle'][0] 136 phi = data['axial_sym_axis_azimuthal_angle'][0] 137 138 # Set up the tensor. 139 spheroid(params=(Dpar, Dper, theta, phi), d_scale=1.0, param_types=3, spheroid_type=type) 140 141 # Ellipsoid. 142 elif data['geometric_shape'] == 'ellipsoid': 143 ellipsoid(params=(Di[0], Di[1], Di[2], alpha, beta, gamma), d_scale=1.0, param_types=3)

144 145

146 -def bmrb_write(star):

147 """Generate the diffusion tensor saveframes for the NMR-STAR dictionary object. 148 149 @param star: The NMR-STAR dictionary object. 150 @type star: NMR_STAR instance 151 """ 152 153 # Get the current data pipe. 154 cdp = pipes.get_pipe() 155 156 # Initialise the spin specific data lists. 157 mol_name_list = [] 158 res_num_list = [] 159 res_name_list = [] 160 atom_name_list = [] 161 isotope_list = [] 162 element_list = [] 163 attached_atom_name_list = [] 164 attached_isotope_list = [] 165 attached_element_list = [] 166 167 # Store the spin specific data in lists for later use. 168 for spin, mol_name, res_num, res_name, spin_id in spin_loop(full_info=True, return_id=True): 169 # Check the data for None (not allowed in BMRB!). 170 if res_num == None: 171 raise RelaxError("For the BMRB, the residue of spin '%s' must be numbered." % spin_id) 172 if res_name == None: 173 raise RelaxError("For the BMRB, the residue of spin '%s' must be named." % spin_id) 174 if spin.name == None: 175 raise RelaxError("For the BMRB, the spin '%s' must be named." % spin_id) 176 if not hasattr(spin, 'isotope') or spin.isotope == None: 177 raise RelaxError("For the BMRB, the spin isotope type of '%s' must be specified." % spin_id) 178 179 # The molecule/residue/spin info. 180 mol_name_list.append(mol_name) 181 res_num_list.append(str(res_num)) 182 res_name_list.append(str(res_name)) 183 atom_name_list.append(str(spin.name)) 184 185 # Interatomic info. 186 interatoms = return_interatom_list(spin_id) 187 if len(interatoms) == 0: 188 raise RelaxError("No interatomic interactions are defined for the spin '%s'." % spin_id) 189 if len(interatoms) > 1: 190 raise RelaxError("The BMRB only handles a signal interatomic interaction for the spin '%s'." % spin_id) 191 192 # Get the attached spin. 193 spin_attached = return_spin(interatoms[0].spin_id1) 194 if id(spin_attached) == id(spin): 195 spin_attached = return_spin(interatoms[0].spin_id2) 196 197 # The attached atom info. 198 if hasattr(spin_attached, 'name'): 199 attached_atom_name_list.append(str(spin_attached.name)) 200 else: 201 attached_atom_name_list.append(None) 202 if hasattr(spin_attached, 'isotope'): 203 attached_element_list.append(element_from_isotope(spin_attached.isotope)) 204 attached_isotope_list.append(str(number_from_isotope(spin_attached.isotope))) 205 else: 206 attached_element_list.append(None) 207 attached_isotope_list.append(None) 208 209 # Other info. 210 isotope_list.append(int(spin.isotope.strip(string.ascii_letters))) 211 element_list.append(spin.element) 212 213 # Convert the molecule names into the entity IDs. 214 entity_ids = zeros(len(mol_name_list), int32) 215 mol_names = get_molecule_names() 216 for i in range(len(mol_name_list)): 217 for j in range(len(mol_names)): 218 if mol_name_list[i] == mol_names[j]: 219 entity_ids[i] = j+1 220 221 # The tensor geometric shape. 222 geometric_shape = cdp.diff_tensor.type 223 if geometric_shape == 'spheroid': 224 geometric_shape = "%s %s" % (cdp.diff_tensor.spheroid_type, geometric_shape) 225 226 # The tensor symmetry. 227 shapes = ['sphere', 'oblate spheroid', 'prolate spheroid', 'ellipsoid'] 228 sym = ['isotropic', 'axial symmetry', 'axial symmetry', 'rhombic'] 229 for i in range(len(shapes)): 230 if geometric_shape == shapes[i]: 231 tensor_symmetry = sym[i] 232 233 # Axial symmetry axis. 234 theta = None 235 phi = None 236 if tensor_symmetry == 'axial symmetry': 237 theta = cdp.diff_tensor.theta 238 phi = cdp.diff_tensor.phi 239 240 # Euler angles. 241 alpha, beta, gamma = None, None, None 242 if tensor_symmetry == 'rhombic': 243 alpha = cdp.diff_tensor.alpha 244 beta = cdp.diff_tensor.beta 245 gamma = cdp.diff_tensor.gamma 246 247 # The tensor eigenvalues. 248 Diso = cdp.diff_tensor.Diso 249 Da = None 250 Dr = None 251 if tensor_symmetry == 'axial symmetry': 252 Da = cdp.diff_tensor.Da 253 elif tensor_symmetry == 'rhombic': 254 Dr = cdp.diff_tensor.Dr 255 256 # The full tensor. 257 tensor_11 = cdp.diff_tensor.tensor[0, 0] 258 tensor_12 = cdp.diff_tensor.tensor[0, 1] 259 tensor_13 = cdp.diff_tensor.tensor[0, 2] 260 tensor_21 = cdp.diff_tensor.tensor[1, 0] 261 tensor_22 = cdp.diff_tensor.tensor[1, 1] 262 tensor_23 = cdp.diff_tensor.tensor[1, 2] 263 tensor_31 = cdp.diff_tensor.tensor[2, 0] 264 tensor_32 = cdp.diff_tensor.tensor[2, 1] 265 tensor_33 = cdp.diff_tensor.tensor[2, 2] 266 267 268 # Add the diffusion tensor. 269 star.tensor.add(tensor_type='diffusion', euler_type='zyz', geometric_shape=geometric_shape, tensor_symmetry=tensor_symmetry, matrix_val_units='s-1', angle_units='rad', iso_val_formula='Diso = 1/(6.tm)', aniso_val_formula='Da = Dpar - Dper', rhomb_val_formula='Dr = (Dy - Dx)/2Da', entity_ids=entity_ids, res_nums=res_num_list, res_names=res_name_list, atom_names=atom_name_list, atom_types=element_list, isotope=isotope_list, axial_sym_axis_polar_angle=theta, axial_sym_axis_azimuthal_angle=phi, iso_val=Diso, aniso_val=Da, rhombic_val=Dr, euler_alpha=alpha, euler_beta=beta, euler_gamma=gamma, tensor_11=tensor_11, tensor_12=tensor_12, tensor_13=tensor_13, tensor_21=tensor_21, tensor_22=tensor_22, tensor_23=tensor_23, tensor_31=tensor_31, tensor_32=tensor_32, tensor_33=tensor_33)

270 271 272

273 -def copy(pipe_from=None, pipe_to=None):

274 """Function for copying diffusion tensor data from one data pipe to another. 275 276 @param pipe_from: The data pipe to copy the diffusion tensor data from. This defaults to the 277 current data pipe. 278 @type pipe_from: str 279 @param pipe_to: The data pipe to copy the diffusion tensor data to. This defaults to the 280 current data pipe. 281 @type pipe_to: str 282 """ 283 284 # Defaults. 285 if pipe_from == None and pipe_to == None: 286 raise RelaxError("The pipe_from and pipe_to arguments cannot both be set to None.") 287 elif pipe_from == None: 288 pipe_from = pipes.cdp_name() 289 elif pipe_to == None: 290 pipe_to = pipes.cdp_name() 291 292 # Test if the pipe_from and pipe_to data pipes exist. 293 pipes.test(pipe_from) 294 pipes.test(pipe_to) 295 296 # Get the data pipes. 297 dp_from = pipes.get_pipe(pipe_from) 298 dp_to = pipes.get_pipe(pipe_to) 299 300 # Test if pipe_from contains diffusion tensor data. 301 if not diff_data_exists(pipe_from): 302 raise RelaxNoTensorError('diffusion') 303 304 # Test if pipe_to contains diffusion tensor data. 305 if diff_data_exists(pipe_to): 306 raise RelaxTensorError('diffusion') 307 308 # Copy the data. 309 dp_to.diff_tensor = deepcopy(dp_from.diff_tensor)

310 311

312 -def data_names():

313 """Function for returning a list of names of data structures associated with the sequence.""" 314 315 names = [ 'diff_type', 316 'diff_params' ] 317 318 return names

319 320

321 -def default_value(param):

322 """Return the default values for the diffusion tensor parameters. 323 324 @param param: The name of the parameter. 325 @type param: str 326 @return: The default value. 327 @rtype: float 328 """ 329 330 # tm. 331 if param == 'tm': 332 return 10.0 * 1e-9 333 334 # Diso, Dx, Dy, Dz, Dpar, Dper. 335 elif param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 336 return 1.666 * 1e7 337 338 # Da, Dr. 339 elif param == 'Da' or param == 'Dr': 340 return 0.0 341 342 # Dratio. 343 elif param == 'Dratio': 344 return 1.0 345 346 # All angles. 347 elif param == 'alpha' or param == 'beta' or param == 'gamma' or param == 'theta' or param == 'phi': 348 return 0.0

349 350 # User function documentation. 351 __default_value_doc__ = Desc_container("Diffusion tensor parameter default values") 352 table = uf_tables.add_table(label="table: diff default values", caption="Diffusion tensor parameter default values.") 353 table.add_headings(["Data type", "Object name", "Value"]) 354 table.add_row(["tm", "'tm'", "10 * 1e-9"]) 355 table.add_row(["Diso", "'Diso'", "1.666 * 1e7"]) 356 table.add_row(["Da", "'Da'", "0.0"]) 357 table.add_row(["Dr", "'Dr'", "0.0"]) 358 table.add_row(["Dx", "'Dx'", "1.666 * 1e7"]) 359 table.add_row(["Dy", "'Dy'", "1.666 * 1e7"]) 360 table.add_row(["Dz", "'Dz'", "1.666 * 1e7"]) 361 table.add_row(["Dpar", "'Dpar'", "1.666 * 1e7"]) 362 table.add_row(["Dper", "'Dper'", "1.666 * 1e7"]) 363 table.add_row(["Dratio", "'Dratio'", "1.0"]) 364 table.add_row(["alpha", "'alpha'", "0.0"]) 365 table.add_row(["beta", "'beta'", "0.0"]) 366 table.add_row(["gamma", "'gamma'", "0.0"]) 367 table.add_row(["theta", "'theta'", "0.0"]) 368 table.add_row(["phi", "'phi'", "0.0"]) 369 __default_value_doc__.add_table(table.label) 370 371

372 -def delete():

373 """Function for deleting diffusion tensor data.""" 374 375 # Test if the current data pipe exists. 376 pipes.test() 377 378 # Test if diffusion tensor data exists. 379 if not diff_data_exists(): 380 raise RelaxNoTensorError('diffusion') 381 382 # Delete the diffusion data. 383 del(cdp.diff_tensor)

384 385

386 -def diff_data_exists(pipe=None):

387 """Function for determining if diffusion data exists in the current data pipe. 388 389 @param pipe: The data pipe to search for data in. 390 @type pipe: str 391 @return: The answer to the question. 392 @rtype: bool 393 """ 394 395 # The data pipe to check. 396 if pipe == None: 397 pipe = pipes.cdp_name() 398 399 # Get the data pipe. 400 dp = pipes.get_pipe(pipe) 401 402 # Test if the data structure exists. 403 if hasattr(dp, 'diff_tensor'): 404 return True 405 else: 406 return False

407 408

409 -def display():

410 """Function for displaying the diffusion tensor.""" 411 412 # Test if the current data pipe exists. 413 pipes.test() 414 415 # Test if diffusion tensor data exists. 416 if not diff_data_exists(): 417 raise RelaxNoTensorError('diffusion') 418 419 # Spherical diffusion. 420 if cdp.diff_tensor.type == 'sphere': 421 # Tensor type. 422 print("Type: Spherical diffusion") 423 424 # Parameters. 425 print("\nParameters {tm}.") 426 print("tm (s): " + repr(cdp.diff_tensor.tm)) 427 428 # Alternate parameters. 429 print("\nAlternate parameters {Diso}.") 430 print("Diso (1/s): " + repr(cdp.diff_tensor.Diso)) 431 432 # Fixed flag. 433 print("\nFixed: " + repr(cdp.diff_tensor.fixed)) 434 435 # Spheroidal diffusion. 436 elif cdp.diff_tensor.type == 'spheroid': 437 # Tensor type. 438 print("Type: Spheroidal diffusion") 439 440 # Parameters. 441 print("\nParameters {tm, Da, theta, phi}.") 442 print("tm (s): " + repr(cdp.diff_tensor.tm)) 443 print("Da (1/s): " + repr(cdp.diff_tensor.Da)) 444 print("theta (rad): " + repr(cdp.diff_tensor.theta)) 445 print("phi (rad): " + repr(cdp.diff_tensor.phi)) 446 447 # Alternate parameters. 448 print("\nAlternate parameters {Diso, Da, theta, phi}.") 449 print("Diso (1/s): " + repr(cdp.diff_tensor.Diso)) 450 print("Da (1/s): " + repr(cdp.diff_tensor.Da)) 451 print("theta (rad): " + repr(cdp.diff_tensor.theta)) 452 print("phi (rad): " + repr(cdp.diff_tensor.phi)) 453 454 # Alternate parameters. 455 print("\nAlternate parameters {Dpar, Dper, theta, phi}.") 456 print("Dpar (1/s): " + repr(cdp.diff_tensor.Dpar)) 457 print("Dper (1/s): " + repr(cdp.diff_tensor.Dper)) 458 print("theta (rad): " + repr(cdp.diff_tensor.theta)) 459 print("phi (rad): " + repr(cdp.diff_tensor.phi)) 460 461 # Alternate parameters. 462 print("\nAlternate parameters {tm, Dratio, theta, phi}.") 463 print("tm (s): " + repr(cdp.diff_tensor.tm)) 464 print("Dratio: " + repr(cdp.diff_tensor.Dratio)) 465 print("theta (rad): " + repr(cdp.diff_tensor.theta)) 466 print("phi (rad): " + repr(cdp.diff_tensor.phi)) 467 468 # Fixed flag. 469 print("\nFixed: " + repr(cdp.diff_tensor.fixed)) 470 471 # Ellipsoidal diffusion. 472 elif cdp.diff_tensor.type == 'ellipsoid': 473 # Tensor type. 474 print("Type: Ellipsoidal diffusion") 475 476 # Parameters. 477 print("\nParameters {tm, Da, Dr, alpha, beta, gamma}.") 478 print("tm (s): " + repr(cdp.diff_tensor.tm)) 479 print("Da (1/s): " + repr(cdp.diff_tensor.Da)) 480 print("Dr: " + repr(cdp.diff_tensor.Dr)) 481 print("alpha (rad): " + repr(cdp.diff_tensor.alpha)) 482 print("beta (rad): " + repr(cdp.diff_tensor.beta)) 483 print("gamma (rad): " + repr(cdp.diff_tensor.gamma)) 484 485 # Alternate parameters. 486 print("\nAlternate parameters {Diso, Da, Dr, alpha, beta, gamma}.") 487 print("Diso (1/s): " + repr(cdp.diff_tensor.Diso)) 488 print("Da (1/s): " + repr(cdp.diff_tensor.Da)) 489 print("Dr: " + repr(cdp.diff_tensor.Dr)) 490 print("alpha (rad): " + repr(cdp.diff_tensor.alpha)) 491 print("beta (rad): " + repr(cdp.diff_tensor.beta)) 492 print("gamma (rad): " + repr(cdp.diff_tensor.gamma)) 493 494 # Alternate parameters. 495 print("\nAlternate parameters {Dx, Dy, Dz, alpha, beta, gamma}.") 496 print("Dx (1/s): " + repr(cdp.diff_tensor.Dx)) 497 print("Dy (1/s): " + repr(cdp.diff_tensor.Dy)) 498 print("Dz (1/s): " + repr(cdp.diff_tensor.Dz)) 499 print("alpha (rad): " + repr(cdp.diff_tensor.alpha)) 500 print("beta (rad): " + repr(cdp.diff_tensor.beta)) 501 print("gamma (rad): " + repr(cdp.diff_tensor.gamma)) 502 503 # Fixed flag. 504 print("\nFixed: " + repr(cdp.diff_tensor.fixed))

505 506

507 -def ellipsoid(params=None, time_scale=None, d_scale=None, angle_units=None, param_types=None):

508 """Function for setting up a ellipsoidal diffusion tensor. 509 510 @param params: The diffusion tensor parameter. 511 @type params: float 512 @param time_scale: The correlation time scaling value. 513 @type time_scale: float 514 @param d_scale: The diffusion tensor eigenvalue scaling value. 515 @type d_scale: float 516 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 517 @type angle_units: str 518 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 519 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 520 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 521 @type param_types: int 522 """ 523 524 # The diffusion type. 525 cdp.diff_tensor.set_type('ellipsoid') 526 527 # (tm, Da, Dr, alpha, beta, gamma). 528 if param_types == 0: 529 # Unpack the tuple. 530 tm, Da, Dr, alpha, beta, gamma = params 531 532 # Scaling. 533 tm = tm * time_scale 534 Da = Da * d_scale 535 536 # Set the parameters. 537 set(value=[tm, Da, Dr], param=['tm', 'Da', 'Dr']) 538 539 # (Diso, Da, Dr, alpha, beta, gamma). 540 elif param_types == 1: 541 # Unpack the tuple. 542 Diso, Da, Dr, alpha, beta, gamma = params 543 544 # Scaling. 545 Diso = Diso * d_scale 546 Da = Da * d_scale 547 548 # Set the parameters. 549 set(value=[Diso, Da, Dr], param=['Diso', 'Da', 'Dr']) 550 551 # (Dx, Dy, Dz, alpha, beta, gamma). 552 elif param_types == 2: 553 # Unpack the tuple. 554 Dx, Dy, Dz, alpha, beta, gamma = params 555 556 # Scaling. 557 Dx = Dx * d_scale 558 Dy = Dy * d_scale 559 Dz = Dz * d_scale 560 561 # Set the parameters. 562 set(value=[Dx, Dy, Dz], param=['Dx', 'Dy', 'Dz']) 563 564 # (Dxx, Dyy, Dzz, Dxy, Dxz, Dyz). 565 elif param_types == 3: 566 # Unpack the tuple. 567 Dxx, Dyy, Dzz, Dxy, Dxz, Dyz = params 568 569 # Build the tensor. 570 tensor = zeros((3, 3), float64) 571 tensor[0, 0] = Dxx 572 tensor[1, 1] = Dyy 573 tensor[2, 2] = Dzz 574 tensor[0, 1] = tensor[1, 0] = Dxy 575 tensor[0, 2] = tensor[2, 0] = Dxz 576 tensor[1, 2] = tensor[2, 1] = Dyz 577 578 # Scaling. 579 tensor = tensor * d_scale 580 581 # Eigenvalues. 582 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 583 584 # Set the parameters. 585 set(value=[Di[0], Di[1], Di[2]], param=['Dx', 'Dy', 'Dz']) 586 587 # Change the angular units. 588 angle_units = 'rad' 589 590 # Unknown parameter combination. 591 else: 592 raise RelaxUnknownParamCombError('param_types', param_types) 593 594 # Convert the angles to radians. 595 if angle_units == 'deg': 596 print("Converting the angles from degrees to radian units.") 597 alpha = (alpha / 360.0) * 2.0 * pi 598 beta = (beta / 360.0) * 2.0 * pi 599 gamma = (gamma / 360.0) * 2.0 * pi 600 601 # Set the orientational parameters. 602 set(value=[alpha, beta, gamma], param=['alpha', 'beta', 'gamma'])

603 604

605 -def fold_angles(sim_index=None):

606 """Wrap the Euler or spherical angles and remove the glide reflection and translational symmetries. 607 608 Wrap the angles such that:: 609 610 0 <= theta <= pi, 611 0 <= phi <= 2pi, 612 613 and:: 614 615 0 <= alpha <= 2pi, 616 0 <= beta <= pi, 617 0 <= gamma <= 2pi. 618 619 620 For the simulated values, the angles are wrapped as:: 621 622 theta - pi/2 <= theta_sim <= theta + pi/2 623 phi - pi <= phi_sim <= phi + pi 624 625 and:: 626 627 alpha - pi <= alpha_sim <= alpha + pi 628 beta - pi/2 <= beta_sim <= beta + pi/2 629 gamma - pi <= gamma_sim <= gamma + pi 630 631 632 @param sim_index: The simulation index. If set to None then the actual values will be folded 633 rather than the simulation values. 634 @type sim_index: int or None 635 """ 636 637 638 # Wrap the angles. 639 ################## 640 641 # Spheroid. 642 if cdp.diff_tensor.type == 'spheroid': 643 # Get the current angles. 644 theta = cdp.diff_tensor.theta 645 phi = cdp.diff_tensor.phi 646 647 # Simulated values. 648 if sim_index != None: 649 theta_sim = cdp.diff_tensor.theta_sim[sim_index] 650 phi_sim = cdp.diff_tensor.phi_sim[sim_index] 651 652 # Normal value. 653 if sim_index == None: 654 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta, 0.0, 2.0*pi)) 655 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi, 0.0, 2.0*pi)) 656 657 # Simulated theta and phi values. 658 else: 659 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta_sim, theta - pi, theta + pi), category='sim', sim_index=sim_index) 660 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi_sim, phi - pi, phi + pi), category='sim', sim_index=sim_index) 661 662 # Ellipsoid. 663 elif cdp.diff_tensor.type == 'ellipsoid': 664 # Get the current angles. 665 alpha = cdp.diff_tensor.alpha 666 beta = cdp.diff_tensor.beta 667 gamma = cdp.diff_tensor.gamma 668 669 # Simulated values. 670 if sim_index != None: 671 alpha_sim = cdp.diff_tensor.alpha_sim[sim_index] 672 beta_sim = cdp.diff_tensor.beta_sim[sim_index] 673 gamma_sim = cdp.diff_tensor.gamma_sim[sim_index] 674 675 # Normal value. 676 if sim_index == None: 677 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha, 0.0, 2.0*pi)) 678 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta, 0.0, 2.0*pi)) 679 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma, 0.0, 2.0*pi)) 680 681 # Simulated alpha, beta, and gamma values. 682 else: 683 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha_sim, alpha - pi, alpha + pi), category='sim', sim_index=sim_index) 684 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta_sim, beta - pi, beta + pi), category='sim', sim_index=sim_index) 685 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma_sim, gamma - pi, gamma + pi), category='sim', sim_index=sim_index) 686 687 688 # Remove the glide reflection and translational symmetries. 689 ########################################################### 690 691 # Spheroid. 692 if cdp.diff_tensor.type == 'spheroid': 693 # Normal value. 694 if sim_index == None: 695 # Fold phi inside 0 and pi. 696 if cdp.diff_tensor.phi >= pi: 697 theta, phi = fold_spherical_angles(cdp.diff_tensor.theta, cdp.diff_tensor.phi) 698 cdp.diff_tensor.set(param='theta', value=theta) 699 cdp.diff_tensor.set(param='phi', value=phi) 700 701 # Simulated theta and phi values. 702 else: 703 # Fold phi_sim inside phi-pi/2 and phi+pi/2. 704 if cdp.diff_tensor.phi_sim[sim_index] >= cdp.diff_tensor.phi + pi/2.0: 705 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 706 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] - pi, category='sim', sim_index=sim_index) 707 elif cdp.diff_tensor.phi_sim[sim_index] <= cdp.diff_tensor.phi - pi/2.0: 708 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 709 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] + pi, category='sim', sim_index=sim_index) 710 711 # Ellipsoid. 712 elif cdp.diff_tensor.type == 'ellipsoid': 713 # Normal value. 714 if sim_index == None: 715 # Fold alpha inside 0 and pi. 716 if cdp.diff_tensor.alpha >= pi: 717 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha - pi) 718 719 # Fold beta inside 0 and pi. 720 if cdp.diff_tensor.beta >= pi: 721 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 722 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta - pi) 723 724 # Fold gamma inside 0 and pi. 725 if cdp.diff_tensor.gamma >= pi: 726 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 727 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta) 728 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma - pi) 729 730 # Simulated theta and phi values. 731 else: 732 # Fold alpha_sim inside alpha-pi/2 and alpha+pi/2. 733 if cdp.diff_tensor.alpha_sim[sim_index] >= cdp.diff_tensor.alpha + pi/2.0: 734 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] - pi, category='sim', sim_index=sim_index) 735 elif cdp.diff_tensor.alpha_sim[sim_index] <= cdp.diff_tensor.alpha - pi/2.0: 736 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] + pi, category='sim', sim_index=sim_index) 737 738 # Fold beta_sim inside beta-pi/2 and beta+pi/2. 739 if cdp.diff_tensor.beta_sim[sim_index] >= cdp.diff_tensor.beta + pi/2.0: 740 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 741 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] - pi, category='sim', sim_index=sim_index) 742 elif cdp.diff_tensor.beta_sim[sim_index] <= cdp.diff_tensor.beta - pi/2.0: 743 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 744 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] + pi, category='sim', sim_index=sim_index) 745 746 # Fold gamma_sim inside gamma-pi/2 and gamma+pi/2. 747 if cdp.diff_tensor.gamma_sim[sim_index] >= cdp.diff_tensor.gamma + pi/2.0: 748 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 749 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 750 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] - pi, category='sim', sim_index=sim_index) 751 elif cdp.diff_tensor.gamma_sim[sim_index] <= cdp.diff_tensor.gamma - pi/2.0: 752 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 753 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 754 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] + pi, category='sim', sim_index=sim_index)

755 756

757 -def init(params=None, time_scale=1.0, d_scale=1.0, angle_units='deg', param_types=0, spheroid_type=None, fixed=1):

758 """Function for initialising the diffusion tensor. 759 760 @param params: The diffusion tensor parameters. 761 @type params: float 762 @param time_scale: The correlation time scaling value. 763 @type time_scale: float 764 @param d_scale: The diffusion tensor eigenvalue scaling value. 765 @type d_scale: float 766 @param angle_units: The units for the angle parameters. 767 @type angle_units: str (either 'deg' or 'rad') 768 @param param_types: The type of parameters supplied. For spherical diffusion, the flag 769 values correspond to 0: tm, 1: Diso. For spheroidal diffusion, 0: {tm, 770 Da, theta, phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 771 3: {Dpar, Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. For 772 ellipsoidal diffusion, 0: {tm, Da, Dr, alpha, beta, gamma}, 1: {Diso, 773 Da, Dr, alpha, beta, gamma}, 2: {Dx, Dy, Dz, alpha, beta, gamma}. 774 @type param_types: int 775 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 776 restrict the tensor to either being 'oblate' or 'prolate'. 777 @type spheroid_type: str 778 @param fixed: A flag specifying whether the diffusion tensor is fixed or can be 779 optimised. 780 @type fixed: bin 781 """ 782 783 # Test if the current data pipe exists. 784 pipes.test() 785 786 # Test if diffusion tensor data already exists. 787 if diff_data_exists(): 788 raise RelaxTensorError('diffusion') 789 790 # Check the validity of the angle_units argument. 791 valid_types = ['deg', 'rad'] 792 if not angle_units in valid_types: 793 raise RelaxError("The diffusion tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.") 794 795 # Add the diff_tensor object to the data pipe. 796 cdp.diff_tensor = DiffTensorData() 797 798 # Set the fixed flag. 799 cdp.diff_tensor.set_fixed(fixed) 800 801 # Spherical diffusion. 802 if isinstance(params, float): 803 num_params = 1 804 sphere(params, time_scale, param_types) 805 806 # Spheroidal diffusion. 807 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 4: 808 num_params = 4 809 spheroid(params, time_scale, d_scale, angle_units, param_types, spheroid_type) 810 811 # Ellipsoidal diffusion. 812 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 6: 813 num_params = 6 814 ellipsoid(params, time_scale, d_scale, angle_units, param_types) 815 816 # Unknown. 817 else: 818 raise RelaxError("The diffusion tensor parameters " + repr(params) + " are of an unknown type.") 819 820 # Test the validity of the parameters. 821 test_params(num_params)

822 823

824 -def map_bounds(param, spin_id=None):

825 """The function for creating bounds for the mapping function. 826 827 @param param: The name of the parameter to return the bounds for. 828 @type param: str 829 @keyword spin_id: The spin identification string. This arg is unused. 830 @type spin_id: None or str 831 @return: The bounds for the parameter. 832 @rtype: list of len 2 of floats 833 """ 834 835 # tm. 836 if param == 'tm': 837 return [0, 10.0 * 1e-9] 838 839 # {Diso, Dx, Dy, Dz, Dpar, Dper}. 840 if param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 841 return [1e6, 1e7] 842 843 # Da. 844 if param == 'Da': 845 return [-3.0/2.0 * 1e7, 3.0 * 1e7] 846 847 # Dr. 848 if param == 'Dr': 849 return [0, 1] 850 851 # Dratio. 852 if param == 'Dratio': 853 return [1.0/3.0, 3.0] 854 855 # theta. 856 if param == 'theta': 857 return [0, pi] 858 859 # phi. 860 if param == 'phi': 861 return [0, 2*pi] 862 863 # alpha. 864 if param == 'alpha': 865 return [0, 2*pi] 866 867 # beta. 868 if param == 'beta': 869 return [0, pi] 870 871 # gamma. 872 if param == 'gamma': 873 return [0, 2*pi]

874 875

876 -def map_labels(index, params, bounds, swap, inc):

877 """Function for creating labels, tick locations, and tick values for an OpenDX map.""" 878 879 # Initialise. 880 labels = "{" 881 tick_locations = [] 882 tick_values = [] 883 n = len(params) 884 axis_incs = 5 885 loc_inc = inc / axis_incs 886 887 # Increment over the model parameters. 888 for i in range(n): 889 # Parameter conversion factors. 890 factor = return_conversion_factor(params[swap[i]]) 891 892 # Parameter units. 893 units = return_units(params[swap[i]]) 894 895 # Labels. 896 if units: 897 labels = labels + "\"" + params[swap[i]] + " (" + units + ")\"" 898 else: 899 labels = labels + "\"" + params[swap[i]] + "\"" 900 901 # Tick values. 902 vals = bounds[swap[i], 0] / factor 903 val_inc = (bounds[swap[i], 1] - bounds[swap[i], 0]) / (axis_incs * factor) 904 905 if i < n - 1: 906 labels = labels + " " 907 else: 908 labels = labels + "}" 909 910 # Tick locations. 911 string = "{" 912 val = 0.0 913 for j in range(axis_incs + 1): 914 string = string + " " + repr(val) 915 val = val + loc_inc 916 string = string + " }" 917 tick_locations.append(string) 918 919 # Tick values. 920 string = "{" 921 for j in range(axis_incs + 1): 922 string = string + "\"" + "%.2f" % vals + "\" " 923 vals = vals + val_inc 924 string = string + "}" 925 tick_values.append(string) 926 927 return labels, tick_locations, tick_values

928 929

930 -def return_conversion_factor(param):

931 """Function for returning the factor of conversion between different parameter units. 932 933 For example, the internal representation of tm is in seconds, whereas the external 934 representation is in nanoseconds, therefore this function will return 1e-9 for tm. 935 936 937 @param param: The name of the parameter to return the conversion factor for. 938 @type param: str 939 @return: The conversion factor. 940 @rtype: float 941 """ 942 943 # Get the object name. 944 object_name = return_data_name(param) 945 946 # tm (nanoseconds). 947 if object_name == 'tm': 948 return 1e-9 949 950 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 951 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 952 return 1e6 953 954 # Angles. 955 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 956 return (2.0*pi) / 360.0 957 958 # No conversion factor. 959 return 1.0

960 961

962 -def return_data_name(name):

963 """Return the parameter name. 964 965 @param name: The name of the parameter to return the name of. 966 @type name: str 967 @return: The parameter name. 968 @rtype: str 969 """ 970 971 # Enforce that the name must be a string. 972 if not isinstance(name, str): 973 raise RelaxStrError('name', name) 974 975 # Local tm. 976 if search('^tm$', name): 977 return 'tm' 978 979 # Diso. 980 if search('[Dd]iso', name): 981 return 'Diso' 982 983 # Da. 984 if search('[Dd]a', name): 985 return 'Da' 986 987 # Dr. 988 if search('[Dd]r$', name): 989 return 'Dr' 990 991 # Dx. 992 if search('[Dd]x', name): 993 return 'Dx' 994 995 # Dy. 996 if search('[Dd]y', name): 997 return 'Dy' 998 999 # Dz. 1000 if search('[Dd]z', name): 1001 return 'Dz' 1002 1003 # Dpar. 1004 if search('[Dd]par', name): 1005 return 'Dpar' 1006 1007 # Dper. 1008 if search('[Dd]per', name): 1009 return 'Dper' 1010 1011 # Dratio. 1012 if search('[Dd]ratio', name): 1013 return 'Dratio' 1014 1015 # alpha. 1016 if search('^a$', name) or search('alpha', name): 1017 return 'alpha' 1018 1019 # beta. 1020 if search('^b$', name) or search('beta', name): 1021 return 'beta' 1022 1023 # gamma. 1024 if search('^g$', name) or search('gamma', name): 1025 return 'gamma' 1026 1027 # theta. 1028 if search('theta', name): 1029 return 'theta' 1030 1031 # phi. 1032 if search('phi', name): 1033 return 'phi'

1034 1035 # User function documentation. 1036 __return_data_name_doc__ = Desc_container("Diffusion tensor parameter string matching patterns") 1037 table = uf_tables.add_table(label="table: diff data type patterns", caption="Diffusion tensor parameter string matching patterns.") 1038 table.add_headings(["Data type", "Object name", "Patterns"]) 1039 table.add_row(["Global correlation time - tm", "'tm'", "'^tm$'"]) 1040 table.add_row(["Isotropic component of the diffusion tensor - Diso", "'Diso'", "'[Dd]iso'"]) 1041 table.add_row(["Anisotropic component of the diffusion tensor - Da", "'Da'", "'[Dd]a'"]) 1042 table.add_row(["Rhombic component of the diffusion tensor - Dr", "'Dr'", "'[Dd]r$'"]) 1043 table.add_row(["Eigenvalue associated with the x-axis of the diffusion tensor - Dx", "'Dx'", "'[Dd]x'"]) 1044 table.add_row(["Eigenvalue associated with the y-axis of the diffusion tensor - Dy", "'Dy'", "'[Dd]y'"]) 1045 table.add_row(["Eigenvalue associated with the z-axis of the diffusion tensor - Dz", "'Dz'", "'[Dd]z'"]) 1046 table.add_row(["Diffusion coefficient parallel to the major axis of the spheroid diffusion tensor - Dpar", "'Dpar'", "'[Dd]par'"]) 1047 table.add_row(["Diffusion coefficient perpendicular to the major axis of the spheroid diffusion tensor - Dper", "'Dper'", "'[Dd]per'"]) 1048 table.add_row(["Ratio of the parallel and perpendicular components of the spheroid diffusion tensor - Dratio", "'Dratio'", "'[Dd]ratio'"]) 1049 table.add_row(["The first Euler angle of the ellipsoid diffusion tensor - alpha", "'alpha'", "'^a$' or 'alpha'"]) 1050 table.add_row(["The second Euler angle of the ellipsoid diffusion tensor - beta", "'beta'", "'^b$' or 'beta'"]) 1051 table.add_row(["The third Euler angle of the ellipsoid diffusion tensor - gamma", "'gamma'", "'^g$' or 'gamma'"]) 1052 table.add_row(["The polar angle defining the major axis of the spheroid diffusion tensor - theta", "'theta'", "'theta'"]) 1053 table.add_row(["The azimuthal angle defining the major axis of the spheroid diffusion tensor - phi", "'phi'", "'phi'"]) 1054 __return_data_name_doc__.add_table(table.label) 1055 1056

1057 -def return_eigenvalues():

1058 """Function for returning Dx, Dy, and Dz.""" 1059 1060 # Reassign the data. 1061 data = cdp.diff_tensor 1062 1063 # Diso. 1064 Diso = 1.0 / (6.0 * data.tm) 1065 1066 # Dx. 1067 Dx = Diso - 1.0/3.0 * data.Da * (1.0 + 3.0 * data.Dr) 1068 1069 # Dy. 1070 Dy = Diso - 1.0/3.0 * data.Da * (1.0 - 3.0 * data.Dr) 1071 1072 # Dz. 1073 Dz = Diso + 2.0/3.0 * data.Da 1074 1075 # Return the eigenvalues. 1076 return Dx, Dy, Dz

1077 1078

1079 -def return_units(param):

1080 """Function for returning a string representing the parameters units. 1081 1082 For example, the internal representation of tm is in seconds, whereas the external 1083 representation is in nanoseconds, therefore this function will return the string 1084 'nanoseconds' for tm. 1085 1086 1087 @param param: The name of the parameter to return the units for. 1088 @type param: str 1089 @return: The parameter units string. 1090 @rtype: str 1091 """ 1092 1093 # Get the object name. 1094 object_name = return_data_name(param) 1095 1096 # tm (nanoseconds). 1097 if object_name == 'tm': 1098 return 'ns' 1099 1100 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 1101 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 1102 return '1e6 1/s' 1103 1104 # Angles. 1105 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1106 return 'deg'

1107 1108

1109 -def set(value=None, param=None):

1110 """Set the diffusion tensor parameters. 1111 1112 @keyword tensor: The diffusion tensor object. 1113 @type tensor: DiffTensorData instance 1114 @keyword value: The list of values to set the parameters to. 1115 @type value: list of float 1116 @keyword param: The list of parameter names. 1117 @type param: list of str 1118 """ 1119 1120 # Set up the diffusion tensor data if it doesn't exist. 1121 if not diff_data_exists(): 1122 raise RelaxNoTensorError('diffusion') 1123 1124 # Initialise. 1125 geo_params = [] 1126 geo_values = [] 1127 orient_params = [] 1128 orient_values = [] 1129 1130 # Loop over the parameters. 1131 for i in range(len(param)): 1132 # Get the object name. 1133 param[i] = return_data_name(param[i]) 1134 1135 # Unknown parameter. 1136 if not param[i]: 1137 raise RelaxUnknownParamError("diffusion tensor", param[i]) 1138 1139 # Default value. 1140 if value[i] == None: 1141 value[i] = default_value(param[i]) 1142 1143 # Geometric parameter. 1144 if param[i] in ['tm', 'Diso', 'Da', 'Dratio', 'Dper', 'Dpar', 'Dr', 'Dx', 'Dy', 'Dz']: 1145 geo_params.append(param[i]) 1146 geo_values.append(value[i]) 1147 1148 # Orientational parameter. 1149 if param[i] in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1150 orient_params.append(param[i]) 1151 orient_values.append(value[i]) 1152 1153 1154 # Spherical diffusion. 1155 ###################### 1156 1157 if cdp.diff_tensor.type == 'sphere': 1158 # Geometric parameters. 1159 ####################### 1160 1161 # A single geometric parameter. 1162 if len(geo_params) == 1: 1163 # The single parameter tm. 1164 if geo_params[0] == 'tm': 1165 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1166 1167 # The single parameter Diso. 1168 elif geo_params[0] == 'Diso': 1169 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1170 1171 # Cannot set the single parameter. 1172 else: 1173 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1174 1175 # More than one geometric parameters. 1176 elif len(geo_params) > 1: 1177 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1178 1179 1180 # Orientational parameters. 1181 ########################### 1182 1183 # ??? 1184 if len(orient_params): 1185 raise RelaxError("For spherical diffusion, the orientation parameters " + repr(orient_params) + " should not exist.") 1186 1187 1188 # Spheroidal diffusion. 1189 ####################### 1190 1191 elif cdp.diff_tensor.type == 'spheroid': 1192 # Geometric parameters. 1193 ####################### 1194 1195 # A single geometric parameter. 1196 if len(geo_params) == 1: 1197 # The single parameter tm. 1198 if geo_params[0] == 'tm': 1199 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1200 1201 # The single parameter Diso. 1202 elif geo_params[0] == 'Diso': 1203 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1204 1205 # The single parameter Da. 1206 elif geo_params[0] == 'Da': 1207 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 1208 1209 # The single parameter Dratio. 1210 elif geo_params[0] == 'Dratio': 1211 Dratio = geo_values[0] 1212 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * cdp.diff_tensor.tm * (Dratio + 2.0))) 1213 1214 # Cannot set the single parameter. 1215 else: 1216 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1217 1218 # Two geometric parameters. 1219 elif len(geo_params) == 2: 1220 # The geometric parameter set {tm, Da}. 1221 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1222 # The parameters. 1223 tm = geo_values[geo_params.index('tm')] 1224 Da = geo_values[geo_params.index('Da')] 1225 1226 # Set the internal parameter values. 1227 cdp.diff_tensor.set(param='tm', value=tm) 1228 cdp.diff_tensor.set(param='Da', value=Da) 1229 1230 # The geometric parameter set {Diso, Da}. 1231 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1232 # The parameters. 1233 Diso = geo_values[geo_params.index('Diso')] 1234 Da = geo_values[geo_params.index('Da')] 1235 1236 # Set the internal parameter values. 1237 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1238 cdp.diff_tensor.set(param='Da', value=Da) 1239 1240 # The geometric parameter set {tm, Dratio}. 1241 elif geo_params.count('tm') == 1 and geo_params.count('Dratio') == 1: 1242 # The parameters. 1243 tm = geo_values[geo_params.index('tm')] 1244 Dratio = geo_values[geo_params.index('Dratio')] 1245 1246 # Set the internal parameter values. 1247 cdp.diff_tensor.set(param='tm', value=tm) 1248 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * tm * (Dratio + 2.0))) 1249 1250 # The geometric parameter set {Dpar, Dper}. 1251 elif geo_params.count('Dpar') == 1 and geo_params.count('Dper') == 1: 1252 # The parameters. 1253 Dpar = geo_values[geo_params.index('Dpar')] 1254 Dper = geo_values[geo_params.index('Dper')] 1255 1256 # Set the internal parameter values. 1257 cdp.diff_tensor.set(param='tm', value=1.0 / (2.0 * (Dpar + 2.0*Dper))) 1258 cdp.diff_tensor.set(param='Da', value=Dpar - Dper) 1259 1260 # The geometric parameter set {Diso, Dratio}. 1261 elif geo_params.count('Diso') == 1 and geo_params.count('Dratio') == 1: 1262 # The parameters. 1263 Diso = geo_values[geo_params.index('Diso')] 1264 Dratio = geo_values[geo_params.index('Dratio')] 1265 1266 # Set the internal parameter values. 1267 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1268 cdp.diff_tensor.set(param='Da', value=3.0 * Diso * (Dratio - 1.0) / (Dratio + 2.0)) 1269 1270 # Unknown parameter combination. 1271 else: 1272 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1273 1274 # More than two geometric parameters. 1275 elif len(geo_params) > 2: 1276 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1277 1278 1279 # Orientational parameters. 1280 ########################### 1281 1282 # A single orientational parameter. 1283 if len(orient_params) == 1: 1284 # The single parameter theta. 1285 if orient_params[0] == 'theta': 1286 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 1287 1288 # The single parameter phi. 1289 elif orient_params[0] == 'phi': 1290 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 1291 1292 # Disallowed parameter. 1293 else: 1294 raise RelaxError("For spheroidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1295 1296 # Two orientational parameters. 1297 elif len(orient_params) == 2: 1298 # The orientational parameter set {theta, phi}. 1299 if orient_params.count('theta') == 1 and orient_params.count('phi') == 1: 1300 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 1301 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 1302 1303 # Unknown parameter combination. 1304 else: 1305 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1306 1307 # More than two orientational parameters. 1308 elif len(orient_params) > 2: 1309 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1310 1311 1312 # Ellipsoidal diffusion. 1313 ######################## 1314 1315 elif cdp.diff_tensor.type == 'ellipsoid': 1316 # Geometric parameters. 1317 ####################### 1318 1319 # A single geometric parameter. 1320 if len(geo_params) == 1: 1321 # The single parameter tm. 1322 if geo_params[0] == 'tm': 1323 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1324 1325 # The single parameter Diso. 1326 elif geo_params[0] == 'Diso': 1327 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1328 1329 # The single parameter Da. 1330 elif geo_params[0] == 'Da': 1331 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 1332 1333 # The single parameter Dr. 1334 elif geo_params[0] == 'Dr': 1335 cdp.diff_tensor.set(param='Dr', value=geo_values[0]) 1336 1337 # Cannot set the single parameter. 1338 else: 1339 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1340 1341 # Two geometric parameters. 1342 elif len(geo_params) == 2: 1343 # The geometric parameter set {tm, Da}. 1344 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1345 # The parameters. 1346 tm = geo_values[geo_params.index('tm')] 1347 Da = geo_values[geo_params.index('Da')] 1348 1349 # Set the internal parameter values. 1350 cdp.diff_tensor.set(param='tm', value=tm) 1351 cdp.diff_tensor.set(param='Da', value=Da) 1352 1353 # The geometric parameter set {tm, Dr}. 1354 elif geo_params.count('tm') == 1 and geo_params.count('Dr') == 1: 1355 # The parameters. 1356 tm = geo_values[geo_params.index('tm')] 1357 Dr = geo_values[geo_params.index('Dr')] 1358 1359 # Set the internal parameter values. 1360 cdp.diff_tensor.set(param='tm', value=tm) 1361 cdp.diff_tensor.set(param='Dr', value=Dr) 1362 1363 # The geometric parameter set {Diso, Da}. 1364 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1365 # The parameters. 1366 Diso = geo_values[geo_params.index('Diso')] 1367 Da = geo_values[geo_params.index('Da')] 1368 1369 # Set the internal parameter values. 1370 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1371 cdp.diff_tensor.set(param='Da', value=Da) 1372 1373 # The geometric parameter set {Diso, Dr}. 1374 elif geo_params.count('Diso') == 1 and geo_params.count('Dr') == 1: 1375 # The parameters. 1376 Diso = geo_values[geo_params.index('Diso')] 1377 Dr = geo_values[geo_params.index('Dr')] 1378 1379 # Set the internal parameter values. 1380 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1381 cdp.diff_tensor.set(param='Dr', value=Dr) 1382 1383 # The geometric parameter set {Da, Dr}. 1384 elif geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1385 # The parameters. 1386 Da = geo_values[geo_params.index('Da')] 1387 Dr = geo_values[geo_params.index('Dr')] 1388 1389 # Set the internal parameter values. 1390 cdp.diff_tensor.set(param='Da', value=Da) 1391 cdp.diff_tensor.set(param='Da', value=Dr) 1392 1393 # Unknown parameter combination. 1394 else: 1395 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1396 1397 # Three geometric parameters. 1398 elif len(geo_params) == 3: 1399 # The geometric parameter set {tm, Da, Dr}. 1400 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1401 # The parameters. 1402 tm = geo_values[geo_params.index('tm')] 1403 Da = geo_values[geo_params.index('Da')] 1404 Dr = geo_values[geo_params.index('Dr')] 1405 1406 # Set the internal parameter values. 1407 cdp.diff_tensor.set(param='tm', value=tm) 1408 cdp.diff_tensor.set(param='Da', value=Da) 1409 cdp.diff_tensor.set(param='Dr', value=Dr) 1410 1411 # The geometric parameter set {Diso, Da, Dr}. 1412 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1413 # The parameters. 1414 Diso = geo_values[geo_params.index('Diso')] 1415 Da = geo_values[geo_params.index('Da')] 1416 Dr = geo_values[geo_params.index('Dr')] 1417 1418 # Set the internal parameter values. 1419 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1420 cdp.diff_tensor.set(param='Da', value=Da) 1421 cdp.diff_tensor.set(param='Dr', value=Dr) 1422 1423 # The geometric parameter set {Dx, Dy, Dz}. 1424 elif geo_params.count('Dx') == 1 and geo_params.count('Dy') == 1 and geo_params.count('Dz') == 1: 1425 # The parameters. 1426 Dx = geo_values[geo_params.index('Dx')] 1427 Dy = geo_values[geo_params.index('Dy')] 1428 Dz = geo_values[geo_params.index('Dz')] 1429 1430 # Set the internal tm value. 1431 if Dx + Dy + Dz == 0.0: 1432 cdp.diff_tensor.set(param='tm', value=1e99) 1433 else: 1434 cdp.diff_tensor.set(param='tm', value=0.5 / (Dx + Dy + Dz)) 1435 1436 # Set the internal Da value. 1437 cdp.diff_tensor.set(param='Da', value=Dz - 0.5*(Dx + Dy)) 1438 1439 # Set the internal Dr value. 1440 if cdp.diff_tensor.Da == 0.0: 1441 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) * 1e99) 1442 else: 1443 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) / (2.0*cdp.diff_tensor.Da)) 1444 1445 # Unknown parameter combination. 1446 else: 1447 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1448 1449 1450 # More than three geometric parameters. 1451 elif len(geo_params) > 3: 1452 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1453 1454 1455 # Orientational parameters. 1456 ########################### 1457 1458 # A single orientational parameter. 1459 if len(orient_params) == 1: 1460 # The single parameter alpha. 1461 if orient_params[0] == 'alpha': 1462 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1463 1464 # The single parameter beta. 1465 elif orient_params[0] == 'beta': 1466 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1467 1468 # The single parameter gamma. 1469 elif orient_params[0] == 'gamma': 1470 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1471 1472 # Disallowed parameter. 1473 else: 1474 raise RelaxError("For ellipsoidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1475 1476 # Two orientational parameters. 1477 elif len(orient_params) == 2: 1478 # The orientational parameter set {alpha, beta}. 1479 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1480 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1481 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1482 1483 # The orientational parameter set {alpha, gamma}. 1484 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1: 1485 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1486 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1487 1488 # The orientational parameter set {beta, gamma}. 1489 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1: 1490 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1491 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1492 1493 # Unknown parameter combination. 1494 else: 1495 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1496 1497 # Three orientational parameters. 1498 elif len(orient_params) == 3: 1499 # The orientational parameter set {alpha, beta, gamma}. 1500 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1501 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1502 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1503 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1504 1505 # Unknown parameter combination. 1506 else: 1507 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1508 1509 # More than three orientational parameters. 1510 elif len(orient_params) > 3: 1511 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1512 1513 1514 # Fold the angles in. 1515 ##################### 1516 1517 if orient_params: 1518 fold_angles()

1519 1520 # User function documentation. 1521 __set_doc__ = Desc_container("Diffusion tensor set details") 1522 __set_doc__.add_paragraph("If the diffusion tensor has not been setup, use the more powerful function 'diffusion_tensor.init' to initialise the tensor parameters. This function cannot be used to initialise a diffusion tensor.") 1523 __set_doc__.add_paragraph("The units of the parameters are:") 1524 __set_doc__.add_list_element("Inverse seconds for tm.") 1525 __set_doc__.add_list_element("Seconds for Diso, Da, Dx, Dy, Dz, Dpar, Dper.") 1526 __set_doc__.add_list_element("Unitless for Dratio and Dr.") 1527 __set_doc__.add_list_element("Radians for all angles (alpha, beta, gamma, theta, phi).") 1528 __set_doc__.add_paragraph("When setting a diffusion tensor parameter, the residue number has no effect. As the internal parameters of spherical diffusion are {tm}, spheroidal diffusion are {tm, Da, theta, phi}, and ellipsoidal diffusion are {tm, Da, Dr, alpha, beta, gamma}, supplying geometric parameters must be done in the following way. If a single geometric parameter is supplied, it must be one of tm, Diso, Da, Dr, or Dratio. For the parameters Dpar, Dper, Dx, Dy, and Dx, it is not possible to determine how to use the currently set values together with the supplied value to calculate the new internal parameters. For spheroidal diffusion, when supplying multiple geometric parameters, the set must belong to one of") 1529 __set_doc__.add_list_element("{tm, Da},") 1530 __set_doc__.add_list_element("{Diso, Da},") 1531 __set_doc__.add_list_element("{tm, Dratio},") 1532 __set_doc__.add_list_element("{Dpar, Dper},") 1533 __set_doc__.add_list_element("{Diso, Dratio},") 1534 __set_doc__.add_paragraph("where either theta, phi, or both orientational parameters can be additionally supplied. For ellipsoidal diffusion, again when supplying multiple geometric parameters, the set must belong to one of") 1535 __set_doc__.add_list_element("{tm, Da, Dr},") 1536 __set_doc__.add_list_element("{Diso, Da, Dr},") 1537 __set_doc__.add_list_element("{Dx, Dy, Dz},") 1538 __set_doc__.add_paragraph("where any number of the orientational parameters, alpha, beta, or gamma can be additionally supplied.") 1539 1540

1541 -def sphere(params=None, time_scale=None, param_types=None):

1542 """Function for setting up a spherical diffusion tensor. 1543 1544 @param params: The diffusion tensor parameter. 1545 @type params: float 1546 @param time_scale: The correlation time scaling value. 1547 @type time_scale: float 1548 @param param_types: The type of parameter supplied. If 0, then the parameter is tm. If 1, then 1549 the parameter is Diso. 1550 @type param_types: int 1551 """ 1552 1553 # The diffusion type. 1554 cdp.diff_tensor.set_type('sphere') 1555 1556 # tm. 1557 if param_types == 0: 1558 # Scaling. 1559 tm = params * time_scale 1560 1561 # Set the parameters. 1562 set(value=[tm], param=['tm']) 1563 1564 # Diso 1565 elif param_types == 1: 1566 # Scaling. 1567 Diso = params * d_scale 1568 1569 # Set the parameters. 1570 set(value=[Diso], param=['Diso']) 1571 1572 # Unknown parameter combination. 1573 else: 1574 raise RelaxUnknownParamCombError('param_types', param_types)

1575 1576

1577 -def spheroid(params=None, time_scale=None, d_scale=None, angle_units=None, param_types=None, spheroid_type=None):

1578 """Function for setting up a spheroidal diffusion tensor. 1579 1580 @param params: The diffusion tensor parameter. 1581 @type params: float 1582 @param time_scale: The correlation time scaling value. 1583 @type time_scale: float 1584 @param d_scale: The diffusion tensor eigenvalue scaling value. 1585 @type d_scale: float 1586 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 1587 @type angle_units: str 1588 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 1589 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 1590 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 1591 @type param_types: int 1592 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 1593 restrict the tensor to either being 'oblate' or 'prolate'. 1594 @type spheroid_type: str 1595 """ 1596 1597 # The diffusion type. 1598 cdp.diff_tensor.set_type('spheroid') 1599 1600 # Spheroid diffusion type. 1601 allowed_types = [None, 'oblate', 'prolate'] 1602 if spheroid_type not in allowed_types: 1603 raise RelaxError("The 'spheroid_type' argument " + repr(spheroid_type) + " should be 'oblate', 'prolate', or None.") 1604 cdp.diff_tensor.set(param='spheroid_type', value=spheroid_type) 1605 1606 # (tm, Da, theta, phi). 1607 if param_types == 0: 1608 # Unpack the tuple. 1609 tm, Da, theta, phi = params 1610 1611 # Scaling. 1612 tm = tm * time_scale 1613 Da = Da * d_scale 1614 1615 # Set the parameters. 1616 set(value=[tm, Da], param=['tm', 'Da']) 1617 1618 # (Diso, Da, theta, phi). 1619 elif param_types == 1: 1620 # Unpack the tuple. 1621 Diso, Da, theta, phi = params 1622 1623 # Scaling. 1624 Diso = Diso * d_scale 1625 Da = Da * d_scale 1626 1627 # Set the parameters. 1628 set(value=[Diso, Da], param=['Diso', 'Da']) 1629 1630 # (tm, Dratio, theta, phi). 1631 elif param_types == 2: 1632 # Unpack the tuple. 1633 tm, Dratio, theta, phi = params 1634 1635 # Scaling. 1636 tm = tm * time_scale 1637 1638 # Set the parameters. 1639 set(value=[tm, Dratio], param=['tm', 'Dratio']) 1640 1641 # (Dpar, Dper, theta, phi). 1642 elif param_types == 3: 1643 # Unpack the tuple. 1644 Dpar, Dper, theta, phi = params 1645 1646 # Scaling. 1647 Dpar = Dpar * d_scale 1648 Dper = Dper * d_scale 1649 1650 # Set the parameters. 1651 set(value=[Dpar, Dper], param=['Dpar', 'Dper']) 1652 1653 # (Diso, Dratio, theta, phi). 1654 elif param_types == 4: 1655 # Unpack the tuple. 1656 Diso, Dratio, theta, phi = params 1657 1658 # Scaling. 1659 Diso = Diso * d_scale 1660 1661 # Set the parameters. 1662 set(value=[Diso, Dratio], param=['Diso', 'Dratio']) 1663 1664 # Unknown parameter combination. 1665 else: 1666 raise RelaxUnknownParamCombError('param_types', param_types) 1667 1668 # Convert the angles to radians. 1669 if angle_units == 'deg': 1670 print("Converting the angles from degrees to radian units.") 1671 theta = (theta / 360.0) * 2.0 * pi 1672 phi = (phi / 360.0) * 2.0 * pi 1673 1674 # Set the orientational parameters. 1675 set(value=[theta, phi], param=['theta', 'phi'])

1676 1677

1678 -def tensor_eigen_system(tensor):

1679 """Determine the eigenvalues and vectors for the tensor, sorting the entries. 1680 1681 @return: The eigenvalues, rotation matrix, and the Euler angles in zyz notation. 1682 @rtype: 3D rank-1 array, 3D rank-2 array, float, float, float 1683 """ 1684 1685 # Eigenvalues. 1686 R, Di, A = svd(tensor) 1687 D_diag = zeros((3, 3), float64) 1688 for i in range(3): 1689 D_diag[i, i] = Di[i] 1690 1691 # Reordering structure. 1692 reorder_data = [] 1693 for i in range(3): 1694 reorder_data.append([Di[i], i]) 1695 reorder_data.sort() 1696 1697 # The indices. 1698 reorder = zeros(3, int) 1699 Di_sort = zeros(3, float) 1700 for i in range(3): 1701 Di_sort[i], reorder[i] = reorder_data[i] 1702 1703 # Reorder columns. 1704 R_new = zeros((3, 3), float64) 1705 for i in range(3): 1706 R_new[:, i] = R[:, reorder[i]] 1707 1708 # Switch from the left handed to right handed universes (if needed). 1709 if norm(cross(R_new[:, 0], R_new[:, 1]) - R_new[:, 2]) > 1e-7: 1710 R_new[:, 2] = -R_new[:, 2] 1711 1712 # Reverse the rotation. 1713 R_new = transpose(R_new) 1714 1715 # Euler angles (reverse rotation in the rotated axis system). 1716 gamma, beta, alpha = R_to_euler_zyz(R_new) 1717 1718 # Collapse the pi axis rotation symmetries. 1719 if alpha >= pi: 1720 alpha = alpha - pi 1721 if gamma >= pi: 1722 alpha = pi - alpha 1723 beta = pi - beta 1724 gamma = gamma - pi 1725 if beta >= pi: 1726 alpha = pi - alpha 1727 beta = beta - pi 1728 1729 # Return the values. 1730 return Di_sort, R_new, alpha, beta, gamma

1731 1732

1733 -def test_params(num_params):

1734 """Function for testing the validity of the input parameters.""" 1735 1736 # An allowable error to account for machine precision, optimisation quality, etc. 1737 error = 1e-4 1738 1739 # tm. 1740 tm = cdp.diff_tensor.tm 1741 if tm <= 0.0 or tm > 1e-6: 1742 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1743 1744 # Spheroid. 1745 if num_params == 4: 1746 # Parameters. 1747 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1748 Da = cdp.diff_tensor.Da 1749 1750 # Da. 1751 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1752 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1753 1754 # Ellipsoid. 1755 if num_params == 6: 1756 # Parameters. 1757 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1758 Da = cdp.diff_tensor.Da 1759 Dr = cdp.diff_tensor.Dr 1760 1761 # Da. 1762 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1763 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1764 1765 # Dr. 1766 if Dr < (0.0 - error) or Dr > (1.0 + error): 1767 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1768 1769

1770 -def unit_axes():

1771 """Function for calculating the unit axes of the diffusion tensor. 1772 1773 Spheroid 1774 ======== 1775 1776 The unit Dpar vector is:: 1777 1778 | sin(theta) * cos(phi) | 1779 Dpar = | sin(theta) * sin(phi) | 1780 | cos(theta) | 1781 1782 1783 Ellipsoid 1784 ========= 1785 1786 The unit Dx vector is:: 1787 1788 | -sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma) | 1789 Dx = | -sin(alpha) * cos(gamma) - cos(alpha) * cos(beta) * sin(gamma) | 1790 | cos(alpha) * sin(beta) | 1791 1792 The unit Dy vector is:: 1793 1794 | cos(alpha) * sin(gamma) + sin(alpha) * cos(beta) * cos(gamma) | 1795 Dy = | cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma) | 1796 | sin(alpha) * sin(beta) | 1797 1798 The unit Dz vector is:: 1799 1800 | -sin(beta) * cos(gamma) | 1801 Dz = | sin(beta) * sin(gamma) | 1802 | cos(beta) | 1803 1804 """ 1805 1806 # Spheroid. 1807 if cdp.diff_tensor.type == 'spheroid': 1808 # Initialise. 1809 Dpar = zeros(3, float64) 1810 1811 # Trig. 1812 sin_theta = sin(cdp.diff_tensor.theta) 1813 cos_theta = cos(cdp.diff_tensor.theta) 1814 sin_phi = sin(cdp.diff_tensor.phi) 1815 cos_phi = cos(cdp.diff_tensor.phi) 1816 1817 # Unit Dpar axis. 1818 Dpar[0] = sin_theta * cos_phi 1819 Dpar[1] = sin_theta * sin_phi 1820 Dpar[2] = cos_theta 1821 1822 # Return the vector. 1823 return Dpar 1824 1825 # Ellipsoid. 1826 if cdp.diff_tensor.type == 'ellipsoid': 1827 # Initialise. 1828 Dx = zeros(3, float64) 1829 Dy = zeros(3, float64) 1830 Dz = zeros(3, float64) 1831 1832 # Trig. 1833 sin_alpha = sin(cdp.diff_tensor.alpha) 1834 cos_alpha = cos(cdp.diff_tensor.alpha) 1835 sin_beta = sin(cdp.diff_tensor.beta) 1836 cos_beta = cos(cdp.diff_tensor.beta) 1837 sin_gamma = sin(cdp.diff_tensor.gamma) 1838 cos_gamma = cos(cdp.diff_tensor.gamma) 1839 1840 # Unit Dx axis. 1841 Dx[0] = -sin_alpha * sin_gamma + cos_alpha * cos_beta * cos_gamma 1842 Dx[1] = -sin_alpha * cos_gamma - cos_alpha * cos_beta * sin_gamma 1843 Dx[2] = cos_alpha * sin_beta 1844 1845 # Unit Dy axis. 1846 Dx[0] = cos_alpha * sin_gamma + sin_alpha * cos_beta * cos_gamma 1847 Dx[1] = cos_alpha * cos_gamma - sin_alpha * cos_beta * sin_gamma 1848 Dx[2] = sin_alpha * sin_beta 1849 1850 # Unit Dz axis. 1851 Dx[0] = -sin_beta * cos_gamma 1852 Dx[1] = sin_beta * sin_gamma 1853 Dx[2] = cos_beta 1854 1855 # Return the vectors. 1856 return Dx, Dy, Dz

1857

Source Code for Module generic_fns.diffusion_tensor