pipe_control.diffusion

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2013 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program is free software: you can redistribute it and/or modify # 8 # it under the terms of the GNU General Public License as published by # 9 # the Free Software Foundation, either version 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program is distributed in the hope that it will be useful, # 13 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 15 # GNU General Public License for more details. # 16 # # 17 # You should have received a copy of the GNU General Public License # 18 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 # Module docstring. 23 """Module for the 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 float64, int32, zeros 29 from re import search 30 import string 31 32 # relax module imports. 33 from data_store.diff_tensor import DiffTensorData 34 from lib.diffusion.main import tensor_eigen_system, tensor_info_table 35 from lib.errors import RelaxError, RelaxNoTensorError, RelaxStrError, RelaxTensorError, RelaxUnknownParamCombError, RelaxUnknownParamError 36 from lib.physical_constants import element_from_isotope, number_from_isotope 37 from pipe_control import pipes 38 from pipe_control.angles import fold_spherical_angles, wrap_angles 39 from pipe_control.interatomic import return_interatom_list 40 from pipe_control.mol_res_spin import get_molecule_names, return_spin, spin_loop 41 from user_functions.data import Uf_tables; uf_tables = Uf_tables() 42 from user_functions.objects import Desc_container 43 44

45 -def bmrb_read(star):

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

142 143

144 -def bmrb_write(star):

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

268 269 270

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

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

308 309

310 -def data_names():

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

317 318

319 -def default_value(param):

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

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

370 -def delete():

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

382 383

384 -def diff_data_exists(pipe=None):

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

405 406

407 -def display():

408 """Function for displaying the diffusion tensor.""" 409 410 # Test if the current data pipe exists. 411 pipes.test() 412 413 # Test if diffusion tensor data exists. 414 if not diff_data_exists(): 415 raise RelaxNoTensorError('diffusion') 416 417 # Printout. 418 if cdp.diff_tensor.type == 'sphere': 419 tensor_info_table(type='sphere', tm=cdp.diff_tensor.tm, Diso=cdp.diff_tensor.Diso, fixed=cdp.diff_tensor.fixed) 420 elif cdp.diff_tensor.type == 'spheroid': 421 tensor_info_table(type='spheroid', tm=cdp.diff_tensor.tm, Diso=cdp.diff_tensor.Diso, Da=cdp.diff_tensor.Da, Dpar=cdp.diff_tensor.Dpar, Dper=cdp.diff_tensor.Dpar, Dratio=cdp.diff_tensor.Dratio, theta=cdp.diff_tensor.theta, phi=cdp.diff_tensor.phi, fixed=cdp.diff_tensor.fixed) 422 elif cdp.diff_tensor.type == 'ellipsoid': 423 tensor_info_table(type='ellipsoid', tm=cdp.diff_tensor.tm, Diso=cdp.diff_tensor.Diso, Da=cdp.diff_tensor.Da, Dr=cdp.diff_tensor.Dr, Dx=cdp.diff_tensor.Dx, Dy=cdp.diff_tensor.Dy, Dz=cdp.diff_tensor.Dz, alpha=cdp.diff_tensor.alpha, beta=cdp.diff_tensor.beta, gamma=cdp.diff_tensor.gamma, fixed=cdp.diff_tensor.fixed)

424 425

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

427 """Function for setting up a ellipsoidal diffusion tensor. 428 429 @param params: The diffusion tensor parameter. 430 @type params: float 431 @param time_scale: The correlation time scaling value. 432 @type time_scale: float 433 @param d_scale: The diffusion tensor eigenvalue scaling value. 434 @type d_scale: float 435 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 436 @type angle_units: str 437 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 438 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 439 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 440 @type param_types: int 441 """ 442 443 # The diffusion type. 444 cdp.diff_tensor.set_type('ellipsoid') 445 446 # (tm, Da, Dr, alpha, beta, gamma). 447 if param_types == 0: 448 # Unpack the tuple. 449 tm, Da, Dr, alpha, beta, gamma = params 450 451 # Scaling. 452 tm = tm * time_scale 453 Da = Da * d_scale 454 455 # Set the parameters. 456 set(value=[tm, Da, Dr], param=['tm', 'Da', 'Dr']) 457 458 # (Diso, Da, Dr, alpha, beta, gamma). 459 elif param_types == 1: 460 # Unpack the tuple. 461 Diso, Da, Dr, alpha, beta, gamma = params 462 463 # Scaling. 464 Diso = Diso * d_scale 465 Da = Da * d_scale 466 467 # Set the parameters. 468 set(value=[Diso, Da, Dr], param=['Diso', 'Da', 'Dr']) 469 470 # (Dx, Dy, Dz, alpha, beta, gamma). 471 elif param_types == 2: 472 # Unpack the tuple. 473 Dx, Dy, Dz, alpha, beta, gamma = params 474 475 # Scaling. 476 Dx = Dx * d_scale 477 Dy = Dy * d_scale 478 Dz = Dz * d_scale 479 480 # Set the parameters. 481 set(value=[Dx, Dy, Dz], param=['Dx', 'Dy', 'Dz']) 482 483 # (Dxx, Dyy, Dzz, Dxy, Dxz, Dyz). 484 elif param_types == 3: 485 # Unpack the tuple. 486 Dxx, Dyy, Dzz, Dxy, Dxz, Dyz = params 487 488 # Build the tensor. 489 tensor = zeros((3, 3), float64) 490 tensor[0, 0] = Dxx 491 tensor[1, 1] = Dyy 492 tensor[2, 2] = Dzz 493 tensor[0, 1] = tensor[1, 0] = Dxy 494 tensor[0, 2] = tensor[2, 0] = Dxz 495 tensor[1, 2] = tensor[2, 1] = Dyz 496 497 # Scaling. 498 tensor = tensor * d_scale 499 500 # Eigenvalues. 501 Di, R, alpha, beta, gamma = tensor_eigen_system(tensor) 502 503 # Set the parameters. 504 set(value=[Di[0], Di[1], Di[2]], param=['Dx', 'Dy', 'Dz']) 505 506 # Change the angular units. 507 angle_units = 'rad' 508 509 # Unknown parameter combination. 510 else: 511 raise RelaxUnknownParamCombError('param_types', param_types) 512 513 # Convert the angles to radians. 514 if angle_units == 'deg': 515 print("Converting the angles from degrees to radian units.") 516 alpha = (alpha / 360.0) * 2.0 * pi 517 beta = (beta / 360.0) * 2.0 * pi 518 gamma = (gamma / 360.0) * 2.0 * pi 519 520 # Set the orientational parameters. 521 set(value=[alpha, beta, gamma], param=['alpha', 'beta', 'gamma'])

522 523

524 -def fold_angles(sim_index=None):

525 """Wrap the Euler or spherical angles and remove the glide reflection and translational symmetries. 526 527 Wrap the angles such that:: 528 529 0 <= theta <= pi, 530 0 <= phi <= 2pi, 531 532 and:: 533 534 0 <= alpha <= 2pi, 535 0 <= beta <= pi, 536 0 <= gamma <= 2pi. 537 538 539 For the simulated values, the angles are wrapped as:: 540 541 theta - pi/2 <= theta_sim <= theta + pi/2 542 phi - pi <= phi_sim <= phi + pi 543 544 and:: 545 546 alpha - pi <= alpha_sim <= alpha + pi 547 beta - pi/2 <= beta_sim <= beta + pi/2 548 gamma - pi <= gamma_sim <= gamma + pi 549 550 551 @param sim_index: The simulation index. If set to None then the actual values will be folded 552 rather than the simulation values. 553 @type sim_index: int or None 554 """ 555 556 557 # Wrap the angles. 558 ################## 559 560 # Spheroid. 561 if cdp.diff_tensor.type == 'spheroid': 562 # Get the current angles. 563 theta = cdp.diff_tensor.theta 564 phi = cdp.diff_tensor.phi 565 566 # Simulated values. 567 if sim_index != None: 568 theta_sim = cdp.diff_tensor.theta_sim[sim_index] 569 phi_sim = cdp.diff_tensor.phi_sim[sim_index] 570 571 # Normal value. 572 if sim_index == None: 573 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta, 0.0, 2.0*pi)) 574 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi, 0.0, 2.0*pi)) 575 576 # Simulated theta and phi values. 577 else: 578 cdp.diff_tensor.set(param='theta', value=wrap_angles(theta_sim, theta - pi, theta + pi), category='sim', sim_index=sim_index) 579 cdp.diff_tensor.set(param='phi', value=wrap_angles(phi_sim, phi - pi, phi + pi), category='sim', sim_index=sim_index) 580 581 # Ellipsoid. 582 elif cdp.diff_tensor.type == 'ellipsoid': 583 # Get the current angles. 584 alpha = cdp.diff_tensor.alpha 585 beta = cdp.diff_tensor.beta 586 gamma = cdp.diff_tensor.gamma 587 588 # Simulated values. 589 if sim_index != None: 590 alpha_sim = cdp.diff_tensor.alpha_sim[sim_index] 591 beta_sim = cdp.diff_tensor.beta_sim[sim_index] 592 gamma_sim = cdp.diff_tensor.gamma_sim[sim_index] 593 594 # Normal value. 595 if sim_index == None: 596 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha, 0.0, 2.0*pi)) 597 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta, 0.0, 2.0*pi)) 598 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma, 0.0, 2.0*pi)) 599 600 # Simulated alpha, beta, and gamma values. 601 else: 602 cdp.diff_tensor.set(param='alpha', value=wrap_angles(alpha_sim, alpha - pi, alpha + pi), category='sim', sim_index=sim_index) 603 cdp.diff_tensor.set(param='beta', value= wrap_angles(beta_sim, beta - pi, beta + pi), category='sim', sim_index=sim_index) 604 cdp.diff_tensor.set(param='gamma', value=wrap_angles(gamma_sim, gamma - pi, gamma + pi), category='sim', sim_index=sim_index) 605 606 607 # Remove the glide reflection and translational symmetries. 608 ########################################################### 609 610 # Spheroid. 611 if cdp.diff_tensor.type == 'spheroid': 612 # Normal value. 613 if sim_index == None: 614 # Fold phi inside 0 and pi. 615 if cdp.diff_tensor.phi >= pi: 616 theta, phi = fold_spherical_angles(cdp.diff_tensor.theta, cdp.diff_tensor.phi) 617 cdp.diff_tensor.set(param='theta', value=theta) 618 cdp.diff_tensor.set(param='phi', value=phi) 619 620 # Simulated theta and phi values. 621 else: 622 # Fold phi_sim inside phi-pi/2 and phi+pi/2. 623 if cdp.diff_tensor.phi_sim[sim_index] >= cdp.diff_tensor.phi + pi/2.0: 624 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 625 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] - pi, category='sim', sim_index=sim_index) 626 elif cdp.diff_tensor.phi_sim[sim_index] <= cdp.diff_tensor.phi - pi/2.0: 627 cdp.diff_tensor.set(param='theta', value=pi - cdp.diff_tensor.theta_sim[sim_index], category='sim', sim_index=sim_index) 628 cdp.diff_tensor.set(param='phi', value=cdp.diff_tensor.phi_sim[sim_index] + pi, category='sim', sim_index=sim_index) 629 630 # Ellipsoid. 631 elif cdp.diff_tensor.type == 'ellipsoid': 632 # Normal value. 633 if sim_index == None: 634 # Fold alpha inside 0 and pi. 635 if cdp.diff_tensor.alpha >= pi: 636 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha - pi) 637 638 # Fold beta inside 0 and pi. 639 if cdp.diff_tensor.beta >= pi: 640 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 641 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta - pi) 642 643 # Fold gamma inside 0 and pi. 644 if cdp.diff_tensor.gamma >= pi: 645 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha) 646 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta) 647 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma - pi) 648 649 # Simulated theta and phi values. 650 else: 651 # Fold alpha_sim inside alpha-pi/2 and alpha+pi/2. 652 if cdp.diff_tensor.alpha_sim[sim_index] >= cdp.diff_tensor.alpha + pi/2.0: 653 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] - pi, category='sim', sim_index=sim_index) 654 elif cdp.diff_tensor.alpha_sim[sim_index] <= cdp.diff_tensor.alpha - pi/2.0: 655 cdp.diff_tensor.set(param='alpha', value=cdp.diff_tensor.alpha_sim[sim_index] + pi, category='sim', sim_index=sim_index) 656 657 # Fold beta_sim inside beta-pi/2 and beta+pi/2. 658 if cdp.diff_tensor.beta_sim[sim_index] >= cdp.diff_tensor.beta + pi/2.0: 659 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 660 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] - pi, category='sim', sim_index=sim_index) 661 elif cdp.diff_tensor.beta_sim[sim_index] <= cdp.diff_tensor.beta - pi/2.0: 662 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 663 cdp.diff_tensor.set(param='beta', value=cdp.diff_tensor.beta_sim[sim_index] + pi, category='sim', sim_index=sim_index) 664 665 # Fold gamma_sim inside gamma-pi/2 and gamma+pi/2. 666 if cdp.diff_tensor.gamma_sim[sim_index] >= cdp.diff_tensor.gamma + pi/2.0: 667 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 668 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 669 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] - pi, category='sim', sim_index=sim_index) 670 elif cdp.diff_tensor.gamma_sim[sim_index] <= cdp.diff_tensor.gamma - pi/2.0: 671 cdp.diff_tensor.set(param='alpha', value=pi - cdp.diff_tensor.alpha_sim[sim_index], category='sim', sim_index=sim_index) 672 cdp.diff_tensor.set(param='beta', value=pi - cdp.diff_tensor.beta_sim[sim_index], category='sim', sim_index=sim_index) 673 cdp.diff_tensor.set(param='gamma', value=cdp.diff_tensor.gamma_sim[sim_index] + pi, category='sim', sim_index=sim_index)

674 675

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

677 """Function for initialising the diffusion tensor. 678 679 @param params: The diffusion tensor parameters. 680 @type params: float 681 @param time_scale: The correlation time scaling value. 682 @type time_scale: float 683 @param d_scale: The diffusion tensor eigenvalue scaling value. 684 @type d_scale: float 685 @param angle_units: The units for the angle parameters. 686 @type angle_units: str (either 'deg' or 'rad') 687 @param param_types: The type of parameters supplied. For spherical diffusion, the flag 688 values correspond to 0: tm, 1: Diso. For spheroidal diffusion, 0: {tm, 689 Da, theta, phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 690 3: {Dpar, Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. For 691 ellipsoidal diffusion, 0: {tm, Da, Dr, alpha, beta, gamma}, 1: {Diso, 692 Da, Dr, alpha, beta, gamma}, 2: {Dx, Dy, Dz, alpha, beta, gamma}. 693 @type param_types: int 694 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 695 restrict the tensor to either being 'oblate' or 'prolate'. 696 @type spheroid_type: str 697 @param fixed: A flag specifying whether the diffusion tensor is fixed or can be 698 optimised. 699 @type fixed: bin 700 """ 701 702 # Test if the current data pipe exists. 703 pipes.test() 704 705 # Test if diffusion tensor data already exists. 706 if diff_data_exists(): 707 raise RelaxTensorError('diffusion') 708 709 # Check the validity of the angle_units argument. 710 valid_types = ['deg', 'rad'] 711 if not angle_units in valid_types: 712 raise RelaxError("The diffusion tensor 'angle_units' argument " + repr(angle_units) + " should be either 'deg' or 'rad'.") 713 714 # Add the diff_tensor object to the data pipe. 715 cdp.diff_tensor = DiffTensorData() 716 717 # Set the fixed flag. 718 cdp.diff_tensor.set_fixed(fixed) 719 720 # Spherical diffusion. 721 if isinstance(params, float): 722 num_params = 1 723 sphere(params, time_scale, param_types) 724 725 # Spheroidal diffusion. 726 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 4: 727 num_params = 4 728 spheroid(params, time_scale, d_scale, angle_units, param_types, spheroid_type) 729 730 # Ellipsoidal diffusion. 731 elif (isinstance(params, tuple) or isinstance(params, list)) and len(params) == 6: 732 num_params = 6 733 ellipsoid(params, time_scale, d_scale, angle_units, param_types) 734 735 # Unknown. 736 else: 737 raise RelaxError("The diffusion tensor parameters " + repr(params) + " are of an unknown type.") 738 739 # Test the validity of the parameters. 740 test_params(num_params)

741 742

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

744 """The function for creating bounds for the mapping function. 745 746 @param param: The name of the parameter to return the bounds for. 747 @type param: str 748 @keyword spin_id: The spin identification string. This arg is unused. 749 @type spin_id: None or str 750 @return: The bounds for the parameter. 751 @rtype: list of len 2 of floats 752 """ 753 754 # tm. 755 if param == 'tm': 756 return [0, 10.0 * 1e-9] 757 758 # {Diso, Dx, Dy, Dz, Dpar, Dper}. 759 if param == 'Diso' or param == 'Dx' or param == 'Dy' or param == 'Dz' or param == 'Dpar' or param == 'Dper': 760 return [1e6, 1e7] 761 762 # Da. 763 if param == 'Da': 764 return [-3.0/2.0 * 1e7, 3.0 * 1e7] 765 766 # Dr. 767 if param == 'Dr': 768 return [0, 1] 769 770 # Dratio. 771 if param == 'Dratio': 772 return [1.0/3.0, 3.0] 773 774 # theta. 775 if param == 'theta': 776 return [0, pi] 777 778 # phi. 779 if param == 'phi': 780 return [0, 2*pi] 781 782 # alpha. 783 if param == 'alpha': 784 return [0, 2*pi] 785 786 # beta. 787 if param == 'beta': 788 return [0, pi] 789 790 # gamma. 791 if param == 'gamma': 792 return [0, 2*pi]

793 794

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

796 """Function for creating labels, tick locations, and tick values for an OpenDX map.""" 797 798 # Initialise. 799 labels = "{" 800 tick_locations = [] 801 tick_values = [] 802 n = len(params) 803 axis_incs = 5 804 loc_inc = inc / axis_incs 805 806 # Increment over the model parameters. 807 for i in range(n): 808 # Parameter conversion factors. 809 factor = return_conversion_factor(params[swap[i]]) 810 811 # Parameter units. 812 units = return_units(params[swap[i]]) 813 814 # Labels. 815 if units: 816 labels = labels + "\"" + params[swap[i]] + " (" + units + ")\"" 817 else: 818 labels = labels + "\"" + params[swap[i]] + "\"" 819 820 # Tick values. 821 vals = bounds[swap[i], 0] / factor 822 val_inc = (bounds[swap[i], 1] - bounds[swap[i], 0]) / (axis_incs * factor) 823 824 if i < n - 1: 825 labels = labels + " " 826 else: 827 labels = labels + "}" 828 829 # Tick locations. 830 string = "{" 831 val = 0.0 832 for j in range(axis_incs + 1): 833 string = string + " " + repr(val) 834 val = val + loc_inc 835 string = string + " }" 836 tick_locations.append(string) 837 838 # Tick values. 839 string = "{" 840 for j in range(axis_incs + 1): 841 string = string + "\"" + "%.2f" % vals + "\" " 842 vals = vals + val_inc 843 string = string + "}" 844 tick_values.append(string) 845 846 return labels, tick_locations, tick_values

847 848

849 -def return_conversion_factor(param):

850 """Function for returning the factor of conversion between different parameter units. 851 852 For example, the internal representation of tm is in seconds, whereas the external 853 representation is in nanoseconds, therefore this function will return 1e-9 for tm. 854 855 856 @param param: The name of the parameter to return the conversion factor for. 857 @type param: str 858 @return: The conversion factor. 859 @rtype: float 860 """ 861 862 # Get the object name. 863 object_name = return_data_name(param) 864 865 # tm (nanoseconds). 866 if object_name == 'tm': 867 return 1e-9 868 869 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 870 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 871 return 1e6 872 873 # Angles. 874 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 875 return (2.0*pi) / 360.0 876 877 # No conversion factor. 878 return 1.0

879 880

881 -def return_data_name(name):

882 """Return the parameter name. 883 884 @param name: The name of the parameter to return the name of. 885 @type name: str 886 @return: The parameter name. 887 @rtype: str 888 """ 889 890 # Enforce that the name must be a string. 891 if not isinstance(name, str): 892 raise RelaxStrError('name', name) 893 894 # Local tm. 895 if search('^tm$', name): 896 return 'tm' 897 898 # Diso. 899 if search('[Dd]iso', name): 900 return 'Diso' 901 902 # Da. 903 if search('[Dd]a', name): 904 return 'Da' 905 906 # Dr. 907 if search('[Dd]r$', name): 908 return 'Dr' 909 910 # Dx. 911 if search('[Dd]x', name): 912 return 'Dx' 913 914 # Dy. 915 if search('[Dd]y', name): 916 return 'Dy' 917 918 # Dz. 919 if search('[Dd]z', name): 920 return 'Dz' 921 922 # Dpar. 923 if search('[Dd]par', name): 924 return 'Dpar' 925 926 # Dper. 927 if search('[Dd]per', name): 928 return 'Dper' 929 930 # Dratio. 931 if search('[Dd]ratio', name): 932 return 'Dratio' 933 934 # alpha. 935 if search('^a$', name) or search('alpha', name): 936 return 'alpha' 937 938 # beta. 939 if search('^b$', name) or search('beta', name): 940 return 'beta' 941 942 # gamma. 943 if search('^g$', name) or search('gamma', name): 944 return 'gamma' 945 946 # theta. 947 if search('theta', name): 948 return 'theta' 949 950 # phi. 951 if search('phi', name): 952 return 'phi'

953 954 # User function documentation. 955 __return_data_name_doc__ = Desc_container("Diffusion tensor parameter string matching patterns") 956 table = uf_tables.add_table(label="table: diff data type patterns", caption="Diffusion tensor parameter string matching patterns.") 957 table.add_headings(["Data type", "Object name", "Patterns"]) 958 table.add_row(["Global correlation time - tm", "'tm'", "'^tm$'"]) 959 table.add_row(["Isotropic component of the diffusion tensor - Diso", "'Diso'", "'[Dd]iso'"]) 960 table.add_row(["Anisotropic component of the diffusion tensor - Da", "'Da'", "'[Dd]a'"]) 961 table.add_row(["Rhombic component of the diffusion tensor - Dr", "'Dr'", "'[Dd]r$'"]) 962 table.add_row(["Eigenvalue associated with the x-axis of the diffusion tensor - Dx", "'Dx'", "'[Dd]x'"]) 963 table.add_row(["Eigenvalue associated with the y-axis of the diffusion tensor - Dy", "'Dy'", "'[Dd]y'"]) 964 table.add_row(["Eigenvalue associated with the z-axis of the diffusion tensor - Dz", "'Dz'", "'[Dd]z'"]) 965 table.add_row(["Diffusion coefficient parallel to the major axis of the spheroid diffusion tensor - Dpar", "'Dpar'", "'[Dd]par'"]) 966 table.add_row(["Diffusion coefficient perpendicular to the major axis of the spheroid diffusion tensor - Dper", "'Dper'", "'[Dd]per'"]) 967 table.add_row(["Ratio of the parallel and perpendicular components of the spheroid diffusion tensor - Dratio", "'Dratio'", "'[Dd]ratio'"]) 968 table.add_row(["The first Euler angle of the ellipsoid diffusion tensor - alpha", "'alpha'", "'^a$' or 'alpha'"]) 969 table.add_row(["The second Euler angle of the ellipsoid diffusion tensor - beta", "'beta'", "'^b$' or 'beta'"]) 970 table.add_row(["The third Euler angle of the ellipsoid diffusion tensor - gamma", "'gamma'", "'^g$' or 'gamma'"]) 971 table.add_row(["The polar angle defining the major axis of the spheroid diffusion tensor - theta", "'theta'", "'theta'"]) 972 table.add_row(["The azimuthal angle defining the major axis of the spheroid diffusion tensor - phi", "'phi'", "'phi'"]) 973 __return_data_name_doc__.add_table(table.label) 974 975

976 -def return_units(param):

977 """Function for returning a string representing the parameters units. 978 979 For example, the internal representation of tm is in seconds, whereas the external 980 representation is in nanoseconds, therefore this function will return the string 981 'nanoseconds' for tm. 982 983 984 @param param: The name of the parameter to return the units for. 985 @type param: str 986 @return: The parameter units string. 987 @rtype: str 988 """ 989 990 # Get the object name. 991 object_name = return_data_name(param) 992 993 # tm (nanoseconds). 994 if object_name == 'tm': 995 return 'ns' 996 997 # Diso, Da, Dx, Dy, Dz, Dpar, Dper. 998 if object_name in ['Diso', 'Da', 'Dx', 'Dy', 'Dz', 'Dpar', 'Dper']: 999 return '1e6 1/s' 1000 1001 # Angles. 1002 if object_name in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1003 return 'deg'

1004 1005

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

1007 """Set the diffusion tensor parameters. 1008 1009 @keyword tensor: The diffusion tensor object. 1010 @type tensor: DiffTensorData instance 1011 @keyword value: The list of values to set the parameters to. 1012 @type value: list of float 1013 @keyword param: The list of parameter names. 1014 @type param: list of str 1015 """ 1016 1017 # Set up the diffusion tensor data if it doesn't exist. 1018 if not diff_data_exists(): 1019 raise RelaxNoTensorError('diffusion') 1020 1021 # Initialise. 1022 geo_params = [] 1023 geo_values = [] 1024 orient_params = [] 1025 orient_values = [] 1026 1027 # Loop over the parameters. 1028 for i in range(len(param)): 1029 # Get the object name. 1030 param[i] = return_data_name(param[i]) 1031 1032 # Unknown parameter. 1033 if not param[i]: 1034 raise RelaxUnknownParamError("diffusion tensor", param[i]) 1035 1036 # Default value. 1037 if value[i] == None: 1038 value[i] = default_value(param[i]) 1039 1040 # Geometric parameter. 1041 if param[i] in ['tm', 'Diso', 'Da', 'Dratio', 'Dper', 'Dpar', 'Dr', 'Dx', 'Dy', 'Dz']: 1042 geo_params.append(param[i]) 1043 geo_values.append(value[i]) 1044 1045 # Orientational parameter. 1046 if param[i] in ['theta', 'phi', 'alpha', 'beta', 'gamma']: 1047 orient_params.append(param[i]) 1048 orient_values.append(value[i]) 1049 1050 1051 # Spherical diffusion. 1052 ###################### 1053 1054 if cdp.diff_tensor.type == 'sphere': 1055 # Geometric parameters. 1056 ####################### 1057 1058 # A single geometric parameter. 1059 if len(geo_params) == 1: 1060 # The single parameter tm. 1061 if geo_params[0] == 'tm': 1062 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1063 1064 # The single parameter Diso. 1065 elif geo_params[0] == 'Diso': 1066 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1067 1068 # Cannot set the single parameter. 1069 else: 1070 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1071 1072 # More than one geometric parameters. 1073 elif len(geo_params) > 1: 1074 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1075 1076 1077 # Orientational parameters. 1078 ########################### 1079 1080 # ??? 1081 if len(orient_params): 1082 raise RelaxError("For spherical diffusion, the orientation parameters " + repr(orient_params) + " should not exist.") 1083 1084 1085 # Spheroidal diffusion. 1086 ####################### 1087 1088 elif cdp.diff_tensor.type == 'spheroid': 1089 # Geometric parameters. 1090 ####################### 1091 1092 # A single geometric parameter. 1093 if len(geo_params) == 1: 1094 # The single parameter tm. 1095 if geo_params[0] == 'tm': 1096 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1097 1098 # The single parameter Diso. 1099 elif geo_params[0] == 'Diso': 1100 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1101 1102 # The single parameter Da. 1103 elif geo_params[0] == 'Da': 1104 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 1105 1106 # The single parameter Dratio. 1107 elif geo_params[0] == 'Dratio': 1108 Dratio = geo_values[0] 1109 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * cdp.diff_tensor.tm * (Dratio + 2.0))) 1110 1111 # Cannot set the single parameter. 1112 else: 1113 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1114 1115 # Two geometric parameters. 1116 elif len(geo_params) == 2: 1117 # The geometric parameter set {tm, Da}. 1118 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1119 # The parameters. 1120 tm = geo_values[geo_params.index('tm')] 1121 Da = geo_values[geo_params.index('Da')] 1122 1123 # Set the internal parameter values. 1124 cdp.diff_tensor.set(param='tm', value=tm) 1125 cdp.diff_tensor.set(param='Da', value=Da) 1126 1127 # The geometric parameter set {Diso, Da}. 1128 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1129 # The parameters. 1130 Diso = geo_values[geo_params.index('Diso')] 1131 Da = geo_values[geo_params.index('Da')] 1132 1133 # Set the internal parameter values. 1134 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1135 cdp.diff_tensor.set(param='Da', value=Da) 1136 1137 # The geometric parameter set {tm, Dratio}. 1138 elif geo_params.count('tm') == 1 and geo_params.count('Dratio') == 1: 1139 # The parameters. 1140 tm = geo_values[geo_params.index('tm')] 1141 Dratio = geo_values[geo_params.index('Dratio')] 1142 1143 # Set the internal parameter values. 1144 cdp.diff_tensor.set(param='tm', value=tm) 1145 cdp.diff_tensor.set(param='Da', value=(Dratio - 1.0) / (2.0 * tm * (Dratio + 2.0))) 1146 1147 # The geometric parameter set {Dpar, Dper}. 1148 elif geo_params.count('Dpar') == 1 and geo_params.count('Dper') == 1: 1149 # The parameters. 1150 Dpar = geo_values[geo_params.index('Dpar')] 1151 Dper = geo_values[geo_params.index('Dper')] 1152 1153 # Set the internal parameter values. 1154 cdp.diff_tensor.set(param='tm', value=1.0 / (2.0 * (Dpar + 2.0*Dper))) 1155 cdp.diff_tensor.set(param='Da', value=Dpar - Dper) 1156 1157 # The geometric parameter set {Diso, Dratio}. 1158 elif geo_params.count('Diso') == 1 and geo_params.count('Dratio') == 1: 1159 # The parameters. 1160 Diso = geo_values[geo_params.index('Diso')] 1161 Dratio = geo_values[geo_params.index('Dratio')] 1162 1163 # Set the internal parameter values. 1164 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1165 cdp.diff_tensor.set(param='Da', value=3.0 * Diso * (Dratio - 1.0) / (Dratio + 2.0)) 1166 1167 # Unknown parameter combination. 1168 else: 1169 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1170 1171 # More than two geometric parameters. 1172 elif len(geo_params) > 2: 1173 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1174 1175 1176 # Orientational parameters. 1177 ########################### 1178 1179 # A single orientational parameter. 1180 if len(orient_params) == 1: 1181 # The single parameter theta. 1182 if orient_params[0] == 'theta': 1183 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 1184 1185 # The single parameter phi. 1186 elif orient_params[0] == 'phi': 1187 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 1188 1189 # Disallowed parameter. 1190 else: 1191 raise RelaxError("For spheroidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1192 1193 # Two orientational parameters. 1194 elif len(orient_params) == 2: 1195 # The orientational parameter set {theta, phi}. 1196 if orient_params.count('theta') == 1 and orient_params.count('phi') == 1: 1197 cdp.diff_tensor.set(param='theta', value=orient_values[orient_params.index('theta')]) 1198 cdp.diff_tensor.set(param='phi', value=orient_values[orient_params.index('phi')]) 1199 1200 # Unknown parameter combination. 1201 else: 1202 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1203 1204 # More than two orientational parameters. 1205 elif len(orient_params) > 2: 1206 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1207 1208 1209 # Ellipsoidal diffusion. 1210 ######################## 1211 1212 elif cdp.diff_tensor.type == 'ellipsoid': 1213 # Geometric parameters. 1214 ####################### 1215 1216 # A single geometric parameter. 1217 if len(geo_params) == 1: 1218 # The single parameter tm. 1219 if geo_params[0] == 'tm': 1220 cdp.diff_tensor.set(param='tm', value=geo_values[0]) 1221 1222 # The single parameter Diso. 1223 elif geo_params[0] == 'Diso': 1224 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * geo_values[0])) 1225 1226 # The single parameter Da. 1227 elif geo_params[0] == 'Da': 1228 cdp.diff_tensor.set(param='Da', value=geo_values[0]) 1229 1230 # The single parameter Dr. 1231 elif geo_params[0] == 'Dr': 1232 cdp.diff_tensor.set(param='Dr', value=geo_values[0]) 1233 1234 # Cannot set the single parameter. 1235 else: 1236 raise RelaxError("The geometric diffusion parameter " + repr(geo_params[0]) + " cannot be set.") 1237 1238 # Two geometric parameters. 1239 elif len(geo_params) == 2: 1240 # The geometric parameter set {tm, Da}. 1241 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1: 1242 # The parameters. 1243 tm = geo_values[geo_params.index('tm')] 1244 Da = geo_values[geo_params.index('Da')] 1245 1246 # Set the internal parameter values. 1247 cdp.diff_tensor.set(param='tm', value=tm) 1248 cdp.diff_tensor.set(param='Da', value=Da) 1249 1250 # The geometric parameter set {tm, Dr}. 1251 elif geo_params.count('tm') == 1 and geo_params.count('Dr') == 1: 1252 # The parameters. 1253 tm = geo_values[geo_params.index('tm')] 1254 Dr = geo_values[geo_params.index('Dr')] 1255 1256 # Set the internal parameter values. 1257 cdp.diff_tensor.set(param='tm', value=tm) 1258 cdp.diff_tensor.set(param='Dr', value=Dr) 1259 1260 # The geometric parameter set {Diso, Da}. 1261 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1: 1262 # The parameters. 1263 Diso = geo_values[geo_params.index('Diso')] 1264 Da = geo_values[geo_params.index('Da')] 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=Da) 1269 1270 # The geometric parameter set {Diso, Dr}. 1271 elif geo_params.count('Diso') == 1 and geo_params.count('Dr') == 1: 1272 # The parameters. 1273 Diso = geo_values[geo_params.index('Diso')] 1274 Dr = geo_values[geo_params.index('Dr')] 1275 1276 # Set the internal parameter values. 1277 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1278 cdp.diff_tensor.set(param='Dr', value=Dr) 1279 1280 # The geometric parameter set {Da, Dr}. 1281 elif geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1282 # The parameters. 1283 Da = geo_values[geo_params.index('Da')] 1284 Dr = geo_values[geo_params.index('Dr')] 1285 1286 # Set the internal parameter values. 1287 cdp.diff_tensor.set(param='Da', value=Da) 1288 cdp.diff_tensor.set(param='Da', value=Dr) 1289 1290 # Unknown parameter combination. 1291 else: 1292 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1293 1294 # Three geometric parameters. 1295 elif len(geo_params) == 3: 1296 # The geometric parameter set {tm, Da, Dr}. 1297 if geo_params.count('tm') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1298 # The parameters. 1299 tm = geo_values[geo_params.index('tm')] 1300 Da = geo_values[geo_params.index('Da')] 1301 Dr = geo_values[geo_params.index('Dr')] 1302 1303 # Set the internal parameter values. 1304 cdp.diff_tensor.set(param='tm', value=tm) 1305 cdp.diff_tensor.set(param='Da', value=Da) 1306 cdp.diff_tensor.set(param='Dr', value=Dr) 1307 1308 # The geometric parameter set {Diso, Da, Dr}. 1309 elif geo_params.count('Diso') == 1 and geo_params.count('Da') == 1 and geo_params.count('Dr') == 1: 1310 # The parameters. 1311 Diso = geo_values[geo_params.index('Diso')] 1312 Da = geo_values[geo_params.index('Da')] 1313 Dr = geo_values[geo_params.index('Dr')] 1314 1315 # Set the internal parameter values. 1316 cdp.diff_tensor.set(param='tm', value=1.0 / (6.0 * Diso)) 1317 cdp.diff_tensor.set(param='Da', value=Da) 1318 cdp.diff_tensor.set(param='Dr', value=Dr) 1319 1320 # The geometric parameter set {Dx, Dy, Dz}. 1321 elif geo_params.count('Dx') == 1 and geo_params.count('Dy') == 1 and geo_params.count('Dz') == 1: 1322 # The parameters. 1323 Dx = geo_values[geo_params.index('Dx')] 1324 Dy = geo_values[geo_params.index('Dy')] 1325 Dz = geo_values[geo_params.index('Dz')] 1326 1327 # Set the internal tm value. 1328 if Dx + Dy + Dz == 0.0: 1329 cdp.diff_tensor.set(param='tm', value=1e99) 1330 else: 1331 cdp.diff_tensor.set(param='tm', value=0.5 / (Dx + Dy + Dz)) 1332 1333 # Set the internal Da value. 1334 cdp.diff_tensor.set(param='Da', value=Dz - 0.5*(Dx + Dy)) 1335 1336 # Set the internal Dr value. 1337 if cdp.diff_tensor.Da == 0.0: 1338 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) * 1e99) 1339 else: 1340 cdp.diff_tensor.set(param='Dr', value=(Dy - Dx) / (2.0*cdp.diff_tensor.Da)) 1341 1342 # Unknown parameter combination. 1343 else: 1344 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1345 1346 1347 # More than three geometric parameters. 1348 elif len(geo_params) > 3: 1349 raise RelaxUnknownParamCombError('geometric parameter set', geo_params) 1350 1351 1352 # Orientational parameters. 1353 ########################### 1354 1355 # A single orientational parameter. 1356 if len(orient_params) == 1: 1357 # The single parameter alpha. 1358 if orient_params[0] == 'alpha': 1359 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1360 1361 # The single parameter beta. 1362 elif orient_params[0] == 'beta': 1363 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1364 1365 # The single parameter gamma. 1366 elif orient_params[0] == 'gamma': 1367 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1368 1369 # Disallowed parameter. 1370 else: 1371 raise RelaxError("For ellipsoidal diffusion, the orientation parameter " + repr(orient_params) + " cannot be set.") 1372 1373 # Two orientational parameters. 1374 elif len(orient_params) == 2: 1375 # The orientational parameter set {alpha, beta}. 1376 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1377 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1378 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1379 1380 # The orientational parameter set {alpha, gamma}. 1381 if orient_params.count('alpha') == 1 and orient_params.count('gamma') == 1: 1382 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1383 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1384 1385 # The orientational parameter set {beta, gamma}. 1386 if orient_params.count('beta') == 1 and orient_params.count('gamma') == 1: 1387 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1388 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1389 1390 # Unknown parameter combination. 1391 else: 1392 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1393 1394 # Three orientational parameters. 1395 elif len(orient_params) == 3: 1396 # The orientational parameter set {alpha, beta, gamma}. 1397 if orient_params.count('alpha') == 1 and orient_params.count('beta') == 1: 1398 cdp.diff_tensor.set(param='alpha', value=orient_values[orient_params.index('alpha')]) 1399 cdp.diff_tensor.set(param='beta', value=orient_values[orient_params.index('beta')]) 1400 cdp.diff_tensor.set(param='gamma', value=orient_values[orient_params.index('gamma')]) 1401 1402 # Unknown parameter combination. 1403 else: 1404 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1405 1406 # More than three orientational parameters. 1407 elif len(orient_params) > 3: 1408 raise RelaxUnknownParamCombError('orientational parameter set', orient_params) 1409 1410 1411 # Fold the angles in. 1412 ##################### 1413 1414 if orient_params: 1415 fold_angles()

1416 1417 # User function documentation. 1418 __set_doc__ = Desc_container("Diffusion tensor set details") 1419 __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.") 1420 __set_doc__.add_paragraph("The units of the parameters are:") 1421 __set_doc__.add_list_element("Inverse seconds for tm.") 1422 __set_doc__.add_list_element("Seconds for Diso, Da, Dx, Dy, Dz, Dpar, Dper.") 1423 __set_doc__.add_list_element("Unitless for Dratio and Dr.") 1424 __set_doc__.add_list_element("Radians for all angles (alpha, beta, gamma, theta, phi).") 1425 __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") 1426 __set_doc__.add_list_element("{tm, Da},") 1427 __set_doc__.add_list_element("{Diso, Da},") 1428 __set_doc__.add_list_element("{tm, Dratio},") 1429 __set_doc__.add_list_element("{Dpar, Dper},") 1430 __set_doc__.add_list_element("{Diso, Dratio},") 1431 __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") 1432 __set_doc__.add_list_element("{tm, Da, Dr},") 1433 __set_doc__.add_list_element("{Diso, Da, Dr},") 1434 __set_doc__.add_list_element("{Dx, Dy, Dz},") 1435 __set_doc__.add_paragraph("where any number of the orientational parameters, alpha, beta, or gamma can be additionally supplied.") 1436 1437

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

1439 """Function for setting up a spherical diffusion tensor. 1440 1441 @param params: The diffusion tensor parameter. 1442 @type params: float 1443 @param time_scale: The correlation time scaling value. 1444 @type time_scale: float 1445 @param param_types: The type of parameter supplied. If 0, then the parameter is tm. If 1, then 1446 the parameter is Diso. 1447 @type param_types: int 1448 """ 1449 1450 # The diffusion type. 1451 cdp.diff_tensor.set_type('sphere') 1452 1453 # tm. 1454 if param_types == 0: 1455 # Scaling. 1456 tm = params * time_scale 1457 1458 # Set the parameters. 1459 set(value=[tm], param=['tm']) 1460 1461 # Diso 1462 elif param_types == 1: 1463 # Scaling. 1464 Diso = params * d_scale 1465 1466 # Set the parameters. 1467 set(value=[Diso], param=['Diso']) 1468 1469 # Unknown parameter combination. 1470 else: 1471 raise RelaxUnknownParamCombError('param_types', param_types)

1472 1473

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

1475 """Function for setting up a spheroidal diffusion tensor. 1476 1477 @param params: The diffusion tensor parameter. 1478 @type params: float 1479 @param time_scale: The correlation time scaling value. 1480 @type time_scale: float 1481 @param d_scale: The diffusion tensor eigenvalue scaling value. 1482 @type d_scale: float 1483 @param angle_units: The units for the angle parameters which can be either 'deg' or 'rad'. 1484 @type angle_units: str 1485 @param param_types: The type of parameters supplied. These correspond to 0: {tm, Da, theta, 1486 phi}, 1: {Diso, Da, theta, phi}, 2: {tm, Dratio, theta, phi}, 3: {Dpar, 1487 Dper, theta, phi}, 4: {Diso, Dratio, theta, phi}. 1488 @type param_types: int 1489 @param spheroid_type: A string which, if supplied together with spheroid parameters, will 1490 restrict the tensor to either being 'oblate' or 'prolate'. 1491 @type spheroid_type: str 1492 """ 1493 1494 # The diffusion type. 1495 cdp.diff_tensor.set_type('spheroid') 1496 1497 # Spheroid diffusion type. 1498 allowed_types = [None, 'oblate', 'prolate'] 1499 if spheroid_type not in allowed_types: 1500 raise RelaxError("The 'spheroid_type' argument " + repr(spheroid_type) + " should be 'oblate', 'prolate', or None.") 1501 cdp.diff_tensor.set(param='spheroid_type', value=spheroid_type) 1502 1503 # (tm, Da, theta, phi). 1504 if param_types == 0: 1505 # Unpack the tuple. 1506 tm, Da, theta, phi = params 1507 1508 # Scaling. 1509 tm = tm * time_scale 1510 Da = Da * d_scale 1511 1512 # Set the parameters. 1513 set(value=[tm, Da], param=['tm', 'Da']) 1514 1515 # (Diso, Da, theta, phi). 1516 elif param_types == 1: 1517 # Unpack the tuple. 1518 Diso, Da, theta, phi = params 1519 1520 # Scaling. 1521 Diso = Diso * d_scale 1522 Da = Da * d_scale 1523 1524 # Set the parameters. 1525 set(value=[Diso, Da], param=['Diso', 'Da']) 1526 1527 # (tm, Dratio, theta, phi). 1528 elif param_types == 2: 1529 # Unpack the tuple. 1530 tm, Dratio, theta, phi = params 1531 1532 # Scaling. 1533 tm = tm * time_scale 1534 1535 # Set the parameters. 1536 set(value=[tm, Dratio], param=['tm', 'Dratio']) 1537 1538 # (Dpar, Dper, theta, phi). 1539 elif param_types == 3: 1540 # Unpack the tuple. 1541 Dpar, Dper, theta, phi = params 1542 1543 # Scaling. 1544 Dpar = Dpar * d_scale 1545 Dper = Dper * d_scale 1546 1547 # Set the parameters. 1548 set(value=[Dpar, Dper], param=['Dpar', 'Dper']) 1549 1550 # (Diso, Dratio, theta, phi). 1551 elif param_types == 4: 1552 # Unpack the tuple. 1553 Diso, Dratio, theta, phi = params 1554 1555 # Scaling. 1556 Diso = Diso * d_scale 1557 1558 # Set the parameters. 1559 set(value=[Diso, Dratio], param=['Diso', 'Dratio']) 1560 1561 # Unknown parameter combination. 1562 else: 1563 raise RelaxUnknownParamCombError('param_types', param_types) 1564 1565 # Convert the angles to radians. 1566 if angle_units == 'deg': 1567 print("Converting the angles from degrees to radian units.") 1568 theta = (theta / 360.0) * 2.0 * pi 1569 phi = (phi / 360.0) * 2.0 * pi 1570 1571 # Set the orientational parameters. 1572 set(value=[theta, phi], param=['theta', 'phi'])

1573 1574

1575 -def test_params(num_params):

1576 """Function for testing the validity of the input parameters.""" 1577 1578 # An allowable error to account for machine precision, optimisation quality, etc. 1579 error = 1e-4 1580 1581 # tm. 1582 tm = cdp.diff_tensor.tm 1583 if tm <= 0.0 or tm > 1e-6: 1584 raise RelaxError("The tm value of " + repr(tm) + " should be between zero and one microsecond.") 1585 1586 # Spheroid. 1587 if num_params == 4: 1588 # Parameters. 1589 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1590 Da = cdp.diff_tensor.Da 1591 1592 # Da. 1593 if Da < (-1.5*Diso - error*Diso) or Da > (3.0*Diso + error*Diso): 1594 raise RelaxError("The Da value of " + repr(Da) + " should be between -3/2 * Diso and 3Diso.") 1595 1596 # Ellipsoid. 1597 if num_params == 6: 1598 # Parameters. 1599 Diso = 1.0 / (6.0 * cdp.diff_tensor.tm) 1600 Da = cdp.diff_tensor.Da 1601 Dr = cdp.diff_tensor.Dr 1602 1603 # Da. 1604 if Da < (0.0 - error*Diso) or Da > (3.0*Diso + error*Diso): 1605 raise RelaxError("The Da value of " + repr(Da) + " should be between zero and 3Diso.") 1606 1607 # Dr. 1608 if Dr < (0.0 - error) or Dr > (1.0 + error): 1609 raise RelaxError("The Dr value of " + repr(Dr) + " should be between zero and one.")

1610 1611

1612 -def unit_axes():

1613 """Function for calculating the unit axes of the diffusion tensor. 1614 1615 Spheroid 1616 ======== 1617 1618 The unit Dpar vector is:: 1619 1620 | sin(theta) * cos(phi) | 1621 Dpar = | sin(theta) * sin(phi) | 1622 | cos(theta) | 1623 1624 1625 Ellipsoid 1626 ========= 1627 1628 The unit Dx vector is:: 1629 1630 | -sin(alpha) * sin(gamma) + cos(alpha) * cos(beta) * cos(gamma) | 1631 Dx = | -sin(alpha) * cos(gamma) - cos(alpha) * cos(beta) * sin(gamma) | 1632 | cos(alpha) * sin(beta) | 1633 1634 The unit Dy vector is:: 1635 1636 | cos(alpha) * sin(gamma) + sin(alpha) * cos(beta) * cos(gamma) | 1637 Dy = | cos(alpha) * cos(gamma) - sin(alpha) * cos(beta) * sin(gamma) | 1638 | sin(alpha) * sin(beta) | 1639 1640 The unit Dz vector is:: 1641 1642 | -sin(beta) * cos(gamma) | 1643 Dz = | sin(beta) * sin(gamma) | 1644 | cos(beta) | 1645 1646 """ 1647 1648 # Spheroid. 1649 if cdp.diff_tensor.type == 'spheroid': 1650 # Initialise. 1651 Dpar = zeros(3, float64) 1652 1653 # Trig. 1654 sin_theta = sin(cdp.diff_tensor.theta) 1655 cos_theta = cos(cdp.diff_tensor.theta) 1656 sin_phi = sin(cdp.diff_tensor.phi) 1657 cos_phi = cos(cdp.diff_tensor.phi) 1658 1659 # Unit Dpar axis. 1660 Dpar[0] = sin_theta * cos_phi 1661 Dpar[1] = sin_theta * sin_phi 1662 Dpar[2] = cos_theta 1663 1664 # Return the vector. 1665 return Dpar 1666 1667 # Ellipsoid. 1668 if cdp.diff_tensor.type == 'ellipsoid': 1669 # Initialise. 1670 Dx = zeros(3, float64) 1671 Dy = zeros(3, float64) 1672 Dz = zeros(3, float64) 1673 1674 # Trig. 1675 sin_alpha = sin(cdp.diff_tensor.alpha) 1676 cos_alpha = cos(cdp.diff_tensor.alpha) 1677 sin_beta = sin(cdp.diff_tensor.beta) 1678 cos_beta = cos(cdp.diff_tensor.beta) 1679 sin_gamma = sin(cdp.diff_tensor.gamma) 1680 cos_gamma = cos(cdp.diff_tensor.gamma) 1681 1682 # Unit Dx axis. 1683 Dx[0] = -sin_alpha * sin_gamma + cos_alpha * cos_beta * cos_gamma 1684 Dx[1] = -sin_alpha * cos_gamma - cos_alpha * cos_beta * sin_gamma 1685 Dx[2] = cos_alpha * sin_beta 1686 1687 # Unit Dy axis. 1688 Dx[0] = cos_alpha * sin_gamma + sin_alpha * cos_beta * cos_gamma 1689 Dx[1] = cos_alpha * cos_gamma - sin_alpha * cos_beta * sin_gamma 1690 Dx[2] = sin_alpha * sin_beta 1691 1692 # Unit Dz axis. 1693 Dx[0] = -sin_beta * cos_gamma 1694 Dx[1] = sin_beta * sin_gamma 1695 Dx[2] = cos_beta 1696 1697 # Return the vectors. 1698 return Dx, Dy, Dz

1699

Source Code for Module pipe_control.diffusion_tensor